Low-carbon remediation decision-making method based on the spatiotemporal differentiation of ecosystem carbon storage.

Abstract

This invention provides a land zoning low-carbon remediation decision-making method based on the spatiotemporal differentiation of ecosystem carbon storage, relating to the field of computer technology. First, it acquires a target geographic area's remote sensing image set and a ground sample plot survey set. Then, it extracts the spatiotemporal differentiation features of carbon storage from the target remote sensing image set and the ground sample plot survey set, obtaining spatial and temporal differentiation maps. Next, it delineates land management zones based on the spatial and temporal carbon storage differentiation maps, obtaining a set of spatial units for each zone. Then, it generates a set of low-carbon remediation instructions for zoning-differentiated land based on driving force analysis and matching remediation items. Finally, it sends the set of low-carbon remediation instructions for zoning-differentiated land to the land management terminal to trigger relevant operations. This invention can accurately achieve low-carbon remediation of land zones and enhance the carbon sink function of land ecosystems.

Description

Technical Field

[0001] This invention relates to the field of computer technology, and more specifically, to a land zoning low-carbon remediation decision-making method based on the spatiotemporal differentiation of ecosystem carbon storage. Background Technology

[0002] As a crucial carbon reservoir, land's spatial and temporal variations in carbon storage have a critical impact on the global carbon cycle. Traditional land remediation decision-making methods are significantly inadequate in addressing changes in ecosystem carbon storage.

[0003] On the one hand, past land remediation decisions have relied heavily on a single data source, such as ground survey data or simple remote sensing data, making it difficult to obtain comprehensive and accurate spatiotemporal information on carbon storage in land ecosystems. While ground surveys can obtain relatively accurate local data, their coverage is limited and cannot reflect changes in carbon storage over a large area; while simple remote sensing data may lack necessary spatial coordinate information and detailed biological and soil parameters, leading to inaccurate estimates of carbon storage.

[0004] On the other hand, existing land remediation decisions lack in-depth analysis of the spatiotemporal differentiation characteristics of carbon storage and fail to implement targeted zoning management based on the carbon storage change characteristics of different regions. When formulating remediation strategies, a "one-size-fits-all" approach is often adopted, without fully considering the needs for carbon sink maintenance, carbon source management, and carbon sink enhancement in different regions. This results in poor remediation outcomes, an inability to effectively enhance the carbon sink function of land ecosystems, and difficulty in meeting the requirements of global carbon emission reduction and ecological sustainable development. Summary of the Invention

[0005] In view of the aforementioned problems, and in conjunction with the first aspect of the present invention, the present invention provides a land zoning low-carbon remediation decision-making method based on the spatiotemporal differentiation of ecosystem carbon storage, the method comprising: Acquire a target remote sensing image set and a ground sample plot survey set for the target geographic area within a historical continuous monitoring period. The target remote sensing image set includes multispectral image gratings and radar image gratings with spatial coordinate labels, and the ground sample plot survey set includes measured records of vegetation biomass and soil organic carbon with spatial coordinate labels. Carbon storage spatiotemporal differentiation features are extracted from the target remote sensing image set and the ground sample plot survey set to obtain a carbon storage spatial differentiation map reflecting the differences in the spatial distribution of carbon storage and a carbon storage temporal differentiation map reflecting the differences in the temporal variation of carbon storage. Land management zones are delineated on the carbon storage spatial differentiation map and the carbon storage temporal differentiation map to obtain a set of zoning spatial units consisting of the boundary coordinates of the carbon sink maintenance zone, the boundary coordinates of the carbon source control zone, and the boundary coordinates of the carbon sink enhancement zone. Based on the analysis of the driving forces of carbon storage changes in different management zones in the spatial unit set, the corresponding restoration entries in the preset land low-carbon restoration strategy library are matched for the boundary coordinate coverage areas of the carbon sink maintenance zone, the boundary coordinate coverage areas of the carbon source control zone, and the boundary coordinate coverage areas of the carbon sink enhancement zone, respectively, and a set of land low-carbon restoration instructions with zoned differences is generated. The set of instructions for low-carbon restoration of land with zoning differences is sent to the land management terminal to trigger vegetation structure adjustment, soil carbon pool protection, and ecological network connectivity enhancement operations for different zones.

[0006] Furthermore, this invention also provides a land zoning low-carbon remediation decision-making system based on the spatiotemporal differentiation of ecosystem carbon storage, comprising: A processor; a machine-readable storage medium for storing machine-executable instructions of the processor; wherein the processor is configured to execute the aforementioned land zoning low-carbon remediation decision-making method based on the spatiotemporal differentiation of ecosystem carbon storage by executing the machine-executable instructions.

[0007] Based on the above, by acquiring target remote sensing image sets and ground sample plot survey sets of the target geographical area within a continuous historical monitoring period, the spatiotemporal differentiation characteristics of carbon storage are extracted from the target remote sensing image sets and ground sample plot survey sets. The generated carbon storage spatial differentiation map and carbon storage temporal differentiation map can reflect the differences in carbon storage changes at different spatial locations and time series. Based on the carbon storage spatiotemporal differentiation map, land management zoning is carried out. The resulting zoning spatial unit set can accurately define the boundaries of carbon sink maintenance areas, carbon source control areas, and carbon sink enhancement areas. According to the analysis of the driving forces of carbon storage changes in different management zones, corresponding restoration items are matched from the preset land low-carbon restoration strategy library. The generated zoning-differentiated land low-carbon restoration instruction set is targeted and operable, which can effectively enhance the carbon sink function of different areas and realize the low-carbon restoration of land ecosystems. The land low-carbon restoration instruction set is sent to the land management terminal to trigger targeted vegetation structure adjustment, soil carbon pool protection, and ecological network connectivity enhancement operations, which helps to build a stable and efficient land ecosystem and improve the carbon sink capacity and stability of the ecosystem. Attached Figure Description

[0008] Figure 1 This is a schematic diagram of the execution process of the land zoning low-carbon remediation decision-making method based on the spatiotemporal differentiation of ecosystem carbon storage provided by the present invention.

[0009] Figure 2 This is a schematic diagram of exemplary hardware and software components of the land zoning low-carbon remediation decision system based on the spatiotemporal differentiation of ecosystem carbon storage provided by the present invention. Detailed Implementation

[0010] Figure 1This is a flowchart illustrating a land zoning low-carbon remediation decision-making method based on the spatiotemporal differentiation of ecosystem carbon storage, provided in one embodiment of the present invention. A detailed description follows.

[0011] The land zoning low-carbon remediation decision-making method based on the spatiotemporal differentiation of ecosystem carbon storage provided in this application can be applied to scenarios where land low-carbon remediation decisions are made for a key ecological functional area spanning administrative regions. For example, the target geographical area includes an upstream water conservation area, a central agricultural cultivation area, and a downstream urban development area, and there are different problems such as forest degradation, farmland carbon loss, and construction land expansion within the area.

[0012] Step S110: Obtain the target remote sensing image set and ground sample plot survey set for the target geographic area during the historical continuous monitoring period. The target remote sensing image set includes multispectral image gratings and radar image gratings with spatial coordinate labels, and the ground sample plot survey set includes measured records of vegetation biomass and soil organic carbon with spatial coordinate labels.

[0013] In this embodiment, firstly, multispectral remote sensing images and radar remote sensing images for a total of T monitoring time phases covering the target geographical area are acquired. After radiometric calibration and atmospheric correction, the multispectral remote sensing images of each monitoring time phase are used to obtain multispectral image gratings containing reflectance values ​​for the blue, green, red, and near-infrared bands. The radar remote sensing images of each monitoring time phase are then subjected to radiometric calibration, speckle filtering, and geocoding to obtain radar image gratings containing backscattering coefficient values. These multispectral image gratings and radar image gratings together constitute the target remote sensing image set.

[0014] Simultaneously, ground plot survey data were acquired for the target geographic area during the same continuous historical monitoring period. This ground plot survey data includes spatial coordinate markers for multiple fixed plots, each recording measured vegetation biomass values ​​(calculated using the quadrat harvesting method or allometric growth equation) and measured soil organic carbon values ​​(determined using the potassium dichromate oxidation method or elemental analysis). All spatial coordinate markers, vegetation biomass records, and soil organic carbon records for all fixed plots constitute the ground plot survey set.

[0015] Step S120: Extract the spatiotemporal differentiation features of carbon storage from the target remote sensing image set and the ground sample plot survey set to obtain a carbon storage spatial differentiation map reflecting the differences in the spatial distribution of carbon storage and a carbon storage temporal differentiation map reflecting the differences in the temporal variation of carbon storage.

[0016] Step S121: Perform vegetation index inversion on the multispectral image raster in the target remote sensing image set to obtain the canopy greenness index raster representing the intensity of vegetation photosynthesis, and interpret the backscattering coefficient of the radar image raster to obtain the canopy height index raster and the canopy water content index raster representing the three-dimensional structure of the ground surface.

[0017] For each monitoring phase of the multispectral image raster, the red band reflectance value R_red and the near-infrared band reflectance value R_nir are extracted, and the normalized vegetation index value is calculated pixel by pixel. The calculation formula is NDVI=(R_nir-R_red) / (R_nir+R_red) to obtain the canopy greenness index raster.

[0018] For each radar image grid at each monitoring time phase, backscattering coefficient values ​​for different polarization modes are extracted. A semi-empirical inversion algorithm based on the radiative transfer model is used to input the backscattering coefficient values ​​for horizontal transmission and vertical reception polarization into the water cloud model. The canopy height index is calculated as a function of the horizontal transmission and vertical reception polarization backscattering coefficient values ​​and the local incident angle, using the formula CHI=a*σ_HV*sqrt(cosθ_inc), where a is a preset empirical coefficient. The canopy water content index is calculated as the product of the polarization backscattering coefficient ratio and an empirical calibration coefficient, using the formula CWC=k_cwc*(σ_HV / σ_VV). The calculated canopy height index and canopy water content index values ​​are then filled into the corresponding pixels to generate canopy height index and canopy water content index grids.

[0019] Step S122: Input the canopy greenness index raster, the canopy height index raster, and the canopy water content index raster into the first feature layer of the carbon storage spatial heterogeneous modeling network. Extract the texture gradient of the canopy greenness index raster, the structural changes of the canopy height index raster, and the water heterogeneity information of the canopy water content index raster in parallel through the multi-scale convolution operator in the first feature layer. Then, stitch them together along the feature channel dimension to generate a target remote sensing fusion feature map. Input the target remote sensing fusion feature map into the spatial correlation layer of the carbon storage spatial heterogeneous modeling network. Calculate the covariance matrix between any points in the target remote sensing fusion feature map through the spatial self-attention operator in the spatial correlation layer. Then, generate a global spatial correlation feature map by weighted summation of the target remote sensing fusion feature map based on the covariance matrix.

[0020] The carbon storage spatial heterogeneous modeling network is a multi-branch deep learning network. Its first feature layer contains three parallel multi-scale convolutional subnetworks, processing canopy greenness index raster, canopy height index raster, and canopy water content index raster, respectively. Each multi-scale convolutional subnetwork uses three dilated convolutional kernels with different dilation rates (r1, r2, and r3) to extract features in parallel. For the canopy greenness index raster, three feature maps are obtained after three dilated convolution operations. These three feature maps are concatenated along the channel dimension and then passed through a 1x1 convolutional layer for channel dimensionality reduction, outputting a canopy greenness texture gradient feature map. Similarly, the canopy height index raster and canopy water content index raster are processed separately to obtain a canopy height structure variation feature map and a canopy water heterogeneity information feature map. These three feature maps are concatenated along the channel dimension to generate a target remote sensing fusion feature map.

[0021] The target remote sensing fusion feature map is input into the spatial association layer. The spatial association layer uses a nonlocal neural network module. First, the target remote sensing fusion feature map is transformed into a query feature map, a key feature map, and a value feature map through three 1x1 convolutional layers. The spatial dimensions of the three feature maps are all HxW, and the number of channels is C. The query feature map and the key feature map are reshaped into two-dimensional matrices of shape (HW, C). The covariance matrix A = softmax(Q*K^T) is calculated, and the shape of A is (HW, HW), where each element represents the spatial association weight between one location and another. The value feature map is reshaped into a two-dimensional matrix of shape (HW, C), and the weighted feature matrix O = A*V is calculated, and the shape of O is (HW, C). O is reshaped back to the spatial dimensions (H, W, C) to obtain the global spatial association feature map.

[0022] Step S123: Align the spatial coordinates of the measured vegetation biomass records and measured soil organic carbon records in the ground sample plot survey set with the coordinate grid of the global spatial association feature map, extract the first feature vector of the spatial coordinates of the corresponding measured vegetation biomass records in the global spatial association feature map, and extract the second feature vector of the corresponding measured soil organic carbon records.

[0023] Obtain the spatial coordinates of each fixed plot in the ground sample plot survey set. Map these coordinates onto the coordinate grid of the global spatial association feature map to determine the cell position of the coordinates in the feature map. Use bilinear interpolation to extract the C-dimensional feature vector at the cell position, which serves as the first feature vector for that plot. Similarly, for the plot coordinates corresponding to the measured soil organic carbon records, extract the corresponding second feature vector.

[0024] Step S124: Construct a vegetation biomass mapping multilayer sensing branch and a soil organic carbon mapping multilayer sensing branch. Input the first feature vector into the vegetation biomass mapping multilayer sensing branch for nonlinear transformation to generate a vegetation biomass spatial distribution map. Input the second feature vector into the soil organic carbon mapping multilayer sensing branch for nonlinear transformation to generate a soil organic carbon spatial distribution map.

[0025] The vegetation biomass mapping multilayer sensing branch is a three-layer fully connected network with an input layer dimension of C, hidden layers dimensions D1 and D2, and an output layer dimension of 1. The first feature vector of each sample plot is input into this branch, undergoes three nonlinear transformations (each followed by a ReLU activation function), and outputs the predicted vegetation biomass value for that sample plot location. Supervised learning is employed, using measured vegetation biomass values ​​from the ground sample plot survey set as labels to train this branch network. After training, the C-dimensional feature vector of each pixel in the global spatial association feature map is input into the trained vegetation biomass mapping multilayer sensing branch, calculating the predicted vegetation biomass value pixel by pixel, generating a spatial distribution map of vegetation biomass covering the entire target geographic area.

[0026] Similarly, the multilayer sensing branch for soil organic carbon mapping employs the same three-layer fully connected network structure. The second feature vector of each sample plot is input into this branch, and training is performed using measured soil organic carbon values ​​as labels. After training, a spatial distribution map of soil organic carbon is generated pixel by pixel.

[0027] Step S125: The spatial distribution map of vegetation biomass and the spatial distribution map of soil organic carbon are spatially superimposed using a grid overlay operator to generate a total spatial map of ecosystem carbon storage, and the gradient of carbon storage change between adjacent grids in the total spatial map of ecosystem carbon storage is calculated to generate a spatial differentiation map of carbon storage.

[0028] Input the spatial distribution maps of vegetation biomass and soil organic carbon into the raster overlay processor. For each pixel, calculate the total ecosystem carbon storage value C_total = B_map * C_factor_B + S_map * C_factor_S, where C_factor_B is the conversion factor for converting vegetation biomass dry weight to carbon mass, and C_factor_S is the conversion factor for converting soil organic carbon mass to carbon mass. Generate the spatial total ecosystem carbon storage map.

[0029] Calculate the spatial gradient of the total spatial map of ecosystem carbon storage. For each cell, calculate the carbon storage difference ΔC_east = C_map(i+1,j) - C_map(i,j) with its eastern neighbor and the difference ΔC_north = C_map(i,j+1) - C_map(i,j) with its northern neighbor. Then calculate the gradient magnitude G_mag = sqrt(ΔC_east^2 + ΔC_north^2). Fill the gradient magnitude of each cell into the corresponding position to generate a spatial differentiation map of carbon storage.

[0030] Step S126: Extract the initial carbon storage spatial total map corresponding to the starting time point and the final carbon storage spatial total map corresponding to the ending time point within the historical continuous monitoring period. Calculate the grid-by-grid difference between the initial carbon storage spatial total map and the final carbon storage spatial total map using a time-series change operator to generate a carbon storage time change raster map.

[0031] From historical continuous monitoring periods, the spatial total carbon storage map of the ecosystem corresponding to the starting monitoring phase t_start and the ending monitoring phase t_end are extracted. Using a time-series change operator, the temporal change in carbon storage ΔC_time = C_map_end(i,j) - C_map_start(i,j) is calculated for each pixel. The change for each pixel is then filled into the corresponding location to generate a raster map of the temporal change in carbon storage.

[0032] Step S127: Determine the direction of change of the carbon storage time change raster, mark the change polarity of each grid in the carbon storage time change raster, the change polarity includes carbon storage accumulation state and carbon storage loss state, and generate a carbon storage time differentiation map based on the combination characteristics of the change polarity of continuous grids.

