A joint prediction method for soil erosion based on Bayesian weights

Abstract

The invention discloses a soil erosion joint prediction method based on Bayesian weights. The invention relates to the technical field of soil erosion prediction. By constructing a multi-source heterogeneous database of a watershed, a soil erosion joint simulation framework is established based on five models, namely, a MUSLE model, a WEPP model, a LISEM model, a EUROSEM model, and a SIMWE model; a parameter sensitivity analysis is performed for each soil erosion model using a FAST algorithm; a parameter calibration and verification is performed for each soil erosion model using a PSO algorithm; a Bayesian sampling is performed for each soil erosion model using a DREAM algorithm to estimate the posterior distribution of parameters; based on the Bayesian sampling results, the uncertainty of the soil erosion model is evaluated, and the Bayesian weights of different soil erosion models are calculated using a BMA algorithm to perform soil erosion joint prediction.

Description

Technical Field

[0001] The invention relates to the technical field of soil erosion prediction, and in particular to a soil erosion joint prediction method based on Bayesian weights. Background Art

[0002] Quantitative analysis of soil erosion is one of the important contents of soil erosion research. There are many methods to study soil erosion. With the increasing development and maturity of remote sensing technology, global positioning system, geographic information system and in-depth research on soil erosion models, soil erosion models play an increasingly important role in the quantitative analysis of soil loss and the formulation of protection strategies. The application of soil erosion models is becoming more and more extensive. Compared with traditional soil erosion research methods, such as field experimental observations, isotope tracing, runoff plot simulation experiments, etc., soil erosion models are not only suitable for small-scale research, but also play an important role in large-scale research that is difficult to reach with traditional methods.

[0003] Soil erosion models can be divided into empirical statistical models and physical cause models according to modeling techniques. The empirical statistical model uses statistical methods and experimental data as research methods to determine the influencing factors of soil erosion, thereby constructing a soil loss calculation formula. The physical process model is based on the physical process of soil erosion, focuses on physical causes, and quantifies the physical process of soil erosion through interactions within the system. The empirical model has low requirements for input data, a relatively simple structure, and is easy to calculate. The physical process model requires more sophisticated spatial data support and a large amount of data calculation. In terms of prediction accuracy, the physical process model can often simulate the soil erosion conditions in the study area more accurately and comprehensively.

[0004] Since the soil erosion process involves the interaction of multiple factors such as climate, topography, soil, and vegetation, the model often has to simplify these complexities, resulting in large parameter uncertainties and structural uncertainties in the soil erosion model. It is necessary to verify and evaluate the internal processes and simulation results of the existing models with the same data set and boundary conditions, quantitatively evaluate the applicability and accuracy of different models, and reduce the uncertainty of erosion models and improve the accuracy of erosion forecasts through scientific and effective methods. Summary of the invention

[0005] The present invention provides a soil erosion joint prediction method based on Bayesian weights to solve the problem of high uncertainty in soil erosion assessment commonly used in the prior art.

[0006] The present invention provides a soil erosion joint prediction method based on Bayesian weights, comprising the following steps:

[0007] Step 1: construct a watershed multi-source heterogeneous database, wherein the data of the watershed multi-source heterogeneous database includes hydrological data, meteorological data, vegetation data, terrain data and soil data; preprocess the data, wherein the preprocessing includes unification of spatial resolution, unification of spatial boundaries and unification of temporal resolution;

[0008] Step 2: According to the multi-source heterogeneous database of the watershed, a soil erosion joint simulation framework is established based on the MUSLE model, WEPP model, LISEM model, EUROSEM model and SIMWE model to carry out soil erosion forecasting;

[0009] Step 3: Use the FAST algorithm to conduct parameter sensitivity analysis for each soil erosion model and determine the sensitive parameter set;

[0010] Step 4: Use the PSO algorithm to calibrate and verify parameters for each soil erosion model, and obtain the optimal parameter set with the best simulation effect as the goal;

[0011] Step 5: Use the DREAM algorithm to perform Bayesian sampling for each soil erosion model and estimate the posterior distribution of parameters;

[0012] Step 6: Based on the Bayesian sampling results, the uncertainty of the soil erosion model is evaluated, and the Bayesian weights of different soil erosion models are calculated using the BMA algorithm to perform joint soil erosion forecasting.

[0013] Further, step 1 comprises:

[0014] Step 1.1, construct a watershed multi-source heterogeneous database, wherein the data of the watershed multi-source heterogeneous database includes hydrological data, meteorological data, vegetation data, terrain data and soil data;

[0015] Step 1.2: If the data time scale is not uniform, data expansion, interpolation, and reduction are used to match the data to ensure the consistency of the data time series;

[0016] Step 1.3: If the spatial scale of the data is not uniform, the inverse distance weighted interpolation method, Kriging interpolation method, and resampling method are used for matching to ensure the consistency of the spatial resolution and spatial range of the raster data.

[0017] Further, step 2 includes:

[0018] Step 2.1, the MUSLE model equation is:

[0019]

[0020] Where: is the sediment yield; is the surface runoff depth; is the peak flow rate; is the area of ​​the hydrological response unit; is the soil erodibility factor; is the terrain factor; are vegetation cover and crop management factors; is the soil and water conservation measures factor; is the coarse debris factor;

[0021] The simulation of sediment transport in the basin is calculated by the runoff sediment entrainment:

[0022]

[0023] Where: It is the sediment-carrying force of runoff; is the roughness factor; is the cumulative value of single-width flow within the time step; is the gradient of the slope or river channel; and are model parameters, which are exponents of flow and slope factors respectively; the runoff sediment entrainment used in the process of sediment transport on the slope is calculated based on the surface runoff, and the runoff sediment entrainment used in the process of sediment transport in the river channel is calculated based on the total runoff;

[0024] Step 2.2: The sediment transport in the WEPP model is controlled by the continuity equation. Soil erosion is divided into inter-rill erosion and rill erosion, and the sediment transport process is calculated based on the magnitude of the runoff sediment carrying capacity. The continuity equation is:

[0025]

[0026] Where: is the sediment load per unit width per unit time; is the distance along the downhill direction; is the rate at which inter-rill eroded sediment is transported to the rills; is the rate of rill erosion or deposition;

[0027] Among them, inter-rill erosion is calculated based on the effective rainfall intensity, rill erosion is calculated based on the shear stress of surface runoff, and the runoff sediment carrying capacity of different sediment particle compositions is obtained using the modified Yalin formula:

[0028]

[0029] Where: is the correction factor for inter-rill erosion that ignores certain factors; is the interrill erodibility coefficient; is the vegetation factor; is the ground cover factor; is the effective rainfall intensity; is the spacing between fine grooves; is the width of the rill;

[0030]

[0031] Where: is a correction factor for rill erosion that ignores certain flow conditions; rill shear stress correction factor to reduce uncertainty in drag coefficient calculations; is the rill erodibility coefficient; is the shear stress of water flow acting on soil particles; is the critical shear stress for soil stripping; It is the sediment-carrying force of rill flow runoff;

[0032] Step 2.3, EUROSEM is a process-based single event model that can simulate the rainfall erosion process in a small watershed based on the physical process of soil erosion, including vegetation interception, rainfall splash erosion, runoff erosion, sediment transport and deposition processes; the rainfall splash erosion is calculated based on the rainfall kinetic energy:

[0033]

[0034] Where: is the rainfall splash erosion; is the soil particle dispersion index; is the density of sediment; is the total kinetic energy of rainfall, including throughfall and canopy intercepted rainfall; is a correction index, ranging from 0.9 to 3.1, depending on soil texture; is the average depth of thin-layer water flow on the surface;

[0035] Runoff scour, sediment transport and deposition are calculated based on unit runoff power:

[0036]

[0037] Where: It is the sediment-carrying force of runoff; is the unit runoff power; is the critical unit runoff power; and The median particle size of sediment Related empirical coefficients;

[0038]

[0039] Where: The runoff flushing capacity; is the runoff scouring efficiency coefficient, which is between 0 and 1 and is related to soil cohesion; For flow width; is the soil particle settling velocity; is the sediment content of runoff;

[0040] Step 2.4, LISEM is a distributed physical model that takes into account the spatial differences of the underlying surface and can simulate the soil erosion and sediment transport process of a single rainfall event, including rainfall splash erosion, runoff erosion, sediment transport and deposition processes; splash erosion is calculated based on the rainfall kinetic energy:

[0041]

[0042] Where: is the sputtering separation amount; is the total aggregate stability of soil; Kinetic energy for rainfall; is the depth of the surface water layer; the amount of rainfall intercepted by the forest canopy; is the total surface area where splash erosion occurs;

[0043] Runoff erosion, sediment transport and deposition are calculated based on the unit runoff power using the EUROSEM model:

[0044]

[0045] Where: It is the sediment-carrying force of runoff; is the density of sediment; is the unit runoff power; is the critical unit runoff power; and The median particle size of sediment Related empirical coefficients;

[0046]

[0047] Where: The runoff flushing capacity; is the efficiency coefficient, between 0 and 1, and is related to soil cohesion; For flow width; is the soil particle settling velocity; is the sediment content of runoff;

[0048] Step 2.5: The SIMWE model simulates the surface runoff and sediment transport process based on the two-dimensional Saint-Venant equations, and simulates the sediment deposition and transport through the sediment continuity equation:

[0049]

[0050] Where: is the net erosion and deposition rate; is the density of sediment; is the sediment content of runoff; is the runoff depth; For time; is the sediment transport rate per single width.

[0051] Further, step 3 includes:

[0052] Step 3.1: Use the FAST algorithm to determine the sensitivity of the parameters to be calibrated; by analyzing the variance of the model output, decompose the overall variance into the contribution of each input parameter, and explain the relative impact of different input parameters on the output results;

[0053] Step 3.2, assign a corresponding characteristic frequency to each input parameter, and design and select the frequency set based on the number and physical characteristics of the input parameters;

[0054] Step 3.3, using the conversion function and characteristic frequency, sampling in the parameter space, and obtaining the corresponding objective function value result based on the model calculation;

[0055] Step 3.4: Perform sensitivity analysis based on the sampling results. Determine the relative impact of each input parameter value on the output result according to the spectral amplitude corresponding to different characteristic frequencies. Explain the contribution of the model input parameters to the output variance. Use the first-order sensitivity index as an evaluation indicator. It is between 0 and 1. The larger the value, the more sensitive the parameter.

[0056] Further, step 4 includes:

[0057] Step 4.1: The sensitive parameter set is calibrated using the PSO algorithm. A group of random particles are initialized in the parameter space for searching. Each particle searches for the optimal solution as an individual extreme value. , the best individual extreme value in the entire particle swarm is the global optimal solution ;

[0058] Step 4.2: Initialize a group of particles in the parameter space, initialize the two attributes of velocity and position, and initialize the individual extreme value of each particle and the global optimal value of the entire particle swarm ;

[0059] Step 4.3, calculate the objective function value according to the current position of the particle;

[0060] Step 4.4: Compare the value of the objective function and update the individual extreme value and the global optimal solution ;

[0061] Step 4.5: According to individual extreme values and the global optimal solution Update the particle's velocity according to the formula and location :

[0062]

[0063]

[0064] Where: and are the learning factors of individuals and groups respectively; and are two random factors between 0 and 1 respectively;

[0065] Step 4.6: Determine whether the end condition is met. If not, return to step 4.3 to continue iterating. Otherwise, output the global optimal solution. , and end the optimization;

[0066] Step 4.7. Select the Kling-Gupta efficiency coefficient and RMS error Two indicators are used as the objective function of the PSO algorithm to compare the simulated sediment transport distribution mean with the measured value:

[0067]

[0068]

[0069] Where: is the Pearson correlation coefficient; is the relative rate of change, which is the ratio of the standard deviation of the simulated value to the standard deviation of the observed value; is the bias ratio, which is the ratio of the sum of simulated values ​​to the sum of observed values; For the sequence observations; For the sequence indivual; is the sequence length.

[0070] It can indicate the overall effect of the model, and its value ranges from negative infinity to 1. The closer it is to 1, the higher the match between the simulated value of the model and the observed value. It indicates the error between the simulated value and the measured value. The value ranges from 0 to positive infinity. The closer it is to 0, the smaller the model simulation error.

[0071] Further, step 5 includes:

[0072] Step 5.1: The DREAM algorithm is used to sample the model parameters and the posterior distribution of the results. The differential evolution algorithm is used to generate the candidate parameter set. The Metropolis-Hastings acceptance criterion is used to determine whether to accept the candidate parameter set. The Gelman-Rubin statistic is used as the judgment condition for convergence, and multiple Markov chains are used to simultaneously explore the posterior distribution of the parameters in the sampling space.

