A method for quantifying terrace pattern evolution and its ecosystem service trade-offs

Abstract

The present invention discloses a method for quantifying terrace pattern evolution and its ecosystem service trade-offs, involving the interdisciplinary field of agricultural ecological technology and geographic information technology. The method constructs a multi-source remote sensing information collaborative framework to obtain multi-source remote sensing data. The multi-source remote sensing data is preprocessed to obtain a standard input dataset. Based on the standard input dataset, a multimodal boundary refinement segmentation algorithm is used to extract quantitative spatiotemporal evolution data, such as the distribution of terraces in typical regions and their interannual expansion / contraction dynamics. Through multimodal deep learning and spatial statistical model linkage, a comprehensive quantitative analysis of the terrace evolution process is achieved. Based on the extracted terrace distribution data, combined with the standard input dataset and ground observation data, an evaluation framework is proposed, and a collaborative analysis of the trade-offs between terrace ecosystem services is performed. The present invention overcomes the problems of insufficient spatiotemporal resolution, static ecosystem service assessment, and weak decision support in the prior art.

Description

Technical Field

[0001] The present invention relates to the field of intersectional technology of agricultural ecological technology and geographic information technology, and more specifically to a method for quantifying terrace pattern evolution and balancing its ecosystem services. Background Art

[0002] Currently, research on terrace pattern monitoring and quantitative analysis techniques falls into two main categories. The first, morphological extraction techniques based on remote sensing image interpretation, are the most widely used. Their main characteristics are the use of high-resolution satellite imagery (such as Sentinel-2 and GF-2) combined with object-oriented classification or other automated classification algorithms to extract terrace boundaries. Morphological algorithms (such as edge detection and texture analysis) are then used to quantify geometric parameters (such as area, slope, and fragmentation). However, these techniques have two main drawbacks. First, they lack spatiotemporal dynamics. These methods focus on single-phase static analysis and lack analysis of long-term evolution and its driving forces, making them difficult to support comprehensive predictive optimization. Second, these techniques have limited applicability at mesoscales. These techniques are often designed for small-scale, high-precision (<10 km²) or large-scale, low-precision (>100 km²), lacking specialized modeling methods for mesoscale terrain units ranging from 1:5000 to 1:50,000. The second is scenario simulation technology based on geographic information systems. The main feature of this technology is that it often uses land change models (such as CLUE-S and CA-Markov) to simulate the expansion / contraction trends of terraces, and combines terrain characteristic indices (such as the terrain wetness index (TWI) and the slope variability index (SVI)) to evaluate the evolution characteristics and laws of spatial patterns. The shortcomings of this type of technical method are reflected in two aspects. First, the driving factors are simplistic. The model relies heavily on natural factors (such as topography and precipitation) and ignores the dynamic interactive effects of socioeconomic factors (such as labor transfer and policy compensation). In addition, the measurement module of service trade-offs is missing. The simulation results are mostly based on area changes. The dynamic response function of ecosystem service value is not coupled, and it cannot reveal the mechanism of the impact of pattern evolution on ecosystem services and their interrelationships.

[0003] Ecosystem service trade-off analysis methods can be categorized into two main groups. The first involves static ecosystem service valuation methods, which often use unit area value equivalents (such as the InVEST model) or emergy analysis (EMA) (Nadalini et al., 2021) to calculate ecosystem service values ​​and prioritize ecosystem services through spatial overlay analysis. This type of assessment is relatively common and widely used, but suffers from timeliness issues. They rely on base-year data and lack a dynamic time series correction mechanism, resulting in results lagging behind actual evolution. Furthermore, they simplify the synergistic trade-offs between ecosystem services, often using linear regression or correlation analysis to characterize service trade-offs, making it difficult to account for nonlinear threshold effects (e.g., when terrace fragmentation exceeds a critical value, its carbon sequestration function suddenly decreases). The second group involves decision-making techniques based on multi-objective optimization, which often employ Pareto frontiers (PFs) or genetic algorithms (GAs) to seek optimal solutions and balance conflicting service objective functions. The main drawback of this method is the lack of spatial visualization. The optimization results are mostly expressed at the statistical level, and there is a lack of implementation paths that match the space positions and terrain information of specific terrace patches. In addition, the dynamic adaptability is poor, and the results do not establish a feedback mechanism with the evolution process of the terrace pattern, making it difficult for the planning scheme to adapt to long-term policy or environmental conditions.

[0004] Overall, focusing on the trade-off and coordinated analysis of terrace pattern evolution and ecosystem services, previous ideas were mostly to combine remote sensing images and land use data to identify and extract target land types, and then to generate static ecosystem service value results based on the equivalent factor method, and to perform simple linear correlation analysis through spatial overlay.

[0005] Therefore, how to solve the high-precision dynamic identification and quantitative analysis of the evolution of terrace spatial pattern from a mesoscale spatial perspective, and to construct a multi-objective collaborative trade-off framework for ecosystem services to overcome the insufficient spatiotemporal resolution, static ecosystem service assessment and weak decision-making support in existing technologies are issues that technical personnel in this field urgently need to solve. Summary of the Invention

[0006] In view of this, the present invention provides a method for quantifying the evolution of terrace patterns and its ecosystem service trade-off to solve the technical problems mentioned in the background technology.

[0007] In order to achieve the above object, the present invention adopts the following technical solutions:

[0008] In one aspect, the present invention discloses a method for quantifying terrace pattern evolution and its ecosystem service trade-off, comprising:

[0009] Build a multi-source remote sensing information coordination framework, identify frequently changing areas, and pre-process multi-source remote sensing data to obtain a standard input dataset through temporal and spatial benchmark unification, key feature enhancement, redundant information compression, and dynamic quality feedback optimization.

[0010] Based on a standard input dataset, a multimodal boundary refinement segmentation algorithm is used to extract quantitative data on the distribution of terraces in typical regions and their interannual expansion / contraction dynamics to characterize their spatiotemporal evolution. By integrating multimodal deep learning with spatial statistical models, a comprehensive quantitative analysis of the terrace evolution process is achieved.

[0011] Based on the extracted terrace distribution data, combined with standard input datasets and ground observation data, an evaluation framework is proposed to quantitatively characterize the terraces in soil conservation services, water conservation services, carbon sequestration services, food supply services and biomass energy supply services, and to conduct a synergistic analysis of the trade-offs among the above typical ecosystem services.

[0012] Preferably, in the above-mentioned method for quantifying terrace pattern evolution and its ecosystem service trade-off, the specific steps of constructing a multi-source remote sensing information collaborative framework are as follows:

[0013] A dynamic weight allocation algorithm is used to achieve a global optimal balance between spatial resolution, coverage efficiency, and time cost through mathematical modeling. The Herfindahl index is also introduced to measure data diversity and redundancy to ensure that weight allocation is comprehensive while avoiding information overload. When necessary, adaptive grid division is performed on the study area to ensure smooth image data integration.

[0014] Construct a "macro-representation-process capture-micro-calibration" framework to carry out complementary fusion of multi-source data;

[0015] A mesoscale adaptive data collection command feedback mechanism is used to improve data collection efficiency and identify frequently changing areas.

[0016] Preferably, in the above-mentioned method for quantifying terrace pattern evolution and its ecosystem service trade-off, the weight calculation model of the dynamic weight allocation algorithm is as follows:

[0017]

[0018] Among them, W i Represents the priority weight of data source i, which determines its collection order in scheduling and is dimensionless; SR i Indicates the spatial resolution of data source i, in meters; CR i represents the coverage efficiency of data source i, that is, the coverage area per unit time, in km2 / day; λ0 represents the time-effectiveness attenuation coefficient; Δt represents the time delay after data collection, in days; HI(D similar) is the Herfindahl index, which measures the redundancy of similar data sources, where s j is the ratio of the coverage area of ​​data source j to the total overlapping area; k represents the total number of data sources.

[0019] Preferably, in the above-mentioned method for quantifying terrace pattern evolution and its ecosystem service trade-off, the specific steps for extracting terrace distribution data for a typical area are as follows:

[0020] Define optical image features, terrain parameter features, and vegetation time series features as input features;

[0021] Constructing a dynamic weighted multimodal segmentation network based on the input features to obtain a segmentation probability map of the terraced area;

[0022] The segmentation probability map is subjected to threshold segmentation, and concave points are detected on the terrace patch contour, the concave points are connected, and the obtained binary mask is converted into a polygonal vector layer, and the final terrace distribution vector map is output.

[0023] Preferably, in the above-mentioned method for quantifying terrace pattern evolution and its ecosystem service trade-off, the specific steps of constructing a dynamic weighted multimodal segmentation network are as follows:

[0024] To build MMF-SegNet, we first extract branch features. In the optical branch, we use a dilated residual module to capture multi-scale textures and output an optical feature map. Secondly, in the terrain branch, we input two-channel data consisting of slope and terrain curvature, encode it using a lightweight U-Net, and output terrain features. Finally, for time series analysis, we input NDRE time series and use a 1D Conv-LSTM to extract time-dependent features, outputting a vector. Finally, for feature module fusion, we introduce terrain-guided attention and dynamically fuse optical and terrain features to obtain fused features.

