Automatic metadata identification and extraction method for remote sensing image and spatial data

Abstract

The invention provides an automatic metadata recognition and extraction method for remote sensing images and spatial data, and belongs to the technical field of remote sensing measurement. Multi-source remote sensing image data are collected and preprocessed to establish a multi-modal image set under a unified geographic coordinate system, and an image pyramid is constructed by using a GPU accelerated block parallel processing architecture; a topographic feature extraction algorithm based on topology durability is used to identify stable topographic markers, a geometric registration feature extraction and matching network is constructed, a structured prediction framework based on a probability graph model is used to carry out accurate registration, and adaptive filtering and multi-modal data fusion are carried out on a registered image set. A semi-supervised classification algorithm based on random walk is used for carrying out ground feature classification, metadata elements are automatically identified according to a ground feature classification result and stable topographic markers, and the technical problem that the registration precision is insufficient due to the fact that topographic feature extraction is prone to noise interference in multi-mode remote sensing image registration is solved.

Description

TECHNICAL FIELD

[0001] The present application belongs to the technical field of remote sensing measurement, and in particular, relates to an automatic metadata recognition and extraction method for remote sensing images and spatial data. BACKGROUND

[0002] Remote sensing image metadata recognition and extraction is a basic work of geographic spatial data management. The traditional method adopts a combination of manual annotation and automatic recognition, and obtains metadata information such as spatial resolution, geographic range, coordinate system through image registration, feature extraction and classification processing. In the current multi-modal remote sensing data processing application, optical images, synthetic aperture radar images and laser radar data need to be accurately registered to realize effective fusion. The registration process relies on stable terrain features as control points. The traditional method uses edge detection, gradient calculation or corner extraction algorithm to identify terrain features. However, the traditional terrain feature extraction method is based on local gray level change calculation, which is highly sensitive to image noise and imaging condition change. In high-resolution remote sensing images, noise fluctuations are easily misidentified as terrain structures, resulting in poor stability of the extracted feature points. At the same time, the traditional registration method independently matches the feature points, ignoring the spatial geometric consistency constraints that should be met between adjacent feature points. In the repeated texture area and seasonal change scene, it is easy to produce false matching, and the registration error is accumulated and propagated to the subsequent data fusion and classification link, affecting the accuracy of metadata recognition. That is, there is a technical problem of insufficient registration accuracy in multi-modal remote sensing image registration due to the noise interference of terrain feature extraction in the prior art. SUMMARY

[0003] Therefore, the present application provides an automatic metadata recognition and extraction method for remote sensing images and spatial data, which can solve the technical problem of insufficient registration accuracy in multi-modal remote sensing image registration due to the noise interference of terrain feature extraction in the prior art.

[0004] The application is implemented in the following manner: the application provides an automatic metadata recognition and extraction method for remote sensing images and spatial data, collects multi-source remote sensing image data and performs preprocessing to establish a multi-modal image set under a unified geographic coordinate system, uses a GPU-accelerated block parallel processing architecture to perform hierarchical processing on the multi-modal image set to construct an image pyramid, applies a topological persistence-based terrain feature extraction algorithm to multi-scale filtering of elevation data in the multi-modal image set and calculation of the birth time and death time of connected components under different gray threshold values to construct a persistence barcode graph to recognize stable terrain markers, constructs a geometric registration feature extraction and matching network, uses a structured prediction framework based on a probabilistic graph model to model the spatial relationship of matching point pairs as a Markov random field, and realizes globally optimal matching by a belief propagation algorithm to accurately register the multi-modal image set, applies a non-local mean and edge-aware adaptive filtering algorithm to the registered image set to suppress speckle noise, uses an iterative registration optimization algorithm based on mutual information maximization to realize multi-modal data fusion in combination with an adaptive weighted network, applies a semi-supervised classification algorithm based on random walk to perform feature classification on the fused image, and forms a metadata information set according to the feature classification result and the automatically recognized metadata elements of the stable terrain markers.

[0005] In the preprocessing step of the multi-source remote sensing image data, the optical image is subjected to radiation calibration and atmospheric correction, the synthetic aperture radar image is subjected to multi-view processing and geographic coding, and the laser radar data is subjected to point cloud filtering and elevation normalization.

[0006] In the GPU-accelerated block parallel processing architecture, a quadtree spatial index is used to construct the image pyramid, an adaptive sampling strategy is applied at each image pyramid level, and the CUDA parallel computing framework is used to distribute computing tasks among multiple GPUs.

[0007] In the topological persistence-based terrain feature extraction algorithm, the persistence value is defined as the difference between the death time and the birth time, and topological features with a persistence value greater than twice the average persistence value are recognized as stable terrain markers.

[0008] In the topological persistence-based terrain feature extraction algorithm, a persistence barcode graph is constructed by calculating the Betti number change under different thresholds, and the connectivity features of the topological space are described by the Betti number sequence.

[0009] In the geometric registration feature extraction and matching network, a double-branch Siamese architecture is used, the two branches share weight parameters, and the similarity between feature vectors is calculated by a cosine similarity measurement layer after feature extraction.

[0010] In the structured prediction framework based on the probabilistic graph model, the node potential function is defined as the feature similarity, the edge potential function is defined as the change amount of the distance and angle between the matching point pairs, and the global optimal matching is realized by maximizing the joint probability distribution.

[0011] In the confidence propagation algorithm, the probability information is propagated between the graph nodes, each node sends a message to the adjacent node, the message content is the confidence of the node under different matching hypotheses, and after multiple iterations, the edge probability of the node converges.

[0012] In the geometric registration feature extraction and matching network, the global elevation model is combined to constrain the matching search space, the bidirectional verification is carried out through the cyclic consistency test, and the RANSAC robust estimation is used to remove the mis-matching points.

[0013] In the adaptive filtering algorithm based on non-local mean and edge perception, the local directionality is calculated through the structure tensor, a large window smoothing is used in the uniform area, and directional filtering is used in the edge area.

[0014] The present application uses the topological persistence based terrain feature extraction algorithm, quantifies the saliency of the terrain feature by using the persistence homology theory in algebraic topology, represents the stability of the terrain structure as the persistence value, automatically filters the short life cycle topological features caused by noise through the persistence threshold, extracts the stable terrain markers with the persistence value greater than twice the average persistence value, and solves the technical defects that the traditional edge detection method is easily disturbed by noise. The present application constructs a geometric registration feature extraction and matching network, uses a structured prediction framework based on a probabilistic graph model, models the spatial relationship of the matching point pairs as a Markov random field, encodes the local feature similarity through the node potential function, encodes the geometric consistency constraint of the adjacent matching point pairs through the edge potential function, realizes the global optimal matching by propagating the probability information between the graph nodes through the confidence propagation algorithm, significantly reduces the mis-matching rate of the repeated texture area, and ensures the registration accuracy. BRIEF DESCRIPTION OF DRAWINGS

[0015] Figure 1 The flowchart of the method of the present application.

[0016] Figure 2 The structure diagram of the geometric registration feature extraction and matching network.

[0017] Figure 3 The energy distribution diagram of multi-scale wavelet decomposition.

[0018] Figure 4 Figure 6 is a curve chart of iterative convergence of random walk arrival probability.

[0019] Figure 5 Figure 7 is a comparison chart of metadata element extraction accuracy. DETAILED DESCRIPTION

[0020] In order to make the purposes, technical solutions and advantages of the embodiments of the present application clearer, the technical solutions in the embodiments of the present application will be clearly and completely described below.

[0021] As Figure 1 shown in Figure 1, which is a flowchart of an automatic metadata recognition and extraction method for remote sensing images and spatial data provided by the present application, the method comprises the following steps:

[0022] S1, collect multi-source remote sensing image data and perform preprocessing, perform radiation calibration and atmospheric correction on optical images, perform multi-view processing and geographic coding on SAR images, perform point cloud filtering and elevation normalization on laser radar data, and establish a multi-modal image set under a unified geographic coordinate system;

[0023] S2, use a GPU-accelerated block parallel processing architecture to perform layered processing on high-resolution images in the multi-modal image set, use a quadtree spatial index to construct an image pyramid, apply an adaptive sampling strategy at each image pyramid level, and distribute computing tasks among multiple GPUs through a CUDA parallel computing framework;

[0024] S3, apply a topological persistence-based terrain feature extraction algorithm to perform multi-scale filtering on elevation data in the multi-modal image set, calculate the birth time and death time of connected components under different gray threshold values, construct a persistence barcode graph and calculate a Betti number sequence, identify topological features with a persistence value greater than twice the average persistence value as stable topographic markers, and extract ridge and valley linear landform structures;

[0025] S4, construct a geometric registration feature extraction and matching network to accurately register the multi-modal image set, the geometric registration feature extraction and matching network uses a Siamese architecture to extract feature descriptors invariant to geometric transformation and radiation variation, combines global elevation model constraints to search the matching space, performs bidirectional verification through loop consistency checking, uses RANSAC robust estimation to remove mismatched points, and outputs the registered image set;

[0026] S5, applying a non-local mean and edge-aware based adaptive filtering algorithm to the registered SAR image set for speckle noise suppression, calculating local directionality through structure tensor, using large window smoothing in homogeneous areas and directional filtering in edge areas, introducing a multi-scale decomposition framework to separate signal components and noise components of different scales, and outputting a filtered image set;

