A land and resource utilization evaluation method and system based on multi-source data fusion

Abstract

The present application relates to the technical field of spatial data processing, in particular to a land resource utilization evaluation method and system based on multi-source data fusion, comprising the following steps: acquiring remote sensing image sequences and collecting ground meteorological parameters of the corresponding area through a multispectral camera, calculating cross-modal differences to generate heterogeneous feature tensors, inputting manifold deep network extraction matrices to construct topological network structures, using a dynamic graph convolution model to correct the forward propagation path to extract multi-scale fusion vectors, and measuring gradient numerical density expectation to analytically calculate output comprehensive indexes. In the present application, the feature distribution heterogeneity is eliminated by introducing a spatiotemporal heterogeneous tensor mechanism, a cross-modal nonlinear interaction constraint is established to guide the iterative update of the projection matrix, and the aggregated weight is dynamically allocated to break the representation bias caused by the dependence on static projection and single scale, so that the resource state under complex terrain can be accurately expressed by directly analyzing the index from the aggregation result.

Description

Technical Field

[0001] This invention relates to the field of spatial data processing technology, and in particular to a method and system for land resource utilization assessment based on multi-source data fusion. Background Technology

[0002] The field of spatial data processing technology involves operations such as acquiring, reducing dimensionality, and restructuring cross-modal geographic information and surface parameters to achieve dynamic monitoring of natural resource elements and spatial quantitative mapping of regional development status. Among these, traditional land resource utilization assessment methods based on multi-source data fusion involve merging the multispectral pixel matrix of satellite remote sensing imagery with humidity values ​​output from ground monitoring stations using a matrix stitching method. The global variance of the fusion matrix is ​​then calculated using a static feature projection function to select representative principal components, and land parcel attribute classification results are output based on a decision tree structure.

[0003] Traditional land resource utilization assessment methods based on multi-source data fusion rely on manually set projection rules for feature alignment. Fixed-dimensional transformation methods are difficult to capture nonlinear correlations between cross-modal features when dealing with heterogeneous remote sensing images and ground sensing vector parameters. Direct matrix merging is prone to conflicting feature dimension information. Preset rigid weight allocation mechanisms cannot be adaptively updated according to the changing patterns of land parcels. When faced with complex terrain boundaries and differences in multi-source modal attributes, the extracted vectors have structural biases. Single-scale classification regression cannot characterize the evolution state of land parcels. Summary of the Invention

[0004] The purpose of this invention is to address the shortcomings of existing technologies by proposing a land resource utilization assessment method and system based on multi-source data fusion.

[0005] To achieve the above objectives, the present invention adopts the following technical solution: a land resource utilization assessment method based on multi-source data fusion, comprising the following steps: S1: Acquire remote sensing image sequences through a spaceborne multispectral camera, collect ground meteorological parameters in conjunction with a micro-monitoring station, calculate the difference Euclidean distance between the remote sensing image sequences and the ground meteorological parameters, expand the feature dimension based on the difference Euclidean distance, and construct a heterogeneous feature tensor. S2: Establish a manifold deep network containing a nonlinear activation layer, input the heterogeneous feature tensor into the manifold deep network to perform subspace projection calculation, extract the high-dimensional manifold matrix, calculate the adjacency connectivity of the high-dimensional manifold matrix based on the geographic coordinate relationship, and construct the topological network structure based on the adjacency connectivity. S3: Establish a dynamic graph convolutional model, input the topological network structure into the dynamic graph convolutional model to perform node information aggregation calculation, generate time decay weights according to the historical evolution trend of nodes, fuse the time decay weights to correct the forward propagation path, and extract multi-scale fusion vectors. S4: Calculate the distribution gradient value of the multi-scale fusion vector under the preset resource category, superimpose and calculate the cumulative magnitude of the distribution gradient value within the discrete spatial grid cell, generate a two-dimensional density expectation, and calculate the evaluation comprehensive index based on the two-dimensional density expectation.

[0006] As a further aspect of the present invention, the specific steps of S1 are as follows: S11: Acquire the remote sensing image sequence and the ground meteorological parameters, perform spatiotemporal registration and outlier removal on the acquired ground meteorological parameters, resample the remote sensing image sequence and the ground meteorological parameters with inconsistent spatial resolution into the same spatial grid system, and extract effective meteorological feature data. S12: Calculate the difference Euclidean distance between the physical quantities of each pixel value in the remote sensing image sequence and the corresponding nodes in the effective meteorological feature data, evaluate the complementary correlation characteristics between different data sources based on the difference Euclidean distance, and generate a feature weight vector. S13: Based on the feature weight vector, the preprocessed data is concatenated according to the feature dimensions, and tensor product operation is performed in the concatenated high-dimensional space to construct the heterogeneous feature tensor.

[0007] As a further aspect of the present invention, the specific steps of S2 are as follows: S21: Establish a manifold deep network containing nonlinear activation layers, input the heterogeneous feature tensor into the manifold deep network to perform subspace projection calculation, use the encoding structure of the autoencoder to perform layer-by-layer dimensionality reduction mapping on the heterogeneous feature tensor, and extract the high-dimensional manifold matrix; S22: Calculate the geographic spatial distance and feature spatial distance between any two nodes in the high-dimensional manifold matrix based on geographic coordinate relationships, and perform weighted fusion calculation on the calculated geographic spatial distance and feature spatial distance to calculate the adjacency connectivity between adjacent nodes. S23: Based on the adjacency connectivity, set a connectivity filtering threshold, remove redundant connection edges with connectivity lower than the connectivity filtering threshold, extract node connection relationships and combine them into a topology network edge set to construct the topology network structure.

[0008] As a further aspect of the present invention, the specific steps of S3 are as follows: S31: Establish a dynamic graph convolution model, input the topological network structure into the dynamic graph convolution model to perform node information aggregation calculation, aggregate the feature representations of the target node and its first-order and second-order neighbor nodes, and extract the node state evolution vector; S32: Analyze the decay change law of the node state evolution vector within a continuous time step using a recurrent neural network unit, calculate the decay coefficient by combining the time interval between the current moment and the historical moment, and generate the time decay weight. S33: The time decay weights are fused to dynamically weight the inter-layer connection weights within the dynamic graph convolutional model, the forward propagation path is corrected, and pooling fusion is performed on the multi-layer convolutional output features along the corrected forward propagation path to extract the multi-scale fusion vector.

