CN103942408A - Annual erosion sediment yield model calculation method of mesoscale drainage basin of loess plateau - Google Patents

Abstract

The invention relates to an annual erosion sediment yield model calculation method of a mesoscale drainage basin of the loess plateau. The method includes the concrete steps that the drainage basin is divided into a plurality of units, and each unit approximately represents one slope surface; through the slope surface soil erosion test observation and research result and knowledge, relations between all factors and the erosion amount, namely unit models, are established based on slope surface soil erosion test observation data; the erosion amount of each unit is calculated through the unit models; with a spatial analysis and spatial statistical method, the overall conditions of soil erosion of the drainage basin are evaluated and analyzed through clustering operation. According to the method, the GIS spatial analysis function is fully utilized, the output results of the method can give out the erosion intensities of all portions in the drainage basin, and thus not only can the requirements of channel water conservancy engineering design be met, but also the demand for arrangement and planning of water and soil conservation measures can be met.

Description

Mesoscale basin, loess plateau year Model of Erosion and Sediment Production computing method

Technical field

The present invention relates to Territorial Soil Erosion field, relate in particular to mesoscale basin, a kind of loess plateau year Model of Erosion and Sediment Production computing method.

Background technology

In prior art, application and improvement for the model system of Water and soil conservation planning and design, water-and-soil conservation measures performance evaluation, a mesoscale flow survey of soil and water loss and monitoring on regional scale are significant to soil test.

Lumped model, sets up factor of soil erosion database based on GIS, at thematic maps, by space statistical analysis, obtains one or several statistical characteristicss of each factor, in order to characterize on Watershed Scale because of subcharacter.Then the sedimentary loading data that all Erosion factor values and hydrologic observation obtained are carried out correlation analysis, set up an empirical model, complete the analog computation of watershed Soil Erosion.The advantage of this thinking is, can to model, carry out calibration and checking with hydrographic data (sedimentary loading), and model is set up and be take statistical study as basis, and the course of work is relatively simple; Model Output rusults can meet the demand of raceway groove Water Conservancy engineering design substantially.

Shortcoming is to embody the requirement of project (the husky empirical model of the annual output of research and development based on GIS); Operation result is mainly the form about basin soil erosion general characteristic, can not meet the demand that water-and-soil conservation measures is laid; Need to take than relatively large hydrographic data as basis, collect more difficult; Simultaneously because not for the field data of certain specific single-factor (other factors is constant), thereby determine that the algorithm of each factor remains at certain difficulty, the transplanting to another one basin is restricted by a basin also therefore to make model.

In view of above-mentioned defect, creator of the present invention has obtained this creation finally through long research and practice.

Summary of the invention

The object of the present invention is to provide mesoscale basin, a kind of loess plateau year Model of Erosion and Sediment Production computing method, in order to overcome above-mentioned technological deficiency.

For achieving the above object, the invention provides mesoscale basin, a kind of loess plateau year Model of Erosion and Sediment Production computing method, its detailed process is:

Step 1, is divided into some unit basin, each unit is approximate represent one domatic;

Step 2, utilizes slope soil erosion test observation achievement in research and understanding, based on slope soil erosion test observation data, sets up relation between each factor and erosion amount, i.e. model of element;

Step 3, completes the calculating to each unit erosion amount with described model of element;

Step 4, utilizes spatial analysis and spatial statistics method, by polymerization computing, completes the evaluate and analyze of watershed soil erosion overall state;

In above-mentioned steps 2, the computation process of obtaining each factor is as follows:

Step 21, calculates website climatic factor R1;

Step 22, calculates climatic factor and calculates interpolation R2;

Step 23, calculates plant measures factor B, with land-use map and vegetation cover degree figure, calculates certain region plant measures factor B;

Step 24, computational engineering measure factor E, utilizes statistics to calculate the engineering measure factor E value in each region;

Step 25, calculates tillage control measure factor T, and tillage control measure factor values T reduces soil losses according to Contour farming under different gradient condition and determines.

Further, in above-mentioned steps 21, calculate the rainfall erosivity that website climatic factor R1 comprises one or more websites, its rainfall amount according to respective site calculates;

Computing formula is as shown in following (1) and (2):

$F=\left(\mathrm{\Σ}{P_i}^2\right)/P---\left(1\right)$

R=αF ^β (2)

In formula, P is average annual rainfall (mm); Pi is the average rainfall (mm) of the i month; R is average rainfall erosivity (MJmmhm for many years ^-2h ^-1a ^-1); α=0.3589, β=1.9462; F exponent m m size is relevant with the season distribution of mean annual rainfall P, and span is between P12-1～P.

Further, the detailed process of calculating plant measures factor B in above-mentioned steps 23 is:

Step 231, obtains land-use map;

Step 232, calculates vegetation cover degree figure;

Step 233, calculates annual vegetation cover degree;

Step 234, calculates plant measures factor B.

Further, in above-mentioned steps 232, utilize NDVI to calculate vegetation cover degree according to following formula (3),

