Natural resource image data analysis method and system combining spatio-temporal features

Abstract

The application relates to a natural resource image data analysis method and system combining space-time features, which comprises: preprocessing multi-temporal remote sensing images to generate an image data cube; extracting and fusing space-time features from the image data cube to generate fused space-time features; learning the features based on a neural network model and outputting a preliminary land feature category probability distribution; performing differentiable compliance measurement on the preliminary distribution and predefined ecological and geographical knowledge rules to generate a knowledge constraint loss signal; and iteratively optimizing model parameters based on the signal through back propagation, updating the features and the probability distribution, and outputting an analysis graph conforming to the field knowledge constraint until the output conforms to the field knowledge constraint. The method embeds field knowledge in the form of differentiability into the model training and reasoning process, solves the problem that existing deep learning models lack prior constraints in natural resource analysis, and improves the reliability and practicality of automatic analysis results.

Description

Technical Field

[0001] This invention relates to the field of electronic digital data processing technology, specifically to a method and system for analyzing natural resource image data that combines spatiotemporal characteristics. Background Technology

[0002] With the rapid development of remote sensing technology, dynamic monitoring of natural resources using multi-temporal remote sensing imagery has become a core tool for environmental assessment, agricultural surveys, and land use planning. In this field, analytical methods combining spatiotemporal characteristics are receiving increasing attention. These methods aim to go beyond the static interpretation of single images, enhancing the depth and accuracy of understanding natural resource types, states, and their changing mechanisms by comprehensively analyzing the spatial distribution patterns and continuous temporal changes of ground features (such as vegetation phenology, water area fluctuations, and urban expansion). Such methods typically involve preprocessing long-term imagery, extracting and fusing spatiotemporal features, and utilizing complex computational models to uncover hidden patterns.

[0003] However, existing deep learning-based spatiotemporal feature analysis methods have a significant drawback: their decision-making process is often a difficult-to-understand "black box." While these models excel at data fitting, they lack the ability to integrate and adhere to prior domain knowledge. The models may make judgments purely based on statistical regularities in the data, resulting in analytical results that clearly violate basic ecological and geographical principles. For example, they might misidentify a snow-covered mountaintop as "evergreen forest" or mistakenly identify a "water body" on a steep ridge where water is impossible. This lack of interpretability and controllability severely undermines the credibility and application value of the analytical results in serious scientific research and precise industry management. Summary of the Invention

[0004] Based on this, the purpose of this invention is to provide a method and system for analyzing natural resource image data that can effectively embed domain knowledge and enable the analysis process and results to have both high accuracy and high reliability by combining spatiotemporal characteristics.

[0005] The objective of this invention is achieved through the following solution:

[0006] In a first aspect, the present invention provides a method for analyzing natural resource image data that combines spatiotemporal characteristics, comprising the following steps:

[0007] S1: Preprocess the acquired multi-temporal remote sensing images of natural resources to eliminate differences in sensor, atmospheric and geometric positions and generate image data cubes with consistent spatiotemporal references.

[0008] S2: Extract and fuse spatiotemporal features of image data cubes. Based on the preset multi-scale convolution and temporal attention model, mine the spectral, texture and phenological change information of the data at different time phases and perform joint characterization to generate fused spatiotemporal features that have both spatial details and temporal evolution laws.

[0009] S3: Based on a pre-set neural network model, learn the high-order nonlinear relationship between the spatial details and temporal evolution of land features in the spatiotemporal features, and output a preliminary probability distribution of land feature categories;

[0010] S4: Based on the preliminary probability distribution of land cover categories, measure the conformity of domain knowledge rules, quantify the predefined ecological and geographical rules into differentiable loss terms, assess the degree of conflict between the probability distribution of land cover categories and domain priors, and generate knowledge constraint loss signals.

[0011] S5: Based on the backpropagation algorithm, the knowledge constraint loss signal is fed back to the neural network model to update the internal weight parameters of the multi-scale convolution and temporal attention model and the neural network model. Based on the updated multi-scale convolution and temporal attention model and the neural network model, the adjusted fused spatiotemporal features and land cover category probability distribution are recalculated. The optimization is iteratively performed until convergence, and an interpretable spatiotemporal change analysis map of natural resources that conforms to the domain knowledge constraints is output.

[0012] In one embodiment, S1 of the natural resource image data analysis method combining spatiotemporal features provided by the present invention specifically includes the following steps:

[0013] S11: Perform joint geometric fine correction processing on multi-source images from different platforms or sensors. By matching corresponding feature points across images and using regional network adjustment techniques, the geographic coordinate framework of all temporal images is unified, generating a geometrically fine-corrected image set with alignment accuracy.

[0014] S12: Perform atmospheric radiative transfer and topographic illumination correction processing on the geometrically finely corrected image set. Based on the atmospheric radiative transfer model and the topographic radiative correction model, eliminate atmospheric scattering and absorption and the topographic shadowing effect caused by slope aspect and slope on a pixel-by-pixel basis to generate a radiation-normalized image with the physical meaning of the true surface reflectance.

[0015] S13: Perform temporal noise filtering and missing data reconstruction on radiometrically normalized images. Combine the continuity characteristics of time series and the adaptive weighted interpolation model based on local spatiotemporal similarity to detect and remove noisy pixels such as clouds, fog and shadows. Then, use the information of its neighboring effective pixels in the temporal and spatial dimensions to perform weighted interpolation and generate a complete and clean temporal image data cube.

[0016] In one embodiment, step S2 of the natural resource image data analysis method combining spatiotemporal features provided by the present invention specifically includes the following steps:

[0017] S21: Extract time-series phenological parameters from the image data cube, decompose the spectral time-series curve of each pixel within a preset time window based on the harmonic model, extract the start and end times of the growing season, peak amplitude and key points of growth rate, and generate a pixel-level phenological feature parameter set.

[0018] S22: Perform multi-scale spatial-spectral joint feature extraction processing on the image data cube, and use a multi-scale convolutional network to aggregate the spectral and texture context information of the pixel neighborhood at different spatial scales to generate a multi-scale spatial-spectral feature set that can characterize the internal structure and spatial distribution pattern of ground objects.

[0019] S23: Adaptive spatiotemporal feature fusion processing is performed on the phenological feature parameter set and the multi-scale spatial-spectral feature set. The importance weights of different features in the time and space dimensions are dynamically evaluated through the attention mechanism and weighted fusion is performed to generate fused spatiotemporal features.

[0020] In one embodiment, step S3 of the natural resource image data analysis method combining spatiotemporal features provided by the present invention specifically includes the following steps:

[0021] S31: Model the long-term and short-term temporal dependencies of the fused spatiotemporal features, capture the forward evolution and backward dependencies of the land cover state in time through a bidirectional gated cyclic unit network, and generate advanced temporal semantic encoding;

[0022] S32: Perform multi-receptive-field spatial context aggregation processing on high-level temporal semantic coding, and use the spatial pyramid pooling module to aggregate feature context information of different spatial ranges in the neighborhood of a pixel to generate context-enhanced spatial features.

[0023] S33: Perform category probability decoding on the context-enhanced spatial features. Map the high-dimensional features to the same dimension as the number of land cover categories through a fully connected layer, and use the Softmax function to normalize and calculate the probability of belonging to each category, generating a preliminary land cover category probability distribution.

[0024] In one embodiment, step S4 of the natural resource image data analysis method combining spatiotemporal features provided by the present invention specifically includes the following steps:

[0025] S41: Perform phenological period coincidence rule measurement on the preliminary probability distribution of land cover categories, compare the prediction results of high probability land cover categories with the legal phenological period calendar of high probability land cover categories obtained from the statistical analysis of pre-stored historical observation data, calculate the degree of conflict of the prediction results appearing in illegal phenological periods, and generate phenological rule loss components.

[0026] S42: Perform spatial distribution compatibility rule measurement on the preliminary probability distribution of land cover categories, check whether the combination of land cover categories between adjacent pixels violates the laws of ecological symbiosis or geographical distribution based on the prior knowledge of the geographic information system, quantify the degree of conflict of unreasonable spatial adjacency relationships, and generate spatial rule loss components.

[0027] S43: Based on the preset importance coefficients of each rule, the loss components of phenological rules and spatial rules are weighted and fused to generate a unified knowledge constraint loss signal.

[0028] In one embodiment, the formula for calculating the knowledge constraint loss signal of a natural resource image data analysis method combining spatiotemporal features provided by the present invention is as follows:

[0029]

[0030] in, This is a signal of knowledge constraint loss. These are the probability values ​​extracted from the preliminary probability distribution of land cover categories, indicating that a pixel is predicted to be category c. Let be the phenological conflict function, representing the degree of deviation between the observed phenological characteristics of the current pixel and the theoretical phenological calendar of category c. Let be the spatial neighborhood conflict function, representing the degree of spatial topological inconsistency between the predicted class of the current pixel and its spatial neighbor pixel n when the current pixel is predicted as class c. , These are the preset trainable weight coefficients.

