Water and soil conservation ecological product value accounting method and system
By splitting the data by slope aspect, calculating the aspect-dependent radiation coupling flux and hydrological connectivity matrix, eliminating topographic bias, and quantifying the cumulative impact of mismatch, this method solves the problem that the aspect difference and the impact of sunlight are not considered in traditional methods, and improves the accuracy and reliability of the value accounting of soil and water conservation ecological products.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- YUNNAN ACAD OF ENVIRONMENTAL SCI
- Filing Date
- 2026-01-15
- Publication Date
- 2026-05-01
AI Technical Summary
Traditional methods for calculating the value of soil and water conservation ecological products fail to effectively consider the impact of slope aspect differences and topographic illumination on vegetation indices, resulting in biased calculation results. They cannot accurately capture the true correlation between illumination, vegetation, and sediment yield, and fail to quantify the cumulative impact of mismatches in multiple upstream and neighboring units, thus reducing the reliability of the calculation results.
By splitting the data by slope aspect, vegetation index, incident light factor and normalized sediment yield intensity sequence are obtained, slope aspect irradiance coupling flux is calculated, hydrological connectivity matrix is established, closed forward propagation is performed, topographic bias is removed, audit weights are mapped, soil retention increment is calculated by combining rain erosion force, and regional weighted summation is performed.
It enables the quantification of slope aspect system deviation and spatial propagation effects, the vegetation index is more in line with the actual cover status, the regional summary results can reflect the actual impact of different units, and provide more realistic data for assessing the value of soil and water conservation ecological products.
Smart Images

