CN113987778A - Water and soil loss analog value space-time weighting correction method based on field station - Google Patents

Water and soil loss analog value space-time weighting correction method based on field station Download PDF

Info

Publication number
CN113987778A
CN113987778A CN202111236015.4A CN202111236015A CN113987778A CN 113987778 A CN113987778 A CN 113987778A CN 202111236015 A CN202111236015 A CN 202111236015A CN 113987778 A CN113987778 A CN 113987778A
Authority
CN
China
Prior art keywords
soil
water
value
factor
data
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202111236015.4A
Other languages
Chinese (zh)
Inventor
李畅
姚佳吟
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Central China Normal University
Original Assignee
Central China Normal University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Central China Normal University filed Critical Central China Normal University
Priority to CN202111236015.4A priority Critical patent/CN113987778A/en
Publication of CN113987778A publication Critical patent/CN113987778A/en
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling

Abstract

The invention discloses a water and soil loss analog value space-time weighting correction method based on a field station, which takes a space-time geographic weighting regression model as an example for correcting a soil erosion modulus of a CSLE model and comprises the following steps: acquiring basic data of water and soil loss influence factors and calculating a regional dynamic water and soil loss simulation value by using a CSLE model; acquiring and processing the actual water and soil loss measurement value of the water and soil loss monitoring station in the area; thirdly, performing regression fitting on the simulated value and the measured value by using a space-time geographic weighting model, and finally obtaining a soil erosion modulus correction model by using a heuristic intelligent optimization algorithm; and fourthly, checking the precision of the correction model. The invention utilizes the actually measured time-space data of the soil and water conservation field station to correct the soil and water loss equation model and the soil and water loss single influence factor value and the soil erosion modulus analog value thereof, and can effectively improve the measuring, calculating and predicting accuracy of the soil and water loss.

Description