[0009] As a further aspect of the present invention, the specific steps of S4 are as follows: S41: Obtain the prior distribution probability of the preset resource category, calculate the feature projection error of the multi-scale fusion vector in the preset resource category space, and use the activation function to perform classification probability mapping calculation based on the feature projection error to obtain the distribution gradient value. S42: Obtain the boundary coordinates of the discrete spatial grid cell, superimpose and calculate the distribution density change trend of the distribution gradient value in the discrete spatial grid cell, perform numerical integration along the latitude and longitude direction of the discrete spatial grid cell, calculate the cumulative size, and generate the two-dimensional density expectation. S43: Based on the two-dimensional density expectation, establish a multi-index comprehensive evaluation system, use the entropy weight method to calculate the objective evaluation weights of the comprehensive evaluation indicators and perform weighted summation and normalization processing to calculate the comprehensive evaluation index; The prior distribution probability is the basic distribution probability value of each resource category in the target area, obtained based on historical surveying data and land planning information statistics.

[0010] As a further aspect of the present invention, the process of calculating the difference Euclidean distance specifically includes: The pixel brightness features of the remote sensing image sequence and the numerical features of the effective meteorological feature data are obtained. The pixel brightness features and numerical features are then dimensionally aligned to generate a feature alignment matrix. The absolute values ​​of the differences between the elements in each dimension of the feature alignment matrix are squared one by one. The squared results of each dimension are weighted and summed using a preset weight allocation vector to calculate the sum of squared distances. Extract the sum of squared distance values, perform a square root smoothing operation on the sum of squared distance values, and adjust the values ​​in combination with the data normalization coefficient to obtain the difference Euclidean distance; The data normalization coefficient is a scaling product parameter used to uniformly map the value range of different modal data to a preset standard interval.

[0011] As a further aspect of the present invention, the process of extracting the high-dimensional manifold matrix specifically includes: The heterogeneous feature tensor is obtained, and a multilayer perceptron is used to perform nonlinear mapping on the channel dimension of the heterogeneous feature tensor to generate a nonlinear feature matrix. The nonlinear feature matrix is ​​input into the manifold dimensionality reduction module to extract intrinsic geometric structural features, calculate the reconstruction error coefficients of the node neighborhood, and obtain the local linear embedding weights. The nonlinear feature matrix is ​​reconstructed by performing dimensionality reduction based on the local linear embedding weights, and the local neighborhood relationships of the data are mapped globally to preserve locality. The reconstructed features are then projected into a continuous low-dimensional manifold space to extract the high-dimensional manifold matrix.

[0012] As a further aspect of the present invention, the process of generating the time decay weight specifically includes: Obtain the node state evolution vector, analyze the feature evolution rate of the node under different time slices, compare the state differences of different historical stages, and calculate the feature change rate ratio between adjacent time nodes; The ratio of the characteristic rate of change is input into the exponential decay function for nonlinear smoothing transformation. The numerical amplitude of the transformation result is limited by a preset decay control factor to generate a basic decay fluctuation sequence. Each element in the basic decay fluctuation sequence is subjected to time-series low-pass filtering, and the changes and fluctuations between each time slice are smoothed through a sliding window mechanism. Then, a standardized scaling operation is performed to generate the time decay weight.

[0013] As a further aspect of the present invention, the process of generating the two-dimensional density expectation specifically includes: Obtain the distribution gradient values ​​and the boundary coordinate sequence of the discrete spatial grid cells, establish the mapping relationship between the distribution gradient values ​​and the discrete spatial grid cells based on spatial location information, and construct a geospatial coordinate system; Under the mapping and association constraints, numerical integration is performed on the distributed gradient values ​​along the horizontal and vertical dimensions of the discrete spatial grid cells to calculate the concentration of elements within the discrete spatial grid cells and to calculate the cumulative size within each discrete spatial grid cell. The cumulative size of all discrete spatial grid cells is obtained, and the numerical smoothing transition calculation of the grid blank area is performed using a spatial two-dimensional bilinear interpolation processing algorithm to construct the density distribution surface feature and generate the two-dimensional density expectation.

[0014] A land and resources utilization assessment system based on multi-source data fusion, the system being used to implement the aforementioned land and resources utilization assessment method based on multi-source data fusion, the system comprising: The heterogeneous tensor construction module is used to acquire remote sensing image sequences through a spaceborne multispectral camera, link with a micro-monitoring station to collect ground meteorological parameters, calculate the difference Euclidean distance between the remote sensing image sequence and the ground meteorological parameters, and expand the feature dimension based on the difference Euclidean distance to construct a heterogeneous feature tensor. The topology network construction module is used to build a manifold deep network containing a nonlinear activation layer. The heterogeneous feature tensor is input into the manifold deep network to perform subspace projection calculation, extract a high-dimensional manifold matrix, calculate the adjacency connectivity of the high-dimensional manifold matrix based on geographic coordinate relationships, and construct a topology network structure based on the adjacency connectivity. The fusion vector extraction module is used to establish a dynamic graph convolution model, input the topological network structure into the dynamic graph convolution model to perform node information aggregation calculation, generate time decay weights based on the historical evolution trend of nodes, fuse the time decay weights to correct the forward propagation path, and extract multi-scale fusion vectors. The comprehensive index evaluation module is used to calculate the distribution gradient value of the multi-scale fusion vector under the preset resource category, superimpose and calculate the cumulative magnitude of the distribution gradient value within the discrete spatial grid cell, generate a two-dimensional density expectation, and calculate the comprehensive index based on the two-dimensional density expectation.