[0027] S6, implementing multi-modal data fusion on the filtered image set using an iterative registration optimization algorithm based on mutual information maximization, improving convergence speed using Powell direction acceleration method, learning the contribution of different modalities in different scenarios through adaptive weighted network, combining a hybrid strategy of feature-level fusion and decision-level fusion, calculating uncertainty quantization indicators of the fused image, and outputting the fused image;

[0028] S7, applying a semi-supervised classification algorithm based on random walk to the fused image for ground feature classification, using a small number of labeled samples as absorbing states of random walk, constructing a sample similarity matrix as transition probability, calculating the arrival probability of unlabeled samples to each absorbing state by solving a linear equation system, introducing a restart random walk mechanism to prevent walk deviation, and outputting ground feature classification results and classification confidence;

[0029] S8, automatically identifying metadata elements according to the ground feature classification results and the stable terrain markers, starting an active learning strategy to select the most representative samples for manual annotation when the classification confidence is lower than the threshold value 0.75, increasing the multi-scale filtering processing level when the persistence value is lower than 60% of the standard persistence value, and extracting the spatial resolution, geographic range, coordinate system, acquisition time and sensor parameters of the fused image to form a metadata information set.

[0030] Wherein, radiometric calibration refers to the process of converting the gray value of remote sensing image to radiance value or apparent reflectance, which eliminates the instrument response difference by applying sensor calibration parameters. Atmospheric correction refers to the process of eliminating the influence of atmospheric scattering and absorption on remote sensing image radiation value, which calculates and corrects atmospheric transmittance using the 6S radiation transfer model. Multi-view processing refers to the average processing of multiple observation data of SAR image, which reduces the coherent speckle noise through non-coherent superposition in spatial domain or frequency domain. Geocoding refers to the process of converting SAR image from slant range projection to geographic coordinate system, which realizes geometric correction through range-doppler positioning model and elevation data. Point cloud filtering refers to the process of separating ground points and non-ground points from raw laser radar point cloud, which gradually identifies ground points through different window sizes of opening operation using progressive morphological filtering algorithm. Elevation normalization refers to the process of converting laser radar point cloud elevation value to height value relative to the ground, which is realized by subtracting the ground point elevation. Multi-modal image set refers to a set of image data obtained by multiple sensors including optical image, SAR image and laser radar elevation data.

[0031] GPU-accelerated block-parallel processing architecture refers to a computing architecture that utilizes the massive parallel computing power of a graphics processor to process images in blocks. The image is divided into multiple non-overlapping block regions, each of which is assigned to a GPU thread group for independent processing. The quadtree spatial index is a hierarchical spatial data organization structure that recursively divides a two-dimensional space into four quadrants. Each node represents a spatial region and stores the image data summary information of the spatial region. The tree structure is used to quickly locate the target region. The image pyramid is a hierarchical data structure composed of the same image at different resolutions. Each layer of the image has a resolution that is half of the previous layer. The bottom layer is the original high-resolution image, and the top layer is the low-resolution overview image. The adaptive sampling strategy is a method that dynamically adjusts the sampling density according to the spatial complexity of the image content. In areas with rapid texture changes, the sampling point density is increased, and in flat areas with uniform gray levels, the sampling point density is reduced to reduce the computational load. The CUDA parallel computing framework is a parallel computing platform and programming model developed by NVIDIA. It allows developers to utilize the massive parallel processing power of GPUs to accelerate compute-intensive tasks by breaking down the computing task into thousands of threads that execute simultaneously on GPU cores.

[0032] The topological persistence-based terrain feature extraction algorithm uses the persistence homology theory in algebraic topology to quantify the significance of terrain features. The topological persistence-based terrain feature extraction algorithm first performs multi-scale filtering on the elevation data using Gaussian filtering, bilateral filtering, and median filtering to generate elevation fields with different smoothing levels. Then, starting from the lowest gray value, the threshold is gradually increased, and the connected components of the super-level set are calculated at each threshold level. The threshold at which each connected component first appears is recorded as the birth time, and the threshold at which the connected component disappears or merges with other connected components is recorded as the death time. The persistence value is defined as the difference between the death time and the birth time, reflecting the significance of the terrain feature. The Betti number sequence is used to describe the connectivity characteristics of the topological space, where the zeroth Betti number represents the number of connected components, and the first Betti number represents the number of independent loops. By calculating the Betti number changes at different thresholds, a persistence barcode graph is constructed, with each bar corresponding to the persistence value of a topological feature. Cluster analysis is performed on the persistence values, and topological features with persistence values greater than twice the average persistence value are identified as stable terrain markers. Ridge extraction is achieved by identifying the persistence maximum value connection, and valley extraction is achieved by identifying the persistence minimum value connection. The reason why the topological persistence-based terrain feature extraction algorithm is naturally robust to noise is that short-lived topological features are usually caused by noise, and the noise effect is automatically eliminated by the persistence threshold filtering.

[0033] The terrain feature extraction algorithm based on topological persistence provides high-precision terrain structure information for the entire scheme, significantly improves the extraction accuracy of key terrain elements in metadata information, and still can stably identify key landform structures under complex terrain conditions, providing reliable spatial constraints for subsequent geometric registration and feature classification, reducing the geometric distortion correction error caused by terrain undulation, and ensuring the consistency of the extraction results between different resolution images through topological invariance. The terrain feature extraction algorithm based on topological persistence quantifies the significance of terrain features into a persistence value through the persistence homology theory, avoids the problem of relying on manual threshold setting in traditional methods, and realizes the automatic extraction of terrain features. The persistence barcode graph intuitively displays the life cycle of terrain features of different scales, and the long bar corresponds to the stable main terrain structure, and the short bar corresponds to noise or secondary undulation. The evolution process of the topological structure is tracked through the change of the Betti number sequence, and the peak value of the zero-order Betti number corresponds to the number of independent mountain peaks, and the peak value of the first-order Betti number corresponds to the number of closed valleys or basins. The terrain feature extraction algorithm based on topological persistence is particularly important in high-resolution remote sensing image processing, because high-resolution images contain rich terrain details, and traditional edge detection or gradient calculation methods are easily disturbed by noise, while persistence analysis filters out truly stable terrain features through multi-scale comprehensive evaluation.

[0034] The structure of the geometric registration feature extraction and matching network is a double-branch Siamese architecture, each branch contains five layers of convolutional layers for feature extraction, the convolution kernel size is 7x7, 5x5, 3x3, 3x3, 1x1 in turn, and the channel number is 64, 128, 256, 256, 512 in turn, and a batch normalization layer and a ReLU activation function are connected after each convolution. The two branches share the weight parameters to ensure the consistency of feature extraction. After feature extraction, the similarity between feature vectors is calculated by a cosine similarity measurement layer, and the feature point pairs with a similarity greater than 0.85 are considered to be matching points. The geometric registration feature extraction and matching network introduces a structured prediction framework based on a probabilistic graph model, modeling the spatial relationship of matching point pairs as a Markov random field. Each node in the graph represents a candidate matching point pair, and the edge connection represents the geometric consistency constraint between adjacent matching point pairs. The node potential function is defined as the feature similarity, and the edge potential function is defined as the change amount of the distance and angle between matching point pairs. The global optimal matching is achieved by maximizing the joint probability distribution. The belief propagation algorithm is used to propagate probability information between graph nodes. Each node sends a message to adjacent nodes, and the message content is the confidence of the node under different matching hypotheses. After multiple iterations, the edge probability of the node converges. The variational inference technique is used to approximate the posterior probability distribution as the product of simple distributions, and the variational parameters are optimized by minimizing the KL divergence. Gibbs sampling is used to generate structured matching results. The state of each matching point pair is sampled from the conditional probability distribution one by one, and the configuration with the highest frequency after multiple sampling is taken as the final matching result. The temperature parameter of the attention mechanism in the geometric registration feature extraction and matching network is determined according to the signal-to-noise ratio, texture complexity and overlap degree of the fused image. The temperature parameter calculation formula is as follows: the temperature parameter is equal to the baseline temperature value multiplied by the signal-to-noise ratio standard value divided by the actual signal-to-noise ratio, and then multiplied by the texture complexity divided by the standard texture complexity, and then multiplied by the overlap degree divided by the standard overlap degree.

[0035] The step of establishing a training data set for the geometric registration feature extraction and matching network includes: collecting multi-modal remote sensing image pairs obtained by different sensors, ensuring that the image pairs contain various combinations of optical and SAR, optical and lidar, and SAR and lidar; manually labeling each image pair and marking at least 50 corresponding control points in the image, the control points being uniformly distributed and covering different ground object types; applying random geometric transformations to the image pair, including a rotation angle of -30 degrees to +30 degrees, a scaling ratio of 0.8 to 1.2, and a translation distance of -50 pixels to +50 pixels, to generate image pairs with different degrees of distortion; applying radiation changes to the image pair, including a brightness adjustment range of -20% to +20%, a contrast adjustment range of 0.8 to 1.2, and a Gaussian noise standard deviation of 5 to 15, to simulate different imaging conditions; combining the processed image pair with the corresponding control point label to form a training sample, each training sample containing a pair of images and the coordinate information of multiple matching point pairs; and dividing the training data set according to a ratio of 70% training set, 15% validation set, and 15% test set.

