A multi-temporal land use change dynamic analysis method based on carbon emission accounting

Abstract

The application discloses a kind of multi-time sequence land use change dynamic analysis method based on carbon emission accounting, and the application relates to carbon emission accounting and land use spatial data processing technical field, construct the graph convolution network model based on space-time attention mechanism, graph structure data and the land use feature sequence of each grid are input into graph convolution network model, the aggregation weight of node on space-time dimension is dynamically adjusted by space-time attention mechanism and time attention mechanism, the deep space-time feature representation that multi-time sequence dynamic evolution characteristics are fused is output, the advantages of the present application are that: the present application breaks through the technical bottleneck of only comparing a few static time nodes in traditional method, by constructing graph structure data and fusing space-time attention mechanism and graph convolution network, the high nonlinearity and collateral effect of land use change in long-time sequence dimension and wide-area spatial topology are fully mined, the accuracy and robustness of evolution feature extraction are greatly improved.

Description

Technical Field

[0001] This invention relates to a multi-temporal dynamic analysis method for land use change based on carbon emission accounting, belonging to the field of land spatial planning. Background Technology

[0002] As the global climate change problem becomes increasingly severe, reducing carbon emissions and achieving carbon neutrality has become a core consensus of the international community. Land use change is one of the most direct manifestations of human activities affecting the natural ecosystems on the Earth's surface, and it is also a key factor causing dynamic fluctuations in regional and even global carbon emissions. Different types of land use, such as urban construction land, agricultural land, forests and wetlands, play completely different roles in the carbon cycle, serving as both carbon sources and carbon sinks. Therefore, accurately calculating the evolution characteristics of carbon emissions brought about by different land use changes is an important foundation for scientifically formulating low-carbon development policies and national spatial planning.

[0003] Currently, there are some existing technologies for land use and carbon budget accounting. For example, Chinese Patent Publication No. CN 118569455 B discloses a method and system for carbon budget accounting before and after land spatial planning. This method acquires relevant data before and after the implementation of land spatial planning, stores the data in an accounting platform and transmits it to a cloud platform, and performs spatiotemporal evolution analysis based on carbon source accounting models and carbon sink accounting models through spatial analysis algorithms. However, through in-depth analysis, it was found that this existing technology mainly relies on the data from two specific time nodes before and after the implementation of the plan for accounting and comparison. It fails to effectively handle the cumulative effect of carbon emissions and the high spatiotemporal heterogeneity caused by frequent changes in land use under multiple time series over a long period of time. At the same time, other existing land use spatial evolution analysis models, such as the CA-Markov model based on Markov chains and cellular automata, are mostly limited to the independent transition probability analysis of a single pixel state. They lack the ability to explore the topological relationship of spatial neighborhoods and the deep dynamic characteristics of long-term evolution, and have the technical defect of being unable to capture complex nonlinear spatiotemporal coupling laws. Summary of the Invention

[0004] The purpose of this invention is to provide a multi-temporal dynamic analysis method for land use change based on carbon emission accounting, so as to solve the above-mentioned problems.

[0005] The present invention adopts the following technical solution:

[0006] A multi-temporal-series land use change dynamic analysis method based on carbon emission accounting, characterized by the following steps:

[0007] Step 1: Acquire multi-temporal remote sensing image data covering the target area and corresponding basic statistical data. Perform radiometric calibration, atmospheric correction and geometric registration preprocessing operations on the multi-temporal remote sensing image data in sequence to obtain a standardized reflectance image sequence.

[0008] Step 2: Spatial gridding is performed on the standardized reflectance image sequence to divide the target area into multiple spatial grids. Each spatial grid is used as a graph node, and an adjacency matrix is ​​constructed based on the geographic spatial distance between grids and the Gaussian decay function to establish graph structure data that represents the spatial topological relationship of land use.

[0009] Step 3: Construct a graph convolutional network model based on spatiotemporal attention mechanism. Input the graph structure data and the land use feature sequences of each spatial grid into the graph convolutional network model. Dynamically adjust the aggregation weights of nodes in the spatiotemporal dimension through spatial attention mechanism and temporal attention mechanism, and output a deep spatiotemporal feature representation that integrates multiple temporal dynamic evolution features.

[0010] Step 4: Establish a multi-dimensional dynamic carbon budget accounting model. Combine the deep spatiotemporal feature representations of each spatial grid with the corresponding land use feature sequences to calculate the agricultural carbon source emissions, building carbon source emissions and forest carbon sink absorption at each time series node, and obtain the multi-time series carbon emission evolution index sequence of the target area.