$F=100\×\frac{\mathrm{float}\left(\mathrm{NDVI}\right)-\mathrm{NDV}{I}_{\min }}{{\mathrm{NDVI}}_{\max }-{\mathrm{NDVI}}_{\min }}---\left(3\right)$

In formula, F is vegetation cover degree, and NDVI is normalized differential vegetation index, and NDVImin and N DVImax are respectively without the N DVI value in vegetation area and the NDVI value of the good covering area of vegetation; Test basis is upper at the scene, in conjunction with the typically sampling of class of remote sensing images, determines the NDVI threshold value of the highest and minimum vegetation cover degree; The result of sampling is carried out to statistical with corresponding NDVI value, and when NDVI value >=NDVImax, cover degree can be considered 100%, and when NDVI value≤NDVImin, cover degree can be considered 0%.

Further, in above-mentioned steps 232, the process that N DVI data are calculated is:

For less basin, utilize TM data to calculate N DVI, refer to formula (4),

$\mathrm{NDVI}=\frac{\mathrm{TM}4-\mathrm{TM}3}{\mathrm{TM}4+\mathrm{TM}3}---\left(4\right)$

In formula, NDVI is required, normalized differential vegetation index; TM3, TM4 are respectively near infrared and the infrared band of remote sensing image;

For compared with large watershed, directly obtain.

Further, in above-mentioned steps 233, while calculating annual vegetation cover degree, by each month flood season vegetation cover degree value and flood season each month rainfall erosivity surface be weighted on average, can obtain the average vegetation cover degree of each month in flood season, computing formula is as shown in the formula (5),

$F_r=\frac{\underset{m=5}{\overset{9}{\mathrm{\Σ}}}{F}_m\×R_m}{\underset{m=5}{\overset{9}{\mathrm{\Σ}}}{R}_m}---\left(5\right)$

In formula: Fm is for pressing the vegetation cover degree of each month in flood season of above-mentioned (3) formula calculating, and Rm is the value on the rainfall erosivity surface of each month in flood season; Fr is cover degree and the weighted mean value of rainfall erosivity in flood season.

Further, in above-mentioned steps 24, utilize statistics to calculate the engineering measure factor E value in each region; Refer to shown in formula (6),

$E=(1-\frac{S_t}{S}\×\mathrm{\α})(1-\frac{S_d}{S}\×\mathrm{\β})(1-\frac{\mathrm{\λ}\×N_{d1}+\mathrm{\ϵ}\×N_{d2}}{A\×S})---\left(6\right)$

In formula, S _tfor terraced fields area, S _dfor silt arrester control area, S is the soil total area, α, and β is respectively the husky coefficient of subtracting of terraced fields and silt arrester, is respectively 0.763 and 1; N _d1, N _d2the quantity that is respectively check dam, check dam, unit is seat, and λ and ε are respectively the sediment trapping quota of check dam and check dam, are respectively 1000t/ seat and 100t/ seat, and A is zone leveling soil erosion modulus, and unit is t/km ².

Further, in above-mentioned steps 2, also obtain ditch eclipse factor G, to have or not rainfall data, comprise two kinds of algorithms,

(1) there is the calculating of rainfall amount data: the relation of ditch eclipse factor and ground inclination and characteristics of rainfall, shown in (9),

G=1+β (9)

In formula, G is shallow ridges coefficient, dimensionless; β is correction coefficient, affected by time rainfall, water catchment area and the gradient;

(2) without the calculating of rainfall data: when without rainfall data in detail, available following formula (11) estimation annual G value:

G=1+1.60sin(θ-15) (11)。

Further,

Above-mentioned while obtaining ditch eclipse factor G, when ground inclination is less than shallow gully erosion generation critical grade, during without shallow gully erosion, β=0; When shallow ridges generation critical grade is 15 °, β calculating formula is following formula (10),

$\mathrm{\β}=\left(\frac{\mathrm{\θ}-15}{15}\right)[1.003{\left({\mathrm{RI}}_{30}\right)}^{0.103}-1]---\left(10\right)$

In formula: θ be ground inclination (°); R is time rainfall amount (mm); I ₃₀for inferior rainfall maximum 30-min rainfall intensity (mm/min).

Further, in above-mentioned steps 4, the analytic process of the overall shape of the watershed soil erosion is as described below:

Average annual erosion intensity At: refer to calculate the average soil erosion intensity in basin or region, concrete computation process is as following formula (12),

$A_t=\left(\underset{c=1}{\overset{n}{\mathrm{\Σ}}}(A_c\×s_c)\right)/\underset{c=1}{\overset{n}{\mathrm{\Σ}}}{s}_c---\left(12\right)$

In formula: A _cyear soil erosion modulus (t/km for computing unit ²a), A _tfor calculating the annual soil erosion modulus (t/km in basin ²a), S _cfor computing unit area;

Produce husky total amount Et per year: refer to calculate the year Erosion and Sediment Production total amount in basin or region, concrete computation process is as following formula (13),

$E_t=A_t\×S_t=\underset{c=1}{\overset{n}{\mathrm{\Σ}}}{A}_c\×s_c---\left(13\right)$