[0036] The step of training the geometric registration feature extraction and matching network includes: initializing the network weights, randomly initializing the convolution layer weights using the He initialization method, and initializing the bias term to zero; defining the loss function as the weighted sum of the contrast loss and the geometric consistency loss, the contrast loss measuring the similarity of the matching point pair features, and the geometric consistency loss measuring the spatial relationship rationality of the matching result, with weight coefficients set to 0.6 and 0.4; updating the parameters using the Adam optimizer, setting the initial learning rate to 0.001, and setting the learning rate decay strategy to decay to 0.5 times the original value every 20 rounds; calculating the loss value by forward propagation on the training set, calculating the gradient by the back propagation algorithm, and updating the network weights; evaluating the network performance on the validation set every 5 rounds of training, recording the matching accuracy and average matching error; stopping training when the validation set performance does not improve for 10 consecutive rounds, and selecting the model parameters with the best validation set performance as the final model; evaluating the generalization ability of the model on the test set, calculating the matching success rate and average re-projection error.

[0037] The structured prediction framework based on probabilistic graphical model converts the isolated point pair matching problem into a global optimization problem by explicitly modeling the spatial dependency between matching point pairs. Traditional matching methods evaluate each candidate matching point pair independently, ignoring the geometric consistency constraints satisfied by neighboring matching point pairs, which is prone to produce false matches in repetitive texture regions. The structured prediction framework based on probabilistic graphical model models the configuration space composed of all candidate matching point pairs as a probability graph, and realizes the propagation and integration of information in the global range through message passing algorithm, so that each matching decision depends not only on the local feature similarity, but also considers the consistency with neighboring matches. The variational inference technique introduces a variational distribution to approximate the complex posterior distribution, decomposes the difficult global inference problem into multiple local optimization sub-problems, greatly reduces the computational complexity while ensuring the quality of the solution. The Gibbs sampling explores the configuration space through a random sampling mechanism to avoid getting stuck in a local optimal solution, and generates a structured matching result that satisfies the global constraints.

[0038] The structured prediction framework based on probabilistic graphical model provides strong global reasoning ability for the geometric registration feature extraction and matching network, significantly reduces the false matching rate, and still maintains high matching accuracy in scenes with repetitive texture, seasonal changes, and shadow occlusion. In the entire metadata extraction scheme, accurate geometric registration is the basis for multi-modal data fusion. The structured prediction framework based on probabilistic graphical model ensures registration accuracy through structured prediction, avoids the influence of false features introduced by registration errors on subsequent feature classification and metadata identification, and provides a basis for fusion weight allocation and uncertainty quantification through probability output to provide matching confidence evaluation. The structured prediction framework based on probabilistic graphical model formulates the matching problem as a maximum a posteriori probability inference on a Markov random field, decomposes the joint probability distribution into the product of node potential functions and edge potential functions, encodes the local evidence of a single matching point pair in the node potential function, and encodes the geometric constraints between neighboring matching point pairs in the edge potential function. The belief propagation algorithm propagates probability information between graph nodes through iterative message passing. Each message represents the confidence evaluation of the sending node on the possible state of the receiving node. The message update rule combines the node potential function and the messages from other neighbors. After multiple iterations, the marginal probability of each node converges to an approximate optimal solution. The variational inference technique introduces a variational distribution to approximate the true posterior distribution, optimizes the variational parameters by minimizing the KL divergence between the variational distribution and the true distribution, and decomposes the complex joint inference problem into independent optimization of each node or small group, greatly reducing the computational complexity. The Gibbs sampling starts from a random initial configuration and updates each matching point pair according to its conditional probability distribution. The conditional probability distribution is determined by the node potential function and the current state of the neighboring nodes. After sufficient sampling iterations, the sample distribution converges to the joint posterior distribution, and the most likely matching configuration is obtained by statistical sampling results.

[0039] Siamese architecture is a kind of two-branch neural network structure, both branches have the same network structure and shared weight parameters, used to learn the similarity measure between pairs of input data. The loop consistency check refers to the matching points from image A to image B and the reverse matching points from image B to image A, if the two matching results point to the same position, it is considered as a reliable verification mechanism. RANSAC robust estimation is the abbreviation of random sample consensus algorithm, by randomly selecting the minimum sample set to fit the model, counting the number of inliers that meet the model, and iteratively selecting the model parameters with the most inliers and removing the outliers. Global elevation model refers to the digital elevation data product covering the global range, providing surface elevation information for geometric correction and terrain analysis. The matching search space refers to the possible position range of the candidate matching points in the image matching process, which reduces the search space through prior constraints to improve the matching efficiency and accuracy.

[0040] The adaptive filtering algorithm based on non-local mean and edge perception combines the filtering strategies of spatial domain and frequency domain. The adaptive filtering algorithm based on non-local mean and edge perception first detects the image edge through Sobel operator and Canny operator, calculates the edge intensity and edge direction of each pixel. For uniform areas with edge intensity lower than threshold 20, non-local mean filtering with a 15x15 window is used, searching for similar pixels in the window to the center pixel gray distribution for weighted average, and the weight is calculated according to the Euclidean distance between the pixel blocks, the smaller the distance, the greater the weight. For edge regions with edge intensity higher than threshold 20, directional filtering along the edge direction is used, the long axis direction of the filtering window is consistent with the edge direction, the long axis length is 11 pixels, and the short axis length is 3 pixels, only the pixels in the window are weighted and averaged. A multi-scale decomposition framework is introduced to decompose the image into three layers of wavelet, and the low-frequency approximation subband and high-frequency detail subband are obtained after decomposition, the low-frequency approximation subband contains the main structure information of the image, and the high-frequency detail subband contains noise and detail information. Large window filtering is applied to the low-frequency approximation subband, and threshold shrinkage is applied to the high-frequency detail subband, coefficients less than three times the noise standard deviation are set to zero, and then wavelet reconstruction is performed to obtain the filtered image set. Structure tensor is a second-order matrix that describes the local image structure, the matrix eigenvalue reflects the main direction and amplitude of local gray change, and the ratio of eigenvalues is used to judge whether the region is a uniform region, an edge region or a corner region.

[0041] An iterative registration optimization algorithm based on mutual information maximization is proposed to achieve accurate registration by maximizing the mutual information between two images. Mutual information measures the statistical dependence between two random variables and is defined as the sum of the entropies of the two variables minus the joint entropy. The larger the mutual information is, the higher the alignment degree of the two images is. The algorithm first calculates the gray level histograms of the reference image and the image to be registered, and then computes the marginal entropies. A two-dimensional joint histogram is constructed to count the frequency of the gray level combinations of the corresponding pixels in the two images and to calculate the joint entropy. The mutual information is equal to the sum of the entropies of the reference image and the image to be registered minus the joint entropy. The Powell direction acceleration method is used to optimize the registration parameters, including translation, rotation and scaling parameters. The Powell direction acceleration method performs one-dimensional search along a set of conjugate directions, and updates one parameter of each direction in each iteration. The convergence is accelerated by the property of the conjugate directions. The adaptive weighted network is a three-layer fully connected neural network, and the input is the feature statistics of different modal images in the current scene, including mean, variance, entropy and texture complexity. The output is the fusion weight of each modality. The adaptive weighted network is trained by supervised learning, and the training samples are the modal features and the optimal weights labeled by artificial annotation in different scenes. Feature-level fusion is to splice the feature vectors of different modal images by weighting, and decision-level fusion is to vote or weighted fusion of the independent classification results of each modality. The uncertainty quantification index is realized by calculating the entropy value of the fused image. The larger the entropy value is, the higher the classification uncertainty is.

[0042] The semi-supervised classification algorithm based on random walk utilizes the similarity relationship between a small number of labeled samples and a large number of unlabeled samples for classification. The semi-supervised classification algorithm based on random walk represents samples as nodes of a graph, and the edge weight between nodes represents the sample similarity, which is calculated by a Gaussian kernel function. The bandwidth parameter of the kernel function is determined according to the median of the distance between samples. The labeled samples are used as absorbing state nodes, and the unlabeled samples are used as transient nodes. Random walk starts from a transient node, moves to an adjacent node according to the transition probability of the normalized edge weight, and is finally absorbed by an absorbing state node. The arrival probability of each transient node to each absorbing state is calculated by solving a linear equation system. The class corresponding to the absorbing state with the maximum arrival probability is the predicted class of the transient node. The coefficient matrix of the linear equation system is the transition probability matrix between transient nodes, and the constant term is the transition probability of the transient node to the absorbing state. The conjugate gradient method or LU decomposition is used for solving. The restart random walk mechanism returns to the starting node with a probability of 0.15 in each step of walking, which prevents the walking from deviating too far from the initial position and improves the local consistency. The sparse similarity matrix is realized by retaining only the nearest neighbor connection of each node. For each node, only the edges of the 10 neighbor nodes with the largest similarity are retained, and the weights of other edges are set to zero, which greatly reduces the matrix storage demand and computational complexity. The active learning strategy selects samples according to the classification uncertainty. The entropy value of each unlabeled sample is calculated, and the sample with the maximum entropy value is selected for manual labeling and added to the labeled sample set. The result of the ground object classification is iteratively updated.

