Regional Evaluation Method of Risk of Wind-Water Composite Erosion Based on Coupling Coordination Degree

Abstract

The present invention discloses a method for evaluating the risk of wind-water combined erosion in a region based on coupling coordination degree, which specifically includes the following steps: First, determine the research area and collect the shp file of the research area boundary, meteorological data, and NDVI data; calculate the monthly wind erosion climate erosivity and rainfall erosivity of the research area; then calculate the three-dimensional coupling coordination degree of wind erosion climate erosivity, rainfall erosivity, and NDVI; according to the probability statistics theory, calculate the theoretical cumulative probability of the three-dimensional coupling coordination degree of wind erosion climate erosivity, rainfall erosivity, and NDVI to reflect the risk of wind-water combined erosion; according to the Mann-Kendall trend test method, the spatial distribution results of the change trend of the risk of wind-water combined erosion in the research area can be identified.

Description

Technical Field

[0001] The present invention belongs to the technical field of soil erosion risk assessment, and relates to a regional wind-water combined erosion risk assessment method based on coupling coordination degree. Background Technique

[0002] Soil erosion is a natural disaster that occurs widely globally and seriously affects the production and life of human society. Wind and water processes have always been considered important driving forces for land degradation and desertification. Water erosion and wind erosion are two of the most common types of soil erosion. Especially in arid or semi-arid regions, the combined erosion process and effect formed by the interaction of the two erosions are different from single wind erosion and water erosion. Wind-water combined erosion refers to the combined action and alternating action of wind and water on the same object (area), resulting in erosion, transportation, and deposition processes different from those of a single agent. Approximately 17.5% of the global land area is affected by wind-water combined erosion. The land degradation caused by severe wind-water combined erosion has caused global environmental problems and is one of the current research hotspots.

[0003] Wind-water combined erosion mainly occurs in semi-arid regions, coasts, riverbanks, and lakeshores. In the past few decades, rich research results have been obtained on wind-water combined erosion, mainly involving the observation of wind-water combined erosion in experimental plots, the analysis of wind-water combined erosion characteristics at different spatio-temporal scales, and the sediment transportation, accumulation, redistribution, and environmental impacts during the wind-water combined erosion process. Although many current research progresses have been made, these results mainly focus on the erosion sediment monitoring and characteristic analysis at the plot scale, and the research on the erosion process in the wind-water combined erosion area from the perspective of wind erosion or water erosion, and then combining the two for analysis, lacking the theory and method for coupling assessment of wind-water combined erosion force at the regional scale. Therefore, it is very important to propose a regional wind-water combined erosion risk assessment method based on coupling coordination degree. Summary of the Invention

[0004] The purpose of the present invention is to provide a regional wind-water combined erosion risk assessment method based on coupling coordination degree, and using this method can quickly and accurately identify the spatial distribution and change trend of regional wind-water combined erosion risk.

[0005] The technical solution adopted by the present invention is a regional wind-water combined erosion risk assessment method based on coupling coordination degree, which is specifically implemented according to the following steps:

[0006] Step 1, determine the research area, and collect the shp file of the research area boundary, the daily meteorological data of the meteorological station, and the monthly NDVI raster data.

[0007] Step 2: Using the daily meteorological data of the meteorological stations in the study area collected in Step 1, calculate the monthly wind erosion climate erosivity of each meteorological station, and obtain the monthly wind erosion climate erosivity raster data of the study area through interpolation and masking extraction;

[0008] Step 3: Using the daily meteorological data of the meteorological stations in the study area collected in Step 1, calculate the monthly rainfall erosivity of each meteorological station, and obtain the monthly rainfall erosivity raster data of the study area through interpolation and masking extraction;

[0009] Step 4: Based on the coupling coordination degree theory, using the monthly NDVI data of the study area collected in Step 1, the monthly wind erosion climate erosivity raster data of the study area calculated in Step 2, and the monthly rainfall erosivity raster data of the study area calculated in Step 3, calculate the three-dimensional coupling coordination degree of wind erosion climate erosivity, rainfall erosivity, and NDVI;

[0010] Step 5: Based on the calculation results of Step 4, according to the probability statistics theory, select several probability distribution functions as backup options, use the linear moment method to estimate the parameters of several probability distribution functions respectively, and then calculate the corresponding theoretical cumulative probabilities using each backup probability distribution function; at the same time, calculate the empirical cumulative probability, and finally select the optimal theoretical cumulative probability of the three-dimensional coupling coordination degree of wind erosion climate erosivity, rainfall erosivity, and NDVI according to the principle of the smallest RMSE between the empirical cumulative probability and the theoretical cumulative probability;

[0011] Step 6: On the basis of the optimal theoretical cumulative probability selected in Step 5, take the average value of the annual theoretical cumulative probabilities to obtain the spatial distribution results of the average wind-water combined erosion risk in the study area over several years;

[0012] Step 7: On the basis of the optimal theoretical cumulative probability selected in Step 5, according to the Mann-Kendall trend test method, the spatial distribution results of the change trend of the wind-water combined erosion risk in the study area can be identified.

[0013] The features of the present invention also lie in:

[0014] The specific process of Step 2 is as follows:

[0015] Step 2.1: According to the meteorological data collected in Step 1, use the FAO Penman-Monteith method to calculate the monthly potential evapotranspiration of each meteorological station. The specific calculation formula of the FAO Penman-Monteith method is as follows:

[0016]

[0017] In the formula: ET0 is the potential evapotranspiration, mm / d; Δ is the slope of the saturated water vapor pressure curve, kPa·°C -1 ; Rn is the net radiation, MJ·m -2 ·d -1 ; G is the soil heat flux density, MJ·m -2 ·d -1 ; γ is the psychrometric constant, kPa·℃ -1 ; T is the average temperature, ℃; e s is the saturation vapor pressure, kPa; e a is the actual vapor pressure, kPa; u2 is the average wind speed at 2m, m / s;

[0018] The calculation formulas for G, γ, u2 and Δ are as follows:

[0019] G = 0.14(T i -T i-1 ) (2)

[0020]

[0021]

[0022]

[0023] In the formula: T i is the average temperature of the i-th month, ℃; T i-1 is the average temperature of the (i - 1)-th month, ℃; ε takes the value of 0.662; cP takes the value of 1.013×10 -5 , MJ·kg -1 ·℃ -1 ; P e is the atmospheric pressure; u z is the wind speed at 10m, m / s; z is the height, taking the value of 10m;

[0024] R n is the difference between the incoming shortwave radiation R ns and the outgoing net longwave radiation R nl , that is:

[0025] R n = R ns - R nl (6)

[0026] R ns = (1 + a)R s (7)

[0027]

[0028] In the formula: a is the albedo, taking the value of 0.23; R S is the received solar radiation, MJ·m -2 ·d -1; σ is the Stefan-Boltzmann constant, with a value of 4.903×10 -9 , MJ·K -4 ·m -2 ·d -1 ; R S0 is the clear sky radiation, MJ·m -2 ·d -1 ; T max,K is the highest temperature of the day, °C; T min,K is the lowest temperature of the day, °C;

[0029] Among them, the calculation formulas of R S and R S0 are as follows:

[0030]

[0031] R S0 = (a s + b s )R a (10)

[0032]

[0033] In the formula: n is the actual sunshine hours, h; N is the maximum possible sunshine hours, h; n / N is the relative sunshine; a s takes a value of 0.25; b s takes a value of 0.5; R a is the extraterrestrial radiation, MJ·m -2 ·d -1 ;

[0034] Among them, the calculation formula of R a is as follows:

[0035]

[0036]

[0037]

[0038]

[0039] In the formula: G sc is the solar constant, with a value of 0.082; d r is the inverse of the mean sun-earth distance; ω s is the hour angle at sunrise, rad; is the latitude, rad; δ is the longitude, rad; J is the day sequence;

[0040]

[0041]

[0042]

[0043]

[0044] In the formula: e 0 min and e 0 max are the instantaneous saturated water vapor pressures corresponding to the minimum air temperature and the maximum air temperature respectively, kPa; T min is the daily minimum air temperature, °C; T max is the daily maximum air temperature, °C; RH is the average relative humidity;

[0045] Step 2.2: According to the monthly potential evapotranspiration results of each meteorological station calculated in Step 2.1, calculate the monthly wind erosion climate erosivity of each meteorological station. The specific formula is as follows:

[0046]

[0047] In the formula: C is the wind erosion climate erosivity, is the monthly average wind speed at 2 m, m / s; ET 0i is the potential evapotranspiration in the i-th month, mm; P i is the precipitation in the i-th month, mm; d i is the number of days in the i-th month, d;

[0048] Step 2.3: Use ArcGIS software to perform Kriging interpolation on the monthly wind erosion climate erosivity of meteorological stations to obtain the raster data of monthly wind erosion climate erosivity, and use the Yellow River Basin boundary shp file to perform clipping in ArcGIS software to obtain the raster file of monthly wind erosion climate erosivity in the Yellow River Basin.

[0049] The specific process of Step 3 is as follows:

[0050] Step 3.1: First, according to the daily meteorological data collected in Step 1, use Zhang's formula to calculate the rainfall erosivity of each meteorological station monthly. The specific formula is as follows:

[0051]

[0052] In the formula: R i is the rainfall erosivity value during the i-th half-month period, MJ·mm·ha -1 ·h -1 ; k is the number of days during the half-month period; P j is the erosive daily rainfall on the j-th day during the half-month period, mm. When the daily rainfall is greater than 12 mm, it is called erosive daily rainfall. Each month from the 1st to the 15th is a half-month, and the remaining time is a half-month;

[0053] α and β are two important parameters in the model and can be calculated by the following formulas:

[0054] β = 0.8363 + 18.144 / P (d12) + 24.455 / P (y12) (22)

[0055] α = 21.586β -7.1897 (23)

[0056] Step 3.2: According to the rainfall erosivity values of each meteorological station every half month calculated in Step 3.1, the monthly rainfall erosivity value can be obtained by summation;

[0057] Step 3.3: Use ArcGIS software to perform Kriging interpolation on the monthly rainfall erosivity of meteorological stations to obtain the monthly rainfall erosivity raster data, and use the Yellow River Basin boundary shp file to perform clipping in ArcGIS software to obtain the monthly Yellow River Basin rainfall erosivity raster file.

[0058] Step 4 is specifically as follows:

[0059] Step 4.1: According to the monthly NDVI raster data of the Yellow River Basin collected in Step 1, the monthly wind erosion climate erosivity raster data of the Yellow River Basin calculated in Step 2, and the monthly rainfall erosivity data of the Yellow River Basin obtained in Step 3, construct an NDVI subsystem, a wind erosion climate erosivity subsystem, and a rainfall erosivity subsystem respectively. The NDVI subsystem is composed of the monthly NDVI indicators of 12 months of each year, the wind erosion climate erosivity subsystem is composed of the monthly wind erosion climate erosivity indicators of 12 months of each year, and the rainfall erosivity subsystem is composed of the monthly rainfall erosivity indicators of 12 months of each year. Each of the above indicators needs to be normalized. The specific method for normalizing the indicators is as follows:

[0060] When u ij is a value where the larger the value, the better for the system, positive normalization:

[0061]

[0062] When u ij is a value where the smaller the value, the better for the system, negative normalization:

[0063]

[0064] In the formula, u ij is the normalized value of the i-th indicator of the subsystem in the j-th year; x ijis the value of the j-th year of the i-th index of the subsystem; since the wind erosion climate erosivity and rainfall erosivity have a positive contribution to erosion, positive normalization is adopted; NDVI has a negative contribution to erosion, so negative normalization is adopted;

[0065] Calculate the contribution degree of each index of each subsystem in each year of the research period:

[0066]

[0067] In the formula, G ij represents the total contribution degree of the index data x ij to the j-th index; m is the number of samples;

[0068] Calculate the information entropy:

[0069]

[0070] In the formula, e ij represents the information entropy of the j-th index; G ij represents the total contribution degree of the index data x ij to the j-th index; m is the number of samples; where K is a constant related to the number of samples, K = 1 / ln m;

[0071] Calculate the weight:

[0072]

[0073] In the formula, w j represents the weight of the j-th index; n is the number of indexes, and in the constructed coupling system of NDVI, wind erosion climate erosivity, and rainfall erosivity, n is 12; e ij represents the information entropy of the j-th index;

[0074] The specific method for calculating the contributions of the NDVI subsystem, wind erosion climate erosivity subsystem, and rainfall erosivity subsystem to the order degree of the total system is as follows:

[0075]

[0076]

[0077] In the formula, u i is the contribution of the i-th subsystem to the order degree of the total system; u ij is the normalized value of the j-th index in the i-th subsystem; w ij is the weight of the j-th index in the i-th subsystem;

[0078] Step 4.2: According to the coupling degree model of the multivariate system, obtain the coupling degree models among the NDVI subsystem, the wind erosion climate erosivity subsystem, and the rainfall erosivity subsystem in the Yellow River Basin, as shown in the following formulas (31) - (33):

[0079]

[0080] D = (C·T) 1 / 2 (32);

[0081] T = au1 + bu2 + cu3 (33);

[0082] In the formula, D is the coupling coordination degree of the three-dimensional system of the NDVI subsystem, the wind erosion climate erosivity subsystem, and the rainfall erosivity subsystem; C is the coupling degree; u1, u2, and u3 are the contributions of the NDVI subsystem, the wind erosion climate erosivity subsystem, and the rainfall erosivity subsystem to the order degree of the total system, respectively; T is the comprehensive harmonic index of the NDVI subsystem, the wind erosion climate erosivity subsystem, and the rainfall erosivity subsystem, which reflects the overall synergy effect or contribution of the NDVI subsystem, the wind erosion climate erosivity subsystem, and the rainfall erosivity subsystem; a, b, and c are undetermined coefficients.

[0083] The specific process of Step 5 is as follows:

[0084] Step 5.1: According to the coupling coordination degree D of the three-dimensional system of the NDVI subsystem, the wind erosion climate erosivity subsystem, and the rainfall erosivity subsystem calculated in Step 4, use the following formula to calculate the empirical cumulative probability of the coupling coordination degree:

[0085]

[0086] where P 经验 is the empirical cumulative probability of each coupling coordination degree D; n is the number of samples, and i is the order of the coupling coordination degree D arranged from large to small;

[0087] Step 5.2: Use the linear moment method to estimate the parameters of several probability distribution functions respectively, specifically:

[0088] D is the coupling coordination degree, the distribution function is F(D), and there exists an inverse function G(F). G(F) is also called the quantile function. For the coupling coordination degrees D1,..., D n , denote the order statistics:

[0089] D 1:n ≤ D 2:n , …, ≤ D n:n , the r-th order linear moment estimator l r is as follows:

[0090]

[0091] Wherein:

[0092] k ≥ 1, r is the order of the linear moment, and n is the number of samples;

[0093] Step 5.3: Using four alternative theoretical probability distribution functions, namely the exponential distribution, the Gumbel distribution, the Pearson type III distribution, and the generalized extreme value distribution, calculate the corresponding theoretical cumulative probabilities respectively:

[0094] P 理论 = F(D; θ) (36);

[0095] In the formula, P 理论 is the cumulative probability of the theoretical distribution function; F(D; θ) is the theoretical probability distribution function; θ is the parameter of the theoretical distribution function; D is the coupling coordination degree;

[0096] Step 5.4: Calculate the root mean square error (RMSE) of the empirical cumulative probability P 经验 and the theoretical cumulative probability P 理论 of each alternative theoretical distribution function respectively to select the optimal theoretical distribution function. The one with the smallest RMSE calculated among the four theoretical distribution functions is the optimal theoretical cumulative probability P 最优理论 . Using the optimal theoretical cumulative probability P 最优理论 can reflect the risk level of the wind-water combined erosion force. The calculation formula of RMSE is as follows:

[0097]

[0098] In the formula, RMSE is the root mean square error; n is the number of samples; P 经验 is the empirical cumulative probability value of the coupling coordination degree; P 理论 is the theoretical cumulative probability value of the coupling coordination degree.

[0099] The specific implementation in Step 7 is as follows:

[0100] Use the Mann-Kendall trend test method to calculate the trend of the annual optimal theoretical cumulative probability P 最优理论 of the coupling coordination degree of each grid point in the study area calculated in Step 5. A positive Z value calculated by the Mann-Kendall trend test method indicates an increasing trend, a negative value indicates a decreasing trend, and an absolute value of the Z value greater than 1.96 indicates that this change is significant at the 0.05 significance level, thus obtaining the spatial distribution of the change trend of the wind-water combined erosion risk in the Yellow River Basin. The specific calculation formula of the Mann-Kendall trend test method is as follows:

[0101]

[0102]

[0103]

[0104] In the formula: S is the Mann-Kendall trend test statistic; Var(S) is the variance of S; Z is the Mann-Kendall trend test statistic after normal standardization; P 最优理论 is the optimal theoretical cumulative probability of the coupling coordination degree.

[0105] The beneficial effect of the present invention is that the regional wind-water combined erosion risk assessment method based on the coupling coordination degree of the present invention has reliable calculation results, simple calculation methods, clear theoretical significance, can fully reflect the spatial distribution of the regional wind-water combined erosion risk size, and identify areas where the risk may increase and decrease, providing scientific support for the formulation and implementation of regional soil and water conservation planning. Description of the Drawings

[0106] Figure 1 is the implementation technical flow chart in the regional wind-water combined erosion risk assessment method based on the coupling coordination degree of the present invention;

[0107] Figure 2 is a schematic diagram of the spatial distribution result of the average annual coupling coordination degree in the Yellow River Basin in the regional wind-water combined erosion risk assessment method based on the coupling coordination degree of the present invention;

[0108] Figure 3 is a schematic diagram of the spatial distribution result of the average annual risk in the Yellow River Basin in the regional wind-water combined erosion risk assessment method based on the coupling coordination degree of the present invention;

[0109] Figure 4 is a schematic diagram of the spatial distribution result of the average annual risk level in the Yellow River Basin in the regional wind-water combined erosion risk assessment method based on the coupling coordination degree of the present invention;

[0110] Figure 5 is a schematic diagram of the risk change trend result in the Yellow River Basin in the regional wind-water combined erosion risk assessment method based on the coupling coordination degree of the present invention. Detailed Embodiments

[0111] The present invention will be described in detail below in conjunction with the drawings and specific embodiments.

[0112] The regional wind-water combined erosion risk assessment method based on the coupling coordination degree of the present invention specifically includes the following steps:

[0113] Step 1, determine the study area, and collect the boundary shp file of the study area, the daily meteorological data of the meteorological station, and the monthly NDVI raster data;

[0114] Step 2: Using the daily meteorological data of the meteorological stations in the study area collected in Step 1, calculate the monthly wind erosion climate erosivity of each meteorological station, and obtain the monthly wind erosion climate erosivity raster data of the study area through interpolation and masking extraction. The specific process of Step 2 is as follows:

[0115] Step 2.1: First, according to the daily meteorological data of the study area collected in Step 1, use the FAO Penman-Monteith method to calculate the monthly potential evapotranspiration of each meteorological station. The specific calculation formula of the FAO Penman-Monteith method is as follows:

[0116]

[0117] In the formula: ET0 is the potential evapotranspiration, mm / d; Δ is the slope of the saturated water vapor pressure curve, kPa·°C -1 ; R n is the net radiation, MJ·m -2 ·d -1 ; G is the soil heat flux density, MJ·m -2 ·d -1 ; γ is the psychrometric constant, kPa·°C -1 ; T is the average temperature, °C; e s is the saturated water vapor pressure, kPa; e a is the actual water vapor pressure, kPa; u2 is the average wind speed at 2m, m / s.

[0118] The calculation formulas of G, γ, u2 and Δ are as follows:

[0119] G = 0.14(T i -T i-1 )

[0120]

[0121]

[0122]

[0123] In the formula: T i is the average temperature of the i-th month, °C; T i-1 is the average temperature of the (i - 1)-th month, °C; ε takes the value of 0.662; cP takes the value of 1.013×10 -5 , MJ·kg -1 ·°C -1 ; P e is the atmospheric pressure; u z is the wind speed at 10m, m / s; z is the height, taking the value of 10m;

[0124] Rn is the shortwave radiation R of income ns and the net longwave radiation R of expenditure nl The difference is, namely:

[0125] R n = R ns - R nl

[0126] R ns = (1 + a)R s

[0127]

[0128] In the formula: a is the albedo, and its value is 0.23; R S is the solar radiation received, MJ·m -2 ·d -1 ; σ is the Stefan-Boltzmann constant, and its value is 4.903×10 -9 , MJ·K -4 ·m -2 ·d -1 ; R S0 is the clear sky radiation, MJ·m -2 ·d -1 ; T max,K is the highest temperature during a day, °C; T min,K is the lowest temperature during a day, °C;

[0129] Among them, the calculation formulas of R S and R S0 are as follows:

[0130]

[0131] R S0 = (a s + b s )R a

[0132]

[0133] In the formula: n is the actual sunshine duration, h; N is the maximum possible sunshine duration, h; n / N is the relative sunshine; a s takes the value of 0.25; b s takes the value of 0.5; R a is the extraterrestrial radiation, MJ·m -2 ·d -1 ;

[0134] Among them, the calculation formula of R a is as follows:

[0135]

[0136]

[0137]

[0138]

[0139] Where: G sc is the solar constant, with a value of 0.082; d r is the inverse of the mean Earth - Sun distance; ω s is the hour angle at sunrise, rad; is the latitude, rad; δ is the longitude, rad; J is the day sequence number;

[0140]

[0141]

[0142]

[0143]

[0144] Where: e 0 min and e 0 max are the instantaneous saturation water vapor pressures corresponding to the minimum temperature and the maximum temperature respectively, kPa; T min is the daily minimum temperature, °C; T max is the daily maximum temperature, °C; RH is the average relative humidity;

[0145] Step 2.2, calculate the monthly wind - erosion climate erosivity of each meteorological station according to the monthly potential evapotranspiration results of each meteorological station calculated in Step 2.1. The specific formula is as follows:

[0146]

[0147] Where: C is the wind - erosion climate erosivity, is the monthly average wind speed at 2m, m / s; ET 0i is the potential evapotranspiration in the i - th month, mm; P i is the precipitation in the i - th month, mm; d i is the number of days in the i - th month, d;

[0148] Use the ArcGIS software to perform Kriging interpolation on the monthly wind - erosion climate erosivity of the meteorological stations to obtain the raster data of the monthly wind - erosion climate erosivity, and use the shp file of the study area boundary to perform clipping in the ArcGIS software to obtain the raster file of the monthly wind - erosion climate erosivity of the study area;

[0149] Step 3: Using the daily meteorological data of the meteorological stations in the study area collected in Step 1, calculate the monthly rainfall erosivity of each meteorological station, and calculate and obtain the monthly rainfall erosivity raster data of the study area through interpolation and masking extraction. The specific steps of Step 3 are as follows:

[0150] Step 3.1: First, according to the daily meteorological data of the study area collected in Step 1, use Zhang's formula to calculate the monthly rainfall erosivity of each meteorological station. The specific formula is as follows:

[0151]

[0152] In the formula: R i is the rainfall erosivity value during the i-th half-month period, MJ·mm·ha -1 ·h -1 ; k is the number of days in the half-month period; P j is the erosive daily rainfall on the j-th day during the half-month period, mm. When the daily rainfall is greater than 12 mm, it is called erosive daily rainfall. Each month from the 1st to the 15th is a half-month, and the remaining time is a half-month.

[0153] α and β are two important parameters in the model and can be calculated through the following formulas:

[0154] β = 0.8363 + 18.144 / P (d12) + 24.455 / P (y12)

[0155] α = 21.586β -7.1897

[0156] Step 3.2: According to the rainfall erosivity values of each meteorological station every half-month calculated in Step 3.1, the monthly rainfall erosivity value can be obtained by summation. Use ArcGIS software to perform Kriging interpolation on the monthly rainfall erosivity of the meteorological stations to obtain the monthly rainfall erosivity raster data, and use the study area boundary shp file to perform clipping in ArcGIS software to obtain the monthly study area rainfall erosivity raster file;

[0157] Step 4: Based on the coupling coordination degree theory, use the monthly NDVI data of the study area collected in Step 1, the monthly wind erosion climate erosivity raster data of the study area calculated in Step 2, and the monthly rainfall erosivity raster data of the study area calculated in Step 3 to calculate the three-dimensional coupling coordination degree of wind erosion climate erosivity, rainfall erosivity, and NDVI;

[0158] The specific steps of Step 4 are as follows: In Step 4.1, first, based on the monthly NDVI raster data of the study area collected in Step 1, the monthly wind erosion climate erosivity raster data of the study area calculated in Step 2, and the monthly rainfall erosivity data of the study area obtained in Step 3, an NDVI subsystem, a wind erosion climate erosivity subsystem, and a rainfall erosivity subsystem are respectively constructed. The NDVI subsystem consists of the monthly NDVI indicators for 12 months of each year, the wind erosion climate erosivity subsystem consists of the monthly wind erosion climate erosivity indicators for 12 months of each year, and the rainfall erosivity subsystem consists of the monthly rainfall erosivity indicators for 12 months of each year. Each of the above indicators needs to be normalized, and the method for normalizing the indicators is as follows:

[0159] When u ij is such that the larger the value, the better for the system (positive normalization): When u ij is such that the smaller the value, the better for the system (negative normalization):

[0160] In the formula, u ij is the normalized value of the j - th year of the i - th indicator of the subsystem; x ij is the value of the j - th year of the i - th indicator of the subsystem. Since wind erosion climate erosivity and rainfall erosivity have a positive contribution to erosion, positive normalization is adopted; NDVI has a negative contribution to erosion, so negative normalization is adopted.

[0161] Calculate the contribution degree of each indicator of each subsystem in each year of the study period:

[0162]

[0163] In the formula, G ij represents the total contribution degree of the indicator data x ij to the j - th indicator; m is the number of samples (the length of the study period in years); Calculate the information entropy:

[0164]

[0165] In the formula, e ij represents the information entropy of the j - th indicator; G ij represents the total contribution degree of the indicator data x ij to the j - th indicator; m is the number of samples (the length of the study period in years); where K is a constant related to the number of samples, K = 1 / ln m;

[0166] Calculate the weight:

[0167] In the formula, w jdenotes the weight of the j-th index; n is the number of indices, and in the constructed coupling system of NDVI, wind erosion climate erosivity, and rainfall erosivity, n is 12; e ij denotes the information entropy of the j-th index;

[0168] The specific method for calculating the contributions of the NDVI subsystem, wind erosion climate erosivity subsystem, and rainfall erosivity subsystem to the order degree of the total system is as follows:

[0169]

[0170] In the formula, u i is the contribution of the i-th subsystem to the order degree of the total system; u ij is the normalized value of the j-th index in the i-th subsystem; w ij is the weight of the j-th index in the i-th subsystem;

[0171] Step 4.2, according to the coupling degree model of the multi-element system, it is easy to obtain the coupling degree model between the NDVI subsystem, wind erosion climate erosivity subsystem, and rainfall erosivity subsystem, which can be expressed in the following form:

[0172] D = (C·T) 1 / 2 ; T = au1 + bu2 + cu3

[0173] In the formula, D is the coupling coordination degree of the three-dimensional system of the NDVI subsystem, wind erosion climate erosivity subsystem, and rainfall erosivity subsystem; C is the coupling degree; u1, u2, and u3 are the contributions of the NDVI subsystem, wind erosion climate erosivity subsystem, and rainfall erosivity subsystem to the order degree of the total system respectively; T is the comprehensive harmony index of the NDVI subsystem, wind erosion climate erosivity subsystem, and rainfall erosivity subsystem, which reflects the overall synergy effect or contribution of the NDVI subsystem, wind erosion climate erosivity subsystem, and rainfall erosivity subsystem; a, b, and c are undetermined coefficients. In practice, it is often considered that the importance of the three subsystems is the same, so a = b = c = 1 / 3.

[0174] Step 5, based on the calculation results of Step 4, according to the probability statistics theory, select several probability distribution functions as alternative options, use the linear moment method to estimate the parameters of several probability distribution functions respectively, and then calculate the corresponding theoretical cumulative probabilities using each alternative probability distribution function; at the same time, calculate the empirical cumulative probability, and finally select the optimal theoretical distribution function and its theoretical cumulative probability of the three-dimensional coupling coordination degree of wind erosion climate erosivity, rainfall erosivity, and NDVI according to the principle of the minimum RMSE between the empirical cumulative probability and the theoretical cumulative probability;

[0175] The specific steps of Step 5 are as follows:

[0176] Among them, P 经验 is the empirical cumulative probability of each coupling coordination degree D; n is the number of samples, i is the order of the coupling coordination degree D arranged from large to small;

[0177] Step 5.2, use the linear moment method to estimate the parameters of several probability distribution functions respectively. Specifically: D is the coupling coordination degree, the distribution function is F(D), and there is an inverse function G(F). G(F) is also called the quantile function. For the coupling coordination degrees D1,..., D n , denote the order statistic: D 1 :n ≤D 2:n , …, ≤D n:n , the r - order linear moment estimator l r is as follows:

[0178]

[0179] Among them: k ≥ 1, r is the order of the linear moment, n is the number of samples;

[0180] Step 5.3, use the four alternative theoretical probability distribution functions, namely the exponential distribution, the Gumbel distribution, the Pearson type III distribution, and the generalized extreme value distribution, to calculate the corresponding theoretical cumulative probabilities respectively:

[0181] P 理论 = F(D; θ) (36);

[0182] In the formula, P 理论 is the cumulative probability of the theoretical distribution function; F(D; θ) is the theoretical probability distribution function; θ is the parameter of the theoretical distribution function; D is the coupling coordination degree;

[0183] Step 5.4, calculate the RMSE of the empirical cumulative probability P 经验 and the theoretical cumulative probability P 理论 of each alternative theoretical distribution function respectively to select the optimal theoretical distribution function. The one with the smallest RMSE calculated among the four theoretical distribution functions is the optimal theoretical cumulative probability P 最优理论 . Using the optimal theoretical cumulative probability P 最优理论 can reflect the risk level of the geomantic - water combined erosion force. The calculation formula of RMSE is as follows:

[0184]

[0185] In the formula, RMSE is the root - mean - square error; n is the number of samples; P 经验 is the empirical cumulative probability value of the coupling coordination degree; P 理论 is the theoretical cumulative probability value of the coupling coordination degree.

[0186] In step 7, the specific implementation is as follows: Use the Mann-Kendall trend test method to calculate the annual optimal theoretical cumulative probability P of the coupling coordination degree for each grid point in the study area obtained in step 5 最优理论 For trend calculation, a positive Z value obtained by the Mann-Kendall trend test method indicates an increasing trend, a negative value indicates a decreasing trend, and an absolute value of the Z value greater than 1.96 indicates that this change is significant at the 0.05 significance level, and thus the spatial distribution of the change trend of the aeolian-water erosion risk in the Yellow River Basin can be obtained. The specific calculation formula of the Mann-Kendall trend test method is as follows:

[0187]

[0188] In the formula: S is the Mann-Kendall trend test statistic; Var(S) is the variance of S; Z is the Mann-Kendall trend test statistic after normal standardization processing; P 最优理论 is the optimal theoretical cumulative probability of the coupling coordination degree.

[0189] Embodiment

[0190] As Figure 1 shown, taking the Yellow River Basin as an example, the technical process is as Figure 1 shown, and the specific implementation is as follows: Step 1, first collect basic data: Collect the daily precipitation, wind speed, maximum temperature, minimum temperature, average relative humidity, and sunshine hours of 289 meteorological stations in the Yellow River Basin and its adjacent areas from 1982 to 2015; the monthly NDVI data from 1982 to 2015; the shp file of the Yellow River Basin boundary;

[0191] Step 2, the implementation is as follows; Step 2.1, first, according to the daily meteorological data of 289 meteorological stations in the Yellow River Basin collected in step 1 from 1982 to 2015, use the FAO Penman-Monteith method to calculate the monthly potential evapotranspiration of each meteorological station. The specific calculation formula of the FAO Penman-Monteith method is as follows:

[0192]

[0193] In the formula: ET0 is the potential evapotranspiration, mm / d; Δ is the slope of the saturated water vapor pressure curve, kPa·℃ -1 ; R n is the net radiation, MJ·m -2 ·d -1 ; G is the soil heat flux density, MJ·m -2 ·d -1 ; γ is the psychrometric constant, kPa·℃ -1 ; T is the average temperature, ℃; es is the saturated water vapor pressure, kPa; e a is the actual water vapor pressure, kPa; u2 is the average wind speed at 2 m, m / s.

[0194] The calculation formulas for G, γ, u2 and Δ are as follows:

[0195] G = 0.14(T i - T i-1 );

[0196]

[0197] In the formula: T i is the average temperature in the i-th month, °C; T i-1 is the average temperature in the (i - 1)-th month, °C; ε takes the value of 0.662; cP takes the value of 1.013×10 -5 , MJ·kg -1 ·°C -1 ; P e is the atmospheric pressure; u z is the wind speed at 10 m, m / s; z is the height, taking the value of 10 m;

[0198] R n is the difference between the incoming shortwave radiation R ns and the outgoing net longwave radiation R nl , that is:

[0199] R n = R ns - R nl ; R ns = (1 + a)R s

[0200]

[0201] In the formula: a is the albedo, taking the value of 0.23; R S is the received solar radiation, MJ·m -2 ·d -1 ; σ is the Stefan - Boltzmann constant, taking the value of 4.903×10 -9 , MJ·K -4 ·m -2 ·d -1 ; R S0 is the clear - sky radiation, MJ·m -2 ·d -1 ; T max,K is the highest temperature in a day, °C; T min,K is the lowest temperature in a day, °C; Among them, the calculation formulas for R S and R S0 are as follows: R S0 = (a s + b s )R a ;

[0202] In the formula: n is the actual sunshine duration, h; N is the maximum possible sunshine duration, h; n / N is the relative sunshine; a s takes the value of 0.25; b s takes the value of 0.5; R a is the extraterrestrial radiation, MJ·m -2 ·d -1 ; Among them, the calculation formula of R a is as follows:

[0203]

[0204]

[0205]

[0206] In the formula: G sc is the solar constant, taking the value of 0.082; d r is the inverse heliocentric mean distance; ω s is the hour angle at sunrise, rad; is the latitude, rad; δ is the longitude, rad; J is the day sequence;

[0207]

[0208]

[0209] In the formula: e 0 min and e 0 max are the instantaneous saturated water vapor pressures corresponding to the minimum temperature and the maximum temperature respectively, kPa; T min is the daily minimum temperature, °C; T max is the daily maximum temperature, °C; RH is the average relative humidity;

[0210] Step 2.2, according to the monthly potential evapotranspiration results of each meteorological station calculated in Step 2.1, calculate the monthly wind erosion climate erosivity of each meteorological station. The specific formula is as follows:

[0211]

[0212] In the formula: C is the wind erosion climate erosivity, is the monthly average wind speed at 2 m, m / s; ET 0i is the potential evapotranspiration in the i-th month, mm; P iis the precipitation in the i-th month, mm; d i is the number of days in the i-th month, d; Kriging interpolation is performed on the monthly wind erosion climate erosivity of meteorological stations using ArcGIS software to obtain the raster data of monthly wind erosion climate erosivity. The raster file of monthly wind erosion climate erosivity in the Yellow River Basin is obtained by clipping in ArcGIS software using the shp file of the Yellow River Basin boundary;