In formula: E _tfor calculating the year sediment production (t/a) in basin, all the other same formula (12);

Erosion grade index Ie: for a parameter of basin or region erosion intensity grade general status; Concrete computation process is as following formula (14)-(16),

I=I ₁×a+I ₂×b (14)

$I_1=\underset{10}{\overset{16}{\mathrm{\Σ}}}s\×c_{\mathrm{wat}}\×{10}^{2\×(C_{\mathrm{wat}}-10)}---\left(15\right)$

$I_2=\underset{20}{\overset{26}{\mathrm{\Σ}}}s\×c_{\mathrm{wib}}\×{10}^{2(C_{\mathrm{win}}-20)}---\left(16\right)$

In formula: I is comprehensive erosion intensity index; I ₁for intensity of water erosion index; I ₂for intensity of wind erosion index; A, b are respectively water erosion and the shared drainage area ratio of wind erosion in basin; S is each erosion intensity grade area and drainage area ratio in basin; C is erosion intensity grade;

Simple eigenvalue St: carry out the most basic statistical indicator and calculate;

Frequency curve Fi: to zoning or each strength grade frequency statistics of basin.

Beneficial effect of the present invention is compared with prior art: the present invention takes full advantage of GIS spatial analysis functions, the result of its output (soil erosion intensity map) can provide the erosion intensity at Nei Ge position, basin, thereby both can meet the requirement of raceway groove Water Conservancy engineering design, also can meet water-and-soil conservation measures and lay the demand with planning; Meanwhile, the model that this method is set up has absorbed the achievement, the particularly model such as USLE, CSLE of slope soil erosion model preferably, than being easier to, promotes.

Accompanying drawing explanation

Fig. 1 is the process flow diagram of the computing method of mesoscale basin, loess plateau of the present invention year Model of Erosion and Sediment Production;

Fig. 2 is the process flow diagram that the present invention obtains the computation process of each factor.

embodiment

Below in conjunction with accompanying drawing, to the present invention is above-mentioned, be described in more detail with other technical characterictic and advantage.

The computing method of mesoscale basin, loess plateau of the present invention year Model of Erosion and Sediment Production, based on distributed year Model of Erosion and Sediment Production, the computation process of this model was:

Step 1, is divided into some unit basin, each unit is approximate represent one domatic;

Step 2, utilizes slope soil erosion test observation achievement in research and understanding, based on slope soil erosion test observation data, sets up relation between each factor and erosion amount---model of element;

Step 3, completes the calculating to each unit erosion amount with model of element;

Step 4, utilizes spatial analysis and spatial statistics method, by polymerization computing, completes the evaluate and analyze of watershed soil erosion overall state.

In the present invention, the factor relating in above-mentioned steps 2 comprises single station climatic factor R1, climatic factor interpolation R2, plant measures B, engineering measure E, tillage control measure T.

In above-mentioned steps 4, watershed Soil Erosion Assessment comprises: potential erosion 1 (A01), potential erosion 1 (A02), basin number of dropouts (A1), basin number of dropouts (A2), water erosion grade (A31), wind erosion class (A32).

In above-mentioned steps 4, the result of statistical study comprises: average annual erosion intensity (At), produce husky total amount (Et) per year, erosion grade index (Ie), simple eigenvalue (St), frequency curve (Fi) and frequency curve (Ff).

Concrete, in above-mentioned steps 2, the computation process of obtaining each factor is as follows:

Step 21, calculates website climatic factor R1, and this website climatic factor comprises the rainfall erosivity of one or more websites, and its rainfall amount according to respective site calculates;

In the present embodiment, rainfall amount data interactive is inputted, also can be imported from external data; The external data call format importing is the ascii text file of CSV symbol, or the .CVS file of Excel derivation; Require wherein the first behavior month, first classifies site name as.Refer to shown in table 1, it is each month precipitation data statistic of each website.

Each month precipitation data statistic of each website of table 1

Website/month	January	February	March	April	May	June	July	August	September	October	November	Dec
													Website 1	1.3	2.2	10.8	16.2	12.6	71.2	48.2	98.8	73	54.2	11.3	0
Website 2	0	5.5	18.3	8.6	18.1	122	58.6	95.3	33.9	11.3	8.2	5.1
													Website 3	3.8	4.5	16.3	12	34.5	21.6	156	24.2	38.6	2.1	23.8	0
Website 4	5.3	5.2	0	2.7	51.6	17.1	61.7	46	16.5	11.8	1.8	1.9
													Website 5	1.2	2.2	11.2	16.2	18.2	71.2	48.2	98.8	56.9	54.2	11.3	0
Website 6	5.4	4.5	3.4	2.7	51.6	17.1	62.5	36	16.5	11.8	1.8	1.9
													Website 7	0.23	0.56	0.56	0.56	0.56	0.56	0.56	0.56	0.56	0.56	0.56	0.56
Website 8	1.3	2.2	10.8	16.2	12.6	71.2	48.2	98.8	71.9	54.2	11.3	0
													Website 9	0	5.5	18.3	8.6	18.1	122	58.6	95.3	33.9	11.3	8.2	5.1