[0073] Step 5.2: Determine the parameters to be estimated based on the characteristics of the soil erosion model:

[0074]

[0075]

[0076] Where: is the simulation result of the erosion model, is the parameter set of the erosion model, is the simulation process function of the model; Represents a parameter set Each parameter in is subject to uniform distribution of their upper and lower bounds;

[0077] Step 5.3, create multiple Markov chains, initialize parameter values ​​and likelihood function values, where the likelihood function uses the logarithmic transformation form of the Gaussian likelihood:

[0078]

[0079] Where: Represents the likelihood function value, which characterizes the possibility of the observed data under the current parameter value, and its range is The larger the value, the greater the possibility, and the more the current simulation results match the observed values; is the length of the data sequence; A number for each data point; is the measurement error of the observed data, which defaults to the observed data point , if it appears The value is from Pick a random number between and issue a warning; is the observation data; For the current parameter set Simulated data under the value;

[0080] Step 5.4, perform correction and mutation operations in the parameter space through the differential evolution algorithm to generate a candidate parameter set;

[0081] Step 5.5, calculate the likelihood function value of the simulation result of the candidate parameter set, and determine whether to accept the candidate parameter set as the new state according to the Metropolis-Hastings acceptance criterion;

[0082] Step 5.6: Adaptively adjust the step size of the sampling update to control the mutation strategy according to the performance of the Markov chain;

[0083] Step 5.7: Check whether the Gelman-Rubin statistic has reached the convergence standard, or the number of runs exceeds the maximum allowed number. If it has not converged, return to step 5.4, where the Gelman-Rubin statistic The calculation is as follows:

[0084]

[0085]

[0086]

[0087]

[0088] Where: is the number of Markov chains; is the number of repetitions; is the total variance within and between chains, For the chain The mean of the variances of the parameters; For chain The variance of the parameters; Represents the first parameters; Representative The mean of the parameters on all chains; Represents the mean of all parameters on all chains;

[0089] Step 5.8: After reaching the convergence criterion, perform a fixed number of extended samplings.

[0090] Further, step 6 comprises:

[0091] Step 6.1: Calculate the weights of each soil erosion model according to the general information criterion to avoid the problem of overfitting. The smaller the value, the better the comprehensive effect of soil erosion. For each data point , The calculation is as follows:

[0092]

[0093]

[0094]

[0095] Where: for The scaling parameters of is the logarithmic point prediction density, which is used to measure how well the model fits the data; is the log-likelihood variance of the posterior sample, which is used to measure the complexity of the model; For the The simulated value is calculated by extracting a set of parameter combinations from the posterior probability of the parameters for each observation; For data points The likelihood value of ;

[0096] Step 6.2: Under the normal distribution assumption, the log-likelihood value matrix is ​​calculated as follows:

[0097]

[0098] Where: is the observed value of the corresponding data point, is the standard deviation of the observation sequence, In the parameter set Model simulation values ​​under conditions;

[0099] Step 6.3: Use WAIC value to replace marginal likelihood function for Bayesian model averaging. The weight calculation formula is:

[0100]

[0101] Where: For the The weight of each model; For the The model is relative to the best model Relative difference; the denominator is the normalization factor.

[0102] The present invention has the following beneficial effects: a soil erosion joint forecasting method based on Bayesian weights of the present invention has sufficient theoretical basis, reliable practical application, and fully considers the deficiencies of existing technologies and methods. It can provide a convenient and feasible joint forecasting approach for soil erosion simulation, and provide a core algorithm for soil erosion forecasting in smart water conservancy, and has good application and promotion value. BRIEF DESCRIPTION OF THE DRAWINGS

[0103] In order to more clearly illustrate the technical solution of the present invention, the drawings required for use in the embodiments are briefly introduced below. Obviously, for ordinary technicians in this field, other drawings can be obtained based on these drawings without creative work.

[0104] Figure 1 The present invention provides a flow chart of a soil erosion joint prediction method based on Bayesian weights according to an embodiment of the present invention. DETAILED DESCRIPTION

[0105] In order to make the purpose, technical solution and advantages of the present invention clearer, the technical solution of the present invention will be clearly and completely described in conjunction with the specific embodiments of the present invention and the corresponding drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention. The technical solutions provided by the embodiments of the present invention are described in detail below in conjunction with the drawings.

[0106] See also Figure 1 The present invention provides a soil erosion joint prediction method based on Bayesian weights, comprising the following steps:

[0107] Step 1: construct a multi-source heterogeneous database of the watershed, wherein the data of the multi-source heterogeneous database of the watershed includes hydrological data, meteorological data, vegetation data, terrain data and soil data; preprocess the data, wherein the preprocessing includes unification of spatial resolution, unification of spatial boundaries and unification of temporal resolution.

[0108] Specifically, step 1 includes:

[0109] Step 1.1, construct a watershed multi-source heterogeneous database, wherein the data of the watershed multi-source heterogeneous database includes hydrological data, meteorological data, vegetation data, terrain data and soil data;

[0110] Step 1.2: If the data time scale is not uniform, data expansion, interpolation, and reduction are used to match the data to ensure the consistency of the data time series;

[0111] Step 1.3: If the spatial scale of the data is not uniform, the inverse distance weighted interpolation method, Kriging interpolation method, and resampling method are used for matching to ensure the consistency of the spatial resolution and spatial range of the raster data.

[0112] Step 2: According to the multi-source heterogeneous database of the watershed, a soil erosion joint simulation framework is established based on five models, namely, the MUSLE model, the WEPP model, the LISEM model, the EUROSEM model, and the SIMWE model, for soil erosion forecasting.

[0113] Specifically, step 2 includes:

[0114] Step 2.1, MUSLE (Modified Universal Soil Loss Equation) model is an empirical soil erosion model improved on the basis of USLE (Universal Soil Loss Equation), and the model equation is:

[0115]

[0116] Where: is the sediment yield, in units of ; is the surface runoff depth, in ; is the peak flow rate, in units of ; is the hydrological response unit area, in units of ; is the soil erodibility factor, in units of ; is the terrain factor; are vegetation cover and crop management factors; is the soil and water conservation measures factor; is the coarse debris factor;