[0025] F fusion =α·F RGB +(1-α)·F DEM ;

[0026] α=σ(RI·W);

[0027]

[0028] Where RI represents the terrain roughness index, which is a composite of slope and terrain curvature; W refers to the learnable parameter, which is optimized by back propagation; F RGB Represents the optical image characteristics, F DEM represents the terrain parameter characteristics; σ is the Sigmoid function; the weight coefficient α is dynamically adjusted by the terrain complexity; θ is the slope angle calculated based on the DEM, reflecting the degree of surface inclination; C represents the terrain curvature;

[0029] Finally, the vegetation temporal feature FNDRE and fusion feature F fusion Splicing, through convolution fusion, outputs the segmentation probability map of the terraced area.

[0030] Preferably, in the above-mentioned method for quantifying terrace pattern evolution and its ecosystem service trade-off, the specific steps for quantifying the spatiotemporal evolution characteristics of terraces and extracting the dynamic representation of interannual expansion / contraction are as follows:

[0031] Calculate the rate of change in terrace area, count the number of expanded / reduced terrace patches, their spatial distribution, and area proportion;

[0032] Evolution of spatial pattern index, the selected spatial pattern indexes include fragmentation index, morphological stability index, center of gravity migration trajectory, and interannual migration distance;

[0033] To generate the spatiotemporal evolution map of terraces, we first generate a density map of the expansion / contraction of terraces year by year in the form of a heat map sequence; then we generate a change trajectory map, and construct a change path network diagram G = (V, E), where the node V is the centroid of the change patch, the edge E connects the centroids of adjacent years, and the weight is the migration distance; finally, we construct a three-dimensional space-time cube model, stacking the two-dimensional masks between years according to the pixels, and display the evolution trend of the vertical time axis through voxel rendering to generate a visual map.

[0034] Preferably, in the above-mentioned method for quantifying terrace pattern evolution and its ecosystem service trade-off, the calculation of terrace area change rate, the statistics of the number of expanded / contracted terrace patches and their spatial distribution and area proportion are carried out in the following specific steps:

[0035] The specific formula for calculating the rate of change of terrace area is as follows:

[0036] ΔA t =A t+1 -A t ;

[0037]

[0038] In the formula, A t represents the total area of ​​terraces in the study area in year t, S k represents the area of ​​the kth terrace patch, K represents the total number of terrace patches; A t+1 represents the total area of ​​terraces in the study area in year t+1, in square meters or hectares; ΔA t Represents the interannual variation of terrace area; if ΔA t > 0 and has been growing for three consecutive years, marked as a stable expansion zone. t <-0.2A t , determined to be a damage event;

[0039] Perform spatial overlay analysis and detect the changed area pixel by pixel. The formula is as follows:

[0040]

[0041] In the formula, C t (i, j) represents the category label of pixel (i, j) from time t to t+1; M t (i, j) represents the binary mask of year t, 0 represents non-terraced fields and 1 represents terraced fields.

[0042] Preferably, in the above-mentioned method for quantifying terrace pattern evolution and its ecosystem service trade-off, the specific steps of the evaluation framework for quantitatively characterizing terrace ecosystem functions are as follows:

[0043] Quantification of soil conservation services: The improved RUSLE model is used, and the specific formula is as follows:

[0044] A actual =R·K·L·S·C terrace ·P terrae ;

[0045] A potential =R·K·L·S;

[0046] Soil Conservation=A potential -A actual ;