[0015] Compared with the prior art, the advantages and positive effects of the present invention are as follows: In this invention, a spatiotemporal heterogeneous tensor construction mechanism is introduced to perform manifold spatial mapping on the features of acquired remote sensing images and ground sensing vector features. Input sources of different modal dimensions are projected into a shared high-dimensional manifold cluster to eliminate differences in feature spatial distribution. Cross-modal feature nonlinear interaction constraints are established to guide the adaptive iterative update of the mapping matrix. Dynamic graph convolutional networks are used to reconstruct the edge strength between nodes based on the spatial proximity topology and temporal evolution sequence of land parcels. Aggregation weights are dynamically allocated to perform graph domain fusion calculation on the multi-scale correlation features of land parcels. This breaks the structural representation bias caused by relying on static projection rules and single-scale features. The land parcel utilization index is directly analyzed from the graph convolutional aggregation results to achieve an accurate expression of resource status under complex terrain. Attached Figure Description

[0016] Figure 1 This is the main flowchart of the land and resources utilization assessment method based on multi-source data fusion of the present invention; Figure 2 A flowchart for constructing heterogeneous feature tensors for this invention; Figure 3 A flowchart illustrating the construction of the topology network structure for this invention; Figure 4 This is a flowchart illustrating the process of extracting multi-scale fusion vectors for this invention. Figure 5 This is a flowchart illustrating the calculation and evaluation of the comprehensive index for this invention. Detailed Implementation

[0017] To make the objectives, technical solutions, and advantages of this invention clearer, the software-based technical solution is described in detail below with reference to system architecture diagrams and embodiments. It should be understood that the specific embodiments described herein are only for explaining the technical solutions of this invention and do not constitute a limitation on the scope of protection.

[0018] In the description of this invention, the system architecture relationships or data processing flows indicated by terms such as "layer," "module," "interface," "data flow," "client," and "server" are all defined based on the architecture diagram or flowchart corresponding to the embodiments. This way of describing is only used to clearly illustrate the logical relationships between the elements in the technical solution, and not to limit the physical deployment form. The term "multiple" includes two or more technical units, including but not limited to multiple data nodes, processing threads, service instances, or functional components and other scalable elements. The specific number is determined according to the actual business scenario and needs to be specifically specified.

[0019] Please see Figure 1 and Figure 2 This invention provides a technical solution: a land and resources utilization assessment method based on multi-source data fusion, comprising the following steps: S1: Acquire remote sensing image sequences through a spaceborne multispectral camera, link with a micro-monitoring station to collect ground meteorological parameters, calculate the difference Euclidean distance between the remote sensing image sequence and the ground meteorological parameters, expand the feature dimension based on the difference Euclidean distance, and construct a heterogeneous feature tensor. The specific steps of S1 are as follows: S11: Acquire remote sensing image sequences and ground meteorological parameters, perform spatiotemporal registration and outlier removal on the acquired ground meteorological parameters, resample remote sensing image sequences and ground meteorological parameters with inconsistent spatial resolutions into the same spatial grid system, and extract effective meteorological feature data. S12: Calculate the difference Euclidean distance between the physical quantities of each pixel value in the remote sensing image sequence and the corresponding nodes in the effective meteorological feature data, evaluate the complementary correlation characteristics between different data sources based on the difference Euclidean distance, and generate a feature weight vector. S13: Based on the feature weight vector, the preprocessed data is concatenated according to the feature dimensions, and tensor product operation is performed in the concatenated high-dimensional space to construct a heterogeneous feature tensor; The process of calculating the difference Euclidean distance specifically includes: The pixel brightness features and numerical features of effective meteorological feature data of remote sensing image sequences are obtained. The pixel brightness features and numerical features are dimensionally aligned to generate a feature alignment matrix. The absolute value of the difference between the elements of each dimension in the feature alignment matrix is ​​squared item by item. The squared results of each dimension are weighted and summed using a preset weight allocation vector to calculate the sum of squared distances. Extract the sum of squared distance values, perform a square root smoothing operation on the sum of squared distance values, and adjust the values ​​in combination with the data normalization coefficient to obtain the difference Euclidean distance; Remote sensing image sequences were acquired using a spaceborne multispectral camera, and ground meteorological parameters were collected in conjunction with micro-monitoring stations. A zero-day period was calculated between the remote sensing image sequences and the ground meteorological parameters to directionally acquire remote sensing image sequences covering the target monitoring area (an agricultural ecological experimental zone of 120 square kilometers). The acquired bands included red (center wavelength 0.66 μm), green (center wavelength 0.55 μm), blue (center wavelength 0.47 μm), and near-infrared (center wavelength 0.83 μm) bands. The initial spatial resolution of the images was 8 meters, and the pixel depth was 16 bits. Simultaneously, 50 IoT micro-monitoring stations deployed within the target area were connected to collect ground meteorological parameters hourly via the MQTT protocol, specifically including Celsius temperature, relative humidity percentage, and rainfall in millimeters.

[0020] Spatiotemporal registration was performed on the collected ground meteorological parameters, and the average meteorological monitoring values ​​within one hour before and after the remote sensing satellite's transit time (set to 10:30 AM daily) were extracted as the registration benchmark. In the outlier removal process, the mean and standard deviation over the past 30 days were calculated for each meteorological parameter. If a temperature value uploaded by a monitoring station exceeded the range of the mean plus or minus three times the standard deviation (e.g., if the mean is 25.4°C and the standard deviation is 1.2, a temperature higher than 29.0°C or lower than 21.8°C), it was identified as an outlier and replaced with the median of the same time from three neighboring normal monitoring stations. Remote sensing image sequences with inconsistent spatial resolution (8-meter resolution) and ground meteorological parameters (discrete point data) were resampled into the same spatial grid system. The specific process involves constructing a standard two-dimensional geographic grid with 1500 rows by 1500 columns, recalculating the pixel grid coordinates of the remote sensing images using bilinear interpolation, and smoothly extending the discrete station data to the entire grid surface using inverse distance weighted interpolation for ground meteorological parameters. The distance power exponent is set to 2, the effective influence radius is set to 5000 meters, and the effective meteorological feature data after grid alignment is extracted.