Website 10	3.8	4.5	16.3	12	34.5	21.6	156	24.2	38.6	2.1	23.8	0
													Website 11	5.3	5.2	0	2.7	51.6	17.1	61.7	46	16.5	11.8	1.8	1.9
Website 12	1.2	2.2	11.2	16.2	18.2	71.2	48.2	98.8	56.9	54.2	11.3	0
													Website 13	5.4	4.5	3.4	2.7	51.6	17.1	62.5	36	16.5	11.8	1.8	1.9

Computing formula is as shown in following (1) and (2):

$F=\left(\mathrm{\Σ}{P_i}^2\right)/P---\left(1\right)$

R=αF ^β (2)

In formula, P is average annual rainfall (mm); Pi is the average rainfall (mm) of the i month; R is average rainfall erosivity (MJmmhm for many years ^-2h ^-1a ^-1); α=0.3589, β=1.9462.F exponent m m size is relevant with the season distribution of mean annual rainfall P, and span is between P12-1～P.

Step 22, calculates climatic factor and calculates interpolation R2; Website climatic factor R1 according to above-mentioned calculating gained, can utilize the methods such as IDW, Kriging and Spline to complete interpolation, and the data that regulation participates in interpolation can not be less than 3 points.

Step 23, calculates plant measures factor B; The present invention, with land-use map and vegetation cover degree figure, calculates certain region plant measures factor B; Detailed process is:

Step 231, obtains land-use map: directly utilize ready-made data;

Step 232, calculates vegetation cover degree figure: if there is ready-made vegetation cover degree figure, directly call; Otherwise, utilize N DVI to calculate according to following formula (3),

$F=100\×\frac{\mathrm{float}\left(\mathrm{NDVI}\right)-\mathrm{NDV}{I}_{\min }}{{\mathrm{NDVI}}_{\max }-{\mathrm{NDVI}}_{\min }}---\left(3\right)$

In formula, F is vegetation cover degree, and NDVI is normalized differential vegetation index, and NDVImin and NDVImax are respectively without the NDVI value in vegetation area and the NDVI value of the good covering area of vegetation; Test basis is upper at the scene, in conjunction with the typically sampling of class (sand ground, the high forest land etc. that covers) of remote sensing images, determines the N DVI threshold value of the highest and minimum vegetation cover degree.The result of sampling is carried out to statistical with corresponding N DVI value, and when NDVI value >=NDVImax, cover degree can be considered 100%, and when NDVI value≤NDVImin, cover degree can be considered 0%.In the present embodiment, the minimum value of NDVI and maximal value get respectively 54 and 230.8.

In above-mentioned formula (3), the process that N DVI data are calculated is:

For less basin (area <=1000km2), can utilize TM data to calculate NDVI, refer to formula (4),

$\mathrm{NDVI}=\frac{\mathrm{TM}4-\mathrm{TM}3}{\mathrm{TM}4+\mathrm{TM}3}---\left(4\right)$

In formula, NDVI is required, normalized differential vegetation index; TM3, TM4 are respectively near infrared and the infrared band of remote sensing image.

For compared with large watershed (area >1000km2), can directly obtain.Download the NDVI data product of MODIS or SPOTVEGETATION.

Step 233, calculates annual vegetation cover degree; By each month flood season cover degree value and flood season each month rainfall erosivity surface be weighted on average, can obtain the average vegetation cover degree of each month in flood season.Computing formula is as shown in the formula (5),

$F_r=\frac{\underset{m=5}{\overset{9}{\mathrm{\&Sigma;}}}{F}_m\&times;R_m}{\underset{m=5}{\overset{9}{\mathrm{\&Sigma;}}}{R}_m}---\left(5\right)$

In formula: Fm is for pressing the vegetation cover degree of each month in flood season of above-mentioned (3) formula calculating, and Rm is the value on the rainfall erosivity surface of each month in flood season; Fr is cover degree and the weighted mean value of rainfall erosivity in flood season.

Step 234, calculates plant measures factor B; According to land-use map and vegetation cover degree, calculate.In the present invention, the code of Land-Use is: 11 arable lands, slope, 12 flatbreaking ground, 20 forest lands, 30 meadows, 40 waters 50 construction lands.Refer to shown in table 2, it is the B value under different land use type of the present invention and different vegetation cover degree.

B value under table 2 different land use type and different vegetation cover degree

Step 24, computational engineering measure factor E, utilizes statistics to calculate the engineering measure factor E value in each region; Refer to shown in formula (6),

$E=(1-\frac{S_t}{S}\&times;\mathrm{\&alpha;})(1-\frac{S_d}{S}\&times;\mathrm{\&beta;})(1-\frac{\mathrm{\&lambda;}\&times;N_{d1}+\mathrm{\&epsiv;}\&times;N_{d2}}{A\&times;S})---\left(6\right)$