[0033] For each pixel in the raster map of carbon storage time changes, determine its change polarity. Set a change discrimination threshold T_change. If ΔC_time(i,j)>T_change, mark the pixel as a carbon storage accumulation state; if ΔC_time(i,j)<-T_change, mark it as a carbon storage loss state; if |ΔC_time(i,j)|≤T_change, mark it as a carbon storage stable state.

[0034] All consecutive adjacent pixels with the same polarity of change are extracted to form connected regions. For each connected region, the average polarity intensity of its internal pixels is calculated. The spatial distribution boundaries of the connected regions for carbon storage accumulation and carbon storage loss are extracted to generate a carbon storage temporal differentiation map, which records the distribution of regions with different polarities of change in the form of planar vector features or raster regions.

[0035] Step S128: Input the carbon storage spatial differentiation map into the spectrum edge enhancement filter for convolution filtering, perform normalized numerical stretching on the carbon storage spatial differentiation map and carbon storage temporal differentiation map after convolution filtering, and map the numerical range of the carbon storage spatial differentiation map and the carbon storage temporal differentiation map after convolution filtering to a preset display dynamic range to obtain the carbon storage spatial differentiation map and the carbon storage temporal differentiation map.

[0036] A Laplacian edge enhancement operator was used as the spectral edge enhancement filter to perform convolution filtering on the carbon storage spatial differentiation map. The Laplacian convolution kernel is a 3x3 matrix with positive values ​​at the center and negative values ​​at the periphery. After the convolution operation, the edge-enhanced carbon storage spatial differentiation map was obtained.

[0037] Normalized numerical stretching is applied to both the edge-enhanced carbon storage spatial differentiation map and the carbon storage temporal differentiation map. The minimum and maximum values ​​of the edge-enhanced carbon storage spatial differentiation map are calculated, and a linear transformation is performed on each pixel value: C_spatial_norm = (C_spatial_enhanced - Min_spatial) / (Max_spatial - Min_spatial) * D_range, where D_range is the upper limit of the preset display dynamic range. Similarly, the same normalization process is applied to the carbon storage temporal differentiation map. The normalized results are then output as the final carbon storage spatial differentiation map and carbon storage temporal differentiation map.

[0038] Step S130: Delineate land management zones on the carbon storage spatial differentiation map and the carbon storage temporal differentiation map to obtain a set of zoning spatial units consisting of the boundary coordinates of the carbon sink maintenance zone, the boundary coordinates of the carbon source control zone, and the boundary coordinates of the carbon sink enhancement zone.

[0039] Step S131: Extract the carbon storage change gradient in the carbon storage spatial differentiation map, and mark the grids with carbon storage change gradients higher than the preset gradient threshold as spatial abrupt change candidate grids. Extract the closed outer contour of the continuously clustered spatial abrupt change candidate grids as carbon storage spatial break zone line vectors. Extract the change polarity in the carbon storage time differentiation map, and extract the patch boundary formed by the continuous grids with change polarity of carbon storage loss state as carbon loss hotspot area boundary line vectors.

[0040] The gradient magnitude of each pixel is extracted from the carbon storage spatial differentiation map. A preset gradient threshold is set. Pixels with gradient magnitudes greater than the threshold are marked as spatial abrupt change candidate grids. An eight-connected region growing algorithm is used to aggregate all spatially adjacent spatial abrupt change candidate grids into connected regions. The outermost pixel boundary of each connected region is extracted and converted into a closed vector line feature, which serves as the spatial fault line vector of carbon storage.

[0041] The polarity labels of carbon storage time differentiation are extracted. All pixels marked as carbon storage loss states are extracted and aggregated into connected patches using an eight-connected region growing algorithm. The outer boundary of each patch is extracted as the boundary vector of the carbon loss hotspot region.

[0042] Step S132: Perform spatial topological superposition on the spatial fault zone vector of carbon reserves and the boundary line vector of carbon loss hotspot area, extract the vector surface unit jointly divided by the spatial fault zone vector of carbon reserves and the boundary line vector of carbon loss hotspot area as the initial partition surface unit, and assign a unique partition identifier to each initial partition surface unit.

[0043] The spatial vector lines of carbon storage fault zones and the boundary vector lines of carbon loss hotspots are spatially merged to obtain a composite linear vector map layer. This composite linear vector map layer is used as a dividing line and spatially overlaid with the overall areal boundary of the target geographic area. A polygon segmentation algorithm is used to divide the entire region into several vector surface units enclosed by the aforementioned dividing lines. Each vector surface unit is an initial partition surface unit. A unique partition identifier is assigned to each initial partition surface unit.

[0044] Step S133: Extract the numerical mean of the carbon storage spatial differentiation map inside each initial partition surface unit as the spatial heterogeneity attribute value of the initial partition surface unit, and extract the proportion of carbon storage loss state grid area in the carbon storage temporal differentiation map inside each initial partition surface unit as the temporal variation risk value of the initial partition surface unit.

[0045] For each initial partition surface cell, obtain all the raster cells contained within it. Extract the values ​​of these cells from the carbon storage spatial differentiation map, calculate their arithmetic mean, and use it as the spatial heterogeneity attribute value of that initial partition surface cell.

[0046] From the carbon storage temporal diversity map, the number of pixels marked as carbon storage loss states (N_loss) within the initial partitioned surface cell, and the total number of pixels within that cell (N_total) are counted. The area ratio of the carbon storage loss state grid is calculated as R_loss = N_loss / N_total, and this ratio is used as the temporal variation risk value of the initial partitioned surface cell.

[0047] Step S134: Construct a two-dimensional partitioning decision feature vector based on the spatial heterogeneous attribute values ​​and temporal change risk values ​​of the initial partitioning surface unit. Input the two-dimensional partitioning decision feature vector into the gradient boosting-based partitioning decision tree ensemble model. Calculate the category probability score of the initial partitioning surface unit belonging to the carbon sink maintenance, carbon source management, or carbon sink enhancement categories through multiple decision subtrees in the partitioning decision tree ensemble model.

[0048] For each initial partition unit, its spatial heterogeneity attribute value and temporal variation risk value are combined into a two-dimensional feature vector. This feature vector is then input into a pre-trained gradient boosting-based partition decision tree ensemble model. This model consists of M decision subtrees. Each subtree traverses downwards according to the splitting rules of the tree nodes based on the feature vector values, eventually reaching a leaf node. This leaf node outputs a three-dimensional vector, representing the original score predicted by the subtree as belonging to the carbon sink maintenance, carbon source management, and carbon sink enhancement categories, respectively. The output vectors of the M subtrees are summed to obtain the total original score vector. Then, a normalized exponential function is applied to convert the total original score vector into a category probability score vector, where each component represents the probability that the unit belongs to the corresponding category.

[0049] Step S135: For each initial partition surface unit, compare its category probability score for carbon sink maintenance, carbon source control, and carbon sink enhancement. Select the category with the highest category probability score as the initial partition label of the initial partition surface unit. Merge the surface features of initial partition surface units with the same initial partition label and spatial adjacency to generate carbon sink maintenance area, carbon source control area, and carbon sink enhancement area with complete geographic entities.

[0050] For each initial partitioned surface cell, compare the three components in its class probability score vector. Select the class corresponding to the component with the largest value as the initial partition label for that cell. After traversing all initial partitioned surface cells, a set of labeled cells is obtained.

[0051] All initial partition surface units labeled as "carbon sink maintenance" are extracted. A spatial fusion algorithm is then used to merge spatially adjacent units into a single surface region, generating the carbon sink maintenance region. Similarly, adjacent units labeled as "carbon source control" and "carbon sink enhancement" are merged to generate the carbon source control region and carbon sink enhancement region, respectively.

[0052] Step S136: Extract the outer boundary coordinates of the carbon sink maintenance area as the carbon sink maintenance area boundary coordinates, extract the outer boundary coordinates of the carbon source control area as the carbon source control area boundary coordinates, and extract the outer boundary coordinates of the carbon sink enhancement area as the carbon sink enhancement area boundary coordinates. Input the boundary coordinates of the carbon sink maintenance area, the carbon source control area, and the carbon sink enhancement area into the partitioned geometric regularization processor. Fill the tiny holes inside each boundary coordinate using the geometric closing operator and remove the small protrusions outside each boundary coordinate using the geometric opening operator.

[0053] The outer boundaries of the carbon sink maintenance zone, carbon source control zone, and carbon sink enhancement zone are extracted separately to obtain the corresponding boundary coordinate sequences. These three boundary coordinate sequences are then input into a partitioned geometric regularization processor. This processor first performs a geometric closing operation: for the region enclosed by the boundary coordinates, an expansion operation is performed followed by an erosion operation, which fills the tiny pores inside the region. Then, a geometric opening operation is performed: an erosion operation is performed followed by an expansion operation, which removes small protrusions outside the region. After the above processing, the regularized boundary coordinates of the carbon sink maintenance zone, carbon source control zone, and carbon sink enhancement zone are obtained.

[0054] Step S137: Perform coordinate thinning on the boundary coordinates of the carbon sink maintenance area, carbon source control area, and carbon sink enhancement area after geometric processing. Combine the coordinate thinned boundary coordinates of the carbon sink maintenance area, carbon source control area, and carbon sink enhancement area to construct a set of partitioned spatial units. Store the association between the boundary coordinates of the carbon sink maintenance area and the carbon sink maintenance category label, the boundary coordinates of the carbon source control area and the carbon source control category label, and the boundary coordinates of the carbon sink enhancement area and the carbon sink enhancement category label in the set of partitioned spatial units.

[0055] The Douglas-Puk algorithm is used to thin the normalized boundary coordinate sequences. A distance threshold is set, and the algorithm recursively deletes points on the boundary line whose distance from the first and last points is less than this threshold, retaining feature points and thus reducing the number of coordinate points while maintaining the basic shape of the boundary. The thinned boundary coordinates of the carbon sink maintenance zone, carbon source control zone, and carbon sink enhancement zone are combined into a geospatial dataset, i.e., a set of partitioned spatial units. In this dataset, each boundary coordinate sequence is associated with its corresponding category label, forming a key-value pair storage structure.

[0056] Step S140: Based on the analysis of the driving force of carbon storage change in different management zones in the spatial unit set, match the corresponding restoration entries in the preset land low-carbon restoration strategy library for the boundary coordinate coverage areas of the carbon sink maintenance zone, the boundary coordinate coverage areas of the carbon source control zone, and the boundary coordinate coverage areas of the carbon sink enhancement zone, and generate a set of land low-carbon restoration instructions for zone differences.

[0057] Step S141: Based on the boundary coordinates of the carbon sink maintenance zone, the carbon storage time differentiation map is cropped, the spatial distribution pattern of the carbon storage accumulation state grid with changing polarity is extracted within the boundary coordinate coverage area of ​​the carbon sink maintenance zone, and the area continuous expansion direction vector and area continuous expansion rate of the carbon storage accumulation state grid are calculated.

[0058] Using the boundary coordinates of the carbon sink sustaining zone as a mask, the carbon storage time differentiation map was cropped to extract pixels within the carbon sink sustaining zone. Pixels exhibiting the carbon storage accumulation state were selected to form an accumulation state pixel set. Spatial clustering analysis was performed on this set to identify multiple accumulation state patches. For each accumulation state patch, its spatial extent across multiple temporal phases within a continuous historical monitoring period was obtained, and the area of ​​the patch in each phase was calculated. Using the monitoring phase number as the independent variable and the area as the dependent variable, a linear regression was performed to obtain the slope of the regression line, which was taken as the continuous area expansion rate of the patch. Simultaneously, the spatial centroid coordinates of each patch in different temporal phases were calculated. The vector from the centroid of the initial phase to the centroid of the final phase represents the continuous area expansion direction vector. The expansion direction vectors of all accumulation state patches were summed to obtain the overall dominant expansion direction of the carbon sink sustaining zone.

[0059] Step S142: Based on the continuous expansion direction vector and continuous expansion rate of the area, retrieve the set of strategy entries associated with the carbon sink maintenance category tag in the preset land low-carbon remediation strategy library, and select carbon sink maintenance strategy entries whose execution intensity is positively correlated with the continuous expansion rate of the area from the strategy entry set. The carbon sink maintenance strategy entries include maintaining the natural succession process and limiting the scope of human disturbance.

[0060] The pre-defined land low-carbon remediation strategy library is a structured database. The strategy entry set associated with the carbon sink maintenance category tag contains multiple entries. Each entry includes a strategy description, an implementation intensity level, and an applicable expansion rate range. Based on the continuous area expansion rate calculated in step S141, the corresponding implementation intensity level is determined. Strategy entries matching the implementation intensity level are retrieved from the strategy entry set, and entries maintaining natural succession processes and limiting the scope of human disturbance are selected.

[0061] Step S143: Spatial range binding of the boundary coordinates of the carbon sink maintenance area and the carbon sink maintenance strategy entry, generating a carbon sink maintenance area repair instruction package carrying spatial range constraints and strategy content description, wherein the spatial range constraints include the complete coordinate sequence of the boundary coordinates of the carbon sink maintenance area.

[0062] Associate the boundary coordinate sequence of the carbon sink maintenance zone with the selected carbon sink maintenance strategy entries. Create a data structure containing a spatial scope field (storing the complete coordinate sequence of the carbon sink maintenance zone boundary coordinates) and a strategy content description field (storing the specific execution requirements of the above strategy entries), thus serving as the carbon sink maintenance zone repair instruction package.

[0063] Step S144: Based on the boundary coordinates of the carbon source control area, the carbon storage time differentiation map is cropped. The spatial aggregation center coordinates of the carbon storage loss state grid with changing polarity are extracted within the boundary coordinate coverage area of ​​the carbon source control area. The spatial influence radiation radius of the spatial aggregation center coordinates to the surrounding area is calculated. A circular carbon source influence core area is constructed with the spatial aggregation center coordinates as the center and the spatial influence radiation radius as the radius. The area located outside the circular carbon source influence core area within the boundary coordinates of the carbon source control area is designated as the carbon source influence buffer diffusion area, generating a two-layer spatial structure description of the carbon source control area.

[0064] Using the boundary coordinates of the carbon source control area as a mask, the carbon storage time differentiation map is cropped to extract pixels within the carbon source control area. Pixels exhibiting carbon storage loss polarity are selected to form a set of loss-state pixels. A mean-shift clustering algorithm is applied to this set to identify the spatial cluster centers of the loss-state pixels, and the spatial coordinates of each cluster center are calculated. For each cluster center, the average distance to all loss-state pixels within that cluster is calculated as the spatial influence radiation radius. A circular carbon source influence core area is constructed with the cluster center coordinates as the center and the radiation radius as the radius. The remaining area outside all circular carbon source influence core areas within the region enclosed by the boundary coordinates of the carbon source control area is designated as the carbon source influence buffer diffusion area. The boundary ranges of the core area and the buffer diffusion area are recorded, forming a two-layer spatial structure description of the carbon source control area.

[0065] Step S145: Based on the description of the dual-layer spatial structure, retrieve the set of strategy entries associated with the carbon source management category tag in the preset land low-carbon remediation strategy library, and select carbon source blocking strategy entries for the core area of ​​the circular carbon source impact and carbon source conversion strategy entries for the buffer diffusion area of ​​the carbon source impact from the set of strategy entries. The carbon source blocking strategy entries include removing carbon emission activities and cutting off carbon loss pathways, and the carbon source conversion strategy entries include replacing the surface cover and breaking the soil sealing layer.

[0066] In the pre-defined low-carbon land remediation strategy database, strategy entries associated with carbon source management tags are retrieved. Based on the two-layer spatial structure description, for the circular carbon source impact core area, entries with the strategy type of carbon source blocking are selected, including removing carbon emission activities and cutting off carbon loss pathways. For the carbon source impact buffer and diffusion area, entries with the strategy type of carbon source transformation are selected, including surface cover replacement and soil sealing layer removal.

[0067] Step S146: Spatially layer and bind the boundary coordinates of the carbon source control area, the carbon source blocking strategy entry, and the carbon source conversion strategy entry to generate a first sub-instruction carrying the core area range constraint and the carbon source blocking strategy, and a second sub-instruction carrying the buffer diffusion area range constraint and the carbon source conversion strategy. Combine the first sub-instruction and the second sub-instruction to form a carbon source control area repair instruction package.

[0068] Create a composite instruction package. The first sub-instruction contains the boundary coordinates of the core area affected by the circular carbon source and the carbon source blocking strategy entry. The second sub-instruction contains the boundary coordinates of the buffer diffusion area affected by the carbon source and the carbon source conversion strategy entry. Combine the first and second sub-instructions to form the carbon source control area repair instruction package.

[0069] Step S147: Based on the boundary coordinates of the carbon sink enhancement zone, crop the carbon storage spatial differentiation map and the carbon storage temporal differentiation map, identify the coordinates of potential patches to be activated within the coverage area of ​​the carbon sink enhancement zone boundary coordinates where the carbon storage spatial differentiation map value is lower than the surrounding area and the carbon storage temporal differentiation map change polarity is in a stable state, and search for a set of strategy entries associated with the carbon sink enhancement category tag in the preset land low-carbon remediation strategy library based on the coordinates of the potential patches to be activated, select carbon sink enhancement strategy entries that include vegetation type optimization configuration and soil exogenous organic supplementation from the strategy entry set, and adaptively match the vegetation function combination features in the vegetation type optimization configuration with the soil organic carbon spatial distribution map values ​​of the potential patches to be activated.

[0070] Using the boundary coordinates of the carbon sink enhancement zone as a mask, the spatial and temporal carbon storage differentiation maps are cropped. Within the carbon sink enhancement zone, pixels that simultaneously meet two conditions are selected: Condition 1, the spatial carbon storage differentiation value of this pixel is lower than the average value of its surrounding neighboring pixels; Condition 2, the polarity of the change in the temporal carbon storage differentiation map of this pixel represents a stable carbon storage state. The connected patches formed by these pixels are extracted as potential patches to be activated, and their boundary coordinates are recorded.

[0071] The strategy entries associated with the carbon sink enhancement category are retrieved from the pre-defined land low-carbon remediation strategy library, and carbon sink enhancement strategy entries are selected, including vegetation type optimization and soil exogenous organic supplementation. For each potential patch to be activated, the spatial distribution map of soil organic carbon at its location is extracted, and based on the level of this value, appropriate vegetation function combination features are selected from the vegetation type optimization entries.

[0072] Step S148: Spatially bind the boundary coordinates of the carbon sink enhancement zone, the carbon sink enhancement strategy entries, and the coordinates of the potential patches to be activated to generate a carbon sink enhancement zone remediation instruction package carrying the coordinate sequence of potential patches to be activated and the carbon sink enhancement strategy entries. Merge the carbon sink maintenance zone remediation instruction package, the carbon source control zone remediation instruction package, and the carbon sink enhancement zone remediation instruction package to form a low-carbon remediation instruction set for zonal differential land.

[0073] The boundary coordinates of the carbon sink enhancement zone, the selected carbon sink enhancement strategy entries, and the coordinate sequences of potential patches to be activated are linked and bound to generate a carbon sink enhancement zone remediation instruction package. Finally, the carbon sink maintenance zone remediation instruction package, the carbon source control zone remediation instruction package, and the carbon sink enhancement zone remediation instruction package are merged to form a complete set of low-carbon remediation instructions for zoning-differentiated land.

[0074] Step S150: Send the set of instructions for low-carbon restoration of land with zoning differences to the land management terminal to trigger vegetation structure adjustment operation, soil carbon pool protection operation and ecological network connectivity enhancement operation for different zones.

[0075] The low-carbon remediation instruction set for zoning-differentiated land generated in step S148 is sent to the land management terminal via wired or wireless network. This land management terminal can be a geographic information system server deployed in the natural resources management department or a mobile patrol device. Upon receiving the instruction set, the terminal parses the spatial scope constraints and strategy descriptions in each instruction packet to generate specific operational tasks. For carbon sink maintenance zones, operations to limit the scope of human disturbance and monitoring tasks to maintain natural succession processes are triggered. For carbon source control zones, engineering operations to remove carbon emission activities are triggered in the core area, and planting operations to replace surface cover are triggered in the buffer diffusion zone. For carbon sink enhancement zones, vegetation structure adjustment operations to plant according to preferred vegetation types and soil carbon pool protection operations to apply exogenous organic replenishment to the soil are triggered at the coordinate sequence of potential patches to be activated. These operations collectively enhance the connectivity of the ecological network.

[0076] Furthermore, the method may also include: step S210: obtaining a continuous monitoring feedback remote sensing set and a continuous monitoring feedback ground set after the execution of the low-carbon remediation instruction set for zonal differential land on the land management terminal, wherein the continuous monitoring feedback remote sensing set includes feedback multispectral image raster and feedback radar image raster with spatial coordinate markers after the remediation is performed, and the continuous monitoring feedback ground set includes feedback vegetation biomass measurement records and feedback soil organic carbon measurement records with spatial coordinate markers after the remediation is performed.

[0077] Within several monitoring cycles following the execution of the zoning-differentiated land low-carbon remediation instruction set on the land management terminal, multispectral and radar remote sensing images covering the target geographic area are acquired again. Following the same preprocessing procedure as in step S110, feedback multispectral image rasters and feedback radar image rasters are generated to form a continuous monitoring feedback remote sensing set. Simultaneously, within the same monitoring cycle, existing fixed sample plots are redeployed or reused to collect feedback vegetation biomass measurement records and feedback soil organic carbon measurement records, forming a continuous monitoring feedback ground set.

[0078] Step S220: Extract the spatiotemporal differentiation features of the feedback carbon storage from the continuous monitoring feedback remote sensing set and the continuous monitoring feedback ground set to obtain a feedback carbon storage spatial differentiation map reflecting the spatial distribution differences of carbon storage after restoration and a feedback carbon storage temporal differentiation map reflecting the temporal changes of carbon storage after restoration.

[0079] Using the same processing method as step S120, the continuous monitoring feedback remote sensing set and the continuous monitoring feedback ground set are taken as input, and all the sub-steps described in steps S121 to S128 are executed in sequence to generate the repaired feedback carbon storage spatial differentiation map and feedback carbon storage temporal differentiation map, respectively.

[0080] Step S230: Based on the boundary coordinates of the carbon sink maintenance area, crop the spatial differentiation map of the feedback carbon storage and the temporal differentiation map of the feedback carbon storage, and calculate the average net accumulation of feedback carbon storage and the slope change of the accumulation rate of feedback carbon storage within the area covered by the boundary coordinates of the carbon sink maintenance area.

[0081] Using the boundary coordinates of the carbon sink maintenance zone as a mask, the spatial and temporal differentiation maps of the feedback carbon reserves were cropped. The temporal variation of carbon reserves for each pixel within the carbon sink maintenance zone was extracted from the temporal differentiation map, and the arithmetic mean of these variations was calculated to obtain the average net accumulation of feedback carbon reserves. From the spatial differentiation maps of feedback carbon reserves across multiple monitoring time phases after restoration, the total carbon reserves within the carbon sink maintenance zone were extracted. Using the monitoring time phase number as the independent variable and the total carbon reserves as the dependent variable, a linear regression was performed to obtain the slope of the regression line, which was used as the slope change of the feedback carbon reserve accumulation rate.

[0082] Step S240: Based on the boundary coordinates of the carbon source control area, crop the feedback carbon storage spatial differentiation map and the feedback carbon storage temporal differentiation map, and extract the area shrinkage ratio of the carbon storage loss state grid with changing polarity within the boundary coordinate coverage area of ​​the carbon source control area, as well as the fragmentation index change of the carbon storage loss state grid at the core boundary.

[0083] Using the boundary coordinates of the carbon source control area as a mask, the feedback carbon storage temporal differentiation map is cropped. The number of pixels marked as carbon storage loss states within the carbon source control area, as well as the total number of pixels within that area, are counted to calculate the percentage of loss states after restoration. The percentage of loss states before restoration is also obtained, and the area shrinkage ratio is calculated.

[0084] Extract pixels within the buffer zone at the boundary of the core area affected by the circular carbon source. Calculate the number and average area of ​​carbon storage loss-prone patches in these pixels, and calculate the fragmentation index. Obtain the fragmentation index before remediation and calculate its change.

[0085] Step S250: Based on the boundary coordinates of the carbon sink enhancement zone, cut the spatial differentiation map of the feedback carbon storage and the temporal differentiation map of the feedback carbon storage, detect the increase in the value of the feedback carbon storage at the coordinates of the potential patch to be activated compared with that before the repair, and analyze whether its change polarity has changed from a stationary state to a carbon storage accumulation state.

[0086] Using the boundary coordinates of the carbon sink enhancement zone as a mask, the spatial and temporal differentiation maps of the feedback carbon reserves are cropped. The pixel positions corresponding to the coordinate sequences of potential patches to be activated are extracted. The feedback carbon reserve values ​​of these pixels are read from the spatial differentiation map of the feedback carbon reserves, and the pre-remediation carbon reserve values ​​at the corresponding positions are read from the pre-remediation spatial total carbon reserve map to calculate the enhancement rate. The polarity of the change at the coordinates of the potential patches to be activated is read from the temporal differentiation map of the feedback carbon reserves to determine whether this polarity represents a carbon reserve accumulation state.

[0087] Step S260: Construct a carbon sink maintenance zone remediation effectiveness deviation analyzer, compare the difference between the slope change of the feedback carbon storage accumulation rate within the boundary coordinate coverage area of ​​the carbon sink maintenance zone and the preset expected rate range of the carbon sink maintenance strategy item, and generate a carbon sink maintenance effectiveness deviation vector.

[0088] The carbon sink maintenance zone remediation effectiveness deviation analyzer is a mathematical calculation unit. It obtains the slope change value of the feedback carbon storage accumulation rate. It reads the preset expected rate range from the carbon sink maintenance strategy entries. It calculates the deviation value. Using the deviation value as components, it constructs a one-dimensional carbon sink maintenance effectiveness deviation vector.

[0089] Step S270: Construct a carbon source control zone remediation effectiveness deviation analyzer. Compare the area shrinkage ratio and fragmentation index changes of the carbon storage loss state grid within the boundary coordinate coverage area of ​​the carbon source control zone with the expected shrinkage interval preset by the carbon source blocking strategy item and the expected fragmentation reduction interval preset by the carbon source conversion strategy item, respectively, to generate a carbon source control effectiveness deviation vector.

[0090] The carbon source control zone remediation effectiveness deviation analyzer is a mathematical calculation unit. It acquires the changes in area shrinkage ratio and fragmentation index. It reads preset expected shrinkage intervals from the carbon source blocking strategy entries and preset expected fragmentation reduction intervals from the carbon source conversion strategy entries. It calculates the deviation values ​​for the area shrinkage ratio and the fragmentation index change, respectively. These two deviation values ​​are combined into a two-dimensional carbon source control effectiveness deviation vector.

[0091] Step S280: Construct a carbon sink enhancement zone restoration effectiveness deviation analyzer, compare the increase in the feedback carbon storage value at the coordinates of the potential patch to be activated with the expected increase range preset by the carbon sink enhancement strategy item, and determine the consistency between the change polarity at the coordinates of the potential patch to be activated and the expected carbon storage accumulation state, generating a carbon sink potential activation effectiveness deviation vector. Adjust the execution intensity of the carbon sink maintenance strategy item according to the carbon sink maintenance effectiveness deviation vector, adjust the boundary scaling of the spatial range of the carbon source blocking strategy item and the carbon source conversion strategy item according to the carbon source control effectiveness deviation vector, and replace or adjust the type of vegetation function combination characteristics in the carbon sink enhancement strategy item according to the carbon sink potential activation effectiveness deviation vector.

[0092] The carbon sink enhancement zone remediation effectiveness deviation analyzer is a mathematical calculation unit. It obtains the magnitude of the increase in feedback carbon storage values ​​at the coordinates of the potential patch to be activated. It reads the preset expected improvement range from the carbon sink enhancement strategy entries and calculates the deviation value of the increase. Simultaneously, it determines whether the polarity of the change at the coordinates of the potential patch to be activated is consistent with the expected carbon storage accumulation state, generating a state consistency identifier. The deviation value and the state consistency identifier are combined to form the carbon sink potential activation effectiveness deviation vector.

[0093] Based on the deviation values ​​in the carbon sink maintenance effectiveness deviation vector, the execution intensity of the carbon sink maintenance strategy items is numerically adjusted. Based on the deviation values ​​of the area shrinkage ratio and fragmentation index change in the carbon source control effectiveness deviation vector, the spatial range boundaries of the carbon source blocking strategy item and the carbon source conversion strategy item are adjusted to shrink or expand, respectively. Based on the deviation values ​​of the improvement magnitude and state consistency indicators in the carbon sink potential activation effectiveness deviation vector, the vegetation function combination characteristics in the carbon sink enhancement strategy items are either replaced by type (if the state is inconsistent) or their proportion is adjusted (if the improvement magnitude deviates).

[0094] Step S290: Repackage the adjusted carbon sink maintenance strategy entries, adjusted carbon source blocking strategy entries, carbon source conversion strategy entries, and adjusted carbon sink enhancement strategy entries into an iterative optimization land low-carbon remediation instruction set based on zoning differences, and send it to the land management terminal to update the land low-carbon remediation operation in progress.

[0095] The carbon sink maintenance strategy entries, carbon source blocking strategy entries, carbon source conversion strategy entries, and carbon sink enhancement strategy entries, adjusted in step S280, are re-generated into carbon sink maintenance zone restoration instruction packages, carbon source control zone restoration instruction packages, and carbon sink enhancement zone restoration instruction packages according to the encapsulation formats described in steps S140 to S148. These three instruction packages are then merged to form an iteratively optimized set of instructions for low-carbon land restoration based on zoning differences. This set of instructions is then sent again to the land management terminal to update ongoing low-carbon land restoration operations.

[0096] Step S310: Obtain the set of climate environmental factors and the set of human activity factors that affect the dynamic changes of carbon storage in the target geographical area. The set of climate environmental factors includes a grid of annual average temperature distribution, an annual average precipitation distribution, and a grid of total solar radiation distribution with spatial coordinate labels. The set of human activity factors includes a grid of land use type, a grid of population activity thermal data, and a grid of transportation network density with spatial coordinate labels.

[0097] In this embodiment, multi-year average annual temperature, average annual precipitation, and total solar radiation values ​​covering the target geographical area are obtained from meteorological observation station data or reanalysis datasets. Using the Kriging interpolation method, these station observations are spatially interpolated onto a grid identical to the carbon storage spatial differentiation map, generating an average annual temperature distribution grid, an average annual precipitation distribution grid, and a total solar radiation distribution grid. These three grids together constitute the climate environmental factor set.

[0098] Simultaneously, land use type vector data for the target geographic area is obtained from natural resource survey data, rasterized, and each raster cell is assigned a corresponding land use type code to generate a land use type raster. Population activity intensity distribution is obtained from census data or mobile signaling data, and after spatial interpolation and normalization, a population activity heat map raster is generated. Road density (road length per unit area) within each raster cell is calculated from road vector data to generate a traffic network density raster. These three raster arrays together constitute a set of human activity factors.

[0099] Step S320: Input the annual average temperature distribution grid, the annual average precipitation distribution grid, and the total solar radiation distribution grid into the climate factor nonlinear coupler. Through the polynomial kernel function in the climate factor nonlinear coupler, map the annual average temperature distribution grid value, the annual average precipitation distribution grid value, and the total solar radiation distribution grid value to the high-dimensional climate space and calculate the climate stress comprehensive index grid for each spatial location.

[0100] The climate factor nonlinear coupler is a computational module based on the kernel method. For each spatial location (i,j), the annual average temperature T_avg(i,j), annual average precipitation P_avg(i,j), and total solar radiation R_solar(i,j) at that location are extracted to form a three-dimensional input vector X(i,j)=[T_avg,P_avg,R_solar]. A second-order polynomial kernel function K(X_p,X_q)=(X_p·X_q+c)^d is used to map the original three-dimensional vector to a high-dimensional feature space. First, the kernel function value K(i,j)=(X(i,j)·X_ref+c)^d between the input vector and the preset climate stress reference vector X_ref is calculated, where c and d are preset kernel function parameters, and X_ref is a regional climate baseline vector statistically obtained from historical climate data. Then, the kernel function value is converted into a comprehensive climate stress index using the nonlinear transformation function f(X)=Σα_k*K(X,X_k), where α_k is the weight coefficient determined through supervised learning or prior knowledge, and X_k is the climate model vector in the training samples. All raster cells are traversed to generate the comprehensive climate stress index raster C_stress.

[0101] Step S330: Input the land use type raster, the population activity heat raster, and the traffic network density raster into the anthropogenic interference intensity quantifier. Use the weighted summation operator in the anthropogenic interference intensity quantifier to perform weighted superposition of the interference coefficient values ​​of different land types, the normalized values ​​of the population activity heat raster, and the normalized values ​​of the traffic network density raster to generate an anthropogenic interference intensity comprehensive index raster.

[0102] The anthropogenic interference intensity quantizer is a linear weighted calculation module. For each spatial location (i,j), the land use type code L_type(i,j) is read from the land use type raster, and converted into the corresponding land use interference coefficient value D_land(i,j) according to a preset mapping table of different land use interference coefficients. The population activity heat map value H_pop(i,j) is read from the population activity heat map raster, and the road network density value D_road(i,j) is read from the road network density raster. Global normalization is performed on H_pop(i,j) and D_road(i,j) respectively to obtain the normalized population heat map value H_norm(i,j) and the normalized road network density value D_norm(i,j). Then, the comprehensive index of human interference intensity, A_intensity(i,j), is calculated using the weighted summation formula: A_intensity(i,j) = w_land * D_land(i,j) + w_pop * H_norm(i,j) + w_road * D_norm(i,j), where w_land, w_pop, and w_road are preset weight coefficients that satisfy w_land + w_pop + w_road = 1. All raster cells are traversed to generate the comprehensive index raster A_intensity for human interference intensity.

[0103] Step S340: Input the climate stress comprehensive index grid and the anthropogenic disturbance intensity comprehensive index grid into the carbon storage change driving force discrimination network, and calculate the first Granger causality between the climate stress comprehensive index grid and the carbon storage time change grid, and the second Granger causality between the anthropogenic disturbance intensity comprehensive index grid and the carbon storage time change grid through the causal inference layer in the carbon storage change driving force discrimination network.

[0104] The carbon storage change driver discrimination network is a computational network based on Granger causality tests. First, the climate stress composite index raster C_stress and the anthropogenic disturbance intensity composite index raster A_intensity are converted into time series data. For each raster cell, the climate stress composite index values ​​for multiple time phases within the historical continuous monitoring period are extracted to form the time series {C_stress(t)}. Similarly, the anthropogenic disturbance intensity composite index time series {A_intensity(t)} is extracted. From the carbon storage time change raster C_time_change, the carbon storage change series {C_change(t)} for the corresponding location and time phase is extracted.

[0105] For each grid cell, two vector autoregressive models are constructed. Model 1: Regression is performed with C_change(t) as the dependent variable and its own lagged p-period values ​​{C_change(t-1),...,C_change(tp)} as independent variables, calculating the residual sum of squares RSS_restricted. Model 2: Based on the independent variables of Model 1, the lagged p-period values ​​{C_stress(t-1),...,C_stress(tp)} of the climate stress composite index series are added, and regression is performed again, calculating the residual sum of squares RSS_unrestricted. The F-statistic F_clim = ((RSS_restricted-RSS_unrestricted) / p) / (RSS_unrestricted / (T-2p-1)), where T is the time series length. This F-statistic is the first Granger causality of the climate stress composite index on carbon storage change. Similarly, a model is constructed with the anthropogenic disturbance intensity composite index series as an additional independent variable, calculating the second Granger causality F_human.

[0106] Step S350: For the grids in the carbon storage time differentiation map where the change polarity is the carbon storage accumulation state, the carbon storage accumulation state driver is classified as either climate-dominated or anthropogenic suppression-weakening driver based on the relative magnitude of the first Granger causality and the second Granger causality, and a carbon sink formation driving force label distribution map is generated.

[0107] Extract all raster cells in the carbon storage time-separation map whose change polarity corresponds to the carbon storage accumulation state. For each cell, obtain its corresponding first Granger causality F_clim and second Granger causality F_human. Set a relative dominance threshold R_dom. If F_clim / F_human > R_dom, the carbon storage accumulation state driver of that cell is classified as climate-dominant; if F_human / F_clim > R_dom, it is classified as anthropogenic suppression-weakening driver; if the ratio is within the threshold range, it is classified as a combined climate-anthropogenic driver. Assign a label code to each driver type to generate a carbon sink formation driver label distribution map.

[0108] Step S360: For the grids in the carbon storage time differentiation map where the change polarity is carbon storage loss state, the carbon storage loss state driver is classified as either climate stress-driven or anthropogenic activity-enhanced driver based on the relative magnitude of the first Granger causality and the second Granger causality, and a carbon source formation driving force label distribution map is generated.

[0109] Extract all raster cells in the carbon storage temporal differentiation map whose change polarity represents carbon storage loss states. For each of these cells, obtain its corresponding first Granger causality F_clim and second Granger causality F_human. Using the same comparison logic as in step S350, if F_clim / F_human > R_dom, then the carbon storage loss state driver of this cell is classified as climate stress-dominated; if F_human / F_clim > R_dom, then it is classified as anthropogenic enhancement-driven; otherwise, it is classified as a combined climate-anthropogenic driver. Assign a label code to each driver type to generate a carbon source formation driving force label distribution map.

[0110] Step S370: Perform spatial overlay analysis on the carbon sink formation driving force label distribution map and the partitioned spatial unit set, and statistically analyze the area ratio of climate-dominant driving grids and anthropogenic suppression weakening driving grids within the boundary coordinate coverage area of ​​the carbon sink maintenance zone, and determine the larger ratio as the main driving force for carbon storage change in the carbon sink maintenance zone; Perform spatial overlay analysis on the carbon source formation driving force label distribution map and the partitioned spatial unit set, and statistically analyze the area ratio of climate stress-dominant driving grids and anthropogenic activity enhancement driving grids within the boundary coordinate coverage area of ​​the carbon source control zone, and determine the larger ratio as the main driving force for carbon storage change in the carbon source control zone; Perform spatial overlay analysis on the carbon sink formation driving force label distribution map and the carbon source formation driving force label distribution map and the partitioned spatial unit set, and statistically analyze the area distribution characteristics of climate-dominant driving grids, anthropogenic suppression weakening driving grids, climate stress-dominant driving grids, and anthropogenic activity enhancement driving grids within the boundary coordinate coverage area of ​​the carbon sink enhancement zone, and generate the main driving force for carbon storage change in the carbon sink enhancement zone.

[0111] Using the boundary coordinates of the carbon sink maintenance zone as a mask, a distribution map of carbon sink driving force labels was cropped. The number of pixels (N_clim_dom) of climate-dominant driving force labels and the number of pixels (N_human_dom) of anthropogenic suppression-weakening driving force labels within the cropped area were statistically analyzed. The area proportions of these two types were calculated as R_clim_maintain = N_clim_dom / (N_clim_dom + N_human_dom) and R_human_maintain = N_human_dom / (N_clim_dom + N_human_dom). The driving force corresponding to the larger proportion was identified as the primary driving force for carbon storage changes in the carbon sink maintenance zone.

[0112] Similarly, by using the boundary coordinates of the carbon source control area to crop the distribution map of carbon source driving force labels, and by counting the number of pixels of climate stress-dominant driving force labels and human activity-enhanced driving force labels, the main driving force of carbon storage changes in the carbon source control area can be determined.

[0113] Using the boundary coordinates of the carbon sink enhancement zone, both the carbon sink formation driving force label distribution map and the carbon source formation driving force label distribution map are cropped. The area proportion of the four driving labels (climate-dominated driving, anthropogenic suppression and weakening driving, climate stress-dominated driving, and anthropogenic activity enhancement driving) is calculated, and the driving type with the highest proportion is taken as the main driving force for the change of carbon storage in the carbon sink enhancement zone.

[0114] Step S380: The main driving force of carbon storage change in the carbon sink maintenance zone, the main driving force of carbon storage change in the carbon source control zone, and the main driving force of carbon storage change in the carbon sink enhancement zone are added to the zonal spatial unit set as input parameters for matching strategy entries in the zonal differential land low-carbon restoration instruction set.

[0115] The primary driving force information of carbon storage changes for each of the three zones determined in step S370 is appended as an attribute field to the corresponding boundary coordinate entries in the spatial unit set of the zones. In the subsequent execution of step S140 for strategy entry matching, the primary driving force information is used as an additional input parameter to filter remediation strategy entries from the preset land low-carbon remediation strategy library that match the driving force type.

[0116] Step S410: Obtain the set of field patrol survey tracks for the target geographic area during the historical continuous monitoring period. The set of field patrol survey tracks includes vegetation community structure records and soil profile morphology records collected along the patrol path with spatial coordinate markers and collection time markers.

[0117] In this embodiment, survey data collected by field patrol personnel along fixed or random patrol routes within a target geographical area during a continuous historical monitoring period are acquired. Each survey record includes spatial coordinate markers, collection time markers, and records of vegetation community structure (textual description, such as "The dominant species in the tree layer is cork oak, the shrub coverage is 65%, and the herbaceous diversity is high") and soil profile morphology (textual description, such as "The humus layer is approximately eight centimeters thick, the soil color is dark brown, and the soil compaction is loose"). All of the above records constitute a set of field patrol survey tracks.

[0118] Step S420: Extract key information from the vegetation community structure record. Extract the dominant tree species names, shrub coverage descriptions, and herb diversity descriptions from the vegetation community structure record using the named entity recognition operator of the natural language processing component. Convert the dominant tree species names, shrub coverage descriptions, and herb diversity descriptions into structured vegetation vertical structure vectors. Also, extract key information from the soil profile morphology record. Extract the humus layer thickness description, soil color description, and soil compaction description from the soil profile morphology record using the named entity recognition operator of the natural language processing component. Convert the humus layer thickness description, soil color description, and soil compaction description into structured soil profile trait vectors.

[0119] A pre-trained natural language processing model (e.g., a transformer-based bidirectional encoder representation model) is used as the named entity recognition operator. The text recording the vegetation community structure is input into this operator, and the model outputs a sequence of entity labels from the text through a conditional random field layer. Entity values ​​labeled as "dominant tree species" (e.g., "Quercus variabilis"), "shrub cover" (e.g., "65%)", and "herbaceous diversity" (e.g., "relatively high") are identified. The names of dominant tree species are converted into numerical codes using a pre-defined native plant coding table; the percentage values ​​in the shrub cover descriptions are extracted; and the herbaceous diversity descriptions (e.g., "relatively low", "moderate", "relatively high") are mapped to rank values ​​(e.g., 1, 2, 3). These three values ​​are combined into a three-dimensional structured vegetation vertical structure vector.

[0120] Similarly, the soil profile morphology record text is input into the named entity recognition operator to extract the thickness value from the humus layer thickness description, the soil color description (converted to numerical code through a preset soil colorimetric card encoding table), and the soil compaction description (e.g., "loose," "compact," and "hardened" mapped to grade values). These three values ​​are then combined into a three-dimensional structured soil profile morphology vector.

[0121] Step S430: Spatially associate the structured vegetation vertical structure vector and the structured soil profile trait vector with the grid in the carbon storage spatial differentiation map according to the spatial coordinate label, to form a field-enhanced carbon storage grid sequence containing the vegetation vertical structure vector and the soil profile trait vector.

[0122] For each field patrol survey record, its spatial coordinates (X_i, Y_i) are obtained. These coordinates are mapped onto the raster grid of the carbon storage spatial differentiation map to determine the cell location. The original carbon storage value, vegetation vertical structure vector, and soil profile trait vector at that cell location are associated and stored to form a field-enhanced carbon storage grid record. The records corresponding to all survey records are summarized to form a field-enhanced carbon storage grid sequence.

[0123] Step S440: Construct a correlation analysis network between vegetation structure complexity and carbon storage stability. Input the structured vegetation vertical structure vector from the field enhanced carbon storage grid sequence into the structure complexity layer of the correlation analysis network. Calculate the species evenness of dominant tree species, the cover continuity described by shrub cover, and the species richness described by herb diversity. Normalize each indicator and weightedly fuse them to generate a comprehensive score map of vegetation structure complexity. Construct a correlation analysis network between soil properties and carbon pool storage persistence. Input the structured soil profile trait vector from the field enhanced carbon storage grid sequence into the carbon pool stability layer of the correlation analysis network. Calculate the accumulated organic matter thickness based on the humus layer thickness description, the degree of organic matter darkening based on the soil color description, and the soil bulk density estimation based on the soil compaction description. Normalize each indicator and weightedly fuse them to generate a comprehensive score map of soil carbon pool persistence.

[0124] The correlation analysis network between vegetation structure complexity and carbon storage stability is a computational network. Its structure complexity layer receives a structured vegetation vertical structure vector. For dominant tree species, the tree species and their frequencies at all sampling points within the region are statistically analyzed, and the Shannon-Wiener diversity index is calculated as the species evenness index E_tree. The cover value in the shrub cover description is directly used as the cover continuity index C_shrub. The rank value in the herbaceous diversity description is used as the species richness index R_herb. E_tree, C_shrub, and R_herb are respectively min-max normalized to obtain normalized values ​​E_norm, C_norm, and R_norm. Then, a weighted fusion formula is used to calculate the comprehensive vegetation structure complexity score V_complex = w_tree * E_norm + w_shrub * C_norm + w_herb * R_norm, where w_tree, w_shrub, and w_herb are preset weight coefficients. The scores of all sampling points are extended to the entire region through spatial interpolation to generate a comprehensive vegetation structure complexity score map.

[0125] The correlation analysis network between soil properties and carbon pool persistence is a computational network. Its carbon pool stability layer receives structured soil profile trait vectors. The thickness value in the humus layer thickness description is directly used as the organic matter accumulation thickness H_humus. The soil color description is converted to a darkening index D_color using a pre-defined organic matter darkening degree mapping table. The soil compaction description is calculated using a pre-defined compaction-bulk weight conversion model to estimate the soil bulk density value B_bulk. H_humus, D_color, and B_bulk are normalized separately, and then a weighted fusion formula is used to calculate the comprehensive soil carbon pool persistence score S_persistence = w_humus*H_norm + w_color*D_norm + w_bulk*B_norm. The scores of all samples are spatially interpolated to extend to the entire region, generating a comprehensive soil carbon pool persistence score map.

[0126] Step S450: Perform spatial raster correlation analysis on the vegetation structure complexity comprehensive score map and the soil carbon pool persistence comprehensive score map with the carbon storage time differentiation map, respectively. Calculate the first spatial correlation coefficient matrix of the carbon storage accumulation rate in the vegetation structure complexity comprehensive score map and the carbon storage time differentiation map, and the second spatial correlation coefficient matrix of the carbon storage accumulation rate in the soil carbon pool persistence comprehensive score map and the carbon storage time differentiation map. Generate a first spatial heat map of the influence of vegetation complexity on carbon storage time differentiation based on the first spatial correlation coefficient matrix, and generate a second spatial heat map of the influence of soil persistence on carbon storage time differentiation based on the second spatial correlation coefficient matrix.

[0127] For each raster cell, extract its vegetation structure complexity score V_complex and carbon storage accumulation rate R_acc (derived from the carbon storage temporal variation) from the carbon storage temporal variation map. Using the Pearson correlation coefficient formula, with V_complex as the independent variable and R_acc as the dependent variable, calculate the correlation coefficient r_local within a local window (e.g., a 3x3 or 5x5 neighborhood). Traverse all raster cells to generate a first spatial correlation coefficient matrix. Map each correlation coefficient value to its corresponding raster cell to generate a first spatial heatmap, where color intensity indicates the strength of the positive or negative correlation between vegetation complexity and carbon storage accumulation rate.

[0128] Similarly, using the same method, the Pearson correlation coefficient between the soil carbon pool persistence comprehensive score S_persistence and the carbon storage accumulation rate R_acc within the local window is calculated, generating a second spatial correlation coefficient matrix and a second spatial heatmap.

[0129] Step S460: Spatially crop the first spatial heat map and the boundary coordinates of the carbon sink enhancement zone to identify the coordinate sequence of the first priority remediation sub-region within the coverage area of ​​the carbon sink enhancement zone boundary coordinates, where the comprehensive score of vegetation structure complexity is positively correlated with the carbon storage accumulation rate; and spatially crop the second spatial heat map and the boundary coordinates of the carbon sink enhancement zone to identify the coordinate sequence of the second priority remediation sub-region within the coverage area of ​​the carbon sink enhancement zone boundary coordinates, where the comprehensive score of soil carbon pool persistence is positively correlated with the carbon storage accumulation rate, and supplement the first and second priority remediation sub-region coordinate sequences to the partitioned spatial unit set to refine the spatial placement of carbon sink enhancement zone strategy items.

[0130] The boundary coordinates of the carbon sink enhancement zone are used as a mask to crop the first spatial heat map. Grid cells with a correlation coefficient greater than a preset positive correlation threshold (e.g., 0.5) are extracted from the cropped heat map. The boundary coordinates of the connected patches formed by these cells are extracted and used as the coordinate sequence of the first priority restoration sub-region. These regions represent areas where increased vegetation structure complexity significantly promotes carbon storage accumulation.

[0131] Similarly, using the boundary coordinates of the carbon sink enhancement zone as a mask, the second spatial heat map is cropped, and grid cells with correlation coefficients greater than the preset positive correlation threshold are extracted to generate the coordinate sequence of the second priority remediation sub-region. The above-mentioned regions represent areas where the improvement of soil carbon pool persistence has a significant promoting effect on carbon storage accumulation.

[0132] The coordinate sequences of the first and second priority remediation sub-regions will be added as new attribute fields to the records associated with the carbon sink enhancement zone in the zonal spatial unit set. When generating subsequent carbon sink enhancement zone remediation instruction packages, these two priority sub-region coordinate sequences will be prioritized as key location areas, and vegetation type optimization and exogenous organic soil replenishment operations in the carbon sink enhancement strategy entries will be configured with priority.

[0133] Step S510: Obtain the current land use planning atlas and ecological protection red line boundary set corresponding to different management zones within the target geographical area. The current land use planning atlas contains legally binding planning land category patch boundaries and planning land category use codes. The ecological protection red line boundary set contains legally binding ecological red line boundaries and ecological red line control levels.

[0134] In this embodiment, vector data of the current land use plan approved and implemented by the natural resources authority is obtained. This data contains multiple planned land use patches, each with a boundary coordinate sequence and a planned land use code (e.g., code "01" represents permanent basic farmland, "02" represents general cultivated land, and "03" represents a permitted construction area). The above data constitutes the current land use plan atlas.

[0135] Simultaneously, vector data of ecological protection red lines delineated by the ecological and environmental authorities are obtained, which includes the coordinate sequence of ecological red line boundaries and the control level corresponding to each red line zone (e.g., "Level 1 Control Zone" and "Level 2 Control Zone"). The above data constitutes the ecological protection red line boundary set.

[0136] Step S520: Perform spatial conflict detection on the boundary coordinates of the planned land use patch in the current land use planning atlas and the boundary coordinates of the carbon sink maintenance area in the zoning spatial unit set, and extract the first set of spatial conflict coordinates that overlap with the boundary coordinates of the carbon sink maintenance area in the boundary of the planned land use patch and whose planned land use code is a permitted development and construction class.

[0137] Spatial overlay analysis is performed on the boundaries of planned land use patches in the current land use planning atlas and the boundaries of carbon sink conservation areas. A polygon intersection algorithm is used to detect whether there are overlapping areas between the planned land use patches and the carbon sink conservation areas. For overlapping patches, it is checked whether their planned land use codes belong to the permitted development and construction category (e.g., code "03" for permitted construction or "04" for mining land). The boundary coordinates of the overlapping areas of patches that meet the conditions of overlap and permitted development and construction are extracted to form the first spatial conflict coordinate set.

[0138] Step S530: Perform spatial conflict detection on the boundary coordinates of the planned land use patch in the current land use planning atlas and the boundary coordinates of the carbon source control area in the zoning spatial unit set, and extract the second set of spatial conflict coordinates in the boundary of the planned land use patch that overlaps with the boundary coordinates of the carbon source control area and whose planned land use code is the carbon emission enhancement class.

[0139] Spatial overlay analysis is performed on the boundaries of planned land use patches in the current land use planning atlas and the boundaries of carbon source control zones. For overlapping patches, it is checked whether their planned land use codes belong to the category allowing carbon emission enhancement (e.g., code "05" for industrial land or "06" for warehousing land). The boundary coordinates of the overlapping areas of patches that meet the criteria of overlapping and belonging to the category allowing carbon emission enhancement are extracted to form a second spatial conflict coordinate set.

[0140] Step S540: Perform spatial conflict detection on the boundary coordinates of the planned land use patch in the current land use planning atlas and the boundary coordinates of the carbon sink enhancement zone in the zoning spatial unit set, and extract the third spatial conflict coordinate set of the planned land use patch boundary that overlaps with the boundary coordinates of the carbon sink enhancement zone and whose planned land use code is prohibited for ecological restoration.

[0141] Spatial overlay analysis is performed on the boundaries of planned land use patches in the current land use planning atlas and the boundaries of carbon sink enhancement zones. For overlapping patches, it is checked whether their planned land use codes belong to the prohibited ecological restoration category (e.g., code "07" for urban and rural road land or "08" for hydraulic engineering construction land). The boundary coordinates of the overlapping areas of patches that meet the conditions of overlapping and being prohibited ecological restoration categories are extracted to form a third spatial conflict coordinate set.

[0142] Step S550: Perform spatial overlay analysis on the ecological red line boundary coordinates of the ecological protection red line boundary set and the carbon sink maintenance area boundary coordinates of the partitioned spatial unit set, extract the coordinates of the first ecological void area that is not covered by the ecological red line boundary within the carbon sink maintenance area boundary coordinates, and calculate the area size and spatial dispersion of the first ecological void area coordinates.

[0143] Spatial overlay analysis is performed between the ecological protection redline boundary set and the carbon sink maintenance zone boundary coordinates. The difference between the area enclosed by the carbon sink maintenance zone boundary coordinates and the ecological redline boundary is calculated, i.e., the area within the carbon sink maintenance zone but not covered by the ecological redline. The boundary coordinates of this area are extracted as the coordinates of the first ecological gap zone. The area value A_gap of this area, the number of patches N_patches within the area, and the ratio of the average patch perimeter to the area are calculated to obtain the spatial dispersion index D_dispersion.

[0144] Step S560: Perform spatial overlay analysis on the ecological red line boundary coordinates of the ecological protection red line boundary set and the carbon source control area boundary coordinates of the partitioned spatial unit set, extract the coordinates of the overlapping area between the carbon source control area boundary coordinates and the ecological red line boundary, and calculate the compatibility index between the carbon source blocking strategy entries and the ecological red line control level within the overlapping area coordinates.

[0145] Spatial overlay analysis is performed on the boundary coordinates of the ecological protection red line and the carbon source control zone boundary to extract the intersection area, which is then used as the coordinates of the overlapping area. For this overlapping area, its ecological red line control level (Level 1 or Level 2) is determined. From the carbon source blocking strategy items selected in step S145, their environmental impact intensity levels are extracted. Using a preset compatibility matrix table, the compatibility degree between the ecological red line control level and the environmental impact intensity level of the carbon source blocking strategy item is queried, and the compatibility index C_compat (ranging from 0 to 1, with higher values ​​indicating greater compatibility) is output.

[0146] Step S570: Perform spatial overlay analysis on the ecological red line boundary coordinates of the ecological protection red line boundary set and the carbon sink enhancement zone boundary coordinates of the partitioned spatial unit set, extract the coordinates of the ecological network connected to the ecological red line boundary outside the carbon sink enhancement zone boundary coordinates, and analyze the minimum spatial distance distribution between the coordinates of the connected to the enhanced zone and the carbon sink enhancement zone boundary coordinates.

[0147] Spatial adjacency analysis is performed between the ecological protection redline boundary set and the carbon sink enhancement zone boundary coordinates. Regions outside the carbon sink enhancement zone boundary coordinates, whose distance from the carbon sink enhancement zone boundary is less than a preset adjacency threshold (e.g., 100 meters), and adjacent to the ecological redline boundary are extracted as the ecological network connectivity enhancement area coordinates. The shortest Euclidean distance from each connected unit in this region to the carbon sink enhancement zone boundary is calculated, generating minimum spatial distance distribution statistics (including minimum, maximum, average, and standard deviation).

[0148] Step S580: Generate a land use code adjustment suggestion for the current land use planning atlas based on the first spatial conflict coordinate set, the second spatial conflict coordinate set, and the third spatial conflict coordinate set. The land use code adjustment suggestion includes changing the permitted development and construction category to the prohibited development and construction category, changing the permitted carbon emission enhancement category to the restricted carbon emission enhancement category, and changing the prohibited ecological restoration category to the permitted ecological restoration category.

[0149] The first, second, and third spatial conflict coordinate sets are aggregated to form a list of conflict areas. For conflict areas located within carbon sink maintenance zones (permitted development and construction category), an adjustment suggestion is generated: change the planned land use code from permitted development and construction to prohibited development and construction category. For conflict areas located within carbon source control zones (permitted carbon emission enhancement category), an adjustment suggestion is generated: change to restricted carbon emission enhancement category. For conflict areas located within carbon sink enhancement zones (prohibited ecological restoration category), an adjustment suggestion is generated: change to permitted ecological restoration category.

[0150] Step S590: Generate a boundary expansion proposal for the ecological protection red line boundary set based on the area size and spatial dispersion of the first ecological void area coordinates. The boundary expansion proposal includes the spatial expansion direction and area size of incorporating the first ecological void area coordinates into the ecological red line boundary.

[0151] Based on the spatial distribution of the coordinates of the first ecological gap area, the direction of the shortest distance between it and the existing ecological red line boundary is calculated, which serves as the direction for spatial expansion. Based on the area size A_gap of the first ecological gap area coordinates, the suggested area size for incorporating it into the ecological red line boundary is determined. A boundary expansion proposal is generated: the area enclosed by the coordinates of the first ecological gap area is included within the ecological red line range, with the expansion direction being towards this area, and the expansion area being approximately A_gap.

[0152] Step S5100: Based on the compatibility index calculation results of the overlapping area coordinates and the minimum spatial distance distribution between the coordinates of the connected area to be strengthened and the boundary coordinates of the carbon sink enhancement area, generate management level optimization suggestions for the ecological protection red line boundary set, and add the land use code adjustment suggestions, boundary expansion suggestions and management level optimization suggestions as auxiliary information to the low-carbon restoration instruction set for zoning differences.

[0153] For the coordinates of the overlapping area, if the compatibility index C_compat is lower than the preset compatibility threshold, an optimization suggestion for the control level is generated: the ecological red line control level of the area is upgraded by one level (for example, from level two control to level one control) to enhance the binding force on the implementation of the carbon source blocking strategy.

[0154] For the coordinates of the connected areas to be strengthened, based on the distribution of their minimum spatial distance from the boundary of the carbon sink enhancement area, the following optimization suggestions for the control level are generated: the coordinates of the connected areas to be strengthened that are close to the boundary of the carbon sink enhancement area (e.g., less than the average distance) are included in the management of the ecological red line buffer zone, and control measures coordinated with the carbon sink enhancement area are implemented.

[0155] The usage code adjustment suggestions generated in step S580, the boundary expansion suggestions generated in step S590, and the control level optimization suggestions generated in this step are added as auxiliary information fields to the low-carbon remediation instruction set for zoning differences, for reference by the land management terminal when making decisions.

[0156] Step S610: Obtain the regional energy consumption structure set and the regional transportation mode set within and around the target geographical area. The regional energy consumption structure set includes the proportion of fossil energy consumption and the proportion of clean energy consumption marked by administrative divisions. The regional transportation mode set includes the motor vehicle travel modal share and the non-motor vehicle travel modal share marked by administrative divisions.

[0157] In this embodiment, the proportions of fossil fuel consumption and clean energy consumption of each administrative unit (e.g., county-level administrative region) within and around the target geographical area are obtained from statistical yearbooks or energy balance sheets to form a regional energy consumption structure set. The motor vehicle travel modal share (the proportion of motor vehicle trips to total trips) and non-motor vehicle travel modal share of each administrative unit are obtained from traffic survey data or mobile phone signaling data to form a regional traffic travel pattern set.

[0158] Step S620: Allocate the proportion of fossil energy consumption and the share of motor vehicle travel to the grid scale according to the spatial weight of nighttime light intensity to generate a gridded fossil carbon emission source intensity map and a gridded traffic carbon emission source intensity map.

[0159] Acquire nighttime light remote sensing imagery of the target geographic area and obtain the nighttime light intensity value NTL(i,j) after radiometric calibration. For each administrative unit, calculate the sum of nighttime light intensity values ​​NTL_sum_admin for all grids within that unit. For each grid within that administrative unit, calculate its weighting coefficient W_ntl(i,j) = NTL(i,j) / NTL_sum_admin. Assign the fossil fuel consumption ratio P_fossil_admin of that administrative unit to each grid according to the weighting coefficient, obtaining the gridded fossil carbon emission source intensity value E_fossil(i,j) = P_fossil_admin * W_ntl(i,j). Similarly, assign the motor vehicle travel share R_vehicle_admin to each grid according to the same weighting coefficient, obtaining the gridded traffic carbon emission source intensity value E_traffic(i,j) = R_vehicle_admin * W_ntl(i,j). Generate the gridded fossil carbon emission source intensity map and the gridded traffic carbon emission source intensity map respectively.

[0160] Step S630: Spatially overlay the gridded fossil carbon emission source intensity map and the gridded transportation carbon emission source intensity map to generate a regional carbon emission source comprehensive intensity map, and extract the grids in the regional carbon emission source comprehensive intensity map whose comprehensive carbon emission source intensity exceeds a preset intensity threshold as a high carbon emission source grid set.

[0161] For each grid cell, calculate the comprehensive carbon emission source intensity E_total(i,j) = E_fossil(i,j) + E_traffic(i,j). Generate a comprehensive carbon emission source intensity map for the region. Set a preset intensity threshold T_emission. Extract grid cells where E_total(i,j) > T_emission to form a high carbon emission source grid set.

[0162] Step S640: Calculate the spatial Euclidean distance matrix between the spatial centroid coordinates of the high carbon emission source grid set and the boundary coordinates of the carbon source control area. Mark the boundary coordinates of the carbon source control area that are less than a preset distance threshold in the spatial Euclidean distance matrix as the boundary coordinates of the near-source control area, and mark the boundary coordinates of the carbon source control area that are greater than the preset distance threshold as the boundary coordinates of the far-source control area.

[0163] Calculate the spatial centroid coordinates of the high-carbon emission source grid set, that is, the arithmetic mean (X_center, Y_center) of all high-carbon emission source grid coordinates. For each boundary point in the carbon source control area boundary coordinate sequence or each control area polygon, calculate the Euclidean distance D_euclidean between it and the spatial centroid coordinates. Set a preset distance threshold T_distance. Mark the control areas with D_euclidean < T_distance as the boundary coordinates of the near-source control area, and mark the control areas with D_euclidean ≥ T_distance as the boundary coordinates of the far-source control area.

[0164] Step S650: Extract the industrial supporting strategy entries for energy structure optimization from the preset land low-carbon restoration strategy library for the area covered by the boundary coordinates of the near-source control area. The industrial supporting strategy entries include regional energy replacement to increase the proportion of clean energy consumption and public transportation induction to reduce the share rate of motor vehicle trips.

[0165] In the preset land low-carbon restoration strategy library, retrieve the set of strategy entries associated with the "near-source control" label. Screen out the industrial supporting strategy entries for energy structure optimization, including: regional energy replacement strategy (replacing fossil energy consumption with clean energy, such as "coal to gas", "coal to electricity") and public transportation induction strategy (optimizing bus routes, increasing bus frequencies, and building non-motor vehicle lanes to reduce the share rate of motor vehicle trips).

[0166] Step S660: Extract the ecological barrier strategy entries for blocking carbon transportation pathways from the preset land low-carbon restoration strategy library for the area covered by the boundary coordinates of the far-source control area. The ecological barrier strategy entries include the windward interception construction of configuring a vegetation belt with high carbon absorption efficiency in the windward direction of the boundary coordinates of the far-source control area and the internal carbon sink enhancement of constructing a multi-layer vegetation three-dimensional configuration inside the boundary coordinates of the far-source control area.

[0167] In the preset land low-carbon restoration strategy library, retrieve the set of strategy entries associated with the "far-source control" label. Screen out the ecological barrier strategy entries for blocking carbon transportation pathways, including: windward interception construction strategy (configuring a tree-shrub-herb complex vegetation belt with high carbon absorption efficiency on the windward main direction side of the boundary coordinates of the far-source control area) and internal carbon sink enhancement strategy (implementing a multi-layer vegetation three-dimensional configuration inside the far-source control area to form a complete carbon absorption structure from the ground cover layer to the canopy layer).

[0168] Step S670: Spatially bind the boundary coordinates of the near-source control area with the industrial supporting strategy entries to generate a regional collaborative near-source restoration sub-instruction; spatially bind the boundary coordinates of the far-source control area with the ecological barrier strategy entries to generate an ecological barrier far-source restoration sub-instruction; and merge the regional collaborative near-source restoration sub-instruction and the ecological barrier far-source restoration sub-instruction into the zonal differential land low-carbon restoration instruction set.

[0169] The coordinate sequence of the near-source pollution control zone boundary is correlated with the industrial support strategy items (regional energy replacement, public transportation guidance) to generate a regional collaborative near-source pollution remediation sub-instruction. The coordinate sequence of the far-source pollution control zone boundary is correlated with the ecological barrier strategy items (windward interception construction, internal carbon sink enhancement) to generate a far-source pollution remediation sub-instruction for ecological barriers. These two sub-instructions are added to the zonal differential land low-carbon remediation instruction set generated in step S148.

[0170] Step S680: Obtain the hydrological connectivity set within the watershed of the target geographical area. The hydrological connectivity set includes a river system linear vector with river hierarchical coding and a sub-watershed area vector with watershed boundary coding. Calculate the riparian vegetation cover integrity index of the river system linear vector and the soil erosion sensitivity index of the sub-watershed area vector.

[0171] Vector data of river systems within the target geographic area's watershed are obtained from water resources departments or hydrological databases. Each river segment has a Strahler classification code. Simultaneously, sub-watershed boundary vector data are acquired, with each sub-watershed surface having a unique watershed boundary code. For linear river system vectors, a fixed-width riparian buffer zone is generated with each river segment as its centerline. Vegetation cover within the buffer zone is calculated using high-resolution remote sensing imagery, and the riparian vegetation cover integrity index C_riparian = (actual vegetation cover area / total buffer zone area) is calculated. For sub-watershed surface vectors, the average slope, soil erodibility factor, and runoff length for each sub-watershed are calculated using a digital elevation model. The soil erosion sensitivity index S_erosion is then calculated using a general soil loss equation.

[0172] Step S690: Perform spatial overlay analysis on the boundary coordinates of the carbon sink maintenance area with the linear vector of the river system and the planar vector of the sub-basin, extract the vector of the riparian vegetation cover integrity index below the preset integrity threshold inside the boundary coordinates of the carbon sink maintenance area, and extract the vector of the slope soil and water conservation area with the soil and water loss sensitivity index above the preset sensitivity threshold inside the boundary coordinates of the carbon sink maintenance area.

[0173] Using the boundary coordinates of the carbon sink maintenance zone as a mask, spatial overlay is performed with the linear vector of the river system. River segments within the carbon sink maintenance zone are selected, and their riparian vegetation cover integrity index (C_riparian) is checked to see if it is lower than the preset integrity threshold (T_riparian). River segments that meet the criteria are extracted and used as vectors for riparian degradation restoration.

[0174] Simultaneously, the boundary coordinates of the carbon sink conservation zone are used as a mask and spatially superimposed with the sub-watershed surface vectors. Sub-watershed surfaces within the carbon sink conservation zone are selected, and their soil erosion sensitivity index S_erosion is checked to see if it exceeds the preset sensitivity threshold T_erosion. Sub-watershed surfaces that meet the conditions are extracted as the slope soil and water conservation zone vectors.

[0175] Step S6100: Extract river buffer zone vegetation reconstruction strategy entries from the preset land low-carbon restoration strategy library for the vector of the waterfront degradation restoration section, and extract slope runoff regulation strategy entries from the preset land low-carbon restoration strategy library for the vector of the slope soil and water conservation area. Bind the river buffer zone vegetation reconstruction strategy entries with the vector space of the waterfront degradation restoration section to generate watershed collaborative waterfront restoration sub-instructions, bind the slope runoff regulation strategy entries with the vector space of the slope soil and water conservation area to generate watershed collaborative slope restoration sub-instructions, and merge the watershed collaborative waterfront restoration sub-instructions and the watershed collaborative slope restoration sub-instructions into the zonal differential land low-carbon restoration instruction set.

[0176] In the pre-defined low-carbon land restoration strategy library, strategy entries associated with "river buffer zone" are retrieved, and vegetation reconstruction strategy entries for the river buffer zone (including planting deep-rooted bank-stabilizing plants and constructing a composite buffer zone of trees, shrubs, and grasses) are selected. The vector of the waterfront degradation restoration section is spatially bound to the strategy entry to generate a watershed collaborative waterfront restoration sub-instruction.

[0177] In the pre-defined low-carbon land restoration strategy library, strategy entries associated with "slope runoff control" are retrieved, and slope runoff control strategy entries (including the construction of horizontal steps, fish-scale pits, and the implementation of contour tillage) are selected. The slope soil and water conservation zone vector is spatially bound to the strategy entry to generate a watershed collaborative slope restoration sub-instruction.

[0178] The two sub-instructions mentioned above are merged into the low-carbon restoration instruction set for zonal differences in land, which will be used to guide ecological restoration operations that are coordinated with hydrological connectivity.

[0179] Step S710: Obtain the multi-period land use change survey vector set of the target geographic area during the historical continuous monitoring period. The multi-period land use change survey vector set includes the land use type patch boundaries and corresponding land use type codes at different survey time points. Based on the land use type codes, extract the forest land to non-forest land patch change vector, grassland to non-grassland patch change vector, and wetland to non-wetland patch change vector.

[0180] In this embodiment, land use change survey vector data of the target geographic area at multiple time points within a continuous historical monitoring period are acquired. The data for each time point includes land use type patch boundaries and corresponding land use type codes (e.g., "031" represents forest land, "032" represents shrubland, "041" represents natural grassland, "043" represents artificial grassland, "110" represents river surface, and "116" represents inland tidal flats). By comparing data from adjacent time points, patches whose codes change from forest land (031, 032) to non-forest land are extracted, generating forest-to-non-forest patch change vectors. Similarly, patches whose codes change from grassland (041, 043) to non-grassland, and wetland (110, 116) to non-wetland, are extracted, generating grassland-to-non-grassland patch change vectors and wetland-to-non-wetland patch change vectors, respectively.

[0181] Step S720: Spatially overlay the forest-to-non-forest patch change vector, the grassland-to-non-grass patch change vector, and the wetland-to-non-wetland patch change vector with the carbon storage time differentiation map, respectively. Calculate the average loss magnitude of the carbon storage time change raster within the forest-to-non-forest patch change vector coverage area, the average loss magnitude of the carbon storage time change raster within the grassland-to-non-grass patch change vector coverage area, and the average loss magnitude of the carbon storage time change raster within the wetland-to-non-wetland patch change vector coverage area. Input the average loss magnitude of the forest-to-non-forest patch change vector, the average loss magnitude of the grassland-to-non-grass patch change vector, and the average loss magnitude of the wetland-to-non-wetland patch change vector into the ecological type conversion carbon loss sensitivity ranking device. Based on the magnitude of the average loss magnitude, sort the sensitivity levels of various patch change vectors to generate an ecological type conversion carbon loss sensitivity ranking list.

[0182] Using the forest-to-non-forest patch change vector as a mask, the carbon storage time-change raster map C_time_change is cropped to extract the carbon storage time-change values ​​ΔC_time(i,j) for all pixels within the changed patch coverage area. The average of these values ​​is calculated as the average loss magnitude A_loss_forest of the forest-to-non-forest patch change vector. Similarly, the average loss magnitude A_loss_grass of the grassland-to-non-grass patch change vector and the average loss magnitude A_loss_wetland of the wetland-to-non-wetland patch change vector are calculated respectively.

[0183] Input A_loss_forest, A_loss_grass, and A_loss_wetland into the ecological type conversion carbon loss sensitivity sorter, which sorts them from largest to smallest according to the average loss magnitude, and generates an ecological type conversion carbon loss sensitivity sort list. For example, the first place in the sort list is wetland to non-wetland (if A_loss_wetland is the largest), the second place is forest to non-forest, and the third place is grassland to non-grassland.

[0184] Step S730: Based on the ecological type conversion carbon loss sensitivity ranking list, select the target ecological conversion carbon loss hot spot patch set that is ranked within the set ranking range by average loss magnitude from the forest-to-non-forest patch change vector, grassland-to-non-grass patch change vector, and wetland-to-non-wetland patch change vector. Extract the spatial boundary coordinate sequence of the target ecological conversion carbon loss hot spot patch set. Perform a spatial consistency comparison between the spatial boundary coordinate sequence and the partitioned spatial unit set. Mark the first type of historical carbon loss hot spot coordinates located within the boundary coordinates of the carbon sink maintenance zone, the second type of historical carbon loss hot spot coordinates located within the boundary coordinates of the carbon source control zone, and the third type of historical carbon loss hot spot coordinates located within the boundary coordinates of the carbon sink enhancement zone in the spatial boundary coordinate sequence of the target ecological conversion carbon loss hot spot patch set.

[0185] Based on the sensitivity ranking list of carbon loss during ecological type transitions, the selection range is set to the top two ranked areas (i.e., the two ecological type transitions with the most severe carbon loss). The spatial boundary coordinate sequences of these patches are extracted from the corresponding patch change vectors to form a set of hotspot patches for carbon loss during the target ecological transition.

[0186] The spatial boundary coordinate sequence of the hotspot patch set is spatially overlaid with the set of partitioned spatial units. For each hotspot patch boundary coordinate, its positional relationship with the carbon sink maintenance zone, carbon source control zone, and carbon sink enhancement zone is determined. If the patch is located inside the boundary coordinates of the carbon sink maintenance zone, it is marked as a first-type historical carbon loss hotspot coordinate; if it is located inside the boundary coordinates of the carbon source control zone, it is marked as a second-type historical carbon loss hotspot coordinate; and if it is located inside the boundary coordinates of the carbon sink enhancement zone, it is marked as a third-type historical carbon loss hotspot coordinate.

[0187] Step S740: For the first type of historical carbon loss hotspot coordinates, extract historical carbon loss ecological restoration strategy entries from the preset land low-carbon restoration strategy library. The historical carbon loss ecological restoration strategy entries include an ecological reverse succession restoration operation that reverses the current land use type code of the area covered by the first type of historical carbon loss hotspot coordinates to the original forest land code, original grassland code, or original wetland code. For the second type of historical carbon loss hotspot coordinates, extract historical carbon loss source blocking strategy entries from the preset land low-carbon restoration strategy library. The historical carbon loss source blocking strategy entries include a source blocking management operation that permanently prohibits development and construction in the area covered by the second type of historical carbon loss hotspot coordinates. For the third type of historical carbon loss hotspot coordinates, extract historical carbon loss potential activation strategy entries from the preset land low-carbon restoration strategy library. The historical carbon loss potential activation strategy entries include a carbon sink function reactivation operation that rapidly increases soil organic matter in the area covered by the third type of historical carbon loss hotspot coordinates.

[0188] In the preset land low-carbon restoration strategy library, search for strategy entries associated with "ecological restoration of historical carbon loss" and select ecological reverse succession restoration operations: reversely restore the current land use type (such as cultivated land or construction land) of the first type of historical carbon loss hotspot coordinate coverage area to the historical forest land, grassland or wetland type through engineering measures (such as returning cultivated land to forest or returning construction land to wetland).

[0189] In the pre-set land low-carbon restoration strategy library, search for strategy entries associated with "historical carbon loss source sealing", and select source sealing management operations: implement permanent prohibition of development and construction in the area covered by the coordinates of the second type of historical carbon loss hotspots, and stop all human activities that may cause carbon loss.

[0190] In the pre-set land low-carbon restoration strategy library, search for strategy entries associated with "activation of historical carbon loss potential" and select carbon sink function reactivation operations: implement soil organic matter rapid improvement measures (such as applying organic fertilizer, planting green manure crops, introducing soil organisms such as earthworms) in the third category of historical carbon loss hotspot coordinate coverage areas.

[0191] Step S750: Bind the coordinates of the first type of historical carbon loss hotspots to the entry space of the historical carbon loss ecological restoration strategy to generate a historical loss restoration and repair sub-instruction; bind the coordinates of the second type of historical carbon loss hotspots to the entry space of the historical carbon loss source containment strategy to generate a historical loss containment and repair sub-instruction; bind the coordinates of the third type of historical carbon loss hotspots to the entry space of the historical carbon loss potential activation strategy to generate a historical loss activation and repair sub-instruction.

[0192] The coordinate sequences of the first type of historical carbon loss hotspots are correlated with ecological reverse succession restoration operations to generate historical loss restoration sub-instructions. The coordinate sequences of the second type of historical carbon loss hotspots are correlated with source control management operations to generate historical loss control restoration sub-instructions. The coordinate sequences of the third type of historical carbon loss hotspots are correlated with carbon sink function reactivation operations to generate historical loss activation restoration sub-instructions.

[0193] Step S760: Merge the historical loss recovery and repair sub-instruction, the historical loss blocking and repair sub-instruction, and the historical loss activation and repair sub-instruction into the zonal differential land low-carbon restoration instruction set, which is used to guide differentiated low-carbon restoration intervention operations for carbon loss hotspots caused by ecological type transformation in historical periods.

[0194] The historical loss recovery and repair sub-instructions, historical loss blocking and repair sub-instructions, and historical loss activation and repair sub-instructions generated in step S750 are added to the zonal differential land low-carbon repair instruction set generated in step S148. The resulting complete instruction set not only includes conventional repair instructions for carbon sink maintenance zones, carbon source control zones, and carbon sink enhancement zones, but also includes special repair instructions for historical carbon loss hotspots, achieving a comprehensive response and differentiated and precise repair of the spatiotemporal differentiation of carbon reserves in the target geographical area.

[0195] In one exemplary embodiment, a land zoning low-carbon remediation decision-making system based on the spatiotemporal differentiation of ecosystem carbon storage is provided. This system can be a terminal, server, etc., and its internal structure diagram can be as follows: Figure 2As shown, the system specifically includes a processor, memory, input / output interface, communication interface, display unit, and input device. The processor, memory, and input / output interface are connected via a system bus, and the communication interface, display unit, and input device are also connected to the system bus via the input / output interface. The processor provides computing and control capabilities. The memory includes non-volatile storage media and internal memory. The non-volatile storage media stores the operating system and computer programs. The internal memory provides the environment for the operation of the operating system and computer programs in the non-volatile storage media. The input / output interface is used for exchanging information between the processor and external devices. The communication interface is used for wired or wireless communication with external terminals; wireless communication can be achieved through Wi-Fi, mobile cellular networks, near-field communication, or other technologies. When the computer program is executed by the processor, it implements a land zoning low-carbon remediation decision-making method based on the spatiotemporal differentiation of ecosystem carbon storage. The display unit is used to form a visually visible image and can be a display screen, projection device, or virtual reality imaging device. The display screen can be an LCD screen or an e-ink screen. The input device can be a touch layer covering the display screen, or a button, trackball, or touchpad set on the shell of the land zoning low-carbon remediation decision system based on the spatiotemporal differentiation of ecosystem carbon storage, or an external keyboard, touchpad, or mouse, etc.

[0196] It should be noted that, in order to simplify the description of the present invention and thus help to understand one or more embodiments of the invention, multiple features may sometimes be grouped into one embodiment, drawing or description thereof in the foregoing description of the embodiments of the present invention.

Claims

1. A land partition low-carbon restoration decision-making method based on the spatiotemporal differentiation of ecosystem carbon storage, characterized in that, The method includes: Acquire a target remote sensing image set and a ground sample plot survey set for the target geographic area within a historical continuous monitoring period. The target remote sensing image set includes multispectral image gratings and radar image gratings with spatial coordinate labels, and the ground sample plot survey set includes measured records of vegetation biomass and soil organic carbon with spatial coordinate labels. Carbon storage spatiotemporal differentiation features are extracted from the target remote sensing image set and the ground sample plot survey set to obtain a carbon storage spatial differentiation map reflecting the differences in the spatial distribution of carbon storage and a carbon storage temporal differentiation map reflecting the differences in the temporal variation of carbon storage. Land management zones are delineated on the carbon storage spatial differentiation map and the carbon storage temporal differentiation map to obtain a set of zoning spatial units consisting of the boundary coordinates of the carbon sink maintenance zone, the boundary coordinates of the carbon source control zone, and the boundary coordinates of the carbon sink enhancement zone. Based on the analysis of the driving forces of carbon storage changes in different management zones in the spatial unit set, the corresponding restoration entries in the preset land low-carbon restoration strategy library are matched for the boundary coordinate coverage areas of the carbon sink maintenance zone, the boundary coordinate coverage areas of the carbon source control zone, and the boundary coordinate coverage areas of the carbon sink enhancement zone, respectively, and a set of land low-carbon restoration instructions with zoned differences is generated. The set of instructions for low-carbon restoration of land with zoning differences is sent to the land management terminal to trigger vegetation structure adjustment, soil carbon pool protection, and ecological network connectivity enhancement operations for different zones.

2. The land zoning low-carbon remediation decision-making method based on the spatiotemporal differentiation of ecosystem carbon storage as described in claim 1, characterized in that, The step of extracting spatiotemporal differentiation features of carbon storage from the target remote sensing image set and the ground sample plot survey set to obtain a carbon storage spatial differentiation map reflecting the spatial distribution differences of carbon storage and a carbon storage temporal differentiation map reflecting the temporal variation differences of carbon storage includes: The vegetation index is inverted on the multispectral image raster of the target remote sensing image set to obtain the canopy greenness index raster that characterizes the intensity of vegetation photosynthesis, and the backscattering coefficient is interpreted on the radar image raster to obtain the canopy height index raster and the canopy water content index raster that characterize the three-dimensional structure of the ground surface. The canopy greenness index raster, the canopy height index raster, and the canopy water content index raster are input into the first feature layer of the carbon storage spatial heterogeneous modeling network. The texture gradient of the canopy greenness index raster, the structural change of the canopy height index raster, and the water heterogeneity information of the canopy water content index raster are extracted in parallel through the multi-scale convolution operator in the first feature layer. The raster is then stitched along the feature channel dimension to generate a target remote sensing fusion feature map. The target remote sensing fusion feature map is input into the spatial correlation layer of the carbon storage spatial heterogeneous modeling network. The spatial self-attention operator in the spatial correlation layer is used to calculate the covariance matrix between any points in the target remote sensing fusion feature map. The target remote sensing fusion feature map is then weighted and summed based on the covariance matrix to generate a global spatial correlation feature map. The spatial coordinates of the measured vegetation biomass records and measured soil organic carbon records in the ground sample plot survey set are aligned with the coordinate grid of the global spatial association feature map. The first feature vector of the corresponding spatial coordinates of the measured vegetation biomass records and the second feature vector of the corresponding spatial coordinates of the measured soil organic carbon records are extracted from the global spatial association feature map. A multi-sensor branch for vegetation biomass mapping and a multi-sensor branch for soil organic carbon mapping are constructed. The first feature vector is input into the multi-sensor branch for vegetation biomass mapping and nonlinear transformation is performed to generate a spatial distribution map of vegetation biomass. The second feature vector is input into the multi-sensor branch for soil organic carbon mapping and nonlinear transformation is performed to generate a spatial distribution map of soil organic carbon. The spatial distribution map of vegetation biomass and the spatial distribution map of soil organic carbon are spatially overlaid using a grid overlay operator to generate a total spatial map of ecosystem carbon storage, and the gradient of carbon storage change between adjacent grids in the total spatial map of ecosystem carbon storage is calculated to generate a spatial differentiation map of carbon storage. Extract the initial carbon storage spatial total map corresponding to the starting time point and the final carbon storage spatial total map corresponding to the ending time point within the historical continuous monitoring period. Then, use the time series change operator to perform grid-by-grid difference calculation on the initial carbon storage spatial total map and the final carbon storage spatial total map to generate a carbon storage time change raster map. The direction of change of the carbon storage over time is determined, and the polarity of change of each grid in the carbon storage over time grid is marked. The polarity of change includes carbon storage accumulation state and carbon storage loss state. A carbon storage time differentiation map is generated based on the combination characteristics of the polarity of change of continuous grids. The carbon storage spatial differentiation map is input into the spectrum edge enhancement filter for convolution filtering. The carbon storage spatial differentiation map and carbon storage temporal differentiation map after convolution filtering are normalized and stretched. The numerical ranges of the carbon storage spatial differentiation map and the carbon storage temporal differentiation map after convolution filtering are mapped to a preset display dynamic range to obtain the carbon storage spatial differentiation map and carbon storage temporal differentiation map.

3. The land zoning low-carbon remediation decision-making method based on the spatiotemporal differentiation of ecosystem carbon storage as described in claim 1, characterized in that, The process of delineating land management zones based on the spatial and temporal differentiation maps of carbon reserves yields a set of spatial units comprising the boundary coordinates of carbon sink maintenance zones, carbon source control zones, and carbon sink enhancement zones, including: Extract the carbon storage change gradient in the carbon storage spatial differentiation map, and mark the grids with carbon storage change gradients higher than a preset gradient threshold as spatial abrupt change candidate grids. Extract the closed outer contour of the continuously clustered spatial abrupt change candidate grids as carbon storage spatial fault zone line vectors. Extract the change polarity in the carbon storage time differentiation map, and extract the patch boundary formed by the continuous grids with change polarity of carbon storage loss state as carbon loss hot spot area boundary line vectors. Spatial topological superposition is performed on the spatial fault zone vector of the carbon storage and the boundary line vector of the carbon loss hotspot area. The vector surface unit jointly divided by the spatial fault zone vector of the carbon storage and the boundary line vector of the carbon loss hotspot area is extracted as the initial partition surface unit, and each initial partition surface unit is assigned a unique partition identifier. The average value of the spatial differentiation map of carbon storage within each initial partition surface unit is extracted as the spatial heterogeneity attribute value of that initial partition surface unit. The proportion of the carbon storage loss state grid area in the temporal differentiation map of carbon storage within each initial partition surface unit is extracted as the temporal variation risk value of that initial partition surface unit. A two-dimensional partition decision feature vector is constructed based on the spatial heterogeneous attribute values ​​and temporal change risk values ​​of the initial partition surface unit. The two-dimensional partition decision feature vector is input into the partition decision tree ensemble model based on gradient boosting. The probability score of the initial partition surface unit belonging to the carbon sink maintenance class, carbon source management class, and carbon sink enhancement class is calculated through multiple decision subtrees in the partition decision tree ensemble model. For each initial partition surface unit, compare the category probability scores of its carbon sink maintenance category, carbon source control category, and carbon sink enhancement category. Select the category with the highest category probability score as the initial partition label of the initial partition surface unit. Merge the surface features of initial partition surface units with the same initial partition label and spatial adjacency to generate carbon sink maintenance area, carbon source control area, and carbon sink enhancement area with complete geographic entities. The outer boundary coordinates of the carbon sink maintenance area are extracted as the boundary coordinates of the carbon sink maintenance area, the outer boundary coordinates of the carbon source control area are extracted as the boundary coordinates of the carbon source control area, and the outer boundary coordinates of the carbon sink enhancement area are extracted as the boundary coordinates of the carbon sink enhancement area. The boundary coordinates of the carbon sink maintenance area, the carbon source control area, and the carbon sink enhancement area are input into the partitioned geometric regularization processor. The geometric closing operation operator fills the tiny holes inside each boundary coordinate and the geometric opening operation operator removes the small protrusions outside each boundary coordinate. The boundary coordinates of the carbon sink maintenance area, carbon source control area, and carbon sink enhancement area after geometric processing are subjected to coordinate thinning. The coordinate thinned boundary coordinates of the carbon sink maintenance area, carbon source control area, and carbon sink enhancement area are combined to construct a set of partitioned spatial units. The association between the boundary coordinates of the carbon sink maintenance area and the carbon sink maintenance category label, the association between the boundary coordinates of the carbon source control area and the carbon source control category label, and the association between the boundary coordinates of the carbon sink enhancement area and the carbon sink enhancement category label are stored in the set of partitioned spatial units.

4. The land zoning low-carbon remediation decision-making method based on the spatiotemporal differentiation of ecosystem carbon storage as described in claim 1, characterized in that, The step involves analyzing the driving forces of carbon storage changes in different management zones within the spatial unit set, matching corresponding restoration entries from a pre-defined land low-carbon restoration strategy library to the boundary coordinate coverage areas of the carbon sink maintenance zone, the carbon source control zone, and the carbon sink enhancement zone, respectively, and generating a set of land low-carbon restoration instructions based on zone differences, including: The carbon storage time differentiation map is cropped based on the boundary coordinates of the carbon sink maintenance area. The spatial distribution pattern of the carbon storage accumulation state grid with changing polarity is extracted within the boundary coordinate coverage area of ​​the carbon sink maintenance area. The area continuous expansion direction vector and area continuous expansion rate of the carbon storage accumulation state grid are calculated. Based on the continuous expansion direction vector and continuous expansion rate of the area, a set of strategy entries associated with carbon sink maintenance category tags is retrieved from the preset land low-carbon remediation strategy library. Carbon sink maintenance strategy entries with execution intensity positively correlated with the continuous expansion rate of the area are selected from the strategy entry set. The carbon sink maintenance strategy entries include maintaining the natural succession process and limiting the scope of human disturbance. The boundary coordinates of the carbon sink maintenance area and the carbon sink maintenance strategy entries are spatially bound together to generate a carbon sink maintenance area repair instruction package carrying spatial range constraints and strategy content descriptions. The spatial range constraints include the entire coordinate sequence of the boundary coordinates of the carbon sink maintenance area. Based on the boundary coordinates of the carbon source control area, the carbon storage time differentiation map is cropped. The coordinates of the spatial aggregation center with the changing polarity of carbon storage loss state grid within the boundary coordinate coverage area of ​​the carbon source control area are extracted. The spatial influence radiation radius of the spatial aggregation center coordinates to the surrounding area is calculated. A circular carbon source influence core area is constructed with the spatial aggregation center coordinates as the center and the spatial influence radiation radius as the radius. The area within the boundary coordinates of the carbon source control area located outside the circular carbon source influence core area is designated as the carbon source influence buffer diffusion area, generating a two-layer spatial structure description of the carbon source control area. Based on the description of the dual-layer spatial structure, a set of strategy entries associated with carbon source management tags is retrieved from the preset land low-carbon remediation strategy library. From the set of strategy entries, carbon source blocking strategy entries for the core area of ​​the circular carbon source impact and carbon source conversion strategy entries for the buffer and diffusion area of ​​the carbon source impact are selected. The carbon source blocking strategy entries include removing carbon emission activities and cutting off carbon loss pathways. The carbon source conversion strategy entries include replacing the surface cover and breaking the soil sealing layer. The boundary coordinates of the carbon source control area, the carbon source blocking strategy entries and the carbon source conversion strategy entries are spatially layered and bound to generate a first sub-instruction carrying the core area range constraints and the carbon source blocking strategy, and a second sub-instruction carrying the buffer diffusion area range constraints and the carbon source conversion strategy. The first sub-instruction and the second sub-instruction are combined to form a carbon source control area repair instruction package. Based on the boundary coordinates of the carbon sink enhancement zone, the spatial differentiation map of carbon storage and the temporal differentiation map of carbon storage are cropped. The coordinates of potential patches to be activated within the boundary coordinates of the carbon sink enhancement zone are identified, where the spatial differentiation map values ​​of carbon storage are lower than those of the surrounding areas and the polarity of the temporal differentiation map of carbon storage is in a stable state. Based on the coordinates of potential patches to be activated, a set of strategy entries associated with carbon sink enhancement category tags is retrieved from a preset land low-carbon remediation strategy library. From the set of strategy entries, carbon sink enhancement strategy entries containing vegetation type optimization configuration and soil exogenous organic supplementation are selected. The vegetation function combination features in the vegetation type optimization configuration are adaptively matched with the soil organic carbon spatial distribution map values ​​of the potential patches to be activated. Spatially bind the boundary coordinates of the carbon sink enhancement zone, the carbon sink enhancement strategy entries, and the coordinates of the potential patches to be activated to generate a carbon sink enhancement zone remediation instruction package carrying the coordinate sequence of potential patches to be activated and the carbon sink enhancement strategy entries. Merge the carbon sink maintenance zone remediation instruction package, the carbon source control zone remediation instruction package, and the carbon sink enhancement zone remediation instruction package to form a low-carbon remediation instruction set for land with zoning differences.

5. The land zoning low-carbon remediation decision-making method based on the spatiotemporal differentiation of ecosystem carbon storage according to claim 1, characterized in that, The method further includes: After the execution of the low-carbon remediation instruction set for zonal differences in land on the land management terminal, a continuous monitoring feedback remote sensing set and a continuous monitoring feedback ground set are obtained. The continuous monitoring feedback remote sensing set includes feedback multispectral image raster and feedback radar image raster with spatial coordinate markers after the remediation is performed. The continuous monitoring feedback ground set includes feedback vegetation biomass measurement records and feedback soil organic carbon measurement records with spatial coordinate markers after the remediation is performed. The spatiotemporal differentiation features of the feedback carbon storage are extracted from the continuous monitoring feedback remote sensing set and the continuous monitoring feedback ground set to obtain a feedback carbon storage spatial differentiation map reflecting the differences in the spatial distribution of carbon storage after restoration and a feedback carbon storage temporal differentiation map reflecting the differences in the temporal change of carbon storage after restoration. Based on the boundary coordinates of the carbon sink maintenance area, the spatial differentiation map and temporal differentiation map of the feedback carbon storage are cropped, and the average net accumulation of feedback carbon storage and the slope change of the accumulation rate of feedback carbon storage within the area covered by the boundary coordinates of the carbon sink maintenance area are calculated. Based on the boundary coordinates of the carbon source control area, the spatial differentiation map of the feedback carbon storage and the temporal differentiation map of the feedback carbon storage are cropped, and the area shrinkage ratio of the carbon storage loss state grid with changing polarity in the boundary coordinate coverage area of ​​the carbon source control area and the fragmentation index change of the carbon storage loss state grid at the boundary of the core area are extracted. Based on the boundary coordinates of the carbon sink enhancement zone, the spatial differentiation map of the feedback carbon storage and the temporal differentiation map of the feedback carbon storage are cropped. The improvement of the feedback carbon storage value at the coordinates of the potential patch to be activated compared with that before the repair is detected, and it is analyzed whether the polarity of its change has changed from a stationary state to a carbon storage accumulation state. A carbon sink maintenance zone remediation effectiveness deviation analyzer is constructed. The difference between the slope change of the feedback carbon storage accumulation rate within the boundary coordinate coverage area of ​​the carbon sink maintenance zone and the expected rate range preset by the carbon sink maintenance strategy item is compared to generate a carbon sink maintenance effectiveness deviation vector. A deviation analyzer for the remediation effectiveness of carbon source control areas is constructed. The changes in the area shrinkage ratio and fragmentation index of the carbon storage loss state grid within the boundary coordinate coverage area of ​​the carbon source control area are compared with the expected shrinkage range preset by the carbon source blocking strategy item and the expected fragmentation reduction range preset by the carbon source conversion strategy item, respectively, to generate a deviation vector for the effectiveness of carbon source control. A carbon sink enhancement zone restoration effectiveness deviation analyzer is constructed. The difference between the increase in the feedback carbon storage value at the coordinates of the potential patch to be activated and the expected increase range preset by the carbon sink enhancement strategy item is compared. The polarity of the change at the coordinates of the potential patch to be activated is judged to be consistent with the expected carbon storage accumulation state. A carbon sink potential activation effectiveness deviation vector is generated. The execution intensity of the carbon sink maintenance strategy item is adjusted according to the carbon sink maintenance effectiveness deviation vector. The spatial range of the carbon source blocking strategy item and the carbon source conversion strategy item is adjusted according to the carbon source control effectiveness deviation vector. The vegetation function combination characteristics in the carbon sink enhancement strategy item are replaced or the proportion is adjusted according to the carbon sink potential activation effectiveness deviation vector. The adjusted carbon sink maintenance strategy items, adjusted carbon source blocking strategy items, carbon source conversion strategy items, and adjusted carbon sink enhancement strategy items are repackaged into an iterative optimization land low-carbon remediation instruction set based on zoning differences, and sent to the land management terminal to update the land low-carbon remediation operations in progress.

6. The land zoning low-carbon remediation decision-making method based on the spatiotemporal differentiation of ecosystem carbon storage according to claim 1, characterized in that, The method further includes: Acquire a set of climate environmental factors and a set of anthropogenic activity factors affecting the dynamic changes of carbon storage within a target geographic region. The climate environmental factor set includes a raster of annual average temperature distribution, an annual average precipitation distribution, and a raster of total solar radiation distribution with spatial coordinate labels. The anthropogenic activity factor set includes a raster of land use type, a raster of population activity thermal data, and a raster of transportation network density with spatial coordinate labels. Input the annual average temperature distribution raster, the annual average precipitation distribution raster, and the raster of total solar radiation distribution into a climate factor nonlinear coupler. The polynomial kernel function maps the annual average temperature distribution raster value, annual average precipitation distribution raster value, and total solar radiation distribution raster value to a high-dimensional climate space and calculates the climate stress comprehensive index raster for each spatial location. The land use type raster, the population activity thermal raster, and the traffic network density raster are input into the anthropogenic disturbance intensity quantizer. The weighted summation operator in the anthropogenic disturbance intensity quantizer performs weighted superposition of the disturbance coefficient values ​​of different land types, the normalized values ​​of the population activity thermal raster, and the normalized values ​​of the traffic network density raster to generate the anthropogenic disturbance intensity comprehensive index raster. The climate stress composite index grid and the anthropogenic disturbance intensity composite index grid are input into the carbon storage change driving force discrimination network. The first Granger causality between the climate stress composite index grid and the carbon storage time change grid, and the second Granger causality between the anthropogenic disturbance intensity composite index grid and the carbon storage time change grid are calculated through the causal inference layer in the carbon storage change driving force discrimination network. For the grids in the carbon storage time differentiation map where the change polarity is the carbon storage accumulation state, the carbon storage accumulation state driver is classified as either climate-dominated or anthropogenic suppression and weakening driver based on the relative magnitude of the first Granger causality and the second Granger causality, thus generating a carbon sink formation driving force label distribution map. For the grids in the carbon storage time differentiation map where the change polarity is carbon storage loss state, the carbon storage loss state driver is classified as either climate stress-driven or anthropogenic activity-driven based on the relative magnitude of the first Granger causality and the second Granger causality, thus generating a carbon source formation driving force label distribution map. The carbon sink formation driving force label distribution map is spatially overlaid with the partitioned spatial unit set to analyze the area ratio of climate-dominant driving grids and anthropogenic suppression weakening driving grids within the boundary coordinate coverage area of ​​the carbon sink maintenance zone, and the grid with the larger ratio is identified as the main driving force for carbon storage changes in the carbon sink maintenance zone. Similarly, the carbon source formation driving force label distribution map is spatially overlaid with the partitioned spatial unit set to analyze the area ratio of climate stress-dominant driving grids and anthropogenic activity enhancement driving grids within the boundary coordinate coverage area of ​​the carbon source control zone, and the grid with the larger ratio is identified as the main driving force for carbon storage changes in the carbon source control zone. Finally, the carbon sink formation driving force label distribution map and the carbon source formation driving force label distribution map are spatially overlaid with the partitioned spatial unit set to analyze the area distribution characteristics of climate-dominant driving grids, anthropogenic suppression weakening driving grids, climate stress-dominant driving grids, and anthropogenic activity enhancement driving grids within the boundary coordinate coverage area of ​​the carbon sink enhancement zone, and the main driving force for carbon storage changes in the carbon sink enhancement zone is generated. The main driving forces of carbon storage changes in carbon sink maintenance areas, carbon storage changes in carbon source control areas, and carbon storage enhancement areas are added to the spatial unit set of the zoning, serving as input parameters for matching strategy entries in the instruction set for low-carbon restoration of land with zoning differences.

7. The land zoning low-carbon remediation decision-making method based on the spatiotemporal differentiation of ecosystem carbon storage according to claim 1, characterized in that, The method further includes: Obtain a set of field patrol and survey tracks for the target geographic area during a continuous historical monitoring period. The set of field patrol and survey tracks includes vegetation community structure records and soil profile morphology records collected along the patrol route, which have spatial coordinate markers and collection time markers. Key information is extracted from the vegetation community structure record. The dominant tree species names, shrub coverage descriptions, and herb diversity descriptions are extracted from the vegetation community structure record using the named entity recognition operator of the natural language processing component. The dominant tree species names, shrub coverage descriptions, and herb diversity descriptions are then converted into structured vegetation vertical structure vectors. Similarly, key information is extracted from the soil profile morphology record. The humus layer thickness description, soil color description, and soil compaction description are extracted from the soil profile morphology record using the named entity recognition operator of the natural language processing component. The humus layer thickness description, soil color description, and soil compaction description are then converted into structured soil profile trait vectors. The structured vegetation vertical structure vector and the structured soil profile trait vector are spatially associated with the grid in the carbon storage spatial differentiation map according to spatial coordinate labels, forming a field-enhanced carbon storage grid sequence containing the vegetation vertical structure vector and the soil profile trait vector. A correlation analysis network was constructed between vegetation structure complexity and carbon storage stability. The structured vegetation vertical structure vectors from the field enhanced carbon storage grid sequence were input into the structure complexity layer of the correlation analysis network. The species evenness of dominant tree species, the cover continuity described by shrub cover, and the species richness described by herb diversity were calculated. Each indicator was normalized and weighted to generate a comprehensive score map of vegetation structure complexity. A correlation analysis network was also constructed between soil properties and carbon pool storage persistence. The structured soil profile trait vectors from the field enhanced carbon storage grid sequence were input into the carbon pool stability layer of the correlation analysis network. The cumulative organic matter thickness was calculated based on the humus layer thickness, the degree of organic matter darkening was calculated based on the soil color, and the soil bulk density was estimated based on the soil compaction. Each indicator was normalized and weighted to generate a comprehensive score map of soil carbon pool persistence. Spatial raster correlation analysis is performed between the vegetation structure complexity comprehensive score map and the soil carbon pool persistence comprehensive score map and the carbon storage time differentiation map, respectively. The first spatial correlation coefficient matrix of carbon storage accumulation rate in the vegetation structure complexity comprehensive score map and the carbon storage time differentiation map, and the second spatial correlation coefficient matrix of carbon storage accumulation rate in the soil carbon pool persistence comprehensive score map and the carbon storage time differentiation map are calculated. A first spatial heat map of the influence of vegetation complexity on carbon storage time differentiation is generated based on the first spatial correlation coefficient matrix, and a second spatial heat map of the influence of soil persistence on carbon storage time differentiation is generated based on the second spatial correlation coefficient matrix. Spatially crop the first spatial heat map and the boundary coordinates of the carbon sink enhancement zone to identify a first priority remediation sub-zone coordinate sequence in which the comprehensive score of vegetation structure complexity and the carbon storage accumulation rate are positively correlated within the coverage area of ​​the carbon sink enhancement zone boundary coordinates; and spatially crop the second spatial heat map and the boundary coordinates of the carbon sink enhancement zone to identify a second priority remediation sub-zone coordinate sequence in which the comprehensive score of soil carbon pool persistence and the carbon storage accumulation rate are positively correlated within the coverage area of ​​the carbon sink enhancement zone boundary coordinates, and supplement the first priority remediation sub-zone coordinate sequence and the second priority remediation sub-zone coordinate sequence into the partitioned spatial unit set to refine the spatial placement of carbon sink enhancement zone strategy items.

8. The land zoning low-carbon remediation decision-making method based on the spatiotemporal differentiation of ecosystem carbon storage according to claim 1, characterized in that, The method further includes: Obtain the current land use planning atlas and ecological protection red line boundary set corresponding to different management zones within the target geographical area. The current land use planning atlas contains legally binding planning land category patch boundaries and planning land category use codes. The ecological protection red line boundary set contains legally binding ecological red line boundaries and ecological red line control levels. Spatial conflict detection is performed between the boundary coordinates of the planned land use patch in the current land use planning atlas and the boundary coordinates of the carbon sink maintenance area in the zoning spatial unit atlas, and the first set of spatial conflict coordinates that overlaps with the boundary coordinates of the carbon sink maintenance area and whose planned land use code is a permitted development and construction class is extracted from the boundary of the planned land use patch. Spatial conflict detection is performed between the planning land use patch boundaries in the current land use planning atlas and the carbon source control area boundary coordinates in the zoning spatial unit atlas. The second set of spatial conflict coordinates is extracted from the planning land use patch boundaries that overlap with the carbon source control area boundary coordinates and whose planning land use code is the carbon emission enhancement category. Spatial conflict detection is performed between the planning land type patch boundaries in the current land use planning atlas and the carbon sink enhancement zone boundary coordinates in the zoning spatial unit atlas. A third set of spatial conflict coordinates is extracted from the planning land type patch boundaries that overlap with the carbon sink enhancement zone boundary coordinates and whose planning land use code is prohibited for ecological restoration. Spatial overlay analysis is performed on the ecological red line boundary of the ecological protection red line boundary set and the carbon sink maintenance area boundary coordinate of the partitioned spatial unit set. The coordinates of the first ecological void area that is not covered by the ecological red line boundary within the carbon sink maintenance area boundary coordinate are extracted, and the area size and spatial dispersion of the first ecological void area coordinate are calculated. Spatial overlay analysis is performed on the ecological red line boundary coordinates of the ecological protection red line boundary set and the carbon source control area boundary coordinates of the partitioned spatial unit set. The coordinates of the overlapping area inside the carbon source control area boundary coordinates and the ecological red line boundary are extracted, and the compatibility index between the carbon source blocking strategy items and the ecological red line control level in the overlapping area coordinates is calculated. Spatial overlay analysis is performed on the ecological red line boundary of the ecological protection red line boundary set and the carbon sink enhancement zone boundary coordinate of the partitioned spatial unit set. The coordinates of the ecological network connected to the ecological red line boundary outside the carbon sink enhancement zone boundary coordinates are extracted, and the minimum spatial distance distribution between the coordinates of the connected to the carbon sink enhancement zone boundary coordinates and the carbon sink enhancement zone boundary coordinates is analyzed. Based on the first spatial conflict coordinate set, the second spatial conflict coordinate set, and the third spatial conflict coordinate set, a land use code adjustment suggestion is generated for the current land use planning atlas. The land use code adjustment suggestion includes changing the permitted development and construction category to the prohibited development and construction category, changing the permitted carbon emission enhancement category to the restricted carbon emission enhancement category, and changing the prohibited ecological restoration category to the permitted ecological restoration category. Based on the area size and spatial dispersion of the first ecological void area coordinates, a boundary expansion proposal is generated for the ecological protection red line boundary set. The boundary expansion proposal includes the spatial expansion direction and area size of incorporating the first ecological void area coordinates into the ecological red line boundary. Based on the compatibility index calculation results of the overlapping area coordinates and the minimum spatial distance distribution between the coordinates of the connected areas to be strengthened and the boundary coordinates of the carbon sink enhancement area, the management level optimization suggestions for the ecological protection red line boundary set are generated. The use code adjustment suggestions, boundary expansion suggestions and management level optimization suggestions are added as auxiliary information to the low-carbon restoration instruction set for zoning differences in land.

9. The land zoning low-carbon remediation decision-making method based on the spatiotemporal differentiation of ecosystem carbon storage according to claim 1, characterized in that, The method further includes: Obtain the regional energy consumption structure set and the regional transportation mode set within and around the target geographical area. The regional energy consumption structure set includes the proportion of fossil energy consumption and the proportion of clean energy consumption marked by administrative divisions. The regional transportation mode set includes the motor vehicle travel modal share and the non-motor vehicle travel modal share marked by administrative divisions. The proportion of fossil energy consumption and the share of motor vehicle travel are allocated to the grid scale according to the spatial weight of nighttime light intensity, generating gridded fossil carbon emission source intensity maps and gridded traffic carbon emission source intensity maps. The gridded fossil carbon emission source intensity map and the gridded transportation carbon emission source intensity map are spatially overlaid to generate a regional carbon emission source comprehensive intensity map, and the grids in the regional carbon emission source comprehensive intensity map whose comprehensive carbon emission source intensity exceeds a preset intensity threshold are extracted as high carbon emission source grid sets. Calculate the spatial Euclidean distance matrix between the spatial centroid coordinates of the high carbon emission source grid set and the boundary coordinates of the carbon source control area. Mark the boundary coordinates of the carbon source control area that are less than a preset distance threshold in the spatial Euclidean distance matrix as the boundary coordinates of the near source control area, and mark the boundary coordinates of the carbon source control area that are greater than the preset distance threshold as the boundary coordinates of the far source control area. For the near-source control zone boundary coordinate coverage area, extract industrial supporting strategy items for energy structure optimization from the preset land low-carbon remediation strategy library. The industrial supporting strategy items include regional energy substitution to increase the proportion of clean energy consumption and public transportation guidance to reduce the share of motor vehicle travel. For the boundary coordinates of the remote source control area, ecological barrier strategy entries targeting the blocking of carbon transport pathways are extracted from the preset land low-carbon restoration strategy library. The ecological barrier strategy entries include the windward interception construction of high carbon absorption efficiency vegetation belts configured in the windward direction of the boundary coordinates of the remote source control area and the internal carbon sink enhancement by constructing multi-layer vegetation three-dimensional configuration inside the boundary coordinates of the remote source control area. Spatially bind the boundary coordinates of the near-source control area with the industrial supporting strategy items to generate regional collaborative near-source restoration sub-instructions; spatially bind the boundary coordinates of the far-source control area with the ecological barrier strategy items to generate ecological barrier far-source restoration sub-instructions; and merge the regional collaborative near-source restoration sub-instructions and the ecological barrier far-source restoration sub-instructions into the zonal differential land low-carbon restoration instruction set. Obtain the hydrological connectivity set within the watershed of the target geographic region. The hydrological connectivity set includes a linear vector of the river system with river hierarchical coding and a sub-watershed area vector with watershed boundary coding. Calculate the riparian vegetation cover integrity index of the linear vector of the river system and the soil erosion sensitivity index of the sub-watershed area vector. Spatial overlay analysis is performed on the boundary coordinates of the carbon sink maintenance area, the linear vector of the river system, and the planar vector of the sub-basin. Vectors of riparian degradation restoration sections with riparian vegetation cover integrity index below a preset integrity threshold are extracted from the boundary coordinates of the carbon sink maintenance area. Vectors of slope soil and water conservation areas with soil and water loss sensitivity index above a preset sensitivity threshold are also extracted from the boundary coordinates of the carbon sink maintenance area. For the vector of the waterfront degradation restoration section, extract the river buffer zone vegetation reconstruction strategy entries from the preset land low-carbon restoration strategy library; for the vector of the slope soil and water conservation area, extract the slope runoff regulation strategy entries from the preset land low-carbon restoration strategy library; bind the river buffer zone vegetation reconstruction strategy entries with the vector space of the waterfront degradation restoration section to generate watershed collaborative waterfront restoration sub-instructions; bind the slope runoff regulation strategy entries with the vector space of the slope soil and water conservation area to generate watershed collaborative slope restoration sub-instructions; and merge the watershed collaborative waterfront restoration sub-instructions and the watershed collaborative slope restoration sub-instructions into the zonal differential land low-carbon restoration instruction set.

10. The land zoning low-carbon remediation decision-making method based on the spatiotemporal differentiation of ecosystem carbon storage according to claim 1, characterized in that, The method further includes: Obtain a multi-period land use change survey vector set for the target geographic area during a continuous historical monitoring period. The multi-period land use change survey vector set includes the boundaries of land use type patches at different survey times and the corresponding land use type codes. Based on the land use type codes, extract the forest land to non-forest land patch change vector, grassland to non-grassland patch change vector, and wetland to non-wetland patch change vector. The forest-to-non-forest patch change vector, the grassland-to-non-grass patch change vector, and the wetland-to-non-wetland patch change vector are spatially overlaid with the carbon storage temporal differentiation map. The average loss magnitude of the carbon storage temporal change raster in the forest-to-non-forest patch change vector coverage area, the average loss magnitude of the carbon storage temporal change raster in the grassland-to-non-grass patch change vector coverage area, and the average loss magnitude of the carbon storage temporal change raster in the wetland-to-non-wetland patch change vector coverage area are calculated. The average loss magnitudes of the forest-to-non-forest patch change vector, the grassland-to-non-grass patch change vector, and the wetland-to-non-wetland patch change vector are then input into an ecotype conversion carbon loss sensitivity ranking algorithm. Based on the magnitude of the average loss magnitude, the sensitivity levels of each type of patch change vector are ranked to generate an ecotype conversion carbon loss sensitivity ranking list. Based on the ecological type conversion carbon loss sensitivity ranking, target ecological conversion carbon loss hotspot patch sets that are ranked within a set ranking range by average loss magnitude are selected from the forest-to-non-forest patch change vector, grassland-to-non-grass patch change vector, and wetland-to-non-wet patch change vector. The spatial boundary coordinate sequence of the target ecological conversion carbon loss hotspot patch sets is extracted, and the spatial boundary coordinate sequence is compared with the spatial unit set of the partition. The spatial boundary coordinate sequence of the target ecological conversion carbon loss hotspot patch sets is marked with the first type of historical carbon loss hotspot coordinates located within the boundary coordinates of the carbon sink maintenance zone, the second type of historical carbon loss hotspot coordinates located within the boundary coordinates of the carbon source control zone, and the third type of historical carbon loss hotspot coordinates located within the boundary coordinates of the carbon sink enhancement zone. For the first type of historical carbon loss hotspot coordinates, historical carbon loss ecological restoration strategy entries are extracted from a preset land low-carbon restoration strategy library. These historical carbon loss ecological restoration strategy entries include an ecological reverse succession restoration operation that reverses the current land use type code of the area covered by the first type of historical carbon loss hotspot coordinates to the original forest land code, original grassland code, or original wetland code. For the second type of historical carbon loss hotspot coordinates, historical carbon loss source control strategy entries are extracted from the preset land low-carbon restoration strategy library. These historical carbon loss source control strategy entries include a source control management operation that permanently prohibits development and construction in the area covered by the second type of historical carbon loss hotspot coordinates. For the third type of historical carbon loss hotspot coordinates, historical carbon loss potential activation strategy entries are extracted from the preset land low-carbon restoration strategy library. These historical carbon loss potential activation strategy entries include a carbon sequestration function reactivation operation that rapidly increases soil organic matter in the area covered by the third type of historical carbon loss hotspot coordinates. The coordinates of the first type of historical carbon loss hotspots are bound to the entry space of the historical carbon loss ecological restoration strategy to generate a historical loss restoration and repair sub-instruction; the coordinates of the second type of historical carbon loss hotspots are bound to the entry space of the historical carbon loss source containment strategy to generate a historical loss containment and repair sub-instruction; and the coordinates of the third type of historical carbon loss hotspots are bound to the entry space of the historical carbon loss potential activation strategy to generate a historical loss activation and repair sub-instruction. The historical loss restoration and repair sub-instruction, the historical loss blocking and repair sub-instruction, and the historical loss activation and repair sub-instruction are merged into the zonal differential land low-carbon restoration instruction set, which is used to guide differentiated low-carbon restoration intervention operations for carbon loss hotspots caused by ecological type transformation in historical periods.

Images