[0043] The random walk based semi-supervised classification algorithm fully utilizes the intrinsic manifold structure of data, realizes efficient classification through propagation of labeled information, and still obtains high classification accuracy in the scene of scarce labeled samples. The random walk based semi-supervised classification algorithm provides efficient ground feature classification capability for the whole scheme, realizes automatic classification of large-scale remote sensing images in the case of only a small number of labeled samples, and significantly reduces the cost of manual labeling. Through the random walk propagation mechanism, the random walk based semi-supervised classification algorithm can capture the local manifold structure of data, so that similar samples obtain consistent classification labels, and the spatial continuity of the ground feature classification result is improved. The restart random walk mechanism and the sparsification similarity matrix strategy greatly improve the calculation efficiency while ensuring the classification accuracy, so that the random walk based semi-supervised classification algorithm is suitable for processing high-resolution large-format remote sensing images. The active learning strategy selects the most representative samples for labeling, so as to obtain the maximum classification performance improvement with the least labeling cost, and realizes the optimal allocation of labeling resources. In the metadata extraction process, accurate ground feature classification result is the basis for identifying land cover type, vegetation type and building distribution metadata elements, and the efficiency and accuracy of the random walk based semi-supervised classification algorithm directly determine the quality and efficiency of metadata extraction. The random walk based semi-supervised classification algorithm models the classification problem as a label propagation process on a graph, constructs a connection graph using the similarity relationship between samples, and simulates the diffusion of label information from labeled samples to unlabeled samples through random walk. The arrival probability describes the association strength between unlabeled samples and labeled samples of each class, and the class with high arrival probability indicates that the unlabeled sample and the labeled sample of the class are close in feature space or are connected through a similar sample chain. The restart random walk mechanism introduces the probability of returning to the initial node, avoids the false association of distant samples caused by infinite diffusion of walk in the graph, and ensures the local consistency and global rationality of the classification result. The sparsification similarity matrix constructs a sparse graph by retaining the most similar neighbor connections, reduces the storage space and calculation time, and at the same time filters weakly related connections to reduce noise interference.

[0044] The standard persistence value refers to the average persistence value of stable terrain features statistically obtained in a large number of typical terrain samples, which is used to determine whether the extracted topological features are significant. The multi-scale filtering processing level refers to the number of times of Gaussian filtering of elevation data at different scales. Increasing the multi-scale filtering processing level extracts more scale terrain features but increases the amount of calculation. The spatial resolution refers to the minimum size of ground objects that can be identified in a remote sensing image, which is usually represented by the size of the image pixel. The geographic range refers to the geographic area covered by the remote sensing image, which is usually represented by the latitude and longitude coordinates of the four corners of the image. The coordinate system refers to the geographic coordinate system or the projection coordinate system used by the remote sensing image, including the geodetic datum, ellipsoid parameters, and projection method information. The acquisition time refers to the date and time when the remote sensing satellite or sensor acquires the image. The sensor parameters include the sensor type, band setting, spectral range, radiation resolution, and swath width technical indicators.

[0045] The signal-to-noise ratio refers to the ratio of signal power to noise power in a remote sensing image, which reflects the image quality. The ratio is estimated by the ratio of the mean value to the standard deviation of the uniform area of the image. The larger the ratio, the higher the signal-to-noise ratio. The texture complexity is calculated by the gray level co-occurrence matrix to obtain the entropy, contrast, and correlation texture features of the image. The comprehensive texture complexity index is obtained by weighted sum of multiple texture features. The overlap degree refers to the ratio of the overlapping area of two images to be registered to the area of the smaller image. The larger the overlap degree, the more information available for registration. The reference temperature value, signal-to-noise ratio standard value, standard texture complexity, and standard overlap degree are all average values of the training data set obtained in the training process of the geometric registration feature extraction and matching network, which are used for normalization calculation. The average persistence value is the arithmetic mean of all persistence values in the training sample. The connected component refers to a set of connected pixels in an image, and the gray level of the pixels in the connected component is greater than or equal to the current threshold. The super level set refers to a set of all pixels with a gray level greater than or equal to a specified threshold. The Markov random field is a probabilistic graph model in which the joint probability distribution of random variables satisfies the Markov property, and each variable is only directly related to its neighbor variables. The node potential function is a function defined on a single node, which measures the likelihood of the node taking a predetermined value. The edge potential function is a function defined on a pair of adjacent nodes, which measures the compatibility of adjacent nodes taking a predetermined value combination. The joint probability distribution is the probability of multiple random variables taking predetermined values simultaneously, which is obtained by normalization of the product of the node potential function and the edge potential function. The conditional probability distribution refers to the probability distribution of a variable given the values of other variables.

[0046] The application optionally further provides a computer-implemented mode to form an automatic metadata recognition and extraction system for remote sensing images and spatial data, wherein a readable storage medium is arranged in the computer, and program instructions are stored in the readable storage medium, and the program instructions perform the above method when running in the computer.

[0047] The specific implementation of the above steps is described in detail below.

[0048] The specific implementation of step S1 is to first receive multi-source remote sensing image data, which includes three types of optical images, SAR images and lidar data. The optical images are subjected to radiation calibration processing to convert the original digital quantization value into a radiation brightness value. The conversion process applies the gain coefficient and offset parameter provided by the sensor to eliminate the nonlinear characteristics of the sensor response. Then, the optical images are subjected to atmospheric correction, and the 6S radiation transfer model is used to calculate the scattering and absorption of atmospheric layer on solar radiation to obtain the true reflectivity of the ground surface. Atmospheric correction requires input of atmospheric parameters such as atmospheric visibility, aerosol optical thickness and water vapor content. The SAR images are subjected to multi-view processing to non-coherently superimpose the complex data of multiple observations to reduce the intensity of coherent speckle noise. The number of views for multi-view processing is usually selected between 4 and 16. Then, the SAR images are subjected to geographic encoding, and the slant range projection coordinates are converted into geographic coordinates using the range-doppler positioning model. The geographic encoding process needs to be combined with orbit parameters and elevation data to achieve geometric correction. The lidar data is subjected to point cloud filtering, and the progressive morphological filtering algorithm is used to separate ground points and non-ground points. The filtering process gradually increases the window size of morphological opening operation from 3x3 pixels to 21x21 pixels. The separated ground points are subjected to elevation normalization, and the elevation values of all points are reduced by the elevation of the ground point at the corresponding position to obtain the height value relative to the ground. Finally, the processed optical images, SAR images and lidar elevation data are unified under the WGS84 geographic coordinate system to complete the coordinate system conversion and spatial registration, forming a multi-modal image set. The purpose of the step is to eliminate the radiation difference and geometric distortion of different sensor data to provide a unified and standardized data basis for subsequent processing.

[0049] The specific implementation of step S2 is GPU accelerated processing of high-resolution images with resolution higher than 1 meter in the multi-modal image set. First, the high-resolution images are divided into blocks with a block size of 512x512 pixels, and a 64-pixel overlap region is set between adjacent blocks to avoid boundary effects. A quadtree spatial index is constructed, and the image space is recursively divided into four quadrants. Each quadrant node stores the minimum bounding rectangle, average gray value, variance, and block file path information of the corresponding region. Based on the quadtree index, an image pyramid is constructed, starting from the original image and down-sampling layer by layer. The number of rows and columns of each layer of the image is half of the previous layer, and the number of pyramid layers is set to 5. At each pyramid level, an adaptive sampling strategy is applied to calculate the gray variance of each block region. Blocks with a variance greater than 30 are marked as high complexity regions, and the sampling density is set to 1 point per 4 pixels. Blocks with a variance less than 30 are marked as low complexity regions, and the sampling density is set to 1 point per 16 pixels. The image blocks are distributed to multiple GPUs using the CUDA parallel computing framework. Each GPU is responsible for processing several consecutive image blocks. The block data is loaded into shared memory within the GPU, a large number of threads are started for parallel processing of the pixels within the block, and the results are written back to global memory after processing is complete. The purpose of the step is to accelerate high-resolution image processing using the large-scale parallel computing capabilities of GPUs, and to reduce computational complexity through block strategies and adaptive sampling.

[0050] The specific implementation of step S3 is to extract the elevation data in the multi-modal image set as input. Gaussian filtering, bilateral filtering and median filtering are applied to the elevation data in sequence for multi-scale filtering processing. The standard deviation parameters of Gaussian filtering are set to 1 pixel, 2 pixels and 4 pixels for three scales, the spatial domain standard deviation of bilateral filtering is set to 3 pixels, the value domain standard deviation is set to 10 meters, and the window size of median filtering is set to 5x5 pixels. The position with the minimum gray value is selected from the filtered elevation field as the initial threshold value, and the threshold value is gradually increased by 1 meter. At each threshold level, the pixels with a gray value greater than or equal to the threshold value form a super-level set. Connected component analysis is performed on the super-level set, and each connected component is labeled using a region growing algorithm. The threshold value at which each connected component first appears is recorded as the birth time, and the threshold value at which the connected component disappears or merges with other connected components is recorded as the death time. The persistence value is calculated as the difference between the death time and the birth time. The arithmetic mean of the persistence values of all connected components is calculated to obtain the average persistence value. A persistence barcode graph is constructed, with the horizontal axis representing the threshold range and the vertical axis representing the connected component number. Each connected component corresponds to a horizontal line segment, with the starting point being the birth time, the ending point being the death time, and the length being the persistence value. The zero-order Betti number sequence is calculated, and the zero-order Betti number is equal to the number of connected components at the current threshold level. The first-order Betti number sequence is calculated, and the first-order Betti number is equal to the number of independent loops. The connected components with a persistence value greater than twice the average persistence value are selected, and the topological feature markers corresponding to the connected components are marked as stable terrain markers. The persistence maxima points are extracted and connected to form a ridge line, and the persistence minima points are extracted and connected to form a valley line. The purpose of the step is to automatically identify significant features in the terrain using topological persistence theory, providing stable control information for geometric registration.

