Method for correcting soil erodibility factor in wind erosion model

Abstract

The application discloses a soil erodibility factor correction calculation method in a wind erosion model, relates to the technical field of wind erosion model correction, and is used for solving the problems that the existing unified empirical formula produces systematic deviation and nonlinear deviation which is difficult to be self-adaptively compensated under different soil texture combinations. The application builds a segmented regression and nonlinear deviation identification mechanism, calculates a linear correction function with the minimum residual sum of squares in each segmented interval based on a scatter diagram of a planar simulation erodibility factor and a point measurement value, further determines an inflection point and a deviation direction, combines leave-one-out cross validation and spatial interpolation confidence marking, generates a segmented correction coefficient table and a reliability diagram, and thus local overcorrection or insufficient correction caused by global unified correction is avoided, and in a complex soil spatial heterogeneity scene, the simulation accuracy of the erodibility factor can be improved, the evaluation confidence is improved, and the robustness of wind erosion prediction is improved.

Description

Technical Field

[0001] This invention relates to the field of wind erosion model correction technology, and more specifically, to a method for correcting the soil erodibility factor in a wind erosion model. Background Technology

[0002] In the field of soil wind erosion assessment, wind erosion models based on areal soil property raster data have been widely used. These models typically extract sand, silt, clay, and organic matter content from soil property databases using computer programs, and calculate the soil erodibility factor value for each raster according to preset texture correction and particle size weighting rules. This value is used to predict regional wind erosion risk and formulate windbreak and sand-fixing measures. Existing technical solutions focus on optimizing the accuracy and computational efficiency of soil component extraction, obtaining areal simulated erodibility factor values ​​through unified empirical formulas. The technical implementation involves the coupled processing of geographic information systems, soil physical parameters, and meteorological factors. Firstly, this type of method can quickly generate erodibility factor distribution maps on a large scale; secondly, standardized processes reduce the cost and time of manual sampling.

[0003] However, existing technologies have significant limitations in simulation accuracy. First, there is a scale difference between the resolution of the areal soil property raster data and the measured sample data, leading to varying degrees of systematic shifts in empirical formulas under different soil texture combinations. Second, the spatial heterogeneity of soil components makes it impossible for the original empirical model to accurately capture nonlinear biases in local areas; that is, the direction and magnitude of biases vary across different numerical intervals, and the inflection point location is difficult to eliminate with a single global correction parameter. Furthermore, due to the lack of a segmented identification and compensation mechanism for the nonlinear relationship between measured values ​​and areal simulation values, existing methods cannot adaptively adjust the correction intensity, often resulting in overcorrection in some areas and undercorrection in others. Ultimately, the accumulation of biases significantly reduces the reliability of the wind erosion model output, weakening the scientific basis for soil loss prediction and control decisions based on erodibility factor maps. Summary of the Invention

[0004] In order to overcome the above-mentioned defects of the prior art, embodiments of the present invention provide a method for correcting the calculation of soil erodibility factor in a wind erosion model to solve the problems mentioned in the background art.

[0005] To achieve the above objectives, the present invention provides the following technical solution: A method for correcting the calculation of soil erodibility factors in a wind erosion model, comprising the following steps: S1. Extract the soil components of the grid where each field survey point is located from the isometric soil property raster data, substitute them into the original wind-eroded soil erosibility factor calculation method to calculate the isometric simulated erosibility factor value, and at the same time calculate the point measured erosibility factor value from the soil sample analysis results of the same geographical location. Pair the isometric simulated erosibility factor value and the point measured erosibility factor value according to the geographical coordinates of the field survey points to form an erosibility factor comparison, and obtain the isometric simulated raster map composed of the isometric simulated erosibility factor value. S2. Based on the corrosivity factor comparison, a scatter plot of the simulated corrosivity factor value of the area and the measured corrosivity factor value of the point is drawn. The nonlinear offset pattern between the simulated corrosivity factor value of the area and the measured corrosivity factor value of the point is identified by piecewise regression, and the deviation direction and inflection point position of the numerical range of the simulated corrosivity factor value of the area are determined. S3. Based on the nonlinear offset mode, the numerical range of the simulated erodibility factor value is divided into different numerical segments. A linear translation correction function is fitted to each numerical segment. The regression coefficient of the correction function for each segment is determined by the least squares method based on the erodibility factor comparison, forming a segmented correction coefficient table and the corresponding numerical range criteria. S4. Apply the segmented correction coefficient table and numerical interval criteria to each grid cell in the areal simulation grid map. That is, first determine the numerical segment to which each grid cell belongs based on the areal simulation erodibility factor value, then look up the segmented correction coefficient table and call the regression coefficient of the correction function of the corresponding numerical segment to calculate and obtain the preliminary corrected spatial distribution map of soil erodibility factor. S5. Using the leave-one-out cross-validation method, each field survey sampling point is used as a validation point in turn. Based on the control of soil erodibility factors after removing validation points, the correction function is re-determined and the correction value of the validation point is generated. The error between the correction value and the measured value is calculated. The errors of all validation points are summarized. Confidence labels are generated through spatial interpolation and superimposed on the preliminary corrected spatial distribution map of soil erodibility factors. The final spatial distribution map of soil erodibility factors is then output.

[0006] In a preferred embodiment, soil components of the grid containing each field survey sampling point are extracted from the areal soil property raster data. These components are then substituted into the original method for calculating the erodibility factor of wind-eroded soil to calculate the simulated areal erodibility factor value. Simultaneously, the measured erodibility factor value at each point is calculated from the analysis results of soil samples from the same geographical location. The simulated areal erodibility factor value and the measured erodibility factor value at each point are paired according to the geographical coordinates of the field survey sampling points to form an erodibility factor comparison. This results in an areal simulation raster map composed of the simulated areal erodibility factor values, specifically including: Soil components are extracted from the raster data of areal soil properties. The extracted soil components are then organized according to the data format specified in the original method for calculating the erodibility factor of wind-eroded soil to form an input parameter set containing sand content, silt content, clay content and organic matter content. Then, according to the preset soil texture correction rules, organic matter content correction rules and particle size weighting rules in the original method for calculating the erodibility factor of wind-eroded soil, the corresponding calculation steps are executed in sequence to obtain the areal simulated erodibility factor value for each raster. Based on the geographic coordinates of the field survey points, the grid where each field survey point is located is located from the areal simulated erodibility factor value corresponding to each grid, and the areal simulated erodibility factor value of the grid is extracted as the areal simulated erodibility factor value corresponding to the field survey point. From the analysis results of soil samples from the same geographical location, the test data of sand content, silt content, clay content and organic matter content are extracted and organized according to the data format specified by the on-site measured erodibility factor calculation method to form an input parameter set containing sand content, silt content, clay content and organic matter content; Based on the preset soil texture correction rules, organic matter content correction rules, and particle size weighting rules in the method for calculating the erodibility factor at the sampling point, the corresponding calculation steps are executed in sequence to obtain the measured erodibility factor value for each field survey sampling point. The simulated erodibility factor value of each field survey point is paired with the measured erodibility factor value at the point according to common geographical coordinates to form a one-to-one erodibility factor comparison. Organize the planar simulated erodibility factor values ​​of all rasters according to their original raster positions to generate a planar simulated raster map composed of the planar simulated erodibility factor values.

[0007] In a preferred embodiment, based on the corrosivity factor comparison, a scatter plot of the simulated corrosivity factor values ​​for the area and the measured corrosivity factor values ​​for the point is plotted, specifically including: The simulated erosibility factor value and the measured erosibility factor value of each field survey point were read from the erosibility factor comparison. The simulated erosibility factor value was plotted on the x-axis and the measured erosibility factor value was plotted on the y-axis in a two-dimensional coordinate system.

[0008] In a preferred embodiment, the nonlinear offset pattern between the simulated erosibility factor value of the areal surface and the measured erosibility factor value at a point is identified by piecewise regression, and the deviation direction and inflection point position of the numerical range of the simulated erosibility factor value of different areal surfaces are determined, specifically including: A candidate segment number sequence is set for the scatter plot. For each candidate segment number in the candidate segment number sequence, the entire numerical range of the simulated erodibility factor value is divided into a corresponding number of continuous numerical intervals. Linear regression is performed independently in each numerical interval to minimize the sum of squared residuals. The total sum of squared residuals of all numerical intervals is calculated. The candidate segment number with the smallest total sum of squared residuals is selected as the final segment number. The boundary points of each numerical interval of the final segment number are determined as inflection points. For each numerical interval, calculate the difference between the mean of the measured erosibility factor values ​​of all field survey points within the numerical interval and the mean of the simulated erosibility factor values ​​of the areal area. If the difference is positive, the deviation direction is positive; if the difference is negative, the deviation direction is negative. Output the range of each numerical interval, the corresponding deviation direction, and the inflection point position to form a non-linear offset mode.

[0009] In a preferred embodiment, based on the nonlinear offset mode, the numerical range of the simulated erodibility factor values ​​for the planar surface is divided into different numerical segments. A linear translation correction function is fitted to each numerical segment, and the regression coefficients of the correction function for each segment are determined by the least squares method based on the erodibility factor comparison, forming a segmented correction coefficient table and corresponding numerical range criteria, specifically including: Using the inflection point positions determined in the nonlinear migration mode as segment boundaries, the entire numerical range of the surface simulation erodibility factor value is divided into different continuous numerical segments, and the deviation direction corresponding to each numerical segment is directly used as the attribute of the numerical segment. For each numerical segment, the areal simulated erodibility factor of all field survey points within the numerical segment is used as the independent variable, and the corresponding measured erodibility factor at the corresponding point is used as the dependent variable. A straight-line correction function is fitted to minimize the sum of squared residuals of all data points within the numerical segment to the fitted straight line. The constant term and linear term coefficients of the straight line are determined by the least squares method and used as the regression coefficients of the numerical segment. The numerical range of each numerical segment, the corresponding regression coefficient, and the deviation direction of that numerical segment are organized in tabular form as a segmented correction coefficient table. At the same time, the position of each inflection point is recorded as a numerical range criterion to determine the numerical segment to which the area simulation erodibility factor value of each grid belongs.

[0010] In a preferred embodiment, the segmented correction coefficient table and numerical interval criteria are applied to each grid cell in the areal simulation raster map. Specifically, the numerical segment to which each grid cell belongs is first determined based on its areal simulation erodibility factor value. Then, the regression coefficient of the correction function for the corresponding numerical segment is retrieved from the segmented correction coefficient table for calculation, resulting in a preliminary corrected spatial distribution map of the soil erodibility factor. This includes: Read the area simulation erodibility factor value corresponding to each grid in the area simulation raster map. Based on the grid, compare the area simulation erodibility factor value of the grid with the numerical interval criteria recorded in the piecewise correction coefficient table one by one to determine the numerical segment to which the area simulation erodibility factor value belongs, and obtain the regression coefficient and deviation direction corresponding to the numerical segment in the piecewise correction coefficient table. Multiply the area simulation erodibility factor value of the raster by the obtained first-order term coefficient to obtain an intermediate product. Then add the intermediate product to the obtained constant term to obtain the correction value of the raster. At the same time, mark the obtained deviation direction as positive or negative and record it as the offset attribute of the raster in the corresponding attribute table. After completing the above operations on all grids, the correction values ​​of all grids are arranged and stored according to their original grid positions to generate a preliminary corrected spatial distribution map of soil erodibility factors.

[0011] In a preferred embodiment, a leave-one-out cross-validation method is used. Each field survey sample point is sequentially used as a validation point. Based on the erosibility factor comparison after removing validation points, a correction function is re-determined and a validation point correction value is generated. The error between the correction value and the measured value is calculated, and the errors of all validation points are summarized. Specifically, this includes: Select one field survey point from the erosibility factor control as a verification point, and combine all the remaining field survey points to form a new erosibility factor control. Based on the new corrosivity factor comparison, a scatter plot of the simulated corrosivity factor value and the measured corrosivity factor value at the point was redrawn. The new nonlinear migration pattern was identified by piecewise regression, the new inflection point position and deviation direction were determined, and the numerical segments were re-divided according to the new nonlinear migration pattern. A linear translation correction function was refitted for each numerical segment, and the new regression coefficients of each segment were determined by the least squares method, forming a new piecewise correction coefficient table and the corresponding numerical interval criteria. Based on the new segmented correction coefficient table and numerical interval criteria, the areal simulated erodibility factor values ​​of the field survey samples selected as verification points are judged to determine the numerical segment to which they belong, and the regression coefficients of the corresponding numerical segment are called to calculate the correction value of the verification point. Calculate the difference between the correction value of the verification point and the corresponding measured erosibility factor value, and record the error of the verification point.

[0012] In a preferred embodiment, confidence markers are generated through spatial interpolation and superimposed onto the preliminarily corrected spatial distribution map of soil erodibility factors to output the final spatial distribution map of soil erodibility factors, specifically including: After iterating through all field survey points, all recorded errors are summarized to obtain an error set. Spatial interpolation is performed on all errors in the error set to generate confidence markers that reflect the reliability of the correction results. The confidence level markers are spatially aligned with the preliminary revised spatial distribution map of soil erodibility factors. The pixel values ​​of the confidence level markers are used as transparency weights and are superimposed on the corresponding pixels of the preliminary revised spatial distribution map of soil erodibility factors using a raster band synthesis method to output the final spatial distribution map of soil erodibility factors.

[0013] Compared with the prior art, the present invention has the following beneficial effects: 1. By identifying the nonlinear shift pattern between simulated and measured values ​​of soil erodibility factors through piecewise regression, this method accurately captures and corrects the deviation direction and inflection point location across different numerical intervals using piecewise linear regression. Based on scatter plot analysis of erodibility factor comparisons and a piecewise correction function constructed using the principle of minimizing the sum of squared residuals, it effectively identifies the nonlinear influence of soil texture and organic matter content on erodibility factors across different intervals. Combined with independently fitted linear translation correction coefficients for each numerical segment, and through the joint application of a piecewise correction coefficient table and numerical interval criteria, systematic shifts in the simulated soil erodibility raster map are compensated grid-by-grid, resulting in a preliminary corrected spatial distribution map of soil erodibility factors. Compared to existing technologies, this method can adaptively adjust the correction intensity under different soil property combinations, reduce nonlinear deviations through piecewise linear regression, significantly improve the consistency between simulated and measured values ​​of soil erodibility factors, and suppress local over-correction or under-correction caused by a single global correction.

[0014] 2. By combining leave-one-out cross-validation with spatial interpolation confidence markers, this method solves the problems of random selection of validation points and lack of visualization of the reliability of correction results in traditional correction methods. It is based on an iterative process of successively removing validation points, refitting the piecewise correction function, and calculating the correction error. The constructed error set can reflect the credibility of the correction result at each spatial location in real time. By synthesizing the confidence markers generated by spatial interpolation with the raster bands of the preliminary correction map, a transparent superposition display of confidence is achieved. This ensures that the final output spatial distribution map of soil erodibility factors is transparent or semi-transparent in low-confidence areas to prompt users to use it with caution, while high-confidence areas are completely opaque to show reliable results. Compared with existing technologies, this method can quantify the uncertainty of correction results in space. By dynamically adjusting the display effect through the transparency weight of the confidence markers, it suppresses the accumulation of model output errors caused by blindly using correction results, and ensures the long-term robustness of wind erosion assessment in complex geographical environments. Attached Figure Description

[0015] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0016] Figure 1 This is a flowchart of a method for correcting the soil erosibility factor in a wind erosion model according to the present invention. Detailed Implementation

[0017] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0018] Example 1: Figure 1 A flowchart of a method for correcting the calculation of soil erodibility factor in a wind erosion model according to the present invention is provided, which includes the following steps: S1. Extract the soil components of the grid where each field survey point is located from the isometric soil property raster data, substitute them into the original wind-eroded soil erosibility factor calculation method to calculate the isometric simulated erosibility factor value, and at the same time calculate the point measured erosibility factor value from the soil sample analysis results of the same geographical location. Pair the isometric simulated erosibility factor value and the point measured erosibility factor value according to the geographical coordinates of the field survey points to form an erosibility factor comparison, and obtain the isometric simulated raster map composed of the isometric simulated erosibility factor value. S2. Based on the corrosivity factor comparison, a scatter plot of the simulated corrosivity factor value of the area and the measured corrosivity factor value of the point is drawn. The nonlinear offset pattern between the simulated corrosivity factor value of the area and the measured corrosivity factor value of the point is identified by piecewise regression, and the deviation direction and inflection point position of the numerical range of the simulated corrosivity factor value of the area are determined. S3. Based on the nonlinear offset mode, the numerical range of the simulated erodibility factor value is divided into different numerical segments. A linear translation correction function is fitted to each numerical segment. The regression coefficient of the correction function for each segment is determined by the least squares method based on the erodibility factor comparison, forming a segmented correction coefficient table and the corresponding numerical range criteria. S4. Apply the segmented correction coefficient table and numerical interval criteria to each grid cell in the areal simulation grid map. That is, first determine the numerical segment to which each grid cell belongs based on the areal simulation erodibility factor value, then look up the segmented correction coefficient table and call the regression coefficient of the correction function of the corresponding numerical segment to calculate and obtain the preliminary corrected spatial distribution map of soil erodibility factor. S5. Using the leave-one-out cross-validation method, each field survey sampling point is used as a validation point in turn. Based on the control of soil erodibility factors after removing validation points, the correction function is re-determined and the correction value of the validation point is generated. The error between the correction value and the measured value is calculated. The errors of all validation points are summarized. Confidence labels are generated through spatial interpolation and superimposed on the preliminary corrected spatial distribution map of soil erodibility factors. The final spatial distribution map of soil erodibility factors is then output.

[0019] Soil components were extracted from the raster data of areal soil properties, and the soil components of each field survey sampling point were substituted into the original method for calculating the areal simulated erosibility factor. Simultaneously, the measured erosibility factor values ​​were calculated from soil samples at the same geographical location. The areal simulated erosibility factor values ​​and the measured erosibility factor values ​​were paired according to the geographical coordinates of the field survey sampling points to form an erosibility factor comparison. A areal simulated raster map composed of the areal simulated erosibility factor values ​​was then obtained. The specific implementation is as follows: Soil components recorded in each raster cell were extracted from the areal soil property raster data. This raster data consisted of pre-prepared raster cells covering the study area with a resolution of 30 meters and different bands. Each band stored the contents of sand, silt, clay, and organic matter, respectively. For each raster cell, the values ​​corresponding to its four bands were extracted and processed according to the data format specified in the original method for calculating the erodibility factor of wind-eroded soil. This resulted in an input parameter set containing the contents of sand, silt, clay, and organic matter. The original method for calculating the erodibility factor of wind-eroded soil adopted an equivalent lookup table simplification method based on the erodibility factor calculation formula in the RWEQ model.

[0020] Then, according to the preset soil texture correction rules in this method, the soil texture factor is first calculated based on the mass percentages of sand, silt, and clay content. This factor is obtained through table lookup or interpolation. Specifically, when the clay content is less than or equal to 15% and the sand content is greater than or equal to 70%, the texture factor is 0.16; when the clay content is greater than 15% and less than or equal to 30%, the texture factor is 0.27. In other cases, the soil texture category is first determined according to the mass percentages of sand, silt, and clay, based on the USDA soil texture triangle, such as loam and sandy clay loam. Each soil texture category corresponds to a pre-determined soil texture factor value. This correspondence is stored in a local lookup table file, for example, taken from the discretization result of the RWEQ original formula under standard texture composition. If the field survey sample points fall on the boundary between two categories, linear interpolation is performed to obtain the soil texture factor.

[0021] Then, according to the organic matter content correction rule, the organic matter content is divided by 100 and multiplied by a constant of 0.2 to obtain the organic matter correction factor. Then, according to the particle size weighting rule, the soil texture factor and the organic matter correction factor are weighted and summed, with weighting coefficients of 0.6 and 0.4, respectively, to obtain the areal simulated erodibility factor value corresponding to each grid. The constants and weighting coefficients involved in the above correction rules are all obtained by equivalent transformation of the standard parameters of the RWEQ model. Those skilled in the art can directly obtain them by consulting the RWEQ model technical documents or commonly used simplified tables.

[0022] The geographic coordinates, such as latitude and longitude or projected coordinates, of each field survey point are obtained from the field survey records. In the raster data of the calculated areal simulated erosibility factor values, the specific raster where each field survey point is located is located according to its geographic coordinates, i.e., the raster where the geographic coordinates fall, and the areal simulated erosibility factor value of that raster is extracted as the areal simulated erosibility factor value corresponding to that field survey point. All field survey points are processed in sequence to obtain the areal simulated erosibility factor value corresponding to each field survey point.

[0023] For soil samples from the same geographical locations as the field survey sampling points, laboratory analysis data on sand content, silt content, clay content, and organic matter content were extracted. These data were then organized into an input parameter set containing sand content, silt content, clay content, and organic matter content, according to the data format specified in the on-site measured erosibility factor calculation method. The soil texture correction rules, organic matter content correction rules, and particle size weighting rules preset in this on-site measured erosibility factor calculation method are exactly the same as those used in the above-mentioned areal simulation erosibility factor calculation. The same lookup rules, correction formulas, and weighting coefficients are used to execute each rule sequentially. Specifically, the soil texture factor is obtained by looking up the table file and interpolating the boundary values ​​based on the content of sand, silt, and clay. Then, the organic matter content is divided by 100 and multiplied by 0.2 to obtain the organic matter correction factor. Finally, the soil texture factor is multiplied by 0.6 and the organic matter correction factor is multiplied by 0.4 to obtain the measured erodibility factor value of the soil sample. The above steps are repeated for the soil samples of each field survey point to obtain the measured erodibility factor value of each field survey point.

[0024] Based on the corrosivity factor comparison, a scatter plot of the simulated corrosivity factor values ​​for the area and the measured corrosivity factor values ​​for the points was plotted. The specific implementation is as follows: Based on the corrosivity factor comparison, the simulated corrosivity factor value of the area and the measured corrosivity factor value of the point are read sequentially for each field survey sample point. The simulated corrosivity factor value of the area is used as the abscissa and the measured corrosivity factor value of the point is used as the ordinate. A data point is plotted for each field survey sample point in a two-dimensional coordinate system, thereby obtaining a scatter plot covering all field survey sample points.

[0025] By identifying the nonlinear offset pattern between the simulated erosibility factor value for a surface and the measured erosibility factor value for a point using piecewise regression, the deviation direction and inflection point location are determined within different numerical ranges of the simulated erosibility factor value for different surfaces. The specific implementation is as follows: A candidate segment number sequence is set for the scatter plot. This candidate segment number sequence consists of a number of pre-defined positive integers, such as an integer sequence from 1 to 10. For each candidate segment number in the candidate segment number sequence, the entire numerical range of the areal simulated erodibility factor value, that is, the minimum to maximum value of the areal simulated erodibility factor value in all field survey samples, is divided into a number of continuous numerical intervals equal to the number of candidate segments. The boundary points of each numerical interval are determined at equal intervals or based on the data distribution.

[0026] Within each numerical interval, linear regression is independently performed on all field survey points within that interval to minimize the sum of squared residuals for all field survey points within that interval—that is, the sum of squared differences between the measured erosibility factor values ​​at each field survey point and the corresponding predicted values ​​on the regression line for that interval. The sum of squared residuals within that interval is then calculated. Finally, the sums of squared residuals from all numerical intervals are summed to obtain the total sum of squared residuals for that candidate number of segments. This total sum of squared residuals is defined as follows:

[0027] in Let j be the isometric simulated erosibility factor value for the j-th field survey sampling point. The corresponding measured corrosion factor value is given. and Let be the intercept and slope of the linear regression within the i-th numerical interval. Perform the above operation sequentially on all candidate segment numbers to obtain the total sum of squared residuals corresponding to each candidate segment number. Select the candidate segment number with the smallest total sum of squared residuals as the final segment number, i.e.:

[0028] in The maximum number of segments is preset, and the boundary points of each numerical interval under this final number of segments are determined as inflection points, that is, the turning points of the erosibility factor value in the planar simulation at the boundaries of different numerical intervals. The above operation is performed on all candidate number of segments in sequence to obtain the total residual sum of squares corresponding to each candidate number of segments. The candidate number of segments with the smallest total residual sum of squares is selected as the final number of segments, and the boundary points of each numerical interval under this final number of segments are determined as inflection points, that is, the turning points of the erosibility factor value in the planar simulation at the boundaries of different numerical intervals.

[0029] After determining the final number of segments and the dividing points (inflection points) of each interval, for each numerical interval, calculate the arithmetic mean of the measured erosibility factor values ​​of all field survey points within that interval, as well as the arithmetic mean of the simulated erosibility factor values ​​of all field survey points. Subtract the arithmetic mean of the simulated erosibility factor values ​​of the area from the arithmetic mean of the measured erosibility factor values ​​to obtain the difference. If the difference is positive, the deviation direction of that numerical interval is determined to be positive; if the difference is negative, the deviation direction is determined to be negative.

[0030] Finally, the range of each numerical interval is output, namely the minimum and maximum values ​​of the surface simulation erodibility factor in that interval, the deviation direction corresponding to that interval (positive or negative), and the positions of all inflection points (i.e., the specific values ​​of the boundary points of each interval). This information together constitutes the nonlinear migration mode.

[0031] Based on the nonlinear migration mode, the numerical range of the simulated erodibility factor values ​​is divided into different numerical segments. A linear translation correction function is fitted to each segment, and the regression coefficients of the correction function for each segment are determined by the least squares method based on the erodibility factor comparison, forming a segmented correction coefficient table and corresponding numerical range criteria. The specific implementation is as follows: In the nonlinear migration mode, the positions of each inflection point have been determined. These inflection point positions are now used as segmentation boundaries to divide the entire numerical range of the simulated erosibility factor values ​​of the areal into different continuous numerical segments. At the same time, the deviation direction corresponding to each numerical segment in the nonlinear migration mode is directly used as the inherent property of that numerical segment. For each numerical segment, the simulated erosibility factor values ​​of the areal of all field survey points within that numerical segment are used as independent variables, and the measured erosibility factor values ​​of the corresponding points are used as dependent variables to fit a straight-line correction function.

[0032] The mathematical form of the correction function is that the dependent variable equals the constant term plus the product of the coefficient of the linear term and the independent variable. The goal of the fitting is to minimize the sum of squared residuals from all data points in the numerical range to the fitted line. The constant term and the coefficient of the linear term are solved by the least squares method. Specifically, the least squares method obtains the coefficient of the linear term by calculating the ratio of the covariance to the variance between the independent and dependent variables of all data points. Then, the constant term is obtained by subtracting the product of the coefficient of the linear term and the mean of the independent variable from the mean of the dependent variable. Thus, the regression coefficient of the numerical range is determined.

[0033] The numerical range of each numerical segment, the corresponding regression coefficients (constant and linear terms), and the deviation direction of the numerical segment are organized in tabular form as a piecewise correction coefficient table. Each row of the table corresponds to a numerical segment, listing the range, regression coefficient, and deviation direction respectively. At the same time, the position of each inflection point is recorded separately as a numerical range criterion. This criterion is used to determine the numerical segment to which the areal simulated erodibility factor value of each raster belongs. That is, when the raster data is subsequently corrected, the relationship between the areal simulated erodibility factor value of the raster and the magnitude of each inflection point is first used to determine which numerical segment it falls into. Then, the regression coefficient and deviation direction corresponding to the numerical segment are extracted from the piecewise correction coefficient table.

[0034] The piecewise correction coefficient table and numerical interval criteria are applied to each grid cell in the areal simulation raster map. Specifically, the numerical range to which each grid cell belongs is first determined based on its areal simulation erodibility factor value. Then, the regression coefficient of the correction function for the corresponding numerical range is retrieved from the piecewise correction coefficient table for calculation, resulting in a preliminary corrected spatial distribution map of the soil erodibility factor. The specific implementation is as follows: First, the simulated erosibility factor value of each grid cell in the simulated area raster is read. For each grid cell, its simulated erosibility factor value is compared one by one with the numerical interval criteria recorded in the piecewise correction coefficient table. The numerical interval criteria are composed of the positions of each inflection point. This criterion defines the boundary value of each numerical segment. For example, if the numerical segment interval is [a, b), then the judgment condition is that the simulated erosibility factor value of the grid cell is greater than or equal to a and less than b. By comparing the simulated erosibility factor value of the grid cell with all inflection points, it is determined which numerical segment it falls into. Once the numerical segment is determined, the regression coefficients corresponding to the numerical segment are obtained from the piecewise correction coefficient table, including the constant term and the linear term coefficients, as well as the deviation direction of the numerical segment. The deviation direction has been marked as positive or negative in the piecewise correction coefficient table.

[0035] After obtaining the regression coefficients and deviation direction, the raster is corrected. First, the raster's simulated erodibility factor value is multiplied by the obtained linear coefficient to obtain an intermediate product. Then, the intermediate product is added to the obtained constant term to obtain the raster's correction value. At the same time, the obtained deviation direction is directly recorded as the raster's offset attribute, marked as positive or negative, and stored in the attribute table corresponding to the raster. The attribute table may contain fields such as raster identifier, original simulated erodibility factor value, correction value, and offset attribute.

[0036] The above operations are performed sequentially on all rasters until each raster in the areal simulation raster map has completed comparison, calculation and attribute recording. Finally, the correction values ​​of all rasters are arranged and stored according to their original raster positions, that is, the spatial coordinates of the rasters remain unchanged. The correction values ​​are used to replace the original areal simulation erodibility factor values ​​to generate a new raster map as a preliminary corrected spatial distribution map of soil erodibility factors.

[0037] The leave-one-out cross-validation method is adopted. Each field survey sampling point is used as a validation point in turn. Based on the erosability factor comparison after removing validation points, the correction function is re-determined and the validation point correction value is generated. The error between the correction value and the measured value is calculated, and the errors of all validation points are summarized. The specific implementation is as follows: Based on all field survey points, a soil erodibility factor control has been established, and a preliminary revised spatial distribution map of soil erodibility factors has been generated. The reliability of the revised results was evaluated using leave-one-out cross-validation. Specifically, one field survey point was selected sequentially from the erodibility factor control as a validation point, and all remaining field survey points were combined to form a new erodibility factor control. Based on the new erodibility factor control, the operations of plotting scatter plots and piecewise regression were re-executed: a new scatter plot was plotted with the areal simulated erodibility factor values ​​of all points in the new control as the x-axis and the measured erodibility factor values ​​at the points as the y-axis. New nonlinear shift patterns were identified through piecewise regression. That is, for each segment number in the candidate segment number sequence, the entire numerical range of the areal simulated erodibility factor values ​​was divided into an equal number of continuous numerical intervals. Linear regression was performed independently within each interval to minimize the sum of squared residuals. The sum of squared residuals of all intervals was accumulated to obtain the total sum of squared residuals. The candidate segment number with the smallest total sum of squared residuals was selected as the final segment number, and the boundary points of each interval were determined as new inflection points.

[0038] For each numerical interval, the arithmetic mean of the measured corrosivity factor values ​​of all sample points within the interval is calculated, minus the arithmetic mean of the simulated corrosivity factor values ​​of the area. If the difference is positive, the deviation direction is determined to be positive; otherwise, it is negative, thus obtaining a new deviation direction. Based on the new inflection point position and deviation direction, the numerical range of the simulated corrosivity factor values ​​of the area is divided into new continuous numerical segments, and each numerical segment inherits its corresponding deviation direction as an inherent attribute.

[0039] For each numerical segment, the simulated erodibility factor of all samples within the segment is used as the independent variable, and the measured erodibility factor of the points is used as the dependent variable. A correction function in the form of a straight line is fitted using the least squares method. First, the covariance between all independent and dependent variables within the segment is calculated and divided by the variance of the independent variables to obtain the coefficient of the first term. Then, the mean of the dependent variable is calculated and the product of the coefficient of the first term and the mean of the independent variables is subtracted to obtain the constant term. Thus, the new regression coefficients for the numerical segment are obtained. The interval range, regression coefficients, and deviation direction of each numerical segment are organized in tabular form as a new piecewise correction coefficient table, and all inflection point positions are recorded as new numerical interval criteria.

[0040] Based on the new segmented correction coefficient table and the corresponding numerical interval criteria, the areal simulated erodibility factor value of the field survey sample points selected as verification points is judged. Specifically, the areal simulated erodibility factor value of the verification point is compared with the position of each inflection point in the new numerical interval criteria to determine which numerical segment it falls into. For example, if the interval is [a, b) and the value is greater than or equal to a and less than b, it is determined to fall into that numerical segment.

[0041] Extract the regression coefficients corresponding to the numerical segment from the new segmented correction coefficient table, including the linear term coefficient and the constant term, as well as the deviation direction of the numerical segment. Then, perform correction calculations on the verification point. First, multiply the simulated erodibility factor value of the verification point by the obtained linear term coefficient to obtain an intermediate product. Then, add the intermediate product to the obtained constant term to obtain the correction value of the verification point. Next, calculate the difference between the correction value of the verification point and the corresponding measured erodibility factor value. Specifically, subtract the measured erodibility factor value of the point from the correction value of the verification point. The result is taken as the error of the verification point, and the error value is recorded. The recorded error may include the geographical coordinates of the verification point.

[0042] Confidence markers are generated through spatial interpolation and superimposed onto the preliminarily corrected spatial distribution map of soil erodibility factors to output the final spatial distribution map of soil erodibility factors. The specific implementation is as follows: In the manner described above, all field survey points in the erosibility factor control are iterated through in an iterative manner. That is, each field survey point is used as a verification point in turn, and the remaining points form a new control. The piecewise correction coefficient table and numerical interval criteria are regenerated, the correction error of the verification point is calculated and recorded, and after all field survey points have been selected once and the error calculation has been completed, all recorded errors are summarized to obtain an error set. This error set contains an error value and its corresponding geographic coordinates at each field survey point.

[0043] Spatial interpolation is performed on all error values ​​in the error set, for example, using inverse distance weighted interpolation or ordinary kriging interpolation. The error value is used as the interpolation field, and the study area is used as the interpolation range to generate a continuous raster layer covering the entire study area. Each pixel value of this raster layer represents the possible error size of the correction result at that location. The error values ​​are converted into confidence labels through normalization or grading. Finally, a confidence label raster map reflecting the reliability of the correction result is obtained. Its pixel value represents the confidence level. The confidence level uses a value between 0 and 1, where 1 represents completely reliable and 0 represents completely unreliable.

[0044] The confidence-marked raster map is spatially aligned with the preliminary corrected spatial distribution map of soil erodibility factors. Then, the pixel value of the confidence mark is used as the transparency weight. Specifically, for each pixel in the preliminary corrected spatial distribution map of soil erodibility factors, if the confidence mark value is 1, the pixel is completely opaque; if the confidence mark value is 0, the pixel is completely transparent; and if it is in between, the transparency is set proportionally. The preliminary corrected spatial distribution map of soil erodibility factors is used as the base layer by raster band synthesis, and the transparency weight is superimposed to output the final spatial distribution map of soil erodibility factors.