[0117] The simulation of sediment transport in the basin is calculated by the runoff sediment entrainment:

[0118]

[0119] Where: is the runoff sediment carrying capacity, in units of ; is the roughness factor; is the cumulative value of the single-width flow within the time step, in units of ; is the gradient of the slope or river channel; and are model parameters, and are exponents of flow and slope factors, respectively. The runoff sediment entrainment used in the process of sediment transport on the slope is calculated based on the surface runoff, and the runoff sediment entrainment used in the process of sediment transport in the river channel is calculated based on the total runoff.

[0120] Step 2.2: The sediment transport in the WEPP (Water Erosion Prediction Project) model is controlled by the continuity equation. Soil erosion is divided into inter-rill erosion and rill erosion, and the sediment transport process is calculated based on the magnitude of the runoff sediment carrying capacity. The continuity equation is:

[0121]

[0122] Where: is the sediment load per unit time and unit width, in units of ; is the distance along the downhill direction, in units of ; is the rate at which eroded sediment between rills is transported to the rills, in ; is the rill erosion or deposition rate, in ;

[0123] Among them, inter-rill erosion is calculated based on the effective rainfall intensity, rill erosion is calculated based on the shear stress of surface runoff, and the runoff sediment carrying capacity of different sediment particle compositions is obtained using the modified Yalin formula:

[0124]

[0125] Where: Some factors are neglected, such as the correction factor for inter-rill erosion in rock cover; is the inter-rill erodibility coefficient, in units of ; is the vegetation factor; is the ground cover factor; is the effective rainfall intensity, in units of ; is the spacing between fine grooves, in ; is the width of the rill, in ;

[0126]

[0127] Where: is a correction factor for rill erosion that ignores certain flow conditions; rill shear stress correction factor to reduce uncertainty in drag coefficient calculations; is the rill erodibility coefficient, in units of ; is the shear stress of water flow acting on soil particles, in units of ; is the critical shear stress of soil peeling, in units of ; is the sediment-carrying capacity of rill flow, in units of ;

[0128] Step 2.3, EUROSEM (The European Soil Erosion Model) is a process-based single event model that can simulate the rainfall erosion process in a small watershed based on the physical process of soil erosion, including vegetation interception, rainfall splash erosion, runoff erosion, sediment transport and deposition processes; the rainfall splash erosion is calculated based on the rainfall kinetic energy:

[0129]

[0130] Where: is the rainfall splash erosion, in units of ; is the soil particle dispersion index, in units of ; is the sediment density, in ; is the total kinetic energy of rainfall, in units of , including throughfall and canopy intercepted rainfall; It is a correction index, ranging from 0.9 to 3.1, depending on soil texture, and is usually taken as 2.0; is the average depth of thin-layer water flow on the surface, in units of ;

[0131] Runoff scour, sediment transport and deposition are calculated based on unit runoff power:

[0132]

[0133] Where: is the runoff sediment carrying capacity, in units of ; is the unit runoff power, unit ; is the critical unit runoff power, in units of ; and The median particle size of sediment Related empirical coefficients;

[0134]

[0135] Where: is the runoff scouring capacity, in units of ; is the runoff scouring efficiency coefficient, which is between 0 and 1 and is related to soil cohesion; is the flow width, in units of ; is the soil particle settling velocity, in units of ; is the sediment content of runoff, in units of ;

[0136] Step 2.4, LISEM (LImburg Soil Erosion Model) is a distributed physical model that takes into account the spatial differences of the underlying surface. It can simulate the soil erosion and sediment transport process of a single rainfall event, mainly including rainfall splash erosion, runoff erosion, sediment transport and deposition processes; splash erosion is calculated based on the rainfall kinetic energy:

[0137]

[0138] Where: is the sputtering separation amount, in units of ; is the total aggregate stability of soil; is the rainfall kinetic energy, in units of ; is the depth of the surface water layer, in ; is the amount of rainfall intercepted by the forest canopy, in ; is the total surface area where splash erosion occurs, in units of ;

[0139] Runoff erosion, sediment transport and deposition are calculated based on the unit runoff power using the EUROSEM model:

[0140]

[0141] Where: is the runoff sediment carrying capacity, in units of ; is the sediment density, in ; is the unit runoff power, in units of ; is the critical unit runoff power, in units of ; and The median particle size of sediment Related empirical coefficients;

[0142]

[0143] Where: is the runoff scouring capacity, in units of ; is the efficiency coefficient, between 0 and 1, and is related to soil cohesion; is the flow width, in units of ; is the soil particle settling velocity, in units of ; is the sediment content of runoff, in units of ;

[0144] Step 2.5, SIMWE (Simulation of Water Erosion Model) simulates the surface runoff and sediment transport process based on the two-dimensional Saint-Venant equations, and simulates the sediment deposition and transport through the sediment continuity equation:

[0145]

[0146] Where: is the net erosion and deposition rate in ; is the sediment density, in ; is the sediment content of runoff, in units of ; is the runoff depth in units of ; is the time, the unit is ; is the sediment transport rate per unit width, in units of .

[0147] Step 3: Use the FAST (Fourier Amplitude Sensitivity Test) algorithm to conduct parameter sensitivity analysis for each soil erosion model and determine the sensitive parameter set.

[0148] Specifically, step 3 includes:

[0149] Step 3.1: Use the FAST algorithm to determine the sensitivity of the parameters to be calibrated; by analyzing the variance of the model output, decompose the overall variance into the contribution of each input parameter, and explain the relative impact of different input parameters on the output results;

[0150] Step 3.2, assign a corresponding characteristic frequency to each input parameter, and design and select the frequency set based on the number and physical characteristics of the input parameters;

[0151] Step 3.3, using the conversion function and characteristic frequency, sampling in the parameter space, and obtaining the corresponding objective function value result based on the model calculation;

[0152] Step 3.4: Perform sensitivity analysis based on the sampling results. Determine the relative impact of each input parameter value on the output result according to the spectral amplitude corresponding to different characteristic frequencies. Explain the contribution of the model input parameters to the output variance. Use the first-order sensitivity index as an evaluation indicator. It is between 0 and 1. The larger the value, the more sensitive the parameter.

[0153] Step 4: Use the PSO algorithm to calibrate and verify parameters for each soil erosion model, and obtain the optimal parameter set with the best simulation effect as the goal.

[0154] Specifically, step 4 includes:

[0155] Step 4.1: The sensitive parameter set is calibrated using the PSO (Particle Swarm Optimization) algorithm. A group of random particles are initialized in the parameter space for searching. Each particle searches for the optimal solution as an individual extreme value. , the best individual extreme value in the entire particle swarm is the global optimal solution ;

[0156] Step 4.2: Initialize a group of particles in the parameter space, initialize the two attributes of velocity and position, and initialize the individual extreme value of each particle and the global optimal value of the entire particle swarm ;

[0157] Step 4.3, calculate the objective function value according to the current position of the particle;

[0158] Step 4.4: Compare the value of the objective function and update the individual extreme value and the global optimal solution ;

[0159] Step 4.5: According to individual extreme values and the global optimal solution Update the particle's velocity according to the formula and location :

[0160]

[0161]

[0162] Where: and As learning factors for individuals and groups respectively; and are two random factors between 0 and 1 respectively;

[0163] Step 4.6: Determine whether the end condition is met. If not, return to step 4.3 to continue iterating. Otherwise, output the global optimal solution. , and end the optimization;

[0164] Step 4.7, Select (Kling-Gupta Efficiency, Kling-Gupta Efficiency) and The two indicators (Root Mean Squared Error) are used as the objective function of the PSO algorithm to compare the simulated sediment transport distribution mean with the measured value:

[0165]

[0166]

[0167] Where: is the Pearson correlation coefficient; is the relative rate of change, which is the ratio of the standard deviation of the simulated value to the standard deviation of the observed value; is the bias ratio, which is the ratio of the sum of simulated values ​​to the sum of observed values; For the sequence observations; For the sequence indivual; is the sequence length.

[0168] It can indicate the overall effect of the model, and its value ranges from negative infinity to 1. The closer it is to 1, the higher the match between the simulated value of the model and the observed value. It indicates the error between the simulated value and the measured value. The value ranges from 0 to positive infinity. The closer it is to 0, the smaller the model simulation error.

[0169] Step 5: Use the DREAM algorithm to perform Bayesian sampling for each soil erosion model and estimate the posterior distribution of parameters.

[0170] Specifically, step 5 includes:

[0171] Step 5.1: The DREAM (DiffeRential Evolution Adaptive Metropolis) algorithm is used to sample the model parameters and the posterior distribution of the results. The Differential Evolution algorithm is used to generate the candidate parameter set. The Metropolis-Hastings acceptance criterion is used to determine whether to accept the candidate parameter set. The Gelman-Rubin statistic is used as the convergence judgment condition, and multiple Markov chains are used to simultaneously explore the posterior distribution of the parameters in the sampling space.

[0172] Step 5.2: Determine the parameters to be estimated based on the characteristics of the soil erosion model:

[0173]

[0174]

[0175] Where: is the simulation result of the erosion model, is the parameter set of the erosion model, is the simulation process function of the model; Represents a parameter set Each parameter in is subject to uniform distribution of their upper and lower bounds;

[0176] Step 5.3, create multiple Markov chains, initialize parameter values ​​and likelihood function values, where the likelihood function uses the logarithmic transformation form of the Gaussian likelihood:

[0177]

[0178] Where: Represents the likelihood function value, which characterizes the possibility of the observed data under the current parameter value, and its range is The larger the value, the greater the possibility, and the more the current simulation results match the observed values; is the length of the data sequence; A number for each data point; is the measurement error of the observed data, which defaults to the observed data point , if it appears The value is from Pick a random number between and issue a warning; is the observation data; For the current parameter set Simulated data under the value;

[0179] Step 5.4, perform correction and mutation operations in the parameter space through the differential evolution algorithm to generate a candidate parameter set;

[0180] Step 5.5, calculate the likelihood function value of the simulation result of the candidate parameter set, and determine whether to accept the candidate parameter set as the new state according to the Metropolis-Hastings acceptance criterion;

[0181] Step 5.6: Adaptively adjust the step size of the sampling update to control the mutation strategy according to the performance of the Markov chain;

[0182] Step 5.7: Check whether the Gelman-Rubin statistic has reached the convergence standard, or the number of runs exceeds the maximum allowed number. If it has not converged, return to step 5.4, where the Gelman-Rubin statistic The calculation is as follows:

[0183]

[0184]

[0185]

[0186]

[0187] Where: is the number of Markov chains; is the number of repetitions; is the total variance within and between chains, For the chain The mean of the variances of the parameters; For chain The variance of the parameters; Represents the first parameters; Representative The mean of the parameters on all chains; Represents the mean of all parameters on all chains;

[0188] Step 5.8: After reaching the convergence criterion, perform a fixed number of extended samplings.

[0189] The above fixed number of times defaults to 100 times.

[0190] Step 6: Based on the Bayesian sampling results, the uncertainty of the soil erosion model is evaluated, and the Bayesian weights of different soil erosion models are calculated using the BMA (Bayesian model average) algorithm to perform joint soil erosion forecasting.

[0191] Specifically, step 6 includes:

[0192] Step 6.1: Calculate the weights of each soil erosion model according to the general information criterion to avoid the problem of overfitting. The smaller the value, the better the comprehensive effect of soil erosion. For each data point , The calculation is as follows:

[0193]

[0194]

[0195]

[0196] Where: for The scaling parameters of is the logarithmic point prediction density, which is used to measure how well the model fits the data; is the log-likelihood variance of the posterior sample, which is used to measure the complexity of the model; For the The simulated value is calculated by extracting a set of parameter combinations from the posterior probability of the parameters for each observation; For data points The likelihood value of ;

[0197] Step 6.2: Under the normal distribution assumption, the log-likelihood value matrix is ​​calculated as follows:

[0198]

[0199] Where: is the observed value of the corresponding data point, is the standard deviation of the observation sequence, In the parameter set Model simulation values ​​under conditions;

[0200] Step 6.3: Use WAIC value to replace marginal likelihood function for Bayesian model averaging. The weight calculation formula is:

[0201]

[0202] Where: For the The weight of each model; For the The model is relative to the best model Relative difference; the denominator is the normalization factor.

[0203] The erosion multi-model joint prediction method of the present invention is described below with reference to specific cases.

[0204] Step 1: construct a watershed multi-source heterogeneous database, wherein the data of the watershed multi-source heterogeneous database includes hydrological data, meteorological data, vegetation data, terrain data and soil data; preprocess the data, wherein the preprocessing includes unification of spatial resolution, unification of spatial boundaries and unification of temporal resolution;

[0205] Taking a typical watershed in a certain area as an example, basic data of the watershed was collected and sorted, and ArcGIS Pro was used to build a multi-source heterogeneous database of the watershed, including data on vegetation cover, soil properties, terrain slope, land use type, meteorological elements, etc. The land use data uses high-resolution remote sensing images, through supervised classification, and field verification. The land use types include six types: forest, grassland, bare land or unused land, cultivated land, water body, and urban land.

[0206] Specifically, ArcGIS Pro software was used to preprocess vegetation data, terrain data, soil data, etc., and the inverse distance weighted method, Kriging interpolation method, and resampling method were used to unify the range of spatial data and the resolution of raster data; data expansion, interpolation, and reduction methods were used to match the time series of precipitation, runoff, and sediment content data.

[0207] Step 2: According to the multi-source heterogeneous database of the watershed, a soil erosion model joint simulation framework is established based on five models, namely, the MUSLE model, the WEPP model, the LISEM model, the EUROSEM model, and the SIMWE model, for soil erosion forecasting.

[0208] Step 3: Use the FAST algorithm to conduct parameter sensitivity analysis for each soil erosion model and determine the sensitive parameter set.

[0209] Specifically, the sensitive parameter of the MUSLE model is the runoff factor index and slope factor index ; The sensitive parameter of the WEPP model is the correction factor for rill shear stress ; The EUROSEM model and the LISEM model are similar, and the sensitive parameters are both Manning coefficients and sediment median particle size ; The sensitive parameter of the SIMWE model is the shear stress coefficient In general, the parameters of each soil erosion model have similar sensitivity characteristics, that is, the parameters of sediment transport process are the most sensitive, followed by the parameters of runoff erosion process, while the parameters of rainfall splash process are relatively less sensitive.

[0210] Step 4: Use the PSO algorithm to calibrate and verify parameters for each soil erosion model, and obtain the optimal parameter set with the best simulation effect as the goal;

[0211] Specifically, the simulation results of the model flood events are shown in Table 1. The simulation results of the WEPP model are the best, with a regular average rate of is 0.620, average The validation period index is 1.897, which is also the best. The LISEM model and the MUSLE model are second, with good overall performance. The EUROSEM model and the SIMWE model perform relatively poorly.

[0212] Table 1 Erosion model simulation results

[0213] Step 5: Use the DREAM algorithm to perform Bayesian sampling for each soil erosion model and estimate the posterior distribution of parameters.

[0214] Specifically, Bayesian modeling was performed on each soil erosion model, and the DREAM algorithm was used to establish five Markov chains to sample four typical processes with regular rates, and the posterior distribution of parameters of the MUSLE model, WEPP model, EUROSEM model, LISEM model, and SIMWE model was obtained. Compared with the results of sensitivity analysis, it can be found that the standard deviation and coefficient of variation of sensitive parameters are relatively small, most of them can quickly converge to a small interval, and the optimal parameter value is slightly different from the mean of the posterior distribution, and the model parameter values ​​are relatively stable; while the standard deviation and coefficient of variation of insensitive parameters are large, and the uncertainty is large.

[0215] Step 6: Based on the Bayesian sampling results, the uncertainty of the soil erosion model is evaluated, and the Bayesian weight of the soil erosion model is calculated using the BMA algorithm to perform joint soil erosion forecasting.

[0216] Specifically, based on the Bayesian sampling results obtained in step 5, the BMA algorithm is used to calculate the weights of each soil erosion model, and the results are: the weight of the MUSLE model is 0.5%, the weight of the WEPP model is 59.8%, the weight of the EUROSEM model is 0.1%, the weight of the LISEM model is 33.1%, and the weight of the SIMWE model is 6.5%. Taking the typical secondary flood event of 1988-07-15 in the basin as an example, the simulation results of the posterior distribution of sediment transport of each soil erosion model are shown, including the simulated value, 95% interval and 50% interval. According to the posterior distribution of the simulation results of the MUSLE model, WEPP model, EUROSEM model, LISEM model, and SIMWE model and the posterior distribution of the simulation results of the Bayesian model average algorithm, it can be concluded that the simulated values ​​of the typical event sediment production and transport process are highly consistent with the observed values, which is mainly reflected in the fitting of the peak time, peak shape and peak value. In addition, the 50% intervals of each model are narrow and close to the mean of the posterior distribution, indicating that the simulation results are relatively accurate and the uncertainty is acceptable. However, the 95% intervals of each model vary greatly. The 95% interval of the posterior distribution of the simulation results calculated by the Bayesian model average (BMA) method is the narrowest and has the smallest uncertainty. The 95% intervals of the WEPP model and the LISEM model are also narrow, with small uncertainties. The 95% interval of the SIMWE model is also narrow, but the upper bound of the interval is high. The 95% intervals of the MUSLE and EUROSEM models are relatively wide, with large uncertainties, which also explains why they are assigned the lowest weights.

[0217] As can be seen from the above, the present invention realizes the joint prediction of multiple soil erosion models based on GIS. Specifically, a multi-source heterogeneous database of the watershed is constructed, and the data of the multi-source heterogeneous database of the watershed includes hydrological, meteorological, vegetation, topographic, soil and other data, and the data is preprocessed, including the unification of spatial resolution, spatial boundary and time resolution; secondly, a soil erosion model is established, and a simulation program under a unified simulation framework is established based on five models, namely, the MUSLE model, the WEPP model, the LISEM model, the EUROSEM model and the SIMWE model, for carrying out soil erosion; finally, the FAST algorithm is used to carry out parameter sensitivity analysis for each soil erosion prediction model, and the PSO algorithm is used to calibrate and verify the physical parameters of each soil erosion prediction model, and the optimal parameter set is obtained with the best simulation effect as the goal, and the DREAM algorithm is used to sample the Bayesian model constructed for different soil erosion prediction models, and the joint simulation result distribution is constructed according to the posterior distribution of the model, and the uncertainty of each model is determined, and the BMA algorithm is used to assign weights, and different soil erosion prediction models are jointly simulated. The present invention performs joint soil erosion forecasting based on Bayesian weights, fully considers the erosion and sediment transport process mechanisms at different scales, can accurately simulate the erosion and sediment transport in large and medium-scale watersheds, and has good promotion and practicability.

[0218] The specific implementation methods and specific parameters in the embodiments described above do not constitute a limitation on the protection scope of the present invention. The patent protection scope of the present invention shall still be based on its claims. Any equivalent structural changes made using the contents of the description and drawings of the present invention should also be included in the protection scope of the present invention.

Claims

1. A soil erosion joint prediction method based on Bayesian weights, characterized in that: The following steps are involved: Step 1: construct a watershed multi-source heterogeneous database, wherein the data of the watershed multi-source heterogeneous database includes hydrological data, meteorological data, vegetation data, terrain data and soil data; preprocess the data, wherein the preprocessing includes unification of spatial resolution, unification of spatial boundaries and unification of temporal resolution; Step 2: According to the multi-source heterogeneous database of the watershed, a soil erosion model joint simulation framework is established based on the MUSLE model, WEPP model, LISEM model, EUROSEM model and SIMWE model to carry out soil erosion forecasting; Step 3: Use the FAST algorithm to conduct parameter sensitivity analysis for each soil erosion model and determine the sensitive parameter set; Step 4: Use the PSO algorithm to calibrate and verify parameters for each soil erosion model, and obtain the optimal parameter set with the best simulation effect as the goal; Step 5: Use the DREAM algorithm to perform Bayesian sampling for each soil erosion model and estimate the posterior distribution of parameters; Step 5 includes: Step 5.1: The DREAM algorithm is used to sample the model parameters and the posterior distribution of the results. The differential evolution algorithm is used to generate the candidate parameter set. The Metropolis-Hastings acceptance criterion is used to determine whether to accept the candidate parameter set. The Gelman-Rubin statistic is used as the judgment condition for convergence, and multiple Markov chains are used to simultaneously explore the posterior distribution of the parameters in the sampling space. Step 5.2: Determine the parameters to be estimated based on the characteristics of the soil erosion model: Where: is the simulation result of the erosion model, is the parameter set of the erosion model, is the simulation process function of the erosion model; Represents a parameter set Each parameter in is subject to uniform distribution of their upper and lower bounds; Step 5.3, create multiple Markov chains, initialize parameter values ​​and likelihood function values, where the likelihood function uses the logarithmic transformation form of the Gaussian likelihood: Where: Represents the likelihood function value, which characterizes the possibility of the observed data under the current parameter value, and its range is The larger the value, the greater the possibility, and the more the current simulation results match the observed values; is the length of the data sequence; A number for each data point; is the measurement error of the observed data, which defaults to the observed data point , if it appears The value is from Pick a random number between and issue a warning; is the observation data; For the current parameter set Simulated data under the value; Step 5.4, perform correction and mutation operations in the parameter space through the differential evolution algorithm to generate a candidate parameter set; Step 5.5, calculate the likelihood function value of the simulation result of the candidate parameter set, and determine whether to accept the candidate parameter set as the new state according to the Metropolis-Hastings acceptance criterion; Step 5.6: Adaptively adjust the step size of the sampling update to control the mutation strategy according to the performance of the Markov chain; Step 5.7: Check whether the Gelman-Rubin statistic has reached the convergence standard, or the number of runs exceeds the maximum allowed number. If it has not converged, return to step 5.4, where the Gelman-Rubin statistic The calculation is as follows: Where: is the number of Markov chains; is the number of repetitions; is the total variance within and between chains, For the chain The mean of the variances of the parameters; For chain The variance of the parameters; Represents the first parameters; Representative The mean of the parameters on all chains; Represents the mean of all parameters on all chains; Step 5.8: After reaching the convergence standard, perform a fixed number of extended samplings; Step 6: Based on the Bayesian sampling results, the uncertainty of the soil erosion model is evaluated, and the Bayesian weights of different soil erosion models are calculated using the BMA algorithm to perform joint soil erosion forecasting.

2. A soil erosion joint prediction method based on Bayesian weights as claimed in claim 1, characterized in that: Step 1 includes: Step 1.1, construct a watershed multi-source heterogeneous database, wherein the data of the watershed multi-source heterogeneous database includes hydrological data, meteorological data, vegetation data, terrain data and soil data; Step 1.2: If the data time scale is not uniform, data expansion, interpolation, and reduction are used to match the data to ensure the consistency of the data time series; Step 1.3: If the spatial scale of the data is not uniform, the inverse distance weighted interpolation method, Kriging interpolation method, and resampling method are used for matching to ensure the consistency of the spatial resolution and spatial range of the raster data.