In formula, S _tfor terraced fields area, S _dfor silt arrester control area, S is the soil total area, α, and β is respectively the husky coefficient of subtracting of terraced fields and silt arrester, is respectively 0.763 and 1; N _d1, N _d2the quantity that is respectively check dam, check dam, unit is seat, and λ and ε are respectively the sediment trapping quota of check dam and check dam, are respectively 1000t/ seat and 100t/ seat, and A is zone leveling soil erosion modulus, and unit is t/km ².

Silt arrester control area calculates according to dissimilar silt arrester control area standard in the < < of upper and middle reaches, the Yellow River management board silt arrester design > >, its middle-size and small-size silt arrester control area <1km ², medium-sized silt arrester control area is 1-3km ², large-scale silt arrester control area is 3-8km ², the present invention gets the intermediate value within the scope of each type silt arrester control area, and small-sized, medium-sized and large-scale silt arrester control area is respectively 0.5,2 and 5.5km ².

Correction about E value: be negative region for result of calculation, we think that he has reached good improvement degree, has artificially defined its E=0.2.

Step 25, calculates tillage control measure factor T: tillage control measure factor values T reduces soil losses according to Contour farming under different gradient condition to determine.

Tillage control measure factor values under table 3 different gradient

In formula: S is slope factor; θ is the gradient; L is slope length factor; λ is length of grade (m); M is length of grade index, according to gradient difference, gets different values.

Step 26, calculates ditch eclipse factor G; Consider that CSLE model only calculates sheetflood, introduce the effect that ditch eclipse factor can be considered shallow gully erosion to a certain extent.To have or not rainfall data, use two kinds of algorithms below.

(1) there is the calculating of rainfall amount data: the relation of ditch eclipse factor and ground inclination and characteristics of rainfall, shown in (9),

G=1+β (9)

In formula, G is shallow ridges coefficient, dimensionless; β is correction coefficient, affected by time rainfall, water catchment area and the gradient.When ground inclination is less than shallow gully erosion generation critical grade, during without shallow gully erosion, β=0; When shallow ridges generation critical grade is 15 °, β calculating formula is following formula (10),

$\mathrm{\&beta;}=\left(\frac{\mathrm{\&theta;}-15}{15}\right)[1.003{\left({\mathrm{RI}}_{30}\right)}^{0.103}-1]---\left(10\right)$

In formula: θ be ground inclination (°); R is time rainfall amount (mm); I ₃₀for inferior rainfall maximum 30-min rainfall intensity (mm/min).

(2) without the calculating of rainfall data: when without rainfall data in detail, available following formula (11) estimation annual G value:

G=1+1.60sin（θ-15) (11)

In above-mentioned steps 4, erosion assessment process is as described below,

Potential erosion 1 (A01): A=RK, represents the erosion of weather and soil;

Potential erosion 1 (A02): A=RKLS, does not consider water guarantor measure, only considers the erosion of weather, soil and landform;

Basin soil losses 1 (A1): A=RKLSBET; Do not consider gullying;

Basin soil losses 2 (A2): A=RKLSBETg; Consider gullying;

Water erosion grade (A31): according to Ministry of Water Resources's Criterion of Classification of Soil Erosion (SL190-2007), complete the evaluation to intensity of water erosion grade.Refer to shown in table 4, it is water erosion (face is invaded) evaluation criterion; Needed input parameter comprises: the gradient, vegetation cover degree, soil utilization.

, there is the situation of deposition in non-corrosive referring to compared with river or the low-lying position of trench bottom (do not distinguish and plough and woods meadow); Mired corrodes and refers to that the gradient is less than or equal to the woods meadow of 5 degree; The realization of upper table, need to call 3 figure layers, i.e. land-use map, vegetation coverage map and slope maps.Refer to shown in table 5, it is soil erosion intensity and gradation standard.

Table 4 water erosion (face is invaded) evaluation criterion

Table 5 soil erosion intensity and gradation standard

Wind erosion class (A32): according to Ministry of Water Resources's Criterion of Classification of Soil Erosion (SL190-2007), complete the evaluation to intensity of wind erosion grade.

Due to the more difficult acquisition of wind erosion thickness data, generally according to vegetation coverage and surface configuration, betray not.Needed input parameter comprises: vegetation cover degree, pedological map.Refer to shown in table 6 and 7, it is respectively wind erosion evaluation criterion and soil erosion intensity evaluation criterion.

Table 6 wind erosion evaluation criterion

Table 7 soil erosion intensity grade evaluation standard

In above-mentioned steps 4, the analytic process of the overall shape of the watershed soil erosion is as described below:

Average annual erosion intensity (At): refer to calculate the average soil erosion intensity in basin or region, concrete computation process is as following formula (12),

$A_t=\left(\underset{c=1}{\overset{n}{\mathrm{\&Sigma;}}}(A_c\&times;s_c)\right)/\underset{c=1}{\overset{n}{\mathrm{\&Sigma;}}}{s}_c---\left(12\right)$

In formula: A _cyear soil erosion modulus (t/km for computing unit ²a), A _tfor calculating the annual soil erosion modulus (t/km in basin ²a), S _cfor computing unit area.

Produce husky total amount (Et) per year: refer to calculate the year Erosion and Sediment Production total amount in basin or region, concrete computation process is as following formula (13),

$E_t=A_t\&times;S_t=\underset{c=1}{\overset{n}{\mathrm{\&Sigma;}}}{A}_c\&times;s_c---\left(13\right)$

In formula: E _tfor calculating the year sediment production (t/a) in basin, all the other same formula (12).

Erosion grade index (Ie): for a parameter of basin or region erosion intensity grade general status.Concrete computation process is as following formula (14)-(16),

I=I ₁×a+I ₂×b (14)

$I_1=\underset{10}{\overset{16}{\mathrm{\&Sigma;}}}s\&times;c_{\mathrm{wat}}\&times;{10}^{2\&times;(C_{\mathrm{wat}}-10)}---\left(15\right)$

$I_2=\underset{20}{\overset{26}{\mathrm{\&Sigma;}}}s\&times;c_{\mathrm{wib}}\&times;{10}^{2(C_{\mathrm{win}}-20)}---\left(16\right)$

In formula: I is comprehensive erosion intensity index; I ₁for intensity of water erosion index; I ₂for intensity of wind erosion index; The b that lies prone is respectively water erosion and the shared drainage area ratio of wind erosion in basin; S is each erosion intensity grade area and drainage area ratio in basin; C is erosion intensity grade (water erosion is corroded violent erosion from mired and is respectively 10,11,12,13,14,15,16, and wind erosion is corroded violent erosion from mired and is respectively 20,21,22,23,24,25,26).I value has considered water erosion and intensity of wind erosion in basin, can be used to water erosion between comparison basin, the comprehensive erosion intensity of eroding size.

Simple eigenvalue (St) is carried out the most basic statistical indicator calculating and is comprised min, max, mean, std, middle number, mode.

Frequency curve (Fi): to zoning or each strength grade frequency statistics of basin; Result is as shown in table 8 below, and intensity of water erosion code is above 10,11 ... 16; Intensity of wind erosion code is 20,21 ..., 26.

Table 8 erosion intensity grade statistical form

In above-mentioned steps 4, the space statistical analysis of space statistical analysis model comprises: conventional statistical study, spatial autocorrelation analysis, regretional analysis, trend analysis.The isoparametric statistical study of average, summation, variance, frequency, coefficient of kurtosis of the main complete paired data set of conventional statistical study.Spatial autocorrelation analysis is understanding spatial distribution characteristic, selects suitable space scale to complete the most frequently used method of spatial analysis.Regretional analysis is for analyzing the correlationship between two or more sets variablees, and common regression analysis equation has: linear regression, index return, logarithm regression, multiple regression etc.By mathematical model, simulate space distribution and the time course of geographic entity, the insufficient section interpolation between the measured data point of geographic element spatial and temporal distributions or prediction out.

Space operational method utilizes the IPixelBlockCursor that AO provides to obtain the value that different figure layers are put, and by the space traversal to a data, when these point value match operations require, difference figure layer point is carried out to computing, realizes like this space statistical study.

The demonstration of space diagram data, first to carry out the division of figure layer to figure, then the pixel in each figure layer is carried out to pel numbering, keep the attribute of each pel corresponding one by one, be that pel numbering is unique, be Property ID, generally comprise the attribute number of line, 3 kinds of the attribute numbers of the attribute number in district and point.In order to articulate like clockwork corresponding attribute, in this kind of line style, do a few thing.Can not change line style, and only change line color according to attribute difference, as can be still black between stratum and stratum, intrusive contact can be red or other colors that you think and other line parameter can distinguish.

The editor of district's (piece) pel is first simple compared with line chart.Formation zone, first will form the curve (segmental arc) in district, and these curves are obtained through " line turns segmental arc " by line chart unit.In the online pel editor of these segmental arcs, passed through " adjustment by method of junction point ", " automatic shearing broken string " etc. and processed early stage, directly line chart unit has been carried out to " line turns segmental arc " operation, then carried out " topology rebuilding " and generate all face pels.Then be that each district is changed to color, and match the look number of standard.

The information attribute of some pel does not generally have identical, with the point of attribute, is generally subgraph, and No. ID of pel, attribute list are corresponding one by one, generate a some pel, there is no space between region-position code and region-position code, and reading out data point is loaded as point diagram layer, then shows.

Raster data is the Method of Data Organization that comes representation space atural object or phenomenon to distribute with the form of two-dimensional matrix.Each matrix unit is called a grid cell (cell).The loading of raster data, demonstration, read grid point data exactly, analyzes the process that reads to internal memory and show after its space coordinates, coordinate, lattice dimensions, ranks number, Wu Zhi district numerical value.

The present invention takes full advantage of GIS spatial analysis functions, the result of its output (soil erosion intensity map) can provide the erosion intensity at Nei Ge position, basin, thereby both can meet the requirement of raceway groove Water Conservancy engineering design, also can meet water-and-soil conservation measures and lay the demand with planning; Meanwhile, the model that this method is set up has absorbed the achievement, the particularly model such as USLE, CSLE of slope soil erosion model preferably, than being easier to, promotes.As for the division of unit, can be irregular polygon (being called plot), can be also regular grid.