[0213] Step 3 is specifically implemented as follows; Step 3.1, first, based on the daily meteorological data of 289 meteorological stations in the Yellow River Basin from 1982 to 2015 collected in Step 1, the rainfall erosivity of each meteorological station per month is calculated using Zhang's formula. The specific formula is as follows:

[0214] In the formula: R i is the rainfall erosivity value during the i-th half-month period, MJ·mm·ha -1 ·h -1 ; k is the number of days in the half-month period; P j is the erosive daily rainfall on the j-th day during the half-month period, mm. When the daily rainfall is greater than 12 mm, it is called erosive daily rainfall. Each month from the 1st to the 15th is half a month, and the remaining time is half a month. α and β are two important parameters in the model and can be calculated through the following formulas:

[0215] β = 0.8363 + 18.144 / P (d12) + 24.455 / P (y12) ; α = 21.586β -7.1897

[0216] Step 3.2, based on the rainfall erosivity value of each meteorological station per half-month calculated in Step 3.1, the monthly rainfall erosivity value can be obtained by summing. Kriging interpolation is performed on the monthly rainfall erosivity of meteorological stations using ArcGIS software to obtain the raster data of monthly rainfall erosivity. The raster file of monthly rainfall erosivity in the Yellow River Basin is obtained by clipping in ArcGIS software using the shp file of the Yellow River Basin boundary;

[0217] Step 4 is specifically implemented as follows: In Step 4.1, first, based on the monthly NDVI raster data of the Yellow River Basin from 1982 to 2015 collected in Step 1, the monthly wind erosion climate erosivity raster data of the Yellow River Basin calculated in Step 2, and the monthly rainfall erosivity data of the Yellow River Basin obtained in Step 3, an NDVI subsystem, a wind erosion climate erosivity subsystem, and a rainfall erosivity subsystem are respectively constructed. Among them, the NDVI subsystem is composed of the monthly NDVI indicators for 12 months of each year, the wind erosion climate erosivity subsystem is composed of the monthly wind erosion climate erosivity indicators for 12 months of each year, and the rainfall erosivity subsystem is composed of the monthly rainfall erosivity indicators for 12 months of each year. Each of the above indicators needs to be normalized, and the specific method for normalizing the indicators is as follows:

[0218] When u ij is such that the larger the value, the better for the system (positive normalization): When u ij is such that the smaller the value, the better for the system (negative normalization): In the formula, u ij is the normalized value of the j-th year of the i-th indicator of the subsystem; x ij is the value of the j-th year of the i-th indicator of the subsystem; Since wind erosion climate erosivity and rainfall erosivity have a positive contribution to erosion, positive normalization is adopted; NDVI has a negative contribution to erosion, so negative normalization is adopted. Calculate the contribution degree of each indicator of each subsystem in each year of the research period: In the formula, G ij represents the total contribution degree of the indicator data x ij to the j-th indicator; m is the number of samples; Calculate the information entropy:

[0219] In the formula, e ij represents the information entropy of the j-th indicator; G ij represents the total contribution degree of the indicator data x ij to the j-th indicator; m is the number of samples; where K is a constant related to the number of samples, K = 1 / ln m; Calculate the weight: In the formula, w j represents the weight of the j-th indicator; n is the number of indicators, and in the constructed NDVI, wind erosion climate erosivity, and rainfall erosivity coupling system, n is 12; e ij represents the information entropy of the j-th indicator; The specific method for calculating the contribution of the NDVI subsystem, the wind erosion climate erosivity subsystem, and the rainfall erosivity subsystem to the order degree of the total system is as follows: In the formula; u i is the contribution of the i-th subsystem to the order degree of the total system; u ij is the normalized value of the j-th indicator in the i-th subsystem; w ijis the weight of the jth index in the ith subsystem;

[0220] Step 4.2, according to the coupling degree model of the multi - system, it is easy to obtain the coupling degree model among the NDVI subsystem, the wind erosion climate erosivity subsystem and the rainfall erosivity subsystem in the Yellow River Basin, which can be expressed in the following form: D=(C·T) 1 / 2 ; T = au1+bu2+cu3; where D is the coupling coordination degree of the three - dimensional system of the NDVI subsystem, the wind erosion climate erosivity subsystem and the rainfall erosivity subsystem; C is the coupling degree; u1, u2, u3 are the contributions of the NDVI subsystem, the wind erosion climate erosivity subsystem and the rainfall erosivity subsystem to the order degree of the total system respectively; T is the comprehensive harmonic index of the NDVI subsystem, the wind erosion climate erosivity subsystem and the rainfall erosivity subsystem, which reflects the overall synergy effect or contribution of the NDVI subsystem, the wind erosion climate erosivity subsystem and the rainfall erosivity subsystem; a, b, c are undetermined coefficients. In practice, it is often considered that the importance of the three subsystems is the same, so a = b = c = 1 / 3. The spatial distribution of the average value of the coupling coordination degree of the three - dimensional system of the NDVI subsystem, the wind erosion climate erosivity subsystem and the rainfall erosivity subsystem in the Yellow River Basin from 1982 to 2015 is as Figure 2 shown.

[0221] Step 5, the specific implementation is as follows; Step 5.1, according to the coupling coordination degree D of the three - dimensional system of the NDVI subsystem, the wind erosion climate erosivity subsystem and the rainfall erosivity subsystem calculated in Step 4, use the following formula to calculate the empirical cumulative probability of the coupling coordination degree: where, P 经验 is the empirical cumulative probability of each coupling coordination degree D; n is the number of samples, and i is the order of the coupling coordination degree D arranged from large to small;

[0222] Step 5.2, use the linear moment method to estimate the parameters of several probability distribution functions respectively. Specifically: D is the coupling coordination degree, the distribution function is F(D), and there is an inverse function G(F), and G(F) is also called the quantile function. For the coupling coordination degrees D1,..., D n , denote the order statistics: D 1:n ≤D 2:n , …, ≤D n:n , the r - order linear moment estimator l r of Hosking is as follows: where: k >= 1, r is the order of the linear moment, and n is the number of samples;

[0223] Step 5.3: Using four alternative theoretical probability distribution functions, namely the exponential distribution, Gumbel distribution, Pearson type III distribution, and generalized extreme value distribution, calculate the corresponding theoretical cumulative probabilities respectively:

[0224] P 理论 = F(D; θ);

[0225] In the formula, P 理论 is the cumulative probability of the theoretical distribution function; F(D; θ) is the theoretical probability distribution function; θ is the parameter of the theoretical distribution function; D is the coupling coordination degree;

[0226] Step 5.4: Calculate the empirical cumulative probability P 经验 and the theoretical cumulative probability P 理论 of each alternative theoretical distribution function respectively, and select the optimal theoretical distribution function by comparing the magnitudes of the RMSE values. The one with the smallest RMSE calculated among the four theoretical distribution functions is the optimal theoretical cumulative probability P 最优理论 . Using the optimal theoretical cumulative probability P 最优理论 can reflect the risk level of the combined water and wind erosion force. The calculation formula of RMSE is as follows:

[0227] In the formula, RMSE is the root mean square error; n is the number of samples; P 经验 is the empirical cumulative probability value of the coupling coordination degree; P 理论 is the theoretical cumulative probability value of the coupling coordination degree.