3. A soil erosion joint prediction method based on Bayesian weights as claimed in claim 1, characterized in that: Step 2 includes: Step 2.1, the MUSLE model equation is: Where: is the sediment yield; is the surface runoff depth; is the peak flow rate; is the area of ​​the hydrological response unit; is the soil erodibility factor; is the terrain factor; are vegetation cover and crop management factors; is the soil and water conservation measures factor; is the coarse debris factor; The simulation of sediment transport in the basin is calculated by the runoff sediment entrainment: Where: It is the sediment-carrying force of runoff; is the roughness factor; is the cumulative value of single-width flow within the time step; is the gradient of the slope or river channel; and are model parameters, which are exponents of flow and slope factors respectively; the runoff sediment entrainment used in the process of sediment transport on the slope is calculated based on the surface runoff, and the runoff sediment entrainment used in the process of sediment transport in the river channel is calculated based on the total runoff; Step 2.2: The sediment transport in the WEPP model is controlled by the continuity equation. Soil erosion is divided into inter-rill erosion and rill erosion, and the sediment transport process is calculated based on the magnitude of the runoff sediment carrying capacity. The continuity equation is: Where: is the sediment load per unit width per unit time; is the distance along the downhill direction; is the rate at which inter-rill eroded sediment is transported to the rills; is the rate of rill erosion or deposition; Among them, inter-rill erosion is calculated based on the effective rainfall intensity, rill erosion is calculated based on the shear stress of surface runoff, and the runoff sediment carrying capacity of different sediment particle compositions is obtained using the modified Yalin formula: Where: is the correction factor for inter-rill erosion that ignores certain factors; is the interrill erodibility coefficient; is the vegetation factor; is the ground cover factor; is the effective rainfall intensity; is the spacing between fine grooves; is the width of the rill; Where: is a correction factor for rill erosion that ignores certain flow conditions; rill shear stress correction factor to reduce uncertainty in drag coefficient calculations; is the rill erodibility coefficient; is the shear stress of water flow acting on soil particles; is the critical shear stress for soil stripping; It is the sediment-carrying force of rill flow runoff; Step 2.3, EUROSEM is a process-based single event model that can simulate the rainfall erosion process in a small watershed based on the physical process of soil erosion, including vegetation interception, rainfall splash erosion, runoff erosion, sediment transport and deposition processes; the rainfall splash erosion is calculated based on the rainfall kinetic energy: Where: is the rainfall splash erosion; is the soil particle dispersion index; is the density of sediment; is the total kinetic energy of rainfall, including throughfall and canopy intercepted rainfall; is a correction index, ranging from 0.9 to 3.1, depending on soil texture; is the average depth of thin-layer water flow on the surface; Runoff scour, sediment transport and deposition are calculated based on unit runoff power: Where: It is the sediment-carrying force of runoff; is the unit runoff power; is the critical unit runoff power; and The median particle size of sediment Related empirical coefficients; Where: The runoff flushing capacity; is the runoff scouring efficiency coefficient, which is between 0 and 1 and is related to soil cohesion; For flow width; is the soil particle settling velocity; is the sediment content of runoff; Step 2.4, LISEM is a distributed physical model that takes into account the spatial differences of the underlying surface and can simulate the soil erosion and sediment transport process of a single rainfall event, including rainfall splash erosion, runoff erosion, sediment transport and deposition processes; splash erosion is calculated based on the rainfall kinetic energy: Where: is the sputtering separation amount; is the total aggregate stability of soil; Kinetic energy for rainfall; is the depth of the surface water layer; the amount of rainfall intercepted by the forest canopy; is the total surface area where splash erosion occurs; Runoff erosion, sediment transport and deposition are calculated based on the unit runoff power using the EUROSEM model: Where: It is the sediment-carrying force of runoff; is the density of sediment; is the unit runoff power; is the critical unit runoff power; and The median particle size of sediment Related empirical coefficients; Where: The runoff flushing capacity; is the efficiency coefficient, between 0 and 1, and is related to soil cohesion; For flow width; is the soil particle settling velocity; is the sediment content of runoff; Step 2.5: The SIMWE model simulates the surface runoff and sediment transport process based on the two-dimensional Saint-Venant equations, and simulates the sediment deposition and transport through the sediment continuity equation: Where: is the net erosion and deposition rate; is the density of sediment; is the sediment content of runoff; is the runoff depth; For time; is the sediment transport rate per single width.

4. A soil erosion joint prediction method based on Bayesian weights as claimed in claim 1, characterized in that: Step 3 includes: Step 3.1: Use the FAST algorithm to determine the sensitivity of the parameters to be calibrated; by analyzing the variance of the model output, decompose the overall variance into the contribution of each input parameter, and explain the relative impact of different input parameters on the output results; Step 3.2, assign a corresponding characteristic frequency to each input parameter, and design and select the frequency set based on the number and physical characteristics of the input parameters; Step 3.3, using the conversion function and characteristic frequency, sampling in the parameter space, and obtaining the corresponding objective function value result based on the model calculation; Step 3.4: Perform sensitivity analysis based on the sampling results. Determine the relative impact of each input parameter value on the output result according to the spectral amplitude corresponding to different characteristic frequencies. Explain the contribution of the model input parameters to the output variance. Use the first-order sensitivity index as an evaluation indicator. It is between 0 and 1. The larger the value, the more sensitive the parameter.

5. The soil erosion joint prediction method based on Bayesian weights according to claim 1, characterized in that: Step 4 includes: Step 4.1: The sensitive parameter set is calibrated using the PSO algorithm. A group of random particles are initialized in the parameter space for searching. Each particle searches for the optimal solution as an individual extreme value. , the best individual extreme value in the entire particle swarm is the global optimal solution ; Step 4.2: Initialize a group of particles in the parameter space, initialize the two attributes of velocity and position, and initialize the individual extreme value of each particle and the global optimal value of the entire particle swarm ; Step 4.3, calculate the objective function value according to the current position of the particle; Step 4.4: Compare the value of the objective function and update the individual extreme value and the global optimal solution ; Step 4.5: According to individual extreme values and the global optimal solution Update the particle's velocity according to the formula and location : Where: and are the learning factors of individuals and groups respectively; and are two random factors between 0 and 1 respectively; Step 4.6: Determine whether the end condition is met. If not, return to step 4.3 to continue iterating. Otherwise, output the global optimal solution. , and end the optimization; Step 4.

7. Select the Kling-Gupta efficiency coefficient and RMS error Two indicators are used as the objective function of the PSO algorithm to compare the simulated sediment transport distribution mean with the measured value: Where: is the Pearson correlation coefficient; is the relative rate of change, which is the ratio of the standard deviation of the simulated value to the standard deviation of the observed value; is the bias ratio, which is the ratio of the sum of simulated values ​​to the sum of observed values; For the sequence observations; For the sequence indivual; is the sequence length; It can indicate the overall effect of the model, and its value ranges from negative infinity to 1. The closer it is to 1, the higher the match between the simulated value of the model and the observed value. It indicates the error between the simulated value and the measured value. The value ranges from 0 to positive infinity. The closer it is to 0, the smaller the model simulation error.

6. A soil erosion joint prediction method based on Bayesian weights as claimed in claim 1, characterized in that: Step 6 includes: Step 6.1: Calculate the weights of each soil erosion model according to the general information criterion to avoid the problem of overfitting. The smaller the value, the better the comprehensive effect of soil erosion. For each data point , The calculation is as follows: Where: for The scaling parameters of is the logarithmic point prediction density, which is used to measure how well the model fits the data; is the log-likelihood variance of the posterior sample, which is used to measure the complexity of the model; For the The simulated value is calculated by extracting a set of parameter combinations from the posterior probability of the parameters for each observation; For data points The likelihood value of ; Step 6.2: Under the normal distribution assumption, the log-likelihood value matrix is ​​calculated as follows: Where: is the observed value of the corresponding data point, is the standard deviation of the observation sequence, In the parameter set Model simulation values ​​under conditions; Step 6.3: Use WAIC value to replace marginal likelihood function for Bayesian model averaging. The weight calculation formula is: Where: For the The weight of each model; For the The model is relative to the best model Relative difference; the denominator is the normalization factor.

Images