The foregoing is only preferred embodiment of the present invention, is only illustrative for invention, and nonrestrictive.Those skilled in the art is understood, and in the spirit and scope that limit, can carry out many changes to it in invention claim, revise, and even equivalence, but all will fall within the scope of protection of the present invention.

Claims (10)

1. loess plateau mesoscale basin year Model of Erosion and Sediment Production computing method, is characterized in that, its detailed process is;

Step 1, is divided into some unit basin, each unit is approximate represent one domatic;

Step 2, utilizes slope soil erosion test observation achievement in research and understanding, based on slope soil erosion test observation data, sets up relation between each factor and erosion amount, i.e. model of element;

Step 3, completes the calculating to each unit erosion amount with described model of element;

Step 4, utilizes spatial analysis and spatial statistics method, by polymerization computing, completes the evaluate and analyze of watershed soil erosion overall state;

In above-mentioned steps 2, the computation process of obtaining each factor is as follows:

Step 21, calculates website climatic factor R1;

Step 22, calculates climatic factor and calculates interpolation R2;

Step 23, calculates plant measures factor B, with land-use map and vegetation cover degree figure, calculates certain region plant measures factor B;

Step 24, computational engineering measure factor E, utilizes statistics to calculate the engineering measure factor E value in each region;

Step 25, calculates tillage control measure factor T, and tillage control measure factor values T reduces soil losses according to Contour farming under different gradient condition and determines.

2. mesoscale basin, loess plateau according to claim 1 year Model of Erosion and Sediment Production computing method, it is characterized in that, in above-mentioned steps 21, calculate the rainfall erosivity that website climatic factor R1 comprises one or more websites, its rainfall amount according to respective site calculates;

Computing formula is as shown in following (1) and (2):

$F=\left(\mathrm{\&Sigma;}{P_i}^2\right)/P---\left(1\right)$

R=αF ^β (2)

In formula, P is average annual rainfall (mm); Pi is the average rainfall (mm) of the i month; R is average rainfall erosivity (MJmmhm for many years ^-2h ^-1a ^-1); α=0.3589, β=1.9462; F exponent m m size is relevant with the season distribution of mean annual rainfall P, and span is between P12-1～P.

3. mesoscale basin, loess plateau according to claim 1 and 2 year Model of Erosion and Sediment Production computing method, is characterized in that, the detailed process of calculating plant measures factor B in above-mentioned steps 23 is:

Step 231, obtains land-use map;

Step 232, calculates vegetation cover degree figure;

Step 233, calculates annual vegetation cover degree;

Step 234, calculates plant measures factor B.

4. mesoscale basin, loess plateau according to claim 3 year Model of Erosion and Sediment Production computing method, is characterized in that, in above-mentioned steps 232, utilize N DVI to calculate vegetation cover degree according to following formula (3),

$F=100\&times;\frac{\mathrm{float}\left(\mathrm{NDVI}\right)-\mathrm{NDV}{I}_{\min }}{{\mathrm{NDVI}}_{\max }-{\mathrm{NDVI}}_{\min }}---\left(3\right)$

In formula, F is vegetation cover degree, and NDVI is normalized differential vegetation index, and NDVImin and NDVImax are respectively without the N DVI value in vegetation area and the NDVI value of the good covering area of vegetation; Test basis is upper at the scene, in conjunction with the typically sampling of class of remote sensing images, determines the NDVI threshold value of the highest and minimum vegetation cover degree; The result of sampling is carried out to statistical with corresponding NDVI value, and when NDVI value > >=NDVImax, cover degree can be considered 100%, and when NDVI value < NDVImin, cover degree can be considered 0%.

5. mesoscale basin, loess plateau according to claim 4 year Model of Erosion and Sediment Production computing method, is characterized in that, in above-mentioned steps 232, the process that NDVI data are calculated is:

For less basin, utilize TM data to calculate N DVI, refer to formula (4),

$\mathrm{NDVI}=\frac{\mathrm{TM}4-\mathrm{TM}3}{\mathrm{TM}4+\mathrm{TM}3}---\left(4\right)$

In formula, NDVI is required, normalized differential vegetation index; TM3, TM4 are respectively near infrared and the infrared band of remote sensing image;

For compared with large watershed, directly obtain.

6. mesoscale basin, loess plateau according to claim 5 year Model of Erosion and Sediment Production computing method, it is characterized in that, in above-mentioned steps 233, while calculating annual vegetation cover degree, by each month flood season vegetation cover degree value and flood season each month rainfall erosivity surface be weighted on average, the average vegetation cover degree that can obtain each month in flood season, computing formula is as shown in the formula (5)

$F_r=\frac{\underset{m=5}{\overset{9}{\mathrm{\&Sigma;}}}{F}_m\&times;R_m}{\underset{m=5}{\overset{9}{\mathrm{\&Sigma;}}}{R}_m}---\left(5\right)$