[0051] The specific implementation of step S4 is to input the image pairs in the multi-modal image set into a geometric registration feature extraction and matching network. The geometric registration feature extraction and matching network includes two branches with the same structure and weight sharing. Each branch sequentially passes through a first convolutional layer with a convolution kernel size of 7x7 and a channel number of 64, a second convolutional layer with a convolution kernel size of 5x5 and a channel number of 128, a third convolutional layer with a convolution kernel size of 3x3 and a channel number of 256, a fourth convolutional layer with a convolution kernel size of 3x3 and a channel number of 256, and a fifth convolutional layer with a convolution kernel size of 1x1 and a channel number of 512 to extract features. After each convolutional layer, a batch normalization layer and a ReLU activation function are connected. The output feature vectors of the two branches are calculated for similarity by a cosine similarity measurement layer. Feature point pairs with a cosine similarity greater than 0.85 are retained as candidate matching point pairs. A Markov random field is constructed, each candidate matching point pair is taken as a graph node, the distance change and the angle change between adjacent candidate matching point pairs are calculated, and an edge is connected between candidate matching point pairs with a distance change less than 5 pixels and an angle change less than 10 degrees. A node potential function is defined to be equal to the cosine similarity of the feature vector, and an edge potential function is defined to be equal to the negative exponential function of the distance change and the angle change. A belief propagation algorithm is used for message passing, each node sends a message to adjacent nodes, the message content is the confidence of the node under the matching or non-matching hypothesis, and the message is iteratively updated until the edge probability of each node converges, and the convergence criterion is that the probability change of two consecutive iterations is less than 0.001. The variational inference technique is used to introduce a variational distribution to approximate the posterior probability distribution, and the variational distribution is assumed to be a Bernoulli distribution independent of each node. The variational parameters are optimized by minimizing the KL divergence. Gibbs sampling is used to generate the matching result. Starting from a random initial configuration, each candidate matching point pair is sampled according to its conditional probability distribution, and the conditional probability is determined by the node potential function and the current state of adjacent nodes. After 1000 sampling iterations, the frequency of each candidate matching point pair being sampled as a matching state is counted, and the candidate matching point pair with a frequency greater than 0.5 is determined as the final matching point pair. RANSAC is used to remove false matching point pairs. Four matching point pairs are randomly selected to calculate an affine transformation matrix, the number of inliers that meet the affine transformation matrix is counted, and the affine transformation matrix with the most inliers is selected after 1000 iterations. Matching point pairs with a re-projection error greater than 2 pixels are removed. The purpose of the step is to realize high-precision multi-modal image registration by using deep learning and probabilistic graph models, and to provide accurate geometric correspondence for data fusion.

[0052] The specific implementation of step S5 is speckle noise suppression processing on the SAR image in the registered post-image set. First, the gradient amplitude of each pixel is calculated using the Sobel operator, and edge detection is performed using the Canny operator, with the high threshold set to 1.5 times the average gradient amplitude and the low threshold set to 0.4 times the high threshold. The edge strength of each pixel is calculated, defined as the number of edge pixels in the pixel and its 8-neighborhood. For uniform regions with edge strength lower than the threshold value 20, non-local mean filtering is performed using a 15x15 pixel window, and the Euclidean distance between the center pixel and the 5x5 neighborhood block of the search pixel is calculated, with smaller distances having larger weights. For edge regions with edge strength higher than the threshold value 20, the local edge direction is calculated using the structure tensor, and an elliptical directional filtering window with a major axis length of 11 pixels and a minor axis length of 3 pixels is used. The image is decomposed into three layers of wavelets, the mean filtering with a 9x9 pixel window is applied to the low-frequency approximation subband, and the threshold shrinkage is applied to the high-frequency detail subband, with coefficients less than three times the noise standard deviation set to zero, and then the wavelet reconstruction is performed to obtain the filtered post-image set. The purpose of the step is to suppress the speckle noise of the SAR image while preserving the edge and detail information.

[0053] The specific implementation of step S6 is multi-modal data fusion on the filtered post-image set. The gray level histogram and entropy value of the reference image and the image to be registered are calculated, a two-dimensional joint histogram is constructed and the joint entropy is calculated, and the mutual information is calculated as the entropy of the reference image plus the entropy of the image to be registered minus the joint entropy. The Powell direction acceleration method is used to optimize the registration parameters, including horizontal translation, vertical translation, rotation angle, horizontal scaling ratio, and vertical scaling ratio, and one-dimensional search is performed in each direction to find the step size that maximizes the mutual information, and iteration is performed until the increment of mutual information is less than the threshold value 0.001. The input of the adaptive weighted network is set to five features of different modal images, including mean, variance, information entropy, contrast, and correlation. The adaptive weighted network includes an input layer, two hidden layers, and an output layer, with the number of hidden layer neurons being 64 and 32 respectively, and the output layer uses the Softmax function to ensure weight normalization. When performing feature-level fusion, the feature vectors of each modal image are weighted and spliced according to the fusion weight, and when performing decision-level fusion, the classification results of each modal are weighted and voted according to the weight. The class probability entropy value of each pixel of the fused image is calculated as the uncertainty quantification index, and the pixels with entropy values greater than the threshold value 1.5 are marked as low confidence regions. The purpose of the step is to achieve accurate registration and adaptive fusion of multi-modal data.

[0054] The specific implementation of step S7 involves classifying ground features in the fused image. Each pixel is represented as a feature vector, including spectral features, texture features, and geometric features. A Gaussian kernel function is used to calculate the sample similarity between pixels, with the kernel function parameter set to the median of the Euclidean distance between the feature vectors. A sample similarity matrix is ​​constructed and sparsified, retaining only the connections between the 10 neighboring pixels with the highest sample similarity for each pixel. Manually labeled pixels are set as absorbing state nodes, and the remaining pixels are set as transient nodes. The sample similarity matrix is ​​normalized to obtain the transition probability matrix, and a system of linear equations is constructed. The conjugate gradient method is used to solve for the arrival probability of each transient node reaching each absorbing state. For each transient node, the category corresponding to the absorbing state with the highest arrival probability is selected as the predicted category. A restarted random walk mechanism is introduced, with each step returning to the starting node with a probability of 0.15. The classification confidence is calculated as the difference between the maximum and second-largest arrival probabilities. The 50 samples with the highest entropy values ​​are selected for manual labeling, and the classification results are iteratively updated. The purpose of these steps is to achieve automatic classification of large-scale remote sensing images using semi-supervised learning.

[0055] The specific implementation of step S8 involves extracting metadata information based on the land cover classification results and the stable terrain markers. The area and proportion of six land cover types—built-up land, farmland, forest land, grassland, water bodies, and bare land—are statistically analyzed and written into the land cover type field of the metadata information set. The total length, average elevation, and maximum elevation of ridgelines and valley lines are statistically analyzed and written into the terrain feature field. The pixel size, corner coordinates, coordinate system definition, acquisition timestamp, and sensor parameters of the fused image are read and written into the corresponding fields. When the proportion of pixels with a classification confidence score below 0.75 exceeds 5%, 100 pixels are selected for manual annotation and reclassification. When the proportion of topological features with a persistence value below 60% of the standard persistence value exceeds 10%, a multi-scale filtering processing layer is added, and terrain features are re-extracted. The purpose of these steps is to automatically generate complete and accurate metadata information by integrating multiple data sources.

[0056] It should be noted that one of the key technical ideas of this invention is a terrain feature extraction algorithm based on topological persistence. This algorithm quantifies the saliency of terrain features into persistence values ​​through persistence homology theory, avoiding the problem of traditional methods relying on manually set thresholds. The persistence barcode image visually displays the lifecycle of terrain features at different scales. Long bars correspond to stable primary terrain structures, while short bars correspond to noise or minor undulations. Noise is automatically filtered out using persistence thresholds, achieving robust extraction of terrain features. The algorithm's advantage in high-resolution remote sensing image processing lies in the fact that high-resolution images contain rich terrain details. Traditional edge detection or gradient calculation methods are easily affected by noise, while persistence analysis, through multi-scale comprehensive evaluation, filters out truly stable terrain features, providing reliable control information for geometric registration and reducing geometric distortion correction errors caused by terrain undulations.

[0057] The second key technical idea of ​​this invention is a structured prediction framework based on probabilistic graphical models, which transforms the isolated point-pair matching problem into a global optimization problem. Traditional matching methods independently evaluate each candidate matching point pair, ignoring the geometric consistency constraints satisfied by adjacent matching points, which easily leads to mismatches in areas with repetitive textures. The proposed framework models the configuration space composed of all candidate matching point pairs as a Markov random field, and uses a message passing algorithm to achieve the propagation and integration of information globally, so that each matching decision not only depends on local feature similarity, but also considers consistency with adjacent matches. Variational inference technology introduces a variational distribution to approximate a complex posterior distribution, decomposing the computationally difficult global inference problem into multiple local optimization sub-problems, significantly reducing computational complexity while ensuring the quality of the solution. Gibbs sampling explores the configuration space through a random sampling mechanism to avoid getting trapped in local optima.

[0058] The third key technical idea of ​​this invention is a semi-supervised classification algorithm based on random walks, which classifies samples by utilizing the similarity relationship between a small number of labeled samples and a large number of unlabeled samples. The algorithm models the classification problem as a label propagation process on a graph, simulating the diffusion of label information from labeled samples to unlabeled samples through random walks. This fully utilizes the inherent manifold structure of the data, ensuring that similar samples receive consistent classification labels. Restarting the random walk mechanism introduces the probability of returning to the initial node, preventing the walk from spreading infinitely across the graph and causing false associations between distant samples, thus ensuring the local consistency and global reasonableness of the classification results. The sparse similarity matrix constructs a sparse graph by retaining the most similar neighbor connections, reducing storage space and computation time. The active learning strategy optimizes the allocation of labeling resources by selecting the most representative sample labels.

[0059] The synergistic effect of the three key technical approaches mentioned above lies in the fact that the stable terrain markers extracted by the topological persistence algorithm provide reliable spatial constraints for the probabilistic graphical model, significantly improving the accuracy and robustness of geometric registration. Accurate geometric registration lays the foundation for multimodal data fusion, avoiding the impact of pseudo-features introduced by registration errors on subsequent processing. The matching confidence output by the probabilistic graphical model provides a basis for fusion weight allocation and uncertainty quantification, improving the quality of the fused image. High-quality fused imagery provides rich feature information for semi-supervised classification algorithms, improving classification accuracy. Accurate land cover classification results are the foundation for identifying metadata elements such as land cover type, vegetation type, and building distribution. The synergistic effect of these three approaches forms a complete technical chain from terrain feature extraction to geometric registration, multimodal fusion, land cover classification, and metadata identification. Compared with the independent processing of each stage in traditional methods, collaborative optimization significantly improves the overall processing accuracy and efficiency.

[0060] It should be noted that this invention also solves the following technical problem: In existing technologies, remote sensing image land cover classification requires a large number of manually labeled samples for supervised learning training, which is costly and time-consuming, especially in high-resolution, large-format remote sensing image processing, where full-coverage labeling is almost infeasible. This invention applies a semi-supervised classification algorithm based on random walks, using the similarity relationship between a small number of labeled samples and a large number of unlabeled samples for classification. Samples are represented as nodes in a graph, and classification labels are determined by calculating the arrival probability of unlabeled samples reaching each labeled sample during a random walk. This fully utilizes the inherent manifold structure of the data, achieving high classification accuracy even in scenarios with scarce labeled samples, significantly reducing the cost of manual labeling. An active learning strategy is introduced to select the most representative samples for manual labeling based on classification uncertainty, achieving maximum classification performance improvement with minimal labeling cost, optimizing the allocation of labeling resources, and solving the technical problem of insufficient labeled samples in remote sensing image land cover classification.

[0061] Furthermore, in existing technologies, the contribution of different modalities in multimodal remote sensing data fusion varies significantly across different scenarios. Fixed-weight fusion strategies struggle to adapt to scene changes, leading to unstable fused image quality. This invention employs an iterative registration optimization algorithm based on maximizing mutual information to achieve multimodal data fusion. An adaptive weighted network learns the contribution of different modalities in different scenarios. The adaptive weighted network inputs feature statistics of different modalities, including mean, variance, entropy, and texture complexity. Supervised learning is used to train and output the fusion weights for each modality. A hybrid strategy combining feature-level fusion and decision-level fusion is combined to calculate the uncertainty quantification index of the fused image to evaluate the fusion quality. This achieves scene-adaptive data fusion, improving the stability and reliability of the fused image in different application scenarios.

[0062] Specifically, the principle of this invention is as follows: The terrain feature extraction algorithm based on topological persistence calculates the birth and death times of connected components under different grayscale thresholds, constructing a persistent barcode map to quantify the lifecycle of terrain features. Topological features with high persistence values ​​correspond to stable primary terrain structures, while those with low persistence values ​​correspond to noise or minor undulations. By setting a persistence threshold, the algorithm automatically distinguishes between real terrain features and noise interference, mathematically ensuring the robustness of the extraction results. The algorithm utilizes the principle of topological invariance to ensure the stability of extracted terrain features under geometric transformations such as image rotation and scaling, providing reliable control points for multimodal image registration. This invention, based on a probabilistic graphical model-based structured prediction framework, transforms the isolated feature point matching problem into a global optimization problem. It explicitly models the spatial dependencies between adjacent matching point pairs using Markov random fields, and uses a confidence propagation algorithm to ensure that each matching decision simultaneously considers local feature similarity and the geometric consistency of adjacent matches. Through multi-round message passing, information is propagated and integrated globally, ensuring that the matching results satisfy global geometric constraints from a probabilistic inference perspective, thus improving registration accuracy.

[0063] The following provides a specific embodiment 1 of the present invention. The specific implementation methods of steps S1, S2, S5 and S8 in this embodiment 1 are the same as those described above, and will not be repeated in detail here. The specific implementation methods of other steps are described in detail below.

[0064] The specific implementation of step S3 is as follows: A topological persistence-based terrain feature extraction algorithm is applied to perform multi-scale filtering on the elevation data in the multimodal image set. First, Gaussian filtering, bilateral filtering, and median filtering are applied sequentially to the elevation data to generate elevation fields with different smoothness levels. Then, starting from the lowest grayscale value, the threshold is gradually increased. At each threshold level, the connected components of the hyperlevel set are calculated, and the threshold at which each connected component first appears is recorded as the birth time. The threshold at which a connected component disappears or merges with other connected components is recorded as the time of disappearance. Durability value The calculation formula is expressed as follows:

[0065] ;

[0066] In the formula, The value is a persistence value, dimensionless; The time of extinction is expressed in gray levels; The birth time is expressed in gray levels. The maximum gray level is typically set to 255, with the unit being gray levels. Then, the zero-order Betti number is calculated for different thresholds. With first-order Betty number Zero-order Betty number The number of connected components is a positive integer; a first-order Betti number. The number of independent loops is represented by a positive integer. A persistent barcode map is constructed by calculating the changes in the Betty number under different thresholds. Cluster analysis is performed on the persistent values. The discriminant formula for the recognition conditions of stable terrain markers is expressed as follows:

[0067] ;

[0068] In the formula, The average persistence value is dimensionless and is obtained by taking the arithmetic mean of all persistence values ​​in the training samples. The normalization coefficient has an empirical value of 1.0 and is dimensionless. Ridge extraction is achieved by connecting persistent maxima, and valley extraction is achieved by connecting persistent minima. When the persistence value is lower than the standard persistence value... When the value reaches 60%, a multi-scale filtering layer is added. The discrimination formula is expressed as follows:

[0069] ;

[0070] In the formula, The standard persistence value is dimensionless and obtained through statistical analysis of a large number of typical terrain samples.

[0071] The specific implementation of step S4 is as follows: a geometric registration feature extraction and matching network is constructed to accurately register the multimodal image set. This network adopts a two-branch Siamese architecture, with each branch containing five convolutional layers. The kernel sizes are 7×7, 5×5, 3×3, 3×3, and 1×1, and the number of channels are 64, 128, 256, 256, and 512, respectively. Each convolutional layer is followed by a batch normalization layer and a modified linear unit activation function. The two branches share weight parameters to ensure the consistency of feature extraction. After feature extraction, the similarity between feature vectors is calculated through a cosine similarity measurement layer. The calculation formula is expressed as follows:

[0072] ;

[0073] In the formula, Cosine similarity, dimensionless; This is the first feature vector, with a dimension of 512; This is the second feature vector, with a dimension of 512; This is a vector dot product operation; For feature vectors The Euclidean norm is calculated using the following formula: ,in For vectors The One component; For feature vectors The Euclidean norm is calculated using the following formula: ,in For vectors The Each component has a similarity greater than 0.85, and feature point pairs are considered matching points. The network introduces a structured prediction framework based on a probabilistic graphical model, modeling the spatial relationship of matching point pairs as a Markov random field, with node potential functions. Defined as feature similarity, edge potential function Defined as the change in distance and angle between matching point pairs, the joint probability distribution The calculation formula is expressed as follows:

[0074] ;

[0075] In the formula, It is a joint probability distribution, dimensionless; Configuration for all matching point pairs; The normalization constant is dimensionless and is obtained by summing the products of potential functions over all possible configurations. The calculation formula is as follows: ; For nodes The potential function is dimensionless and its value is equal to the feature similarity of the matching point pair corresponding to that node. For nodes With nodes The boundary potential function between them is dimensionless, and its calculation formula is as follows: ,in For matching point pairs and The change in distance between them, in pixels. For matching point pairs and The change in angle between them, in radians. is the standard deviation parameter of the edge potential function, with an empirical value of 10, in pixels; For nodes The state indicates the matching hypothesis selected by that node; For nodes The state of a node represents the matching hypothesis it has chosen. A confidence propagation algorithm is used to propagate probability information between graph nodes. Each node sends a message to its neighbors, containing its confidence level under different matching hypotheses. After multiple iterations, the marginal probabilities of the nodes converge. Variational inference is used to approximate the posterior probability distribution as a product of simple distributions. The variational parameters are optimized by minimizing the Körbeklebler divergence. Gibbs sampling is used to generate structured matching results. The state of each matching point pair is sampled sequentially from the conditional probability distribution. After multiple samplings, the configuration with the highest statistical frequency is used as the final matching result. The matching search space is constrained by a global elevation model. Bidirectional verification is performed through a cyclic consistency test. Mismatched points are eliminated using robust estimation of random sampling consistency. The registered image set is output, including the temperature parameter of the attention mechanism. The calculation formula is expressed as follows:

[0076] ;

[0077] In the formula, This is a temperature parameter, dimensionless. The reference temperature value is 1.0, which is dimensionless. The signal-to-noise ratio standard value is dimensionless and is obtained through statistics from the training dataset. The actual signal-to-noise ratio is dimensionless and is estimated by the ratio of the mean to the standard deviation of the uniform region of the image. The texture complexity is dimensionless and is obtained by weighted summation of the image's entropy, contrast, and correlation texture features calculated using the gray-level co-occurrence matrix. The standard texture complexity is dimensionless and is obtained through statistics from the training dataset. The overlap is dimensionless and is the ratio of the overlapping area of ​​the two images to be registered to the area of ​​the smaller image. The standard overlap is dimensionless and is obtained through statistics from the training dataset.

[0078] The specific implementation of step S6 is to use an iterative registration optimization algorithm based on maximizing mutual information to achieve multimodal data fusion of the filtered image set, and to achieve mutual information fusion. The calculation formula is expressed as follows:

[0079] ;

[0080] In the formula, Mutual information, dimensionless; The entropy of the reference image, in bits; The entropy of the image to be registered is expressed in bits. Joint entropy, in bits; The maximum possible entropy value, in bits, is calculated using the following formula: ,in This represents the number of gray levels, typically set to 256. Among these, entropy... The calculation formula is expressed as follows:

[0081] ;

[0082] In the formula, The entropy of the reference image, in bits; This represents the number of gray levels, typically set to 256. gray levels in the reference image The probability of occurrence is dimensionless, determined by statistical analysis of gray levels in the image. The number of pixels is obtained by dividing the total number of pixels; This is a grayscale index, with values ​​ranging from 0 to 255. Joint entropy. The calculation formula is expressed as follows:

[0083] ;

[0084] In the formula, Joint entropy, in bits; The reference image has a gray level of And the gray level of the image to be registered is The joint probability, dimensionless, is obtained by statistically analyzing the gray-level combinations of corresponding pixels in two images. The number of pixel pairs is obtained by dividing the total number of pixels; The grayscale index of the image to be registered ranges from 0 to 255. The Powell orientation acceleration method is used to optimize the registration parameters, including translation, rotation, and scaling parameters. An adaptive weighted network is used to learn the contribution of different modalities in different scenes. A hybrid strategy combining feature-level fusion and decision-level fusion is employed to calculate the uncertainty quantification index of the fused image, and the fused image is output.

[0085] The specific implementation of step S7 is as follows: a semi-supervised classification algorithm based on random walk is applied to classify ground features in the fused image. A small number of labeled samples are used as the absorption states of the random walk, and a sample similarity matrix is ​​constructed as the transition probability. The calculation formula, derived using the Gaussian kernel function, is as follows:

[0086] ;

[0087] In the formula, For the sample With sample The similarity between them is dimensionless; For the sample The feature vectors, whose dimensions are determined by the sum of the number of bands and the number of texture features in the fused image; For the sample Feature vectors, dimensions and same; This is the Euclidean distance operator, and the calculation formula is: ,in For the feature vector dimension, For the sample The eigenvector of the first One portion, For the sample The eigenvector of the first One component; The kernel function bandwidth parameter is dimensionless and is determined based on the median distance between samples. The maximum Euclidean distance between samples, with the same dimensions as the eigenvector components, is used for normalization. The arrival probability of unlabeled samples reaching each absorption state is calculated by solving a system of linear equations. The class corresponding to the absorption state with the highest arrival probability is the predicted class of the unlabeled sample. The linear equation system is in the form of:

[0088] ;

[0089] In the formula, It is an identity matrix, dimensionless; The transition probability matrix between transient nodes is dimensionless, and its elements are... Indicates from transient node Transfer to transient node The probability of; This is the arrival probability vector of the transient node reaching each absorption state, which is dimensionless; Let be the dimensionless vector representing the transition probability from the transient node to the absorption state. The vector elements are... Indicates from transient node The probability of directly transitioning to the absorption state. A restart random walk mechanism is introduced, returning to the starting node with a probability of 0.15 in each step of the walk to prevent the walk from deviating too far from the initial position and improve local consistency. The sparsity of the similarity matrix is ​​achieved by retaining only the nearest neighbor connections of each node. For each node, only the edges of the 10 neighboring nodes with the highest similarity are retained, and the weights of other edges are reset to zero. The land cover classification results and classification confidence are output.

[0090] To better understand and implement this invention, a specific application scenario of the invention is provided below as Example 2: To verify the metadata extraction capability of this invention in complex mountainous environments, technicians selected a typical mountainous area as the test region. This area has drastically undulating terrain, including various landform types such as ridges, valleys, and steep slopes. The test dataset includes Landsat-8 optical imagery, Sentinel-1 SAR imagery, and airborne lidar point cloud data. Technicians first performed radiometric calibration on the optical imagery, converting the digital quantization values ​​into apparent reflectance. Then, they used the 6S radiative transfer model for atmospheric correction, setting the aerosol optical thickness to 0.15 and the water vapor content to 2.5. The SAR imagery was processed using a four-view method, and then the slant range projection was converted to a geographic coordinate system using a range-Doppler positioning model. A progressive morphological filtering algorithm was used to extract ground points from the lidar point cloud data. The morphological filtering window size was gradually increased from 3m to 21m, with a window size growth factor set to 1.5 and a maximum window size of 21m. The elevation difference threshold was gradually increased from 0.2m to 3.5m. This resulted in the extraction of 3.28 million ground points and 1.42 million non-ground points.

[0091] Technicians input the preprocessed multimodal image set into a GPU-accelerated block-parallel processing architecture, dividing the image into 256×256 pixel blocks, with each block assigned to a GPU thread group for processing. A quadtree spatial index was used to construct an image pyramid, containing five layers: layer 0 (original image resolution 10m), layer 1 (20m), layer 2 (40m), layer 3 (80m), and layer 4 (160m). In ridge areas with dramatic texture changes, the sampling density was set to 1000 sampling points per square kilometer, while in flat areas with uniform grayscale, the sampling density was reduced to 200 sampling points per square kilometer. Technicians applied a topologically persistent terrain feature extraction algorithm to process the elevation data, first performing multi-scale processing including Gaussian filtering, bilateral filtering, and median filtering. Starting from a minimum elevation value of 850m, the threshold was gradually increased to a maximum elevation value of 2650m, with a threshold step size of 10m. At each threshold level, the connected components of the hyperlevel set were calculated. The average persistence value of all topological features was calculated to be 85m. Features with persistence values ​​greater than 170m were identified as major ridge and valley structures. Technicians extracted 15 major ridge lines with a total length of approximately 78 kilometers and 8 major valley lines with a total length of approximately 52 kilometers.

[0092] Technicians constructed a geometric registration feature extraction and matching network to accurately register multimodal image sets. This network adopts the Siamese architecture, such as... Figure 2As shown, the dataset is divided into three parts: 70% training set, 15% validation set, and 15% test set. The training set contains 2450 pairs, the validation set 525 pairs, and the test set 525 pairs. Table 1 shows the sample distribution of different modality combinations in the training dataset, indicating that the optical and SAR combinations have the highest proportion. Technicians initialized the convolutional layer weights using the He initialization method, setting the weight coefficients for contrastive loss and geometric consistency loss to 0.6 and 0.4, respectively. The Adam optimizer was used for parameter updates, with an initial learning rate of 0.001, which decayed to 0.5 times its original value every 20 rounds. After 80 rounds of training, the validation set performance reached its optimal level, with a matching accuracy of 92.7% and an average matching error of 1.8 pixels. Technicians calculated the temperature parameters of the attention mechanism. The actual signal-to-noise ratio was 18.5, the standard signal-to-noise ratio was 20.0, the texture complexity was 0.82, the standard texture complexity was 0.75, the overlap was 0.88, the standard overlap was 0.85, and the baseline temperature was 1.0. The calculated temperature parameter was 1.18. The network output a total of 1847 matching point pairs. After the cycle consistency test, 1623 were retained. After RANSAC robust estimation to remove false matches, 1456 reliable matching point pairs were finally retained.