[0021] The pixel brightness features and numerical features of effective meteorological data from remote sensing image sequences are acquired. Dimensional alignment is performed on the pixel brightness features and numerical features. Specifically, maximum-minimum normalization is performed, mapping pixel brightness features (0 to 65535) to the 0-1 interval, and Celsius temperature features are mapped to the 0-1 interval according to historical extreme values ​​(lower limit -20.0, upper limit 45.0), generating a feature alignment matrix containing 7 channel dimensions. The absolute values ​​of the differences between elements in each dimension of the feature alignment matrix are squared item by item. The squared results of each dimension are then weighted and summed using a preset weight allocation vector to calculate the sum of squared distances. The preset weight allocation vector is specifically set as follows: near-infrared band weight 0.35, red band weight 0.15, green band weight 0.10, blue band weight 0.10, Celsius temperature weight 0.15, relative humidity weight 0.10, and rainfall weight 0.05.

[0022] The sum of squared distances is extracted, and a square root smoothing operation is performed on the sum of squared distances. This is then adjusted using a data normalization coefficient to obtain the difference in Euclidean distance. The data normalization coefficient is set to 0.85 to uniformly map the value range of different modal data to a preset standard interval. Experimental verification shows that the stability of feature fusion is highest when the normalization coefficient is in the range of 0.8 to 0.9; in this embodiment, 0.85 is selected. The calculation process is as follows: the obtained feature alignment matrix is... Squared difference terms in each dimension With weight Multiply and sum to get the result. Then, the square root operation is performed on it and multiplied by 0.85 to obtain the final difference Euclidean distance. Based on the feature weight vector, the preprocessed data is concatenated according to the feature dimensions, and tensor product operation is performed in the concatenated high-dimensional space to construct a heterogeneous feature tensor.

[0023] The formula for calculating the difference Euclidean distance is: ; in, Represents the difference in Euclidean distance. The normalization coefficient represents the data, with a value of 0.85. This represents the total number of feature dimensions, with a value of 7. Represents the first element in the preset weight allocation vector. Weight values ​​for each dimension, This represents the normalized brightness value of the corresponding pixel in the remote sensing image sequence. This represents the normalized value of the corresponding node in the valid meteorological feature data.

[0024] The aforementioned spaceborne multispectral camera refers to an imaging device installed on an artificial satellite that can record electromagnetic radiation information of ground objects within multiple specific wavelength ranges.

[0025] The aforementioned inverse distance weighted interpolation method is an algorithm for numerically estimating unknown points based on spatial proximity. The closer the known monitoring point is, the greater its influence on the unknown point.

[0026] The tensor product operation mentioned above refers to a multiplication operation between two vectors or tensors, which is used to linearly expand low-dimensional features to a high-dimensional manifold space to capture deep nonlinear correlations between multi-source data.

[0027] Table 1. Heterogeneous Feature Tensor Mesh Alignment Parameters ; See Table 1, which lists the initial physical quantity ranges of each dimension of data, the baseline denominator used for normalization, and the final weight allocation ratios for the weighted calculation during the construction of the heterogeneous feature tensor. Scaling product parameters.

[0028] Please see Figure 1 and Figure 3 S2: Establish a manifold deep network with nonlinear activation layers, input heterogeneous feature tensors into the manifold deep network to perform subspace projection calculation, extract high-dimensional manifold matrix, calculate the adjacency connectivity of the high-dimensional manifold matrix based on geographic coordinate relationships, and construct the topological network structure based on the adjacency connectivity. The specific steps of S2 are as follows: S21: Establish a manifold deep network containing nonlinear activation layers, input heterogeneous feature tensors into the manifold deep network to perform subspace projection calculations, and use the encoding structure of an autoencoder to perform layer-by-layer dimensionality reduction mapping on the heterogeneous feature tensors to extract high-dimensional manifold matrices. S22: Calculate the geographic spatial distance and feature spatial distance between any two nodes in the high-dimensional manifold matrix based on geographic coordinate relationships, and perform weighted fusion calculation on the calculated geographic spatial distance and feature spatial distance to calculate the adjacency connectivity between adjacent nodes. S23: Set a connectivity filtering threshold based on adjacency connectivity, remove redundant connection edges with connectivity lower than the connectivity filtering threshold, extract node connection relationships and combine them into a set of topology network edges to construct the topology network structure; The process of extracting high-dimensional manifold matrices specifically includes: Obtain heterogeneous feature tensors, and use a multilayer perceptron to perform nonlinear mapping on the channel dimensions of the heterogeneous feature tensors to generate nonlinear feature matrices. The nonlinear feature matrix is ​​input into the manifold dimensionality reduction module to extract intrinsic geometric structural features, calculate the reconstruction error coefficients of the node neighborhood, and obtain the local linear embedding weights. The nonlinear feature matrix is ​​reconstructed by performing dimensionality reduction based on local linear embedding weights, and the local neighborhood relationships of the data are mapped globally to preserve locality. The reconstructed features are then projected into a continuous low-dimensional manifold space to extract the high-dimensional manifold matrix.

[0029] A manifold deep network with nonlinear activation layers was constructed, consisting of one input layer, three hidden layers, and one bottleneck layer. The number of neurons in the hidden layers was set to 1024, 512, and 256, respectively, and the bottleneck layer (i.e., the output layer) was set to 128. Heterogeneous feature tensors were input into the manifold deep network in spatial batches (batch size set to 64) to perform subspace projection computation. The encoding structure of the autoencoder was used to perform layer-by-layer dimensionality reduction mapping of the heterogeneous feature tensors. LeakyReLU was used as a nonlinear activation layer between the hidden layers (the negative half-axis slope parameter was set to 0.01) to extract the high-dimensional manifold matrix.

[0030] Heterogeneous feature tensors are obtained, and a multilayer perceptron (MLP) is used to perform nonlinear mapping on the channel dimensions of the heterogeneous feature tensors to generate a nonlinear feature matrix. This nonlinear feature matrix is ​​then input into a manifold dimensionality reduction module to extract intrinsic geometric structural features, calculate the reconstruction error coefficients of the node neighborhood, and obtain the local linear embedding weights. During the calculation process, for each feature node, its 15 nearest neighbors in the feature space are searched. By solving the constrained optimization problem, each node can be optimally represented by a linear combination of its neighboring nodes, thus obtaining weight coefficients that reflect the local topology.

[0031] Dimensionality reduction and reconstruction of the nonlinear feature matrix are performed based on locally linear embedding weights. Local neighborhood relationships are mapped globally to preserve locality, and the reconstructed features are projected onto a continuous low-dimensional manifold space, extracting a 128-dimensional high-dimensional manifold matrix. The geographic distance and feature space distance between any two nodes in the high-dimensional manifold matrix are calculated based on geographic coordinate relationships. The latitude and longitude coordinates of the nodes are extracted and substituted into the spherical semi-versus equation to calculate the actual physical distance, where the Earth's radius is taken as 6,371,000 meters. The results are then mapped to a normalized interval of 0 to 1 to obtain the geographic distance. Simultaneously, the cosine similarity of the feature vectors is calculated in the 128-dimensional manifold space, and its inverse value is taken as the feature space distance.

[0032] The calculated geospatial distance and feature spatial distance are weighted and fused to calculate the adjacency connectivity between adjacent nodes. The geospatial distance is weighted at 0.65, and the feature spatial distance at 0.35. Experimental comparisons show that when the spatial weight is higher than 0.6, the topological network improves the accuracy of identifying land resource boundaries by 12.5% ​​compared to the equal-weight model. A connectivity filtering threshold is set based on adjacency connectivity, with a baseline threshold of 0.78. This threshold is selected based on the connectivity distribution of 200 known sample points within the target area, effectively filtering out pseudo-connections caused by sensor noise.

[0033] Redundant edges with connectivity below the connectivity threshold of 0.78 are removed. Node connections are extracted and combined into a set of topological network edges, constructing a topological network structure that includes node features and connection weights. This structure can accurately reflect the dual constraints of physical adjacency and feature similarity in land resource utilization.

[0034] The aforementioned autoencoder refers to a type of unsupervised learning neural network that extracts the most essential features from the data by compressing the input data into a hidden representation (encoding) and then attempting to reconstruct the input from the hidden representation (decoding).

[0035] The aforementioned local linear embedding refers to a nonlinear dimensionality reduction algorithm. Its core idea is to assume that the manifold is locally linear, and each point can be linearly reconstructed from its neighboring points, while maintaining the reconstruction weights unchanged after dimensionality reduction.

[0036] The connectivity screening threshold mentioned above refers to the numerical limit used to determine whether there is a substantial connection between two nodes when constructing a graph structure. Only connections with a value greater than this threshold are retained in the topology network.

[0037] Table 2 Comparison of Topology Network Connectivity Screening Experiments ; See Table 2, which records the edge retention, network clustering coefficient, and accuracy comparison of the final identification of land resource utilization boundaries under different connectivity filtering threshold settings. Experiments show that the performance is optimal when the threshold is set to 0.78.

[0038] Please see Figure 1 and Figure 4 S3: Establish a dynamic graph convolution model, input the topology network structure into the dynamic graph convolution model to perform node information aggregation calculation, generate time decay weights based on the historical evolution trend of nodes, fuse the time decay weights to correct the forward propagation path, and extract multi-scale fusion vectors. The specific steps for S3 are as follows: S31: Establish a dynamic graph convolution model, input the topological network structure into the dynamic graph convolution model to perform node information aggregation calculation, aggregate the feature representations of the target node and its first-order and second-order neighbor nodes, and extract the node state evolution vector. S32: Analyze the decay change law of the node state evolution vector within a continuous time step using recurrent neural network units, and calculate the decay coefficient by combining the time interval between the current moment and the historical moment to generate time decay weight. S33: The time decay weights of the fusion model are dynamically weighted to calculate the inter-layer connection weights within the dynamic graph convolution model, correcting the forward propagation path. Pooling and fusion operations are then performed on the multi-layer convolution output features along the corrected forward propagation path to extract the multi-scale fusion vector. The process of generating time decay weights specifically includes: Obtain the node state evolution vector, analyze the feature evolution rate of the node under different time slices, compare the state differences of different historical stages, and calculate the feature change rate ratio between adjacent time nodes. The ratio of the characteristic rate of change is input into the exponential decay function for nonlinear smoothing transformation. The numerical amplitude of the transformation result is limited by a preset decay control factor to generate a basic decay fluctuation sequence. Each element in the basic decay fluctuation sequence is subjected to time-series low-pass filtering, and the changes and fluctuations between time slices are smoothed through a sliding window mechanism. Then, a standardized scaling operation is performed to generate time decay weights.

[0039] A dynamic graph convolutional model is established, consisting of two graph convolutional layers (GCNLayer) and one temporal recurrent unit layer. The first graph convolutional layer has an output dimension of 64, and the second layer has an output dimension of 32. The topological network structure is input into the dynamic graph convolutional model to perform node information aggregation calculations, aggregating the feature representations of the target node and its first-order and second-order neighbor nodes. During the calculation process, adjacency matrix normalization is used to ensure the stability of numerical propagation, and node state evolution vectors are extracted.

[0040] A recurrent neural network (LSTM) with 32 hidden nodes was used to analyze the decay pattern of node state evolution vectors over a continuous time step (set as the past 12 monitoring periods). The node state evolution vectors were obtained, and the feature evolution rate of nodes at different time slices was analyzed. State differences across different historical stages were compared, and the ratio of feature change rates between adjacent time nodes was calculated. For example, if the feature value at the current time is 0.85 and the previous time was 0.80, then the ratio of change rates is 0.0625.

[0041] The ratio of characteristic rates of change is input into an exponential decay function for nonlinear smoothing transformation. A preset decay control factor is used to limit the numerical amplitude of the transformation result, generating a basic decaying fluctuation sequence. The preset decay control factor is set to 0.65. Each element in the basic decaying fluctuation sequence undergoes time-series low-pass filtering, and a sliding window mechanism (window width set to 3 time periods) is used to smooth the fluctuations between time slices. The smoothed result is then standardized and multiplied by a base gain coefficient of 1.15 to generate time decay weights.

[0042] The time decay weights are dynamically weighted to calculate the inter-layer connection weights within the dynamic graph convolutional model, thus correcting the forward propagation path. During the computation of each graph convolutional layer, the original adjacency matrix weights are multiplied by the corresponding time decay weights. Average pooling is then performed on the multi-layer convolutional output features along the corrected forward propagation path to extract a multi-scale fusion vector. This vector has a fixed dimension of 32 and integrates composite features of spatial topology and temporal dynamic evolution.

[0043] The aforementioned dynamic graph convolution model refers to a deep learning model that can perform convolution operations on graph structure data and dynamically update graph connection weights or node features over time.

[0044] The aforementioned time decay weight refers to a numerical coefficient used to characterize the degree of influence of historical information on the current state. The further back in time or the more drastic the historical changes, the lower the contribution weight of the feature to the current assessment.

[0045] The aforementioned forward propagation path refers to the logic and weight combination that allows information to travel from the input layer through the hidden layer to the output layer in a neural network.

[0046] Table 3. Parameter Configuration Table for Dynamic Graph Convolution Model ; See Table 3, which details the technical basis for setting and selecting key hyperparameters in the implementation of the dynamic graph convolution model, ensuring the model's accurate modeling of the dynamic evolution of land resources.

[0047] Please see Figure 1 and Figure 5 S4: Calculate the distribution gradient values ​​of the multi-scale fusion vector under the preset resource category, superimpose and measure the cumulative magnitude of the distribution gradient values ​​within the discrete spatial grid cell, generate the two-dimensional density expectation, and calculate the comprehensive evaluation index based on the two-dimensional density expectation. The specific steps of S4 are as follows: S41: Obtain the prior distribution probability of the preset resource category, calculate the feature projection error of the multi-scale fusion vector in the preset resource category space, combine the feature projection error with the activation function to perform classification probability mapping calculation, and obtain the distribution gradient value. S42: Obtain the boundary coordinates of the discrete spatial grid cell, overlay and calculate the distribution density change trend of the distribution gradient value in the discrete spatial grid cell, perform numerical integration along the latitude and longitude direction of the discrete spatial grid cell, calculate the cumulative magnitude, and generate the two-dimensional density expectation. S43: Based on the two-dimensional density expectation, a multi-index comprehensive evaluation system is established. The entropy weight method is used to calculate the objective evaluation weights of the comprehensive evaluation indicators and perform weighted summation and normalization to calculate the comprehensive evaluation index. The prior distribution probability is the basic distribution probability value of each resource category in the target area, obtained based on historical surveying data and land planning information statistics. The process of generating the two-dimensional density expectation specifically includes: Obtain the numerical distribution gradient and the boundary coordinate sequence of discrete spatial grid cells, establish the mapping relationship between the numerical distribution gradient and discrete spatial grid cells based on spatial location information, and construct a geospatial coordinate system; Under the mapping and association constraints, numerical integration is performed on the distributed gradient values ​​along the horizontal and vertical dimensions of the discrete spatial grid cells to statistically analyze the concentration of elements within the discrete spatial grid cells and calculate the cumulative size within each discrete spatial grid cell. The cumulative size of all discrete spatial grid cells is obtained, and the numerical smoothing transition calculation of the grid blank area is performed using a spatial two-dimensional bilinear interpolation processing algorithm to construct the density distribution surface feature and generate the two-dimensional density expectation.

[0048] Obtain the prior distribution probability of the preset resource categories. Based on the historical mapping data of the target area over the past 5 years, the basic distribution probability of each resource category is statistically obtained: agricultural planting land is 0.42, ecological forest land is 0.28, industrial construction land is 0.18, and natural water bodies are 0.12. Calculate the feature projection error of the multi-scale fusion vector in the preset resource category space.

[0049] By combining feature projection error and using activation functions, classification probability mapping is calculated to obtain distribution gradient values. Specifically, the Softmax function is used to map the 32-dimensional multi-scale fusion vector into a 4-dimensional probability distribution vector, where each dimension represents the probability that the grid belongs to the corresponding resource category. This probability is then differiated from the prior distribution probability to extract the distribution gradient values ​​reflecting the degree of resource utilization anomalies. The boundary coordinates of discrete spatial grid cells are obtained, and a mapping relationship is established between the distribution gradient values ​​and discrete spatial grid cells, constructing a 1500x1500 geospatial coordinate system.

[0050] Under the constraints of mapping association, numerical integration of the distribution gradient values ​​is performed along the horizontal and vertical dimensions of the discrete spatial grid cells. Simpson's integral method is used to calculate the concentration of features within each grid cell, and the cumulative magnitude within the grid is calculated. The cumulative magnitudes of all discrete spatial grid cells are obtained, and a two-dimensional bilinear interpolation algorithm is used to perform numerical smoothing calculations on the blank areas of the grid, constructing a continuous density distribution surface feature and generating a two-dimensional density expectation.

[0051] A multi-indicator comprehensive evaluation system is established based on the two-dimensional density expectation. The objective evaluation weights of the comprehensive evaluation indicators are calculated using the entropy weight method. First, the three indicators—land utilization rate, resource coverage density, and development intensity—are standardized. The entropy value and difference coefficient of each indicator are calculated, and the final weight allocation is determined: land utilization rate 0.38, resource coverage density 0.32, and development intensity 0.30. The two-dimensional density expectation values ​​are then substituted into the evaluation system for weighted summation and normalization to calculate the comprehensive evaluation index, with a value range of 0 to 100.

[0052] Experimental data show that when using this method to assess land resource utilization, the accuracy rate for identifying newly added construction land reaches 94.2%, which is 16.8% higher than the traditional method based on a single remote sensing image. Furthermore, the coefficient of variation of the comprehensive assessment index is reduced by 11.5%, demonstrating the stability and accuracy of the assessment results.

[0053] The aforementioned distribution gradient value refers to the rate of change between the actual monitored resource distribution probability and the expected planning probability, used to characterize the changes in the direction and intensity of land resource utilization.

[0054] The aforementioned two-dimensional density expectation refers to the expected intensity of resource distribution within a geospatial grid, taking into account both the characteristic distribution probability and the spatial integral accumulation. It reflects the degree of resource intensification per unit area.

[0055] The aforementioned entropy weighting method refers to an objective weighting method that determines the weight of an indicator based on the amount of information provided by the observed values ​​of each indicator. The greater the amount of information (i.e., the smaller the entropy), the higher the weight.

[0056] The results show that by fusing multi-source data and using dynamic graph convolution modeling, it is possible to accurately quantify the temporal and spatial utilization characteristics of land resources, providing scientific data support for land resource protection and rational planning.

[0057] A land and resources utilization assessment system based on multi-source data fusion, the system being used to execute the aforementioned land and resources utilization assessment method based on multi-source data fusion, the system comprising: The heterogeneous tensor construction module is used to acquire remote sensing image sequences through a spaceborne multispectral camera, link with a micro-monitoring station to collect ground meteorological parameters, calculate the difference Euclidean distance between the remote sensing image sequence and the ground meteorological parameters, and expand the feature dimension based on the difference Euclidean distance to construct a heterogeneous feature tensor. The topology network construction module is used to build a manifold deep network containing nonlinear activation layers. It inputs heterogeneous feature tensors into the manifold deep network to perform subspace projection calculations, extracts high-dimensional manifold matrices, calculates the adjacency connectivity of the high-dimensional manifold matrix based on geographic coordinate relationships, and constructs the topology network structure based on the adjacency connectivity. The fusion vector extraction module is used to build a dynamic graph convolution model. The topological network structure is input into the dynamic graph convolution model to perform node information aggregation calculation. Time decay weights are generated based on the historical evolution trend of nodes. The time decay weights are fused to correct the forward propagation path and extract multi-scale fusion vectors. The comprehensive index evaluation module is used to calculate the distribution gradient values ​​of the multi-scale fusion vector under the preset resource category, superimpose and calculate the cumulative magnitude of the distribution gradient values ​​within the discrete spatial grid cell, generate the two-dimensional density expectation, and calculate the comprehensive index based on the two-dimensional density expectation.

[0058] The above embodiments illustrate preferred embodiments of the present invention. Any equivalent adjustments to the technical solution based on software engineering methods are within the scope of protection, including but not limited to: implementing algorithm logic using different programming languages, refactoring functional modules into services, adjusting data interaction protocols, and optimizing resource scheduling strategies. Any implementation scheme derived from reasonable modifications to the data processing flow, service call chain, or system architecture layer without departing from the core technology of the present invention should be considered within the protection scope defined by the technical solution of the present invention.

Claims

1. A land and resources utilization assessment method based on multi-source data fusion, characterized in that, Includes the following steps: S1: Acquire remote sensing image sequences through a spaceborne multispectral camera, collect ground meteorological parameters in conjunction with a micro-monitoring station, calculate the difference Euclidean distance between the remote sensing image sequences and the ground meteorological parameters, expand the feature dimension based on the difference Euclidean distance, and construct a heterogeneous feature tensor. S2: Establish a manifold deep network containing a nonlinear activation layer, input the heterogeneous feature tensor into the manifold deep network to perform subspace projection calculation, extract the high-dimensional manifold matrix, calculate the adjacency connectivity of the high-dimensional manifold matrix based on the geographic coordinate relationship, and construct the topological network structure based on the adjacency connectivity. S3: Establish a dynamic graph convolutional model, input the topological network structure into the dynamic graph convolutional model to perform node information aggregation calculation, generate time decay weights according to the historical evolution trend of nodes, fuse the time decay weights to correct the forward propagation path, and extract multi-scale fusion vectors. S4: Calculate the distribution gradient value of the multi-scale fusion vector under the preset resource category, superimpose and calculate the cumulative magnitude of the distribution gradient value within the discrete spatial grid cell, generate a two-dimensional density expectation, and calculate the evaluation comprehensive index based on the two-dimensional density expectation.

2. The land and resources utilization assessment method based on multi-source data fusion according to claim 1, characterized in that, The specific steps of S1 are as follows: S11: Acquire the remote sensing image sequence and the ground meteorological parameters, perform spatiotemporal registration and outlier removal on the acquired ground meteorological parameters, resample the remote sensing image sequence and the ground meteorological parameters with inconsistent spatial resolution into the same spatial grid system, and extract effective meteorological feature data. S12: Calculate the difference Euclidean distance between the physical quantities of each pixel value in the remote sensing image sequence and the corresponding nodes in the effective meteorological feature data, evaluate the complementary correlation characteristics between different data sources based on the difference Euclidean distance, and generate a feature weight vector. S13: Based on the feature weight vector, the preprocessed data is concatenated according to the feature dimensions, and tensor product operation is performed in the concatenated high-dimensional space to construct the heterogeneous feature tensor.

3. The land and resources utilization assessment method based on multi-source data fusion according to claim 1, characterized in that, The specific steps of S2 are as follows: S21: Establish a manifold deep network containing nonlinear activation layers, input the heterogeneous feature tensor into the manifold deep network to perform subspace projection calculation, use the encoding structure of the autoencoder to perform layer-by-layer dimensionality reduction mapping on the heterogeneous feature tensor, and extract the high-dimensional manifold matrix; S22: Calculate the geographic spatial distance and feature spatial distance between any two nodes in the high-dimensional manifold matrix based on geographic coordinate relationships, and perform weighted fusion calculation on the calculated geographic spatial distance and feature spatial distance to calculate the adjacency connectivity between adjacent nodes. S23: Based on the adjacency connectivity, set a connectivity filtering threshold, remove redundant connection edges with connectivity lower than the connectivity filtering threshold, extract node connection relationships and combine them into a topology network edge set to construct the topology network structure.

4. The land and resources utilization assessment method based on multi-source data fusion according to claim 1, characterized in that, The specific steps of S3 are as follows: S31: Establish a dynamic graph convolution model, input the topological network structure into the dynamic graph convolution model to perform node information aggregation calculation, aggregate the feature representations of the target node and its first-order and second-order neighbor nodes, and extract the node state evolution vector; S32: Analyze the decay change law of the node state evolution vector within a continuous time step using a recurrent neural network unit, calculate the decay coefficient by combining the time interval between the current moment and the historical moment, and generate the time decay weight. S33: The time decay weights are fused to dynamically weight the inter-layer connection weights within the dynamic graph convolutional model, the forward propagation path is corrected, and pooling fusion is performed on the multi-layer convolutional output features along the corrected forward propagation path to extract the multi-scale fusion vector.

5. The land and resources utilization assessment method based on multi-source data fusion according to claim 1, characterized in that, The specific steps of S4 are as follows: S41: Obtain the prior distribution probability of the preset resource category, calculate the feature projection error of the multi-scale fusion vector in the preset resource category space, and use the activation function to perform classification probability mapping calculation based on the feature projection error to obtain the distribution gradient value. S42: Obtain the boundary coordinates of the discrete spatial grid cell, superimpose and calculate the distribution density change trend of the distribution gradient value in the discrete spatial grid cell, perform numerical integration along the latitude and longitude direction of the discrete spatial grid cell, calculate the cumulative size, and generate the two-dimensional density expectation. S43: Based on the two-dimensional density expectation, establish a multi-index comprehensive evaluation system, use the entropy weight method to calculate the objective evaluation weights of the comprehensive evaluation indicators and perform weighted summation and normalization processing to calculate the comprehensive evaluation index; The prior distribution probability is the basic distribution probability value of each resource category in the target area, obtained based on historical surveying data and land planning information statistics.

6. The land and resources utilization assessment method based on multi-source data fusion according to claim 2, characterized in that, The process of calculating the difference Euclidean distance specifically includes: The pixel brightness features of the remote sensing image sequence and the numerical features of the effective meteorological feature data are obtained. The pixel brightness features and numerical features are then dimensionally aligned to generate a feature alignment matrix. The absolute values ​​of the differences between the elements in each dimension of the feature alignment matrix are squared one by one. The squared results of each dimension are weighted and summed using a preset weight allocation vector to calculate the sum of squared distances. Extract the sum of squared distance values, perform a square root smoothing operation on the sum of squared distance values, and adjust the values ​​in combination with the data normalization coefficient to obtain the difference Euclidean distance; The data normalization coefficient is a scaling product parameter used to uniformly map the value range of different modal data to a preset standard interval.

7. The land and resources utilization assessment method based on multi-source data fusion according to claim 3, characterized in that, The process of extracting the high-dimensional manifold matrix specifically includes: The heterogeneous feature tensor is obtained, and a multilayer perceptron is used to perform nonlinear mapping on the channel dimension of the heterogeneous feature tensor to generate a nonlinear feature matrix. The nonlinear feature matrix is ​​input into the manifold dimensionality reduction module to extract intrinsic geometric structural features, calculate the reconstruction error coefficients of the node neighborhood, and obtain the local linear embedding weights. The nonlinear feature matrix is ​​reconstructed by performing dimensionality reduction based on the local linear embedding weights, and the local neighborhood relationships of the data are mapped globally to preserve locality. The reconstructed features are then projected into a continuous low-dimensional manifold space to extract the high-dimensional manifold matrix.

8. The land and resources utilization assessment method based on multi-source data fusion according to claim 4, characterized in that, The process of generating the time decay weight specifically includes: Obtain the node state evolution vector, analyze the feature evolution rate of the node under different time slices, compare the state differences of different historical stages, and calculate the feature change rate ratio between adjacent time nodes; The ratio of the characteristic rate of change is input into the exponential decay function for nonlinear smoothing transformation. The numerical amplitude of the transformation result is limited by a preset decay control factor to generate a basic decay fluctuation sequence. Each element in the basic decay fluctuation sequence is subjected to time-series low-pass filtering, and the changes and fluctuations between each time slice are smoothed through a sliding window mechanism. Then, a standardized scaling operation is performed to generate the time decay weight.

9. The land and resources utilization assessment method based on multi-source data fusion according to claim 5, characterized in that, The process of generating the two-dimensional density expectation specifically includes: Obtain the distribution gradient values ​​and the boundary coordinate sequence of the discrete spatial grid cells, establish the mapping relationship between the distribution gradient values ​​and the discrete spatial grid cells based on spatial location information, and construct a geospatial coordinate system; Under the mapping and association constraints, numerical integration is performed on the distributed gradient values ​​along the horizontal and vertical dimensions of the discrete spatial grid cells to calculate the concentration of elements within the discrete spatial grid cells and to calculate the cumulative size within each discrete spatial grid cell. The cumulative size of all discrete spatial grid cells is obtained, and the numerical smoothing transition calculation of the grid blank area is performed using a spatial two-dimensional bilinear interpolation processing algorithm to construct the density distribution surface feature and generate the two-dimensional density expectation.

10. A land and resources utilization assessment system based on multi-source data fusion, characterized in that, The system is used to implement the land and resources utilization assessment method based on multi-source data fusion as described in any one of claims 1-9, and the system comprises: The heterogeneous tensor construction module is used to acquire remote sensing image sequences through a spaceborne multispectral camera, link with a micro-monitoring station to collect ground meteorological parameters, calculate the difference Euclidean distance between the remote sensing image sequence and the ground meteorological parameters, and expand the feature dimension based on the difference Euclidean distance to construct a heterogeneous feature tensor. The topology network construction module is used to build a manifold deep network containing a nonlinear activation layer. The heterogeneous feature tensor is input into the manifold deep network to perform subspace projection calculation, extract a high-dimensional manifold matrix, calculate the adjacency connectivity of the high-dimensional manifold matrix based on geographic coordinate relationships, and construct a topology network structure based on the adjacency connectivity. The fusion vector extraction module is used to establish a dynamic graph convolution model, input the topological network structure into the dynamic graph convolution model to perform node information aggregation calculation, generate time decay weights based on the historical evolution trend of nodes, fuse the time decay weights to correct the forward propagation path, and extract multi-scale fusion vectors. The comprehensive index evaluation module is used to calculate the distribution gradient value of the multi-scale fusion vector under the preset resource category, superimpose and calculate the cumulative magnitude of the distribution gradient value within the discrete spatial grid cell, generate a two-dimensional density expectation, and calculate the comprehensive index based on the two-dimensional density expectation.

Images