In formula: Fm is for pressing the vegetation cover degree of each month in flood season of above-mentioned (3) formula calculating, and Rm is the value on the rainfall erosivity surface of each month in flood season; Fr is cover degree and the weighted mean value of rainfall erosivity in flood season.

7. mesoscale basin, loess plateau according to claim 1 year Model of Erosion and Sediment Production computing method, is characterized in that, in above-mentioned steps 24, utilize statistics to calculate the engineering measure factor E value in each region; Refer to shown in formula (6),

$E=(1-\frac{S_t}{S}\&times;\mathrm{\&alpha;})(1-\frac{S_d}{S}\&times;\mathrm{\&beta;})(1-\frac{\mathrm{\&lambda;}\&times;N_{d1}+\mathrm{\&epsiv;}\&times;N_{d2}}{A\&times;S})---\left(6\right)$

In formula, S _tfor terraced fields area, S _dfor silt arrester control area, S is the soil total area, α, and β is respectively the husky coefficient of subtracting of terraced fields and silt arrester, is respectively 0.763 and 1; N _d1, N _d2the quantity that is respectively check dam, check dam, unit is seat, and λ and ε are respectively the sediment trapping quota of check dam and check dam, are respectively 1000t/ seat and 100t/ seat, and A is zone leveling soil erosion modulus, and unit is t/km ².

8. mesoscale basin, loess plateau according to claim 2 year Model of Erosion and Sediment Production computing method, is characterized in that, also obtain ditch eclipse factor G in above-mentioned steps 2, to have or not rainfall data, comprise two kinds of algorithms,

(1) there is the calculating of rainfall amount data: the relation of ditch eclipse factor and ground inclination and characteristics of rainfall, shown in (9),

G=1+β (9)

In formula, G is shallow ridges coefficient, dimensionless; β is correction coefficient, affected by time rainfall, water catchment area and the gradient;

(2) without the calculating of rainfall data: when without rainfall data in detail, available following formula (11) estimation annual G value:

G=1+1.60sin(θ-15) (11)。

9. mesoscale basin, loess plateau according to claim 8 year Model of Erosion and Sediment Production computing method, is characterized in that,

Above-mentioned while obtaining ditch eclipse factor G, when ground inclination is less than shallow gully erosion generation critical grade, during without shallow gully erosion, β=0; When shallow ridges generation critical grade is 15 °, β calculating formula is following formula (10),

$\mathrm{\&beta;}=\left(\frac{\mathrm{\&theta;}-15}{15}\right)[1.003{\left({\mathrm{RI}}_{30}\right)}^{0.103}-1]---\left(10\right)$

In formula: θ be ground inclination (°); R is time rainfall amount (mm); I ₃₀for inferior rainfall maximum 30-min rainfall intensity (mm/min).

10. mesoscale basin, loess plateau according to claim 8 year Model of Erosion and Sediment Production computing method, is characterized in that, in above-mentioned steps 4, the analytic process of the overall shape of the watershed soil erosion is as described below:

Average annual erosion intensity At: refer to calculate the average soil erosion intensity in basin or region, concrete computation process is as following formula (12),

$A_t=\left(\underset{c=1}{\overset{n}{\mathrm{\&Sigma;}}}(A_c\&times;s_c)\right)/\underset{c=1}{\overset{n}{\mathrm{\&Sigma;}}}{s}_c---\left(12\right)$

In formula: A _cyear soil erosion modulus (t/km for computing unit ²a), A _tfor calculating the annual soil erosion modulus (t/km in basin ²a), S _cfor computing unit area;

Produce husky total amount Et per year: refer to calculate the year Erosion and Sediment Production total amount in basin or region, concrete computation process is as following formula (13),

$E_t=A_t\&times;S_t=\underset{c=1}{\overset{n}{\mathrm{\&Sigma;}}}{A}_c\&times;s_c---\left(13\right)$

In formula: E _tfor calculating the year sediment production (t/a) in basin, all the other same formula (12);

Erosion grade index Ie: for a parameter of basin or region erosion intensity grade general status; Concrete computation process is as following formula (14)-(16),

I=I ₁×a+I ₂×b (14)

$I_1=\underset{10}{\overset{16}{\mathrm{\&Sigma;}}}s\&times;c_{\mathrm{wat}}\&times;{10}^{2\&times;(C_{\mathrm{wat}}-10)}---\left(15\right)$

$I_2=\underset{20}{\overset{26}{\mathrm{\&Sigma;}}}s\&times;c_{\mathrm{wib}}\&times;{10}^{2(C_{\mathrm{win}}-20)}---\left(16\right)$

In formula: I is comprehensive erosion intensity index; I ₁for intensity of water erosion index; I ₂for intensity of wind erosion index; A, b are respectively water erosion and the shared drainage area ratio of wind erosion in basin; S is each erosion intensity grade area and drainage area ratio in basin; C is erosion intensity grade;

Simple eigenvalue St: carry out the most basic statistical indicator and calculate;

Frequency curve Fi: to zoning or each strength grade frequency statistics of basin.