[0031] In one embodiment, step S5 of the natural resource image data analysis method combining spatiotemporal features provided by the present invention specifically includes the following steps:

[0032] S51: Calculate and update the gradient of model parameters based on the knowledge-constrained loss signal. Calculate the gradient of the loss signal with respect to the trainable parameters of each layer in the neural network model through the backpropagation algorithm. Update the parameters along the negative gradient direction using the gradient descent algorithm to generate a primary neural network model with optimized parameters.

[0033] S52: Call the optimized primary neural network model to perform forward propagation recalculation on the fused spatiotemporal features, input the features back into the model, obtain the updated land cover category probability distribution guided by knowledge constraints, and generate the knowledge-guided intermediate probability distribution result.

[0034] S53: Iteratively optimize and process the intermediate probability distribution results guided by knowledge until the loss signal converges or the preset number of iterations is reached. Then, convert the final stable probability distribution results into hard classification labels and render them into thematic maps to generate spatiotemporal change analysis maps of natural resources that conform to domain knowledge constraints.

[0035] Secondly, the present invention provides a natural resource image data analysis system that combines spatiotemporal features, the system being configured with the following modules:

[0036] The remote sensing image preprocessing module is used to preprocess the acquired multi-temporal remote sensing images of natural resources, eliminate differences in sensor, atmospheric and geometric positions, and generate image data cubes with consistent spatiotemporal references.

[0037] The spatiotemporal feature fusion module is used to extract and fuse spatiotemporal features of image data cubes. Based on the preset multi-scale convolution and temporal attention model, it mines the spectral, texture and phenological change information of the data at different time phases and performs joint characterization to generate fused spatiotemporal features that have both spatial details and temporal evolution laws.

[0038] The ground feature probability prediction module is used to learn the high-order nonlinear relationship between the spatial details and temporal evolution of ground features in the spatiotemporal features based on a preset neural network model, and output a preliminary probability distribution of ground feature categories.

[0039] The knowledge constraint loss generation module is used to measure the conformity of domain knowledge rules based on the preliminary probability distribution of land cover categories. It quantifies the predefined ecological and geographical rules into differentiable loss terms, assesses the degree of conflict between the probability distribution of land cover categories and domain priors, and generates knowledge constraint loss signals.

[0040] The model iterative optimization output module is used to backpropagate the knowledge constraint loss signal to the neural network model based on the backpropagation algorithm, update the internal weight parameters of the multi-scale convolution and temporal attention model and the neural network model, recalculate the adjusted fused spatiotemporal features and land cover category probability distribution based on the updated multi-scale convolution and temporal attention model and the neural network model, iterate and optimize until convergence, and output an interpretable spatiotemporal change analysis map of natural resources that conforms to the domain knowledge constraints.

[0041] Thirdly, this application provides a computer device, including a memory and a processor, wherein the memory stores a computer program, and the processor executes the computer program to implement any of the above-mentioned methods for analyzing natural resource image data combining spatiotemporal characteristics.

[0042] Fourthly, this application provides a computer-readable storage medium storing a computer program thereon, which, when executed by a processor, implements any of the above-mentioned methods for analyzing natural resource image data incorporating spatiotemporal features.

[0043] In summary, the natural resource image data analysis method provided in this application, which combines spatiotemporal features, can construct a high-quality data foundation with unified spatiotemporal benchmarks for subsequent analysis through rigorous preprocessing of multi-temporal images, thereby eliminating noise and ensuring temporal comparability. Utilizing the synergy of multi-scale convolution and temporal attention mechanisms, it can achieve efficient and deep joint representation of the complex spatial texture and dynamic temporal evolution of ground features, thus capturing subtle change patterns that are difficult to extract using traditional methods. Furthermore, a mechanism is designed to automatically measure the domain knowledge compliance of the preliminary classification results. By transforming ecological and geographical rules into differentiable loss functions, it can quantitatively and in real-time assess the degree of conflict between the model output and physical common sense and professional logic, and generate accurate correction signals. Through feedback optimization based on backpropagation, this signal can be effectively fed back to the feature extraction and classification network to dynamically adjust model parameters, driving the entire system to iteratively evolve in a direction that both fits the image data and strictly adheres to domain rules. This process enables the model to learn autonomously and produce predictions that are not only statistically accurate but also reasonable and interpretable from a professional perspective. For example, it can effectively avoid common-sense errors such as misclassifying snow on a mountaintop as a permanent glacier or misclassifying a seasonal pond as a constant body of water.

[0044] To better understand and implement this invention, the following detailed description is provided in conjunction with the accompanying drawings. Attached Figure Description

[0045] Figure 1 A flowchart illustrating a method for analyzing natural resource image data that combines spatiotemporal features, provided as an embodiment of this application;

[0046] Figure 2 A schematic diagram illustrating the process of generating a unified knowledge constraint loss signal provided in an embodiment of this application;

[0047] Figure 3 This is a schematic diagram of a natural resource image data analysis device that combines spatiotemporal features, provided as another embodiment of this application. Detailed Implementation

[0048] To facilitate understanding of the present invention, a more complete description will be given below with reference to the accompanying drawings. Preferred embodiments of the invention are shown in the drawings. However, the invention can be implemented in many different forms and is not limited to the embodiments described herein. Rather, these embodiments are provided to provide a thorough and complete understanding of the disclosure of the invention.

[0049] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention pertains. The terminology used herein in the description of the invention is for the purpose of describing particular embodiments only and is not intended to be limiting of the invention. The term "and / or" as used herein includes any and all combinations of one or more of the associated listed items.

[0050] In one embodiment, such as Figure 1 As shown, a method for analyzing natural resource imagery data combining spatiotemporal features is provided. This embodiment illustrates the application of this method to a terminal. It is understood that this method can also be applied to a server, and further to a system including both a terminal and a server, and is implemented through interaction between the terminal and the server. In this embodiment, the method includes the following steps:

[0051] S1: Preprocess the acquired multi-temporal remote sensing images of natural resources to eliminate differences in sensor, atmospheric and geometric positions, and generate image data cubes with consistent spatiotemporal references.

[0052] Specifically, the system receives and acquires multi-temporal remote sensing images of natural resources, eliminating differences in sensor response, atmospheric transmission, and geometric positioning to construct a standardized data carrier with consistent spatiotemporal references. When performing sensor radiometric correction, the system first obtains correlation coefficients based on the radiometric calibration file provided by the satellite launcher, converting the digital quantization values ​​of the original images into apparent radiance. Then, using a correction model based on the blackbody radiation law, it performs response consistency calibration on images from different sources, ensuring that comparable radiance data for the same ground feature are generated in images acquired by different sensors.

[0053] Furthermore, the system incorporates a radiative transfer model, taking atmospheric parameters from the image formation process as input. These parameters are derived from corresponding atmospheric products or ground-based measured data. By simulating the atmospheric scattering and absorption of solar radiation, the system retrieves the true reflectance of ground features. For cloud-covered areas in the image, the system employs a spatiotemporal adaptive interpolation method. Based on the spectral-temporal characteristics of cloudless pixels surrounding the cloud area, a relevant regression model is used to fill in the reflectance data of the cloud-masked region, ensuring the continuity of the time series.

[0054] During geometric correction, the system uses a unified geodetic coordinate system as a reference, employs a rational function model, and selects stable features such as road intersections, bridges, and mountaintops as ground control points. The unknown parameters of the model are solved using the least squares method to correct image geometric distortion. Then, resampling technology is used to unify all temporal images to the same pixel resolution and projected coordinate system, achieving spatial reference consistency. After preprocessing, the system organizes the multi-temporal images according to a three-dimensional spatial-spectral-temporal structure, forming an image data cube. The spatial dimension corresponds to the row and column indices of the images, the spectral dimension contains multi-band information of the images, and the temporal dimension corresponds to the multi-temporal sequence. For missing temporal phases, the system uses linear interpolation to supplement them, ensuring uniform time series intervals.

[0055] S2: Extract and fuse spatiotemporal features of image data cubes. Based on the preset multi-scale convolution and temporal attention model, mine the spectral, texture and phenological change information of the data in different time phases and perform joint characterization to generate fused spatiotemporal features that have both spatial details and temporal evolution laws.

[0056] Specifically, the system extracts and fuses spatiotemporal features from the generated image data cubes. Through multi-scale spatial feature mining and temporal attention mechanisms, it achieves deep fusion of spectral, texture, and phenological information. During spatial feature extraction, the system employs dilated convolution configurations with different kernel sizes and dilation rates. Through parallel computation, it captures ground feature features at different spatial levels. Combinations of smaller kernels and corresponding dilation rates are used to extract fine-grained spatial information such as ground feature boundaries and textures. Combinations of medium-sized kernels and corresponding dilation rates are used to capture patch-level features, while combinations of larger kernels and corresponding dilation rates are used to extract regional spatial features.