[0047] In the formula, A potential and A actual They represent the potential soil erosion and actual soil erosion in the region, and the difference between them represents the amount of soil retained by terrace engineering in the region; R represents the rainfall erosivity (MJ·mm / (hm 2 h·yr), the total annual rainstorm erosivity is calculated based on local rain gauges or CHIRPS remote sensing rainfall data; K represents soil erodibility in t·hm 2 ·h / (hm 2 ·MJ·mm), which was parameterized by the EPIC model using the soil texture map; L·S represents the topographic factor, and the slope and slope length are generated based on the satellite elevation data of the terrace distribution area. The formula is:

[0048]

[0049] S=10.8sinθ+0.03, θ<5°;

[0050] S=16.8sinθ-0.5,θ≥5°;

[0051] In the formula, λ is the slope length, unit: meter. After extracting the watershed watershed line, the actual slope length is calculated by dividing the fields into terraces; m is the slope length factor index; θ is the slope angle of the slope; C terrace Represents the terrace vegetation cover management factor, dimensionless, based on the dynamic inversion of the red edge band:

[0052]

[0053] NDRE max represents the maximum NDRE during the crop growing season; P terrace represents the terrace engineering factor, dimensionless;

[0054] Water conservation service assessment: A distributed water conservation model for terraces was constructed, coupling the spatial heterogeneity of terrace structure and soil hydraulic parameters. The specific formula is as follows:

[0055]

[0056] ET i =ET natural ·(1-γ·NDWI wet );

[0057]

[0058]

[0059] In the formula, P i represents the annual precipitation of grid cell i, in mm, obtained through interpolation of weather station data or GSMaP satellite precipitation products; ET i represents the actual evapotranspiration of the terrace in grid unit i, in mm, corrected by the water-saving effect of the terrace; R i Indicates the surface runoff depth of the terrace in grid cell i, unit: mm; A i represents the area of ​​grid cell i; η i represents the soil infiltration enhancement coefficient, dimensionless, dynamically corrected based on soil type and terraced tillage layer thickness; ETnatural represents the evapotranspiration of natural vegetation, mm; where γ represents the interception coefficient of terraced fields in the wet season, calibrated by ground infiltration experiments, and NDWI wet is the canopy moisture index of crops in the wet season; it is calculated using the SCS-CN model: CN value represents the runoff curve number, dimensionless; P represents annual precipitation, mm; S max Indicates the potential maximum water retention, mm;

[0060] Carbon sequestration service assessment: Construct a hierarchical carbon sink model. First, the vegetation carbon sequestration amount. The specific formula is as follows:

[0061] NPP terrace =APAR·ε·f crop ;

[0062] In the formula, APAR represents absorbed photosynthetically active radiation, MJ / m 2 ;ε represents the light energy utilization rate, gC / MJ; f crop is the crop rotation coefficient; the second is the soil carbon sequestration amount, and the specific calculation formula is as follows:

[0063]

[0064] In the formula, C input represents the carbon input of crop litter, kg / ha / yr; SOC initial represents the initial soil organic carbon content, %; α and β represent the organic carbon conversion and mineralization rates, which are assigned according to the tillage method; BD represents the soil bulk density, g / cm 3 ; D represents the depth of the tillage layer, m;

[0065] Food supply service assessment: Using the mixed yield model, first calculate the potential yield. The specific formula is as follows:

[0066] Y potential =NPP·HI·CF;

[0067] In the formula, NPP represents net primary productivity, gC / m2 / year; HI represents harvest index, dimensionless; CF represents the carbon content-yield conversion factor, kg / gC; combined with actual yield correction, the specific formula is as follows:

[0068] Y actual =Y potential ·(1-δ drought -δ erosion );

[0069] In the formula, δ drought represents the drought stress loss rate, which is inverted by the vegetation temperature condition index; δ erosion It represents the yield reduction rate caused by soil erosion and is negatively correlated with soil conservation;

[0070] Biomass energy supply service assessment: Utilize the multi-source biomass inversion model to calculate the energy content of relevant crop straw. The specific formula is as follows:

[0071]

[0072] In the formula, Yi represents the yield of the i-th type of crop, kg / ha, and n represents the number of crop types; RPR i represents the grass-to-grain ratio; η i represents the energy conversion efficiency; then the biodiesel potential of oil crops is calculated, the specific formula is as follows:

[0073] Ebiodiesel =ρ·A·Y oil ·η transester ;

[0074] In the formula, ρ represents the ratio of crop planting area to terrace area, dimensionless; A is the total area of ​​terrace, ha; Y oil Indicates oil yield per unit area, L / ha; η transester It represents the conversion rate of transesterification reaction.

[0075] Preferably, in the above-mentioned method for quantifying terrace pattern evolution and its ecosystem service trade-off, the specific steps for conducting a collaborative analysis of trade-offs between ecosystems are as follows:

[0076] Construct a framework for analyzing trade-off synergy relationships: construct a foundational layer and conduct statistical correlation analysis to analyze multi-scale trade-off synergy relationships;

[0077] Core layer construction, machine learning causal identification, and a two-stage causal-mechanism decoupling model are proposed: Random forest feature importance ranking is used to screen the main controlling factors that affect the trade-offs between services; based on a dual machine learning framework, the causal effects of natural factors and human activities are separated;

[0078] Output layer construction, trade-off synergy index calculation, and definition of spatiotemporal dynamic trade-off index;

[0079] The application layer is constructed, spatial explicit collaborative optimization is performed, and a multi-objective particle swarm optimization model is built. The objective function is to maximize the total synergistic gain of five ecosystem services. The constraints include the non-decreasing threshold of ecosystem services and the upper limit of terrace engineering transformation.

[0080] Through the above technical solutions, it can be seen that compared with the existing technology, the present invention provides a method for quantifying the evolution of terrace patterns and weighing their ecosystem services, which overcomes the problems of insufficient spatiotemporal resolution, static ecosystem service assessment and weak decision-making support in the existing technology; the present invention proposes a dynamic balance model of spatial resolution-coverage efficiency-timeliness-redundancy to solve the nonlinear contradiction of multispectral / high-resolution data collaborative scheduling, and the scheduling efficiency is improved by 40% compared with traditional methods; the MMF-SegNet multimodal segmentation network driven by the terrain-guided attention mechanism (TGA) is developed, and the terrace boundary recognition accuracy (ES) is improved from 0.74 to 0.91 by fusing optical texture temporal features with LiDAR terrain curvature information; breaking through the limitations of traditional static models, the red edge band remote sensing inversion parameters are introduced. Dynamic assignment technology (random forest dynamic modeling of vegetation management factor C) constructs a spatial heterogeneity adaptive assessment framework for soil conservation and water conservation services, reducing verification errors by 15%-22%; independently designs spatiotemporal explicit trade-off collaborative analysis technology (STTI), integrating dual machine learning causal inference and MOPSO multi-objective optimization algorithm to analyze the natural-human driving relationship in terrace multi-ecosystem services and support the generation of differentiated regulatory rules; through the full-chain logical closed loop of "dynamic data coupling-intelligent feature mining-heterogeneity assessment-collaborative optimization", it not only overcomes the shortcomings of traditional methods in mesoscale monitoring, such as model rigidity and service association fragmentation, but also significantly reduces R&D risks and deployment thresholds with the help of compatibility adaptation of mature technologies, providing a high-precision and scalable toolset for the sustainable management of terrace ecosystems. BRIEF DESCRIPTION OF THE DRAWINGS

[0081] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following briefly introduces the drawings required for use in the embodiments or the description of the prior art. Obviously, the drawings described below are merely embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on the provided drawings without paying any creative work.

[0082] Figure 1 The accompanying drawing is a flow chart of the method provided by the present invention. DETAILED DESCRIPTION

[0083] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying 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 making creative efforts are within the scope of protection of the present invention.

[0084] like Figure 1As shown, the embodiment of the present invention discloses a method for quantifying the evolution of terraced fields and its ecosystem service trade-off, including:

[0085] S101 builds a collaborative framework for multi-source remote sensing information, identifies frequently changing areas, and pre-processes multi-source remote sensing data to obtain a standard input dataset through spatial and temporal benchmark unification, key feature enhancement, redundant information compression, and quality dynamic feedback optimization.

[0086] Based on a standard input dataset, S102 uses a multimodal boundary refinement segmentation algorithm to extract quantitative data on the distribution of terraces in typical regions and their interannual expansion / contraction dynamics to characterize their spatiotemporal evolution. By integrating multimodal deep learning with spatial statistical models, it achieves a comprehensive quantitative analysis of the terrace evolution process.

[0087] S103 proposes an evaluation framework based on the extracted terrace distribution data, combined with standard input data sets and ground observation data, to quantitatively characterize the terraces in soil conservation services, water conservation services, carbon sequestration services, food supply services and biomass energy supply services, and conducts a synergistic analysis of the trade-offs among the above typical ecosystem services.

[0088] Specifically, in S101, during the data collection phase, in order to obtain comprehensive and accurate mesoscale terrace distribution information, the embodiment of the present invention constructs a "sky-air-ground" multi-source remote sensing information collaboration framework. The specific solution is as follows:

[0089] S1011 Dynamic data weight allocation algorithm based on mesoscale sensitivity: In mesoscale terrace monitoring, the existing data collection has the contradiction between "high-resolution local coverage" and "low-resolution global observation". For example, the coverage of high-resolution UAV data is limited (1km 2 / flights), while the satellite revisit cycle for the entire region is long (e.g., 7 days), resulting in low monitoring efficiency. To address this problem, this embodiment proposes a dynamic weight allocation algorithm. The goal is to achieve a global optimal balance between spatial resolution, coverage efficiency, and time cost through mathematical modeling. At the same time, the Herfindahl Index (HI) is introduced to measure data diversity and redundancy to ensure that the weight allocation is comprehensive while avoiding information overload. The weight calculation model is as follows:

[0090]

[0091] Among them, W i The priority weight of data source i determines its collection order in scheduling. It is dimensionless and ranges from [0,1]. i Indicates the spatial resolution of data source i. The smaller the value, the higher the accuracy. The unit is meter. i Indicates the coverage efficiency of data source i (i.e., coverage area per unit time), in km2 / day; λ0 represents the timeliness attenuation coefficient, which is used to control the negative impact of the newness of the data on the weight. This parameter can be dynamically adjusted according to the specific conditions of the region, and its gradient is defined for different terrain complexities (based on surface roughness). If the surface roughness in the region is lower than 0.05 (plain area), λ0 = 0.15; if the surface roughness is between 0.05 and 0.3 (hilly area), λ0 = 0.2; if it is higher than 0.3 (mountainous area), λ0 = 0.25; Δt represents the time delay after data collection, in days; HI(D similar ) is the Herfindahl index, which measures the redundancy of similar data sources and has a value range of [0,1]. j is the ratio of the coverage area of ​​data source j to the total overlapping area.

[0092] In addition, the image data provided by the ground-based UAV equipment should also be prepared to be connected with the space-based satellite image data. Therefore, the study area is adaptively gridded, and the grid size is divided according to the following formula:

[0093]

[0094] Where G represents the adaptive grid size, which is the minimum unit for controlling data fusion and is expressed in km; CR min Indicates the minimum data source coverage efficiency requirement, in km 2 .

[0095] S1012 Complementary fusion of multi-source data: Different types of sensors have blind spots in feature expression in terraced field recognition. For example, optical satellites are susceptible to cloud and fog interference, drones lack terrain penetration, and LiDAR is expensive. The present invention uses a multi-level fusion strategy to build a data chain of "macro-representation-process capture-micro-calibration" to overcome the perception limitations of a single data source. See Table 1 for details:

[0096] Table 1 Parameters of different data sources, collection indicators and data fusion algorithms

[0097]

[0098] The slope error correction between lidar data and satellite elevation data adopts the residual least square method. The specific formula is as follows:

[0099]

[0100] Where θ is the slope and σ is the standard deviation.

[0101] The selection logic of the above technologies is:

[0102] Space-based WorldView-3 red edge band: Utilizing the sensitivity of the band near 750nm to vegetation biochemical parameters, it enhances the spectral distinction between terraced ridges and crops, overcoming the fuzzy ridges in traditional RGB images of terraced fields.

[0103] Airborne multispectral drone time series interpolation: To address the gaps in cloud cover caused by the satellite revisit cycle (for example, there are only 5 to 7 days of valid data per month during the rainy season), vegetation growth curves are fitted based on Fourier transform to ensure monitoring continuity.

[0104] Ground-based LiDAR point cloud correction: Because satellite optical images have shadow distortion in steep slope areas (>25°), LiDAR point clouds are introduced to generate sub-meter DEM through surface reconstruction to eliminate slope inversion errors.

[0105] S1013 mesoscale adaptive data acquisition command feedback mechanism:

[0106] The dynamic changes in terraced fields (such as heavy rains destroying ridges and tillage leading to the merging of fields) are spatiotemporally heterogeneous. Traditional regular inspections result in a large amount of ineffective data collection (approximately 60% of inspections show no significant changes). By building an event-driven feedback mechanism, a closed-loop control system of "data collection-change identification-task triggering" is implemented, focusing limited resources on areas of high change. The fragmentation index used here reflects the degree of fragmentation of terraced fields. The specific formula is as follows:

[0107]

[0108] Where Area(P i ) is the area of ​​a single terraced field in square meters; P total is the total area of ​​the region, in square meters; when FI is higher than 0.7, it can be determined as a highly fragmented terraced area.

[0109] The priority scheduling of the drone supplementary mission triggering mechanism uses the Hungarian algorithm to match high-priority areas with idle drones. The specific formula is as follows:

[0110] Cost ij =α·ΔFI j +β·d ij (6);

[0111] Where ΔFI j represents the rate of change of the fragmentation index of region j, which reflects the urgency of monitoring and has a value range of [0,1); d ij represents the flight distance from the current position of UAV i to region j, in kilometers; α and β represent the weight coefficients of the cost function, satisfying α + β = 1, and their value range is [0, 1].

[0112] The main purpose of S101 is to resolve the contradiction between "high resolution but low coverage" and "low resolution but wide coverage" in mesoscale terrace monitoring, to maximize data acquisition efficiency (balance of coverage, timeliness and quality) and complementary use of multi-source data (breaking the limitations of a single sensor), and to provide a high-quality input basis for subsequent analysis.

[0113] In S101, multi-source remote sensing data is dynamically preprocessed:

[0114] To meet the demand for efficient collection and integration of multi-source data from "sky, air, and ground", a heterogeneous data pre-standardization technology for dynamic monitoring of terrace patterns is proposed. The goal is to resolve modeling errors caused by data source complexity and spatiotemporal heterogeneity through spatiotemporal consistency correction, redundant information compression, and unified representation of gradient terrain features. The specific content is as follows:

[0115] (1) Unification of time and space benchmarks and geometric correction:

[0116] Because the spatial references and resolutions of multi-source data collection platforms (satellites, drones, and LiDAR) vary significantly, they require processing before use. The first step is standardization of data coordinates and resolution. All data must be unified to the same coordinate system (WGS84 UTM) and spatial resolution (mesoscale core resolution 0.5m to 5m). A dynamic resolution compensation algorithm is used to avoid information loss. The specific formula is as follows:

[0117]

[0118] Where R i is the resolution of the data, W i The mixed-resolution data is weighted and fused according to its acquisition weight. Secondly, the data is geometrically corrected. For drone imagery, SIFT feature matching combined with affine transformation (see Table 1) is used to eliminate stitching offsets caused by unstable flight attitude. For LiDAR and satellite elevation data, a residual correction model is established to eliminate local slope inversion errors.

[0119] (2) Multi-source data feature enhancement:

[0120] Focusing on issues such as blurred terrace boundaries and terrain obstruction, key information features are enhanced by combining data source characteristics. First, spectral cloud removal and radiometric normalization are performed. For satellite optical imagery (such as WorldView-3), a multi-temporal cloud mask synthesis method is used to dynamically remove invalid pixels based on the Normalized Difference Cloud Index (NDCI). The specific formula is as follows:

[0121]

[0122] Where ρ NIRRepresents the near-infrared band reflectance (dimensionless, ranging from [0,1]), which is used to detect the difference between vegetation and cloud; ρ NIR NDCI represents the shortwave infrared reflectance (dimensionless, ranging from [0,1]) and is sensitive to clouds and fog. If NDCI > 0.2, the pixel is considered to be effectively cloud-free (empirical threshold). Secondly, the UAV multispectral data is corrected for radiance using a solar altitude model to eliminate illumination variations. Finally, to enhance terrain features, LiDAR point cloud processing uses the Poisson surface reconstruction algorithm to generate a 0.1m precision DEM (see Table 1). Shadow compensation is performed on steep slopes (>25°). After slope correction using Equation 4, rasterized terrain parameters (slope, roughness) are output. Superpixel segmentation (SLIC) combined with Canny edge detection is used to enhance the geometric features of terraced ridges and mitigate interference from crop cover.

[0123] (3) Redundant information compression and dynamic quality assessment:

[0124] Redundancy elimination was previously performed during the data collection phase, but for the redundancy problem of repeated coverage areas of multi-source data, an optimization strategy needs to be proposed. The first is dynamic clipping based on the Herfindahl Index (HI). Redundant data is automatically eliminated based on the redundancy index defined in Formula 2. The threshold is set to eliminate similar low-weight data when the HI index is higher than 0.6, reducing storage and subsequent calculations. Secondly, adaptive grid quality control is used. The grid size in Formula 3 is used to check the data integrity (coverage efficiency CR) of each unit in blocks. i ≥CR min ), triggering drone re-collection instructions for insufficient areas (see Formula 6). Finally, the data priority is dynamically updated. Based on the pre-processed data quality feedback (such as registration error and signal-to-noise ratio), the weight parameter λ0 in Formula 1 is reversely optimized, forming a "collection-preprocessing-optimization" closed loop (for example, when the registration error in the plain area exceeds the threshold, λ0 is increased from 0.15 to 0.18 to accelerate the elimination of old data).

[0125] The above technical solution aims to solve the modeling error problem caused by complex sources, resolution differences, terrain interference and redundant coverage of multi-source remote sensing data through the unification of spatiotemporal benchmarks, enhancement of key features, compression of redundant information and optimization of dynamic quality feedback, and generate a standard input data set that is consistent in time and space, has enhanced features and controllable redundancy, laying a high-precision, low-noise data foundation for the dynamic extraction and evolution analysis of terrace patterns.

[0126] In S102, the multimodal boundary refinement segmentation algorithm is used to extract terrace distribution data for a typical area. The specific steps are as follows:

[0127] By integrating optical texture, terrain curvature, and temporal vegetation characteristics, we can solve the problems of weak edges (crop cover) and fragmentation (shadow interference) of ridges and generate continuous vector boundaries of terrace patches, as shown below:

[0128] Step 1: Input feature definition and channel processing;

[0129] The input data consists of three types of features. The first is optical image features, including the RGB bands and the red edge band. The time phase is selected during the period of low crop cover (such as early spring when no crops are cultivated) to maximize the exposure of terraced field ridge structures. The second is terrain parameter features, including slope and terrain curvature. The specific formula for slope is as follows:

[0130]

[0131] Where θ is the slope angle calculated based on the DEM, reflecting the degree of surface inclination, and the range is [0°, 90°]; z represents the elevation value (unit: meter), which is derived from the digital elevation model; and are the elevation gradients of the DEM in the x and y directions respectively; the terrain curvature is used to enhance the contour features of the terraced ridges. The specific formula is as follows:

[0132]

[0133] Where C represents the terrain curvature (unit: curvature per unit length, usually 1 / meter). Positive curvature corresponds to raised terraced ridges, while negative curvature corresponds to valleys. and is the second-order derivative of DEM in the x and y directions, The mixed second-order derivative of DEM is calculated by 3×3 window difference. Finally, the vegetation time series characteristics are the interannual NDRE (Normalized Difference Red Edge) time series, with a weekly interval of 16 days. The specific calculation formula is as follows:

[0134]

[0135] Where ρ NIR Represents the reflectivity in the near-infrared band (wavelength of about 770-900nm). Strong reflection areas usually correspond to vegetation canopies; ρ Red Edge Represents the red-edge reflectance (wavelength approximately 700-740 nm), is sensitive to chlorophyll content, and is used to distinguish vegetation types; the NDRE value range is [-1, 1]. Higher values ​​indicate denser vegetation cover or higher chlorophyll content.

[0136] Step 2: Dynamic weight multimodal segmentation network;

[0137] Design MMF-SegNet (Multi-Model Fusion Segmentation Network), first extract branch features, in the optical branch, use the dilated residual block to capture multi-scale texture, the convolution kernel dilation rate is set to 1, 2, 4, the output feature map F RGB Secondly, the terrain branch inputs the 2-channel data of slope and terrain curvature, uses lightweight U-Net encoding, and outputs terrain features F DEM Finally, for time series analysis, we input the NDRE time series (time step T = 24), use 1DConv-LSTM to extract time-dependent features, and output the vector FNDRE. Next, we perform feature module fusion, introducing Terrain-Guided Attention (TGA) to dynamically fuse optical and terrain features:

[0138] F fusion =α·F RGB +(1-α)·F DEM (12);

[0139] α=σ(RI·W) (13);

[0140]

[0141] Where RI represents the terrain roughness index, which is a composite of slope and terrain curvature; W refers to a learnable parameter, initialized to 0.8 and optimized through backpropagation; σ represents the Sigmoid function, which constrains the weight to [0, 1]; the weight coefficient α is dynamically adjusted by the terrain complexity; high RI areas (steep and broken terrain) reduce the optical weight (α) to enhance terrain features and suppress shadow interference. Finally, the temporal feature is inserted to convert the NDRE feature F NDRE and fusion feature F fusion Splicing, through 3×3 convolution fusion, outputs the segmentation probability map P∈[0,1] of the terraced area.

[0142] Step 3: Post-processing and vectorization;

[0143] This step includes two parts. The first is morphological optimization, which performs threshold segmentation on the probability map P with a threshold of τ = 0.7. The post-processing process first performs an opening operation and a structural operation B1 = 3×3 rectangle to eliminate small noise; the second is concave point connection, which detects concave points on the terrace patch contour. If the concave point depth d ≥ 2 pixels, B-spline interpolation is used to smoothly connect the broken boundaries. Vectorization requires converting the binary mask into a polygonal vector layer, filtering out pseudo-fields with an area less than 500m2 (excluding the actual minimum size of the terrace), and outputting the final terrace distribution vector map M.t (Year t).

[0144] In S102, the temporal and spatial evolution characteristics of terraces are quantified to extract dynamic representations of interannual expansion / contraction. The goal is to extract dynamic indicators such as the spatial expansion / contraction trend, morphological fragmentation rate, and stability of terraces based on the interannual terrace vector sequence. The specific indicators are as follows:

[0145] Step 1: Detection of interannual changes and event classification;

[0146] First, the terrace area change rate is calculated. The specific formula is as follows:

[0147] ΔA t =A t+1 -A t (15);

[0148]

[0149] In the formula, A t represents the total area of ​​terraces in the study area in year t, S k represents the area of ​​the kth terrace patch, K represents the total number of terrace patches; A t+1 represents the total area of ​​terraces in the study area in year t+1 (unit: square meters or hectares); ΔA t Represents the inter-annual change in the area of ​​terraced fields. t > 0 and has been growing for three consecutive years, marked as a stable expansion zone. t <-0.2A t (i.e., a reduction of more than 20% in a single year) is considered a damage event. Then, a spatial overlay analysis is performed to detect the changed area pixel by pixel. The formula is as follows:

[0150]

[0151] In the formula, C t (i, j) represents the category label of pixel (i, j) from time t to t+1; M t (i, j) represents the binary mask of year t, 0 represents non-terraced fields and 1 represents terraced fields.

[0152] Step 2: Analysis of spatial pattern index evolution;

[0153] First, select the appropriate spatial pattern index, as follows:

[0154] Fragmentation Index (FI):

[0155]

[0156] N in the formula tIndicates the total number of terrace patches; represents the average plaque area; FI t The larger the value, the more fragmented the terraced landscape is and the greater the degree of disturbance.

[0157] Shape Stability Index (SSI):

[0158]

[0159] Perimeter in the formula k Area represents the perimeter of the kth terrace patch (unit: meter); k represents the area of ​​the kth terrace patch (unit: square meters); the denominator is the perimeter of the patch when it is assumed to be circular (idealized shape); the closer the SSI value is to 1, the closer the patch shape is to a circle (regular edges); if it is greater than 1, the shape is complex (irregular or elongated).

[0160] Center of gravity migration trajectory:

[0161] Coordinates of the centroid of terrace distribution in each year (Xt, Yt):

[0162]

[0163] In the formula, S k represents the area of ​​the kth terrace patch; (x k ,y k ) represents the centroid coordinates of the kth patch (based on the geometric center of the vector polygon); (X t , Y t ) represents the coordinates of the total centroid of the terrace distribution in year t.

[0164] Interannual migration distance:

[0165] D in the formula t Indicates the migration distance of the center of gravity from year t to year t+1 (the unit is consistent with the coordinates, meters or kilometers).

[0166] Step 3: Generation of spatiotemporal evolution map;

[0167] This step primarily generates visualization maps. First, a density map of the terrace expansion and contraction over time is generated as a series of heat maps. A kernel density estimation bandwidth of 1 km is used to depict spatial clustering. Next, a change trajectory map is constructed, constructing a change path network (G = (V, E)). Node V represents the centroid of a change patch, and edges E connect the centroids of adjacent years, with weights corresponding to migration distances. Finally, a three-dimensional space-time cube model is constructed. An inter-annual two-dimensional mask, Mt, is stacked pixel by pixel, and voxel rendering is used to visualize the evolutionary trend along the vertical time axis.

[0168] In addition, technical verification and accuracy evaluation are mainly performed on the above technologies and generated results. First, the accuracy of the segmentation results is evaluated. Here, the confusion matrix is ​​calculated based on the verification sample, and the F1-Score is used to balance the accuracy (Precision = 0.89, ). In addition, the edge similarity (ES) indicator is used to measure the spatial alignment accuracy between the predicted boundary and the true boundary. The closer the value is to 1, the more accurate the edge positioning is. The specific formula is as follows:

[0169]

[0170] In the formula, EDT(P) represents the Euclidean distance transform map of the predicted result P, and each pixel value represents the distance to the nearest predicted edge; EDT(G) represents the Euclidean distance transform map of the true annotation G.

[0171] Next is the verification of the spatiotemporal evolution analysis results. This part mainly verifies the accuracy of change event detection through historical Google Earth image sampling, and uses an error matrix to reflect the accuracy of the results.

[0172] The above technical solution mainly focuses on the precise extraction of terrace distribution and the characterization of its spatiotemporal dynamics. By linking multimodal deep learning with spatial statistical models, it achieves a comprehensive quantitative analysis of the terrace evolution process.

[0173] In S103, based on previously extracted terrace distribution data (vector boundaries, spatiotemporal change trajectories), combined with multi-source remote sensing and ground observation data, an assessment framework of "spatiotemporal heterogeneity modeling - key process coupling - multi-dimensional service collaboration" was proposed to quantitatively characterize the terraces' ability to quantify soil conservation services, assess water conservation services, assess carbon sequestration services, assess food supply services, and assess biomass energy supply services. This framework emphasizes the dynamicization of model parameters, spatial heterogeneity responses, and multi-service collaborative optimization. The specific content is as follows:

[0174] (1) Quantification of soil conservation services;

[0175] The goal is to evaluate the ability of terraces to control soil erosion and quantify their contribution to reducing sediment loss. The specific technical solution is to use an improved RUSLE model (revised soil and water conservation equation). Compared with traditional methods, this model uses the red edge band to dynamically invert vegetation cover and terrace engineering factors, improving the spatial resolution of the parameters. The specific formula is as follows:

[0176] A actual =R·K·L·S·C terrace ·P terrace ;

[0177] A potential=R·K·L·S (24);

[0178] Soil Conservation=A potential -A actual (25);

[0179] A potential and A actual They represent the potential soil erosion and actual soil erosion in the region, and the difference between them represents the amount of soil retained by terrace engineering in the region; R represents the rainfall erosivity (MJ·mm / (hm 2 h·yr), the total annual rainstorm erosivity is calculated based on local rain gauges or CHIRPS remote sensing rainfall data; K represents soil erodibility in t·hm 2 ·h / (hm 2 ·MJ·mm), which was parameterized by the EPIC model using the soil texture map; L·S represents the topographic factor, and the slope and slope length are generated based on the satellite elevation data of the terrace distribution area. The formula is:

[0180]

[0181] S=10.8sinθ+0.03, θ<5° (27);

[0182] S=16.8sinθ-0.5, θ≥5° (28);

[0183] In the formula, λ is the slope length (unit: meter). The actual slope length is calculated by extracting the watershed line from the DEM and dividing the watershed into terraced fields; m is the slope length factor index; θ is the slope angle of the slope; C terrace Represents the terrace vegetation cover management factor (dimensionless), based on the dynamic inversion of the red edge band (linked with time series NDRE):

[0184]

[0185] NDRE max It represents the maximum NDRE value during the crop growing season, with a value range of [0,1], reflecting the protective effect of vegetation cover on the soil; P terrace It represents the terrace engineering factor (dimensionless), and the assignment rule is: 0.1 for horizontal terraces, 0.15 for reverse slope terraces, and 0.2 to 0.3 for earthen ridge terraces (assigned based on ridge height through drone image classification).

[0186] (2) Water conservation service assessment;

[0187] The goal is to assess the ability of terraces to regulate runoff and enhance groundwater recharge. The technical solution is to construct a terraced distributed water conservation model (TWHM), coupling the spatial heterogeneity of terrace structure and soil hydraulic parameters. The specific formula is as follows:

[0188]

[0189] In the formula, P i represents the annual precipitation of grid cell i (mm), which is obtained through interpolation of weather station data or GSMaP satellite precipitation products; A i represents the area of ​​grid cell i; ET i represents the actual evapotranspiration of the terrace in grid cell i, (mm), based on the inversion of the SEBAL model or MODIS ET product (MOD16), combined with the correction of the water-saving effect of the terrace:

[0190]

[0191] In the formula, γ represents the interception coefficient of terraced fields in the wet season (calibrated by ground seepage test), NDWI wet is the crop canopy moisture index in the wet season; R i represents the surface runoff depth of the terrace within grid cell i (unit: mm), calculated using the SCS-CN model:

[0192]

[0193] CN value adjustment: The CN value of the terraced area is reduced to 75-85 (depending on the water storage capacity of the terraced ridges); η i It represents the soil infiltration enhancement coefficient, which is dynamically corrected based on soil type (e.g., η = 1.2 for clay loam and η = 0.8 for sandy soil) and the thickness of the terraced tillage layer (derived from LiDAR inversion).

[0194] (3) Carbon sequestration service assessment;

[0195] The goal is to measure the amount of photosynthetic carbon sequestration and soil organic carbon accumulation in terraced crops. The specific technical solution is to construct a hierarchical carbon sink model. The first part is the amount of carbon sequestration by vegetation. The specific formula is as follows:

[0196] NPP terrace =APAR·ε·f crop (34);

[0197] In the formula, APAR (MJ / m 2 ) represents absorbed photosynthetically active radiation, which is calculated by inverting FPAR (vegetation coverage) from the red edge band + solar radiation data; ε represents light energy utilization efficiency (gC / MJ), which is assigned a value based on the crop type (1.8 for rice, 2.1 for corn, and 1.5 for wheat); f crop is the crop rotation coefficient (e.g. double-season rice f = 1.5). The second part is the soil carbon sequestration capacity, which is calculated as follows:

[0198]

[0199] In the formula, C input represents the carbon input of crop litter (kg / ha / yr), which is inverted by biomass (NPP) and root-to-shoot ratio; SOC initial represents the initial soil organic carbon content (%), based on the inversion of the Sentinel-2 shortwave infrared band; α and β represent the organic carbon transformation and mineralization rates, which are assigned according to the tillage method (e.g., no-tillage α = 0.6); BD represents the soil bulk density (g / cm 3 ), which is spatially interpolated from measured soil profile data; D represents the depth of the tillage layer (m), which is obtained by combining DEM with remote sensing classification results of tillage activities.

[0200] (4) Food supply service assessment;

[0201] The goal is to measure the grain productivity per unit area of ​​terraced fields and its spatial differentiation. The technical solution adopted is to use a mixed yield model to first calculate the potential yield. The specific formula is as follows:

[0202] Y potential =NPP·HI·CF (36);

[0203] In the formula, HI represents the harvest index (crop economic allocation ratio, such as rice HI = 0.45); CF represents the carbon content-yield conversion coefficient (rice CF = 0.4). The second part is to correct the actual yield. The specific formula is as follows:

[0204] Y actual =Y potential ·(1-δ drought -δ erosion ) (37);

[0205] In the formula, δ drought represents the drought stress loss rate, which is inverted by the vegetation temperature condition index; δ erosion It represents the yield reduction rate caused by soil erosion and is negatively correlated with soil conservation.

[0206] (5) Biomass energy supply service assessment;

[0207] The goal is to evaluate the potential for utilization of biomass energy such as terraced rice straw and oil crops. The technical solution is to use a multi-source biomass inversion model to first calculate the energy content of the relevant crop straw. The specific formula is as follows:

[0208]

[0209] Y in the formula i represents the yield of crop type i (kg / ha), and n represents the number of crop types; RPR irepresents the grass-to-grain ratio; η i represents the energy conversion efficiency. Then the biodiesel potential of oil crops is calculated. The specific formula is as follows:

[0210] E biodiesel =ρ·A·Y oil ·η transester (39);

[0211] In the formula, ρ represents the ratio of crop planting area to terrace area; Y oil represents the oil yield per unit area (L / ha), referring to the FAO crop parameter database; η transester It represents the conversion rate of transesterification reaction (typical value 0.95).

[0212] This technical solution, based primarily on spatial and temporal data from terraced fields, proposes an ecosystem service assessment framework based on "dynamic parameter modeling, multi-process coupling, and multi-dimensional collaborative optimization." This framework focuses on five core functions: quantifying soil conservation services, assessing water conservation services, assessing carbon sequestration services, assessing food supply services, and assessing biomass energy supply services. By leveraging high-resolution remote sensing data to dynamically modify model parameters, this approach transcends the limitations of traditional static parameters and provides key technical support for the ecological protection and sustainable management of terraced fields.

[0213] In S103, to reveal the interactions (trade-offs / synergies) between terraced ecosystem services, a comprehensive technical approach of "spatiotemporal heterogeneity analysis - driver factor decoupling - multi-objective optimization mapping" was proposed. By coupling geostatistical models, machine learning causal inference, and spatially explicit optimization algorithms, the nonlinear relationships between soil conservation service quantification, water conservation service assessment, carbon sequestration service assessment, food supply service assessment, and biomass energy supply service assessment were analyzed at multiple scales. The specific method is as follows:

[0214] (1) Framework for analyzing trade-off synergy relationships;

[0215] The goal is to quantify the trade-off intensity and synergy gain between different terraced ecosystem services and identify the primary driving factors. The technical approach involves three aspects: first, data input, based on the spatial distribution maps of the five ecosystem services generated in the previous chapter, a vector map of terrace distribution, multi-temporal environmental variables (rainfall, temperature, topography), and an index of human activity interference intensity; second, relationship analysis, using a multi-model cascade analysis (correlation statistics, causal inference, and spatial clustering); and third, output, generating a trade-off synergy intensity matrix, a spatial heterogeneity heat map, and a set of multi-objective optimization solutions.

[0216] (2) Constructing a framework for analyzing trade-off synergy relationships;

[0217] Step 1: Statistical correlation analysis (base layer);

[0218] Calculate the Pearson correlation coefficient and partial correlation coefficient between services based on grid cells:

[0219]

[0220] R in the formula jk represents the correlation coefficient between service j and service k. Positive values ​​represent synergistic relationships, while negative values ​​represent trade-off relationships. In addition, significance correction was performed using Bonferroni correction to eliminate false positives from multiple testing, retaining only strong correlations with p < 0.01.

[0221] Step 2: Machine learning causal identification (core layer);

[0222] A two-stage causal mechanism decoupling model (TCMD) is proposed:

[0223] Phase 1 (driving factor screening): Random forest feature importance ranking is used to screen the main controlling factors that affect the trade-off between services (such as slope, vegetation coverage, and terrace ridge density).

[0224] Phase 2 (Causal Effect Estimation): Based on the Double ML framework, the causal effects of natural factors and human activities are separated. The specific formula is as follows:

[0225] Trade off jk =α·X natural +β·X human +∈ (41);

[0226] X in the formula natural represents natural factors (climate, soil thickness, etc.), and instrumental variables (such as altitude) are added to control endogeneity; X human It represents human activities (fertilizer application amount, terrace maintenance frequency), which are inverted through night light data and farmer survey data.

[0227] Step 3: Calculation of trade-off synergy index (output layer);

[0228] Define the Spatiotemporal Trade-off Index (STTI):

[0229]

[0230] STTI in the formula jk (x, t) represents the trade-off index between services j and k at coordinate x and time t (negative value = trade-off, positive value = synergy); represents the marginal change rate between services based on Geographically Weighted Regression (GWR); D max represents the maximum possible trade-off distance in the study area (taking the 99% quantile of all D(x,t)); D(x,t) represents the trade-off distance, reflecting the Euclidean distance (degree of deviation) between the current state and the theoretical optimal state, which can be specifically expressed by the following formula:

[0231]

[0232] In the formula It represents the expected value of the mth service under the theoretical optimal state (determined by the production possibility frontier PPF).

[0233] Step 4: Spatially explicit collaborative optimization (application layer);

[0234] A multi-objective particle swarm optimization model (MOPSP) is constructed. The objective function is to maximize the total synergistic gain of the five ecosystem services. The specific formula is as follows:

[0235]

[0236] In the formula, Z represents the total synergy gain objective function; Synergy jk represents the synergistic gain of services j and k (if R jk ≥0, it is the actual value; otherwise it is 0); w jk The constraints are twofold: the ecosystem service non-decline threshold, such as soil conservation ≥ 80% of the current level; and the upper limit of terrace engineering transformation, such as slope ≤ 25°.

[0237] The above steps yield a synergistic trade-off map for the five ecosystem services provided by terraces, further enabling the formation of spatial hotspot maps (e.g., zones of conflict between high grain production and low water conservation) and priority areas for spatial optimization (regions with low STTI values). Furthermore, by summarizing typical ecosystem service trade-offs across terraces in different climate regions (e.g., the strong carbon sequestration-grain trade-off in dryland terraces in northern China), this can provide a basis for cross-regional policy adjustments.

[0238] This embodiment achieves accurate quantification of the multi-ecosystem service relationships of terraces through the deep integration of causal inference and spatial explicit optimization, providing a scientific and practical tool for multi-objective collaborative management, and filling the gaps in traditional methods in spatiotemporal dynamic analysis and refined management.

[0239] In summary, this embodiment, through system integration and innovative transformation, at the public technology level, adopts SIFT feature matching of multi-source remote sensing data and Poisson reconstruction of LiDAR point clouds to achieve geometric alignment, optimizes the open source fragmentation index (FI) and stability index (SSI) extraction modules for landscape ecology, and constructs a basic evaluation framework based on classic ecological models (RUSLE soil erosion equation, SCS-CN runoff model, and light energy utilization rate (NPP) model) to ensure the scientific nature and verifiability of the method system. In terms of innovation, a "four-dimensional dynamic weight allocation algorithm" (spatial resolution-coverage efficiency-timeliness-redundancy dynamic balance model) was proposed for the first time to solve the nonlinear contradiction in the coordinated scheduling of multispectral / high-resolution data, and the scheduling efficiency was improved by 40% compared with traditional methods; the MMF-SegNet multimodal segmentation network driven by the terrain-guided attention mechanism (TGA) was developed, and the terrace boundary recognition accuracy (ES) was improved from 0.74 to 0.91 by fusing the optical texture temporal characteristics and LiDAR terrain curvature information; breaking through the limitations of traditional static models, the dynamic assignment technology of red edge band remote sensing inversion parameters (dynamic modeling of vegetation management factor C random forest) was introduced to construct a spatial heterogeneity adaptive assessment framework for soil conservation and water conservation services, and the verification error was reduced by 15%-22%; the independent design of the spatiotemporal explicit trade-off collaborative analysis technology (STTI) integrated dual machine learning causal inference and the MOPSO multi-objective optimization algorithm to analyze the natural-human driving relationship in terrace multi-ecosystem services and support the generation of differentiated regulatory rules. Through the full-chain logical closed loop of "dynamic data coupling-intelligent feature mining-heterogeneity assessment-cooperative optimization", the solution not only overcomes the shortcomings of traditional methods in mesoscale monitoring, such as model rigidity and service association fragmentation, but also significantly reduces R&D risks and deployment thresholds through compatibility adaptation of mature technologies, providing a high-precision and scalable toolset for the sustainable management of terraced ecosystems.

[0240] The technical effects are as follows:

[0241] (1) Dynamic optimization acquisition technology solves the contradiction between mesoscale monitoring efficiency and coverage; a spatial resolution-coverage efficiency-timeliness-redundancy balance model is constructed through a "four-dimensional dynamic weight allocation algorithm", and an event-driven adaptive grid division and feedback mechanism is combined to achieve intelligent scheduling of multi-source data. Data redundancy is dynamically optimized based on the Herfindahl index, the priority weight of drone supplementation is adjusted according to the complexity of the terrain, and invalid acquisition is eliminated through adaptive grids. Compared with the traditional fixed weight method, the scheduling efficiency is improved by 40%, and the full area coverage cycle is shortened by 50% within the 1:5000 to 1:50000 mesoscale unit. At the same time, the amount of redundant data is reduced by 30%, effectively solving the contradiction between "high-resolution local coverage" and "lack of timeliness of wide-area observation" in the background technology.

[0242] (2) Heterogeneous data preprocessing enhances the recognizability of terrace features; using spatiotemporal benchmark unification, NDCI dynamic cloud removal, Poisson surface reconstruction, and residual slope correction techniques to eliminate geometric and spectral errors in multi-source data. Through radiometric normalization (Table 1) and terrain feature enhancement (SLIC superpixel segmentation combined with Canny edge detection), the boundary blurring problem caused by shadow distortion and cloud occlusion in traditional preprocessing is overcome. The terrain parameter inversion error is reduced to ±0.5°, and the effective pixel retention rate of NDCI cloud removal is increased to 92%, providing high-precision input for subsequent terrace segmentation and solving the core defect of "poor mesoscale applicability" in the background technology.

[0243] (3) A multimodal segmentation network improves the accuracy of terrace boundary recognition. A terrain-guided attention mechanism (TGA) is designed to drive the MMF-SegNet network, which integrates optical texture, terrain curvature, and temporal NDRE features to construct a dynamic weighted multimodal segmentation model. The TGA mechanism dynamically suppresses shadow interference through the terrain roughness index (RI). In steep slope areas (RI>15), the optical weight is reduced to 0.3, significantly enhancing the contour features of the ridges. The terrace boundary recognition accuracy (ES) is improved from 0.74 of the traditional single-modal model to 0.91, and the missed detection rate of fragmented terraces (FI>0.7) is reduced by 45%, solving the pain point of "single-phase static analysis leads to insufficient representation of dynamic evolution" in the background technology.

[0244] (4) Dynamic parameter modeling is used to evaluate the heterogeneity of ecosystem services. The red edge band is introduced into classic models such as RUSLE and SCS-CN to invert the vegetation management factor C and the dynamic assignment of the CN value of the ridge water storage capacity classification to build a spatial adaptive evaluation framework. The model parameters are dynamically corrected based on the time series NDRE and LiDAR inversion of the tillage layer thickness, breaking through the data lag defect of the traditional static equivalent factor method. The soil conservation verification error is reduced by 22%, and the water conservation simulation accuracy (R 2 ) was improved from 0.68 to 0.83, resolving the key limitation of “static ecosystem service assessment” in the background technology.

[0245] (5) Causal optimization coupling technology enhances the effectiveness of collaborative service decision-making; a spatiotemporal explicit trade-off index (STTI) and a dual machine learning causal inference-MOPSO optimization cascade model are proposed to analyze nonlinear relationships between services and generate spatial control rules. STTI combines marginal change rate and trade-off distance to quantify spatiotemporal heterogeneity, and the MOPSO optimization algorithm generates a multi-objective Pareto solution set of slope, vegetation cover, and terrace ridge density through particle swarm optimization. In typical conflict areas (such as "high grain production and low water conservation" areas), the collaborative gain is increased by 35%, and the efficiency of policy adaptation plan generation is increased by 50%, breaking through the dual bottlenecks of "statistical correlation cannot analyze the interaction between natural and human factors" and "poor spatial implementation of optimization results" in the background technology.

[0246] In summary, the system systematically solves the problems of data fragmentation, model rigidity, and extensive service analysis in mesoscale terrace monitoring. Its core indicators (identification accuracy, assessment error, and synergistic gain) are improved by 35%-50% compared with existing technologies, providing a technical tool that is both scientific and practical for the sustainable management of terrace ecosystems.

[0247] The various embodiments in this specification are described in a progressive manner, with each embodiment focusing on the differences from other embodiments. Reference can be made to the common and similar parts between the various embodiments. For the devices disclosed in the embodiments, since they correspond to the methods disclosed in the embodiments, the description is relatively simple, and the relevant parts can be referred to the method description.

[0248] The above description of the disclosed embodiments is intended to enable one skilled in the art to implement or use the present invention. Various modifications to these embodiments will be readily apparent to one skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the present invention. Therefore, the present invention is not limited to the embodiments shown herein but is intended to conform to the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A method for quantifying terrace pattern evolution and its ecosystem service trade-off, characterized by: include: A multi-source remote sensing information collaborative framework was constructed to identify frequently changing areas. Multi-source remote sensing data was pre-processed to obtain a standard input dataset through temporal and spatial benchmark unification, key feature enhancement, redundant information compression, and quality dynamic feedback optimization. The specific steps for constructing the multi-source remote sensing information collaborative framework are as follows: A dynamic weight allocation algorithm is used to achieve a global optimal balance between spatial resolution, coverage efficiency, and time cost through mathematical modeling. The Herfindahl index is also introduced to measure data diversity and redundancy to ensure that weight allocation is comprehensive while avoiding information overload. When necessary, adaptive grid division is performed on the study area to ensure smooth image data integration. Construct a "macro-representation-process capture-micro-calibration" framework to carry out complementary fusion of multi-source data; Utilize a mesoscale adaptive data collection command feedback mechanism to improve data collection efficiency and identify frequently changing areas; Based on a standard input dataset, a multimodal boundary refinement segmentation algorithm is used to extract quantitative data on the distribution of terraces in typical regions and their interannual expansion / contraction dynamics, representing their spatiotemporal evolution. By integrating multimodal deep learning with spatial statistical models, a comprehensive quantitative analysis of the terrace evolution process is achieved. A terrain-guided attention mechanism (TGA) is designed, driven by the MMF-SegNet network, which integrates optical texture, terrain curvature, and temporal NDRE features to construct a dynamic weighted multimodal segmentation model. Based on the extracted terrace distribution data, combined with the standard input dataset and ground observation data, an evaluation framework is proposed to quantitatively characterize the terraces in soil conservation services, water conservation services, carbon sequestration services, food supply services and biomass energy supply services, and conduct a synergistic analysis of the trade-offs between typical ecosystem services; a spatiotemporal explicit trade-off index (STTI) and a dual machine learning causal inference-MOPSO optimization cascade model are proposed to analyze the nonlinear relationship between services and generate spatial regulation rules.

2. The method for quantifying terrace pattern evolution and its ecosystem service trade-off according to claim 1 is characterized in that: The weight calculation model of the dynamic weight allocation algorithm is as follows: Among them, W i Represents the priority weight of data source i, which determines its collection order in scheduling and is dimensionless; SR i Indicates the spatial resolution of data source i, in meters; CR i Indicates the coverage efficiency of data source i, that is, the coverage area per unit time, in km 2 / day; λ0 represents the time-dependent attenuation coefficient; Δt represents the time delay after data collection, in days; HI(D similar ) is the Herfindahl index, which measures the redundancy of similar data sources, where s j is the ratio of the coverage area of ​​data source j to the total overlapping area; k represents the total number of data sources.

3. The method for quantifying terrace pattern evolution and its ecosystem service trade-off according to claim 1 is characterized in that: The specific steps for extracting terrace distribution data for typical areas are as follows: Define optical image features, terrain parameter features, and vegetation time series features as input features; Constructing a dynamic weighted multimodal segmentation network based on the input features to obtain a segmentation probability map of the terraced area; The segmentation probability map is subjected to threshold segmentation, and concave points are detected on the terrace patch contour, the concave points are connected, and the obtained binary mask is converted into a polygonal vector layer, and the final terrace distribution vector map is output.

4. The method for quantifying terrace pattern evolution and its ecosystem service trade-off according to claim 3 is characterized in that: The specific steps to build a dynamic weight multimodal segmentation network are as follows: To build MMF-SegNet, we first extract branch features. In the optical branch, we use a dilated residual module to capture multi-scale textures and output an optical feature map. Secondly, in the terrain branch, we input two-channel data consisting of slope and terrain curvature, encode it using a lightweight U-Net, and output terrain features. Finally, for time series analysis, we input NDRE time series and use a 1D Conv-LSTM to extract time-dependent features, outputting a vector. Finally, for feature module fusion, we introduce terrain-guided attention and dynamically fuse optical and terrain features to obtain fused features. F fusion =α·F RGB +(1-α)·F DEM ; α=σ(RI·W); Where RI represents the terrain roughness index, which is a composite of slope and terrain curvature; W refers to the learnable parameter, which is optimized by back propagation; F RGB Represents the optical image characteristics, F DEM represents the terrain parameter characteristics; σ is the Sigmoid function; the weight coefficient α is dynamically adjusted by the terrain complexity; θ is the slope angle calculated based on the DEM, reflecting the degree of surface inclination; C represents the terrain curvature; Finally, the vegetation temporal feature F NDRE and fusion feature F fusion Splicing, through convolution fusion, outputs the segmentation probability map of the terraced area.

5. The method for quantifying terrace pattern evolution and its ecosystem service trade-off according to claim 1 is characterized in that: The specific steps for quantifying the spatiotemporal evolution of terraces and extracting the dynamic representation of interannual expansion / contraction are as follows: Calculate the rate of change in terrace area, count the number of expanded / reduced terrace patches, their spatial distribution, and area proportion; Evolution of spatial pattern index, the selected spatial pattern indexes include fragmentation index, morphological stability index, center of gravity migration trajectory, and interannual migration distance; To generate the spatiotemporal evolution map of terraces, we first generate a density map of the expansion / contraction of terraces year by year in the form of a heat map sequence; then we generate a change trajectory map, and construct a change path network diagram G = (V, E), where the node V is the centroid of the change patch, the edge E connects the centroids of adjacent years, and the weight is the migration distance; finally, we construct a three-dimensional space-time cube model, stacking the two-dimensional masks between years according to the pixels, and display the evolution trend of the vertical time axis through voxel rendering to generate a visual map.

6. The method for quantifying terrace pattern evolution and its ecosystem service trade-off according to claim 5 is characterized in that: To calculate the rate of change in terrace area, the number of expanded / reduced terrace patches, their spatial distribution, and area proportion were counted. The specific steps are as follows: The specific formula for calculating the rate of change of terrace area is as follows: ΔA t = Yes t+1 -IN t ; In the formula, A t represents the total area of ​​terraces in the study area in year t, S k represents the area of ​​the kth terrace patch, K represents the total number of terrace patches; A t+1 represents the total area of ​​terraces in the study area in year t+1, in square meters or hectares; ΔA t Represents the interannual variation of terrace area; if ΔA t > 0 and has been growing for three consecutive years, marked as a stable expansion zone. t <-0.2A t , determined to be a damage event; Perform spatial overlay analysis and detect the changed area pixel by pixel. The formula is as follows: In the formula, C t (i, j) represents the category label of pixel (i, j) from time t to t+1; M t (i, j) represents the binary mask of year t, 0 represents non-terraced fields and 1 represents terraced fields.

7. The method for quantifying terrace pattern evolution and its ecosystem service trade-off according to claim 1 is characterized in that: The specific steps of the assessment framework to quantitatively characterize the functions of terraced ecosystems are as follows: Quantification of soil conservation services: The improved RUSLE model is used, and the specific formula is as follows: A actual =R·K·L·S·C terrace ·P terrace ; A potential =R·K·L·S; SOil Conservation=A potential -A actual ; In the formula, A potential and A actual They represent the potential soil erosion and actual soil erosion in the region, respectively. The difference between the two represents the amount of soil retained by terrace projects in the region. R represents rainfall erosivity, which is calculated based on the total annual rainstorm erosivity based on local rain gauges or CHIRPS remote sensing rainfall data. K represents soil erodibility, which is parameterized by the EPIC model using soil texture maps. L·S represents the topographic factor, which generates slope and slope length based on satellite elevation data of the terrace distribution area. The formula is: S=10.8sinθ+0.03, θ<5°; S=16.8sinθ-0.5,θ≥5°; In the formula, λ is the slope length, unit: meter. After extracting the watershed watershed, the actual slope length is calculated by dividing the fields into terraces. m is the slope length factor index; θ is the slope angle of the slope; C terrace Represents the terrace vegetation cover management factor, dimensionless, based on the dynamic inversion of the red edge band: NDRE max represents the maximum NDRE during the crop growing season; P terrace represents the terrace engineering factor, dimensionless; Water conservation service assessment: A distributed water conservation model for terraces was constructed, coupling the spatial heterogeneity of terrace structure and soil hydraulic parameters. The specific formula is as follows: AND i =ET natural ·(1-γ·NDWI wet ); In the formula, P i represents the annual precipitation of grid cell i, obtained through interpolation of weather stations or GSMaP satellite precipitation products; ET i represents the actual evapotranspiration of the terraced field in grid cell i; R i A represents the surface runoff depth of the terrace in grid cell i; i represents the area of ​​grid cell i; η i represents the soil infiltration enhancement coefficient, dimensionless, and dynamically corrected based on soil type and terraced tillage layer thickness; ET natural represents the evapotranspiration of natural vegetation; where γ represents the interception coefficient of terraced fields in the wet season, which is calibrated by ground leakage experiments, and NDWI wet is the canopy moisture index of crops in the wet season; it is calculated using the SCS-CN model: CN value represents the runoff curve number, dimensionless; P represents annual precipitation; S max Indicates the maximum potential retention of water; Carbon sequestration service assessment: Construct a hierarchical carbon sink model. First, the vegetation carbon sequestration amount. The specific formula is as follows: NPP terrace =APAR·e·f crop ; In the formula, APAR represents absorbed photosynthetically active radiation; ε represents light energy utilization efficiency; f crop is the crop rotation coefficient; the second is the soil carbon sequestration amount, and the specific calculation formula is as follows: In the formula, C input Represents the carbon input of crop litter; SOC initial represents the initial soil organic carbon content; α and β represent the organic carbon transformation and mineralization rates, which are assigned according to the tillage method; BD represents the soil bulk density; D represents the tillage layer depth; Food supply service assessment: Using the mixed yield model, first calculate the potential yield. The specific formula is as follows: Y potential =NPP·HI·CF; In the formula, NPP stands for net primary productivity; HI stands for harvest index, which is dimensionless; and CF stands for carbon content-yield conversion coefficient. After adjusting for actual yield, the specific formula is as follows: Y actual =Y potential ·(1-d drought -d erosion ); In the formula, δ drought represents the drought stress loss rate, which is inverted by the vegetation temperature condition index; δ erosion It represents the yield reduction rate caused by soil erosion and is negatively correlated with soil conservation; Biomass energy supply service assessment: Utilize the multi-source biomass inversion model to calculate the energy content of relevant crop straw. The specific formula is as follows: Y in the formula i represents the yield of crop type i, and n represents the number of crop types; RPR i represents the grass-to-grain ratio; η i represents the energy conversion efficiency; then the biodiesel potential of oil crops is calculated, the specific formula is as follows: E biodiesel =ρ·A·Y oil ·or transester ; In the formula, ρ represents the ratio of crop planting area to terrace area, dimensionless; A is the total area of ​​terrace; Y oil represents the oil production per unit area; η transester It represents the conversion rate of transesterification reaction.

8. The method for quantifying terrace pattern evolution and its ecosystem service trade-off according to claim 1 is characterized in that: The specific steps for conducting a collaborative analysis of trade-offs between ecosystems are as follows: Construct a framework for analyzing trade-off relationships: First, data input, coordinating generated data; second, relationship analysis, using multiple models to couple relevant data information; third, result output, generating a trade-off relationship matrix, spatial heterogeneity heat map, and a set of multi-objective optimization solutions; The multi-scale trade-off synergy analysis includes the following four steps: The foundation layer was constructed and statistical correlation analysis was performed to conduct multi-scale terrace ecosystem service correlation analysis; Core layer construction, machine learning causal identification, and a two-stage causal-mechanism decoupling model are proposed: Random forest feature importance ranking is used to screen the main controlling factors that affect the trade-offs between services; based on a dual machine learning framework, the causal effects of natural factors and human activities are separated; Output layer construction, trade-off synergy index calculation, and definition of spatiotemporal dynamic trade-off index; The application layer is constructed, spatial explicit collaborative optimization is performed, and a multi-objective particle swarm optimization model is built. The objective function is to maximize the total synergistic gain of five ecosystem services. The constraints include the non-decreasing threshold of ecosystem services and the upper limit of terrace engineering transformation.

Images