Water and soil loss analog value space-time weighting correction method based on field station
Technical Field
The invention belongs to the technical field of remote sensing and geographic process simulation, and relates to a water and soil loss analog value space-time weighting correction method based on a wild appearance site.
Background
Soil erosion has been a serious ecological problem in the world. The accuracy of the computer simulation value of the regional water and soil loss is improved, the rapid and efficient monitoring and control of the water and soil loss are facilitated, and the method has profound significance on Chinese ecological civilization construction, "Yangtze river protection" and social and economic sustainable development. At present, based on GIS and RS technologies, the calculation of water and soil loss by using a water and soil loss equation is a main mode for quantitative evaluation of water and soil loss in China.
General Soil Loss equations (USLE), modified Soil Loss equations (ruisle), and Chinese Soil Loss Equations (CSLE) are widely used in academic research, data statistics, and production and construction projects. Although CSLE is considered to be closer to the actual situation in china, the model was originally proposed mainly based on the loess plateau water and soil loss data summary, and the calculated soil erosion modulus is only a computer simulation value, and has a certain deviation from the actual water and soil loss in the field in other areas.
Disclosure of Invention
The invention aims to provide a method for solving the problem of deviation between a water loss and soil loss analog value and an actual value and correcting the soil erosion modulus of a water loss and soil loss equation by utilizing field station actual measurement data. The method is suitable for correcting the water and soil loss equation model and the water and soil loss single influence factor value and the soil erosion modulus thereof.
The technical scheme adopted by the invention is as follows: a water and soil loss analog value space-time weighting correction method based on field stations comprises the following steps:
step 1: acquiring basic data of the water and soil loss influence factors and calculating a regional dynamic water and soil loss simulation value by using a CSLE model;
the CSLE model is as follows:
A＝R×K×L×S×B×E×T；
in the formula: a represents the annual average soil erosion amount in the unit of t.hm-2·a-1(ii) a R represents rainfall erosion force, unit MJ.mm.hm-2·h-1·a-1(ii) a K is soil erodability factor with unit t.hm2·h·hm-2·MJ-1·mm-1(ii) a L, S are topographic factors, also called slope length and gradient factors, all of which are dimensionless; b is a vegetation cover and biological measure factor without dimension; e is an engineering measure factor and is dimensionless; t is a factor of cultivation measure and has no dimension；
Step 2: acquiring and processing the actual water and soil loss measurement value of the area through a water and soil loss monitoring station;
and step 3: performing regression fitting on the simulated value and the measured value by adopting a space-time weighting method to obtain a soil erosion modulus correction model;
and 4, step 4: and predicting the future water and soil loss by using the soil erosion modulus correction model.
Compared with the prior art, the invention has the beneficial effects that:
(1) at present, relevant research at home and abroad is mostly correction on a certain influence factor, regionality is strong, but universality is weak, and meanwhile, the research on correction of the soil erosion modulus is insufficient. The time-space weighting correction method for the water and soil loss analog value based on the field station can simultaneously realize the correction of the influence factors and the soil erosion modulus, is suitable for constructing a correction model under each water and soil loss equation frame, and has strong compatibility and universality.
(2) The Chinese water and soil loss equation and other Chinese regional water and soil loss equations are provided by performing regression fitting according to observation data of a corresponding regional water and soil conservation observation station on the basis of the US USLE model and the RUSLE model, the time and related researches are provided more than 10 years ago, and the corrected technical method needs to be further updated, optimized and expanded. The water and soil loss analog value space-time weighting correction method based on the field station fully considers the space-time diversity rule, introduces a GTWR model for the first time to perform local regression fitting of a computer analog value and a field station actual observed value, ensures the strong regionality of the correction method, optimizes GTWR by adopting an intelligent search algorithm, and further improves the measurement and calculation accuracy of water and soil loss.
(3) Compared with other domestic regional water and soil loss equations, the method has the advantages that the time-space characteristics are considered, and the regional water and soil loss can be predicted by a time-space weighting correction method based on the water and soil loss analog value of the field site; the method is suitable for correcting the soil erosion modulus of each area; can be quickly and efficiently butted with the prior art of the water conservancy department for use. The invention utilizes the measured data of the soil and water conservation field station to correct the soil and water loss equation model and the soil and water loss single influence factor value and the soil erosion modulus analog value thereof, thereby effectively improving the measuring, calculating and predicting accuracy of the soil and water loss.
Drawings
FIG. 1 is a flow chart of an embodiment of the present invention;
FIG. 2 is a flowchart of a spatio-temporal geoweighted regression model (GTWR) construction correction model according to an embodiment of the present invention.
Detailed Description
In order to facilitate the understanding and implementation of the present invention for those of ordinary skill in the art, the present invention is further described in detail with reference to the accompanying drawings and examples, it is to be understood that the embodiments described herein are merely illustrative and explanatory of the present invention and are not restrictive thereof.
Referring to fig. 1, the invention provides a water and soil loss analog value space-time weighting correction method based on a field site, which comprises the following steps:
step 1: acquiring basic data of the water and soil loss influence factors and calculating a regional dynamic water and soil loss simulation value by using a CSLE model;
the CSLE model is:
A＝R×K×L×S×B×E×T；
in the formula: a represents the annual average soil erosion amount in the unit of t.hm-2·a-1(ii) a R represents rainfall erosion force, unit MJ.mm.hm-2·h-1·a-1(ii) a K is soil erodability factor with unit t.hm2·h·hm-2·MJ-1·mm-1(ii) a L, S are topographic factors, also called slope length and gradient factors, all of which are dimensionless; b is a vegetation cover and biological measure factor without dimension; e is an engineering measure factor and is dimensionless; t is a cultivation measure factor and is dimensionless;
the basic data of the soil erosion influence factors of the embodiment include:
(1) satellite photos or aerial photos;
(2) precipitation data: rainfall day by day;
(3) soil data: soil type, soil physicochemical property data;
(4) basic geographic data: digital line Drawing (DLG), Digital Elevation Model (DEM), topography;
(5) land utilization data: national land utilization annual change survey data and forest land grassland resource survey data
Water and soil conservation key engineering data: project type, area of implementation, distribution, number or area of major water and soil conservation measures;
(6) artificial water and soil loss data: and (5) maintaining monitoring information of water and soil of the production construction project.
In this embodiment, the remote sensing image data is preprocessed, and all images are preprocessed by radiation correction, orthorectification, fusion and mosaic.
The implementation obtains R, K, L, S, B, E, T seven impact factor values based on basic data and relevant national industry standards;
the R factor of this example is calculated as follows:
Figure BDA0003317690820000031
Figure BDA0003317690820000032
Figure BDA0003317690820000041
in the formula (I), the compound is shown in the specification,
Figure BDA0003317690820000042
is the average annual rainfall erosive power of many years, MJ.mm.hm-2·h-1·a-1(ii) a k takes 1,2, a..., 24, which means that a year is divided into 24 semimonths;
Figure BDA0003317690820000043
is the rainfall erosive power of the kth half-month, MJ.mm-2.h-1(ii) a Get i1, 2.... cndot.n; n refers to the time sequence of 1986-2015, and the subsequent updating is carried out according to the five-year sequence; j takes 0,1, a. m is the number of erosive rainfall days (the daily rainfall is more than or equal to 10mm) in the kth half and a month in the ith year; pi，j，kIn the ith year, the jth erosive rainfall in the kth half-and-half, mm; if there is no erosive rainfall in a certain half-month of a year, i.e. j is 0, let Pi，0，k0; alpha is a parameter, the alpha is 0.3937 in the warm season (5-9 months), and 0.3101 in the cold season (10-12 months, 1-4 months);
Figure BDA0003317690820000044
average rainfall erosive power for the kth half-month
Figure BDA0003317690820000045
The average annual rainfall erosive power of many years
Figure BDA0003317690820000046
The ratio of (A) to (B);
the K factor obtaining method of this embodiment is as follows:
updating and calculating the soil erodibility factor based on the collected runoff plot observation data and the first national water conservancy general survey water and soil conservation condition general survey soil erodibility factor calculation method; or directly adopting the first water conservancy general survey water and soil conservation condition general survey soil erodibility factor achievement or updating based on the observation data of the standard runoff plot, wherein the standard runoff plot calculates the soil erodibility factor K formula as follows:
Figure BDA0003317690820000047
in the formula, A represents the slope length of 22.13m, the slope of 9 percent (5 degrees), and the soil erosion modulus (t.hm) of the multi-year average (generally requiring more than 12 years of continuous observation and properly reducing the southern observation period) observed in the clear-ploughing leisure runoff plot-2·a-1) (ii) a R represents the average annual rainfall erosion force over many years (MJ mm. hm) corresponding to the plot soil erosion observation-2.h-1.a-1) (ii) a Through resampling, 10m spatial resolution (corresponding topographic map) is generatedScale 1: 10000) or 30m spatial resolution (corresponding to a terrain map scale 1: 50000);
the present embodiment calculates the L factor and the S factor by using DEM data, and the L, S factor is calculated as follows:
Figure BDA0003317690820000048
Figure BDA0003317690820000049
Figure BDA00033176908200000410
wherein λ isi，λi-1Represents the slope length of the ith and (i-1) th slope segments in m; m is a slope length index which changes along with the slope, and theta is the slope and is in unit degree; the resolution of the generated L, S raster data is 10m (corresponding to a terrain map scale of 1: 10000) or 30m (corresponding to a terrain map scale of 1: 50000);
the B factor obtaining method of the embodiment is as follows:
utilizing an MODIS normalized vegetation index NDVI product and a TM multispectral image (comprising 4 blue, green, red and near infrared wave bands), obtaining vegetation coverage FVC with 24 half-moon spatial resolution and 30m spatial resolution by adopting a fusion calculation method or a parameter revision method, and calculating a B factor by combining with 24 half-moon rainfall erosion force factor ratios; b factor grid data with 10m spatial resolution (corresponding to a topographic map scale of 1: 10000) or 30m spatial resolution (corresponding to a topographic map scale of 1: 50000) is generated through resampling;
the cultivated land B factor of the embodiment is obtained by referring to an industry experience value (table 1) provided by a water conservation and soil conservation monitoring center of the China Water conservancy department;
TABLE 1
Figure BDA0003317690820000051
Figure BDA0003317690820000061
Garden, forest, grass B are calculated as follows:
Figure BDA0003317690820000062
in the formula, WRiWhen the R factor is calculated, the erosion force of the ith half-moon rainfall accounts for the erosion force proportion of the whole year, and the value range is 0-1; SLRiThe soil loss proportion of the ith semilunar field, the forest land and the grassland is dimensionless, and the value range is 0-1;
SLR for tea garden and shrub forest landi：
Figure BDA0003317690820000063
SLR for orchard, woodland and other woodlandi：
SLRi＝0.44468×e(-3.20096×GD)-0.04099×e(FVC-FVC×GD)+0.025；
Grassland SLRi：
Figure BDA0003317690820000064
In the formula, the FVC is vegetation coverage calculated based on NDVI, and the value range is 0-1; the under-forest coverage of the GD arbor forest is in a value range of 0-1, comprises the under-forest coverage formed by all vegetation except the crown layer of the arbor forest, and is subjected to on-site investigation or empirical value;
the FVC calculation method is as follows:
Figure BDA0003317690820000065
Figure BDA0003317690820000066
wherein FVC is vegetation coverage, NDVI is the pixel NDVI value, NDVImax、NDVIminThe conversion coefficient of the ground class where the pixel is located; k is a nonlinear coefficient, the pixels where the NDVI maximum values and the bare soil NDVI minimum values of different vegetation are located are determined in the same climate type, and the average value of the NDVI in the pixels is taken as a conversion coefficient; NIR is the reflectance in the near infrared band; r is the reflectivity of the visible red wave band;
the embodiment carries out indoor interpretation and field rechecking on the remote sensing image, and extracts an E factor and a T factor by referring to an empirical value;
(1) the remote sensing image map is interpreted indoors to obtain a land utilization current situation map (vector data), and the value is assigned according to the classification of the land utilization current situation (GB/T210102017);
(2) rechecking and correcting the interpretation result by utilizing national land utilization annual change survey data, forest land grassland resource survey data, water and soil conservation key engineering data and artificial water and soil loss data;
(3) and assigning values to the interpretation graph according to the national rotation district name and code of the Classification of Water and soil conservation measures to obtain corresponding E factors and T factors, and finally converting the E factors and T factors into factor grid data.
In this embodiment, the product operation is performed on the seven influence factors to obtain a soil erosion modulus, namely a water and soil loss simulation value a under the CSLE model;
step 2: and acquiring the actual water loss and soil erosion measured values of the long-time sequence in the region through a water loss and soil erosion monitoring station, and sorting and counting the actual water loss and soil erosion measured values into data of the same time unit.
And step 3: and performing regression fitting on the simulated value and the measured value by using a space-time weighting method (such as a space-time geographic weighting regression model, a geographic space-time weighting neural network weighting regression model and the like) to obtain a soil erosion modulus correction model.
The invention takes a space-time geographic weighted regression model (GTWR) as an example to simulate the space-time change of water and soil loss.
The GTWR model is as follows:
yi＝β0(ui，vi，ti)+∑kβk(ui，vi，ti)xik+εi，i＝1，2，...；
wherein i represents the ith point in the n groups of observed values; (u)i，vi，ti) Representing the spatio-temporal coordinates (i ═ 1, 2.., n) and time at the regression point i; yi represents the dependent variable at the regression point i, xikRepresents the kth independent variable at the regression point i; beta is ak(ui，vi，ti) Representing the regression coefficient at the regression point i; epsiloniIs the error term for sample point i, obeys ∈i～N(0，σ2) And Cov (. epsilon.)i，εj)＝0(i≠j)。
In this embodiment, applying GTWR to the water and soil loss simulation value based on the wild appearance site, the soil erosion modulus calibration model is constructed as follows:
WASLi＝β0(ui，vi，ti)+∑kβk(ui，vi，ti)Aik+εi，i＝1，2，…n；
in the formula, AikIs an independent variable representing a computer-simulated value of a modulus of erosion due to soil erosion under k CSLE models of the ith observation point, namely a in the formula of step 1 in claim 1; i represents the ith observation point in the soil erosion monitoring station; WASLi(Water and soil loss) is a dependent variable, which represents the actual measured value of the soil and Water loss at the point i and is a dependent variable; (u)i，vi，ti) Representing geographic coordinates (latitude and longitude) and time coordinates; beta is ak(ui，vi，ti) Regression coefficients representing the kth independent variable of a sample point i; beta is a0(ui，vi，ti) Is a constant term of the equation; epsiloniIs a random error term obeying epsiloni～N(0，σ2) And Cov (. epsilon.)i，εj)＝0(i≠j)。
The above model algorithm is simplified as follows:
WASLi＝Aikβ+ε；
Figure BDA0003317690820000081
in this embodiment, the geographic coordinates and time coordinates of i observation points and the corresponding soil erosion modulus A thereofiActual measurement value WASL of water loss and soil erosioniInputting the GTWR into a model constructed by the GTWR;
in this embodiment, the ideal value of the regression coefficient
Figure BDA0003317690820000089
The calculation is as follows:
Figure BDA0003317690820000082
in the formula (I), the compound is shown in the specification,
Figure BDA0003317690820000083
the ideal value of the regression parameter of the ith observation point is represented; a. theTDenotes the independent variable AikThe transposed matrix of (2); w (u)i，vi，ti)＝diag(wi1，wi2，…，win) Is (u)i，vi，ti) A spatiotemporal weight matrix of (c);
measured value WASL of ith pointiIdeal value of
Figure BDA0003317690820000084
Comprises the following steps:
Figure BDA0003317690820000085
in the formula, SWASLThe hat matrix representing the GTWR model.
In this embodiment, the spatio-temporal weight matrix may be trained using a kernel computation or a geospatial weighted neural network (GNNWR).
The present invention takes kernel function calculation as an example.
(1) Gaussian kernel functions (Gaussian kernel functions):
spatio-temporal weight matrix W (u)i，vi，ti) Comprises the following steps:
Figure BDA0003317690820000086
Figure BDA0003317690820000087
Figure BDA0003317690820000088
Figure BDA0003317690820000091
Figure BDA0003317690820000092
Figure BDA0003317690820000093
Figure BDA0003317690820000094
Figure BDA0003317690820000095
in the formula, WijRepresenting spatio-temporal weight matrix elements;
Figure BDA0003317690820000096
representing the distance from the observation point i to the regression point j; τ represents a space-time ratio parameter; b represents bandwidth, which is used to describe spatial weight andparameters of the non-negative attenuation relation of spatial distance: the larger the bandwidth value is, the slower the weight decays with the increase of the space-time distance; the smaller the bandwidth value, the faster the weight decays with increasing spatiotemporal distance.
(2) Near Gaussian kernel function method (Bi-square kernel function): when the data are very discrete, a near-Gaussian kernel function model is selected, points with no influence or little influence are removed, the calculation efficiency is improved, and a large amount of calculation caused by the long tail effect of the Gaussian kernel function is avoided.
Spatio-temporal weight matrix W (u)i，vi，ti) Comprises the following steps:
Figure BDA0003317690820000097
in the formula, WijRepresenting spatio-temporal weight matrix elements;
Figure BDA0003317690820000098
representing the distance from the observation point i to the regression point j; b represents the bandwidth, the regression point j is in the range of the bandwidth, the weight of the data point is calculated, and the weight of the part exceeding the bandwidth is totally marked as 0.
The selection of the bandwidth value of the kernel function has an important influence on the fitting result of the correction model constructed by the GTWR, and if the bandwidth value is too small, the correction result is over-fitted; if the bandwidth value is too large, the accuracy of the correction result is reduced. Currently, cross validation and akachi pool information criteria are the mainstream bandwidth selection methods.
(1) Cross validation (Cross validation, CV)
Firstly, the space-time distance range of the data is normalized, the value range of the bandwidth b is given, and then the optimal bandwidth is obtained by traversing and searching in the value range. The algorithm is as follows:
Figure BDA0003317690820000099
in the formula (I), the compound is shown in the specification,
Figure BDA0003317690820000101
and (3) representing an ideal value of water and soil loss in the bandwidth of the ith observation point, and when the regression parameters are estimated, in order to prevent the overfitting phenomenon, not including the i observation point, performing regression parameter calculation only according to the data around the i observation point. And drawing the different bandwidths and the different CVs into a trend line, and finding out the corresponding optimal bandwidth when the CV value is minimum.
(2) Chichi information criterion (AIC)
Figure BDA0003317690820000102
Figure BDA0003317690820000103
In the formula, n represents the number of sample points;
Figure BDA0003317690820000104
representing a GTWR model maximum likelihood variance estimate; tr (S)WASL) Representing hat matrix SWASLThe trace of (2); RSS represents the sum of the squares of the residuals.
The embodiment introduces an intelligent optimization algorithm on the basis of the GTWR model, for example: and a genetic algorithm or a particle swarm algorithm, wherein after the initial value of the coefficient of the correction model is obtained through the GTWR model, a space threshold value is set, and random space search is carried out on the initial value of the coefficient to obtain an optimal solution.
Taking particle swarm optimization as an example:
Figure BDA0003317690820000105
Figure BDA0003317690820000106
Figure BDA0003317690820000107
Figure BDA0003317690820000108
Figure BDA0003317690820000109
Figure BDA00033176908200001010
wherein, i is 1,2, …, m, which represents the ith group of particles, i.e. the ith group of coefficient initial values each includes n regression coefficients β, a bandwidth parameter b and a space-time ratio parameter τ in this embodiment; d is 1,2, …, D represents the D-th dimension, and n +2 coefficients are optimized in this embodiment, so D is n + 2;
Figure BDA00033176908200001011
a position representing an initial value of the ith set of coefficients;
Figure BDA00033176908200001012
indicating the location where the ith set of coefficient initial values has experienced the best;
Figure BDA00033176908200001013
represents the best position of the whole particle group;
Figure BDA00033176908200001014
a speed representing an initial value of the ith set of coefficients; w is a non-negative number, called the inertia factor; r is1And r2Is [0,1 ]]Random numbers transformed within a range; alpha is a constraint factor and controls the weight of the speed; c. C1And c2Is a learning factor, i.e., an acceleration constant.
And 4, step 4: and predicting the future water and soil loss by using the soil erosion modulus correction model.
The embodiment also provides a verifying partyThe method is used for verifying a time-space weighting correction method based on water and soil loss analog values of field stations and constructing a soil erosion modulus correction model of a certain area for m years; the soil erosion modulus A under the CSLE model of the (m + 1) th yearmInputting a correction equation to obtain a corrected water and soil loss ideal value of the (m + 1) th year
Figure BDA0003317690820000111
And comparing with the actual measured value WASL of the soil erosion of the place in the m +1 yearm+1Carrying out comparison and inspection; if it is
Figure BDA0003317690820000112
And WASLm+1And (4) approximating, predicting the future water and soil loss amount of the ground by using the soil erosion modulus correction model.
The method combines field observation and indoor computer simulation, improves the simulation accuracy of the water and soil loss equation, is suitable for correcting the current mainstream water and soil loss equation, and can correct the soil erosion modulus and the water and soil loss influence factor. The water and soil conservation principle and method based on the spatial information technology are enriched, and the practicability of the water and soil conservation principle and method is effectively promoted.
It should be understood that the above description of the preferred embodiments is given for clarity and not for any purpose of limitation, and that various changes, substitutions and alterations can be made herein without departing from the spirit and scope of the invention as defined by the appended claims.