[0057] After each convolutional operation, the system performs batch normalization and activation function processing to accelerate model convergence and introduce nonlinearity. Then, feature concatenation is used to fuse spatial features from different levels, forming a multi-scale spatial feature map. In the temporal feature extraction stage, the system constructs a temporal processing flow based on a self-attention mechanism. First, the spectral time series of each pixel is compressed to a set dimension through one-dimensional convolution. Then, temporal sequence information is injected through positional encoding. Subsequently, a multi-head self-attention mechanism is used to calculate the correlation weights between each temporal phase and other temporal phases. Temporal correlation features are obtained by parallel computation and concatenation using multiple attention heads. A recurrent unit is used to perform sequence modeling on the temporal correlation features, capturing long-term dependencies in the time series and outputting the temporal feature vector for each pixel.

[0058] In the spatiotemporal feature fusion stage, the system adopts an attention-guided fusion strategy. First, it unifies the dimensions of spatial features and temporal features through one-dimensional convolution. Then, it calculates the correlation matrix of the two types of features to obtain attention weights that reflect their mutual contribution. Feature fusion is achieved through weighted summation. The fused spatiotemporal features contain both spatial details of ground features and temporal evolution patterns.

[0059] S3: Based on a pre-defined neural network model, learns the high-order nonlinear relationship between the spatial details and temporal evolution of land features in the spatiotemporal features, and outputs a preliminary probability distribution of land feature categories.

[0060] Specifically, the system can employ an improved network structure as the core classification network. This network balances feature extraction and spatial localization accuracy through an encoder-decoder structure, making it suitable for semantic segmentation tasks of remote sensing images. The encoder consists of multiple convolutional blocks, each containing consecutive convolutional layers, batch normalization layers, and activation functions. Pooling layers are used to downsample between convolutional blocks, progressively compressing the fused spatiotemporal features into low-resolution, high-dimensional feature maps, thereby capturing the high-order nonlinear relationship between spatial details and temporal evolution. The decoder consists of multiple deconvolutional blocks, each upsampled through deconvolutional layers and connected to the corresponding feature maps of the encoder layers via skip connections to supplement spatial detail information. A bottleneck layer is set between the encoder and decoder, containing multiple convolutional layers and a global average pooling layer to enhance global feature extraction capabilities. The network output layer maps the decoder output features to the dimension corresponding to the number of land cover categories through one-dimensional convolution, and then converts the output into a probability distribution through an activation function.

[0061] During model initialization, the system can initialize the weight parameters of convolutional and fully connected layers using a specific distribution, and initialize the bias parameters to zero to ensure gradient stability in the early stages of training. The system selects the appropriate optimizer, sets the initial learning rate and decay strategy, and uses cross-entropy loss as the initial training loss. During training, the system inputs fused spatiotemporal features into the network according to a set batch size. After feature transformation and nonlinear mapping by the encoder-decoder, it outputs a probability distribution matrix of land cover categories for each pixel. This probability distribution matrix is ​​generated solely based on data-driven statistical laws and does not incorporate domain prior knowledge. It may contain predictions that conflict with basic ecological and geographical common sense, requiring correction through subsequent domain knowledge constraints. The probability distribution matrix is ​​stored in tensor form, with the probability values ​​for each category of each pixel falling within a set range and satisfying the normalization condition that the sum of the probabilities of a single pixel is one.

[0062] S4: Based on the preliminary probability distribution of land cover categories, measure the conformity of domain knowledge rules, quantify the predefined ecological and geographical rules into differentiable loss terms, assess the degree of conflict between the probability distribution of land cover categories and domain priors, and generate knowledge constraint loss signals.

[0063] Specifically, the system predefines three core domain knowledge rules, covering three dimensions: topographic constraints, phenological matching, and ecological association. All rules are based on authoritative geographical and ecological research findings and industry standards. Topographic constraint rules clarify the adaptation relationship between land cover categories and topographic factors; phenological matching rules establish the necessary correlation between land cover categories and time-series spectral characteristics; and ecological association rules constrain the spatial distribution correlation between land cover categories. A quantitative model is designed for each type of rule, transforming qualitative descriptions into computable mathematical expressions. During topographic constraint quantification, the system extracts topographic factors such as slope and altitude based on a digital elevation model, constructs a topographic adaptation function, and the function's output reflects the degree of topographic adaptation. During phenological matching quantification, the system calculates phenological indicators for the time-series spectral characteristics of each pixel, including peak time, valley value, and annual range, constructs a phenological matching function, and the function's output reflects the degree of phenological matching.

[0064] In the ecological association quantification stage, the system can use spatial neighborhood analysis to calculate the association degree between land cover categories, constructing an ecological association function. The function input is the category probability distribution of neighboring pixels, and the output value reflects the degree of ecological association conformity. Based on the quantified rule adaptation function, the system calculates the degree of conflict between the preliminary probability distribution and domain knowledge, constructing a comprehensive knowledge constraint loss term. The terrain constraint loss is obtained by summing the logarithmic products of pixel category probabilities and the output values ​​of the terrain adaptation function; the phenological matching loss is obtained by summing the logarithmic products of pixel category probabilities and the output values ​​of the phenological matching function; and the ecological association loss is obtained by summing the logarithmic products of pixel category probabilities and the output values ​​of the ecological association function. The system weights and sums the three types of losses according to set weight coefficients to obtain the comprehensive knowledge constraint loss term. The system integrates the comprehensive knowledge constraint loss term with the cross-entropy loss to obtain the total loss function. The total loss function includes both the constraint of data fitting accuracy and the rationality requirement of domain knowledge, forming a knowledge constraint loss signal to guide model optimization.

[0065] S5: Based on the backpropagation algorithm, the knowledge constraint loss signal is fed back to the neural network model to update the internal weight parameters of the multi-scale convolution and temporal attention model and the neural network model. Based on the updated multi-scale convolution and temporal attention model and the neural network model, the adjusted fused spatiotemporal features and land cover category probability distribution are recalculated. The optimization is iteratively performed until convergence, and an interpretable spatiotemporal change analysis map of natural resources that conforms to the domain knowledge constraints is output.

[0066] Specifically, based on the chain rule, the system calculates the gradient of the total loss function with respect to the weight parameters of each layer of the model, and the gradient of the total loss with respect to the output layer's logit value. This gradient is obtained by combining the gradient related to data fitting and the gradient related to knowledge constraints. The gradient related to knowledge constraints is obtained by differentiating the losses related to terrain, phenology, and ecology separately and then summing them according to their weights. Subsequently, the system backpropagates the gradients from the output layer to the decoder, bottleneck layer, and encoder, sequentially calculating the weight gradients and bias gradients of the deconvolution layer and convolution layer. The gradients are further propagated to the relevant computational units in the spatiotemporal feature extraction stage, ensuring that all parameters related to feature extraction and classification are adjusted according to knowledge constraints.

[0067] During the parameter update phase, the system uses an optimizer to update all trainable parameters of the model based on the calculated gradients. Weights and biases are adjusted according to their corresponding update formulas. During the update process, spatial feature extraction parameters are adjusted to enhance terrain-adaptive spatial features, while temporal feature extraction parameters prioritize key time phases that match phenological rules. Classification network parameters reduce their reliance on statistical features that violate domain knowledge. The system iteratively optimizes the process of spatiotemporal feature extraction, probability distribution generation, loss calculation, gradient propagation, and parameter updates, reprocessing data based on the updated model in each iteration.

[0068] During the convergence judgment phase, the system sets dual convergence criteria: first, the change in the total loss function over multiple consecutive iterations reaches a set threshold; second, the knowledge constraint loss reaches a set range. Meeting either criterion constitutes convergence. If the number of iterations reaches a preset maximum but the convergence criteria are not met, the system automatically terminates the iteration and outputs the current optimal model result. After iterative convergence, the system outputs the final land cover category probability distribution based on the updated model, uses the maximum probability method to determine the final category of each pixel, and generates a natural resource spatiotemporal change analysis map. The analysis map is stored in a standard geographic information format, containing metadata such as spatial reference information, pixel resolution, and time series identifiers. Output content includes a single-phase land cover classification map, a spatiotemporal change detection map, and optionally a category probability heatmap. The output results strictly adhere to domain knowledge rules, and the classification results can be backtested using quantitative indicators such as terrain adaptability and phenological matching degree.

[0069] In summary, the natural resource image data analysis method provided in this application, which combines spatiotemporal features, can construct a high-quality data foundation with unified spatiotemporal benchmarks for subsequent analysis through rigorous preprocessing of multi-temporal images, thereby eliminating noise and ensuring temporal comparability. Utilizing the synergy of multi-scale convolution and temporal attention mechanisms, it can achieve efficient and deep joint representation of the complex spatial texture and dynamic temporal evolution of ground features, thus capturing subtle change patterns that are difficult to extract using traditional methods. Furthermore, a mechanism is designed to automatically measure the domain knowledge compliance of the preliminary classification results. By transforming ecological and geographical rules into differentiable loss functions, it can quantitatively and in real-time assess the degree of conflict between the model output and physical common sense and professional logic, and generate accurate correction signals. Through feedback optimization based on backpropagation, this signal can be effectively fed back to the feature extraction and classification network to dynamically adjust model parameters, driving the entire system to iteratively evolve in a direction that both fits the image data and strictly adheres to domain rules. This process enables the model to learn autonomously and produce predictions that are not only statistically accurate but also reasonable and interpretable from a professional perspective. For example, it can effectively avoid common-sense errors such as misclassifying snow on a mountaintop as a permanent glacier or misclassifying a seasonal pond as a constant body of water.

[0070] The method provided in this application systematically solves the problem of "black box" decision-making and unreliable results caused by the lack of domain knowledge guidance in the analysis of remote sensing images of natural resources by deep learning models in the existing technology. It constructs a closed-loop framework that deeply embeds domain prior knowledge in a differentiable form into the entire process of model training and inference, which enables the final analysis results to have both high accuracy driven by data and high credibility guided by knowledge. This significantly improves the professional credibility and practical value of natural resource classification and change detection results, while enhancing the acceptability and interpretability of artificial intelligence models in professional fields. It provides an innovative technical approach for high-precision and high-reliability automated natural resource supervision and assessment.

[0071] In one embodiment, S1 of the natural resource image data analysis method combining spatiotemporal features provided by the present invention specifically includes the following steps:

[0072] S11: Perform joint geometric fine correction processing on multi-source images from different platforms or sensors. By matching corresponding feature points across images and using regional network adjustment techniques, the geographic coordinate framework of all temporal images is unified, generating a geometrically fine-corrected image set with alignment accuracy.

[0073] Specifically, the system receives multi-source, multi-temporal remote sensing images of natural resources from different platforms or sensors. Through cross-image homonymous feature point matching and regional network adjustment techniques, it unifies the geographic coordinate framework of all temporal images, generating a geometrically calibrated image set with compliant alignment accuracy. Preferably, the system extracts feature points from all input images, selecting stable ground features in the images as the source of homonymous feature points. These ground features have consistent spatial location and morphological characteristics across different platforms and temporal images. After feature point extraction, the system performs cross-image homonymous feature point matching. By calculating the descriptor similarity of feature points, it establishes the correspondence between feature points in different images, eliminating erroneous correspondences generated during the matching process, and ensuring the accuracy of homonymous feature point matching.

[0074] Furthermore, the system constructs a regional network adjustment model, incorporating all imagery into a unified adjustment calculation framework. Using a pre-defined geographic coordinate system as a reference, the system solves for the exterior orientation elements and ground coordinates of feature points in the imagery using the least squares method. This eliminates local and systematic errors present in the geometric correction of individual images, achieving uniformity across all temporal images within the geographic coordinate framework. After the adjustment calculation is completed, the system performs geometric distortion correction on each image, mapping image pixels to a unified coordinate grid using resampling technology. This ensures complete spatial alignment of images from different platforms and sensors, generating a geometrically precisely corrected image set.

[0075] S12: Perform atmospheric radiative transfer and topographic illumination correction processing on the geometrically finely corrected image set. Based on the atmospheric radiative transfer model and the topographic radiative correction model, eliminate atmospheric scattering and absorption and topographic shadowing effects caused by slope aspect and slope on a pixel-by-pixel basis to generate a radiation-normalized image with the physical meaning of the true surface reflectance.

[0076] Specifically, the system performs atmospheric radiative transfer and topographic illumination correction processing on the geometrically fine-calibrated image set. Through the synergistic effect of two models, it eliminates interference factors pixel-by-pixel, generating a radiation-normalized image with the true physical meaning of surface reflectance. The system calls the atmospheric radiative transfer model, inputting the atmospheric parameters used in image imaging. These parameters are obtained through corresponding atmospheric observation products or synchronous ground observation data. The atmospheric radiative transfer model simulates the transmission process of solar radiation in the atmosphere, including the scattering effect of atmospheric molecules and aerosols, and the absorption effect of atmospheric components. The model calculates and separates atmospheric path radiation and ground-reflected radiation, achieving pixel-by-pixel elimination of the effects of atmospheric scattering and absorption.

[0077] Building upon atmospheric radiometric correction, the system utilizes a topographic radiometric correction model. Based on topographic data provided by the digital elevation model, it extracts the aspect and slope information for each pixel. The model calculates the impact of topographic shadows on pixel radiance based on the geometric relationship between the pixel's aspect, slope, solar altitude angle, and azimuth angle. By adjusting the radiance values, it eliminates radiometric distortion caused by topographic shadows, ensuring that pixel radiance under different topographic conditions, such as slopes and shaded areas, accurately reflects the reflectivity of the ground features. The system employs a pixel-by-pixel processing mode, fusing the calculation results from the atmospheric radiative transfer model and the topographic radiometric correction model. It standardizes and corrects the radiance value of each pixel, outputting a radiometrically standardized image. This ensures that the radiance values ​​of images from different time phases and topographic regions have a unified physical meaning and comparability.

[0078] S13: Perform temporal noise filtering and missing data reconstruction on radiometrically normalized images. Combine the continuity characteristics of time series and the adaptive weighted interpolation model based on local spatiotemporal similarity to detect and remove noisy pixels such as clouds, fog and shadows. Then, use the information of its neighboring effective pixels in the temporal and spatial dimensions to perform weighted interpolation and generate a complete and clean temporal image data cube.

[0079] Specifically, the system performs temporal noise filtering and missing data reconstruction on radiometrically normalized images. Combining the continuity characteristics of the time series with an adaptive weighted interpolation model based on local spatiotemporal similarity, it generates a complete and clean temporal image data cube. The system also detects noisy pixels in the radiometrically normalized images, identifying pixels that deviate from normal temporal variation patterns based on the continuity characteristics of the time series. These pixels are mainly caused by interference factors such as clouds, fog, and shadows. Simultaneously, the system uses spectral feature analysis to assist in the identification of noisy pixels, incorporating pixels with typical noise spectral characteristics into the noise set, thus achieving accurate identification and labeling of noisy pixels.

[0080] After removing noisy pixels, the system initiates an adaptive weighted interpolation model based on local spatiotemporal similarity for missing data regions in the image. This model first searches for neighboring valid pixels in the spatial dimension, which share similar land cover types and spatial association characteristics with the missing pixel. Simultaneously, it mines valid pixels from the same or adjacent time periods in the temporal dimension, utilizing the continuity of the time series to obtain temporal variation reference information for the missing pixel. The system assigns weights to valid pixels in the spatial and temporal dimensions, with weight values ​​determined based on the degree of local spatiotemporal similarity; higher similarity results in a larger weight. The system then calculates the imputation value for the missing pixel through weighted summation.

[0081] During the interpolation process, the system performs pixel-by-pixel completion for each missing region, ensuring that the completed data conforms to the continuity of the time series and the consistency of spatial distribution. After all noisy pixels are removed and missing data is completed, the system organizes the multi-temporal radiometric normalized images according to the spatial-spectral-temporal three-dimensional dimensions to form a complete and clean time-series image data cube.

[0082] In one embodiment, step S2 of the natural resource image data analysis method combining spatiotemporal features provided by the present invention specifically includes the following steps:

[0083] S21: Extract time-series phenological parameters from the image data cube, decompose the spectral time-series curve of each pixel within a preset time window based on the harmonic model, extract the start and end times of the growing season, peak amplitude and key points of growth rate, and generate a pixel-level phenological feature parameter set.

[0084] Specifically, the system extracts temporal phenological parameters from the image data cube, decomposes the spectral temporal curve using a harmonic model, and obtains a pixel-level set of phenological feature parameters. The system defines a preset time window for each pixel in the image data cube, covering the temporal image segment corresponding to the complete growth cycle of the land cover, ensuring that the complete spectral change process from the start to the end of growth can be captured. Further, the system calls the harmonic model, using the spectral temporal curve of each pixel within the time window as model input. The model decomposes the spectral temporal data, and the decomposition process is based on the principle of trigonometric function superposition, separating the periodic and non-periodic components in the spectral curve, where the periodic component corresponds to the regular changes during the land cover growth process.

[0085] Preferably, through fitting calculations using a harmonic model, the system extracts key phenological parameters from the decomposed periodic components, including the start and end times of the growing season, peak amplitude, and key points of growth rate. The start and end times of the growing season are determined by identifying the starting node where the spectral curve begins to deviate from the baseline and continues to rise, and the ending node where the curve falls back to near the baseline; the peak amplitude is obtained by calculating the difference between the peak value of the spectral curve and the baseline; the key points of growth rate are obtained by solving the extreme points of the first derivative of the spectral curve, reflecting the key stages of rate change during the growth of ground features. After all phenological parameters are extracted, the system integrates the parameters corresponding to each pixel to generate a pixel-level phenological feature parameter set, which is stored in matrix form, consistent with the spatial dimension of the image data cube.

[0086] S22: Perform multi-scale spatial-spectral joint feature extraction processing on the image data cube, and use a multi-scale convolutional network to aggregate the spectral and texture context information of the pixel neighborhood at different spatial scales to generate a multi-scale spatial-spectral feature set that can characterize the internal structure and spatial distribution pattern of ground objects.

[0087] Specifically, the system performs multi-scale spatial-spectral joint feature extraction on the image data cube. It aggregates spectral and texture context information at different spatial scales through a multi-scale convolutional network to generate a multi-scale spatial-spectral feature set. The system constructs a multi-scale convolutional network architecture containing multiple parallel convolutional branches, each corresponding to a different spatial scale. By setting convolutional kernels of different sizes, it achieves coverage of neighborhood information within different ranges. During the operation of each convolutional branch, the system uses the spectral dimension information of the image data cube as a basis, combined with the spatial distribution information of the pixel neighborhood, to perform convolution operations on the spectral data of each pixel and its surrounding neighboring pixels, while simultaneously aggregating texture information within the neighborhood.

[0088] Preferably, the texture information aggregation in this application is achieved by analyzing the gray value variation patterns and spatial arrangement patterns among neighboring pixels, covering features such as texture coarseness and density distribution. Different convolutional branches process neighborhood ranges at different spatial scales. Branches corresponding to smaller convolutional kernels focus on the small neighborhood around the pixel, extracting the spatial-spectral features corresponding to the fine structure inside the ground feature; branches corresponding to larger convolutional kernels cover the large neighborhood around the pixel, capturing the spatial-spectral features related to the overall spatial distribution pattern of the ground feature. After each convolutional branch is completed, the system standardizes the output features to eliminate the dimensional differences between features at different scales. Then, the output features of all branches are concatenated and integrated to form a multi-scale spatial-spectral feature set that can characterize the internal structure and spatial distribution pattern of the ground feature. This feature set corresponds one-to-one with the spatial pixels of the image data cube.

[0089] S23: Adaptive spatiotemporal feature fusion processing is performed on the phenological feature parameter set and the multi-scale spatial-spectral feature set. The importance weights of different features in the time and space dimensions are dynamically evaluated through the attention mechanism and weighted fusion is performed to generate fused spatiotemporal features.

[0090] Specifically, the system performs adaptive spatiotemporal feature fusion processing on the phenological feature parameter set and the multi-scale spatial-spectral feature set. It dynamically evaluates the importance weights of features through an attention mechanism and completes weighted fusion to generate fused spatiotemporal features. Furthermore, the system unifies the dimensions of the phenological feature parameter set and the multi-scale spatial-spectral feature set by mapping both types of features to the same feature dimension through one-dimensional convolution operations, ensuring the feasibility of the fusion process. Further, the system introduces an attention mechanism to construct a spatiotemporal attention evaluation model. This model evaluates the importance of phenological features in the time dimension and multi-scale spatial-spectral features in the spatial dimension. In the time dimension, the attention model analyzes the correlation between each parameter in the phenological feature parameter set and the determination of land cover categories, assigning corresponding weights based on the influence of different phenological stages on the classification results. In the spatial dimension, the attention model evaluates the contribution of features at different spatial scales in the multi-scale spatial-spectral feature set, assigning higher weights to spatial scale features that accurately characterize land cover types, taking into account the actual characteristics of land cover spatial distribution.

[0091] Preferably, the system automatically generates features based on the statistical correlation between features and land cover categories through model calculations, requiring no manual intervention and achieving adaptive adjustment. After the weights are determined, the system performs a weighted summation of the phenological feature parameter set and the multi-scale spatial-spectral feature set according to their corresponding weights, completing the deep fusion of spatiotemporal features. During the fusion process, the system ensures that the phenological patterns in the time dimension and the spatial-spectral features in the spatial dimension complement and synergize with each other. The generated fused spatiotemporal features not only include the temporal variation patterns in the land cover growth process but also cover the spectral and textural features at different spatial scales, comprehensively characterizing the spatiotemporal attributes of land covers.

[0092] In one embodiment, step S3 of the natural resource image data analysis method combining spatiotemporal features provided by the present invention specifically includes the following steps:

[0093] S31: Model the long-term and short-term temporal dependencies of the fused spatiotemporal features, capture the forward evolution and backward dependencies of the land cover state in time through a bidirectional gated cyclic unit network, and generate advanced temporal semantic encoding.

[0094] Specifically, after receiving the fused spatiotemporal features, the system models long-term and short-term temporal dependencies. It captures the temporal evolution patterns and dependencies of land cover states through a bidirectional gated recurrent unit (ROU) network, generating high-level temporal semantic codes. The system reshapes the fused spatiotemporal features along the time dimension, forming continuous sequence data for each spatial pixel. This sequence data encompasses the spatiotemporal attribute changes of land cover within the entire time window. Further, the system inputs the reshaped temporal sequence data into the bidirectional gated recurrent unit (ROU) network, which consists of forward-gated and backward-gated recurrent units operating in parallel. The forward-gated recurrent units traverse the temporal sequence in chronological order, while the backward-gated recurrent units traverse it in reverse chronological order. Both networks dynamically adjust information accumulation and forgetting through gating mechanisms.

[0095] Preferably, the system calculates the gating state and the hidden layer output using the following formula:

[0096]

[0097]

[0098]

[0099]

[0100] in, To update the door status, To reset the door status, In the candidate hidden state, Hide the current state. Input features at time t, The state was hidden in the previous moment. , , , , , These are the network weight parameters. It is the sigmoid activation function. The hyperbolic tangent activation function is used. This involves element-wise multiplication. The system concatenates and integrates the forward and backward hidden states, maps them to a preset dimension through a linear transformation, eliminates redundant dimensional interference, and finally generates an advanced temporal semantic code, which corresponds one-to-one with spatial pixels.

[0101] S32: Perform multi-receptive-field spatial context aggregation processing on high-level temporal semantic coding, and use the spatial pyramid pooling module to aggregate feature context information of different spatial ranges within the pixel neighborhood to generate context-enhanced spatial features.

[0102] Specifically,

[0103] S33: Perform category probability decoding on the context-enhanced spatial features. Map the high-dimensional features to the same dimension as the number of land cover categories through a fully connected layer, and use the Softmax function to normalize and calculate the probability of belonging to each category, generating a preliminary land cover category probability distribution.

[0104] Specifically, the system performs multi-receptive-field spatial context aggregation processing on the high-level temporal semantic encoding, using a spatial pyramid pooling mechanism to aggregate feature context information from different spatial ranges, generating context-enhanced spatial features. The system converts the high-level temporal semantic encoding into a feature map format suitable for spatial feature processing, ensuring that the spatial dimension of the feature map remains consistent with the original image data cube, and embedding temporal semantic information into the channel dimension of the feature map. Further, the system performs pooling operations on the feature map using multiple sliding windows with different spatial ranges, each sliding window corresponding to a spatial receptive field, and aggregates the feature values ​​within the window. The system integrates multi-scale pooled features using the following formula:

[0105]

[0106] in, This is the result of multi-scale pooling feature concatenation. These are the pooling features corresponding to different receptive fields. The feature concatenation operation integrates spatial context information at different scales. Pooling operations with smaller windows aggregate local spatial context information of ground features, while pooling operations with larger windows capture global spatial context information of the area where the ground feature is located. After each pooling process, the system performs dimensionality compression and standardization on each pooled feature to eliminate dimensional differences and redundancy between features from different receptive fields. Further feature extraction and optimization are then performed on the concatenated features through convolutional layers to enhance their discriminative ability and generate context-enhanced spatial features. These features contain both temporal semantic information of ground features and multi-scale spatial context information, providing comprehensive feature support for category probability decoding.

[0107] In one embodiment, such as Figure 2 As shown, S4 of the natural resource image data analysis method combining spatiotemporal features provided by the present invention specifically includes the following steps:

[0108] S41: Perform phenological period coincidence rule measurement on the preliminary probability distribution of land cover categories, compare the prediction results of high-probability land cover categories with the legal phenological period calendar of high-probability land cover categories obtained from pre-stored historical observation data, calculate the degree of conflict of the prediction results appearing in illegal phenological periods, and generate phenological rule loss components.

[0109] Specifically, the system initiates phenological period consistency rule measurement processing on the preliminary probability distribution of land cover categories. By comparing the prediction results with the legal phenological period calendar, the degree of conflict is calculated to generate a phenological rule loss component. The system calls a pre-stored high-probability legal phenological period calendar for land cover categories, which is formed based on historical observation data statistics. This calendar covers the legal growth and existence periods of each land cover category in the time dimension.

[0110] Furthermore, the system extracts the high-probability predicted category for each pixel in the preliminary land cover category probability distribution. It then compares the prediction result of this high-probability category with the corresponding legal phenological calendar along the time dimension to determine whether the observation period corresponding to the prediction result falls within the legal phenological period for that category. The system uses a phenological conflict function... Quantitative comparison results This indicates the degree of deviation between the observed phenological characteristics of the current pixel and the theoretical phenological calendar for category c. A greater deviation results in a higher function output value. The system combines the probability value extracted from the preliminary land cover category probability distribution to predict the pixel as category c. ,Will and Perform a product operation to quantify the degree of conflict in the prediction results during illegal phenological periods, and generate the phenological rule loss component corresponding to each pixel.

[0111] S42: Perform spatial distribution compatibility rule measurement on the preliminary probability distribution of land cover categories. Based on the prior knowledge of the geographic information system, check whether the combination of land cover categories between adjacent pixels violates the laws of ecological symbiosis or geographical distribution. Quantify the degree of conflict of unreasonable spatial adjacency relationships and generate spatial rule loss components.

[0112] Specifically, the system performs spatial distribution compatibility rule measurement on the preliminary land cover category probability distribution. By checking the rationality of adjacent cell category combinations, it quantifies the degree of conflict to generate a spatial rule loss component. The system calls upon pre-stored prior knowledge of the geographic information system, which covers the ecological symbiotic relationships and geographical distribution patterns among land cover categories, clarifying the spatially permissible adjacent category combinations for each land cover category. Further, the system extracts the predicted category of each cell in the preliminary land cover category probability distribution, and simultaneously obtains the set of spatially neighboring cells of that cell and its corresponding predicted category, forming adjacent cell category combination pairs. The system checks whether each category combination pair conforms to the pre-stored ecological symbiosis or geographical distribution laws, determining whether there are any unreasonable spatial adjacency relationships.

[0113] Preferably, the system uses a spatial neighborhood conflict function. Quantify the inspection results, This indicates the degree of spatial topological inconsistency between the predicted class of the current pixel and its spatial neighbor pixels n when the current pixel is predicted as class c. A higher degree of inconsistency results in a larger function output value. The system combines the probability value of the pixel being predicted as class c. ,Will and Perform a product operation to quantify the degree of conflict of unreasonable spatial adjacency relationships, and generate a spatial rule loss component corresponding to each pixel. This component is consistent with the spatial dimension of the phenological rule loss component.

[0114] S43: Based on the preset importance coefficients of each rule, the loss components of phenological rules and spatial rules are weighted and fused to generate a unified knowledge constraint loss signal.

[0115] Specifically, based on preset importance coefficients for each rule, the system weights and fuses the loss components of phenological rules and spatial rules to generate a unified knowledge constraint loss signal. The system obtains preset trainable weight coefficients. , This coefficient is used to balance the importance of phenological rule constraints and spatial rule constraints, and can be adaptively adjusted according to actual data characteristics during model training. Preferably, the system calls the phenological rule loss component and the spatial rule loss component, and combines the phenological rule loss component with the weight coefficient. Perform a product operation to obtain the weighted phenological rule loss; then combine the spatial rule loss components with the weighting coefficients. Perform a product operation to obtain the weighted spatial rule loss. Preferably, the formula for calculating the knowledge constraint loss signal is:

[0116]

[0117] in, This is a signal of knowledge constraint loss. These are the probability values ​​extracted from the preliminary probability distribution of land cover categories, indicating that a pixel is predicted to be category c. Let be the phenological conflict function, representing the degree of deviation between the observed phenological characteristics of the current pixel and the theoretical phenological calendar of category c. Let be the spatial neighborhood conflict function, representing the degree of spatial topological inconsistency between the predicted class of the current pixel and its spatial neighbor pixel n when the current pixel is predicted as class c. , These are the preset trainable weight coefficients. The system uses this formula to sum the weighted losses of the two classes and combines them with the probability values. Once the integration is complete, a unified knowledge constraint loss signal is generated, which can comprehensively reflect the degree of conflict between the prediction results and the domain knowledge rules.

[0118] In one embodiment, step S5 of the natural resource image data analysis method combining spatiotemporal features provided by the present invention specifically includes the following steps:

[0119] S51: Calculate and update the gradient of model parameters based on the knowledge-constrained loss signal. Calculate the gradient of the loss signal with respect to the trainable parameters of each layer in the neural network model using the backpropagation algorithm. Update the parameters along the negative gradient direction using the gradient descent algorithm to generate a primary neural network model with optimized parameters.

[0120] Specifically, based on the knowledge constraint loss signal, the system calculates and updates the gradients of the model parameters, generating a primary neural network model with optimized parameters. The system calls the error backpropagation algorithm and performs gradient calculation according to the chain rule. This process starts from the model output layer and propagates the loss signal backward layer by layer to the feature extraction layer, temporal modeling layer, spatial aggregation layer, etc., solving for the partial derivatives of the knowledge constraint loss signal with respect to the trainable parameters of each layer. The gradient calculation follows the formula below:

[0121]

[0122] in, For the loss signal to parameters gradient vector, This represents the probability distribution output by the model. The category mapping values ​​output by the fully connected layer. This encompasses all trainable parameters in the model, including convolutional layer weights, fully connected layer weights, and bias terms. Preferably, after constructing the complete gradient matrix using this formula, the system can perform parameter updates using the gradient descent algorithm, with the update formula being:

[0123]

[0124] in, For the updated parameters, For the current parameter, As the learning rate, the system adjusts all trainable parameters along the negative gradient direction to reduce the value of the loss signal, thereby optimizing the model parameters in a direction that conforms to the domain knowledge constraints and generating a primary neural network model with optimized parameters.

[0125] S52: Call the optimized primary neural network model to perform forward propagation recalculation on the fused spatiotemporal features, input the features back into the model, obtain the updated land cover category probability distribution guided by knowledge constraints, and generate the knowledge-guided intermediate probability distribution result.

[0126] Specifically, the system calls the optimized primary neural network model to perform forward propagation recalculation on the fused spatiotemporal features, generating knowledge-guided intermediate probability distribution results. The system re-inputs the fused spatiotemporal features into the optimized model. The features first enter the feature extraction layer, where spatial and spectral features are extracted again using updated convolution parameters. They are then passed to the temporal modeling layer, where optimized recurrent unit parameters capture long- and short-term temporal dependencies. Next, the features undergo pooling operations in the spatial aggregation layer to integrate multi-scale spatial context information, and finally enter the decoding layer. Preferably, the decoding layer performs feature mapping through a fully connected layer, with the mapping formula as follows:

[0127]

[0128] in, For category mapping values, This is the aggregated high-dimensional feature vector. This is the updated weight matrix of the fully connected layer. This is the updated bias vector. After feature mapping, the system performs normalization calculations using the Softmax function:

[0129]

[0130] in, Let be the probability that a pixel belongs to category c. The mapping value for category c. This represents the total number of land cover categories. The forward propagation process uses updated parameters throughout, correcting feature responses that conflict with domain knowledge, strengthening feature representations that conform to rules, and outputting an updated probability distribution of land cover categories, thus forming a knowledge-guided intermediate probability distribution result.

[0131] S53: Iteratively optimize and process the intermediate probability distribution results guided by knowledge until the loss signal converges or the preset number of iterations is reached. Then, convert the final stable probability distribution results into hard classification labels and render them into thematic maps to generate spatiotemporal change analysis maps of natural resources that conform to domain knowledge constraints.

[0132] Specifically, the system iteratively optimizes and processes the intermediate probability distribution results guided by knowledge until the convergence condition is met or the preset number of iterations is reached. The system feeds back the intermediate probability distribution results to the previous step, recalculates the phenological rule loss component and the spatial rule loss component, and fuses them to obtain the updated knowledge constraint loss signal. Then, the gradient calculation and parameter update in step S51 and the forward propagation recalculation in step S52 are repeated to form a closed-loop iterative cycle. Preferably, the convergence judgment follows the following formula:

[0133]

[0134] in, This is the loss signal for the current iteration. This is the loss signal from the previous iteration. The convergence threshold is used to determine convergence when the difference in the loss signal of consecutive iterations is lower than this threshold. If the number of iterations reaches the preset upper limit, the process will automatically terminate.

[0135] After the iteration terminates, the system determines the hard classification label using the maximum probability decision criterion:

[0136]

[0137] in, For the final hard classification label, This represents the category probability in the final probability distribution. The system performs spatial rendering on the hard classification labels of all pixels, generates thematic maps according to the preset category-color mapping relationship, integrates spatial positioning information and time series identifiers, and generates a natural resource spatiotemporal change analysis map that conforms to domain knowledge constraints, for use in dynamic monitoring and management decision-making of natural resources.

[0138] It should be understood that although the steps in the flowcharts of the embodiments described above are shown sequentially according to the arrows, these steps are not necessarily executed in the order indicated by the arrows. Unless explicitly stated herein, there is no strict order restriction on the execution of these steps, and they can be executed in other orders. Moreover, at least some steps in the flowcharts of the embodiments described above may include multiple steps or multiple stages. These steps or stages are not necessarily completed at the same time, but can be executed at different times. The execution order of these steps or stages is not necessarily sequential, but can be performed alternately or in turn with other steps or at least some of the steps or stages of other steps.

[0139] Based on the same inventive concept, this application also provides a spatiotemporal feature-based natural resource image data analysis device for implementing the above-mentioned spatiotemporal feature-based natural resource image data analysis method. The solution provided by this device is similar to the solution described in the above method. Therefore, the specific limitations of one or more spatiotemporal feature-based natural resource image data analysis device embodiments provided below can be found in the above-described limitations of the spatiotemporal feature-based natural resource image data analysis method, and will not be repeated here.

[0140] In one embodiment, such as Figure 3 As shown, the present invention provides a natural resource image data analysis system 600 that combines spatiotemporal features. This system is configured with the following modules:

[0141] The remote sensing image preprocessing module 610 is used to preprocess the acquired multi-temporal natural resource remote sensing images, eliminate differences in sensor, atmosphere and geometric location, and generate image data cubes with consistent spatiotemporal reference.

[0142] The spatiotemporal feature fusion module 620 is used to extract and fuse spatiotemporal features of image data cubes. Based on the preset multi-scale convolution and temporal attention model, it mines the spectral, texture and phenological change information of the data at different time phases and performs joint characterization to generate fused spatiotemporal features that have both spatial details and temporal evolution laws.

[0143] The ground feature probability prediction module 630 is used to learn the high-order nonlinear relationship between the spatial details and temporal evolution of ground features in the spatiotemporal features based on a preset neural network model, and output a preliminary probability distribution of ground feature categories.

[0144] The knowledge constraint loss generation module 640 is used to measure the conformity of domain knowledge rules based on the preliminary probability distribution of land cover categories, quantify the predefined ecological and geographical rules into differentiable loss terms, evaluate the degree of conflict between the probability distribution of land cover categories and the domain prior, and generate a knowledge constraint loss signal.

[0145] The model iterative optimization output module 650 is used to backpropagate the knowledge constraint loss signal to the neural network model based on the backpropagation algorithm, update the internal weight parameters of the multi-scale convolution and temporal attention model and the neural network model, recalculate the adjusted fused spatiotemporal features and land cover category probability distribution based on the updated multi-scale convolution and temporal attention model and the neural network model, iterate and optimize until convergence, and output an interpretable spatiotemporal change analysis map of natural resources that conforms to the domain knowledge constraints.

[0146] Preferably, the remote sensing image preprocessing module 610 provided in this application is configured with the following units:

[0147] The geometric fine correction processing unit is used to perform joint geometric fine correction processing on multi-source images from different platforms or sensors. It unifies the geographic coordinate framework of all temporal images through cross-image homonymous feature point matching and regional network adjustment technology, and generates a geometric fine correction image set with compliant alignment accuracy.

[0148] The radiation normalization correction unit is used to perform atmospheric radiation transfer and topographic illumination correction processing on the geometrically finely corrected image set. Based on the atmospheric radiation transfer model and the topographic radiation correction model, it eliminates atmospheric scattering and absorption and the topographic shadowing effect caused by slope aspect and slope on a pixel-by-pixel basis, generating radiation normalized images with the physical meaning of the true surface reflectance.

[0149] The temporal image restoration and completion unit is used to perform temporal noise filtering and missing data reconstruction on radiometrically normalized images. It combines the continuity characteristics of time series and an adaptive weighted interpolation model based on local spatiotemporal similarity to detect and remove noisy pixels such as clouds, fog, and shadows. It then uses the information of its neighboring valid pixels in the temporal and spatial dimensions to perform weighted interpolation, generating a complete and clean temporal image data cube.

[0150] Preferably, the spatiotemporal feature fusion module 620 provided in this application is configured with the following units:

[0151] The temporal phenological parameter extraction unit is used to extract temporal phenological parameters from the image data cube. Based on the harmonic model, it decomposes the spectral temporal curve of each pixel within a preset time window, extracts the start and end times of the growing season, peak amplitude and key points of growth rate, and generates a pixel-level phenological feature parameter set.

[0152] The multi-scale spatial-spectral feature extraction unit is used to perform multi-scale spatial-spectral joint feature extraction processing on image data cubes. It uses a multi-scale convolutional network to aggregate spectral and texture context information of pixel neighborhoods at different spatial scales to generate a multi-scale spatial-spectral feature set that can characterize the internal structure and spatial distribution pattern of ground objects.

[0153] The adaptive spatiotemporal feature fusion unit is used to perform adaptive spatiotemporal feature fusion processing on the phenological feature parameter set and the multi-scale spatial-spectral feature set. It dynamically evaluates the importance weights of different features in the time and space dimensions through an attention mechanism and performs weighted fusion to generate fused spatiotemporal features.

[0154] Preferably, the ground feature probability prediction module 630 provided in this application is configured with the following units:

[0155] The temporal dependency modeling unit is used to model the long-term and short-term temporal dependencies of fused spatiotemporal features. It captures the forward evolution and backward dependency of land cover states in time through a bidirectional gated recurrent unit network, and generates advanced temporal semantic encoding.

[0156] The spatial context aggregation unit is used to perform multi-receptive field spatial context aggregation processing on high-level temporal semantic coding. It uses the spatial pyramid pooling module to aggregate feature context information of different spatial ranges within the pixel neighborhood to generate context-enhanced spatial features.

[0157] The category probability decoding unit is used to perform category probability decoding on the context-enhanced spatial features. It maps high-dimensional features to the same dimension as the number of land cover categories through a fully connected layer, and uses the Softmax function to normalize and calculate the probability of belonging to each category, generating a preliminary land cover category probability distribution.

[0158] Preferably, the knowledge constraint loss generation module 640 provided in this application is configured with the following units:

[0159] The phenological rule measurement unit is used to perform phenological period consistency rule measurement processing on the preliminary probability distribution of land cover categories. It compares the prediction results of high-probability land cover categories with the legal phenological period calendar of high-probability land cover categories obtained from pre-stored historical observation data, calculates the degree of conflict of the prediction results appearing in illegal phenological periods, and generates phenological rule loss components.

[0160] The spatial rule measurement unit is used to perform spatial distribution compatibility rule measurement on the preliminary probability distribution of land cover categories. Based on the prior knowledge of the geographic information system, it checks whether the combination of land cover categories between adjacent pixels violates the ecological symbiosis or geographical distribution law, quantifies the degree of conflict of unreasonable spatial adjacency relationships, and generates spatial rule loss components.

[0161] The knowledge loss fusion unit is used to weight and fuse the phenological rule loss component and the spatial rule loss component based on the preset importance coefficients of various rules, so as to generate a unified knowledge constraint loss signal.

[0162] Preferably, the model iterative optimization output module 650 provided in this application is configured with the following units:

[0163] The model parameter gradient update unit is used to calculate and update the model parameter gradient based on the knowledge constraint loss signal. It calculates the gradient of the loss signal with respect to the trainable parameters of each layer in the neural network model through the error backpropagation algorithm, and updates the parameters along the negative gradient direction using the gradient descent algorithm to generate a primary neural network model with optimized parameters.

[0164] The feature forward recalculation unit is used to call the parameter-optimized primary neural network model to perform forward propagation recalculation on the fused spatiotemporal features, input the features back into the model, obtain the updated land cover category probability distribution guided by knowledge constraints, and generate the knowledge-guided intermediate probability distribution result.

[0165] The iterative optimization mapping unit is used to iteratively optimize and process the intermediate probability distribution results guided by knowledge until the loss signal converges or the preset number of iterations is reached. The final stable probability distribution results are converted into hard classification labels and rendered into thematic maps to generate spatiotemporal change analysis maps of natural resources that conform to domain knowledge constraints.

[0166] In one embodiment, this application also provides a computer device, including a memory and a processor, wherein the memory stores a computer program, and the processor executes the computer program to implement the above-described method for analyzing natural resource image data incorporating spatiotemporal features.

[0167] In one embodiment, this application also provides a computer-readable storage medium storing a computer program thereon, which, when executed by a processor, implements the above-described method for analyzing natural resource image data incorporating spatiotemporal features.

[0168] In the description of this specification, references to terms such as "one embodiment," "some embodiments," "example," "specific example," or "some examples," etc., indicate that a specific feature, structure, material, or characteristic described in connection with that embodiment or example is included in at least one embodiment or example of this application. Furthermore, the specific features, structures, materials, or characteristics described may be combined in any suitable manner in one or more embodiments or examples. Moreover, without contradiction, those skilled in the art can combine and integrate the different embodiments or examples described in this specification, as well as the features of those different embodiments or examples.

[0169] For the device embodiments, since they basically correspond to the method embodiments, the relevant parts can be referred to in the description of the method embodiments. The device embodiments described above are merely illustrative. The components described as separate parts may or may not be physically separate, and the components shown as units may or may not be physical units, that is, they may be located in one place or distributed across multiple network units. Some or all of the modules can be selected to achieve the purpose of this disclosure according to actual needs. Those skilled in the art can understand and implement this without creative effort.

[0170] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any person skilled in the art can easily conceive of various variations or substitutions within the technical scope disclosed in this application, and these should all be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

Claims

1. A method for analyzing natural resource image data combining spatiotemporal characteristics, characterized in that, Includes the following steps: S1: Preprocess the acquired multi-temporal remote sensing images of natural resources to eliminate differences in sensor, atmospheric and geometric positions and generate image data cubes with consistent spatiotemporal references. S2: Perform spatiotemporal feature extraction and fusion on the image data cube. Based on the preset multi-scale convolution and temporal attention model, mine its spectral, texture and phenological change information at different time phases and perform joint characterization to generate fused spatiotemporal features that have both spatial details and temporal evolution laws. S3: Based on the preset neural network model, learn the high-order nonlinear relationship between the spatial details and temporal evolution of land features in the fused spatiotemporal features, and output the preliminary probability distribution of land feature categories; S4: Based on the preliminary probability distribution of the land cover categories, measure the conformity of the domain knowledge rules, quantify the predefined ecological and geographical rules into differentiable loss terms, evaluate the degree of conflict between the probability distribution of the land cover categories and the domain prior, and generate a knowledge constraint loss signal. S5: Based on the backpropagation algorithm, the knowledge constraint loss signal is fed back to the neural network model, the internal weight parameters of the multi-scale convolution and temporal attention model and the neural network model are updated, and the adjusted fused spatiotemporal features and land cover category probability distribution are recalculated based on the updated multi-scale convolution and temporal attention model and the neural network model. The optimization is iteratively performed until convergence, and an interpretable spatiotemporal change analysis map of natural resources that conforms to the domain knowledge constraints is output.

2. The method according to claim 1, characterized in that, S1 includes: S11: Perform joint geometric fine correction processing on multi-source images from different platforms or sensors. By matching corresponding feature points across images and using regional network adjustment techniques, the geographic coordinate framework of all temporal images is unified, generating a geometrically fine-corrected image set with alignment accuracy. S12: Perform atmospheric radiative transfer and topographic illumination correction processing on the geometrically corrected image set. Based on the atmospheric radiative transfer model and the topographic radiative correction model, eliminate atmospheric scattering and absorption and the topographic shadowing effect caused by slope aspect and slope pixel by pixel to generate a radiation-normalized image with the physical meaning of the true surface reflectance. S13: Perform temporal noise filtering and missing data reconstruction on the radiometrically normalized image. Combine the continuity characteristics of the time series and the adaptive weighted interpolation model based on local spatiotemporal similarity to detect and remove noisy pixels such as clouds, fog and shadows. Then, use the information of its neighboring effective pixels in the time and space dimensions to perform weighted interpolation to generate a complete and clean temporal image data cube.

3. The method according to claim 1, characterized in that, S2 includes: S21: Extract time-series phenological parameters from the image data cube, decompose the spectral time-series curve of each pixel within a preset time window based on the harmonic model, extract the start and end times of the growing season, peak amplitude and key points of growth rate, and generate a pixel-level phenological feature parameter set. S22: Perform multi-scale spatial-spectral joint feature extraction processing on the image data cube, and use a multi-scale convolutional network to aggregate the spectral and texture context information of the pixel neighborhood at different spatial scales to generate a multi-scale spatial-spectral feature set that can characterize the internal structure and spatial distribution pattern of ground objects. S23: Adaptive spatiotemporal feature fusion processing is performed on the phenological feature parameter set and the multi-scale spatial-spectral feature set. The importance weights of different features in the time and space dimensions are dynamically evaluated through an attention mechanism and weighted fusion is performed to generate fused spatiotemporal features.

4. The method according to claim 1, characterized in that, S3 includes: S31: Model the long-term and short-term temporal dependencies of the fused spatiotemporal features, capture the forward evolution and backward dependency of the land cover state in time through a bidirectional gated cyclic unit network, and generate advanced temporal semantic encoding; S32: Perform multi-receptive-field spatial context aggregation processing on the advanced temporal semantic coding, and use the spatial pyramid pooling module to aggregate feature context information of different spatial ranges in the neighborhood of a pixel to generate context-enhanced spatial features. S33: Perform category probability decoding processing on the context-enhanced spatial features, map the high-dimensional features to the same dimension as the number of land cover categories through a fully connected layer, and use the Softmax function to normalize and calculate the probability of belonging to each category, generating a preliminary land cover category probability distribution.

5. The method according to claim 1, characterized in that, S4 includes: S41: Perform phenological period coincidence rule measurement processing on the preliminary probability distribution of the land cover categories, compare the prediction results of high probability land cover categories with the legal phenological period calendar of high probability land cover categories obtained from the statistical analysis of pre-stored historical observation data, calculate the degree of conflict of the prediction results appearing in illegal phenological periods, and generate phenological rule loss components. S42: Perform spatial distribution compatibility rule measurement processing on the preliminary probability distribution of land cover categories, check whether the combination of land cover categories between adjacent pixels violates the ecological symbiosis or geographical distribution law based on the prior knowledge of the geographic information system, quantify the degree of conflict of unreasonable spatial adjacency relationships, and generate spatial rule loss components. S43: Based on the preset importance coefficients of each rule, the loss components of the phenological rules and the spatial rules are weighted and fused to generate a unified knowledge constraint loss signal.

6. The method according to claim 5, characterized in that, The formula for calculating the knowledge constraint loss signal is as follows: in, This is a signal of knowledge constraint loss. These are the probability values ​​extracted from the preliminary probability distribution of land cover categories, indicating that a pixel is predicted to be category c. Let be the phenological conflict function, representing the degree of deviation between the observed phenological characteristics of the current pixel and the theoretical phenological calendar of category c. Let be the spatial neighborhood conflict function, representing the degree of spatial topological inconsistency between the predicted class of the current pixel and its spatial neighbor pixel n when the current pixel is predicted as class c. , These are the preset trainable weight coefficients.

7. The method according to any one of claims 1-6, characterized in that, S5 includes: S51: Calculate and update the gradient of the model parameters based on the knowledge-constrained loss signal, calculate the gradient of the loss signal with respect to the trainable parameters of each layer in the neural network model through the backpropagation algorithm, update the parameters along the negative gradient direction using the gradient descent algorithm, and generate a primary neural network model with optimized parameters. S52: Call the optimized primary neural network model to perform forward propagation recalculation on the fused spatiotemporal features, input the features back into the model, obtain the updated land cover category probability distribution guided by knowledge constraints, and generate the knowledge-guided intermediate probability distribution result. S53: Iteratively optimize and process the intermediate probability distribution results guided by knowledge until the loss signal converges or the preset number of iterations is reached. Then, convert the final stable probability distribution results into hard classification labels and render them into thematic maps to generate spatiotemporal change analysis maps of natural resources that conform to domain knowledge constraints.

8. A natural resource image data analysis system combining spatiotemporal characteristics, characterized in that, The system includes: The remote sensing image preprocessing module is used to preprocess the acquired multi-temporal remote sensing images of natural resources, eliminate differences in sensor, atmospheric and geometric positions, and generate image data cubes with consistent spatiotemporal references. The spatiotemporal feature fusion module is used to extract and fuse spatiotemporal features of the image data cube. Based on the preset multi-scale convolution and temporal attention model, it mines the spectral, texture and phenological change information of the data at different time phases and performs joint characterization to generate fused spatiotemporal features that have both spatial details and temporal evolution laws. The ground feature probability prediction module is used to learn the high-order nonlinear relationship between the spatial details and temporal evolution of ground features in the fused spatiotemporal features based on a preset neural network model, and output a preliminary probability distribution of ground feature categories. The knowledge constraint loss generation module is used to measure the conformity of domain knowledge rules based on the preliminary probability distribution of the land cover categories, quantify the predefined ecological and geographical rules into differentiable loss terms, evaluate the degree of conflict between the probability distribution of the land cover categories and the domain prior, and generate a knowledge constraint loss signal. The model iterative optimization output module is used to backpropagate the knowledge constraint loss signal to the neural network model based on the backpropagation algorithm, update the internal weight parameters of the multi-scale convolution and temporal attention model and the neural network model, recalculate the adjusted fused spatiotemporal features and land cover category probability distribution based on the updated multi-scale convolution and temporal attention model and the neural network model, iterate and optimize until convergence, and output an interpretable spatiotemporal change analysis map of natural resources that conforms to domain knowledge constraints.

9. A computer device comprising a memory and a processor, wherein the memory stores a computer program, characterized in that, When the processor executes the computer program, it implements the method of any one of claims 1 to 7.

10. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by a processor, it implements the method of any one of claims 1 to 7.

Images