A natural resource data analysis method and system based on artificial intelligence

Abstract

The application discloses a kind of natural resource data analysis method and system based on artificial intelligence.The method learns the causal graph structure between environmental factor and land class by time series-cause fusion graph neural network, determines the environmental response term coefficient in differential equation using causal graph, establishes continuous state vector evolution model to describe the gradual transition process of land class state;Introduce spectral constraint, time series constraint and spatial constraint to carry out multi-dimensional physical constraint optimization to state vector;Utilize causal path to enhance state representation training strategy network to output optimal management strategy;Change time is backtracked from current state by Bayes inversion and driving factor is traced back;Intelligent analysis report is generated.The application fuses time series law and causal mechanism to improve identification accuracy, expresses continuous evolution process of land class, introduces multi-dimensional physical constraint to ensure result rationality, realizes causal mechanism identification and intervention strategy optimization, and improves change time tracing accuracy.

Description

Technical Field

[0001] This invention belongs to the field of natural resource management technology, specifically relating to a natural resource data analysis method and system based on artificial intelligence. Background Technology

[0002] Natural resource data analysis is a crucial foundation for land spatial planning, resource management, and ecological protection. With the development of monitoring technologies, utilizing multi-source data for dynamic monitoring and analysis of natural resources has become a common technique, primarily encompassing key aspects such as land cover identification, state evolution, risk assessment, and change tracing. However, existing technologies face the following challenges in practical applications: First, existing land use identification methods mainly rely on statistical correlation to establish a mapping relationship between observed features and land use types, without fully considering the temporal variation patterns of land use types. The same land use type may show significant differences at different times, resulting in the identification accuracy being greatly affected by time factors, making it difficult to accurately distinguish land use types with similar features but different variation patterns.

[0003] Secondly, existing land use change detection methods employ discrete time-point comparative analysis, which can only identify the final result of land use change and cannot depict the continuous evolution process of land use. The transformation of natural resource land use is usually a gradual process, and traditional methods are unable to express the intermediate transitional states, making it difficult to detect changes in the early stages in a timely manner.

[0004] Third, existing data analysis methods are not adaptable enough when dealing with mixed states. They lack comprehensive constraints from multi-dimensional physical laws such as spectral, temporal, and spatial dimensions, which may lead to analysis results that do not conform to the actual situation and make it difficult to meet the demand for high-precision data in refined management.

[0005] Fourth, existing change analysis methods are mostly based on correlation analysis, lack the ability to identify the causal relationship of driving factors, make it difficult to accurately assess the actual effect of different management measures, and lack the ability to predict future change trends, thus failing to support proactive risk management.

[0006] Fifth, existing technologies are insufficient in tracing changes over time, relying on direct comparison of multi-temporal observation data. The detection accuracy is limited by the observation time interval, making it difficult to accurately deduce the historical evolution process from the current state.

[0007] Therefore, there is a need for a natural resource data analysis method and system that can combine temporal patterns and causal mechanisms, express gradual transition states, introduce physical constraints, quantify the effects of management measures, and improve the accuracy of time detection. Summary of the Invention

[0008] In view of the shortcomings of existing technologies, this invention aims to solve the technical problems of insufficient utilization of temporal patterns, lack of expression of continuous evolution of land types, insufficient introduction of physical constraints, weak ability to identify causal mechanisms, and low accuracy of tracing change time.

[0009] In a first aspect, the present invention provides a natural resource data analysis method based on artificial intelligence, comprising the following steps: S1. Multi-source data acquisition and temporal feature extraction: Acquire multi-temporal remote sensing image data, environmental auxiliary data, historical land cover data and management data of the study area, preprocess the data and extract temporal feature variables; S2. Construction of Temporal-Causal Fusion Graph Neural Network: Construct a temporal-causal fusion graph neural network, which includes a temporal encoder, a causal graph structure learner, and a physical constraint layer. The temporal encoder encodes temporal features to obtain temporal representation vectors. The causal graph structure learner learns the causal graph structure between environmental factors, temporal feature parameters, and land use labels. The physical constraint layer constrains the causal graph structure based on agronomic and ecological mechanisms. By jointly optimizing prediction accuracy, causal graph topological validity, and physical rationality, the causal graph structure and posterior probability of land use are output. S3. Continuous State Vector Evolution Modeling: Establish a continuous state vector evolution model, using continuous vectors to represent the land use state intensity of spatial units, describing the evolution process of the state vector based on differential equations, the coefficients of the environmental response term of the differential equations are determined according to the cause-effect graph structure, and output the time evolution trajectory of the state vector. S4. Multidimensional physical constraint optimization: The state vector is optimized by introducing spectral constraints, temporal constraints and spatial constraints. The spectral constraints ensure that the state vector is consistent with the observed spectrum, the temporal constraints ensure that the state evolution conforms to the temporal law, and the spatial constraints ensure that the states of adjacent units are continuous. The optimized state vector and land use label are output. S5. Causal intervention strategy optimization: Construct a causal intervention strategy optimization framework based on reinforcement learning, use the causal graph to identify the causal path from the management measure node to the risk node and calculate the path weight, concatenate the path weight with the state vector to form an enhanced state representation, train the strategy network based on the enhanced state representation, and output the optimal management strategy. S6. Change Detection and Time Tracing: Based on the optimized state vector, change detection is performed, and the change time is inferred from the current state through Bayesian inversion. The driving factors are traced using the causal graph. S7. Intelligent Analysis Report Generation: Integrate the analysis results from steps S1 to S6 to generate an intelligent analysis report that includes the current status of land types, evolution trends, risk assessment, and management recommendations.

[0010] Furthermore, the temporal encoder uses a bidirectional long short-term memory network (BiLSTM) to process temporal feature sequences, captures temporal dependencies through forward and backward propagation, and outputs a fixed-dimensional temporal representation vector; the causal graph structure learner maps the node representation vectors to an adjacency matrix through a multilayer perceptron, where the elements of the adjacency matrix represent the probability of causal edges between corresponding nodes, and introduces a directed acyclic graph constraint to penalize loop structures by calculating the trace of the exponential function of the adjacency matrix.

[0011] Furthermore, the physical constraint layer constructs a physical law library based on agronomic and ecological mechanisms. The physical law library defines allowed causal relationships and prohibited causal relationships. For each element in the adjacency matrix, its consistency score with the physical law library is calculated. Edges that conform to allowed relationships receive positive weights, while edges that violate prohibited relationships receive negative penalty weights. The consistency scores are then subjected to nonlinear transformation and regularization to form physical constraint terms. These terms are then weighted and summed with the prediction loss and the directed acyclic graph constraint loss to form a total loss function for optimization.

[0012] Furthermore, each component of the K-dimensional continuous vector represents the state intensity of the corresponding land type, with values ​​ranging from zero to one. The sum of all components equals one to satisfy probability constraints, thus expressing the transitional and mixed states of land types. The coefficient of the environmental response term is determined by traversing all directed edges from environmental factor nodes to land type nodes in the causal graph. If a causal edge exists from an environmental factor to a land type, the response coefficient of that environmental factor to that land type is set to a non-zero value, and its magnitude is determined according to the strength of the causal edge. If no causal edge exists, the response coefficient is set to zero. The coefficient of the competitive transformation term is calculated based on the transformation cost between land types and the distribution of land types in the spatial neighborhood.

[0013] Furthermore, the causal path weights are calculated as follows: for each directed path from the management measure node to the risk node in the causal graph, the path weight is calculated as the product of the causal strengths of all edges on the path. The shorter the path and the greater the causal strength of the edges, the higher the path weight. The enhanced state representation is obtained by concatenating the path weights of each causal path as an additional dimension with the current state vector, enabling the policy network to simultaneously perceive the current state and the causal influence mechanism.

[0014] Furthermore, the Bayesian inversion is achieved by constructing a likelihood function and a prior distribution. The likelihood function represents the probability of observing the current state under a given time of change. Based on the evolutionary model, it is integrated forward from the hypothetical time of change to the current time and the degree of agreement with the observed state is calculated. The prior distribution is set based on the statistical regularity of historical changes. The posterior probability distribution of the time of change is calculated using the Bayesian formula. The peak of the posterior distribution corresponds to the most likely time of change.