Claims (10)

1. A water and soil loss analog value space-time weighting correction method based on field stations is characterized by comprising the following steps:
step 1: acquiring basic data of the water and soil loss influence factors and calculating a regional dynamic water and soil loss simulation value by using a CSLE model;
the CSLE model is as follows:
A＝R×K×L×S×B×E×T；
in the formula: a represents the annual average soil erosion amount in the unit of t.hm-2·a-1(ii) a R represents rainfall erosion force, unit MJ.mm.hm-2·h-1·a-1(ii) a K is soil erodability factor with unit t.hm2·h·hm-2·MJ-1·mm-1(ii) a L, S are topographic factors, also called slope length and gradient factors, all of which are dimensionless; b is a vegetation cover and biological measure factor without dimension; e is an engineering measure factor and is dimensionless; t is a cultivation measure factor and is dimensionless;
step 2: acquiring and processing the actual water and soil loss measurement value of the area through a water and soil loss monitoring station;
and step 3: performing regression fitting on the simulated value and the measured value by adopting a space-time weighting method to obtain a soil erosion modulus correction model;
and 4, step 4: and predicting the future water and soil loss by using the soil erosion modulus correction model.
2. The method for correcting the water and soil loss analog value space-time weighting based on the field station according to claim 1, is characterized in that: in the step 1, acquiring basic data of water and soil loss influence factors, wherein the basic data comprises remote sensing image data, specifically comprising satellite photos and aerial photos; precipitation data, specifically including daily rainfall; soil data, specifically including soil type and soil physicochemical property data; basic geographic data, specifically including a digital line drawing, a digital elevation model and a topographic map; the land utilization data specifically comprises national land utilization annual change survey data and forest land grassland resource survey data; water and soil conservation key engineering data specifically comprising engineering types, implementation areas, and distribution, quantity or area of water and soil conservation measures; the artificial water and soil loss data specifically comprises water and soil conservation monitoring information of a production construction project.
3. The field-site-based water and soil loss analog value space-time weighting correction method according to claim 2, characterized in that: the remote sensing image data is preprocessed, and all images are preprocessed through radiation correction, orthorectification, fusion and mosaic.
4. The method for correcting the water and soil loss analog value space-time weighting based on the field station according to claim 1, is characterized in that: in step 1, acquiring R, K, L, S, B, E, T seven influence factor values based on basic data and relevant national industry standards;
the R factor is calculated as follows:
Figure FDA0003317690810000021
Figure FDA0003317690810000022
Figure FDA0003317690810000023
in the formula (I), the compound is shown in the specification,
Figure FDA0003317690810000024
is the average annual rainfall erosive power of many years, MJ.mm.hm-2·h-1·a-1(ii) a k takes 1,2, … …,24, which means that a year is divided into 24 semimonths;
Figure FDA0003317690810000025
the rainfall erosive power of the kth half-moon, MJ.mm.hm-2·h-1(ii) a Taking 1,2, … …, N; n refers to the time sequence of 1986-2015, and the subsequent updating is carried out according to the five-year sequence; j is 0,1, … …, m; m is the number of erosive rainfall days (the daily rainfall is more than or equal to 10mm) in the kth half and a month in the ith year; pi,j,kIn the ith year, the jth erosive rainfall in the kth half-and-half, mm; if there is no erosive rainfall in a certain half-month of a year, i.e. j is 0, let Pi,0,k0; alpha is a parameter, the alpha is 0.3937 in the warm season, and 0.3101 in the cold season;
Figure FDA0003317690810000026
average rainfall erosive power for the kth half-month
Figure FDA0003317690810000027
The average annual rainfall erosive power of many years
Figure FDA0003317690810000028
The ratio of (A) to (B);
the K factor acquisition method comprises the following steps:
updating and calculating the soil erodibility factor based on the collected runoff plot observation data and the first national water conservancy general survey water and soil conservation condition general survey soil erodibility factor calculation method; or directly adopting the first water conservancy general survey water and soil conservation condition general survey soil erodibility factor achievement or updating based on the observation data of the standard runoff plot, wherein the standard runoff plot calculates the soil erodibility factor K formula as follows:
Figure FDA0003317690810000029
wherein A represents the slope length of 22.13m and the slope of 9 percent, and the average soil erosion modulus (t.hm) of the clear-ploughing leisure runoff plot observed for many years-2·a-1) (ii) a R represents the average annual rainfall erosion force (MJ mm hm) for many years corresponding to the observation of soil erosion in the community-2·h-1·a-1) (ii) a Generating K-factor grid data with 10m spatial resolution or 30m spatial resolution through resampling;
and calculating an L factor and an S factor by using DEM data, wherein the L, S factor is calculated as follows:
Figure FDA00033176908100000210
Figure FDA00033176908100000211
Figure FDA0003317690810000031
wherein λ isi，λi-1Represents the slope length of the ith and (i-1) th slope segments in m; m is a slope length index which changes along with the slope, and theta is the slope and is in unit degree; the generated L, S raster data has the resolution of 10m or 30 m;
the B factor acquisition method comprises the following steps:
utilizing an MODIS normalized vegetation index NDVI product and a TM multispectral image, obtaining vegetation coverage FVC with spatial resolution of 24 half-and-a-half and 30m by adopting a fusion calculation method or a parameter revision method, and calculating a B factor by combining with a rainfall erosive power factor ratio of 24 half-and-a-half; generating B factor grid data with 10m spatial resolution or 30m spatial resolution through resampling;
the cultivated land B factor is obtained by referring to an industry experience value provided by a water conservation and soil conservation monitoring center of the China Water conservancy department;
garden, forest, grass B are calculated as follows:
Figure FDA0003317690810000032
in the formula, WRiWhen the R factor is calculated, the erosion force of the ith half-moon rainfall accounts for the erosion force proportion of the whole year, and the value range is 0-1; SLRiThe soil loss proportion of the ith semilunar field, the forest land and the grassland is dimensionless, and the value range is 0-1;
SLR for tea garden and shrub forest landi：
Figure FDA0003317690810000033
SLR for orchard, woodland and other woodlandi：SLRi＝0.44468×e(-3.20096×GD)-0.04099×e(FVC-FVC×GD)+0.025；
Grassland SLRi：
Figure FDA0003317690810000034
In the formula, the FVC is vegetation coverage calculated based on NDVI, and the value range is 0-1; the under-forest coverage of the GD arbor forest is in a value range of 0-1, comprises the under-forest coverage formed by all vegetation except the crown layer of the arbor forest, and is subjected to on-site investigation or empirical value;
the FVC calculation method is as follows:
Figure FDA0003317690810000035
Figure FDA0003317690810000036
wherein FVC is vegetation coverage, NDVI is the pixel NDVI value, NDVImax、NDVIminThe conversion coefficient of the ground class where the pixel is located; k is a nonlinear coefficient, the pixels where the NDVI maximum values and the bare soil NDVI minimum values of different vegetation are located are determined in the same climate type, and the average value of the NDVI in the pixels is taken as a conversion coefficient; NIR is the reflectance in the near infrared band; r is the reflectivity of the visible red wave band;
E. the T factor acquisition method comprises the following steps:
indoor interpretation and field rechecking are carried out on the remote sensing image, and an E factor and a T factor are extracted by referring to the empirical value.
5. The method for correcting the water and soil loss analog value space-time weighting based on the field station according to claim 4, is characterized in that: the remote sensing image map is interpreted indoors to obtain a current land utilization map, and the value is assigned according to the current land utilization classification (GB/T21010-2017); rechecking and correcting the interpretation result by utilizing national land utilization annual change survey data, forest land grassland resource survey data, water and soil conservation key engineering data and artificial water and soil loss data; and assigning values to the interpretation graph according to the national rotation district name and code of the Classification of Water and soil conservation measures to obtain corresponding E factors and T factors, and finally converting the E factors and T factors into factor grid data.
6. The method for correcting the space-time weighting of the water and soil loss analog value based on the field station according to any one of claims 1 to 5, wherein: in step 3, the soil erosion modulus correction model is as follows:
WASLi＝β0(ui,vi,ti)+∑kβk(ui,vi,ti)Aik+εi,i＝1,2,…n；
in the formula, AikIs an independent variable and represents a computer simulation value of a soil erosion modulus under k CSLE models of the ith observation point; i represents the ith observation point in the soil erosion monitoring station; WASLiThe dependent variable represents the actual measured value of the soil erosion at the ith point and is a dependent variable; (u)i,vi,ti) Representing geographic latitude and longitude coordinates and time coordinates; beta is ak(ui,vi,ti) Regression coefficients representing the kth independent variable of a sample point i; beta is a0(ui,vi,ti) Is a constant term of the equation; epsiloniIs a random error term obeying epsiloni～N(0,σ2) And Cov (. epsilon.)i,εj) When i ≠ j, 0;
the above model is simplified as follows:
WASLi＝Aikβ+ε；
Figure FDA0003317690810000041
ideal value of regression coefficient
Figure FDA0003317690810000042
The calculation is as follows:
Figure FDA0003317690810000043
in the formula (I), the compound is shown in the specification,
Figure FDA0003317690810000044
the ideal value of the regression parameter of the ith observation point is represented; a. theTDenotes the independent variable AikThe transposed matrix of (2); w (u)i,vi,ti)＝diag(wi1,wi2,…,win) Is (u)i,vi,ti) A spatiotemporal weight matrix of (c);
measured value WASL of ith pointiIdeal value of
Figure FDA0003317690810000051
Comprises the following steps:
Figure FDA0003317690810000052
in the formula, SWASLThe hat matrix representing the GTWR model.
7. The method for correcting the water and soil loss analog value space-time weighting based on the field station according to claim 6, is characterized in that: the space-time weight matrix W (u)i,vi,ti) Comprises the following steps:
Figure FDA0003317690810000053
Figure FDA0003317690810000054
Figure FDA0003317690810000055
Figure FDA0003317690810000056
Figure FDA0003317690810000057
Figure FDA0003317690810000058
Figure FDA0003317690810000059
in the formula, WijRepresenting spatio-temporal weight matrix elements;
Figure FDA00033176908100000510
representing the distance from the observation point i to the regression point j; τ represents a space-time ratio parameter; b represents the bandwidth, a parameter describing the non-negative attenuation relationship of spatial weight to spatial distance: the larger the bandwidth value is, the slower the weight decays with the increase of the space-time distance; the smaller the bandwidth value is, the faster the weight is attenuated along with the increase of the space-time distance;
the space-time weight matrix W (u)i,vi,ti) Comprises the following steps:
Figure FDA0003317690810000061
in the formula, WijRepresenting spatio-temporal weight matrix elements;
Figure FDA0003317690810000062
representing the distance from the observation point i to the regression point j; b represents the bandwidth, the regression point j is in the range of the bandwidth, the weight of the data point is calculated, and the weight of the part exceeding the bandwidth is totally marked as 0.
8. The method for correcting the water and soil loss analog value space-time weighting based on the field station according to claim 7, is characterized in that: the selection of the bandwidth value adopts cross validation or Chichi information criterion;
the cross validation is that firstly, the space-time distance range of the data is normalized, the value range of the bandwidth b is given, and then the value range is searched in a traversing way to obtain the optimal bandwidth;
Figure FDA0003317690810000063
in the formula (I), the compound is shown in the specification,
Figure FDA0003317690810000064
representing an ideal value of water and soil loss in the bandwidth of the ith observation point, and only performing regression parameter calculation according to data around the observation point without including the i observation point in order to prevent an overfitting phenomenon when estimating the regression parameters; drawing different bandwidths and different CVs into a trend line, and finding out the corresponding optimal bandwidth when the CV value is minimum;
the Red pool information criterion AICc(b)：
Figure FDA0003317690810000065
Figure FDA0003317690810000066
In the formula, n represents the number of sample points;
Figure FDA0003317690810000067
representing a GTWR model maximum likelihood variance estimate; tr (S)WASL) Representing hat matrix SWASLThe trace of (2); RSS represents the sum of the squares of the residuals.
9. The method according to claims 1-8, obtaining initial values of coefficients of a soil erosion modulus correction model; introducing a particle swarm algorithm on the basis of the GTWR model, setting a space threshold value, and performing random space search near the initial value of the coefficient of the correction model to obtain the optimal solution of the coefficient of the correction model;
the particle swarm algorithm is as follows:
Figure FDA0003317690810000068
Figure FDA0003317690810000069
Figure FDA00033176908100000610
Figure FDA00033176908100000611
Figure FDA0003317690810000071
Figure FDA0003317690810000072
wherein, i is 1,2, …, m, which represents the ith group of particles, i.e. the ith group of coefficient initial values each includes n regression coefficients β, a bandwidth parameter b and a space-time ratio parameter τ in this embodiment; d is 1,2, …, D, which represents the D-th dimension;
Figure FDA0003317690810000073
bits representing initial values of the ith set of coefficientsPlacing;
Figure FDA0003317690810000074
indicating the location where the ith set of coefficient initial values has experienced the best;
Figure FDA0003317690810000075
represents the best position of the whole particle group;
Figure FDA0003317690810000076
a speed representing an initial value of the ith set of coefficients; w is a non-negative number, called the inertia factor; r is1And r2Is [0,1 ]]Random numbers transformed within a range; alpha is a constraint factor and controls the weight of the speed; c. C1And c2Is a learning factor, i.e., an acceleration constant.
10. A method of authenticating the method of claims 1-9, wherein: constructing a soil erosion modulus correction model of a certain area for m years; the soil erosion modulus A under the CSLE model of the (m + 1) th yearmInputting a correction equation to obtain a corrected water and soil loss ideal value of the (m + 1) th year
Figure FDA0003317690810000077
And comparing with the actual measured value WASL of the soil erosion of the place in the m +1 yearm+1Carrying out comparison and inspection; if it is
Figure FDA0003317690810000078
And WASLm+1And (4) approximating, predicting the future water and soil loss amount of the ground by using the soil erosion modulus correction model.
CN202111236015.4A 2021-10-22 2021-10-22 Water and soil loss analog value space-time weighting correction method based on field station Pending CN113987778A (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111236015.4A CN113987778A (en) 2021-10-22 2021-10-22 Water and soil loss analog value space-time weighting correction method based on field station

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111236015.4A CN113987778A (en) 2021-10-22 2021-10-22 Water and soil loss analog value space-time weighting correction method based on field station

Publications (1)

Publication Number Publication Date
CN113987778A true CN113987778A (en) 2022-01-28

Family

ID=79740603

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111236015.4A Pending CN113987778A (en) 2021-10-22 2021-10-22 Water and soil loss analog value space-time weighting correction method based on field station

Country Status (1)

Country Link
CN (1) CN113987778A (en)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115906717A (en) * 2023-03-09 2023-04-04 深圳市睿拓新科技有限公司 Erosion calculation method and system for water and soil conservation

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115906717A (en) * 2023-03-09 2023-04-04 深圳市睿拓新科技有限公司 Erosion calculation method and system for water and soil conservation

Similar Documents

Publication Publication Date Title
CN110751094B (en) Crop yield estimation method based on GEE comprehensive remote sensing image and deep learning method
CN108764255B (en) Method for extracting winter wheat planting information
CN110610054B (en) Method and system for constructing cuboid inversion model of soil humidity
CN111368736B (en) Rice refined estimation method based on SAR and optical remote sensing data
CN107491758B (en) Yangtze river basin water body information extraction and space coding method thereof
Montzka et al. Modelling the water balance of a mesoscale catchment basin using remotely sensed land cover data
CN110927120B (en) Early warning method for vegetation coverage
CN106199627A (en) Grassland vegetation parameter acquiring method
CN115204691B (en) Urban artificial heat emission estimation method based on machine learning and remote sensing technology
CN111582575A (en) Method for identifying urban thermal environment formation development leading factors under multiple space-time scales
Qing-Ling et al. Topographical effects of climate data and their impacts on the estimation of net primary productivity in complex terrain: A case study in Wuling mountainous area, China
Li et al. Mapping spatiotemporal decisions for sustainable productivity of bamboo forest land
Lou et al. An effective method for canopy chlorophyll content estimation of marsh vegetation based on multiscale remote sensing data
CN114254802B (en) Prediction method for vegetation coverage space-time change under climate change drive
CN105930633A (en) Method for forecasting urban heat island effect
CN113987778A (en) Water and soil loss analog value space-time weighting correction method based on field station
Chen et al. Unravelling the multilevel and multi-dimensional impacts of building and tree on surface urban heat islands
CN113297904A (en) Alpine grassland biomass estimation method and system based on satellite driving model
CN107576399A (en) Towards bright the temperature Forecasting Methodology and system of MODIS forest fire detections
Wan et al. Contribution degree of different surface factors in urban interior to urban thermal environment
CN104899464B (en) A kind of sampling study machine remote sensing quantitative inversion method under adaptation noise conditions
CN111723525A (en) PM2.5 inversion method based on multi-source data and neural network model
CN110705010A (en) Remote sensing-based next-day night surface heat island simulation method
CN115527108A (en) Method for rapidly identifying water and soil loss artificial disturbance plots based on multi-temporal Sentinel-2
CN110287915B (en) City impervious bed extraction method based on Landsat remote sensing image

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination