Terraced field pattern evolution quantification and ecological system service tradeoff method

Abstract

The invention discloses a terrace pattern evolution quantification and ecological system service tradeoff method, and relates to the technical field of crossing of an agricultural ecological technology and a geographic information technology, and the method comprises the steps: constructing a multi-source remote sensing information cooperation framework, obtaining multi-source remote sensing data, carrying out the preprocessing of the multi-source remote sensing data, and obtaining a standard input data set; on the basis of a standard input data set, spatio-temporal evolution quantitative data such as terrace distribution of a typical area and interannual expansion / contraction dynamic representation of the terrace are extracted by using a multi-modal boundary refined segmentation algorithm, and comprehensive quantitative analysis of the terrace evolution process is realized through linkage of multi-modal deep learning and a spatial statistical model; and based on the extracted terrace distribution data, combining a standard input data set and ground observation data, proposing an evaluation framework, and carrying out trade-off collaborative analysis among terrace ecosystem services. According to the method, the problems that in the prior art, the temporal-spatial resolution is insufficient, and the provided ecological system service evaluation is static and the decision support is weak are solved.

Description

Technical Field

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

[0002] At present, the current research results on terrace pattern monitoring and quantitative analysis technology are mainly divided into two categories. One is the morphological extraction technology based on remote sensing image interpretation. This type of technology is the most widely used. Its main feature is to use high-resolution satellite images (such as Sentinel-2, GF-2), combined with object-oriented classification or some other automatic classification algorithms to extract terrace boundaries, supplemented by morphological algorithms (such as edge detection, texture analysis) to quantify geometric parameters (area, slope, field fragmentation, etc.). However, the defects of this type of technology mainly include two aspects. First, the spatiotemporal dynamics are insufficient. The method focuses on single-phase static analysis, lacks the evolution of long-term time series and its driving force analysis, and is difficult to support comprehensive optimization of prediction; second, the technology is poorly applicable in the mesoscale. This technology is mostly designed for small-scale high-precision (<10 km2) or large-scale low-precision (>100km2), and lacks a dedicated modeling method for 1:5000 to 1:50000 mesoscale terrain units. The second is the scenario simulation technology based on geographic information system. The main feature of this technology is to use land change models (such as CLUE-S, CA-Markov) to simulate the expansion / contraction trend of terraces, and combine terrain characteristic indexes (such as terrain wetness index TWI, slope variability index SVI) to evaluate the evolution characteristics and laws of spatial pattern. The defects of this type of technical method are reflected in two aspects. First, the driving factors are simplistic. The model relies on natural factors (topography and precipitation, etc.) and ignores the dynamic interaction of socio-economic factors (such as labor transfer and policy compensation); in addition, the measurement module of service trade-off is missing, and the simulation results are mainly based on area changes. The dynamic response function of ecosystem service value is not coupled, and it is impossible to reveal the mechanism of pattern evolution on ecosystem services and their interrelationships.

[0003] The methods of ecosystem service trade-off analysis can be divided into two categories. The first is the static ecosystem service value assessment method. This type of method often uses unit area value equivalent (such as InVEST model) or Emergy Analysis (EMA) (Nadalini et al., 2021) to calculate the value of ecosystem services, and analyzes the priority of ecosystem services through spatial superposition. This type of assessment method is relatively common and widely used, but it has the defect of being out of touch with time. It relies on base year data and does not introduce a dynamic correction mechanism for time series, which causes the results to lag behind the actual evolution process. In addition, the synergistic trade-off relationship between ecosystem services is simplified, and the service trade-off relationship is often described by linear regression or correlation analysis, which makes it difficult to analyze nonlinear threshold effects (such as when the fragmentation of terraces exceeds the critical value, its carbon sink function suddenly decreases). The second is a decision-making technology based on multi-objective optimization, which is mostly based on Pareto Frontier (PF) or Genetic Algorithm (GA) to seek the optimal solution set and balance the service conflict objective function. 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, which makes it difficult for the planning scheme to adapt to changes in long-term policies or environmental conditions.

[0004] In general, focusing on the synergistic analysis of the trade-offs between terrace pattern evolution and ecosystem services, previous ideas were mostly to identify and extract target land types by combining remote sensing images with land use data, 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, so as to overcome the insufficient spatiotemporal resolution in existing technologies, the static nature of ecosystem service assessments, and weak decision-making support, is an urgent problem that technicians in this field 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 solution: In one aspect, the present invention discloses a method for quantifying the evolution of terrace pattern and its ecosystem service trade-off, comprising: Construct a multi-source remote sensing information coordination framework, identify frequently changing areas, pre-process multi-source remote sensing data, and obtain a standard input data set through the unification of spatiotemporal benchmarks, key feature enhancement, redundant information compression, and quality dynamic feedback optimization; Based on the standard input data set, a multimodal boundary refinement segmentation algorithm is used to extract the distribution of terraces in typical areas and their interannual expansion / contraction dynamics to characterize the quantitative temporal and spatial evolution. Through the linkage of multimodal deep learning and spatial statistical models, a comprehensive quantitative analysis of the terrace evolution process is achieved. Based on the extracted terrace distribution data, combined with the standard input data set and ground observation data, an evaluation framework is proposed to quantitatively characterize the terraces in soil conservation service quantification, water conservation service assessment, carbon sequestration service assessment, food supply service assessment and biomass energy supply service assessment, and to conduct a synergistic analysis of the trade-offs among the above typical ecosystem services.

[0008] Preferably, in the above-mentioned method for quantifying the evolution of terrace pattern and its ecosystem service trade-off, the specific steps of constructing a multi-source remote sensing information coordination framework are as follows: Using a dynamic weight allocation algorithm, mathematical modeling is used to achieve a global optimal balance between spatial resolution, coverage efficiency, and time cost. The Herfindahl index is introduced to measure data diversity redundancy to ensure that the weight allocation is comprehensive and avoids information overload. When necessary, the study area is adaptively gridded to ensure good image data connection. Construct a "macro-representation-process capture-micro-calibration" framework to carry out complementary fusion of multi-source data; The mesoscale adaptive data collection command feedback mechanism is used to improve data collection efficiency and identify frequently changing areas.

[0009] Preferably, in the above-mentioned method for quantifying the evolution of terrace pattern and its ecosystem service trade-off, the weight calculation model of the dynamic weight allocation algorithm is as follows: ; ; in, It represents the priority weight of data source i, which determines its collection order in scheduling and is dimensionless; Indicates the spatial resolution of data source i, in meters; Represents the data source The coverage efficiency is the coverage area per unit time, in km2 / day; represents the time-dependent attenuation coefficient; Indicates the time delay after data collection, in days; is the Herfindahl index, which measures the redundancy of similar data sources. is the ratio of the coverage area of ​​data source j to the total overlapping area; k represents the total number of data sources.

[0010] Preferably, in the above-mentioned method for quantifying the evolution of terrace pattern and its ecosystem service trade-off, the specific steps of extracting terrace distribution data for a typical area are as follows: Define optical image features, terrain parameter features, and vegetation time series features as input features; Building a dynamic weight multimodal segmentation network based on the input features to obtain a segmentation probability map of the terraced field area; The segmentation probability map is subjected to threshold segmentation, and concave points are detected on the contour of the terraced field patches. The concave points are connected, and the obtained binary mask is converted into a polygonal vector layer, and the final terraced field distribution vector layer is output.

[0011] Preferably, in the above-mentioned method for quantifying the evolution of terrace pattern and its ecosystem service trade-off, the specific steps of constructing a dynamic weight multimodal segmentation network are as follows: To build MMF-SegNet, we first extract branch features. In the optical branch, we use the dilated residual module to capture multi-scale textures and output an optical feature map. Secondly, in the terrain branch, we input the 2-channel data of slope and terrain curvature splicing, use lightweight U-Net encoding, and output terrain features. Finally, for time series analysis, we input NDRE time series, use 1D Conv-LSTM to extract time-dependent features, and output vectors. In the fusion of feature modules, we introduce terrain-guided attention and dynamically fuse optical and terrain features to obtain fused features: ; ; ; In the formula, Represents the terrain roughness index, which is a composite of slope and terrain curvature; Refers to learnable parameters, optimized by back-propagation; Represents the optical image characteristics, Represents terrain parameter characteristics; is the Sigmoid function; weight coefficient Dynamically adjusted by terrain complexity; It is the slope angle calculated based on DEM, reflecting the degree of surface inclination; C represents the terrain curvature; Finally, the vegetation temporal characteristics With fusion features Splicing, through convolution fusion, output the segmentation probability map of the terrace area.

[0012] Preferably, in the above-mentioned method for quantifying the evolution of terrace pattern and its ecosystem service trade-off, the specific steps of quantifying the spatial and temporal evolution characteristics of terraces and extracting the dynamic representation of interannual expansion / contraction are as follows: Calculate the rate of change of terrace area, count the number of expanded / reduced terrace patches, their spatial distribution and area proportion; The evolution of spatial pattern index, the selected spatial pattern index, specifically includes fragmentation index, morphological stability index, centroid 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; secondly, we generate a change trajectory map and construct a change path network map (G=(V, E)), where node V is the centroid of the change patch, 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, stack 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.

[0013] Preferably, in the above-mentioned method for quantifying the evolution of terrace pattern and its ecosystem service trade-off, the calculation of terrace area change rate, the statistical number of expanded / reduced terrace patches and their spatial distribution and area proportion are performed, and the specific steps are as follows: The specific formula for calculating the change rate of terrace area is as follows: ; ; in the formula represents the total area of ​​terraces in the study area in year t, represents the area of ​​the kth terrace patch, and K represents the total number of terrace patches; represents the total area of ​​terraces in the study area in year t+1, in square meters or hectares; represents the annual change in the area of ​​terraced fields; if And it has been growing for three consecutive years, marked as a stable expansion area. , determined as a damage event; Perform spatial overlay analysis and detect the change area pixel by pixel. The formula is as follows: ; in the formula Represents the category label of pixel (i, j) from time t to t+1; Represents the binary mask of year t, 0 represents non-terraced fields and 1 represents terraced fields.

[0014] Preferably, in the above-mentioned method for quantifying the evolution of terrace pattern and its ecosystem service trade-off, the specific steps of the evaluation framework for quantitatively describing the function of terrace ecosystem are as follows: Quantification of soil conservation services: The improved RUSLE model is used, and the specific formula is as follows: ; ; ; in the formula and They represent the potential soil erosion and actual soil erosion in the region, respectively, and the difference between the two represents the amount of soil retained in the region by terrace engineering; represents rainfall erosivity MJ·mm / (hm²·h·yr), which is based on local rain gauges or CHIRPS remote sensing rainfall data to calculate the total annual rainstorm erosivity; K represents soil erodibility t·hm²·h / (hm²·MJ·mm), which is parameterized by the EPIC model using soil texture maps; Represents the terrain factor, and generates the slope and slope length based on the satellite elevation data of the terrace distribution area. The formula is: ; ; ; in the formula is the slope length, unit: meter. After extracting the watershed watershed line, the actual slope length is calculated by dividing the terraced fields. m is the slope length factor index. is the slope angle of the slope; Represents the terrace vegetation cover management factor, dimensionless, based on the dynamic inversion of the red edge band: ; It indicates the maximum NDRE value during the crop growing season; represents the terrace engineering factor, dimensionless; Water conservation service evaluation: A distributed water conservation model for terraces is constructed, coupling the spatial heterogeneity of terrace structure and soil hydraulic parameters. The specific formula is as follows: ; ; ; ; in the formula represents the annual precipitation of grid cell i, in mm, obtained through interpolation from weather stations or GSMaP satellite precipitation products; represents the actual evapotranspiration of the terrace in grid unit i, mm, corrected by the water-saving effect of the terrace; Represents the surface runoff depth of the terrace in grid unit i, unit: mm; represents the area of ​​grid cell i; represents the soil infiltration enhancement coefficient, dimensionless, dynamically corrected based on soil type and terraced field cultivation layer thickness; ETnatural represents natural vegetation evapotranspiration, mm; where It represents the interception coefficient of terraced fields in the wet season, which is calibrated by ground leakage experiment. It is the canopy moisture index of crops in the wet season; it is calculated using the SCS-CN model: the CN value represents the runoff curve number, which is dimensionless; represents annual precipitation, mm; Indicates the potential maximum retention of water, mm; Carbon sequestration service assessment: Construct a hierarchical carbon sink model, first of all, vegetation sequestration, the specific formula is as follows: ; in the formula Indicates absorbed photosynthetically active radiation, MJ / m²; It represents the light energy utilization rate, gC / MJ; is the crop rotation coefficient; the second is soil carbon sequestration, the specific calculation formula is as follows: ; in the formula represents the carbon input of crop litter, kg / ha / yr; represents the initial soil organic carbon content, %; , It represents the rate of organic carbon transformation and mineralization, and is assigned according to the tillage method; BD represents soil bulk density, g / cm³; D represents the depth of the tillage layer, m; Food supply service assessment: Using the mixed yield model, first calculate the potential yield. The specific formula is as follows: ; in the formula represents net primary productivity, gC / m2 / year; represents the harvest index, dimensionless; Indicates the carbon content-yield conversion coefficient, kg / gC; combined with the actual yield correction, the specific formula is as follows: ; in the formula represents the drought stress loss rate, which is inverted by the vegetation temperature condition index; It indicates the yield reduction rate caused by soil erosion, which is negatively correlated with soil conservation; Biomass energy supply service assessment: Use the multi-source biomass inversion model to calculate the energy content of relevant crop straw. The specific formula is as follows: ; In the formula, Yi represents the yield of the i-th crop, kg / ha, and n represents the number of crop types; represents the grass-to-grain ratio; It represents the energy conversion efficiency; then the biodiesel potential of oil crops is calculated, and the specific formula is as follows: ; in the formula It represents the proportion of crop planting area to terrace area, dimensionless; A is the total area of ​​terrace, ha; It indicates the oil yield per unit area, L / ha; It represents the conversion rate of transesterification reaction.

[0015] Preferably, in the above-mentioned terrace pattern evolution quantification and ecosystem service trade-off method, the specific steps of conducting the trade-off collaborative analysis among ecosystems are as follows: Construct a framework for analyzing trade-off synergy relationships: construct a basic layer and conduct multi-scale trade-off synergy relationship analysis through statistical 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 control 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 constructed. 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 the transformation of terrace projects.

[0016] Through the above technical solutions, it can be seen that compared with the prior art, the present invention provides a method for quantifying the evolution of terrace patterns and weighing their ecosystem services, overcoming the problems of insufficient spatiotemporal resolution, static ecosystem service assessment and weak decision-making support in the prior art; 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 the traditional method; 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 the optical texture temporal features with the LiDAR terrain curvature information; breaking through the limitations of the traditional static model, introducing the red edge band remote sensing inversion parameters 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 relies on the compatibility adaptation of mature technologies to significantly reduce R&D risks and deployment thresholds, providing a high-precision and scalable toolset for the sustainable management of terrace ecosystems. BRIEF DESCRIPTION OF THE DRAWINGS

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

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

[0019] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0020] like Figure 1As shown, the embodiment of the present invention discloses a method for quantifying the evolution of terrace pattern and its ecosystem service trade-off, including: S101 builds a multi-source remote sensing information coordination framework, identifies frequently changing areas, and pre-processes multi-source remote sensing data to obtain a standard input data set through the unification of spatiotemporal benchmarks, key feature enhancement, redundant information compression, and quality dynamic feedback optimization; S102 uses a multimodal boundary refinement segmentation algorithm based on a standard input data set to extract the distribution of terraces in typical regions and their interannual expansion / contraction dynamics to characterize the spatial and temporal evolution of quantitative data. Through the linkage of multimodal deep learning and spatial statistical models, a comprehensive quantitative analysis of the terrace evolution process is achieved. 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.

[0021] Specifically, in S101, in the data collection stage, 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 coordination framework, and the specific scheme is as follows: S1011 Dynamic data weight allocation algorithm based on mesoscale sensitivity: In mesoscale terrace monitoring, the existing data collection has a contradiction between "high-resolution local coverage" and "low-resolution global observation". For example, the coverage of high-resolution drone data is limited (1km² / flight), while the satellite revisit cycle of the entire area is long (such as 7 days), resulting in low monitoring efficiency. In response to this problem, this embodiment proposes a dynamic weight allocation algorithm, the goal of which is to achieve a global optimal balance between spatial resolution, coverage efficiency and time cost through mathematical modeling, and at the same time introduce the Herfindahl Index (HI) to measure data diversity redundancy to ensure that the weight allocation is comprehensive and avoids information overload. The weight calculation model is as follows: (1); (2); in, Represents the priority weight of data source i, which determines its collection order in scheduling. It is dimensionless and ranges from [0,1]; Indicates the spatial resolution of data source i. The smaller the value, the higher the accuracy. The unit is meter. Represents the data source The coverage efficiency (i.e., coverage area per unit time) is expressed in km2 / day; It represents the timeliness attenuation coefficient, which is used to control the negative impact of the newness and oldness 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.15, if the surface roughness is between 0.05 and 0.3 (hilly area), =0.2, if higher than 0.3 (mountainous area), then =0.25; Indicates the time delay after data collection, in days; is the Herfindahl index, which measures the redundancy of similar data sources and ranges from [0,1]. For data source The ratio of coverage area to total overlapping area.

[0022] In addition, the image data provided by the ground-based UAV equipment should also be prepared to connect with the satellite image data in the space. Therefore, the study area is adaptively gridded, and the grid size is divided according to the following formula: (3); Where G represents the adaptive grid size, which controls the minimum unit of data fusion and is in km; Indicates the minimum data source coverage efficiency requirement, in km2.

[0023] S1012 Complementary fusion of multi-source data: Different types of sensors have blind spots in feature expression in terrace identification, such as optical satellites are easily interfered by clouds and fog, drones lack terrain penetration, and laser radar (LiDAR) is expensive. The embodiment of the present invention uses a multi-level fusion strategy to build a data chain of "macro-representation-process capture-micro-calibration" to solve the perception limitations of a single data source. See Table 1 for details: Table 1 Parameters of different data sources, collection indicators and data fusion algorithms The slope error correction of lidar data and satellite elevation data adopts the residual least square method. The specific formula is as follows: (4); In the formula is the slope, is the standard deviation.

[0024] The selection logic of the above technologies is: Space-based WorldView-3 red edge band: Utilizes the sensitivity of the band near 750nm to vegetation biochemical parameters to enhance the spectral distinction between terraced field ridges and crops, and overcomes the fuzzy field ridge problem of traditional RGB images.

[0025] Airborne multispectral UAV time series interpolation: In order to address the gaps in cloud and fog caused by the satellite revisit cycle (for example, there are only 5 to 7 days of valid data per month in the rainy season), the vegetation growth curve is fitted based on Fourier transform to ensure monitoring continuity.

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

[0027] S1013 mesoscale adaptive data acquisition command feedback mechanism: The dynamic changes of terraces (such as heavy rains destroying ridges and farming leading to the merging of fields) are spatiotemporally heterogeneous, and traditional regular inspections result in a large amount of invalid data collection (about 60% of the inspections have no significant changes). By building an event-driven feedback mechanism, a closed-loop control of "data collection-change identification-task triggering" is achieved, and limited resources are focused on high-change areas. The fragmentation index is quoted here to reflect the degree of fragmentation of terraced fields. The specific formula is as follows: (5); In the formula It is expressed as the area of ​​a single terraced field in square meters; is the total area of ​​the region, in square meters; when FI is higher than 0.7, it can be determined as a highly broken terrace area.

[0028] The priority scheduling of the UAV supplementary mission trigger mechanism adopts the Hungarian algorithm to match the high-priority areas with idle UAVs. The specific formula is as follows: (6); In the formula It represents the change rate of the fragmentation index of region j. This indicator reflects the urgency of monitoring and its value range is [0,1); Indicates the flight distance from the current position of drone i to region j, in kilometers; , Represents the weight coefficient of the cost function, satisfying , and its value range is [0,1].

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

[0030] In S101, dynamic preprocessing of multi-source remote sensing data: To meet the requirements of efficient acquisition and fusion of "space-air-ground" multi-source data, a heterogeneous data pre-standardization technology for dynamic monitoring of terrace patterns is proposed. The goal is to solve the modeling errors caused by the complexity of data sources and spatio-temporal heterogeneity through spatio-temporal consistency correction, redundant information compression, and unified characterization of gradient terrain features. The specific content is as follows: (1) Unification of spatio-temporal reference and geometric correction: Due to the significant differences in the spatial reference and resolution of multi-source data acquisition platforms (satellites, drones, LiDAR), they need to be processed before use. First, standardize the data coordinates and resolution. All data need to be unified to the same coordinate system (WGS84 UTM) and spatial resolution (mesoscale core resolution 0.5m - 5m). Based on the dynamic resolution compensation algorithm, information loss can be avoided. The specific formula is as follows: (7); In the formula is the resolution of the data, corresponds to its acquisition weight, and weighted fusion of mixed-resolution data is performed. Secondly, perform geometric correction on the data. For drone images, use SIFT feature matching combined with affine transformation (refer to Table 1) to eliminate the stitching offset caused by unstable flight attitude; for LiDAR and satellite elevation data, establish a residual correction model to eliminate the local slope inversion error.

[0031] (2) Feature enhancement of multi-source data: Focus on problems such as blurred terrace boundaries and terrain occlusion, and enhance key information features in combination with the characteristics of data sources. First is spectral cloud removal and radiation normalization. For satellite optical images (such as WorldView-3), use the multi-temporal cloud mask synthesis method to dynamically eliminate invalid pixels based on the normalized cloud index (NDCI). The specific formula is as follows: (8); In the formula represents the reflectance in the near-infrared band (dimensionless, with a value range of [0,1]), which is used to detect the difference between vegetation and clouds; It represents the reflectivity of the short-wave infrared band (dimensionless, ranging from [0,1]), which is sensitive to clouds and fog; if NDCI>0.2, it is judged as an effective cloud-free pixel (empirical threshold). Secondly, the multispectral data of the drone is corrected for the radiant brightness through the solar altitude angle model to eliminate the difference in illumination. Finally, the terrain features are enhanced. The LiDAR point cloud processing uses the Poisson surface reconstruction algorithm to generate a 0.1m accuracy DEM (refer to Table 1), and the steep slope area (>25°) is shadow compensated. After the slope is corrected in Formula 4, the rasterized terrain parameters (slope, roughness) are output; superpixel segmentation (SLIC algorithm) combined with Canny edge detection is applied to enhance the geometric features of terraced ridges and suppress interference from crop cover.

[0032] (3) Redundant information compression and dynamic quality assessment: Redundancy was eliminated in the previous data collection phase, but the system needs to propose optimization strategies for the redundancy problem of repeated coverage areas of multi-source data. The first is dynamic clipping based on the Herfindahl Index (HI). According to the redundancy index defined in Formula 2, redundant data is automatically eliminated. 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, it is adaptive grid quality control. The grid size in Formula 3 is used to check the data integrity (coverage efficiency) in each unit in blocks. ), triggering the drone supplementary sampling command for the insufficient area (refer to formula 6). Finally, the data priority is dynamically updated, and the weight parameters in formula 1 are reversely optimized based on the pre-processed data quality feedback (such as registration error, signal-to-noise ratio). , forming a closed loop of "collection-preprocessing-optimization" (for example, when the registration error in the plain area exceeds the threshold, Improved from 0.15 to 0.18 to accelerate the elimination of old data).

[0033] 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 feedback of quality, 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.

[0034] In S102, the specific steps of extracting terrace distribution data from a typical area using a multimodal boundary refinement segmentation algorithm are as follows: By integrating optical texture, terrain curvature and temporal vegetation characteristics, the weak edge (crop cover) and fragmentation (shadow interference) problems of ridges are solved, and continuous vector boundaries of terrace patches are generated, as follows: Step 1: Input feature definition and channel processing; The input data includes three types of features. The first is the optical image feature, which includes the RGB three bands and the red edge band. The phase selects the low crop coverage period (such as early spring when there is no cultivation) to maximize the exposure of the terrace ridge structure; the second is the terrain parameter feature, which includes slope and terrain curvature. The specific formula for slope is as follows: (9); In the formula It is the slope angle calculated based on DEM, reflecting the degree of surface inclination, ranging from [0°, 90°]; z represents the elevation value (unit: meter), which comes from the digital elevation model; and are the elevation gradients of DEM in the x and y directions respectively; the terrain curvature is used to enhance the contour features of the terraced fields. The specific formula is as follows: (10); Where C represents the terrain curvature (unit: curvature per unit length, usually 1 / meter). Positive curvature corresponds to raised terraces, and 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, the interannual NDRE (Normalized Difference RedEdge) time series, are calculated every 16 days. The specific calculation formula is as follows: (11); In the formula Indicates the reflectivity in the near-infrared band (wavelength of about 770-900 nm). Strong reflection areas usually correspond to vegetation canopies. It represents the reflectance of the red edge band (wavelength of about 700-740 nm), which is sensitive to chlorophyll content and is used to distinguish vegetation types. The NDRE value range is [-1,1]. The higher the value, the denser the vegetation coverage or the higher the chlorophyll content.

[0035] Step 2: Dynamic weight multimodal segmentation network; MMF-SegNet (Multi-Model Fusion Segmentation Network) is designed. First, branch features are extracted. In the optical branch, a dilated residual block is used to capture multi-scale textures. The convolution kernel dilation rate is set to 1, 2, and 4, and the feature map is output. Secondly, the terrain branch inputs the 2-channel data of slope and terrain curvature, uses lightweight U-Net encoding, and outputs terrain features Finally, for time series analysis, we input the NDRE time series (time step T=24), use 1D Conv-LSTM to extract time-dependent features, and output the vector FNDRE. Then, we integrate the feature modules and introduce Terrain-Guided Attention (TGA) to dynamically integrate optical and terrain features: (12); (13); (14); In the formula Represents the terrain roughness index (Roughness Index), which is a composite of slope and terrain curvature; Refers to the learnable parameters, initialized to 0.8, and optimized by back-propagation; Represents the Sigmoid function, constraining the weight to [0,1]; weight coefficient Dynamically adjusted by terrain complexity; high RI areas (steep and broken terrain) have reduced optical weight ( ), enhancing terrain features to suppress shadow interference. Finally, the temporal features are inserted to transform the NDRE features With fusion features Splicing, through 3×3 convolution fusion, output the segmentation probability map P ∈ [0,1] of the terrace area.

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

[0037] In S102, the spatial-temporal evolution characteristics of terraces are quantitatively extracted to dynamically characterize the 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, as follows: Step 1: Interannual change detection and event classification; First, the terrace area change rate is calculated. The specific formula is as follows: (15); (16); in the formula represents the total area of ​​terraces in the study area in year t, represents the area of ​​the kth terrace patch, and K represents the total number of terrace patches; represents the total area of ​​terraces in the study area in year t+1 (unit: square meters or hectares); represents the annual change in the area of ​​terraced fields. And it has been growing for three consecutive years, marked as a stable expansion area. (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: (17); in the formula Represents the category label of pixel (i, j) from time t to t+1; Represents the binary mask of year t, 0 represents non-terraced fields and 1 represents terraced fields.

[0038] Step 2: Analysis of spatial pattern index evolution; First, select the appropriate spatial pattern index, as follows: Fragmentation Index (FI): (18); in the formula Indicates the total number of terrace patches; represents the average plaque area; The larger it is, the more fragmented the terrace landscape is and the greater the degree of disturbance.

[0039] Shape Stability Index (SSI): (19); in the formula represents the perimeter of the kth terrace patch (unit: meter); 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).

[0040] Center of gravity migration trajectory: The centroid coordinates of terrace distribution in each year (Xt, Yt): (20); (twenty one); in the formula represents the area of ​​the kth terrace patch; ( , ) represents the centroid coordinates of the kth patch (based on the geometric center of the vector polygon); , ) represents the total centroid coordinates of the terrace distribution in year t.

[0041] Annual migration distance: (twenty two); in the formula Represents 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).

[0042] Step 3: Generation of spatiotemporal evolution map; This step mainly generates visualization maps. First, a density map of terrace expansion / contraction is generated year by year in the form of a heat map sequence. The kernel density estimation bandwidth is 1km to express spatial aggregation. The second is the change trajectory map, which constructs a change path network map (G=(V, E)). Node V is the centroid of the change patch, edge E connects the centroids of adjacent years, and the weight is the migration distance. Finally, a three-dimensional space-time cube model is constructed. According to the pixel stacking inter-annual two-dimensional mask Mt, the evolution trend of the vertical time axis is displayed through voxel rendering.

[0043] In addition, technical verification and accuracy evaluation: mainly verify and evaluate the accuracy of 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, Recall=0.92➡F1=0.904). In addition, the edge similarity (ES) indicator is used to measure the spatial alignment accuracy of 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: (twenty three); In the formula, EDT(P) represents the Euclidean distance transform map of the prediction 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.

[0044] The next step is to verify the results of the spatiotemporal evolution analysis. 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.

[0045] The above technical solution mainly focuses on the accurate extraction of terrace distribution and the characterization of its spatiotemporal dynamics. Through the linkage of multimodal deep learning and spatial statistical models, a comprehensive quantitative analysis of the terrace evolution process can be achieved.

[0046] In S103, based on the previously extracted terrace distribution data (vector boundaries, spatiotemporal change trajectories), combined with multi-source remote sensing and ground observation data, an evaluation framework of "spatiotemporal heterogeneity modeling-key process coupling-multidimensional service coordination" is proposed to quantitatively characterize the terraces' ability to quantify soil conservation services, evaluate water conservation services, evaluate carbon sequestration services, evaluate food supply services, and evaluate biomass energy supply services. Emphasis is placed on the dynamic model parameters, spatial heterogeneity response, and multi-service coordinated optimization. The specific contents are as follows: (1) Quantification of soil conservation services; 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 the improved RUSLE model (revised soil and water conservation equation). Compared with the traditional method, the spatial resolution of the parameters is improved by dynamically inverting vegetation coverage and terrace engineering factors through the red edge band. The specific formula is as follows: ; (twenty four); (25); and They represent the potential soil erosion and actual soil erosion in the region, respectively, and the difference between the two represents the amount of soil retained in the region by terrace engineering; represents rainfall erosivity MJ·mm / (hm²·h·yr), which is based on local rain gauges or CHIRPS remote sensing rainfall data to calculate the total annual rainstorm erosivity; K represents soil erodibility t·hm²·h / (hm²·MJ·mm), which is parameterized by the EPIC model using soil texture maps; Represents the terrain factor, and generates the slope and slope length based on the satellite elevation data of the terrace distribution area. The formula is: (26); (27); (28); 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 it into terraced fields; m is the slope length factor index; is the slope angle of the slope; Represents the terrace vegetation cover management factor (dimensionless), based on the dynamic inversion of the red edge band (linked with time series NDRE): (29); 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; It represents the terrace engineering factor (dimensionless), and the assignment rule is: 0.1 for horizontal terrace, 0.15 for reverse slope terrace, 0.2 to 0.3 for earthen ridge terrace (assigned according to the ridge height through drone image classification).

[0047] (2) Water conservation service assessment; The goal is to evaluate the ability of terraces to regulate runoff and enhance groundwater recharge. The technical solution is to construct a terrace distributed water conservation model (TWHM), coupling the terrace structure with the spatial heterogeneity of soil hydraulic parameters. The specific formula is as follows: (30); in the formula represents the annual precipitation of grid cell i (mm), obtained through interpolation of meteorological stations or GSMaP satellite precipitation products; represents the area of ​​grid cell i; It represents the actual evapotranspiration of the terrace in grid cell i, (mm), which is inverted based on the SEBAL model or MODIS ET product (MOD16), combined with the correction of the water-saving effect of the terrace: (31); in the formula represents the interception coefficient of terraced fields in the wet season (calibrated by ground leakage experiment), It is the crop canopy moisture index in the wet season; represents the surface runoff depth of the terrace in grid cell i (unit: mm), calculated using the SCS-CN model: (32); (33); CN value adjustment: CN value of terraced area is reduced to 75-85 (depending on the water storage capacity of terraced field ridges); It represents the soil infiltration enhancement coefficient, which is dynamically corrected based on soil type (such as clay loam η=1.2, sandy soil η=0.8) and terrace tillage layer thickness (LiDAR inversion).

[0048] (3) Carbon sequestration service assessment; The goal is to measure the photosynthetic carbon fixation of terraced crops and the accumulation of soil organic carbon. The specific technical solution is to build a hierarchical carbon sink model. The first part is vegetation carbon fixation (NPP). The specific formula is as follows: (34); in the formula (MJ / m²) represents absorbed photosynthetically active radiation, which is calculated by inverting FPAR (vegetation coverage) + solar radiation data through the red edge band; represents the light energy utilization efficiency (gC / MJ), which is assigned values ​​according to the crop type (1.8 for rice, 2.1 for corn, and 1.5 for wheat); is the crop rotation coefficient (such as double-season rice f=1.5). The second part is soil carbon sequestration, and the specific calculation formula is as follows: (35); in the formula represents the carbon input of crop litter (kg / ha / yr), which is inverted by biomass (NPP) and root-to-shoot ratio; represents the initial soil organic carbon content (%), based on the inversion of the Sentinel-2 shortwave infrared band; , Indicates the rate of organic carbon transformation and mineralization, depending on the tillage method (such as no-till =0.6); BD represents soil bulk density (g / cm³), which is spatially interpolated from measured data of soil profiles; D represents the depth of cultivated layer (m), which is obtained by combining DEM with remote sensing classification results of cultivated activities.

[0049] (4) Food supply service assessment; 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: (36); in the formula It represents the harvest index (crop economic allocation ratio, such as rice HI=0.45); Indicates 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: (37); in the formula represents the drought stress loss rate, which is inverted by the vegetation temperature condition index; It indicates the yield reduction rate caused by soil erosion and is negatively correlated with soil conservation.

[0050] (5) Evaluation of biomass energy supply services; The goal is to evaluate the potential for utilization of biomass energy such as terraced field straw and oil crops. The technical solution adopted 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: (38); in the formula represents the yield of the i-th crop (kg / ha), and n represents the number of crop types; represents the grass-to-grain ratio; It represents the energy conversion efficiency. Then the biodiesel potential of oil crops is calculated. The specific formula is as follows: (39); in the formula It indicates the proportion of crop planting area to terrace area; Indicates oil yield per unit area (L / ha), referring to FAO crop parameter database; It represents the conversion rate of transesterification reaction (typical value 0.95).

[0051] The above technical solution is mainly based on the temporal and spatial distribution data of terraces, and proposes an ecosystem service assessment framework of "dynamic parameter modeling-multi-process coupling-multi-dimensional collaborative optimization", focusing on the five core functions of soil conservation service quantification, water conservation service assessment, carbon fixation service assessment, food supply service assessment and biomass energy supply service assessment. The high-resolution remote sensing data drives the dynamicization of model parameters, breaking through the limitations of traditional static parameters, and providing key technical support for the ecological protection and sustainable management of terraces.

[0052] In S103, in order to reveal the interaction (trade-off / synergy) between the service functions of terraced field ecosystems, a comprehensive technical solution of "spatiotemporal heterogeneity analysis-driving factor decoupling-multi-objective optimization mapping" is proposed. By coupling geostatistical models, machine learning causal inference and spatial explicit optimization algorithms, the nonlinear relationship between soil conservation service quantification, water conservation service assessment, carbon sequestration service assessment, food supply service assessment and biomass energy supply service assessment is analyzed from multiple scales. The specific methods are as follows: (1) Framework for analyzing trade-off synergy relationships; The goal is to quantify the trade-off intensity and synergy gain between different ecosystem services of terraces and to identify the main driving factors. The technical route includes three aspects: the first is data input, based on the spatial distribution map of the five ecosystem service capabilities generated in the previous chapters, the vector map of terrace distribution, multi-temporal environmental variables (rainfall, temperature, topography), and the human activity interference intensity index; the second is relationship analysis, using multi-model cascade analysis (correlation statistics-causal inference-spatial clustering); the third is result output, generating a trade-off synergy intensity matrix, a spatial heterogeneity heat map, and a multi-objective optimization solution set.

[0053] (2) Constructing a framework for analyzing trade-off synergy relationships; Step 1: Statistical correlation analysis (base layer); Calculate the Pearson correlation coefficient and partial correlation coefficient between services based on grid cells: (40); in the formula Represents the correlation coefficient between service j and service k. A positive value represents a synergistic relationship, while a negative value represents a trade-off relationship. In addition, significance correction is performed through Bonferroni correction to eliminate false positives from multiple tests, and only strong correlations with p < 0.01 are retained.

[0054] Step 2: Machine learning causal identification (core layer); A two-stage causal mechanism decoupling model (TCMD) is proposed: Phase 1 (driving factor screening): Random forest feature importance ranking is used to screen the main controlling factors that affect the trade-offs between services (such as slope, vegetation coverage, and terrace ridge density).

[0055] 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: (41); in the formula represents natural factors (climate, soil thickness, etc.), and instrumental variables (such as altitude) are added to control endogeneity; It represents human activities (fertilizer application, terrace maintenance frequency), which are inverted through night light data and farmer survey data.

[0056] Step 3: Calculation of trade-off synergy index (output layer); Define the Spatio-Temporal Trade-off Index (STTI): (42); in the formula represents the trade-off index between services j and k at coordinate x at time t (negative value = trade-off, positive value = synergy); represents the marginal change rate between services based on Geographically Weighted Regression (GWR); represents the maximum possible trade-off distance in the study area (taking all 99% quantile); It represents the trade-off distance, which reflects the Euclidean distance (degree of deviation) between the current state and the theoretical optimal state. It can be specifically expressed by the following formula: (43); In the formula represents the expected value of the mth service under the theoretical optimal state (determined by the production possibility frontier PPF).

[0057] Step 4: Spatially explicit collaborative optimization (application layer); A multi-objective particle swarm optimization model (MOPSP) was constructed. The objective function was to maximize the total synergistic gain of the five ecosystem services. The specific formula is as follows: (44); In the formula, Z represents the total synergy gain objective function; represents the synergy gain between services j and k (if , then it is the actual value; otherwise it is 0); represents the synergistic weight of service jk (reflecting the management priority). The constraints include two aspects: one is the threshold of ecosystem services that cannot be reduced, such as soil conservation ≥ 80% of the current status; the other is the upper limit of the transformation of terrace projects, such as slope ≤ 25°.

[0058] After the above steps, we can obtain the trade-off synergy map of the five ecosystem services provided by terraces, and further form a spatial hotspot map (such as the "high grain-low water conservation" conflict area) and spatial optimization priority areas (low STTI value areas). In addition, by summarizing the typical trade-off patterns of ecosystem services of terraces in different climate regions (such as the strong trade-off of "carbon sequestration-grain" in northern dryland terraces), we can provide relevant basis for the adjustment of cross-regional policies.

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

[0060] In summary, this embodiment, through system integration and innovative transformation, at the public technology level, adopts the SIFT feature matching of multi-source remote sensing data and the Poisson reconstruction of LiDAR point cloud to achieve geometric alignment, optimizes the open source fragmentation (FI) and stability index (SSI) extraction modules of landscape ecology, and constructs a basic evaluation framework based on classic ecological models (RUSLE soil erosion equation, SCS-CN runoff model, light energy utilization rate NPP model) to ensure the scientificity and verifiability of the method system. In terms of innovation, the "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 of multi-spectral / high-resolution data collaborative scheduling, 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 spatiotemporal explicit trade-off collaborative analysis technology (STTI) was independently designed, integrating dual machine learning causal inference and MOPSO multi-objective optimization algorithm, analyzing the natural-human driving relationship in terrace multi-ecosystem services, and supporting the generation of differentiated regulation 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.

[0061] The technical effects are as follows: (1) Dynamic optimization of acquisition technology solves the contradiction between mesoscale monitoring efficiency and coverage; through the "four-dimensional dynamic weight allocation algorithm", a balance model of spatial resolution-coverage efficiency-timeliness-redundancy is constructed, combined with event-driven adaptive grid division and feedback mechanism to realize intelligent scheduling of multi-source data. Based on the Herfindahl index, data redundancy is dynamically optimized, the priority weight of drone supplementation is adjusted according to the complexity of the terrain, and invalid acquisition is eliminated through adaptive grid. 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%, which effectively solves the contradiction between "high-resolution local coverage" and "lack of timeliness of wide-area observation" in the background technology.

[0062] (2) Heterogeneous data preprocessing enhances the recognizability of terrace features; the geometric and spectral errors of multi-source data are eliminated by using the unification of spatiotemporal benchmarks, NDCI dynamic declouding, Poisson surface reconstruction and residual slope correction technology. Through radiation 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 NDCI declouding effective pixel retention rate 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.

[0063] (3) Multimodal segmentation network improves the accuracy of terrace boundary recognition; the terrain guided attention mechanism TGA is designed, and the driven MMF-SegNet network integrates optical texture, terrain curvature and time series NDRE features to construct a dynamic weight multimodal segmentation model. The TGA mechanism dynamically suppresses shadow interference through the terrain roughness index (RI), reduces the optical weight to 0.3 in steep slope areas (RI>15), and significantly enhances the contour features of the ridge. 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.

[0064] (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 dynamically assign the vegetation management factor C and the CN value of the ridge water storage capacity classification, and a spatial adaptive evaluation framework is constructed. 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²) is improved from 0.68 to 0.83, solving the key limitation of "static ecosystem service evaluation" in the background technology.

[0065] (5) Causal optimization coupling technology enhances the effectiveness of collaborative decision-making in services; proposes a spatiotemporal explicit trade-off index (STTI) and a dual machine learning causal inference-MOPSO optimization cascade model 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 coverage, and terrace ridge density through particle swarm optimization. In typical conflict areas (such as "high grain-low water conservation" areas), the collaborative gain is increased by 35%, and the efficiency of generating policy adaptation solutions 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.

[0066] In summary, the system has systematically solved 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.

[0067] In this specification, each embodiment is described in a progressive manner, and each embodiment focuses on the differences from other embodiments. The same or similar parts between the embodiments can be referred to each other. For the device disclosed in the embodiment, since it corresponds to the method disclosed in the embodiment, the description is relatively simple, and the relevant parts can be referred to the method part.

[0068] The above description of the disclosed embodiments enables one skilled in the art to implement or use the present invention. Various modifications to these embodiments will be 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 will not be limited to the embodiments shown herein, but rather to the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A method for quantifying the evolution of terrace patterns and its ecosystem service trade-off, characterized in that: include: Construct a multi-source remote sensing information coordination framework, identify frequently changing areas, pre-process multi-source remote sensing data, and obtain a standard input data set through the unification of spatiotemporal benchmarks, key feature enhancement, redundant information compression, and quality dynamic feedback optimization; Based on the standard input data set, a multimodal boundary refinement segmentation algorithm is used to extract the distribution of terraces in typical areas and their interannual expansion / contraction dynamics to characterize the quantitative temporal and spatial evolution. Through the linkage of multimodal deep learning and spatial statistical models, a comprehensive quantitative analysis of the terrace evolution process is achieved. Based on the extracted terrace distribution data, combined with the standard input data set and ground observation data, an evaluation framework is proposed to quantitatively characterize the terraces in soil conservation service quantification, water conservation service assessment, carbon sequestration service assessment, food supply service assessment and biomass energy supply service assessment, and to conduct a synergistic analysis of the trade-offs among the above typical ecosystem services.

2. A method for quantifying the evolution of terrace patterns and its ecosystem service trade-off according to claim 1, characterized in that: The specific steps to build a multi-source remote sensing information collaborative framework are as follows: Using a dynamic weight allocation algorithm, mathematical modeling is used to achieve a global optimal balance between spatial resolution, coverage efficiency, and time cost. The Herfindahl index is introduced to measure data diversity redundancy to ensure that the weight allocation is comprehensive and avoids information overload. When necessary, the study area is adaptively gridded to ensure good image data connection. Construct a "macro-representation-process capture-micro-calibration" framework to carry out complementary fusion of multi-source data; The mesoscale adaptive data collection command feedback mechanism is used to improve data collection efficiency and identify frequently changing areas.

3. A method for quantifying the evolution of terrace patterns and weighing their ecosystem services according to claim 2, characterized in that: The weight calculation model of the dynamic weight allocation algorithm is as follows: ； ； in, It represents the priority weight of data source i, which determines its collection order in scheduling and is dimensionless; Indicates the spatial resolution of data source i, in meters; Represents the data source The coverage efficiency is the coverage area per unit time, in km2 / day; represents the time-dependent attenuation coefficient; Indicates the time delay after data collection, in days; is the Herfindahl index, which measures the redundancy of similar data sources. is the ratio of the coverage area of ​​data source j to the total overlapping area; k represents the total number of data sources.

4. The method for quantifying the evolution of terrace pattern and its ecosystem service trade-off according to claim 1, 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; Building a dynamic weight multimodal segmentation network based on the input features to obtain a segmentation probability map of the terraced field area; The segmentation probability map is subjected to threshold segmentation, and concave points are detected on the contour of the terraced field patches. The concave points are connected, and the obtained binary mask is converted into a polygonal vector layer, and the final terraced field distribution vector layer is output.

5. A method for quantifying the evolution of terrace patterns and weighing its ecosystem services according to claim 4, 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 the dilation residual module to capture multi-scale textures and output the optical feature map. Secondly, in the terrain branch, we input the 2-channel data of slope and terrain curvature splicing, use lightweight U-Net encoding, and output terrain features. Finally, for time series analysis, we input the NDRE time series, use 1D Conv-LSTM to extract time-dependent features, and output a vector. For the fusion of feature modules, we introduce terrain to guide attention, and dynamically fuse optical and terrain features to obtain fused features: ； ； ； In the formula, Represents the terrain roughness index, which is a composite of slope and terrain curvature; Refers to learnable parameters, optimized by back-propagation; Represents the optical image characteristics, Represents terrain parameter characteristics; is the Sigmoid function; weight coefficient Dynamically adjusted by terrain complexity; It is the slope angle calculated based on DEM, reflecting the degree of surface inclination; C represents the terrain curvature; Finally, the vegetation temporal characteristics With fusion features Splicing, through convolution fusion, output the segmentation probability map of the terrace area.

6. The method for quantifying the evolution of terrace patterns and weighing its ecosystem services according to claim 1, characterized in that: The specific steps for quantifying the spatiotemporal evolution characteristics of terraces and extracting the dynamic representation of interannual expansion / contraction are as follows: Calculate the rate of change of terrace area, count the number of expanded / reduced terrace patches, their spatial distribution and area proportion; The evolution of spatial pattern index, the selected spatial pattern index, specifically includes fragmentation index, morphological stability index, centroid migration trajectory, and interannual migration distance; To generate the spatiotemporal evolution map of terraces, we first generate a density map of terrace expansion / contraction year by year in the form of a heat map sequence; The second is the change trajectory diagram, which constructs the change path network diagram ( ), node V is the centroid of the change patch, edge E connects the centroids of adjacent years, and the weight is the migration distance; finally, the three-dimensional space-time cube model is constructed, and the inter-year two-dimensional mask is stacked according to the pixels. The evolution trend of the vertical time axis is displayed through voxel rendering to generate a visual atlas.

7. A method for quantifying the evolution of terrace patterns and weighing its ecosystem services according to claim 6, characterized in that: The calculation of the terrace area change rate, statistics of the number of expanded / reduced terrace patches and their spatial distribution and area proportion, the specific steps are as follows: The specific formula for calculating the change rate of terrace area is as follows: ； ； in the formula represents the total area of ​​terraces in the study area in year t, represents the area of ​​the kth terrace patch, and K represents the total number of terrace patches; represents the total area of ​​terraces in the study area in year t+1, in square meters or hectares; represents the annual change in the area of ​​terraced fields; if And it has been growing for three consecutive years, marked as a stable expansion area. , determined as a damage event; Perform spatial overlay analysis and detect the change area pixel by pixel. The formula is as follows: ； in the formula Represents the category label of pixel (i, j) from time t to t+1; Represents the binary mask of year t, 0 represents non-terraced fields and 1 represents terraced fields.

8. The method for quantifying the evolution of terrace pattern 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 field ecosystems are as follows: Quantification of soil conservation services: The improved RUSLE model is used, and the specific formula is as follows: ； ； ； in the formula and They represent the potential soil erosion and actual soil erosion in the region, respectively, and the difference between the two represents the amount of soil retained in the region by terrace engineering; represents rainfall erosivity MJ·mm / (hm²·h·yr), which is based on local rain gauges or CHIRPS remote sensing rainfall data to calculate the total annual rainstorm erosivity; K represents soil erodibility t·hm²·h / (hm²·MJ·mm), which is parameterized by the EPIC model using soil texture maps; Represents the terrain factor, and generates the slope and slope length based on the satellite elevation data of the terrace distribution area. The formula is: ； ； ； in the formula is the slope length, unit: meter. After extracting the watershed watershed, the actual slope length is calculated by dividing the terraced fields. m is the slope length factor index; is the slope angle of the slope; Represents the terrace vegetation cover management factor, dimensionless, based on the dynamic inversion of the red edge band: ； It indicates the maximum NDRE value during the crop growing season; represents the terrace engineering factor, dimensionless; Water conservation service evaluation: A distributed water conservation model for terraces is constructed, coupling the spatial heterogeneity of terrace structure and soil hydraulic parameters. The specific formula is as follows: ； ； ； ； in the formula represents the annual precipitation of grid cell i, in mm, obtained through interpolation from weather stations or GSMaP satellite precipitation products; represents the actual evapotranspiration of the terrace in grid unit i, mm; Represents the surface runoff depth of the terrace in grid unit i, unit: mm; represents the area of ​​grid cell i; It represents the soil infiltration enhancement coefficient, dimensionless, and is dynamically corrected based on soil type and terraced field cultivation layer thickness; ETnatural represents the evapotranspiration of natural vegetation, mm; It represents the interception coefficient of terraced fields in the wet season, which is calibrated by ground leakage experiment. It is the canopy moisture index of crops in the wet season; it is calculated using the SCS-CN model: the CN value represents the runoff curve number, which is dimensionless; represents annual precipitation, mm; Indicates the potential maximum retention of water, mm; Carbon sequestration service assessment: Construct a hierarchical carbon sink model, first of all, vegetation sequestration, the specific formula is as follows: ； in the formula Indicates absorbed photosynthetically active radiation, MJ / m²; It represents the light energy utilization rate, gC / MJ; is the crop rotation coefficient; the second is soil carbon sequestration, the specific calculation formula is as follows: ； in the formula represents the carbon input of crop litter, kg / ha / yr; represents the initial soil organic carbon content, %; , It represents the rate of organic carbon transformation and mineralization, and is assigned according to the tillage method; BD represents soil bulk density, g / cm³; D represents the depth of the tillage layer, m; Food supply service assessment: Using the mixed yield model, first calculate the potential yield. The specific formula is as follows: ； in the formula represents net primary productivity, gC / m2 / year; represents the harvest index, dimensionless; Indicates the carbon content-yield conversion coefficient, kg / gC; combined with the actual yield correction, the specific formula is as follows: ； in the formula represents the drought stress loss rate, which is inverted by the vegetation temperature condition index; It indicates the yield reduction rate caused by soil erosion, which is negatively correlated with soil conservation; Biomass energy supply service assessment: Use the multi-source biomass inversion model to calculate the energy content of relevant crop straw. The specific formula is as follows: ； in the formula represents the yield of the i-th crop, kg / ha, and n represents the number of crop types; represents the grass-to-grain ratio; It represents the energy conversion efficiency; then the biodiesel potential of oil crops is calculated, and the specific formula is as follows: ； in the formula It represents the proportion of crop planting area to terrace area, dimensionless; A is the total area of ​​terrace, ha; It indicates the oil yield per unit area, L / ha; It represents the conversion rate of transesterification reaction.

9. The method for quantifying the evolution of terrace pattern and its ecosystem service trade-off according to claim 1, 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 the generated data; second, relationship analysis, using multiple models to couple related data information; third, result output, generating a trade-off relationship matrix, a 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 base 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 control 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 constructed. 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 the transformation of terrace projects.

Images