[0093] Table 1. Sample distribution of different modality combinations in the training dataset

[0094]

[0095] Technicians applied an adaptive filtering algorithm based on nonlocal means and edge awareness to suppress speckle noise in the registered SAR image. First, image edges were detected using the Sobel and Canny operators, with an edge intensity threshold set to 20. For uniform regions, a 15×15 window nonlocal means filter was used, with a search window size of 31×31, a similarity window size of 7×7, and a smoothing parameter set to 0.8. For edge regions, directional filtering was used, with a filter window length of 11 pixels on the major axis and 3 pixels on the minor axis. The image was then subjected to three-level wavelet decomposition using the Daubechies wavelet basis, yielding low-frequency approximate subbands and high-frequency detail subbands. Figure 3 As shown, technicians plotted the energy distribution of wavelet coefficients at different scales. The first high-frequency subband accounted for 28% of the energy, the second high-frequency subband accounted for 15%, the third high-frequency subband accounted for 8%, and the low-frequency approximate subband accounted for 49%. Threshold shrinkage was applied to the high-frequency detail subband, with a noise standard deviation estimate of 6.2. Coefficients with a standard deviation less than 18.6 were set to zero, and then wavelet reconstruction was performed to obtain the filtered image.

[0096] Technicians employed an iterative registration optimization algorithm based on maximizing mutual information to achieve multimodal data fusion. The grayscale histograms of the reference image and the image to be registered were calculated, with 256 grayscale levels. The reference image entropy was 7.32, the image to be registered entropy was 6.95, the initial joint entropy was 12.18, and the initial mutual information was 2.09. The Powell orientation acceleration method was used to optimize the registration parameters. The translation parameter search range was -100 to +100 pixels, the rotation parameter search range was -10 degrees to +10 degrees, and the scaling parameter search range was 0.95 to 1.05. After 15 iterations, the mutual information converged to 2.87, and the registration parameters were: a horizontal translation of 5.3 pixels and a vertical translation of -2.7 pixels, a rotation angle of 1.2 degrees, and a scaling ratio of 1.01. Technicians constructed an adaptive weighted network to learn the fusion weights for different modalities. The network input consisted of optical imagery with a mean of 0.35, variance of 0.082, entropy of 7.32, and texture complexity of 0.82; SAR imagery with a mean of 0.28, variance of 0.095, entropy of 6.95, and texture complexity of 0.76; and lidar elevation data with a mean of 1450m and a variance of 158000. The entropy is 5.63 and the texture complexity is 0.45. As shown in Table 2, the fusion weight distribution of each modality in different scenarios shows that the LiDAR data has a higher weight in mountainous scenarios. Technicians calculated the uncertainty quantification index of the fused image, and the entropy value of the fused image is 3.25, indicating that the classification uncertainty is at a moderate level.

[0097] Table 2. Fusion weights of each modality in different scenarios

[0098]

[0099] Technicians applied a semi-supervised classification algorithm based on random walks to classify land features in the fused imagery. 420 labeled samples were selected, including 120 forests, 85 grasslands, 95 bare lands, 65 water bodies, and 55 buildings. A sample similarity matrix was constructed, and a Gaussian kernel function was used to calculate the similarity between samples. The kernel function bandwidth parameter was determined to be 35.8 based on the median distance between samples. The samples were represented as nodes in a graph, containing a total of 328,000 nodes, of which 420 were absorbing nodes and 327,580 were transient nodes. For example... Figure 4 As shown, the technicians plotted the arrival probability of the random walk as a function of the number of iterations. The arrival probability converged after 12 iterations. The conjugate gradient method was used to solve the linear equation system, and the iteration error threshold was set to... After 45 iterations, the solution was completed. The random walk probability was restarted with a value of 0.15, and the sparsified similarity matrix retained 10 nearest neighbors for each node, achieving a sparsity of 99.997%. The classification results showed that forest area accounted for 38.5%, grassland area for 27.2%, bare land area for 18.6%, water area for 9.8%, and building area for 5.9%. Technicians calculated the classification confidence scores for each category: forest 0.87, grassland 0.82, bare land 0.73, water body 0.91, and building 0.68. Since the classification confidence scores for bare land and buildings were below the threshold of 0.75, the technicians initiated an active learning strategy, manually labeling the 50 samples with the highest entropy values. The labeled samples were then added to the labeled sample set, updating the land cover classification results. The confidence score for bare land increased to 0.79, and the confidence score for buildings increased to 0.77.

[0100] like Figure 5 As shown, technicians created a comparison chart of metadata feature extraction accuracy, demonstrating the extraction accuracy of different metadata features. Sensor parameters included: Landsat-8 with 11 bands, a spectral range of 0.43 to 12.51 micrometers, a radiometric resolution of 12 bits, and a swath width of 185 kilometers; Sentinel-1 with a center frequency of 5.405 GHz, polarization modes VV and VH, a radiometric resolution of 16 bits, and a swath width of 250 kilometers; and the airborne lidar with a laser wavelength of 1064 nm, a scanning frequency of 70 Hz, a point cloud density of 15 points per square meter, and a vertical accuracy of 0.15 m. Because the persistence value for some ridge areas was 135 m, which is 60% lower than the standard persistence value of 170 m (102 m), technicians added multi-scale filtering layers, adding Gaussian filters at scales 8.0 and 16.0, and re-extracted topological features, improving the persistence value to 148 m, successfully identifying topographic landmarks in the area.

[0101] It should be noted that the variables involved in this invention are explained in detail in Table 3.

[0102] Table 3. Variable Explanation Table

[0103]

[0104] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in the present invention should be included within the scope of protection of the present invention.

Claims

1. An automated metadata identification and extraction method for remote sensing images and spatial data, characterized in that, Multi-source remote sensing image data was collected and preprocessed to establish a multimodal image set under a unified geographic coordinate system. A GPU-accelerated block-parallel processing architecture was used to perform hierarchical processing on the multimodal image set to construct an image pyramid. A topological persistence-based terrain feature extraction algorithm was applied to perform multi-scale filtering on the elevation data in the multimodal image set and to calculate the birth and death times of connected components under different grayscale thresholds to construct a persistent barcode image for identifying stable terrain markers. A geometric registration feature extraction and matching network was constructed, and a structured prediction framework based on a probabilistic graphical model was used to model the spatial relationship of matching point pairs as a Markov random field. A confidence propagation algorithm was used to achieve globally optimal matching for accurate registration of the multimodal image set. After registration, an adaptive filtering algorithm based on nonlocal mean and edge perception was applied to the image set to suppress speckle noise. An iterative registration optimization algorithm based on maximizing mutual information, combined with an adaptive weighted network, was used to achieve multimodal data fusion. A semi-supervised classification algorithm based on random walk was applied to classify ground features in the fused image. Based on the ground feature classification results and stable terrain markers, metadata elements were automatically identified to form a metadata information set.

2. The method according to claim 1, characterized in that, In the preprocessing steps of multi-source remote sensing image data, radiometric calibration and atmospheric correction are performed on optical images, multi-view processing and geocoding are performed on synthetic aperture radar images, and point cloud filtering and elevation normalization are performed on lidar data.

3. The method according to claim 2, characterized in that, The GPU-accelerated block-parallel processing architecture uses a quadtree spatial index to build an image pyramid, applies an adaptive sampling strategy at each level of the image pyramid, and distributes computing tasks among multiple GPUs through the CUDA parallel computing framework.

4. The method according to claim 3, characterized in that, In the topological persistence-based terrain feature extraction algorithm, the persistence value is defined as the difference between the time of extinction and the time of birth. Topological features with persistence values ​​greater than twice the average persistence value are identified as stable terrain markers.

5. The method according to claim 4, characterized in that, The topological persistence-based terrain feature extraction algorithm constructs a persistent barcode map by calculating the changes in Betti number under different thresholds, and describes the connectivity features of the topological space through the Betti number sequence.

6. The method according to claim 5, characterized in that, The geometric registration feature extraction and matching network adopts a two-branch Siamese architecture, with the two branches sharing weight parameters. After feature extraction, the similarity between feature vectors is calculated through a cosine similarity measurement layer.

7. The method according to claim 6, characterized in that, In the structured prediction framework based on probabilistic graphical models, the node potential function is defined as the feature similarity, and the edge potential function is defined as the change in distance and angle between matching point pairs. The global optimal matching is achieved by maximizing the joint probability distribution.

8. The method according to claim 7, characterized in that, The confidence propagation algorithm propagates probability information between graph nodes. Each node sends a message to its neighboring nodes, and the message content is the confidence of the node under different matching assumptions. After multiple iterations, the marginal probabilities of the nodes converge.

9. The method according to claim 8, characterized in that, Geometric registration feature extraction and matching network are combined with global elevation model to constrain the matching search space. Bidirectional verification is performed through cycle consistency test, and RANSAC robust estimation is used to eliminate mismatched points.

10. The method according to claim 9, characterized in that, The adaptive filtering algorithm based on nonlocal mean and edge awareness calculates local directionality using a structure tensor, employs large-window smoothing in uniform regions, and directional filtering in edge regions.

Images