[0045] The above embodiments can be implemented, in whole or in part, by software, hardware, firmware, or any other combination thereof. When implemented using software, the above embodiments can be implemented, in whole or in part, in the form of a computer program product.

[0046] Those skilled in the art will recognize that the modules and algorithm steps of the various examples described in conjunction with the embodiments disclosed herein can be implemented in electronic hardware, computer software, and a combination of electronic hardware and software. Whether these functions are implemented in hardware or software depends on the specific application and design constraints of the technical solution. Those skilled in the art can use different methods to implement the described functions for each specific application, but such implementation should not be considered beyond the scope of this application.

[0047] In addition, the functional modules in the various embodiments of this application can be integrated into one processing module, or each module can exist physically separately, or two or more modules can be integrated into one module.

[0048] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

[0049] In conclusion, the above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for correcting the calculation of soil erodibility factors in a wind erosion model, characterized in that, Includes the following steps: S1. Extract the soil components of the grid where each field survey point is located from the isometric soil property raster data, substitute them into the original wind-eroded soil erosibility factor calculation method to calculate the isometric simulated erosibility factor value, and at the same time calculate the point measured erosibility factor value from the soil sample analysis results of the same geographical location. Pair the isometric simulated erosibility factor value and the point measured erosibility factor value according to the geographical coordinates of the field survey points to form an erosibility factor comparison, and obtain the isometric simulated raster map composed of the isometric simulated erosibility factor value. S2. Based on the corrosivity factor comparison, a scatter plot of the simulated corrosivity factor value of the area and the measured corrosivity factor value of the point is drawn. The nonlinear offset pattern between the simulated corrosivity factor value of the area and the measured corrosivity factor value of the point is identified by piecewise regression, and the deviation direction and inflection point position of the numerical range of the simulated corrosivity factor value of the area are determined. S3. Based on the nonlinear offset mode, the numerical range of the simulated erodibility factor value is divided into different numerical segments. A linear translation correction function is fitted to each numerical segment. The regression coefficient of the correction function for each segment is determined by the least squares method based on the erodibility factor comparison, forming a segmented correction coefficient table and the corresponding numerical range criteria. S4. Apply the segmented correction coefficient table and numerical interval criteria to each grid cell in the areal simulation grid map. That is, first determine the numerical segment to which each grid cell belongs based on the areal simulation erodibility factor value, then look up the segmented correction coefficient table and call the regression coefficient of the correction function of the corresponding numerical segment to calculate and obtain the preliminary corrected spatial distribution map of soil erodibility factor. S5. Using the leave-one-out cross-validation method, each field survey sampling point is used as a validation point in turn. Based on the control of soil erodibility factors after removing validation points, the correction function is re-determined and the correction value of the validation point is generated. The error between the correction value and the measured value is calculated. The errors of all validation points are summarized. Confidence labels are generated through spatial interpolation and superimposed on the preliminary corrected spatial distribution map of soil erodibility factors. The final spatial distribution map of soil erodibility factors is then output.

2. The method for correcting the calculation of soil erodibility factor in a wind erosion model according to claim 1, characterized in that: Soil components were extracted from the raster data of areal soil properties, and the soil components of each field survey sampling point were substituted into the original method for calculating the areal simulated erosibility factor. Simultaneously, the measured erosibility factor values ​​were calculated from soil samples at the same geographical location. The areal simulated erosibility factor values ​​and the measured erosibility factor values ​​were paired according to the geographical coordinates of the field survey sampling points to form an erosibility factor comparison. A areal simulated raster map composed of the areal simulated erosibility factor values ​​was then obtained, specifically including: Soil components are extracted from the raster data of areal soil properties. The extracted soil components are then organized according to the data format specified in the original method for calculating the erodibility factor of wind-eroded soil to form an input parameter set containing sand content, silt content, clay content and organic matter content. Then, according to the preset soil texture correction rules, organic matter content correction rules and particle size weighting rules in the original method for calculating the erodibility factor of wind-eroded soil, the corresponding calculation steps are executed in sequence to obtain the areal simulated erodibility factor value for each raster. Based on the geographic coordinates of the field survey points, the grid where each field survey point is located is located from the areal simulated erodibility factor value corresponding to each grid, and the areal simulated erodibility factor value of the grid is extracted as the areal simulated erodibility factor value corresponding to the field survey point. From the analysis results of soil samples from the same geographical location, the test data of sand content, silt content, clay content and organic matter content are extracted and organized according to the data format specified by the on-site measured erodibility factor calculation method to form an input parameter set containing sand content, silt content, clay content and organic matter content; Based on the preset soil texture correction rules, organic matter content correction rules, and particle size weighting rules in the method for calculating the erodibility factor at the sampling point, the corresponding calculation steps are executed in sequence to obtain the measured erodibility factor value for each field survey sampling point. The simulated erodibility factor value of each field survey point is paired with the measured erodibility factor value at the point according to common geographical coordinates to form a one-to-one erodibility factor comparison. Organize the planar simulated erodibility factor values ​​of all rasters according to their original raster positions to generate a planar simulated raster map composed of the planar simulated erodibility factor values.

3. The method for correcting the calculation of soil erodibility factor in a wind erosion model according to claim 1, characterized in that: Based on the corrosivity factor comparison, a scatter plot of the simulated corrosivity factor values ​​for the area and the measured corrosivity factor values ​​for the points was plotted, specifically including: The simulated erosibility factor value and the measured erosibility factor value of each field survey point were read from the erosibility factor comparison. The simulated erosibility factor value was plotted on the x-axis and the measured erosibility factor value was plotted on the y-axis in a two-dimensional coordinate system.