[0015] Secondly, the present invention also provides an artificial intelligence-based natural resource data analysis system, comprising: The data acquisition module is used to acquire multi-temporal remote sensing image data, environmental auxiliary data, historical land cover data and management data of the study area, and to perform quality control, preprocessing and temporal feature extraction. The temporal-causal graph neural network module is used to learn temporal representations and causal graph structures through graph neural networks, and to perform constraint optimization on the causal graph using a physical law library. The state evolution modeling module is used to establish a continuous state vector evolution model and calculate the evolution trajectory of the state over time based on differential equations. The physical constraint optimization module is used to optimize the state vector by introducing spectral constraints, temporal constraints, and spatial constraints. The intervention strategy optimization module is used to optimize management decisions based on reinforcement learning and causal inference, enhance state representation through causal paths, and output the optimal intervention plan. The change detection and tracing module is used to detect land use changes, trace the time of change through Bayesian inversion, and identify driving factors using causal graphs. The report generation module integrates the analysis results from various modules and automatically generates intelligent analysis reports. The modules are connected via a data interface. The output of the data acquisition module serves as the input of the temporal-causal graph neural network module. The causal graph structure output by the temporal-causal graph neural network module is used to guide the state evolution modeling module in determining evolution parameters, to guide the intervention strategy optimization module in constructing enhanced state representations, and to guide the change detection and tracing module in identifying driving factors. The output of the state evolution modeling module is processed by the physical constraint optimization module and then used by the intervention strategy optimization module and the change detection and tracing module.

[0016] Beneficial effects: 1. Improving recognition accuracy by integrating temporal patterns and causal mechanisms: By using a graph neural network that integrates temporal and causal relationships, the temporal variation patterns of land types and the causal influence mechanisms of environmental factors are captured simultaneously. This enables the differentiation of land types with similar spectral characteristics but different temporal variation patterns, thus solving the problem of insufficient recognition accuracy caused by not considering temporal patterns.

[0017] 2. Expressing the continuous evolution process of land use types: By replacing discrete land use type labels with continuous state vectors, the evolution process of state vectors is described by differential equations, which can characterize the transitional state and mixed pixel state of land use type changes, overcoming the shortcomings of discrete detection in being unable to express continuous transition processes.

[0018] 3. Introducing multidimensional physical constraints to ensure the rationality of results: The state vector is optimized by spectral constraints, temporal constraints and spatial constraints, so that the analysis results conform to the actual observation law, temporal change law and spatial continuity, avoiding unreasonable results caused by lack of physical constraints.

[0019] 4. Achieve causal mechanism identification and intervention strategy optimization: Utilize causal graphs to learn the causal relationship between environmental factors and land types, identify causal paths from management measures to risks, and integrate causal path information into a reinforcement learning framework to optimize intervention strategies, overcoming the shortcomings of correlation analysis in lacking causal identification and effect assessment capabilities.

[0020] 5. Improve the accuracy of time tracking of changes: By combining the Bayesian inversion method with the continuous state evolution model, the time of change and its uncertainty range can be deduced from the current state, thus overcoming the problem that the accuracy of time tracking of changes is limited by the observation time interval. Attached Figure Description

[0021] Figure 1 A flowchart illustrating the steps of the method described in this invention is shown. Figure 2 A schematic diagram of the system architecture described in this invention is shown. Detailed Implementation

[0022] Exemplary embodiments of the present invention will now be described in detail with reference to the accompanying drawings.

[0023] Combination Figure 1 In a first aspect, the present invention provides a natural resource data analysis method based on artificial intelligence, comprising the following steps: S1. Multi-source data acquisition and temporal feature extraction Acquire multi-temporal remote sensing image data, environmental auxiliary data, historical land cover data and management data of the study area, preprocess the data and extract temporal feature variables.

[0024] The remote sensing imagery data includes seven spectral bands: blue, green, red, near-infrared, shortwave infrared 1, shortwave infrared 2, and thermal infrared, with a spatial resolution of 30 meters and a temporal resolution of 16 days. Four consecutive years of imagery data for the study area were acquired, totaling 92 temporal phases, covering the complete seasonal variation cycle.

[0025] Environmental ancillary data includes digital elevation models, slope, aspect, soil type distribution maps, annual average temperature data, annual cumulative precipitation data, and annual cumulative evapotranspiration data. This data provides information on the driving factors of land cover distribution, laying the foundation for causal relationship learning.

[0026] Historical land use data includes land use labeling for known pixels, covering six categories: cultivated land, forest land, grassland, water area, construction land, and unused land. This historical land use data is used to train the supervised learning model.

[0027] The management data includes policy texts and vector data of implementation areas related to natural resource management. These policies include agricultural subsidy policies, policies on returning farmland to forest, water resource management measures, and policies on construction land control. The management data provides a realistic basis for optimizing intervention strategies.

[0028] Preprocessing is performed on remote sensing image data. Preprocessing includes radiometric calibration, atmospheric correction, and geometric correction. Radiometric calibration converts the image's digital encoded values ​​into the physical quantity of radiance; atmospheric correction converts radiance into surface reflectance; and geometric correction registers the image to a unified geographic coordinate system with a registration accuracy controlled within one pixel. The preprocessed remote sensing image has a unified radiometric and spatial reference, which can be used for subsequent feature extraction and analysis.

[0029] Extract temporal feature variables. Temporal feature variables include three categories: spectral temporal features, vegetation index temporal features, and texture features. These three types of features characterize the spectral and spatial characteristics of land cover from different perspectives.

[0030] Spectral temporal features reflect the variation of the spectral response of land cover types in different bands over time. The mean, standard deviation, maximum value, and minimum value of the time series were calculated for each of the seven bands, yielding four statistical features for each band, resulting in a total of 28 spectral temporal features for the seven bands.

[0031] The temporal characteristics of vegetation indices reflect the temporal changes in vegetation growth. The Normalized Difference Vegetation Index (NDVI), Enhanced Vegetation Index (EVI), and Soil-Regulated Vegetation Index (SAVI) were calculated. For each vegetation index, the mean, standard deviation, maximum, and minimum values ​​of the time series were calculated, resulting in 12 temporal characteristics for the three indices. The calculation of vegetation indices involves division of band reflectance; to avoid the denominator being zero, a small constant is added to the denominator. .

[0032] Texture features reflect the spatial texture patterns of land cover types. Texture features are calculated based on a gray-level co-occurrence matrix (GLCM), which statistically analyzes the co-occurrence frequency of gray values ​​for pixel pairs at specific directions and distances. In this embodiment, the calculation direction of the GLCM is set to 0 degrees (horizontal), and the distance is set to 1 pixel. Four texture features are extracted from the GLCM: contrast, correlation, energy, and homogeneity.

[0033] Contrast ratio reflects the degree of drastic change in grayscale, and the calculation formula is: Where: C represents contrast; is the number of gray levels, with a value of 256; i is the gray level value in the row direction, ranging from 1 to... j represents the column-direction grayscale value, ranging from 1 to... ; For elements of the gray-level co-occurrence matrix, represent the normalized frequency of co-occurrence of gray-level values ​​i and j in a given direction and distance.

[0034] Correlation reflects the degree of linear correlation between gray levels of pixels, and the calculation formula is: Where: R represents correlation; Let be the mean of the gray-level co-occurrence matrix in the i-th direction, calculated as follows: ; Let be the mean of the gray-level co-occurrence matrix in the j-direction, calculated as follows: ; The standard deviation in the i-direction is calculated as follows: ; Let j be the standard deviation in the j-direction, calculated as follows: ; For small constant values This is used to prevent the denominator from being zero.

[0035] Energy reflects the coarseness of the texture, and the calculation formula is: Where E is energy, also known as the second angular moment.

[0036] Homogeneity reflects the uniformity of texture, and the calculation formula is: , where: H represents homogeneity, also known as inverse difference moment; This represents the absolute value of the grayscale difference.

[0037] The extracted 28 spectral temporal features, 12 vegetation index temporal features, and 4 texture features were concatenated to form a 44-dimensional feature vector. This feature vector comprehensively characterizes the spectral properties, vegetation status, and spatial texture of pixels, providing rich information for subsequent land cover identification.

[0038] Environmental auxiliary data processing yields environmental factor feature vectors. Six environmental factors—slope, aspect, soil type, annual average temperature, annual cumulative precipitation, and annual cumulative evapotranspiration—are extracted from the environmental auxiliary data. For continuous environmental factors, including slope, temperature, precipitation, and evapotranspiration, normalization is performed by subtracting the minimum value of the factor across all pixels from the original value and then dividing by the difference between the maximum and minimum values, normalizing the numerical range to between 0 and 1. For categorical environmental factors, such as soil type, one-heat encoding is performed, converting the categorical values ​​into binary vectors. The processed six environmental factors form a 6-dimensional environmental factor feature vector.

[0039] Historical land use data is processed into land use label vectors. One-hot encoding is used to convert land use categories into 6-dimensional vectors, where only one element is 1 to indicate that the pixel belongs to the corresponding land use category, and the remaining elements are 0. For example, the one-hot encoded vector for cultivated land is... Woodland Corresponding Grassland corresponds Water area corresponding Construction land corresponding Unused land correspond .

[0040] After data acquisition and temporal feature extraction are completed, a 44-dimensional feature vector, a 6-dimensional environmental factor feature vector, and a 6-dimensional land category label vector are output as inputs for the construction step of the temporal-causal fusion graph neural network.

[0041] S2. Construction of a Temporal-Causal Fusion Graph Neural Network A temporal-causal fusion graph neural network is constructed, comprising a temporal encoder, a causal graph structure learner, and a physical constraint layer. The temporal encoder encodes temporal features to obtain temporal representation vectors. The causal graph structure learner learns the causal graph structure between environmental factors, temporal feature parameters, and land use labels. The physical constraint layer constrains the causal graph structure based on agronomic and ecological mechanisms. By jointly optimizing prediction accuracy, causal graph topological validity, and physical rationality, the causal graph structure and posterior probabilities of land use types are output.

[0042] The land cover status of natural resources is driven by both environmental factors and human management measures. This driving relationship is essentially causal rather than a simple statistical correlation. Traditional methods identify land cover based on spectral features, neglecting the causal mechanism between land cover and environmental factors. This leads to difficulties in distinguishing land cover types with similar spectral features but different evolutionary patterns. By learning the causal graph structure, the driving effect of environmental factors on land cover can be identified, thereby more accurately predicting the evolutionary trend of land cover. This approach unifies temporal patterns and causal mechanisms within an end-to-end optimization framework, discovering the causal graph structure while learning temporal representations.

[0043] The temporal encoder uses a bidirectional long short-term memory (BiLSTM) network to encode temporal feature sequences. BiLSTM includes both forward and backward LSTM directions, enabling it to simultaneously capture the forward and backward dependencies of temporal features. The hidden layer has 128 neurons.

[0044] For a 44-dimensional input feature vector The BiLSTM processes the features sequentially from the first to the 44th feature dimension. At each dimension d, the forward LSTM and backward LSTM output the hidden state, respectively. and The dimensions are all 128. The forward hidden state of the last dimension and the backward hidden state of the first dimension are concatenated to obtain the temporal representation vector: in: This is a time series representation vector with a dimension of 256; the semicolon indicates a vector concatenation operation. This represents the hidden state of the forward LSTM in the last dimension, which is 128. This represents the hidden state of the backward LSTM in the first dimension, which is 128.

[0045] The dimensionality of the 256-dimensional representation vector is reduced to 128-dimensional using a single fully connected network: in: This is the final time series representation vector, with a dimension of 128; This is the projection layer weight matrix, with dimensions of 128 x 256; This is the projection layer bias vector, with a dimension of 128; To correct the linear unit activation function, the temporal representation vector contains the most discriminative information for land use identification from the temporal features, serving as input for causal graph learning.

[0046] The causal graph structure learner maps temporal representation vectors, environmental factor feature vectors, and land category label vectors to adjacency matrices of a causal graph. Nodes in the causal graph represent variables, and directed edges represent causal relationships. Three types of nodes are defined: temporal nodes, environmental nodes, and land category nodes, corresponding to temporal representations, environmental factors, and land category labels, respectively.

[0047] For each pixel i, its temporal representation vector , Environmental factor eigenvectors Historical place category label vector These serve as the three node representations of the causal graph. A unified dimension of 128 is achieved using zero-padding, forming the node representation matrix corresponding to pixel i. The dimensions are 3 by 128.

[0048] The node representation matrix is ​​mapped to an adjacency matrix using a two-layer fully connected network. The first layer calculation process is as follows: ,in: The hidden layer output for pixel i has a dimension of 64 by 3; This is the first-layer weight matrix, with dimensions of 64 x 128; This is the first-level bias vector, with a dimension of 64; Let be the transpose of the node representation matrix of pixel i, with dimensions 128 x 3.

[0049] The second layer uses the Sigmoid function to map the hidden layer output to probability values ​​between 0 and 1: in: Let be the adjacency matrix of pixel i, with dimensions 3 by 3; This is the second-layer weight matrix, with dimensions of 3 x 64; This is the second-layer bias vector, with a dimension of 3; It is an exponential function. The elements of the adjacency matrix. This represents the probability that there is a causal edge from node j to node k, where j is the row index ranging from 1 to 3 and k is the column index ranging from 1 to 3.

[0050] Cause-effect graphs need to satisfy the topological constraints of directed acyclic graphs. Constraint terms are introduced into the loss function to penalize adjacency matrix configurations that produce loops. These constraint terms are calculated based on the trace of an exponential function of the adjacency matrix. ,in: For directed acyclic graph constraints; This is the trace operation of a matrix, which is the sum of the elements on the main diagonal; It is a 3x3 identity matrix; 3 is the element-wise square of the adjacency matrix of pixel i; 3 is the number of nodes. This constraint term is a second-order Taylor expansion approximation of the matrix exponent, which is relevant when loops exist in the causal graph. Taking a positive value, the loop can be eliminated by minimizing this constraint term.

[0051] The physical constraint layer integrates mechanistic knowledge from agronomy and ecology into the causal graph learning process. A physical law library is constructed based on domain knowledge, explicitly defining permitted and prohibited causal relationships.

[0052] Permissible causal relationships include: causal edges from environmental nodes to time-series nodes, causal edges from environmental nodes to land type nodes (including the effects of precipitation on grassland growth, temperature on farmland growth, slope on forest distribution, soil type on arable land distribution, and evapotranspiration on vegetation growth), causal edges from time-series nodes to land type nodes, and self-feedback edges from time-series nodes. Each permissible causal relationship is assigned a positive weight value, ranging from 0.3 to 0.9, with higher-confidence causal relationships assigned larger weights.

[0053] Prohibited causal relationships include: causal edges from land type nodes to environmental nodes, causal edges from land type nodes to time series nodes, causal edges from temperature to water distribution, causal edges from precipitation to construction land, and causal edges from slope aspect to water distribution. Each prohibited causal relationship is assigned a negative weight value, ranging from -1.0 to -0.4. Prohibited relationships with severe violations are subject to greater penalties.

[0054] For each element of the adjacency matrix Based on the types of node j and node k, the corresponding causal relationship rules are searched in the physical law database, and the consistency score is calculated. If the causal edge from node j to node k is an allowed causal relationship, then Take the corresponding positive weight value; if it belongs to a prohibited causal relationship, then... Take the corresponding negative weight value; if the rule base does not explicitly specify this, then... The value is 0. The physical constraint term is calculated as follows: in: For physical constraints; The consistency score is calculated from node j to node k. This is the nonlinear enhancement coefficient, with a value of 0.5. It is the hyperbolic tangent function; This is the saturation coefficient, with a value of 3.0; This is the L2 regularization coefficient, with a value of 0.01; This is the penalty coefficient for bidirectional edges, with a value of 0.02. The term calculates the product of the probabilities of causal edges from node j to node k and from node k to node j. The first term is a consistency score weighting term with nonlinear enhancement, and the negative sign indicates that causal edges that conform to physical mechanisms are encouraged; the second term is an L2 regularization term to prevent the adjacency matrix elements from being too large; the third term is a bidirectional edge penalty term to ensure the unidirectionality of causal relationships.

[0055] The temporal-causal fusion graph neural network achieves coordinated optimization of its three modules through a joint loss function. The total loss function is: ,in: Total loss; To predict losses, measure the accuracy of land use predictions; The weight coefficient for the directed acyclic graph constraint is 0.1. This is the weighting coefficient for the physical constraints, with a value of 0.05.

[0056] The prediction loss uses the cross-entropy loss function: Where: N is the number of training pixels; K is the number of land use categories, with a value of 6; i is the pixel index; c is the land use category index; This is the true label for pixel i belonging to land category c, and its value is 0 or 1. Let be the predicted posterior probability that pixel i belongs to land class c; This is the operation for the natural logarithm.

[0057] The posterior probability of land use is calculated using a graph convolutional network based on the adjacency matrix and node representations. The graph convolutional network consists of two layers. The calculation process of the first layer of graph convolution is as follows: in: Let i be the initial node representation matrix, i.e., the node representation matrix. The dimensions are 3 by 128; This is the node representation matrix after the first layer of graph convolution, with a dimension of 3 x 128; The normalized adjacency matrix is ​​calculated as follows: ; This is the weight matrix for the first layer, with dimensions of 128 by 128.

[0058] The calculation process for the second layer of graph convolution is as follows: in: This is the node representation matrix after the second layer of graph convolution, with a dimension of 3 x 128; This is the weight matrix for the second layer, with dimensions of 128 by 128.

[0059] The posterior probability of a land class is calculated from its representation vector, and the Softmax function is used to normalize the multi-class probability. in: Let be the node representation value of land class c corresponding to pixel i, from Extract from the land category node row, i.e., the 3rd row; For pixel i, the corresponding land class The node representation value; For summation index.

[0060] The Adam optimizer is used to minimize the total loss function via gradient descent. The learning rate is set to 0.001, the batch size is 256, and the training iterations are 5000. During training, the weight parameters of the BiLSTM, the projection layer, the causal graph learner, and the graph convolutional network are updated synchronously to achieve end-to-end joint optimization.

[0061] After training, the optimal causal graph structure is obtained. The posterior probability of land class is calculated. Edges with an element greater than 0.5 in the causal graph structure are considered significant causal edges. These causal edges and their corresponding nodes are extracted to form the final causal graph. The causal graph clearly depicts the driving effect of environmental nodes on temporal nodes and the influence of temporal nodes on land class nodes, providing guidance on causal mechanisms for subsequent continuous state evolution modeling and intervention strategy optimization. The posterior probability of land class is output as the result of land class identification and also serves as the initialization basis for the continuous state vector.

[0062] S3, Continuous State Vector Evolution Modeling A continuous state vector evolution model is established, which uses continuous vectors to represent the land use state intensity of spatial units. The evolution process of the state vector is described based on differential equations. The coefficients of the environmental response term of the differential equations are determined according to the causal graph structure, and the time evolution trajectory of the state vector is output.

[0063] Traditional land cover identification methods employ discrete classification, assigning each pixel to a single land cover category. This approach fails to capture the transitional and mixed states of land cover. The evolution of natural resources through land cover is a continuous and gradual process. For example, the transformation from degraded grassland to sandy land involves a transitional phase of gradually decreasing grassland cover and gradually increasing bare soil area. Discrete classification cannot accurately depict this mixed state, while continuous state vectors can quantitatively represent the degree to which a pixel simultaneously belongs to multiple land cover categories.

[0064] The continuous state vector uses a K-dimensional vector to represent the land use state intensity distribution of a pixel, where K is the number of land use categories, and in this embodiment, it is set to 6. The state vector of pixel i at time t is denoted as... The vector dimension is 6. The c-th component of the state vector is denoted as... , represents the state intensity of pixel i at time t corresponding to land type c, where c is the land type index, with values ​​from 1 to 6. The state intensity value ranges from 0 to 1, where 0 indicates that the pixel does not belong to the land type at all, and 1 indicates that the pixel belongs to the land type at all.

[0065] The sum of all components of the state vector equals 1, satisfying the probability normalization constraint: ,in: The summation symbol is used; K represents the number of land use categories, with a value of 6.

[0066] The initial values ​​of the state vector are set based on the posterior probability of land type obtained in S2. For pixel i at the initial time... State vector: ,in: For pixel i at the initial time The state intensity corresponding to land category c; This is the initial time. Let S2 be the posterior probability that pixel i belongs to land class c, predicted by the graph neural network. Through this initialization method, the state vector inherits the temporal features and causal relationship information learned by the graph neural network.

[0067] The evolution of the state vector is described by ordinary differential equations. The rate of change consists of three parts: an intrinsic growth term, an environmental response term, and a competition transition term. The complete form of the differential equation is: in: Let be the derivative of the state vector with respect to time t, and its dimension is 6; Let i be the state vector of pixel i at time t; Let i be the environmental factor feature vector of pixel i, with a dimension of 6; The adjacency matrix of the causal graph learned in S2 has a dimension of 3 by 3; Let i be the set of spatial neighborhood state vectors of pixel i, containing the state vectors of 8 neighboring pixels; This is an intrinsic growth term, reflecting the natural growth trend of land use status; This is an environmental response item, reflecting the driving effect of environmental factors on land cover status through causal mechanisms; As a competitive transformation term, it reflects the competitive and transformation relationships between different land types and the impact of spatial neighborhood.

[0068] The intrinsic growth function adopts the logistic growth model to reflect the natural growth pattern of land use status without external disturbances. The intrinsic growth term for land use c is calculated as follows: ,in: This is the intrinsic growth term for land category c; This represents the intrinsic growth rate of land type c, expressed as negative first power in years, ranging from 0 to 0.5. Larger values ​​are used for rapidly succeeding grasslands and unused land, such as 0.4 to 0.5; smaller values ​​are used for stable forest land and water areas, such as 0.1 to 0.2; and medium values ​​are used for cultivated land and construction land, such as 0.2 to 0.3. The current state intensity of land type c; This is a saturation term, ensuring that the state strength does not exceed 1.

[0069] The vector form of the intrinsic growth function is: Where: the superscript T indicates transpose; arrive These are the intrinsic growth items for cultivated land, forest land, grassland, water areas, construction land, and unused land, respectively.

[0070] The environmental response function characterizes the influence of environmental factors on land cover status through causal mechanisms. A response coefficient matrix from environmental factors to land cover is established. The matrix is ​​6x6 in dimension, with row indices corresponding to 6 environmental factors and column indices corresponding to 6 land use categories. Matrix elements. This represents the response coefficient of environmental factor m to land type c. m is the index of environmental factor, from 1 to 6, corresponding to slope, aspect, soil type, annual average temperature, annual cumulative precipitation, and annual cumulative evapotranspiration, respectively.

[0071] The response coefficient is determined as follows: if a causal edge exists in the causal graph from an environmental factor to a land type, then the response coefficient of that environmental factor to that land type is non-zero, and its magnitude is determined by the strength of the causal edge; otherwise, the response coefficient is zero. For environmental factor m and land type c, the adjacency matrix of the causal graph learned by S2 is checked. Elements from environment nodes to land category nodes. Elements in the adjacency matrix. As the basic weight, the final response coefficient is determined by combining the ecological relationship between environmental factor m and land type c.

[0072] The response coefficient is set to a range of -1 to +1, with a larger absolute value indicating a more significant impact of the environmental factor on the land type. Specific values ​​are determined by analyzing the correlation between historical land type change data and environmental factors in the study area. The calculation method involves linear regression of the change in the state intensity of each land type with the normalized values ​​of each environmental factor. The regression coefficient serves as the initial value for the response coefficient, which is then adjusted using knowledge from agronomy and ecology. The adjustment principle ensures that the sign of the response coefficient aligns with the known causal relationship, and that the numerical value is proportional to the strength of the corresponding edge in the adjacency matrix of the causal graph.

[0073] The calculation process of the environmental response function for land category c is as follows: in: For land category c, the environmental response item is M; M is the number of environmental factors, with a value of 6. The response coefficient of environmental factor m to land type c is read from the response coefficient matrix; is the normalized value of the m-th environmental factor of pixel i, and its value ranges from 0 to 1; It is a saturated term.

[0074] The vector form of the environmental response function is: The competitive transformation function characterizes the competition and transformation relationships between land use types and the influence of spatial neighborhood. The calculation process of the competitive transformation function for land use type c includes two terms: the first is the inter-land use transformation term, and the second is the spatial diffusion term. in: This is a competing transformation term for land category c; For other land use categories, index from 1 to K, but not equal to c; To be classified by land The conversion rate to land category c, expressed in negative first power of years; From land category C to land category The conversion rate; The term represents the interaction strength when two land types coexist; This is the spatial neighborhood influence coefficient, with a value of 0.1. Let be the set of 8 neighboring pixels of pixel i; j is the index of the neighboring pixels. The spatial weight of pixel j to pixel i is calculated as follows: And normalization is required. ; Let Euclidean distance be the distance between the centers of pixel i and pixel j, for neighboring pixels in orthogonal directions. Meters, for neighboring pixels in the diagonal direction A meter is approximately equal to 42.43 meters; The term calculates the difference in state intensity between neighboring pixels and the central pixel on land category c.

[0075] Land Class Conversion Rate Determined based on the land use conversion cost matrix. Land use conversion cost matrix It is a 6x6 matrix, matrix elements This indicates a conversion from land category C to land category D. The cost of conversion. The numerical range of the conversion cost is set from 1 to 5, where 1 indicates easy conversion and 5 indicates difficult conversion.

[0076] The conversion rate is inversely proportional to the conversion cost, and this mapping relationship can be established using an exponential function: in: From land category C to land category The conversion rate is expressed in negative one-th power of years; This is the conversion rate scaling factor, with a value of 0.5; It is an exponential function, with the base being the base e of the natural logarithm; From land category C to land category The conversion cost.

[0077] The vector form of the competition transition function is: The evolution of the state vector is achieved through numerically solving the differential equations. A fourth-order Runge-Kutta method is used to numerically integrate the differential equations; this method has fourth-order accuracy and maintains high computational accuracy even with large time steps. The time step is set to... A year is 36.5 days. From the initial moment Initially, the state vector is iteratively calculated at each time step until the target time. .

[0078] The computation of the fourth-order Runge-Kutta method involves a fourth function evaluation and a weighted average. The first evaluation is performed at the current time step. conduct: in: The state change obtained in the first evaluation has a dimension of 6. Let be the set of neighborhood state vectors of pixel i at time t.

[0079] The second assessment at time Status is proceed at: in: The state change obtained from the second evaluation has a dimension of 6.

[0080] The third assessment at the time Status is proceed at: in: The state change obtained from the third evaluation has a dimension of 6.

[0081] The fourth assessment at time Status is proceed at: in: The state change obtained from the fourth evaluation has a dimension of 6.

[0082] The state vector at the next time step is calculated by weighted averaging: in: For pixel i at time... State vector; weights , , , These weights, corresponding to the contributions of the four evaluations, enable the Runge-Kutta method to achieve fourth-order accuracy.

[0083] To ensure that the state vector always satisfies the probability normalization constraint during evolution, the state vector is normalized after each update: in: This represents the assignment operation; the numerator is the state intensity before normalization; the denominator is the sum of the state intensities of all land types.

[0084] By using a continuous state vector evolution model, the land use status distribution of a pixel at any given time can be obtained, and the temporal evolution trajectory of the state vector can be output. The evolution trajectory is represented as a time series. It reflects the gradual process of land use status change, can identify early signals of land use change, provide a basis for timely management measures, and also serve as input for subsequent multidimensional physical constraint optimization steps.

[0085] S4, Multidimensional Physical Constraint Optimization The state vector is optimized by introducing spectral constraints, temporal constraints, and spatial constraints, and the optimized state vector and land class labels are output.

[0086] The state vector evolution trajectory obtained in S3 is based on a differential equation model prediction. Although it considers mechanisms such as intrinsic growth, environmental response, and competitive transition, it does not directly utilize remote sensing observation data for correction. Remote sensing imagery provides observational information about the true state of the Earth's surface. By matching the state vector with the observational data, the bias in the model prediction can be corrected, improving the accuracy of the state vector. Physical constraint optimization utilizes three constraints—spectral mixture model, temporal continuity, and spatial continuity—to incorporate information from the observational data into the state vector.

[0087] The spectral constraints are based on a linear mixture model, assuming that the observed spectrum of a pixel is the result of a state-intensity weighted mixture of the pure spectra of each land type. The calculation process for the spectral constraint terms is as follows: in: This is a spectral constraint term, dimensionless; The band indexes, from 1 to 7, correspond to blue light, green light, red light, near-infrared, short-wave infrared 1, short-wave infrared 2, and thermal infrared bands, respectively. For pixel i in band The observed reflectance, obtained from the preprocessed data of S1, is dimensionless; For land type c in band The reference spectral reflectance is dimensionless and is obtained by averaging the reflectance of all pixels belonging to the same land type in each band in the training samples, selecting pixels with pure state (i.e., single land type state intensity greater than 0.9). This is the c-th component of the state vector; Bands for state vector reconstruction Reflectivity.

[0088] Temporal constraints ensure that the evolution of the state vector conforms to the laws of temporal change. Changes in land use status should be smooth and continuous; there should be no drastic jumps in the state vector between adjacent time points, because the evolution of natural processes is usually gradual. The calculation process for the temporal constraint term is as follows: in: This is a time-series constraint term, dimensionless; The L2 norm is used to calculate the Euclidean length of a vector, defined as follows: ; For first-order temporal constraints, the difference in state vectors between adjacent time steps is penalized; This is the second-order temporal constraint weight coefficient, with a value of 0.5; As a second-order temporal constraint, it penalizes the second-order difference of the state vector. This term is the second-order derivative in discrete form, ensuring that the curvature of the evolution trajectory is not too large.

[0089] Spatial constraints ensure the spatial continuity of state vectors between adjacent cells. Land use types typically exhibit clustering characteristics, meaning that the same land use type tends to be distributed in contiguous areas. The calculation process for the spatial constraint term is as follows: in: This is a spatial constraint term, dimensionless. Let be the set of 8 neighboring pixels of pixel i; j is the index of the neighboring pixels. For zero-order spatial constraints, penalize the difference in state vectors between the center pixel and its neighboring pixels; This is the first-order spatial constraint weight coefficient, with a value of 0.3; Let be the spatial gradient of the state vector of pixel i at time t, with a dimension of 6; This is a first-order spatial constraint that penalizes the difference in the gradient of the state vector space.

[0090] Spatial gradient of the state vector The central difference method is used for calculation, which has second-order accuracy. The calculation process is as follows: in: It is the partial derivative of the state vector in the x-direction, i.e., the east-west direction, with a dimension of 6; It is the partial derivative of the state vector in the y-direction, i.e., the north-south direction, with a dimension of 6; , , , These are the orthogonal neighbor pixel indices of pixel i in the east, west, north, and south directions, respectively; The pixel spacing in the east-west direction is set to 30 meters. The pixel spacing in the north-south direction is set to 30 meters. This is the element-wise square root operation.

[0091] The spectral constraints, temporal constraints, and spatial constraints are integrated into a unified physical constraint optimization objective function: in: This refers to the total physical constraints. This is the spectral constraint weighting coefficient, with a value of 1.0; This is the weighting coefficient for the time-series constraints, with a value of 0.5. This is the spatial constraint weight coefficient, with a value of 0.8. The specific value of the weight coefficient is determined by performing a grid search on the validation set, with a search range of 0.1 to 2.0 and a step size of 0.1.

[0092] The state vector is optimized by minimizing the total physical constraints. The gradient descent method is used to iteratively optimize the state vector. in: This is the state vector at the k-th iteration, with a dimension of 6; This is the state vector at the k+1th iteration; k is the iteration index, starting from 0. To optimize the learning rate, a value of 0.01 was set. The gradient of the physical constraint term with respect to the state vector has a dimension of 6 and is obtained through automatic differentiation.

[0093] The gradient of the physical constraint term with respect to the state vector consists of three parts. The gradient of the spectral constraint is: in: is the partial derivative of the spectral constraint term with respect to the c-th component of the state vector; the coefficient 2 comes from the derivative of the squared term; To sum the index, traverse all land classes.

[0094] The gradient of a temporal constraint involves the state vectors of adjacent time steps. For a first-order temporal constraint: For second-order timing constraints: The gradient of a spatial constraint involves the state vectors of neighboring pixels. For a zeroth-order spatial constraint: The optimization iteration continues until the change in the physical constraints is less than the convergence threshold or the maximum number of iterations is reached. The convergence condition is: ,in: This represents the value of the physical constraint term at the k-th iteration. The convergence threshold is set to a value of [value to be filled in]. The maximum number of iterations is set to 100.

[0095] After each iteration, the state vector is normalized and boundary constraints are applied to ensure that the probability constraints are always satisfied. The normalization process is as follows: in: Ensure that the state strength is non-negative; normalization ensures that the sum of all components is 1.

[0096] After optimization, output the optimized state vector. Land category labels. Land category labels are determined by selecting the land category corresponding to the largest component in the state vector: in: The optimized land class label for pixel i; This is an operation to retrieve the index corresponding to the maximum value; Let c be the c-th component of the optimized state vector.

[0097] Multidimensional physical constraint optimization ensures that the state vector satisfies the data-driven evolution law while conforming to the physical constraints of spectral observation, temporal continuity, and spatial clustering, thus improving the reliability of land cover identification and evolution prediction. The optimized state vector and land cover labels serve as inputs for intervention strategy optimization and change detection steps.

[0098] S5, Optimization of Causal Intervention Strategies A causal intervention strategy optimization framework based on reinforcement learning is constructed. Causal graphs are used to identify causal paths and calculate path weights. The path weights are concatenated with state vectors to form an enhanced state representation. The policy network is then trained to output the optimal management strategy.

[0099] Natural resource management requires human intervention to influence land cover evolution processes in order to achieve the goals of ecological protection and resource utilization. The formulation of intervention strategies necessitates a quantitative assessment of the causal effects of different management measures on land cover changes and the optimization of multi-measure combinations under budget constraints. Based on the causal graph structure learned in S2, a reinforcement learning framework is used to optimize causal intervention strategies.

[0100] Causal intervention strategy optimization models the management decision-making process as a Markov decision process. The state is defined as the state vector of pixel i at time t. With environmental factor eigenvectors splicing: in: The state represents the Markov decision process, with a dimension of 6 plus 6 equaling 12; the semicolon indicates a vector concatenation operation.

[0101] An action is defined as a subset of the set of optional management measures. These measures include six categories: returning farmland to forest, returning farmland to grassland, water resource allocation, agricultural subsidies, vegetation restoration, and construction land control. (Action vector) Given a 6-dimensional binary vector, the m-th component... A value of 1 indicates that the m-th measure is implemented, and a value of 0 indicates that it is not implemented. m is the index of the management measure, ranging from 1 to 6.

[0102] The reward function measures the improvement in state after taking an action, minus management costs. in: To take an action at time t The reward received later; This is the state vector at the next time step obtained after taking an action and evolving through one time step. The target state vector has a dimension of 6; The squared Euclidean distance between the state vector and the target state vector is given. This is the cost weighting coefficient, with a value of 1.0; Let m be the unit cost of the management measure, expressed in yuan. This is the highest cost among all management measures; The normalized cost; Let m be the m-th component of the action vector.

[0103] The impact of management measures on the evolutionary process is achieved by modifying the response coefficients in the environmental response function: in: The environmental response coefficient after implementing management measures; The original environmental response coefficient; The influence intensity coefficient of management measure m ranges from 0.1 to 0.5. This indicates the impact of management measure m on land category c, and its value is either 0 or 1. Let m be the m-th component of the action vector.

[0104] The causal path enhancement state representation utilizes the causal graph learned in S2 to identify causal paths from management measures to target land classes. For management measure m and target land class c, all directed paths from the management measure node to the land class c node are found in the causal graph. Path weights are calculated as the product of the probabilities of all causal edges on the path: in: The causal path weight from management measure m to land category c, with a value ranging from 0 to 1; Let represent a directed path from management measure node m to land category node c; e is an edge on the path; Let e ​​be the causal edge probability corresponding to edge e; This is a series multiplication operation.

[0105] The causal path weights of all management measures to all land types are organized into a causal path matrix. The dimension is 6 by 6. Multiplying the causal path matrix by the action vector yields the comprehensive causal impact vector of management measures on various land types: ,in: The comprehensive causal impact vector of management measures has 6 dimensions; This is a causal path matrix with a dimension of 6 by 6; This is the action vector, with a dimension of 6.

[0106] The causal influence vector is concatenated with the state vector and environmental factor vector to form an enhanced state representation: ,in: To enhance state representation, the dimension is 6 plus 6 plus 6 equals 18. This enhanced state representation includes causal impact information of management actions, enabling the policy network to simultaneously perceive the current state and the effects of management actions through causal mechanisms.

[0107] The policy network outputs the probability distribution of the optimal action based on the enhanced state representation. The policy network consists of three fully connected layers: the first layer has 64 neurons, the second layer has 32 neurons, and the output layer has 64 neurons, corresponding to 64 possible action combinations. The forward computation process of the policy network is as follows: in: This is the output of the first layer, with a dimension of 64; This is the first-layer weight matrix, with dimensions of 64 x 18; This is the first-level bias vector, with a dimension of 64; To modify the activation function of the linear unit; This is the output of the second layer, with a dimension of 32; This is the second-layer weight matrix, with dimensions of 32 x 64; This is the second-layer bias vector, with a dimension of 32; Take an action given an enhanced state representation The probability of; The output layer weight matrix has a dimension of 64 x 32; This is the output layer bias vector, with a dimension of 64; It is an exponential function; It represents the set of all possible actions, containing 64 actions. For summation index.

[0108] The policy network is trained using the proximal policy optimization algorithm. The objective function of the proximal policy optimization algorithm is: in: Optimize the objective function for the near-end strategy; These are the parameters of the policy network; For the expectation calculation of time step t; The strategy probability ratio is calculated as follows: ; This is the estimated value of the dominance function; The truncation function is defined as follows: ; The cutoff threshold is set to 0.2. This is for calculating the minimum value.

[0109] The advantage function is calculated using the generalized advantage estimation method: in: This is the discount factor, with a value of 0.99; The parameter for estimating generalized dominance is 0.95. This is the time step offset; For a moment The timing difference error is calculated as follows: ; For a moment The rewards received; It is a value function.

[0110] The calculation process of the value network is as follows: in: State value, scalar; , , These are the weight matrices for the three layers of the value network, with dimensions of 64 x 18, 32 x 64, and 1 x 32, respectively. , , These are the bias vectors for the three layers, with dimensions of 64, 32, and 1 respectively.

[0111] Value networks are trained by minimizing the error of the value function: in: Loss is a value function. The parameters of the value network; The output of the value network; For cumulative discount rewards, the calculation is as follows: .

[0112] The policy network and value network are updated alternately, with a batch size of 128, a learning rate of 0.0003, and 10,000 training iterations. Training data is collected by executing the current policy in the environment, collecting empirical data for 2048 time steps each time.

[0113] After training, the policy network can output the optimal management strategy based on the enhanced state representation. The output optimal management strategy includes the combination of management measures selected at each time step, the expected land use evolution trajectory, the expected ecological benefits, and the management costs, providing support for management decisions.

[0114] S6. Change Detection and Time Tracing Change detection is performed based on the optimized state vector, the change time is inferred through Bayesian inversion, and the driving factors are traced using a causal graph.

[0115] Change detection is performed based on the optimized state vector in S4. For pixel i, its dominant land class at different times is calculated: in: The dominant geographic class of pixel i at time t; This is an operation to retrieve the index corresponding to the maximum value; is the c-th component of the state vector after physical constraint optimization; K is the number of land categories, with a value of 6.

[0116] A change in land use type is considered a land use type change when the dominant land use type changes. The criteria for determining a change are: in: This is a marker for detecting changes in pixel i at time t, where 1 indicates a change and 0 indicates no change. For reference time; The state strength threshold is set to 0.6. Duration of dominant land type; The minimum duration threshold is 0.5 years.

[0117] The change time tracing uses the Bayesian inversion method. The form of Bayes' theorem is: in: Given the current state vector, this represents the posterior probability distribution of the time of change. The time of variation is defined as the range from the reference time. up to the current moment ; This is the observation state vector at the current moment, with a dimension of 6; It is the likelihood function; For the prior probability distribution of the time-varying variables; The marginal probability of the observed state.

[0118] The likelihood function is calculated by forward integration of the evolutionary model and takes the Gaussian likelihood form: in: The standard deviation of the observed noise is set to 0.1. This is the normalization constant; It is an exponential function; The squared Euclidean distance between the predicted state vector and the observed state vector; is the variance parameter of the Gaussian distribution.

[0119] Predicted state vector The calculation process is as follows: during the time of change The state vector undergoes a sudden change, from the changed state vector Initially, the evolutionary model of S3 is numerically integrated to the current time step using the fourth-order Runge-Kutta method. : in: This represents the numerical integration operation of the fourth-order Runge-Kutta method; The time step is 0.1 years. , and These are the intrinsic growth function, environmental response function, and competition transition function defined in S3, respectively.

[0120] The prior probability distribution adopts a hybrid form of uniform and exponential distribution: in: The weight is mixed, and the value is 0.7. The density function is a uniform distribution. The rate parameter is for the exponential distribution, and its value is 2.0 years to the power of negative one. This is the normalization factor.

[0121] The calculation process for the posterior probability distribution is as follows: in: To obtain the sum index, iterate through all possible time changes.

[0122] The peak of the posterior probability distribution corresponding to the most likely time of change: in: This is the most likely time estimate of change.

[0123] The range of uncertainty in time is represented by the confidence interval of the posterior probability distribution. Calculate the 95% confidence interval, which is the smallest time interval containing 95% of the posterior probability quality.

[0124] Driving factor tracing utilizes the causal graph learned in S2 to identify the driving factors leading to land cover change. The changes in environmental factors before and after the change are calculated. in: This represents the change in environmental factor m; m is the environmental factor index, from 1 to 6. For moments of change The value of environmental factor m; For reference time The value of environmental factor m.

[0125] Calculate the contribution of each environmental factor to land cover change: in: denoted as m, representing the contribution of environmental factor m to land cover change; c represents the changed land cover. This represents the absolute value of the causal edge weight; It represents the absolute value of the change in environmental factors; The sign function of the response coefficient; denoted as the response coefficient of environmental factor m to land type c.

[0126] The contributions of all environmental factors were ranked, and the top three environmental factors with the largest contributions were selected as the main driving factors.

[0127] The output includes the results of change detection and time tracking, including the change occurrence identifier, the most likely change time, the change time confidence interval, the land type before and after the change, and the main driving factors, providing a basis for management decisions and report generation.

[0128] S7, Intelligent Analysis Report Generation Integrating the analysis results from steps S1 to S6, an intelligent analysis report is generated that includes the current status of land types, evolution trends, risk assessment, and management recommendations. Report generation involves three stages: data statistics, chart creation, and text generation.

[0129] The data statistics process calculates the area and proportion of each land type: The area of ​​land category C is calculated as follows: ,in: , where is the area of ​​land category c, in square meters; N is the total number of pixels in the study area; Let be the optimized state intensity of class c corresponding to pixel i at time t; The area of ​​a single pixel is calculated as 30 meters by 30 meters, which equals 900 square meters.

[0130] The land category proportion is calculated as follows: ,in: The area percentage of land category c is dimensionless. For summation index.

[0131] The number of changed pixels is calculated as follows: ,in: The number of pixels that changed; For indicator functions; This is a change detection identifier for pixel i.

[0132] The elements of the land category transformation matrix are calculated as follows: ,in: To convert from land category C to land category The number of pixels; For pixel i, the land class before the change; The land class after the change of pixel i; For logical AND operator.

[0133] The area of ​​land class c at the future time is calculated as follows: ,in: Let be the predicted area of ​​land category c at future times; For pixel i at future time The state intensity corresponding to land class c.

[0134] High-risk area identification. The high-risk indicator is calculated as follows: ,in: 6 represents the risk identifier for pixel i; 6 represents the land use index for unused land. The risk threshold is set to 0.15.

[0135] Management effectiveness evaluation. The change in area of ​​land category C after the implementation of management measures is calculated as follows: ,in: This represents the change in area of ​​land category C after the implementation of management measures; To determine the area of ​​land category c under the intervention scenario; This represents the area of ​​land class c under the baseline scenario.

[0136] The charting process generates a current land use distribution map, a land use change detection map, a land use conversion chord diagram, an evolution trend curve, and a risk area distribution map. The current land use distribution map is a spatial raster map, with each cell assigned a different color based on the dominant land use type. The land use change detection map labels cells that have changed. The land use conversion chord diagram shows the conversion relationships between land use types, with the line thickness indicating the number of cells involved in the conversion. The evolution trend curve shows the change in area of ​​each land use type over time. The risk area distribution map labels high-risk cells.

[0137] The text generation stage generates the report text based on statistical data and predefined templates. The report consists of four parts: title, abstract, main body, and conclusion. The title concisely summarizes the report's theme. The abstract summarizes the main findings, including the total area of ​​the study region, the distribution of major land use types, the area of ​​change, the main types of change, future trend predictions, and the area of ​​high-risk areas. The main body is divided into five sections: current status of land use types, changes in land use types, evolutionary trends, risk assessment, and management recommendations. The conclusion summarizes the core findings and key recommendations.

[0138] The report outputs in a structured document format, supporting PDF, Word, and HTML formats. It offers interactive query functionality, allowing management decision-makers to click on any cell on the map to view detailed analysis results, including the state vector evolution trajectory, temporal probability distribution of changes, and the contribution of driving factors.

[0139] Combination Figure 2 Secondly, the present invention provides a natural resource intelligent analysis system based on time-series-causal fusion and counterfactual reasoning, comprising: The data acquisition module is used to acquire multi-time series remote sensing image data, environmental auxiliary data, historical land cover data and management data. It performs radiometric calibration, atmospheric correction and geometric correction preprocessing on remote sensing image data, extracts spectral time series features, vegetation index time series features and texture features, processes environmental auxiliary data and historical land cover data, and outputs 44-dimensional feature vector, 6-dimensional environmental factor feature vector and 6-dimensional land cover label vector.

[0140] The temporal-causal graph neural network module receives the feature vector and label vector output by the data acquisition module. It encodes the temporal features through a bidirectional long short-term memory network to obtain the temporal representation vector. It learns the causal graph structure among environmental factors, temporal feature parameters, and land use labels through a causal graph structure learner. It constrains the causal graph structure based on agronomic and ecological mechanisms through a physical constraint layer. It jointly optimizes the prediction accuracy, causal graph topological validity, and physical rationality, and outputs the causal graph structure and the posterior probability of land use.

[0141] The state evolution modeling module receives the causal graph structure and land use posterior probabilities output by the time-series causal graph neural network module. It initializes the state vector based on the land use posterior probabilities and describes the evolution process of the state vector through ordinary differential equations. The ordinary differential equations include intrinsic growth terms, environmental response terms, and competing transition terms. The coefficients of the environmental response terms are determined according to the causal graph structure. The differential equations are numerically solved using the fourth-order Runge-Kutta method, and the time evolution trajectory of the state vector is output.

[0142] The physical constraint optimization module receives the temporal evolution trajectory of the state vector output by the state evolution modeling module. It introduces spectral constraints, temporal constraints, and spatial constraints to optimize the state vector. Spectral constraints ensure that the state vector is consistent with the observed spectrum, temporal constraints ensure that the evolution conforms to the temporal pattern, and spatial constraints ensure that the states of adjacent pixels are continuous. The total physical constraint term is minimized through the gradient descent method, and the optimized state vector and land cover label are output.

[0143] The intervention strategy optimization module receives the causal graph structure output by the temporal-causal graph neural network module and the optimized state vector output by the physical constraint optimization module. It models the management decision-making process as a Markov decision process, uses the causal graph to identify the causal path from the management measures to the target land type and calculates the path weights. The path weights are concatenated with the state vector to form an enhanced state representation. The policy network outputs the probability distribution of the optimal action based on the enhanced state representation. The policy network is trained using a proximal policy optimization algorithm to output the optimal management strategy.

[0144] The change detection and tracing module receives the optimized state vector from the physical constraint optimization module and the causal graph structure from the temporal-causal graph neural network module. It calculates the dominant land use type of a pixel at different times and determines that a change has occurred. It uses the Bayesian inversion method to infer the change time, calculates the likelihood function through forward integration of the evolutionary model, calculates the posterior probability distribution by combining the prior probability distribution, and uses the causal graph to identify the driving factors that cause land use type changes. It outputs the change occurrence identifier, the most likely change time, the change time confidence interval, the land use type before and after the change, and the main driving factors.

[0145] The report generation module receives the optimized state vector and land category labels output by the physical constraint optimization module, the optimal management strategy output by the intervention strategy optimization module, and the change detection results output by the change detection and tracing module. It calculates the area and proportion of each land category, the number of changed pixels, the land category transformation matrix, the predicted area of ​​each land category at future times, the area of ​​high-risk areas, and the management effectiveness assessment. It generates a land category status distribution map, a land category change detection map, a land category transformation chord diagram, an evolution trend curve, and a risk area distribution map. Based on statistical data and predefined templates, it generates report text and outputs an intelligent analysis report in a structured document format.

[0146] The modules are interconnected via data flow. The output of the data acquisition module serves as the input to the temporal-causal graph neural network module. The output of the temporal-causal graph neural network module includes the causal graph structure and the posterior probability of land use types. The causal graph structure is passed to the state evolution modeling module, the intervention strategy optimization module, and the change detection and tracing module, while the posterior probability of land use types is passed to the state evolution modeling module. The output of the state evolution modeling module is passed to the physical constraint optimization module. The output of the physical constraint optimization module includes the optimized state vector and land use type labels. The optimized state vector is passed to the intervention strategy optimization module and the change detection and tracing module, while the land use type labels are passed to the report generation module. The output of the intervention strategy optimization module is passed to the report generation module, and the output of the change detection and tracing module is passed to the report generation module.

[0147] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any changes or substitutions made by those skilled in the art based on the technical solutions of the present invention should be covered within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be determined by the scope defined in the claims.

Claims

1. A natural resource data analysis method based on artificial intelligence, characterized in that, Includes the following steps: S1. Multi-source data acquisition and temporal feature extraction: Acquire multi-temporal remote sensing image data, environmental auxiliary data, historical land cover data and management data of the study area, preprocess the data and extract temporal feature variables; S2. Construction of Temporal-Causal Fusion Graph Neural Network: Construct a temporal-causal fusion graph neural network, which includes a temporal encoder, a causal graph structure learner, and a physical constraint layer. The temporal encoder encodes temporal features to obtain temporal representation vectors. The causal graph structure learner learns the causal graph structure between environmental factors, temporal feature parameters, and land use labels. The physical constraint layer constrains the causal graph structure based on agronomic and ecological mechanisms. By jointly optimizing prediction accuracy, causal graph topological validity, and physical rationality, the causal graph structure and posterior probability of land use are output. S3. Continuous State Vector Evolution Modeling: Establish a continuous state vector evolution model, using continuous vectors to represent the land use state intensity of spatial units, describing the evolution process of the state vector based on differential equations, the coefficients of the environmental response term of the differential equations are determined according to the cause-effect graph structure, and output the time evolution trajectory of the state vector. S4. Multidimensional physical constraint optimization: The state vector is optimized by introducing spectral constraints, temporal constraints and spatial constraints. The spectral constraints ensure that the state vector is consistent with the observed spectrum, the temporal constraints ensure that the state evolution conforms to the temporal law, and the spatial constraints ensure that the states of adjacent units are continuous. The optimized state vector and land use label are output. S5. Causal intervention strategy optimization: Construct a causal intervention strategy optimization framework based on reinforcement learning, use the causal graph to identify the causal path from the management measure node to the risk node and calculate the path weight, concatenate the path weight with the state vector to form an enhanced state representation, train the strategy network based on the enhanced state representation, and output the optimal management strategy. S6. Change Detection and Time Tracing: Based on the optimized state vector, change detection is performed, and the change time is inferred from the current state through Bayesian inversion. The driving factors are traced using the causal graph. S7. Intelligent Analysis Report Generation: Integrate the analysis results from steps S1 to S6 to generate an intelligent analysis report that includes the current status of land use, evolution trends, risk assessment, and management recommendations. In step S2, the time encoder uses a bidirectional long short-term memory network to encode the time-series feature sequence; The causal graph structure learner maps the temporal representation vector into an adjacency matrix through a neural network. The elements of the adjacency matrix represent the probability that there is a causal edge between the corresponding nodes. The introduction of directed acyclic graph constraints ensures the topological validity of the causal graph. In S2, the physical constraint layer constructs a physical law library, which defines allowed causal relationships and prohibited causal relationships. For each element in the adjacency matrix, its consistency score with the physical law library is calculated. Causal edges that conform to allowed relationships receive positive weights, while causal edges that violate prohibited relationships receive negative penalty weights. The consistency scores are then subjected to nonlinear transformation and regularization to form physical constraint terms, which are jointly optimized with the prediction loss and the directed acyclic graph constraint loss.

2. The method according to claim 1, characterized in that, In S3, the continuous vector is a K-dimensional vector, where K is the number of land types, each component represents the state intensity of the corresponding land type, and the value ranges from zero to one. The sum of all components is equal to one. The differential equation includes an intrinsic growth term, an environmental response term, and a competition transformation term.

3. The method according to claim 2, characterized in that, In S3, the environmental response coefficient is determined as follows: if there is a causal edge from the environmental factor to the land type in the causal graph, the response coefficient of the environmental factor to the land type is non-zero, and its value is determined according to the strength of the causal edge; otherwise, the response coefficient is zero.

4. The method according to claim 1, characterized in that, In S4, the spectral constraint is achieved by establishing a linear mixing relationship between the state vector and the observed spectrum, requiring that the deviation between the pure spectrum of each land type weighted by the state vector and the actual observed spectrum be minimized; the temporal constraint requires that the state evolution trajectory conforms to the temporal change law; the spatial constraint achieves spatial continuity by penalizing the difference in the state vectors of adjacent spatial units.

5. The method according to claim 1, characterized in that, In S5, the causal path weight is calculated as follows: for each directed path from the management measure node to the risk node in the causal graph, the path weight is calculated as the product of the causal strengths of all edges on the path; the enhanced state representation is obtained by concatenating the path weights of each causal path as an additional dimension with the current state vector.

6. The method according to claim 1, characterized in that, In S6, the Bayesian inversion is achieved by constructing a likelihood function and a prior distribution. The likelihood function is based on the evolutionary model, which is integrated forward from the hypothetical change time point to the current time and the degree of agreement with the observed state is calculated. The prior distribution is set based on the statistical regularity of historical changes and the posterior probability distribution of the change time is calculated using the Bayesian formula.

7. An artificial intelligence-based natural resource data analysis system, applied to the method as described in any one of claims 1-6, characterized in that, include: The data acquisition module is used to acquire multi-temporal remote sensing image data, environmental auxiliary data, historical land cover data and management data of the study area, preprocess the data and extract temporal feature variables; The temporal-causal graph neural network module is used to construct a graph neural network that integrates temporal and causal relationships. It jointly optimizes the learning of causal graph structure and posterior probability of land use through a temporal encoder, a causal graph structure learner, and a physical constraint layer. The state evolution modeling module is used to establish a continuous state vector evolution model and calculate the time evolution trajectory of the state vector based on differential equations. The coefficients of the environmental response term in the differential equations are determined by the causal graph structure output by the time-series-causal graph neural network module. The physical constraint optimization module is used to optimize the state vector output by the state evolution modeling module by introducing spectral constraints, temporal constraints and spatial constraints, and output the optimized state vector and land type label; The intervention strategy optimization module is used to construct a causal intervention strategy optimization framework based on reinforcement learning. It uses the causal graph output by the time-series-causal graph neural network module to identify causal paths and calculate path weights. The path weights are concatenated with the state vector to form an enhanced state representation. The strategy network is trained to output the optimal management strategy. The change detection and tracing module is used to detect changes based on the state vector output by the physical constraint optimization module, infer the change time through Bayesian inversion, and identify the driving factors using the causal graph output by the time-series-causal graph neural network module. The report generation module is used to integrate the analysis results from various modules to generate intelligent analysis reports.

Images