[0228] Step 6: Based on the theoretical cumulative probability of the optimal theoretical distribution function selected in Step 5, use the ArcGIS software to average the annual theoretical cumulative probabilities to obtain the spatial distribution results of the average combined water and wind erosion risk in the Yellow River Basin from 1982 to 2015 for 34 years; the spatial distribution of the combined water and wind erosion risk in the Yellow River Basin from 1982 to 2015 is as Figure 3 shown. By using the natural breaks method in the ArcGIS software to divide it into 3 categories, the spatial distribution of the combined water and wind erosion risk level in the Yellow River Basin can be obtained, as Figure 4 shown.

[0229] Step 7: Use the Mann-Kendall trend test method to calculate the trend of the annual optimal theoretical cumulative probability P 最优理论 of the coupling coordination degree at each grid point in the Yellow River Basin calculated in Step 5. A positive Z value calculated by the Mann-Kendall trend test method indicates an increasing trend, a negative value indicates a decreasing trend, and an absolute value of the Z value greater than 1.96 indicates that this change is significant at the 0.05 significance level. Thus, the spatial distribution of the change trend of the combined water and wind erosion risk in the Yellow River Basin can be obtained. The specific calculation formula of the Mann-Kendall trend test method is as follows:

[0230] In the formula: S is the Mann-Kendall trend test statistic; Var(S) is the variance of S; Z is the Mann-Kendall trend test statistic after normal standardization; P 最优理论 is the optimal theoretical cumulative probability of the coupling coordination degree. The spatial distribution of the change trend of the risk of wind-water combined erosion in the Yellow River Basin from 1982 to 2015 is as Figure 5 shown.

Claims

1. A regional water and soil erosion risk assessment method based on coupling coordination degree, characterized in that The implementation is carried out according to the following steps: Step 1: Determine the research area, and collect the shp file of the research area boundary, the daily meteorological data of meteorological stations, and the monthly NDVI raster data. Step 2: Use the daily meteorological data of the meteorological stations in the research area collected in Step 1 to calculate the monthly wind erosion climate erosivity of each meteorological station, and obtain the monthly wind erosion climate erosivity raster data of the research area through interpolation and masking extraction. The specific process of Step 2 is as follows: Step 2.1: According to the meteorological data collected in Step 1, use the FAO Penman-Monteith method to calculate the monthly potential evapotranspiration of each meteorological station. The specific calculation formula of the FAO Penman-Monteith method is as follows: （1） where: ET0 is the potential evapotranspiration, mm / d; Δ is the slope of the saturated water vapor pressure curve, kPa·°C -1 ; R n is the net radiation, MJ·m -2 ·d -1 ; G is the soil heat flux density, MJ·m -2 ·d -1 ; γ is the psychrometric constant, kPa·°C -1 ; T is the average air temperature, °C; e s is the saturated water vapor pressure, kPa; e a is the actual water vapor pressure, kPa; u2 is the average wind speed at 2 m, m / s; The calculation formulas of G, γ, u2, and Δ are as follows: （2） （4） （5） where: T i is the average temperature in the i-th month, °C; T i-1 is the average temperature in the (i - 1)-th month, °C; The value of ε is 0.662; the value of cP is 1.013×10 -5 , MJ·kg -1 ·°C -1 ; P e is the atmospheric pressure; u z is the wind speed at 10 m, m / s; z is the height, with a value of 10 m; R n is the shortwave radiation R of the income ns and the net longwave radiation R of the expenditure nl The difference, that is: （6） （7） （8） In the formula: a is the albedo, with a value of 0.23; R S is the received solar radiation, MJ·m -2 ·d -1 ; σ is the Stefan-Boltzmann constant, with a value of 4.903×10 -9 , MJ·K -4 ·m -2 ·d -1 ; R S0 is the clear-sky radiation, MJ·m -2 ·d -1 ; T max,K is the highest temperature during a day, °C; T min,K is the lowest temperature during a day, °C; Among them, R S and R S0 are calculated as follows: （9） （10） （11） Where: n is the actual sunshine duration, h; N is the maximum sunshine duration, h; n / N is the relative sunshine; a s The value is 0.25; b s The value is 0.5; R a is the extraterrestrial radiation, MJ·m -2 ·d -1 ; Among them, R a is calculated as follows: （12） （13） （14） （15） Where: G sc is the solar constant, with a value of 0.082; d r is the inverse mean Earth-sun distance; ω s is the hour angle at sunrise, in rad; φ is the latitude, in rad; δ is the longitude, in rad; J is the day sequence; （16） （17） （18） （19） where: e 0 min and e 0 max are the instantaneous saturation water vapor pressures corresponding to the lowest and highest air temperatures, respectively, in kPa; T min is the daily minimum air temperature, in °C; T max is the daily maximum air temperature, in °C; RH is the average relative humidity; Step 2.2: According to the monthly potential evapotranspiration results of each meteorological station calculated in Step 2.1, calculate the monthly wind erosion climate erosivity of each meteorological station. The specific formula is as follows: （20） Where: C is the wind erosion climate erosivity, is the monthly average wind speed at 2 m, m / s; ET 0i is the potential evapotranspiration in the i-th month, mm; P i is the precipitation in the i-th month, mm; d i is the number of days in the i-th month, d; Step 2.3: Use ArcGIS software to perform Kriging interpolation on the monthly wind erosion climate erosivity of meteorological stations to obtain the monthly wind erosion climate erosivity raster data, and use the shp file of the Yellow River Basin boundary to perform clipping in ArcGIS software to obtain the monthly wind erosion climate erosivity raster file of the Yellow River Basin. Step 3: Use the daily meteorological data of the meteorological stations in the research area collected in Step 1 to calculate the monthly rainfall erosivity of each meteorological station, and obtain the monthly rainfall erosivity raster data of the research area through interpolation and masking extraction. Step 4: Based on the coupling coordination degree theory, use the monthly NDVI data of the research area collected in Step 1, the monthly wind erosion climate erosivity raster data of the research area calculated in Step 2, and the monthly rainfall erosivity raster data of the research area calculated in Step 3 to calculate the three-dimensional coupling coordination degree of wind erosion climate erosivity, rainfall erosivity, and NDVI. Step 5, based on the calculation results of Step 4, according to the probability statistics theory, select several probability distribution functions as alternative choices, use the linear moment method to estimate the parameters of several probability distribution functions respectively, and then calculate the corresponding theoretical cumulative probabilities using each alternative probability distribution function; at the same time, calculate the empirical cumulative probability, and finally, according to the principle of the smallest RMSE between the empirical cumulative probability P 经验 and the theoretical cumulative probability P 理论 , select the optimal theoretical cumulative probability P 最优理论 of the three-dimensional coupling coordination degree D of wind erosion climate erosivity, rainfall erosivity, and NDVI; Step 6, based on the optimal theoretical cumulative probability P 最优理论 selected in Step 5, take the average value of P 最优理论 for each grid year by year to obtain the spatial distribution results of the average fengshui combined erosion risk in the study area for several years; Step 7, based on the optimal theoretical cumulative probability P 最优理论 selected in Step 5, the spatial distribution results of the risk change trend of the combined wind and water erosion in the study area can be identified according to the Mann-Kendall trend test method.