4. The method for correcting the calculation of soil erodibility factor in a wind erosion model according to claim 3, characterized in that: By identifying the nonlinear offset pattern between simulated erosibility factor values ​​for areal surfaces and measured erosibility factor values ​​for specific locations through piecewise regression, the direction of deviation and inflection point location are determined within different numerical ranges of simulated erosibility factor values ​​for different areal surfaces. Specifically, this includes: A candidate segment number sequence is set for the scatter plot. For each candidate segment number in the candidate segment number sequence, the entire numerical range of the simulated erodibility factor value is divided into a corresponding number of continuous numerical intervals. Linear regression is performed independently in each numerical interval to minimize the sum of squared residuals. The total sum of squared residuals of all numerical intervals is calculated. The candidate segment number with the smallest total sum of squared residuals is selected as the final segment number. The boundary points of each numerical interval of the final segment number are determined as inflection points. For each numerical interval, calculate the difference between the mean of the measured erosibility factor values ​​of all field survey points within the numerical interval and the mean of the simulated erosibility factor values ​​of the areal area. If the difference is positive, the deviation direction is positive; if the difference is negative, the deviation direction is negative. Output the range of each numerical interval, the corresponding deviation direction, and the inflection point position to form a non-linear offset mode.

5. The method for correcting the calculation of soil erodibility factor in a wind erosion model according to claim 1, characterized in that: Based on the nonlinear migration mode, the numerical range of the simulated erodibility factor values ​​for the areal region is divided into different numerical segments. A linear translation correction function is fitted to each segment, and the regression coefficients of the correction function for each segment are determined by the least squares method based on the erodibility factor comparison, forming a segmented correction coefficient table and corresponding numerical range criteria, specifically including: Using the inflection point positions determined in the nonlinear migration mode as segment boundaries, the entire numerical range of the surface simulation erodibility factor value is divided into different continuous numerical segments, and the deviation direction corresponding to each numerical segment is directly used as the attribute of the numerical segment. For each numerical segment, the areal simulated erodibility factor of all field survey points within the numerical segment is used as the independent variable, and the corresponding measured erodibility factor at the corresponding point is used as the dependent variable. A straight-line correction function is fitted to minimize the sum of squared residuals of all data points within the numerical segment to the fitted straight line. The constant term and linear term coefficients of the straight line are determined by the least squares method and used as the regression coefficients of the numerical segment. The numerical range of each numerical segment, the corresponding regression coefficient, and the deviation direction of that numerical segment are organized in tabular form as a segmented correction coefficient table. At the same time, the position of each inflection point is recorded as a numerical range criterion to determine the numerical segment to which the area simulation erodibility factor value of each grid belongs.

6. The method for correcting the calculation of soil erodibility factor in a wind erosion model according to claim 1, characterized in that: The piecewise correction coefficient table and numerical interval criteria are applied to each grid cell in the areal simulation raster map. First, the numerical range to which each grid cell belongs is determined based on its areal simulation erodibility factor value. Then, the regression coefficient of the correction function for the corresponding numerical range is retrieved from the piecewise correction coefficient table for calculation, resulting in a preliminary corrected spatial distribution map of the soil erodibility factor. Specifically, this includes: Read the area simulation erodibility factor value corresponding to each grid in the area simulation raster map. Based on the grid, compare the area simulation erodibility factor value of the grid with the numerical interval criteria recorded in the piecewise correction coefficient table one by one to determine the numerical segment to which the area simulation erodibility factor value belongs, and obtain the regression coefficient and deviation direction corresponding to the numerical segment in the piecewise correction coefficient table. Multiply the area simulation erodibility factor value of the raster by the obtained first-order term coefficient to obtain an intermediate product. Then add the intermediate product to the obtained constant term to obtain the correction value of the raster. At the same time, mark the obtained deviation direction as positive or negative and record it as the offset attribute of the raster in the corresponding attribute table. After completing the above operations on all grids, the correction values ​​of all grids are arranged and stored according to their original grid positions to generate a preliminary corrected spatial distribution map of soil erodibility factors.

7. The method for correcting the calculation of soil erodibility factor in a wind erosion model according to claim 1, characterized in that: The leave-one-out cross-validation method was adopted, with each field survey sampling point serving as a validation point in turn. Based on the erosability factor comparison after removing validation points, the correction function was re-determined and the validation point correction value was generated. The error between the correction value and the measured value was calculated, and the errors of all validation points were summarized, specifically including: Select one field survey point from the erosibility factor control as a verification point, and combine all the remaining field survey points to form a new erosibility factor control. Based on the new corrosivity factor comparison, a scatter plot of the simulated corrosivity factor value and the measured corrosivity factor value at the point was redrawn. The new nonlinear migration pattern was identified by piecewise regression, the new inflection point position and deviation direction were determined, and the numerical segments were re-divided according to the new nonlinear migration pattern. A linear translation correction function was refitted for each numerical segment, and the new regression coefficients of each segment were determined by the least squares method, forming a new piecewise correction coefficient table and the corresponding numerical interval criteria. Based on the new segmented correction coefficient table and numerical interval criteria, the areal simulated erodibility factor values ​​of the field survey samples selected as verification points are judged to determine the numerical segment to which they belong, and the regression coefficients of the corresponding numerical segment are called to calculate the correction value of the verification point. Calculate the difference between the correction value of the verification point and the corresponding measured erosibility factor value, and record the error of the verification point.

8. The method for correcting the calculation of soil erodibility factor in a wind erosion model according to claim 7, characterized in that: Confidence markers are generated through spatial interpolation and superimposed onto the initially corrected spatial distribution map of soil erodibility factors, outputting the final spatial distribution map of soil erodibility factors, which specifically includes: After iterating through all field survey points, all recorded errors are summarized to obtain an error set. Spatial interpolation is performed on all errors in the error set to generate confidence markers that reflect the reliability of the correction results. The confidence level markers are spatially aligned with the preliminary revised spatial distribution map of soil erodibility factors. The pixel values ​​of the confidence level markers are used as transparency weights and are superimposed on the corresponding pixels of the preliminary revised spatial distribution map of soil erodibility factors using a raster band synthesis method to output the final spatial distribution map of soil erodibility factors.

Images