Figure CN121961650A_ABST
Abstract
Description
A method and system for calculating the value of soil and water conservation ecological products Technical Field
[0001] This invention relates to the field of ecological value accounting technology, and more specifically, to a method and system for accounting for the value of soil and water conservation ecological products. Background Technology
[0002] The accounting of the ecological product value of soil and water conservation provides key data support for ecological protection decision-making and management, and is of great significance in the field of ecological value assessment. However, most of the current mainstream accounting methods are based on the analysis of regional overall data, which do not fully consider the systematic interference of topographic features on the accounting process, resulting in deviations between the accounting results and the actual situation, making it difficult to meet the needs of precise management.
[0003] Slope aspect is a core topographic factor influencing surface processes, directly determining the amount of solar radiation received by the surface. Sunny and shady slopes inherently differ in the duration and intensity of sunlight exposure, leading to significant variations in vegetation growth rates and soil moisture content, thus altering soil erosion and sediment yield. Traditional methods fail to separate data by slope aspect, mixing vegetation, sunlight, and sediment yield data from sunny and shady slopes for analysis. This masks the systematic bias caused by slope aspect within the overall data, failing to accurately capture the true correlation between sunlight, vegetation, and sediment yield. Within a watershed, slope units form an organic whole through hydrological connectivity, and the effects of mismatch propagate and accumulate along confluence paths. Traditional methods do not consider this spatial propagation characteristic, calculating the accounting indicators of individual slope units in isolation, ignoring the cumulative effects of mismatch from upstream and neighboring units, resulting in inaccurate quantification of overall regional mismatch risk.
[0004] Furthermore, traditional methods fail to effectively eliminate the linear bias of topographic illumination on vegetation indices, making it difficult for vegetation indices to reflect the true cover status. Simultaneously, the weighting design does not incorporate differences in mismatch risk and does not adaptively adjust the contributions of units prone to mismatch, resulting in regional aggregation results that fail to reflect the actual impact of different units. This further reduces the reliability of the calculation results and makes it difficult to support accurate assessment and management decisions regarding the value of soil and water conservation ecological products. Summary of the Invention
[0005] This invention provides a method and system for calculating the value of soil and water conservation ecological products, solving the technical problems mentioned in the background.
[0006] This invention provides a method for calculating the value of soil and water conservation ecological products, comprising the following steps: Step S101, within the pre-flood season window, the slope unit is divided into two parts according to slope aspect, and the vegetation index sequence, incident light factor sequence, and normalized sediment yield intensity sequence are obtained. The unit area and rain erosion force are obtained, and the difference between the vegetation index and the difference between the incident light factor are calculated respectively, and the difference between the normalized sediment yield intensity is calculated; Step S102, the local regression sensitivity difference between the difference between the vegetation index and the difference between the incident light factor is calculated in the low coverage and high coverage intervals respectively, and the sensitivity difference is weighted by combining the difference between the normalized sediment yield intensity to generate the aspect irradiance coupling flux; Step S103, a hydrological connectivity matrix is established, and the aspect irradiance coupling flux is calculated according to the symbol... Step S104: Using the aspect ratio irradiance coupling flux as the initial value of the node, closed forward propagation is performed through the propagation graph weight matrix and the convergence coefficient based on the spectral radius to calculate the mismatch potential energy. Step S105: The topographic bias is removed by removing the local slope of the incident light factor sequence from the vegetation index sequence to obtain the debiased vegetation index. The mismatch potential energy is mapped to the audit weight based on the regional mismatch potential energy mean. Step S106: The interannual variation of the debiased vegetation index is calculated. The soil retention increment is calculated by combining the normalized sediment yield intensity sequence and the rain erosion force. The interannual variation and the soil retention increment are regionally weighted and summarized through the audit weight to obtain the accounting result.
[0007] This invention provides a system for calculating the value of soil and water conservation ecological products, comprising: a data processing module, which divides slope units into two parts according to slope aspect within the pre-flood season window, obtains vegetation index sequences, incident light factor sequences, and normalized sediment yield intensity sequences, obtains unit area and rain erosion capacity, calculates the difference in vegetation indices and incident light factor respectively, and calculates the difference in normalized sediment yield intensity; a coupling flux generation module, which calculates the local regression sensitivity difference of vegetation index differences to incident light factor differences in low and high coverage intervals respectively, and weights the sensitivity difference with the normalized sediment yield intensity difference to generate slope aspect irradiance coupling flux; and a propagation graph weight matrix calculation module, which establishes a hydrological connectivity matrix and calculates the propagation graph weight matrix according to the sign of the slope aspect irradiance coupling flux. The system employs several modules: a directional screening and edge assignment module, a normalized edge weight matrix for propagation graph weighting, a mismatch potential energy calculation module (using aspect-radiation coupling flux as initial node values and closed-loop forward propagation based on the propagation graph weight matrix and convergence coefficients based on spectral radius), an audit weight mapping module (removing topographic bias from the local slope of the incident light factor sequence using the vegetation index sequence to obtain the debiased vegetation index, mapping the mismatch potential energy to audit weights based on the regional mean), and a calculation result module (calculating the interannual variation of the debiased vegetation index, calculating the soil retention increment by combining the normalized sediment yield intensity sequence and rain erosion force, and performing regional weighted summation of the interannual variation and soil retention increment using audit weights to obtain the calculation result).
[0008] The beneficial effects of this invention are as follows: Based on the influence of slope aspect on surface processes and the hydrological connectivity characteristics of watersheds, this invention effectively eliminates the linear bias of topographic illumination on vegetation indices by splitting data according to slope aspect, capturing the correlation differences between sunlight, vegetation, and sediment yield, and combining closed forward propagation to quantify the cumulative impact of mismatches. Simultaneously, it adapts the mismatch risk differences of different units with audit weights. This allows the calculation process to fully incorporate slope aspect system bias and spatial propagation effects, avoiding biases caused by isolated calculations and data mixing, making the vegetation index more closely reflect the actual coverage state, and the regional summary results can reflect the actual impact of different units. The final generated total increase in regional vegetation cover and total increase in soil retention can provide more realistic basic data for the assessment of the value of soil and water conservation ecological products, providing reliable support for related protection decisions and management work. Attached Figure Description
[0009] Figure 1 is a flowchart of a method for calculating the value of soil and water conservation ecological products according to the present invention; Figure 2 is a schematic diagram of a system for calculating the value of soil and water conservation ecological products according to the present invention.
[0010] In the diagram: Data processing module 201, Coupling flux generation module 202, Propagation graph weight matrix calculation module 203, Mismatch potential energy calculation module 204, Audit weight mapping module 205, Accounting result calculation module 206. Detailed Implementation
[0011] The subject matter described herein will now be discussed with reference to exemplary embodiments. It should be understood that these embodiments are discussed only to enable those skilled in the art to better understand and implement the subject matter described herein, and changes may be made to the function and arrangement of the elements discussed without departing from the scope of this specification. Various processes or components may be omitted, substituted, or added as needed in the examples. Furthermore, features described in some examples may be combined in other examples.
[0012] It should be noted that, unless otherwise defined, the technical or scientific terms used in one or more embodiments of the present invention should have the ordinary meaning understood by one of ordinary skill in the art to which this invention pertains. The terms "first," "second," and similar terms used in one or more embodiments of the present invention do not indicate any order, quantity, or importance, but are merely used to distinguish different components. Terms such as "comprising" or "including" indicate that the element or object preceding the term encompasses the elements or objects listed following the term and their equivalents, without excluding other elements or objects. Terms such as "connected" or "linked" are not limited to physical or mechanical connections, but can include electrical connections, whether direct or indirect. Terms such as "upper," "lower," "left," and "right" are used only to indicate relative positional relationships; when the absolute position of the described object changes, the relative positional relationship may also change accordingly.
[0013] As shown in Figures 1 and 2, a method for calculating the value of soil and water conservation ecological products includes the following steps: Step S101: Within the pre-flood season window, the slope unit is divided into two parts according to slope aspect, and the vegetation index sequence, incident light factor sequence, and normalized sediment yield intensity sequence are obtained. The unit area and rain erosion force are obtained, and the difference between vegetation index and incident light factor are calculated respectively, and the difference in normalized sediment yield intensity is calculated; Step S102: The local regression sensitivity difference between the difference between vegetation index and the difference between incident light factor is calculated in the low coverage and high coverage intervals respectively. The sensitivity difference is weighted by combining the difference in normalized sediment yield intensity to generate the aspect irradiance coupling flux; Step S103: A hydrological connectivity matrix is established, and the aspect irradiance coupling flux is calculated based on the... Step S104: Using the aspect ratio irradiance coupling flux as the initial value of the nodes, closed forward propagation is performed through the propagation graph weight matrix and the convergence coefficient based on the spectral radius to calculate the mismatch potential energy. Step S105: The topographic bias is removed by removing the local slope of the incident light factor sequence from the vegetation index sequence to obtain the debiased vegetation index. The mismatch potential energy is mapped to the audit weight based on the regional mismatch potential energy mean. Step S106: The interannual variation of the debiased vegetation index is calculated. The soil retention increment is calculated by combining the normalized sediment yield intensity sequence and the rain erosion force. The interannual variation and the soil retention increment are regionally weighted and summarized by the audit weight to obtain the accounting result.
[0014] In one embodiment of the present invention, within the pre-flood season window, slope units are divided into two parts according to slope aspect, and vegetation index sequences, incident light factor sequences, and normalized sediment yield intensity sequences are obtained. Unit area and rain erosion capacity are obtained, and the differences in vegetation indices and incident light factors are calculated, along with the difference in normalized sediment yield intensity. This includes: within the pre-flood season window... Inside, slope unit Construct a cell set containing all cells. Based on slope aspect data, the pixel set is... Divided into mutually exclusive sets of sunny slopes Collection with shady slopes ,satisfy and Obtain the area of the slope unit. Time series of rain erosion forces Obtain pixel-level vegetation index sequences With pixel-level incident illumination factor sequence For the set of pixels The arithmetic mean was calculated to obtain the vegetation index sequence for the slope unit. With slope unit incident light factor sequence : ; ;in For pixel position, Pre-flood season window Inner moment, For a set of pixels Number of internal pixels; Obtain the sediment yield intensity sequence per unit area normalized to rainfall erosion force. For the sunny slope set Collection with shady slopes Perform arithmetic mean calculation to obtain the vegetation index mean sequence for sunny slopes. Vegetation index mean sequence for shady slopes Incident light factor mean sequence of sunny slope Incident light factor mean sequence of shaded slope Normalized sediment yield intensity mean sequence for sunny slopes and normalized sediment yield intensity mean sequence for shady slopes : ; ; ; ; ; ;in Set on the sunny slope Number of internal pixels Collection of shady slopes Number of internal pixels; Calculation of vegetation index difference Incident light factor difference and normalized sediment yield intensity difference : ; ; .
[0015] It should be noted that slope units represent the basic spatial units for accounting for the value of soil and water conservation ecological products, reflecting sloping areas with relatively consistent hydrodynamic and topographic attributes. Data collection is based on Digital Elevation Model (DEM) data, using GIS software hydrological analysis tools to extract watersheds and contour lines, dividing the area according to river flow direction and isogradability lines. After division, it is ensured that the fluctuation range of attributes such as slope, slope length, and confluence path within the unit does not exceed 10%. Slope aspect data can be extracted from 30-meter resolution DEM data at the pixel level, with azimuth values ranging from 0 to 360 degrees, calculated using GIS software surface analysis tools. Rainfall erosivity time series can be obtained using a daily rainfall erosivity calculation model, utilizing daily rainfall data within the pre-flood season window, calculating daily rainfall erosivity values using formulas, and arranging them in chronological order to form a series. Rainfall data originates from measured data from meteorological stations within the study area or from a meteorological data sharing platform. Pixel-level vegetation index sequences can be generated from Landsat or Sentinel series multispectral remote sensing images. After preprocessing such as atmospheric and geometric corrections, pixel-level data is retrieved using the vegetation index calculation formula and stitched together chronologically to form a sequence. The image temporal resolution must match the monitoring frequency within the pre-flood season window, ensuring at least one image every 16 days. Pixel-level incident light factor sequences can be generated based on solar altitude angle, solar azimuth angle, pixel aspect, and slope data. The cosine of the angle between solar radiation and the surface normal is calculated using the cosine theorem, and combined with the remote sensing image imaging time to calculate data at each moment, forming a sequence. Pixel-level normalized sediment yield intensity sequences can be obtained from monitoring stations. After each rainfall event, the sediment yield in the control area is measured, and the sediment yield per unit area is calculated and divided by the concurrent rain erosion force value to obtain the monitoring station's pixel data. For pixels without monitoring stations, the inverse distance weighted interpolation method is used, based on data from at least three surrounding monitoring stations, with the interpolation results calculated by spatial distance weighting to form a complete sequence. The above-mentioned calculation models for daily rainfall erosivity, pixel-level vegetation index, pixel-level incident light factor, and sediment yield intensity are all existing technical methods and will not be elaborated here.
[0016] It should be noted that the cell set represents the collection of all remote sensing image cells within a single slope unit. Slope aspect data represents the orientation of the land surface towards solar radiation. The sunny slope cell set represents the subset of cells in the cell set whose topographic orientation makes them more susceptible to receiving solar radiation. The shady slope cell set represents the subset of cells in the cell set whose topographic orientation makes them less susceptible to receiving solar radiation. The slope unit area value represents the actual spatial area of a single slope unit. The pre-flood season window represents the main rainfall period in South China from April to June, reflecting the key time range for active soil sedimentation and erosion. The rain erosivity time series represents the data sequence of the erosive capacity of rainfall on the soil at various times within the pre-flood season window. The cell-level vegetation index series represents the data sequence related to the vegetation cover of each cell in the remote sensing image at different times, reflecting the temporal changes in the vegetation growth and cover status of each cell. The cell-level incident light factor series represents the data sequence related to the angle of solar radiation received by each cell at different times, reflecting the temporal changes in the light conditions of each cell. The pixel-level normalized sediment yield intensity sequence represents the sediment yield data sequence per unit area of each pixel after being normalized by rain erosion force at different times, reflecting the temporal variation of sediment yield intensity of each pixel.
[0017] It should be noted that the slope unit vegetation index sequence represents the data sequence of vegetation indices for all pixels within a single slope unit. The slope unit incident light factor sequence represents the data sequence of incident light factors for all pixels within a single slope unit. The vegetation index sunny slope mean sequence represents the average data sequence of vegetation indices for all pixels within the sunny slope pixel set. The vegetation index shady slope mean sequence represents the average data sequence of vegetation indices for all pixels within the shady slope pixel set. The incident light factor sunny slope mean sequence represents the average data sequence of incident light factors for all pixels within the sunny slope pixel set. The incident light factor shady slope mean sequence represents the average data sequence of incident light factors for all pixels within the shady slope pixel set. The normalized sediment yield intensity sunny slope mean sequence represents the average data sequence of normalized sediment yield intensity for all pixels within the sunny slope pixel set. The normalized sediment yield intensity shady slope mean sequence represents the average data sequence of normalized sediment yield intensity for all pixels within the shady slope pixel set. The vegetation index difference sequence represents the data sequence of the difference between the mean vegetation index of sunny and shady slopes, reflecting the difference in vegetation cover status caused by slope aspect. The incident light factor difference sequence represents the data sequence of the difference between the mean incident light factor of sunny and shady slopes, reflecting the difference in light reception conditions caused by slope aspect. The normalized sediment yield intensity difference sequence represents the data sequence of the difference between the mean normalized sediment yield intensity of sunny and shady slopes, reflecting the difference in sediment yield intensity caused by slope aspect.
[0018] It should be noted that the slope aspect classification adopts the industry-standard azimuth angle classification, with 0 degrees as true north. The azimuth angle is calculated by rotating clockwise. 337.5 degrees to 22.5 degrees is the north slope, 22.5 degrees to 67.5 degrees is the northeast slope, 67.5 degrees to 112.5 degrees is the east slope, 112.5 degrees to 157.5 degrees is the southeast slope, 157.5 degrees to 202.5 degrees is the south slope, 202.5 degrees to 247.5 degrees is the southwest slope, 247.5 degrees to 292.5 degrees is the west slope, and 292.5 degrees to 337.5 degrees is the northwest slope. Among them, the south, southeast, and southwest slopes are classified as sunny slope pixel sets, while the north, northeast, and northwest slopes are classified as shady slope pixel sets, ensuring that the two sets are mutually exclusive and completely cover the original pixel sets. For example, if a pixel has an azimuth angle of 180 degrees, it belongs to the south slope and is classified into the sunny slope pixel set; if a pixel has an azimuth angle of 0 degrees, it belongs to the north slope and is classified into the shady slope pixel set.
[0019] In one embodiment of the present invention, the local regression sensitivity difference between the vegetation index difference and the incident light factor difference in the low-coverage and high-coverage intervals is calculated respectively. The sensitivity difference is then weighted by combining the normalized sediment yield intensity difference to generate the slope aspect irradiance coupling flux, including: slope units. Vegetation index series Calculate the 20th percentile value 40th place 70th place and the 90th place value Define a set of low-coverage intervals. High coverage interval set : ; ;in For the pre-flood season window; clustered in low-coverage areas. Inside, based on the difference in incident light factor The independent variable is the difference in vegetation indices. Perform linear regression calculations on the dependent variable to obtain the local regression slope for low coverage. It satisfies the minimum objective function: ;in For the regression intercept, The slope to be solved; set in the high coverage interval. Inside, based on the difference in incident light factor The independent variable is the difference in vegetation indices. Perform linear regression calculations on the dependent variable to obtain the high-coverage local regression slope. It satisfies the minimum objective function: ;in The regression intercept is used to calculate the difference in local regression sensitivity. : ; for the pre-flood season window Normalized sediment yield intensity difference within The normalized sediment yield intensity weighted value is obtained by performing an arithmetic mean calculation. : ;in The number of time steps within the pre-flood season window; utilizing the difference in local regression sensitivity. Multiply by normalized sediment yield intensity weighted value Generate slope aspect irradiation coupling flux : .
[0020] It should be noted that the 20th, 40th, 70th, and 90th decimal places represent the cumulative percentages of the slope unit vegetation index series when arranged in ascending order, reaching 20%, 40%, 70%, and 90%, respectively. The low-coverage interval set represents the set of times when the slope unit vegetation index falls between the 20th and 40th decimal places. The high-coverage interval set represents the set of times when the slope unit vegetation index falls between the 70th and 90th decimal places. The incident light factor difference series reflects the differences in light conditions caused by slope aspect differences. The vegetation index difference series reflects the differences in vegetation cover status caused by slope aspect differences. The local regression sensitivity difference represents the difference in the local regression slope between low and high cover, reflecting the difference in the response intensity of slope aspect vegetation differences to slope aspect light differences under different vegetation cover levels. The normalized sediment yield intensity difference series reflects the differences in sediment yield intensity caused by slope aspect differences. The normalized sediment yield intensity weighted average reflects the average level of slope aspect sediment yield differences during the pre-flood season. The aspect-coupled irradiance flux represents the product of the local regression sensitivity difference and the normalized weighted average sediment yield intensity, reflecting the coupling relationship between vegetation cover and light response differences and aspect-coupled sediment yield differences.
[0021] It should be noted that the 20th to 40th percentiles correspond to vegetation cover of 30% to 50%, at which point the proportion of bare land is relatively high, and soil iron oxides and bare land have a more significant impact on remote sensing imagery. The 70th to 90th percentiles correspond to vegetation cover of 60% to 80%, at which point the vegetation canopy is dominant, and the influence of bare land and iron oxides is weaker. If the 10th to 30th percentiles are chosen, the low cover range is too narrow, and outliers have a greater impact. If the 80th to 100th percentiles are chosen, the amount of high cover data is insufficient, and the stability of the regression slope is poor. The 20th to 40th percentiles and the 70th to 90th percentiles can achieve a balance between data volume and representativeness, accurately distinguishing the differences in surface response under different cover levels. For example, in a certain region, the 20th percentile of the slope unit vegetation index sequence corresponds to a cover of 32%, the 40th percentile to 48%, the 70th percentile to 63%, and the 90th percentile to 79%. This range can clearly separate the two periods of bare land dominance and canopy dominance.
[0022] It should be noted that, given the differences in surface response to light and soil erosion processes under different vegetation cover levels, we divided the period into low and high cover periods using quantile values, separately calculated the slope of the response to light differences based on vegetation differences, and then combined this with the average level of aspect-dependent sediment yield differences to generate aspect-dependent irradiance coupled flux. This process separates the surface response characteristics under different cover levels, quantifies the correlation between light, vegetation, and sediment yield, avoids the ambiguity of response characteristics caused by mixing data from different cover levels, and ensures that the coupled flux accurately reflects the correlation strength of key surface processes, thereby improving the accuracy of subsequent quantitative parameter calculations.
[0023] In one embodiment of the present invention, a hydrological connectivity matrix is established, and connecting edges are selected and assigned values based on the sign-homogeneity of slope aspect irradiance coupling flux. The edge weights are then normalized to obtain a propagation graph weight matrix, including: setting up slope units... Defined as a set of nodes, a hydrological connectivity matrix is established. For any pair of slope elements and If slope unit With slope unit There exist upstream-to-downstream confluence relationships or adjacent confluence relationships, making the hydrological connectivity matrix elements Otherwise, let the elements of the hydrological connectivity matrix... ; Calculate temporary edge weights : ;in and Slope unit With slope unit The corresponding slope aspect irradiance coupling flux, This indicates selecting the larger value between zero and the result of the expression within the parentheses; for temporary edge weights Perform normalization to obtain the elements of the propagation graph weight matrix. Elements of the propagation graph weight matrix Constructing the propagation graph weight matrix: ;in Slope unit The sum of all temporary edge weights within the corresponding row. It is a numerically stable constant and .
[0024] It should be noted that the node set represents the collective term for the nodes in the graph structure, mapping the set of slope units. The upstream-to-downstream confluence relationship indicates a unidirectional flow between two slope units, reflecting the natural conduction path of the water flow. The adjacent confluence relationship indicates that two slope units are adjacent at their boundaries and share a confluence path. The hydrological connectivity matrix represents a matrix characterizing the confluence association state between slope units, used to limit the set of edges allowed for connection in the graph structure. The temporary edge weight represents the initial edge weight after filtering by hydrological connectivity state and coupling flux homogeneity, reflecting the initial strength of the propagation association between nodes. The numerically stable constant represents a small positive number (e.g., 10 to the power of -4) used to avoid zero denominators in calculations. The propagation graph weight matrix elements represent the normalized propagation weights between nodes, reflecting the relative strength of mutual influence between nodes. The propagation graph weight matrix represents a matrix composed of the elements of the propagation graph weight matrix, used for subsequent closed-loop forward propagation calculations of node states.
[0025] It should be noted that the upstream-to-downstream confluence relationship can be determined using GIS hydrological analysis tools. Specifically, 30-meter resolution DEM data of the study area is acquired, and preprocessed by filling depressions and smoothing. Then, water flow direction grids are extracted, and the D8 algorithm is used to determine the water flow direction of each cell. Next, the confluence path is extracted based on the water flow direction grids to determine the outlet location of each slope unit. Finally, for any two slope units, if the water flow from the outlet of the first slope unit eventually flows into the inlet of the second slope unit, and there are no other slope units intervening in the confluence path, then the first slope unit and the second slope unit are determined to have an upstream-to-downstream confluence relationship. For example, if the water flow from the outlet of slope unit A directly flows into the inlet of slope unit B, and there are no other slope units in between, then A and B are determined to have an upstream-to-downstream confluence relationship.
[0026] It should be noted that the determination of adjacent confluence relationships requires two conditions: 1. Spatial adjacency: the vector boundaries of the two slope units have overlapping line segments, and the overlap length is not less than 50 meters; 2. Hydrological correlation: the confluence paths of the two slope units intersect in the boundary overlap area, and the water flow direction is consistent. For example, the boundary overlap length of slope unit C and slope unit D is 80 meters, and their confluence paths converge in the overlap area and flow to the same downstream unit, thus determining that C and D have an adjacent confluence relationship. In addition, based on the propagation law of soil and water conservation mismatch effects within the watershed, the sign of the aspect irradiance coupling flux reflects the consistency of the mismatch type. Same sign indicates the same mismatch type, and opposite sign indicates the opposite mismatch type. The mismatch effect can only propagate to units with the same mismatch type within the watershed through hydrological paths. The mismatch effects between units with opposite signs will cancel each other out and cannot form an effective propagation. Therefore, this invention only retains the connection edges between units with the same sign to ensure that the propagation weight matrix can accurately depict the actual propagation path of the mismatch effect.
[0027] In one embodiment of the present invention, using the slope aspect irradiance coupling flux as the initial value of the nodes, closed forward propagation is performed through the propagation graph weight matrix and the convergence coefficient based on the spectral radius to calculate the mismatch potential energy, including: dividing the slope element... Corresponding slope aspect irradiance coupling flux Arrange them in order and construct the initial value vector of the nodes. : ;in The number of slope units. to This represents the slope aspect irradiance coupling flux corresponding to each slope unit. Represents vector transpose operation; calculates the propagation graph weight matrix. spectral radius : ;in For the propagation graph weight matrix, For the propagation graph weight matrix The set of all eigenvalues, For eigenvalues, This indicates selecting the maximum absolute value among the eigenvalues; based on the spectral radius. Calculate the convergence coefficient : ; build identity matrix of order Perform closed-loop forward propagation calculations to obtain the mismatch potential energy vector. : ;in Represents a matrix Perform the inversion operation; extract the mismatch potential energy vector. The amount As a slope unit The mismatch potential energy, of which .
[0028] It should be noted that the spectral radius represents the maximum absolute value among all eigenvalues of the propagation graph weight matrix. The convergence coefficient represents the coefficients used to ensure the convergence of the propagation series. The identity matrix represents a square matrix with diagonal elements of 1 and other elements of 0. The decay matrix represents the product of the convergence coefficient and the propagation graph weight matrix, reflecting the decay characteristics of the propagation intensity between nodes. The difference matrix represents the difference between the identity matrix and the decay matrix, reflecting the baseline deviation matrix after the propagation intensity decays. Matrix inversion represents the method of solving the inverse matrix based on matrix operation rules, used to transform the difference matrix into an invertible inverse matrix. The inverse matrix represents the result of the inverse operation of the difference matrix. The propagation operator matrix represents the difference between the inverse matrix and the identity matrix. The mismatch potential energy vector represents the product of the propagation operator matrix and the node initial value vector. The mismatch potential energy represents a single component in the mismatch potential energy vector, reflecting the cumulative effect of mismatches from multiple upstream or neighboring orders in a single slope element.
[0029] It should be noted that the convergence coefficient is used to ensure the convergence of the closed forward propagation series, and the key condition for series convergence is that the spectral radius after multiplying the gamma by the propagation graph weight matrix is less than one. The spectral radius is the maximum absolute value of all eigenvalues of the propagation graph weight matrix, and its value is usually between zero and one. Choosing 0.9 avoids both excessive attenuation of the propagation effect due to a value that is too small, thus missing the role of higher-order indirect propagation, and excessive value that would cause the spectral radius after multiplying the gamma by the propagation graph weight matrix to be close to or greater than one, thereby causing series divergence. It achieves a balance between preserving the complete propagation effect and ensuring computational convergence, ensuring that the calculation results are stable and consistent with reality. In addition, based on the multi-order propagation characteristics of the graph structure, using the slope aspect irradiance coupling flux as the initial value of the nodes, the convergence coefficient is determined by the spectral radius to ensure series convergence. By constructing the propagation operator matrix through matrix operations, the closed forward propagation and accumulation of mismatch effects are realized. This allows for the quantification of the multi-order upstream and neighborhood mismatch effects on each slope unit without the need for training or iteration. The use of the identity matrix and matrix inversion ensures the rationality of the propagation operator matrix, and the derivation of the convergence coefficient avoids divergence of results. The obtained mismatch potential energy can accurately reflect the spatial cumulative effect of mismatch influence, providing a reliable basis for the construction of subsequent audit weights.
[0030] In one embodiment of the present invention, a debiased vegetation index is obtained by removing topographic bias from the local slope of the incident light factor sequence using a vegetation index sequence. The mismatch potential energy is then mapped to audit weights based on the mean regional mismatch potential energy, including: during the pre-flood season window... Inside, with slope unit Incident light factor sequence The independent variable is a vegetation index series. Perform linear regression calculations on the dependent variable to obtain the local first-order slope. The local first-order slope satisfies the minimization condition: ;in For the regression intercept, The slope to be solved; for the incident illumination factor sequence Perform an arithmetic mean calculation to obtain the mean incident light factor. : ;in This represents the number of time steps within the pre-flood season window; the vegetation index sequence Subtracting the terrain bias term yields the debiased vegetation index. : Obtain the mismatch potential energy vector composed of the mismatch potential energy of all slope elements. Calculate the mean of the mismatch potential energy in the region. : ;in The number of slope units. For the first Mismatch potential energy of each slope unit The absolute value of the mismatch potential energy; using the mean of the mismatch potential energy in the region. Slope unit Mismatch potential energy Mapping is performed to obtain audit weights. : ;in It is a numerically stable constant and , Slope unit The absolute value of the mismatch potential energy.
[0031] It should be noted that the local first-order slope represents the linear response coefficient of the vegetation index to the incident light factor within the pre-flood season window, reflecting the intensity of the first-order change in the vegetation index affected by light. The light fluctuation sequence represents the difference sequence between the incident light factor sequence and its actual mean, reflecting the degree of fluctuation of light conditions at each time point relative to the average level. The topographic bias sequence represents the product sequence of the light fluctuation sequence and the local first-order slope, reflecting the degree of linear interference of topographic light fluctuation on the vegetation index. The debiased vegetation index represents the vegetation index sequence after removing the linear interference of topographic light. The regional mismatch potential energy mean represents the arithmetic mean of the absolute values of mismatch potential energy of all slope units, reflecting the average level of overall mismatch risk in the study area. The relative mismatch ratio represents the ratio of the mismatch potential energy of a single slope unit to the regional average level, reflecting the relative intensity of unit mismatch risk. The audit weight represents the accounting weighting parameter obtained based on the relative mismatch ratio, reflecting the weight of the impact of unit mismatch risk on the regional accounting results.
[0032] It should be noted that the calculation period for the local first-order slope must be completely consistent with the pre-flood season window. That is, only the incident light factor sequence and vegetation index sequence data within the pre-flood season window should be used for regression, and no data outside the window should be included. This ensures that the slope can accurately reflect the response characteristics of vegetation indices to light during the key sand-producing season. For example, if the pre-flood season window is from April to June, only the sequence data from this period should be used to calculate the local first-order slope to avoid the deviation of response coefficients caused by data from non-key seasons. Adding 1 to the relative mismatch ratio to obtain the audit weight ensures that the audit weight is not less than 1, preventing the contribution of units prone to mismatch from being weakened due to a weight less than 1. The relative mismatch ratio reflects the multiple of a unit's mismatch risk relative to the regional average level. After adding 1, the weight corresponding to the regional average mismatch risk is 2. Units with mismatch risk higher than the average level have a weight greater than 2, while units with mismatch risk lower than the average level have a weight between 1 and 2. This preserves the weight differences between units while ensuring that the contribution of all units is not excessively suppressed. For example, if a unit has a relative mismatch ratio of 0.8, the corresponding audit weight is 1.8; if another unit has a relative mismatch ratio of 1.5, the corresponding audit weight is 2.5, thus distinguishing the weights of units with different mismatch risks.
[0033] It should be noted that, based on the linear interference characteristics of topographic illumination on vegetation indices, the local first-order slope is solved by linear regression to eliminate topographic bias caused by illumination. At the same time, the mismatch potential energy of a single unit is mapped to an audit weight based on the regional mismatch potential energy mean. This results in bias-free data that reflects the true state of vegetation, as well as accounting weights that are adapted to mismatch risks, avoiding interference of topographic illumination bias on vegetation information and making the vegetation index more consistent with the actual coverage state.
[0034] In one embodiment of the present invention, the interannual variation of the debiased vegetation index is calculated, and the soil retention increment is calculated by combining the normalized sediment yield intensity sequence and rain erosion capacity. The interannual variation and soil retention increment are then regionally weighted and aggregated using audit weights to obtain the accounting results, including: [details of the calculation for the current year are missing from the original text]. Compared with the base year Pre-flood season window Within the de-biased vegetation index and The average vegetation index for the current year is calculated by averaging over the execution time. Compared with the baseline year mean vegetation index : ; ;in This refers to the number of time steps within the pre-flood season window. For time; calculate slope unit interannual variation : Get the current year Compared with the base year Normalized sediment yield intensity sequence and Execute the slope unit cell set respectively Spatial average calculation within and pre-flood season window The average sediment yield intensity for the current year is obtained by averaging over a period of time. Compared with the average sediment yield intensity of the baseline year : ; ;in For slope unit cell set Number of internal pixels For pixel location; calculate the current year's soil erosion. Soil loss in the baseline year And calculate the soil retention increment. : ; ; ;in and These represent the pre-flood season erosion forces for the current year and the baseline year, respectively. For slope unit area; using audit weights The total increase in regional vegetation cover was obtained by performing a weighted aggregation of all slope units. Maintaining the total increase in regional soil : ; ;in This represents the number of slope units.
[0035] It should be noted that the mean adjusted vegetation index for the baseline year represents the arithmetic mean of the adjusted vegetation index series within the pre-flood season window of the baseline year, reflecting the average level of actual vegetation cover in the baseline year. The interannual variation of the adjusted vegetation index represents the difference between the mean adjusted vegetation index for the current year and the baseline year, reflecting the degree of interannual variation in vegetation cover status. The mean sediment yield intensity for the baseline year represents the result of spatially and temporally averaging the normalized sediment yield intensity series for the baseline year, reflecting the average sediment yield intensity of slope units in the baseline year. The pre-flood season erosivity for the baseline year represents the cumulative erosive capacity of rainfall on the soil within the pre-flood season window of the baseline year, reflecting the total erosion intensity of rainfall in the baseline year. The soil loss for the baseline year represents the total soil loss caused by erosion in the slope unit in the baseline year, reflecting the actual scale of soil erosion in the baseline year. The increase in soil retention represents the difference between soil loss in the baseline year and the current year. The total increase in regional vegetation cover represents the cumulative interannual change in vegetation cover of all slope units after audit weighting and area weighting, reflecting the overall effect of improving vegetation cover in the region. The total increase in regional soil conservation represents the cumulative increase in physical soil conservation across all slope units after audit weighting, reflecting the overall effectiveness of soil conservation in the region.
[0036] It should be noted that the selection of the base year must meet three conditions: first, it must be 3 to 10 years away from the current year to ensure that the time span is sufficient to reflect the changes in the value of ecological products; second, there must be no major soil and water conservation projects implemented in the region during the base year, and no extreme weather events such as torrential rain or severe drought, to ensure the benchmark nature of the data; and third, the data sources of the base year and the current year must be consistent, and the resolution, data format, and distribution of monitoring stations of remote sensing images must be consistent to avoid accounting deviations caused by data differences.
[0037] It should be noted that the area differences of different slope units will lead to different contributions of their vegetation cover changes to the overall region. The larger the area of the slope unit, the more significant the impact of its vegetation cover changes on the total regional impact. Introducing area weighting can make the total increase in regional vegetation cover truly reflect the comprehensive contribution of slope units of different sizes, avoiding the unreasonable situation of small area units having the same contribution as large area units. For example, slope unit A has an area of 10 hectares and an interannual variation of 0.1; slope unit B has an area of 20 hectares and an interannual variation of 0.1. After area weighting, the contribution of A is 1.0 and the contribution of B is 2.0, which is more in line with the actual situation of the region. When the increase in soil retention volume is negative, it indicates that the soil loss in the current year is greater than that in the baseline year, meaning that the soil and water conservation effect has deteriorated. In this case, the negative value should be retained for regional aggregation, without being zeroed out or taken as an absolute value, to ensure that the accounting results truly reflect the positive and negative trends of changes in the value of ecological products. For example, if the soil loss in a certain slope unit is 500 tons in the baseline year and 600 tons in the current year, the increase in soil retention volume is -100 tons. This negative value should be accurately included in the regional total to reflect the degradation of the unit.
[0038] In one embodiment of the present invention, as shown in Figure 2, a soil and water conservation ecological product value accounting system includes: a data processing module 201, which divides slope units into two parts according to slope aspect within the pre-flood season window, obtains vegetation index sequences, incident light factor sequences, and normalized sediment yield intensity sequences, obtains unit area and rain erosion force, calculates the difference in vegetation index and the difference in incident light factor respectively, and calculates the difference in normalized sediment yield intensity; a coupling flux generation module 202, which calculates the local regression sensitivity difference of vegetation index difference to incident light factor difference in low coverage and high coverage intervals respectively, and weights the sensitivity difference with the normalized sediment yield intensity difference to generate slope aspect irradiance coupling flux; and a propagation graph weight matrix calculation module 203, which establishes a hydrological connectivity matrix and calculates the propagation graph weight matrix based on the slope aspect irradiance coupling flux. The flux is filtered for connection edges based on sign and direction and assigned values. The edge weights are normalized to obtain the propagation graph weight matrix. The mismatch potential energy calculation module 204 uses the aspect-radiation coupling flux as the initial value of the nodes and performs closed forward propagation through the propagation graph weight matrix and the convergence coefficient based on the spectral radius to calculate the mismatch potential energy. The audit weight mapping module 205 obtains the debiased vegetation index by removing the topographic bias from the local slope of the incident light factor sequence using the vegetation index sequence. The mismatch potential energy is mapped to the audit weight based on the regional mismatch potential energy mean. The accounting result calculation module 206 calculates the interannual variation of the debiased vegetation index, calculates the soil retention increment by combining the normalized sediment yield intensity sequence and rain erosion force, and performs regional weighted summation of the interannual variation and soil retention increment through the audit weight to obtain the accounting result.
[0039] It should be noted that the interval and threshold sizes are set for ease of comparison. The size of the threshold depends on the amount of sample data and the base number set by those skilled in the art for each set of sample data, as long as it does not affect the proportional relationship between the parameter and the quantized value. Furthermore, the above formulas are all dimensionless calculations, and the formulas are derived from software simulations using a large amount of collected data to obtain the most recent real-world results. The preset parameters in the formulas are set by those skilled in the art according to the actual situation.
[0040] The embodiments of this example have been described above. However, this example is not limited to the specific implementation methods described above. The specific implementation methods described above are merely illustrative and not restrictive. Those skilled in the art can make many other forms based on the guidance of this example, and all of them are within the protection scope of this example.
Claims
1. A method for calculating the value of soil and water conservation ecological products, characterized in that, Includes the following steps: Step S101: Within the pre-flood season window, divide the slope unit into two parts according to slope aspect, obtain the vegetation index sequence, incident light factor sequence and normalized sediment yield intensity sequence, obtain the unit area and rain erosion force, calculate the difference in vegetation index and the difference in incident light factor respectively, and calculate the difference in normalized sediment yield intensity. Step S102: Calculate the local regression sensitivity difference between the vegetation index difference and the incident light factor difference in the low and high coverage intervals, respectively. Weight the sensitivity difference using the normalized sediment yield intensity difference to generate the aspect radiation coupling flux. Step S103: Establish a hydrological connectivity matrix. Based on the sign-homogeneity of the aspect radiation coupling flux, filter and assign values to connecting edges. Normalize the edge weights to obtain the propagation graph weight matrix. Step S104: Using the aspect radiation coupling flux as the initial node value, and through the propagation graph weight matrix and... Closed forward propagation is performed based on the convergence coefficient of the spectral radius to calculate the mismatch potential energy; in step S105, the topographic bias is removed by the local slope of the incident light factor sequence through the vegetation index sequence to obtain the debiased vegetation index, and the mismatch potential energy is mapped to the audit weight according to the regional mismatch potential energy mean; in step S106, the interannual variation of the debiased vegetation index is calculated, and the soil retention increment is calculated by combining the normalized sediment yield intensity sequence and rain erosion force. The interannual variation and soil retention increment are regionally weighted and summarized by the audit weight to obtain the accounting result.
2. The method for calculating the value of soil and water conservation ecological products according to claim 1, characterized in that, Within a slope unit, a pixel set containing all pixels is constructed. Based on the slope aspect data, the pixel set is divided into mutually exclusive sunny slope pixel sets and shady slope pixel sets. The slope unit area values and the rain erosion force time series within the pre-flood season window are extracted. Pixel-level vegetation index sequences, pixel-level incident light factor sequences, and pixel-level normalized sediment yield sequences are extracted. The arithmetic mean of the pixel set is calculated to obtain the slope unit vegetation index sequence and the slope unit incident light factor sequence.
3. The method for calculating the value of soil and water conservation ecological products according to claim 2, characterized in that, Arithmetic mean calculations were performed on the pixel-level vegetation index sequences, pixel-level incident light factor sequences, and pixel-level normalized sediment yield sequences within the pixel sets of sunny and shady slopes, respectively, to obtain the vegetation index mean sequence for sunny slopes, the vegetation index mean sequence for shady slopes, the incident light factor mean sequence for sunny slopes, the incident light factor mean sequence for shady slopes, the normalized sediment yield mean sequence for sunny slopes, and the normalized sediment yield mean sequence for shady slopes. The vegetation index difference sequence was obtained by subtracting the vegetation index mean sequence for shady slopes from the vegetation index mean sequence for sunny slopes, the incident light factor difference sequence was obtained by subtracting the incident light factor mean sequence for shady slopes from the incident light factor mean sequence for sunny slopes, and the normalized sediment yield difference sequence was obtained by subtracting the normalized sediment yield mean sequence for shady slopes from the normalized sediment yield mean sequence for sunny slopes.
4. The method for calculating the value of soil and water conservation ecological products according to claim 1, characterized in that, The 20th, 40th, 70th, and 90th percentile values of the vegetation index sequence for the slope unit are calculated. The times in the vegetation index sequence where the value is between the 20th and 40th percentile values are classified into the low coverage interval set, and the times in the vegetation index sequence where the value is between the 70th and 90th percentile values are classified into the high coverage interval set.
5. The method for calculating the value of soil and water conservation ecological products according to claim 4, characterized in that, Within the low-coverage interval set, linear regression was performed with the incident light factor difference sequence as the independent variable and the vegetation index difference sequence as the dependent variable, and the regression slope was extracted as the local regression slope for low-coverage intervals. Within the high-coverage interval set, linear regression was performed with the incident light factor difference sequence as the independent variable and the vegetation index difference sequence as the dependent variable, and the regression slope was extracted as the local regression slope for high-coverage intervals. The local regression sensitivity difference was obtained by subtracting the local regression slope for high-coverage intervals from the local regression slope for low-coverage intervals. The normalized sediment yield intensity difference sequence within the pre-flood season window is calculated by performing an arithmetic mean to obtain the normalized sediment yield intensity weighted quantity; Multiply the difference in local regression sensitivity by the normalized sediment yield intensity weighted value to obtain the aspect-coupled irradiance flux.
6. The method for calculating the value of soil and water conservation ecological products according to claim 1, characterized in that, The slope unit set is defined as a node set. For any pair of slope units, if there is an upstream-to-downstream confluence relationship or an adjacent confluence relationship between the slope units, the corresponding element in the hydrological connectivity matrix is set to one; otherwise, the corresponding element in the hydrological connectivity matrix is set to zero. The product of the slope aspect irradiance coupling fluxes corresponding to the two slope units is calculated. The larger value between the product and zero is selected and multiplied by the corresponding element in the hydrological connectivity matrix to obtain the temporary edge weights. The sum of all temporary edge weights in the corresponding row of the slope unit is calculated. The sum plus the numerically stable constant is used as the denominator, and the temporary edge weights are used as the numerator. Division is performed to obtain the elements of the propagation graph weight matrix. The propagation graph weight matrix is constructed from the elements of the propagation graph weight matrix.
7. The method for calculating the value of soil and water conservation ecological products according to claim 1, characterized in that, Arrange the slope aspect irradiance coupling fluxes corresponding to all slope elements in order to construct the node initial value vector; calculate all eigenvalues of the propagation graph weight matrix, select the maximum absolute value of the eigenvalues as the spectral radius, and divide the numerical value of 0.9 by the spectral radius to obtain the convergence coefficient. Construct an identity matrix with the same order as the propagation graph weight matrix. Multiply the convergence coefficients by the propagation graph weight matrix to obtain the decay matrix. Subtract the decay matrix from the identity matrix to obtain the difference matrix. Perform matrix inversion on the difference matrix to obtain the inverse matrix. Subtract the identity matrix from the inverse matrix to obtain the propagation operator matrix. Multiply the propagation operator matrix by the node initial value vector to obtain the mismatch potential energy vector. Extract the components in the mismatch potential energy vector as the mismatch potential energy.
8. The method for calculating the value of soil and water conservation ecological products according to claim 1, characterized in that, Within the pre-flood season window, linear regression was performed with the incident light factor sequence as the independent variable and the vegetation index sequence as the dependent variable, and the regression slope was extracted as the local first-order slope. The arithmetic mean of the incident light factor sequence was calculated to obtain the mean of the incident light factor. The incident light factor sequence was subtracted from the incident light factor mean to obtain the light fluctuation sequence. The light fluctuation sequence was multiplied by the local first-order slope to obtain the topographic bias sequence. The vegetation index sequence was subtracted from the topographic bias sequence to obtain the debiased vegetation index. The arithmetic mean of the absolute values of mismatch potential energy corresponding to all slope units was calculated as the regional mismatch potential energy mean. The absolute value of mismatch potential energy of a single slope unit was used as the numerator, and the regional mismatch potential energy mean plus a numerically stable constant was used as the denominator. The division was performed to obtain the relative mismatch ratio. The relative mismatch ratio was added to a value of one to obtain the audit weight.
9. The method for calculating the value of soil and water conservation ecological products according to claim 1, characterized in that, The arithmetic mean of the corrected vegetation index series within the current year's pre-flood season window and the baseline year's pre-flood season window is calculated to obtain the mean of the corrected vegetation index for the current year and the mean of the corrected vegetation index for the baseline year. The interannual variation of the corrected vegetation index is obtained by subtracting the mean of the corrected vegetation index for the baseline year from the mean of the corrected vegetation index for the current year. The spatial average within the slope unit pixel set and the temporal average within the pre-flood season window are calculated for the normalized sediment yield intensity series for the current year and the baseline year, respectively, to obtain the mean of the sediment yield intensity for the current year and the mean of the sediment yield intensity for the baseline year. The pre-flood season rainfall erosion intensity of the current year is multiplied by the current year's... The soil loss for the current year is obtained by multiplying the average sediment yield intensity by the slope unit area. The soil loss for the baseline year is obtained by multiplying the pre-flood season rain erosion intensity by the average sediment yield intensity for the baseline year and then by the slope unit area. The soil loss for the current year is obtained by subtracting the soil loss for the current year from the soil loss for the baseline year. The total increase in soil retention is obtained by multiplying the audit weight by the interannual change in the debiased vegetation index and then by the slope unit area, and summing the results of the product of all slope units. The total increase in regional vegetation cover is obtained by multiplying the audit weight by the increase in soil retention and summing the results of the product of all slope units. The total increase in regional soil retention is obtained by multiplying the audit weight by the increase in soil retention.
10. A system for calculating the value of soil and water conservation ecological products, characterized in that, The method for calculating the value of soil and water conservation ecological products as described in any one of claims 1 to 9 includes: a data processing module, which divides slope units into two parts according to slope aspect within the pre-flood season window, obtains vegetation index sequences, incident light factor sequences, and normalized sediment yield intensity sequences, obtains unit area and rain erosion capacity, calculates the difference in vegetation index and the difference in incident light factor respectively, and calculates the difference in normalized sediment yield intensity; a coupling flux generation module, which calculates the local regression sensitivity difference of vegetation index difference to incident light factor difference in low and high coverage intervals respectively, and weights the sensitivity difference with the normalized sediment yield intensity difference to generate slope aspect radiation coupling flux; and a propagation graph weight matrix calculation module, which establishes a hydrological connectivity matrix and calculates the propagation graph weight matrix based on the slope aspect radiation coupling flux. The system employs several methods: First, it filters and assigns values to connecting edges based on the sign-direction property of the quantities, then normalizes the edge weights to obtain the propagation graph weight matrix. Second, it calculates mismatch potential energy by using the aspect-radiative coupling flux as the initial node value and performing closed-loop forward propagation using the propagation graph weight matrix and convergence coefficients based on spectral radius. Third, it maps the mismatch potential energy to audit weights by removing topographic bias from the local slope of the incident light factor sequence using the vegetation index sequence and mapping the mismatch potential energy to audit weights based on the regional mismatch potential energy mean. Finally, it calculates the interannual variation of the debiased vegetation index, calculates the soil retention increment by combining the normalized sediment yield intensity sequence and rain erosion force, and performs regional weighted summation of the interannual variation and soil retention increment using audit weights to obtain the final accounting result.