2. The method for evaluating the risk of regional water and wind erosion based on the coupling coordination degree according to claim 1, wherein The specific process of Step 3 is as follows: Step 3.1: First, according to the daily meteorological data collected in Step 1, use Zhang's formula to calculate the monthly rainfall erosivity of each meteorological station. The specific formula is as follows: （21） Where: R i is the rainfall erosivity value during the i-th half-month period, MJ·mm·ha -1 ·h -1 ; k is the number of days in the half-month period; P j is the erosive daily rainfall on the j-th day during the half-month period, mm. When the daily rainfall is greater than 12 mm, it is called erosive daily rainfall. Each month from the 1st to the 15th is a half-month, and the remaining time is a half-month; α and β are two important parameters in the model and can be calculated by the following formula: （22） （23） Step 3.2: According to the rainfall erosivity values of each meteorological station every half month calculated in Step 3.1, sum them up to obtain the monthly rainfall erosivity value. Step 3.3: Use ArcGIS software to perform Kriging interpolation on the monthly rainfall erosivity of meteorological stations to obtain the monthly rainfall erosivity raster data, and use the shp file of the Yellow River Basin boundary to perform clipping in ArcGIS software to obtain the monthly rainfall erosivity raster file of the Yellow River Basin.

3. The regional water and wind erosion risk assessment method based on coupling coordination degree according to claim 2, wherein Step 4 is specifically as follows: Step 4.1: Construct an NDVI subsystem, a wind erosion climate erosivity subsystem, and a rainfall erosivity subsystem based on the monthly NDVI raster data of the Yellow River Basin collected in Step 1, the monthly wind erosion climate erosivity raster data calculated in Step 2, and the monthly rainfall erosivity data obtained in Step 3. The NDVI subsystem consists of monthly NDVI indicators for 12 months of each year, the wind erosion climate erosivity subsystem consists of monthly wind erosion climate erosivity indicators for 12 months of each year, and the rainfall erosivity subsystem consists of monthly rainfall erosivity indicators for 12 months of each year. Normalize each indicator. The specific method for normalizing the indicators is as follows: When u ij is a value where the larger the value is, the better it is for the system, positive normalization: （24） When u ij is such that the smaller the value, the better for the system, negative normalization: （25） where u ij is the normalized value of the j-th year of the i-th index of the subsystem; x ij is the value of the j-th year of the i-th index of the subsystem; since the wind erosion climate erosivity and rainfall erosivity have a positive contribution to erosion, positive normalization is adopted; NDVI has a negative contribution to erosion, so negative normalization is adopted; Calculate the contribution degree of each indicator of each subsystem in each year during the research period: （26） Where G ij represents the total contribution degree of the index data x ij to the j-th index; m is the number of samples; Calculate the information entropy: （27） where, e ij represents the information entropy of the j-th index; G ij represents the total contribution degree of the index data x ij to the j-th index; m is the number of samples; where K is a constant related to the number of samples, ; Calculate the weight: （28） where w j represents the weight of the j-th index; n is the number of indices, and in the constructed coupled system of NDVI, wind erosion climate erosivity, and rainfall erosivity, n is 12; e ij represents the information entropy of the j-th index; The specific method for calculating the contributions of the NDVI subsystem, the wind erosion climate erosivity subsystem, and the rainfall erosivity subsystem to the order degree of the total system is as follows: （29） （30） where, u i is the contribution of the i-th subsystem to the order degree of the total system; u ij is the normalized value of the j-th index in the i-th subsystem; w ij is the weight of the j-th index in the i-th subsystem; Step 4.2: According to the coupling degree model of the multi - element system, obtain the coupling degree model among the NDVI subsystem, the wind erosion climate erosivity subsystem, and the rainfall erosivity subsystem of the Yellow River Basin, as shown in the following formulas (31) - (33): In the formula, D is the coupling coordination degree of the three - dimensional system of the NDVI subsystem, the wind erosion climate erosivity subsystem, and the rainfall erosivity subsystem; C is the coupling degree; u1, u2, and u3 are the contributions of the NDVI subsystem, the wind erosion climate erosivity subsystem, and the rainfall erosivity subsystem to the order degree of the total system respectively; T is the comprehensive coordination index of the NDVI subsystem, the wind erosion climate erosivity subsystem, and the rainfall erosivity subsystem, which reflects the overall synergy effect or contribution of the NDVI subsystem, the wind erosion climate erosivity subsystem, and the rainfall erosivity subsystem; a, b, and c are undetermined coefficients.

4. The regional geomantic water erosion risk assessment method based on coupling coordination degree according to claim 3, characterized in that The specific process of Step 5 is as follows: Step 5.1: According to the coupling coordination degree D of the three - dimensional system of the NDVI subsystem, the wind erosion climate erosivity subsystem, and the rainfall erosivity subsystem calculated in Step 4, use the following formula to calculate the empirical cumulative probability of the coupling coordination degree: （34） where P 经验 is the empirical cumulative probability of each coupling coordination degree D; n is the number of samples, and i is the order of the coupling coordination degree D arranged from large to small; Step 5.2, use the linear moment method to estimate the parameters of several probability distribution functions respectively, specifically: D is the coupling coordination degree, the distribution function is , and there exists an inverse function , G ( F ) is also called the quantile function. For the coupling coordination degree D 1,..., D n , record the order statistic: , the r - order linear moment estimator of Hosking is as follows: （35）； Wherein: ; ; , k >=1, where r is the order of the linear moment and n is the number of samples; Step 5.3: Use four alternative theoretical probability distribution functions, namely the exponential distribution, the Gumbel distribution, the Pearson type III distribution, and the generalized extreme value distribution, to calculate the corresponding theoretical cumulative probabilities respectively: （36）； where P 理论 is the cumulative probability of the theoretical distribution function; F ( D; θ ) is the theoretical probability distribution function; is the parameter of the theoretical distribution function; D is the coupling coordination degree; Step 5.4, calculate the empirical cumulative probability P 经验 and the theoretical cumulative probability P 理论 of each alternative theoretical distribution function, and select the optimal theoretical distribution function by the magnitude of RMSE. The one with the smallest RMSE calculated among the four theoretical distribution functions is the optimal theoretical distribution function. The theoretical cumulative probability calculated by the optimal theoretical distribution function can reflect the risk level of the water and wind combined erosion force. The calculation formula of RMSE is as follows: Wherein, RMSE is the root mean square error; n is the number of samples; is the empirical cumulative probability value of the coupling coordination degree; is the theoretical cumulative probability value of the coupling coordination degree.

5. The regional geomantic water erosion risk assessment method based on coupling coordination degree according to claim 4, characterized in that The specific implementation in Step 7 is as follows: Using the Mann-Kendall trend test method, the annual optimal theoretical cumulative probability of the coupling coordination degree of each grid point in the study area calculated in step 5 is used for trend calculation. A positive Z value obtained by the Mann-Kendall trend test method indicates an increasing trend, and a negative value indicates a decreasing trend. An absolute value of the Z value greater than 1.96 indicates that this change is significant at the 0.05 significance level, and thus the spatial distribution of the change trend of the risk of combined erosion by wind and water in the study area can be obtained. The specific calculation formula of the Mann-Kendall trend test method is as follows: Where: S is the Mann-Kendall trend test statistic; Var(S) is the variance of S; Z is the Mann-Kendall trend test statistic after normal standardization; is the optimal theoretical cumulative probability of the coupling coordination degree.

Images