[0011] Step 5: Based on the deep spatiotemporal feature representation and the multi-time series carbon emission evolution index sequence, a time series state transition prediction mechanism for the joint evolution of land use and carbon emissions is constructed using a hidden Markov model. The land use distribution and carbon emission change trends of the target area at future time nodes are calculated and predicted, and the analysis and prediction results are visualized and output.

[0012] Furthermore, the multi-temporal land use change dynamic analysis method based on carbon emission accounting of the present invention also has the following feature: In step one, an adjacency matrix is ​​constructed based on the geospatial distance between grids and the Gaussian decay function to establish graph structure data representing the spatial topological relationship of land use, including the following steps:

[0013] Define the spatial land use topology map structure as follows: ,in The set of all spatial grid nodes within the target area. This is the set of edges representing the topological connections between different grid nodes;

[0014] For any two grid nodes and Calculate its Euclidean geospatial distance:

[0015] ;

[0016] in , Grid nodes The central x-coordinate and y-coordinate values ​​in the geographic coordinate system , Grid nodes The central x-coordinate and y-coordinate values ​​in the geographic coordinate system;

[0017] Calculate the adjacency weights between nodes based on the Gaussian decay function:

[0018] ;

[0019] in, The standard deviation is used as a control parameter to control the magnitude of spatial correlation decay.

[0020] Construct an adjacency matrix based on the adjacency weights. Define the degree matrix Given a diagonal matrix, the elements on its diagonal are the sum of all elements in the corresponding row of the adjacent matrix. Then, the normalized Laplace matrix is ​​calculated:

[0021] ;

[0022] in, It is a unit diagonal matrix whose dimensions match the number of grid nodes.

[0023] Furthermore, the multi-temporal land use change dynamic analysis method based on carbon emission accounting of the present invention also has the following feature: In step three, the spatial attention mechanism is implemented as follows:

[0024] For the input node feature matrix formed by multiple time series Calculate the spatial attention score matrix:

[0025] ;

[0026] in, This is the weight transformation parameter matrix for the spatial attention layer. These are the bias adjustment parameters used in spatial attention calculations. For spatial attention, reduce the dimensionality of the mapping vector. This is the transpose of the feature matrix of the input nodes. This is the spatial attention matrix between nodes.

[0027] Furthermore, the multi-temporal land use change dynamic analysis method based on carbon emission accounting of the present invention also has the following feature: In step three, the implementation of the time attention mechanism is as follows:

[0028] For the input node feature matrix formed by multiple time series Calculate the time attention matrix:

[0029] ;

[0030] in, Transform the weight matrix for time attention features. Here is the temporal attention bias parameter matrix. The output of temporal attention is mapped to a weight vector. This is the dependency weight matrix between multiple time steps.

[0031] Furthermore, the multi-temporal land use change dynamic analysis method based on carbon emission accounting of the present invention also has the following feature: in step three, the graph convolutional network model's... The layer forward propagation calculation method is as follows:

[0032] ;

[0033] in, Represents the Hadamard product of matrices. It is an adjacency matrix. The spatial attention score matrix is... The improvement matrix after adding self-loop connections, For the first The spatiotemporal features of the input to layer graph convolution, For the first The model parameter matrix optimized by the backpropagation algorithm for each layer. It is a non-linear activation function;

[0034] Initial input From the input node feature matrix With the time attention matrix The Hadama product is determined. This is a high-dimensional feature representation of the fused spatial structure and temporal attributes output after computation by this layer.

[0035] Furthermore, the multi-temporal land use change dynamic analysis method based on carbon emission accounting of the present invention also has the following feature: In step four, the multi-dimensional dynamic carbon budget accounting model specifically includes:

[0036] Agricultural carbon emissions The calculation method is as follows: the actual coverage area of ​​agricultural land in the target grid is... Combined coefficient of agricultural carbon emissions per unit area based on historical statistical fitting Multiplication, that is ;

[0037] Building carbon emissions The calculation method is as follows: the total energy consumption of the building area is calculated by inverting the nighttime light index and the density of points of interest within the statistical period. Energy consumption carbon emission factor constant per unit standard coal equivalent within the region Multiplication, that is ;

[0038] Forest carbon sequestration The calculation method is as follows: the forest patch coverage area extracted through remote sensing is used as the basis for the calculation. Compared with the annual standard carbon sequestration rate per unit area of ​​this forest type Multiplication, that is ;

[0039] The carbon source emissions and carbon sink absorption under each grid are fused and superimposed to obtain a dynamic rolling carbon emission index sequence under continuous time-series nodes.

[0040] Furthermore, the multi-temporal land use change dynamic analysis method based on carbon emission accounting of the present invention also has the following feature: In step five, a temporal state transition prediction mechanism for the joint evolution of land use and carbon emissions is constructed using a hidden Markov model, including the following steps:

[0041] The quintuples for constructing a Hidden Markov Model include:

[0042] Observation sequence composed of remote sensing features and carbon emission indicators Hidden state sequences characterizing the intrinsic state of true land cover transformation and carbon emission trends State transition probability matrix Observation probability matrix and the initial state probability distribution ;

[0043] The forward-backward algorithm is used to iteratively optimize the model parameters, and the forward probability variable is defined. Indicates in Historical sequences are observed at all times and the current hidden state is The joint probability, defining the backward probability variable. Indicates in The state at time is The conditional probability of future observed sequences occurring under given conditions is iteratively updated using an expectation-maximization approach. , and Until the model converges;

[0044] The Viterbi algorithm is used to calculate the optimal hidden state path and infer the direction of land use type transformation and carbon emission evolution of each grid in the target area at future time nodes.

[0045] Furthermore, the multi-temporal land use change dynamic analysis method based on carbon emission accounting of the present invention also has the following feature: after obtaining the standardized reflectance image sequence in step one and before constructing the map structure data in step two, it also includes an outlier detection and cleaning step based on isolated forests.

[0046] Using the Isolation Forest model to analyze each data point in the temporal spatial sample point set Calculate the anomaly score:

[0047] ;

[0048] in, For data points The average path height required to successfully isolate a single isolated tree. The expected value of the path height of this data point calculated for all isolated trees. The normalization constant for the average path length of a binary search tree with n sample points;

[0049] When the abnormal score When the value exceeds the preset anomaly detection threshold, the feature data of the corresponding grid is marked as abnormal noise and removed, and spatial neighborhood mean interpolation is performed for repair.

[0050] Furthermore, the multi-temporal land use change dynamic analysis method based on carbon emission accounting of the present invention also has the following feature: in step four, before performing the multi-dimensional dynamic carbon budget accounting, a time-series smoothing step based on Kalman filtering is also included.

[0051] Based on the system prior error covariance matrix Observation transformation matrix and measurement noise covariance matrix Calculate the Kalman gain parameters:

[0052] ;

[0053] Based on the Kalman gain parameter, the actual measured value is... By incorporating the state equations, the a posteriori update and correction of the system's internal state estimates are completed:

[0054] ;

[0055] in, The predicted prior state value is calculated based on the system state evolution equation at the current moment. This is the posterior precise state estimate after fusion and correction of observational data;

[0056] By iterating the Kalman filter update process step by step, the temporal land use parameters and carbon source carbon summation containing uncertain perturbation components are smoothed.

[0057] Furthermore, the multi-temporal land use change dynamic analysis method based on carbon emission accounting of the present invention also has the following feature: In step two, the division of multiple spatial grids adopts a multi-temporal analysis execution engine based on a distributed parallel computing architecture, including the following steps:

[0058] A spatial grid dynamic cutting mechanism based on quadtree depth is adopted to allocate continuous geographic space in blocks. The termination threshold condition for the spatial recursive segmentation depth is set as the total number of pixels under the current data block is less than the memory carrying capacity limit of a single computing node in the cluster.

[0059] When generating boundaries by cutting the spatial graph, an overlapping halo expansion technique is implemented, that is, each computational subgraph expands outward from its absolute spatial edge with a layer of thickness. Graph node information.

[0060] Beneficial effects of the invention:

[0061] (1) This invention breaks through the technical bottleneck of traditional methods that only compare a few static time nodes. By constructing graph structure data and integrating spatiotemporal attention mechanism and graph convolutional network, it fully explores the high nonlinearity and cascading effects of land use change in long time series dimension and wide-area spatial topology, and greatly improves the accuracy and robustness of evolution feature extraction.

[0062] (2) This invention realizes the deep spatiotemporal coupling of land use dynamic evolution and carbon budget accounting. By constructing a joint driving mechanism of multidimensional dynamic carbon budget accounting model and hidden Markov prediction model, it not only realizes the quantitative analysis of historical carbon emission evolution law, but also has the ability to accurately predict future land use and carbon emission evolution trends, providing strong and reliable technical support for intelligent and low-carbon modern land space refined governance.

[0063] Based on a refined big data and spatial information processing technology framework, multi-temporal remote sensing images are transformed into topological graph structures. By combining spatiotemporal attention mechanisms and graph convolutional networks (STA-GCN), land use evolution trajectories are dynamically extracted from global topology and deep nonlinear spatiotemporal dimensions.

[0064] It is deeply integrated with multidimensional carbon emission dynamic accounting models and hidden Markov models to achieve high-precision mining and prediction of land use and carbon emission evolution in multiple time series, providing intelligent and precise decision support for low-carbon management of national land space planning. Attached Figure Description

[0065] Figure 1This is a schematic diagram of the overall process logic of the multi-temporal land use change dynamic analysis method based on carbon emission accounting in an embodiment of the present invention;

[0066] Figure 2 This is a flowchart illustrating the execution of multi-temporal joint prediction using a hidden Markov model in an embodiment of the present invention. Detailed Implementation

[0067] To enable those skilled in the art to fully understand the content of this invention and to realize its various functions without creative effort, the following is a comprehensive and clear description of the specific implementation methods of this invention, in conjunction with specific embodiments and detailed formula principles.

[0068] Throughout the entire description of this invention, all formula variables and their corresponding representative symbols are kept unique, and the same symbol is not used to refer to different parameters, so as to ensure the rigor and unambiguity of the technical solution.

[0069] Please see the appendix Figure 1 and attached Figure 2 This invention provides a dynamic analysis method for multi-temporal land use change based on carbon emission accounting.

[0070] Example 1

[0071] This embodiment details the complete processing flow of a multi-temporal-series land use change dynamic analysis method based on carbon emission accounting. For example... Figure 1 As shown:

[0072] I. Multi-source data acquisition and preprocessing module

[0073] For the target area, long-term continuous multi-time series remote sensing image datasets are obtained from multi-source satellite remote sensing platforms, such as Landsat and Sentinel series, as the basic input source. In order to eliminate the interference of the characteristics of the satellite sensors themselves and external objective factors such as sunlight and atmospheric scattering on the true reflectivity of the surface, a strict preprocessing procedure must be performed on the original remote sensing images.

[0074] First, the original remote sensing images To perform radiometric calibration calculations, the dimensionless grayscale values ​​recorded by the sensor are converted into apparent radiance. The radiometric calibration calculation formula is as follows:

[0075] ;

[0076] in, Represents the apparent radiance image matrix after radiometric calibration;

[0077] This represents the radiation calibration gain coefficient of the sensor;

[0078] This represents a constant indicating the radiometric calibration offset of the sensor.

[0079] Next, the apparent radiance image matrix was analyzed using a radiative transfer model. Atmospheric correction is performed to eliminate the path effects of atmospheric scattering and absorption, thereby obtaining a true surface reflectance image. .

[0080] The conversion formula is expressed as follows:

[0081] ;

[0082] in, This represents the amount of atmospheric path radiation interference.

[0083] This is the total atmospheric transmittance parameter in the direction of remote sensing observation;

[0084] This represents the average solar spectral irradiance in the upper atmosphere;

[0085] The solar zenith angle at the moment of imaging.

[0086] After radiometric calibration and atmospheric correction, geometric registration is performed using control point matching technology to eliminate image distortion and finally output a standardized multi-temporal reflectance image sequence.

[0087] II. Multi-temporal land use topology map structure reconstruction module

[0088] In order to characterize the changes in land use in the target area from the perspective of overall spatial connectivity, this invention abandons the traditional approach of processing individual pixels in isolation and proposes to perform regular grid segmentation on the preprocessed reflectance image, for example, dividing it into spatial units of 100 meters × 100 meters, and constructing a spatial topology map based on this.

[0089] Define the spatial land use topology map structure as follows: ,in The set of all spatial grid nodes within the target area, denoted as . ;

[0090] The set of edges used to characterize the topological connections between different grid nodes.

[0091] To accurately reflect the correlation strength of different grids in spatial location, an adjacency matrix with a Gaussian decay function is used to construct a graph structure. .

[0092] For any two grid nodes in the graph and Its corresponding adjacency weight The calculation formula is:

[0093] ;

[0094] ;

[0095] In the formula, , and , Representing grid nodes respectively and grid nodes The center coordinates in the geographic coordinate system;

[0096] For the calculated nodes and Euclidean geospatial distance;

[0097] The standard deviation is a control parameter used to control the magnitude of spatial correlation decay.

[0098] according to Construct the complete adjacency matrix Then, define the degree matrix. It is a diagonal matrix, where the elements on the diagonal are the sum of all elements in the adjacent matrix of the corresponding row. The normalized Laplace matrix can then be calculated from this. :

[0099] ;

[0100] in, The dimension is The unit diagonal matrix, the Laplace matrix It contains the Laplacian smoothness between the entire spatial grid and is the core computational factor for subsequent graph convolution operations.

[0101] III. Feature Extraction Module Based on Spatiotemporal Attention Graph Convolutional Network STA-GCN

[0102] For the input node feature matrix formed by multiple time series The matrix contains the spectral, index and other feature dimensions of each grid in multiple consecutive time periods. This module first executes a spatial attention mechanism to mine global long-distance dependencies.

[0103] In the spatial attention layer, a spatial attention score matrix is ​​defined. as follows:

[0104] ;

[0105] in, This represents the weight transformation parameter matrix trained on the spatial attention layer;

[0106] This represents the bias adjustment parameter in spatial attention calculation;

[0107] Dimensionality reduction mapping vector for spatial attention;

[0108] After mapping and activation, the mutual attention matrix between nodes is obtained. .

[0109] Meanwhile, a time attention mechanism is used to dynamically measure the differences in the contribution of data from different time series stages to the current evolution trend:

[0110] ;

[0111] in, The weight matrix represents the transformation of temporal attention features;

[0112] Represents the temporal attention bias parameter matrix;

[0113] The output of temporal attention is mapped to a weight vector;

[0114] The matrix describes the dependency weights between multiple time steps.

[0115] The aforementioned spatial and temporal attention weights are applied to the graph structure to perform spectral domain graph convolution.

[0116] For the The forward propagation calculation formula for a layered graph convolutional neural network is as follows:

[0117] ;

[0118] in, Represents the Hadamard product (element-by-element product) of matrices.

[0119] For the first The spatiotemporal features of the input to the layer graph convolution (initial time) );

[0120] The improvement matrix after adding self-loop connections;

[0121] For the first The layer needs to be optimized using the backpropagation algorithm to obtain the model parameter matrix;

[0122] It is a non-linear activation function (such as the ReLU function);

[0123] This refers to the high-dimensional feature representation output after calculation by this layer, which contains spatial structure and temporal fusion attributes.

[0124] After multi-layer stacking computation, the network finally outputs a deep comprehensive feature representation that describes the evolution trajectory of land use changes in time and space.

[0125] IV. Multidimensional Dynamic Carbon Budget Accounting Module

[0126] Based on the accurate acquisition of land use classification information for each spatial grid, carbon emissions and carbon sinks are measured for each grid unit to calculate the total net carbon emissions for the entire target area. It is necessary to comprehensively consider agricultural carbon source emissions. Carbon emissions from buildings And deduct forest land carbon sequestration amount .

[0127] 1. Agricultural carbon source emission accounting:

[0128] Agricultural land use is largely derived from fertilizer application and agricultural machinery consumption; the calculation formula is as follows:

[0129] ;

[0130] in, This represents the total actual agricultural land coverage area determined within the target grid.

[0131] This represents the comprehensive coefficient of agricultural carbon emissions per unit area, based on statistical fitting from historical yearbooks.

[0132] 2. Carbon emission accounting for buildings and transportation:

[0133] Carbon emissions from building sites are primarily based on electricity consumption and fossil fuel consumption, calculated using the following formula:

[0134] ;

[0135] in, The total standard energy consumption of the building area during the statistical period (obtained by inversion of nighttime light index and point of interest density).

[0136] This represents the energy consumption carbon emission factor constant per unit standard coal equivalent within a region.

[0137] 3. Carbon sequestration accounting for forest land and ecological land:

[0138] Forests possess remarkable photosynthetic carbon sequestration capabilities, calculated using the following formula:

[0139] ;

[0140] in, This represents the forest patch coverage area extracted through remote sensing;

[0141] The annual standard carbon sequestration rate per unit area for this forest type (evergreen broadleaf, deciduous conifer, etc.).

[0142] Finally, the carbon source and carbon sink indicators under this grid are fused, superimposed, and accumulated to obtain a dynamic rolling carbon emission indicator sequence under continuous time-series nodes.

[0143] V. Hidden Markov Multi-Time Evolution Prediction Module

[0144] By combining the spatiotemporal features extracted from deep graph convolution with the continuously updated carbon budget sequence, a Hidden Markov Model (HMM) is used to predict future time steps. Figure 2 As shown:

[0145] Constructing HMM quintuples: observation sequence Remote sensing features and carbon emission indicators; hidden state sequences The intrinsic state of actual land cover transformation and carbon emission trends; state transition probability matrix. Observation probability matrix and the initial state probability distribution .

[0146] The model training process uses the forward-backward algorithm, namely the Baum-Welch algorithm, to iteratively optimize the model parameters.

[0147] Define forward probability variables , indicating in Historical sequences are observed at all times and the current hidden state is The joint probability;

[0148] Define backward probability variables , indicating in The state at time is Under certain conditions, the conditional probability of future observed sequences occurring.

[0149] Based on the above definition, time is calculated. In state probability distribution and the joint distribution of state pairs at adjacent time points :

[0150] ;

[0151] ;

[0152] Update using the iterative expectation maximization (EM) method , and Once the model converges, the Viterbi algorithm is applied to calculate the optimal hidden state path. This allows for the inference and reconstruction of the land use type changes and corresponding expected carbon emission peaks and evolution directions for each grid in the target area at the next moment. The dynamic simulation results serve as the core basis for comprehensive land space management and emission reduction plan formulation.

[0153] Example 2

[0154] To address the pain points of real-world business scenarios where multi-source heterogeneous data often contains significant observation noise and exhibits high uncertainty, this embodiment expands upon the first embodiment with technical extensions and data flow optimization, with particular emphasis on anomaly cleaning, time-series filtering and smoothing, and computational strategies under a distributed parallel architecture.

[0155] In the process of integrating and connecting multi-source basic data with statistical yearbooks, local data extreme values ​​are often missing or incorrectly entered. This embodiment introduces an unsupervised outlier detection and cleaning mechanism based on Isolation Forest.

[0156] For a temporal spatial sample point set, Isolation Forest constructs multiple isolated trees by randomly selecting data features and randomly determining cutoff values, defining each data point. Outlier scores in isolated forest models The mathematical expression is:

[0157] ;

[0158] in, Representative data points The average path height experienced when successfully isolated in a single tree;

[0159] The expected value of the path height of this data point is calculated for all generated isolated trees;

[0160] The normalization constant for the average path length of a binary search tree consisting of a set of n sample points is given if and only if the calculated path length is given by the normalization constant. When the value significantly approaches 1.0, a certain feature data of the grid is marked as strong anomalous noise and removed or repaired by spatial neighborhood mean interpolation.

[0161] After anomaly cleaning, this embodiment constructs a Kalman Filter sequence smoothing subsystem to smooth the time-series carbon emissions and characteristic index curve fluctuations caused by short-term sensor failures or meteorological interference. This subsystem is used to perform optimal state estimation smoothing of time-series land use parameters and total carbon source carbon containing uncertain disturbance components.

[0162] The state update process of Kalman filtering involves two steps: prior estimation and posterior update. First, the Kalman gain parameter, which measures the prediction error at the current time with respect to the subsequent estimation confidence, is calculated. :

[0163] ;

[0164] In the formula, This is the system prior error covariance matrix derived from the previous time step;

[0165] This is the observation transformation matrix that maps the real internal state space to the external measurement space;

[0166] The covariance matrix represents the inherent measurement noise of the measurement system;

[0167] This is the transpose of the observation transformation matrix.

[0168] In obtaining Kalman gain Then, the actual measurement values ​​acquired by the sensing instrument are... By incorporating state equations, the internal state estimates of a smooth system can be obtained. and the posterior error covariance matrix Synchronous update fixes:

[0169] ;

[0170] ;

[0171] in, This represents the system state evolution equation at the current moment, including the system evolution matrix. The calculated prior state prediction value;

[0172] This is the posterior precise state estimate after smoothing and fusion correction of the observed data. Meanwhile, the iterative preparation formula for the prior prediction variance of the next time step is: ,in The external environmental process disturbance noise covariance.

[0173] This Kalman filter sequence processing module effectively offsets the carbon emission calculation jump errors caused by factors such as remote sensing cloud cover.

[0174] For massive images of ultra-large-scale land space and grid processing tasks involving hundreds of millions of pixels, ordinary computing nodes cannot meet the timeliness requirements of the calculation. Therefore, this embodiment innovatively designs a multi-time-series analysis execution engine based on a distributed parallel computing architecture, such as the Hadoop / Spark distributed cluster ecosystem, at the hardware and software deployment level.

[0175] In the data input stage, to ensure that the spatiotemporal attention-based graph convolutional network can effectively process spatial graph data without boundary effects, this invention adopts a spatial grid dynamic segmentation mechanism based on quadtree depth to perform block-like allocation of continuous geographic space and set a spatial recursive segmentation depth. The termination threshold condition is that the total number of pixels in the current data block is less than the maximum memory capacity of a single compute node in the cluster. When generating boundaries by cutting the spatial graph, the HaloOverlap expansion technique is implemented, that is, each compute subgraph expands outward from its absolute spatial edge with a layer of thickness. Graph node information to ensure that boundary nodes are processed when executing the Laplace matrix. No truncation information loss occurs during convolutional aggregation calculations.

[0176] After all distributed computing nodes have completed the STA-GCN feature extraction and multidimensional carbon accounting tasks, the redundant states in overlapping areas are removed using the Reduce aggregation mechanism, and the overall global prediction map of the multi-temporal evolution of land space carbon balance is spliced ​​out.

[0177] This parallel design fully leverages the advantages of big data framework clusters, achieving industrial-grade acceleration and ultra-large-scale adaptability for carbon emission evolution projection while meeting rigorous spatiotemporal dependency computation rules.

[0178] The above-described embodiments are merely for illustrating the specific technical implementation methods adopted by the present invention. In actual implementation, the core functional modules, rigorous formula variable derivation, and their derived big data distributed optimization improvements are described. It should be noted that for those skilled in the art, any minor improvements, equivalent parameter substitutions, or hardware adaptation adjustments made without departing from the fundamental creative concept and rigorous logical algorithm principle of the present invention should be included and fall within the scope of protection that the present invention intends to seek through legal claims.

Claims

1. A multi-temporal-series land use change dynamic analysis method based on carbon emission accounting, characterized in that, Includes the following steps: Step 1: Acquire multi-temporal remote sensing image data covering the target area and corresponding basic statistical data. Perform radiometric calibration, atmospheric correction and geometric registration preprocessing operations on the multi-temporal remote sensing image data in sequence to obtain a standardized reflectance image sequence. Step 2: Spatial gridding is performed on the standardized reflectance image sequence to divide the target area into multiple spatial grids. Each spatial grid is used as a graph node, and an adjacency matrix is ​​constructed based on the geographic spatial distance between grids and the Gaussian decay function to establish graph structure data that represents the spatial topological relationship of land use. Step 3: Construct a graph convolutional network model based on spatiotemporal attention mechanism. Input the graph structure data and the land use feature sequences of each spatial grid into the graph convolutional network model. Dynamically adjust the aggregation weights of nodes in the spatiotemporal dimension through spatial attention mechanism and temporal attention mechanism, and output a deep spatiotemporal feature representation that integrates multiple temporal dynamic evolution features. Step 4: Establish a multi-dimensional dynamic carbon budget accounting model. Combine the deep spatiotemporal feature representations of each spatial grid with the corresponding land use feature sequences to calculate the agricultural carbon source emissions, building carbon source emissions and forest carbon sink absorption at each time series node, and obtain the multi-time series carbon emission evolution index sequence of the target area. Step 5: Based on the deep spatiotemporal feature representation and the multi-time series carbon emission evolution index sequence, a time series state transition prediction mechanism for the joint evolution of land use and carbon emissions is constructed using a hidden Markov model. The land use distribution and carbon emission change trends of the target area at future time nodes are calculated and predicted, and the analysis and prediction results are visualized and output.

2. The multi-temporal land use change dynamic analysis method based on carbon emission accounting according to claim 1, characterized in that, In step one, an adjacency matrix is ​​constructed based on the geospatial distance between grids and the Gaussian decay function to establish graph-structured data representing the spatial topological relationships of land use. This includes the following steps: Define the spatial land use topology map structure as follows: ,in The set of all spatial grid nodes within the target area. This is the set of edges representing the topological connections between different grid nodes; For any two grid nodes and Calculate its Euclidean geospatial distance: ； in , Grid nodes The central x-coordinate and y-coordinate values ​​in the geographic coordinate system , Grid nodes The central x-coordinate and y-coordinate values ​​in the geographic coordinate system; Calculate the adjacency weights between nodes based on the Gaussian decay function: ； in, The standard deviation is used as a control parameter to control the magnitude of spatial correlation decay. Construct an adjacency matrix based on the adjacency weights. Define the degree matrix Given a diagonal matrix, the elements on its diagonal are the sum of all elements in the corresponding row of the adjacent matrix. Then, the normalized Laplace matrix is ​​calculated: ； in, It is a unit diagonal matrix whose dimensions match the number of grid nodes.

3. The multi-temporal land use change dynamic analysis method based on carbon emission accounting according to claim 1, characterized in that, In step three, the spatial attention mechanism is implemented as follows: For the input node feature matrix formed by multiple time series Calculate the spatial attention score matrix: ； in, This is the weight transformation parameter matrix for the spatial attention layer. These are the bias adjustment parameters used in spatial attention calculations. To reduce the dimensionality of spatial attention mapping vectors, This is the transpose of the feature matrix of the input nodes. This is the spatial attention matrix between nodes.

4. The multi-temporal land use change dynamic analysis method based on carbon emission accounting according to claim 1, characterized in that, In step three, the temporal attention mechanism is implemented as follows: For the input node feature matrix formed by multiple time series Calculate the time attention matrix: ； in, Transform the weight matrix for time attention features. Here is the temporal attention bias parameter matrix. The output mapping weight vector for time attention. This is the dependency weight matrix between multiple time steps.

5. The multi-temporal land use change dynamic analysis method based on carbon emission accounting according to claim 1, characterized in that, In step three, the graph convolutional network model's first... The layer forward propagation calculation method is as follows: ； in, Represents the Hadamard product of matrices. It is an adjacency matrix. The spatial attention score matrix is... The improvement matrix after adding self-loop connections, For the first The spatiotemporal features of the input to layer graph convolution, For the first The model parameter matrix optimized by the backpropagation algorithm for each layer. It is a non-linear activation function; Initial input From the input node feature matrix With the time attention matrix The Hadama product is determined. This is a high-dimensional feature representation of the fused spatial structure and temporal attributes output after computation by this layer.

6. The multi-temporal land use change dynamic analysis method based on carbon emission accounting according to claim 5, characterized in that, In step four, the multidimensional dynamic carbon budget accounting model specifically includes: Agricultural carbon emissions The calculation method is as follows: the actual coverage area of ​​agricultural land in the target grid is... Combined coefficient of agricultural carbon emissions per unit area based on historical statistical fitting Multiplication, that is ; Building carbon emissions The calculation method is as follows: the total energy consumption of the building area is calculated by inverting the nighttime light index and the density of points of interest within the statistical period. Energy consumption carbon emission factor constant per unit standard coal equivalent within the region Multiplication, that is ; Forest carbon sequestration The calculation method is as follows: the forest patch coverage area extracted through remote sensing is used as the basis for the calculation. Compared with the annual standard carbon sequestration rate per unit area of ​​this forest type Multiplication, that is ; The carbon source emissions and carbon sink absorption under each grid are fused and superimposed to obtain a dynamic rolling carbon emission index sequence under continuous time-series nodes.

7. The multi-temporal land use change dynamic analysis method based on carbon emission accounting according to claim 6, characterized in that, Step five involves constructing a time-series state transition prediction mechanism for the joint evolution of land use and carbon emissions using a hidden Markov model, including the following steps: The quintuples for constructing a Hidden Markov Model include: Observation sequence composed of remote sensing features and carbon emission indicators Hidden state sequences characterizing the intrinsic state of true land cover transformation and carbon emission trends State transition probability matrix Observation probability matrix and the initial state probability distribution ; The forward-backward algorithm is used to iteratively optimize the model parameters, and the forward probability variable is defined. Indicates in Historical sequences are observed at all times and the current hidden state is The joint probability, defining the backward probability variable. Indicates in The state at time is The conditional probability of future observed sequences occurring under given conditions is iteratively updated using an expectation-maximization approach. , and Until the model converges; The Viterbi algorithm is used to calculate the optimal hidden state path and infer the direction of land use type transformation and carbon emission evolution of each grid in the target area at future time nodes.

8. The multi-temporal land use change dynamic analysis method based on carbon emission accounting according to claim 1, characterized in that, After obtaining the standardized reflectance image sequence in step one and before constructing the graph structure data in step two, an outlier detection and cleaning step based on isolated forest is also included: Using the Isolation Forest model to analyze each data point in the temporal spatial sample point set Calculate the anomaly score: ； in, For data points The average path height required to successfully isolate a single isolated tree. The expected value of the path height of this data point calculated for all isolated trees. The normalization constant for the average path length of a binary search tree with n sample points; When the abnormal score When the anomaly threshold is exceeded, the feature data of the corresponding grid is marked as abnormal noise and removed, and spatial neighborhood mean interpolation is performed for repair.

9. The multi-temporal land use change dynamic analysis method based on carbon emission accounting according to claim 1, characterized in that, Step four includes a time-series smoothing step based on Kalman filtering, prior to the multidimensional dynamic carbon budget calculation. Based on the system prior error covariance matrix Observation transformation matrix and measurement noise covariance matrix Calculate the Kalman gain parameters: ； Based on the Kalman gain parameter, the actual measured value is... By incorporating the state equations, the a posteriori update and correction of the system's internal state estimates are completed: ； in, The predicted prior state value is calculated based on the system state evolution equation at the current moment. This is the posterior precise state estimate after observation data fusion and correction; By iterating the Kalman filter update process step by step, the temporal land use parameters and carbon source carbon summation containing uncertain perturbation components are smoothed.

10. The multi-temporal land use change dynamic analysis method based on carbon emission accounting according to claim 1, characterized in that, In step two, the division of multiple spatial grids employs a multi-temporal analysis execution engine based on a distributed parallel computing architecture, including the following steps: A spatial grid dynamic cutting mechanism based on quadtree depth is adopted to allocate continuous geographic space in blocks. The termination threshold condition for the spatial recursive segmentation depth is set as the total number of pixels under the current data block is less than the memory carrying capacity limit of a single computing node in the cluster. When generating boundaries by cutting the spatial graph, an overlapping halo expansion technique is implemented, that is, each computational subgraph expands outward from its absolute spatial edge with a layer of thickness. Graph node information.

Images