A multimodal remote sensing image matching method based on coupled phase structure and depth features

Abstract

The present invention relates to the technical field of remote sensing image matching, and in particular to a multimodal remote sensing image matching method that couples phase structure and depth features. The method uses a phase structure feature detection algorithm to detect feature points of a multimodal remote sensing image; a multi-level feature interaction and aggregation network is constructed to extract deep feature descriptors, and the multi-level feature interaction and aggregation network includes a dual-branch backbone network, a multi-level feature interaction and aggregation module based on a coordinate attention mechanism, and a dense feature fusion module based on a coordinate attention mechanism. A data set is constructed to train a network model, and the trained network model is used to extract deep feature descriptors from an image; then the deep feature descriptors are matched, and a multi-level mismatch elimination strategy is used to eliminate mismatches to obtain the final matching result. The present invention can combine local and global semantic information, highlight the expression of key features, and improve the number of correct matching points and the positioning accuracy of matching points.

Description

Technical Field

[0001] The present invention relates to the technical field of remote sensing image matching, and in particular to a multi-modal remote sensing image matching method of coupled phase structure and depth features. Background Art

[0002] With the rapid development of remote sensing technology, obtaining the required information from remote sensing images has become a very important means of information acquisition. Earth observation images from multiple types of sensors such as visible light, near infrared, and synthetic aperture radar (SAR) are becoming increasingly abundant. Heterogeneous remote sensing images obtained by sensors on different platforms have certain complementarity, providing a massive data source for in-depth mining of remote sensing information and analysis of big data. The combination of remote sensing images is widely used in tasks such as change detection, image registration, image fusion, image stitching, and 3D reconstruction, and image matching is the prerequisite for these applications. Image matching is a key link in the preprocessing of remote sensing information, and the matching accuracy has a very important impact on many subsequent remote sensing applications. Due to the differences in imaging methods, viewing angles, resolutions, and phases, multimodal remote sensing images have complex geometric distortions and significant nonlinear radiation distortions on the images, making multimodal remote sensing image matching still very challenging.

[0003] To solve the above problems, researchers have proposed many image matching methods. Multimodal image matching methods are mainly divided into two categories: traditional matching methods and deep learning-based matching methods. Traditional matching methods mainly use the intensity and gradient of images for matching. Intensity-based methods directly match images at the pixel level. These methods use the pixel values ​​of images for direct calculation. When the intensity values ​​of images differ greatly, these methods cannot achieve the expected results. They are also extremely sensitive to illumination and noise, and have poor robustness to nonlinear radiation differences. Gradient-based matching methods are also called feature-based matching methods, but artificially designed feature descriptors have certain limitations, insufficient representation capabilities, and lack robustness to geometric distortion and significant nonlinear radiation differences of images. Convolutional neural networks are introduced into image matching as high-level feature extractors of images. Compared with traditional feature extraction methods, deep convolutional neural networks can extract high-level semantic information of images. The deep learning method is applied to multimodal remote sensing image matching tasks in order to learn the semantic information of high-level features and achieve feature robustness matching. Existing deep learning methods mainly use deep convolutional neural networks to construct feature descriptors. Compared with manually designed descriptors, feature descriptors based on learning can learn more valuable information hidden in the data through convolutional neural networks. Existing multimodal remote sensing image matching methods based on deep learning have poor adaptability to complex geometric distortions and significant nonlinear radiation distortions. The existing single network model is difficult to extract robust features of multimodal remote sensing images. It is difficult to design a network to extract local and global features of the image and realize feature interaction and aggregation. At the same time, there is a lack of public multimodal remote sensing image training datasets, which limits the feature representation and performance of the model. It is crucial to obtain diverse training data for training robust and generalizable network models. How to effectively use deep learning methods to achieve high-quality and high-precision image matching has become a key and valuable issue in the field of image matching.

[0004] The invention patent application with publication number CN118072053A discloses a multimodal remote sensing image matching method, which performs step-by-step Gaussian filtering and downsampling on the input raw image data to construct a scale space; obtains multi-directional convolution results based on the linear synthesis of directional adjustable filtering; introduces a multi-scale feature response aggregation model, aggregates multi-directional feature information based on the moment analysis equation, generates a maximum moment value map, and detects significant feature points; the invention method improves the robustness of the descriptor, but also increases the complexity of the calculation. The invention patent application with publication number CN116597183A discloses a multimodal image feature matching method based on spatial and channel dual-dimensional attention. The method designs an attention mechanism that includes spatial and channel dual-dimensional attention, which can simultaneously extract features from spatial and channel dual dimensions to improve the model's information interaction ability between channels, help the model pay more attention to the positioning of important features, and improve matching efficiency. However, the invention method has poor matching effect in the face of significant nonlinear radiation differences and geometric distortion. The invention patent with the publication number CN112149728A discloses a fast multimodal image template matching method, which first uses dense gradient features to construct a multidimensional structural feature descriptor, and then performs feature dimensionality reduction on the descriptor to form an effective low-dimensional feature descriptor, namely the L-CMOG descriptor, which can effectively characterize the common attributes between multimodal images, such as structure, contour and shape features, etc., and can adapt to the complex nonlinear grayscale differences between images. Then, the L-CMOG descriptor is converted to frequency space using fast Fourier transform, and the matching similarity measure is constructed using cross-correlation technology, and the template matching strategy is used to identify the same-name points. The invention method uses the template matching strategy for matching, which is not robust to nonlinear radiation differences. The invention patent application with publication number CN117649541A discloses a multimodal image feature matching method based on window local and global attention. The invention uses the FPN architecture to preliminarily extract a set of features after image data enhancement, uses window attention to interact local features, selects important features within each window, and interacts global information for all features to complete the final feature extraction and achieve feature matching of multimodal images. The invention method only considers the global information of the image, but does not consider the local detail information of the image, and is not robust to nonlinear radiation differences. The invention patent application with publication number CN116012615A discloses a multimodal image feature point detection and matching method. The invention collects multimodal images; extracts feature vectors of multimodal images; calculates the inner product matrix between feature vectors of multimodal images to obtain the similarity of pixel pairs between different modal images; performs Softmax normalization on feature vectors of different modal images to obtain pixel pairs with the highest similarity; obtains similar pixel pairs by screening through a preset similarity threshold to achieve feature point detection and matching. This method uses similarity metrics for pixel matching, which is computationally intensive and the model's generalization ability is insufficient.The invention patent with publication number CN117710711A discloses an optical and SAR image matching method based on a lightweight deep convolutional network. The invention method first uses the Fast algorithm for feature detection, then uses a single deep inverse residual full convolutional neural network to extract the depth descriptor of the image, and finally uses a strategy of refining the coarse to remove the fine to eliminate false matches. The invention method uses the Fast feature detection algorithm to extract feature points. When facing multimodal remote sensing images, the feature points extracted by this algorithm have more noise points. At the same time, a single convolutional neural network is used to obtain the depth descriptor of the image. The extracted feature information has certain limitations, which is not conducive to multimodal image matching tasks. The invention method is mainly used for matching optical and SAR images, and the matching effect is poor for other types of remote sensing images, such as optical and optical, optical and near-infrared, etc. Summary of the invention

[0005] The purpose of the present invention is to provide a multimodal remote sensing image matching method of coupled phase structure and deep features. The method uses a multi-level feature interaction and aggregation network to extract local and global semantic features of the image. The multi-level feature interaction and aggregation module based on the coordinate attention mechanism is designed to realize the interaction and aggregation of features at different stages. The dense feature fusion module based on the coordinate attention mechanism is designed to highlight the expression of key features, and the key features are fused into the deep feature descriptor for image matching. The multi-level mismatch elimination strategy is used to improve the number and accuracy of matching points.

[0006] To achieve the above purpose, the technical solution adopted by the present invention is as follows.

[0007] A multimodal remote sensing image matching method for coupling phase structure and depth features comprises the following steps:

[0008] S1: Detect feature points of multimodal remote sensing images using phase structure feature detection algorithm;

[0009] The phase structure feature detection algorithm includes a point shape texture filtering algorithm, a phase consistency calculation and a Brisk algorithm; the point shape texture filtering algorithm is used to obtain the structural features of the multimodal remote sensing image, the phase consistency calculation is used to enhance the edge information of the multimodal remote sensing image, and the Brisk algorithm is used to perform key point detection to obtain stable key feature points;

[0010] S2: construct a multi-level feature interaction and aggregation network, wherein the multi-level feature interaction and aggregation network includes a dual-branch backbone network, a multi-level feature interaction and aggregation module based on a coordinate attention mechanism, and a dense feature fusion module based on a coordinate attention mechanism;

[0011] The dual-branch backbone network includes an upper branch network and a lower branch network, the upper branch network is a convolutional neural network, which is used to obtain local features of the image, and the lower branch network is a Transformer network, which is used to extract global features of the image;

[0012] The multi-level feature interaction and aggregation module based on the coordinate attention mechanism includes a channel feature interaction module and a spatial feature interaction module. The channel feature interaction module is used to interact and aggregate the features of different stages of the dual-branch backbone network in the channel dimension. The spatial feature interaction module is used to interact and aggregate the features of different stages of the dual-branch backbone network in the spatial dimension.

[0013] The dense feature fusion module based on the coordinate attention mechanism includes a feature aggregation module and a feature fusion module. The feature aggregation module is used to aggregate the local and global features extracted by the dual-branch backbone network. The feature fusion module is used to fuse the aggregated features to obtain a deep feature descriptor and filter redundant features.

[0014] S3: Construct a multimodal remote sensing image dataset and train a multi-level feature interaction and aggregation network. Use the hybrid similarity metric triplet loss function to calculate the minimum distance between the matching deep feature descriptor and the non-matching deep feature descriptor from the deep feature descriptors extracted by the multi-level feature interaction and aggregation network. Calculate the value of the hybrid similarity metric triplet loss function based on the distance and perform back propagation to optimize the network model, prevent gradient disappearance and overfitting, optimize the network model, and complete network training;

[0015] S4: using the trained multi-level feature interaction and aggregation network to extract deep feature descriptors according to the key feature points detected in step S1, and using the nearest neighbor matching algorithm to match the deep feature descriptors to obtain an initial matching result;

[0016] S5: Utilize a multi-level mismatch elimination strategy to eliminate mismatches from the initial matching result to obtain a final matching result; the multi-level mismatch elimination strategy includes adaptive threshold constraints, rough screening based on a grid motion statistics algorithm, and fine matching based on a marginalized sample consensus algorithm.

[0017] Furthermore, the specific method of step S1 is:

[0018] S101: Image preprocessing of multimodal remote sensing images, including geometric correction, atmospheric correction and radiometric calibration;

[0019] S102: calculating an edge confidence map of the preprocessed multimodal remote sensing image, and performing point shape texture filtering on the preprocessed multimodal remote sensing image according to the edge confidence map to obtain structural features of the multimodal remote sensing image;

[0020] S103: performing phase consistency calculation on the filtered multimodal remote sensing image to obtain a maximum moment diagram of the multimodal remote sensing image and highlight edge information;

[0021] S104: Use the Brisk algorithm to perform feature point detection on the maximum moment graph, obtain key feature point information, and obtain a stable feature graph;

[0022] S105: With the key feature points detected in step S104 as the center, a number of image blocks of a size of 224×224 are cropped from the original multimodal remote sensing image to complete image block extraction.

[0023] Furthermore, the convolutional neural network includes five stages, the first stage is a 3×3 convolutional layer, the second stage, the third stage, the fourth stage and the fifth stage are formed by stacking multiple inverted residual structures, and the second stage, the third stage, the fourth stage and the fifth stage stack the inverted residual structure 3 times, 3 times, 7 times and 4 times respectively; wherein, the second stage is formed by stacking an inverted residual structure with a step size of 1, an inverted residual structure with a step size of 2 and an inverted residual structure with a step size of 1, the third stage is formed by stacking an inverted residual structure with a step size of 2 and two inverted residual structures with a step size of 1 in sequence, the fourth stage is formed by stacking an inverted residual structure with a step size of 2 and six inverted residual structures with a step size of 1 in sequence, and the fifth stage is formed by stacking an inverted residual structure with a step size of 2 and three inverted residual structures with a step size of 1;

[0024] The Transformer network comprises five stages, the first stage is a 3×3 convolutional layer, the second stage is four inverted residual structures, and the third, fourth and fifth stages are each composed of an inverted residual structure and a MobileViT structure; wherein, the second stage is composed of three inverted residual structures with a step size of 1 and one inverted residual structure with a step size of 2 stacked together, and the step sizes of the inverted residual structures in the third, fourth and fifth stages are all 2.

[0025] Furthermore, the MobileViT structure includes three stages. The first stage is a 3×3 convolution layer and a 1×1 convolution layer. The second stage first flattens the feature map output by the first stage, and then extracts features through the Transformer module. The Transformer module is divided into two layers. The first layer normalizes the flattened features, and then performs self-attention mechanism calculation on the normalized features. The initial input features of the first layer are added to the features output by the self-attention mechanism to obtain the output of the first layer. The second layer normalizes the features output by the first layer, and then performs multi-layer perceptron operation on the normalized features. Finally, the output features of the first layer are added to the output features of the multi-layer perceptron operation to obtain the output feature map of the final Transformer module. Then, the output feature map of the Transformer module is folded to the size of the input feature map of the second stage. The third stage first uses a 1×1 convolution layer to restore the feature map output by the second stage to the size of the input feature map of the first stage, and then splices the feature map input by the first stage with the feature map output by the 1×1 convolution on the channel, and finally uses a 3×3 convolution layer to fuse to obtain the output feature map.

[0026] Furthermore, the channel feature interaction module includes a coordinate attention mechanism, a global average pooling operation, a global maximum pooling operation, a channel splicing operation, a multi-layer perceptron operation and a Sigmoid activation function, and the spatial feature interaction module includes a channel splicing operation, a multi-layer perceptron operation and a Sigmoid activation function; the feature map obtained by the upper branch network of the dual-branch backbone network is , the feature map obtained by the lower branch network of the dual-branch backbone network is The channel feature interaction module first uses the coordinate attention mechanism to input feature maps and feature map Extract key features, and then perform global maximum pooling and global average pooling operations on the feature map along the channel to obtain four feature maps , , and ; Then, the four feature maps are concatenated to obtain the feature map , , Represents a 3D feature map; finally, the feature map is operated using a multi-layer perceptron Perform deep interaction, then use the Sigmoid activation function to calculate the weight of the feature map after deep interaction and divide it into two weight feature maps of the same size and , , , and then according to the weight feature map And weight feature map Get the feature map after channel correction and ; The calculation formula of the feature map extracted by the channel feature interaction module is expressed as:

[0027] ,

[0028] ,

[0029] ,

[0030] In the formula, Represents the function of the coordinate attention mechanism to calculate the feature map, Indicates a global pooling operation on the feature map. Indicates the global maximum pooling operation on the feature map. Indicates the channel splicing operation on the feature map. represents the multi-layer perceptron operation, Represents the Sigmoid activation function of the channel feature interaction module, represents the channel-by-channel splitting function, Represents the feature map of the upper branch network after channel correction, Represents the feature map of the lower branch network after channel correction, Indicates channel multiplication;

[0031] The spatial feature interaction module first performs the input feature map and feature map Perform channel splicing operation, and then use multi-layer perceptron operation to realize feature interaction and obtain feature map , , and then use the Sigmoid activation function to calculate the feature map The weights of , and divide it into two weight feature maps of the same size and , , , and then according to the weight feature map And weight feature map Get the feature map after spatial correction and ; The calculation formula for extracting feature maps by the spatial feature interaction module is:

[0032] ,

[0033] ,

[0034] ,

[0035] In the formula, The Sigmoid activation function representing the spatial feature interaction module, Represents the feature map of the upper branch network after spatial correction, Represents the feature map of the lower branch network after spatial correction;

[0036] After the feature maps obtained by different branches of the dual-branch backbone network are interacted and corrected by the channel feature interaction module and the spatial feature interaction module, the formula for obtaining the corrected feature map is expressed as:

[0037] ,

[0038] In the formula, It represents the feature map obtained by the upper branch network after the multi-level feature interaction and aggregation module based on the coordinate attention mechanism outputs the feature map. It represents the feature map obtained by the lower branch network after the multi-level feature interaction and aggregation module based on the coordinate attention mechanism outputs the feature map.

[0039] Furthermore, the coordinate attention mechanism includes a parallel horizontal average pooling layer and a vertical average pooling layer, the outputs of which are subjected to a channel splicing operation and then passed through a 1×1 convolution layer, a normalization layer and a nonlinear layer, the nonlinear layer is followed by two parallel 1×1 convolution layers and two parallel Sigmoid activation functions, and the outputs of the two Sigmoid activation functions are respectively element-wise added to the input of the coordinate attention mechanism as the output of the coordinate attention mechanism.

[0040] Furthermore, the dense feature fusion module based on the coordinate attention mechanism includes a feature aggregation module and a feature fusion module, the feature aggregation module includes four 1×1 convolution layers, a 3×3 DW convolution layer, an activation function Relu, a batch normalization layer BN, a Softmax function and a Split function, and the feature fusion module includes a coordinate attention mechanism, a 7×7 convolution layer, a batch normalization layer BN and Normalization;

[0041] The feature aggregation module first performs multi-level feature interaction based on the coordinate attention mechanism and the feature map output by the aggregation module. and feature map Perform channel splicing operations, and use four 1×1 convolutional layers, a 3×3 DW convolutional layer, an activation function Relu, a batch normalization layer BN and a Softmax function to achieve deep aggregation of deep features; secondly, use the Softmax function and Split function to obtain two weight maps and , and then use the weight map and weight graph Get the aggregated dense feature map and ; The calculation formula of the feature aggregation module aggregation feature map is expressed as:

[0042] ,

[0043] ,

[0044] ,

[0045] In the formula, They represent the feature maps output by two 1×1 convolutional layers, represents the Softmax function, Represents a composite function consisting of a 1×1 convolution layer, a 3×3 DW convolution layer, a batch normalization layer BN, and an activation function Relu. represents the dense feature map of the upper branch network output by the feature aggregation module, Represents the dense feature map of the branch network output under the feature aggregation module output;

[0046] The feature fusion module is used to fuse the aggregated dense features to obtain a deep feature descriptor for matching. First, the dense feature map output by the feature aggregation module is and dense feature maps Perform channel splicing operations and use the coordinate attention mechanism to enhance the expression of key features and suppress irrelevant features; secondly, use a 7×7 convolution layer and a batch normalization layer BN to obtain a 256-dimensional feature vector , and the eigenvector conduct Normalization operation obtains deep feature descriptor ; The expression of the feature fusion module is:

[0047] ,

[0048] ,

[0049] In the formula, Represents a composite function consisting of a 7×7 convolutional layer and a batch normalization layer BN. represents the deep feature descriptor, express Normalization operation.

[0050] Further, in step S3, the hybrid similarity metric triplet loss function is given by norm normalization and triplet loss, The calculation formula for norm normalization is:

[0051] ,

[0052] In the formula, Descriptor pair and conduct The result after norm normalization is and yes A pair of matched deep feature descriptors before norm normalization, Indicates the number of deep feature descriptors;

[0053] The triplet loss is calculated as:

[0054] ,

[0055] ,

[0056] In the formula, represents the triplet loss value, represents the angle between the matching deep feature descriptor pairs, represents the angle between non-matching deep feature descriptor pairs, represents a hybrid similarity measure for matching pairs of deep feature descriptors, Representing a hybrid similarity measure for pairs of non-matching deep feature descriptors; represents a mixed similarity measure, represents the cosine value of the angle between the deep feature descriptors, , A function representing the angle of the deep feature descriptor, , is a scalar used to adjust the ratio, ; is the maximum magnitude gradient normalization factor, is the angle between the deep feature descriptors, Represents the margin, which is a constant with a value of 1;

[0057] Finally, the hybrid similarity metric triplet loss function is defined as:

[0058] ,

[0059] In the formula, represents the hybrid similarity metric triplet loss function, is the regularization parameter.

[0060] Furthermore, the initial matching result obtained by the nearest neighbor matching algorithm contains the nearest matching point pair and the second nearest matching point pair in Euclidean distance. If the distance between the nearest matching point pair and the depth feature descriptor is Much smaller than the distance between the next nearest matching point and the deep feature descriptor , then the corresponding initial matching point pair is considered to have good quality. The adaptive threshold constraint calculates the mean of the Euclidean distances between the nearest matching point pair and the second nearest matching point pair in all initial matching results as the basis for judging whether to retain the corresponding matching point pair. The calculation formula for the mean of the Euclidean distance is:

[0061] ;

[0062] in, is the number of key feature points on the reference image, is the mean of the Euclidean distance, is the distance between the next neighboring matching point and the depth feature descriptor, is the distance between the nearest matching point pair and the depth feature descriptor. For each set of nearest matching point pairs and second nearest matching point pairs, when When , the corresponding nearest matching point pair is retained, otherwise the matching point pair is removed to filter out the initial correct matching point pair result;

[0063] Then, the initial correct matching point pairs obtained after adaptive threshold constraint screening are subjected to motion smoothing constraint statistics based on the grid motion statistics algorithm to obtain the matching point pairs after rough screening; then, the marginalized sample consensus algorithm is used to solve the estimated value using iterative weighted least squares to achieve robust estimation, realize precise matching, and obtain the final correct matching points.

[0064] Furthermore, the rough screening based on the grid motion statistics algorithm is to find the matching points that meet the conditions from the initial correct matching point pairs obtained after the adaptive threshold constraint screening, find the nearest neighborhood and the close neighborhood with the coordinates of the input initial matching point as the center, count the number of matches in the neighborhood, and then count the probabilities of correct matches and incorrect matches to obtain the probabilities of correct matches and incorrect matches in each neighborhood, determine the correct and incorrect matches according to the probability threshold, and the more matching points there are in the image, the higher the distinction between correct and incorrect matches; in the process of implementing the grid motion statistics algorithm, the grid is used instead of the neighborhood to realize the fast and efficient operation of the algorithm.

[0065] Furthermore, the marginalized sample consensus algorithm is an iterative algorithm for correctly estimating mathematical model parameters from a set of data containing erroneous matching points. The probability of a random point being a correct matching point is solved by marginalization, and the probability is used as a weight for weighted least squares optimization, thereby optimizing the model. The model with the largest number of correct matching points is found through iteration, and the final matching points are obtained.

[0066] The beneficial technical effects of the present invention are:

[0067] (1) The present invention constructs a multimodal remote sensing image matching method that couples phase structure and depth features. The method proposes a phase structure feature detection algorithm, uses a point shape texture filtering algorithm to obtain the structural features of the image, then calculates the phase consistency of the filtered image, and finally uses the Brisk algorithm to extract stable key feature points. The phase structure feature detection algorithm can effectively resist the influence of significant nonlinear radiation differences and noise on the detection of key feature points of the image. The method uses a dual-branch backbone network to realize the extraction of local and global features of the image. The multi-level feature interaction and aggregation module based on the coordinate attention mechanism is designed to interact and aggregate features at different stages to promote the expression of key features. The dense feature fusion module based on the coordinate attention mechanism focuses on the expression of key features, fuses the features extracted by the dual-branch backbone network into the deep feature descriptor, filters redundant features, and solves the problem that the existing network model structure is single and difficult to extract local and global features of the image at the same time, and limits the interaction and aggregation of features; the network is trained using a self-made multimodal remote sensing image dataset to solve the problem that the lack of comprehensive training data limits the model's feature expression and performance.

[0068] (2) In order to construct a deep feature descriptor that is robust to significant nonlinear radiation differences and complex geometric distortions, a dense feature fusion module based on the coordinate attention mechanism is used to fuse dense features into the deep feature descriptor, effectively fusing the local features and global features extracted by the network model; a multi-level false matching elimination strategy is adopted in the false matching elimination stage to improve the number and accuracy of correct matching points; the problem that multimodal remote sensing images are difficult to match due to large nonlinear radiation differences and complex geometric distortions is solved; the multimodal remote sensing image matching method of the coupled phase structure and deep features of the present invention achieves good results in terms of matching quantity, accuracy and robustness, can meet the image matching requirements between multimodal remote sensing images due to complex geometric distortions and large nonlinear radiation differences, and has good applicability to multimodal remote sensing image matching. BRIEF DESCRIPTION OF THE DRAWINGS

[0069] Figure 1 is a flow chart of a multimodal remote sensing image matching method for coupling phase structure and depth features in an embodiment of the present invention;

[0070] Figure 2 It is a flow chart of a phase structure feature detection method of a multimodal remote sensing image matching method coupling phase structure and depth features in an embodiment of the present invention;

[0071] Figure 3Schematic diagram of the network model structure of the multimodal remote sensing image matching method of coupling phase structure and depth features in an embodiment of the present invention;

[0072] Figure 4 It is a structural schematic diagram of a multi-level feature interaction and aggregation module based on a coordinate attention mechanism of a multimodal remote sensing image matching method for coupling phase structure and depth features in an embodiment of the present invention;

[0073] Figure 5 Schematic diagram of the structure of a dense feature fusion module based on a coordinate attention mechanism in a multimodal remote sensing image matching method coupling phase structure and depth features in an embodiment of the present invention;

[0074] Figure 6 Schematic diagram of the inverse residual structure of the multimodal remote sensing image matching method coupling phase structure and depth features in an embodiment of the present invention;

[0075] Figure 7 Schematic diagram of the MobileViT structure of the multimodal remote sensing image matching method coupling phase structure and depth features in an embodiment of the present invention;

[0076] Figure 8 is a schematic diagram of three groups of test image pairs in an embodiment of the present invention;

[0077] Fig. 9 Schematic diagram of the matching result of an optical image pair in an embodiment of the present invention;

[0078] Fig.10 Schematic diagram of matching results of optical and infrared image pairs in an embodiment of the present invention;

[0079] Fig.11 Schematic diagram of the matching results of the optical and SAR image pairs in an embodiment of the present invention. DETAILED DESCRIPTION

[0080] In order to more clearly understand the above-mentioned purpose, features and advantages of the present invention, the present invention is further described in detail below in conjunction with the accompanying drawings and specific embodiments. In the following description, a lot of specific details are set forth to facilitate a full understanding of the present invention, but the present invention can also be implemented in other ways different from those described herein, and therefore, the present invention is not limited to the limitation of the specific embodiments disclosed below. Unless otherwise defined, the technical terms or scientific terms used herein should be the usual meanings understood by people with ordinary skills in the field described in the application.

[0081] like Figure 1 As shown, a multimodal remote sensing image matching method for coupling phase structure and depth features includes the following steps:

[0082] S1: Use phase structure feature detection algorithm to detect feature points of multimodal remote sensing images, such as Figure 2 As shown; the specific method is:

[0083] S101: Image preprocessing of multimodal remote sensing images, including geometric correction, atmospheric correction and radiation calibration;

[0084] S102: calculating an edge confidence map of the preprocessed multimodal remote sensing image, and performing point shape texture filtering on the preprocessed multimodal remote sensing image according to the edge confidence map to obtain structural features of the multimodal remote sensing image;

[0085] S103: performing phase consistency calculation on the filtered multimodal remote sensing image to obtain a maximum moment diagram of the multimodal remote sensing image and highlight edge information;

[0086] S104: Use the Brisk algorithm to perform feature point detection on the maximum moment graph, obtain key feature point information, and obtain a stable feature graph;

[0087] S105: With the key feature points detected in step S104 as the center, a number of image blocks of a size of 224×224 are cropped from the original multimodal remote sensing image to complete image block extraction.

[0088] S2: Construct a multi-level feature interaction and aggregation network, such as Figure 3 As shown, the multi-level feature interaction and aggregation network includes a dual-branch backbone network, a multi-level feature interaction and aggregation module based on a coordinate attention mechanism, and a dense feature fusion module based on a coordinate attention mechanism. The dual-branch backbone network includes an upper branch network and a lower branch network. The upper branch network is a convolutional neural network (CNN), which is mainly used to extract local features of the image, and the lower branch network is a Transformer network, which is mainly used to extract global features of the image. The multi-level feature interaction and aggregation module based on the coordinate attention mechanism is used to interact and aggregate the local and global features of the dual-branch backbone network at different stages, and the dense feature fusion module based on the coordinate attention mechanism is used to fuse the local and global features extracted by the dual-branch backbone network to obtain a deep feature descriptor. The image block with R, G, and B channel information is subjected to a multi-level feature interaction and aggregation network to obtain a 256-dimensional deep feature descriptor.

[0089] The upper branch network of the dual-branch backbone network, i.e., the convolutional neural network, includes five stages. The first stage is a 3×3 convolutional layer. The second, third, fourth and fifth stages are composed of multiple inverted residual structures stacked together. The second, third, fourth and fifth stages stack the inverted residual structures 3 times, 3 times, 7 times and 4 times respectively. Among them, the second stage is composed of an inverted residual structure with a step size of 1, an inverted residual structure with a step size of 2 and an inverted residual structure with a step size of 1 stacked together. The third stage is composed of an inverted residual structure with a step size of 2 and two inverted residual structures with a step size of 1 stacked together in sequence. The fourth stage is composed of an inverted residual structure with a step size of 2 and six inverted residual structures with a step size of 1 stacked together in sequence. The fifth stage is composed of an inverted residual structure with a step size of 2 and three inverted residual structures with a step size of 1 stacked together. The sizes of the feature maps output by each stage of the upper branch network are respectively , , , , ,in, and are the length and width of the image input to the network respectively.

[0090] The lower branch network of the dual-branch backbone network, namely the Transformer network, contains five stages. The first stage is a 3×3 convolutional layer, the second stage is four inverted residual structures, and the third, fourth and fifth stages are composed of an inverted residual structure and a MobileViT structure. Among them, the second stage is composed of three inverted residual structures with a step size of 1 and an inverted residual structure with a step size of 2. The step sizes of the inverted residual structures in the third, fourth and fifth stages are all 2. The size of the feature map output by each stage of the lower branch network is respectively , , , , ,in, and are the length and width of the image input to the network respectively.

[0091] like Figure 6 As shown, the inverted residual structure is divided into two structures according to the different step sizes. Figure 6 (a) in the figure is an inverse residual structure with a step size of 1. Figure 6(b) in the figure is an inverted residual structure with a stride of 2; the inverted residual structure with a stride of 1 is divided into four stages. The first stage uses a 1×1 convolution layer for dimensionality increase, followed by a batch normalization layer BN and an activation function Relu6. The second stage is a 3×3 DW convolution layer, a batch normalization layer BN and an activation function Relu6. The third stage is a 1×1 convolution layer, a linear layer and a batch normalization layer BN. The fourth stage concatenates the input features and the features output by the third stage on the feature channel. The inverted residual structure with a stride of 2 It is divided into three stages. The first stage is a 1×1 convolution layer, a batch normalization layer BN and an activation function Relu6. The second stage is a 3×3 DW convolution layer, a batch normalization layer BN and an activation function Relu6. The third stage is a 1×1 convolution layer, a linear layer and a batch normalization layer BN. The DW convolution layer is a convolution kernel responsible for one channel, and one channel is only convolved by one convolution kernel. The DW convolution layer can greatly reduce the model parameters of the network and make the model lightweight. The calculation formula for extracting the feature map of the inverted residual structure is as follows:

[0092] ,

[0093] or,

[0094] ;

[0095] in, For the The output of the inverted residual structure, For the The output of the inverted residual structure, It is a composite function composed of convolutional layer, activation function Relu6 and batch normalization layer BN. For the The inverted residual structure outputs a 3D tensor of feature maps, where are the dimensions of the 3D tensor corresponding to the feature map, and correspond to the length, width, and number of channels of the 3D tensor of the feature map, respectively. For the The inverted residual structure outputs a 3D tensor of feature maps, where are the dimensions of the 3D tensor corresponding to the feature map, and correspond to the length, width, and number of channels of the 3D tensor of the feature map, respectively. is the step size of the DW convolution in the inverted residual structure.

[0096] The input size of the inverted residual structure with a step size of 1 is The size of the feature map output by the first stage 1×1 convolutional layer is ,in is the expansion factor, is the number of channels of the feature map. In the embodiment of the present invention, the expansion factor is set to 0.6. Then, the size of the feature map output by the second stage DW convolution layer is , and then the size of the feature map output by the third stage 1×1 convolutional layer and linear layer is , the size of the feature map output by the fourth stage is , represents the number of channels of the feature map finally output by the inverted residual structure with a step size of 1; the sizes of the feature maps output by each stage of the inverted residual structure with a step size of 2 are , and , Represents the number of channels of the feature map finally output by the inverted residual structure with a stride of 2.

[0097] like Figure 7 As shown in the figure, the MobileViT structure is mainly divided into three stages. The first stage is a 3×3 convolution layer and a 1×1 convolution layer. The second stage first flattens the feature map output by the first stage, and then extracts features through the Transformer module. The Transformer module is mainly divided into two layers. The first layer normalizes the flattened features, and then performs a self-attention mechanism (Self-Attention) calculation on the normalized features. The initial input features of the first layer are added to the features output by the self-attention mechanism to obtain the output of the first layer. The second layer normalizes the features output by the first layer, and then performs a multi-layer perceptron operation (i.e., MLP) on the normalized features. Finally, the output features of the first layer are added to the output features of the multi-layer perceptron operation to obtain the output feature map of the final Transformer module. The output feature map of the Transformer module is then folded to the size of the input feature map of the second stage. In the third stage, a 1×1 convolution layer is first used to restore the feature map output by the second stage to the size of the input feature map of the first stage. The feature map input by the first stage and the feature map output by the 1×1 convolution are then concatenated on the channel. Finally, a 3×3 convolution layer is used to fuse them to obtain the output feature map.

[0098] like Figure 4 As shown, the multi-level feature interaction and aggregation module based on the coordinate attention mechanism includes a channel feature interaction module and a spatial feature interaction module. The channel feature interaction module includes a coordinate attention mechanism (i.e. CA), a global average pooling operation (i.e. AP), a global maximum pooling operation (i.e. MP), a channel concatenation operation (i.e. Concatenate operation), a multi-layer perceptron operation and a Sigmoid activation function. The spatial feature interaction module includes a channel concatenation operation, a multi-layer perceptron operation and a Sigmoid activation function. The feature map obtained by the upper branch network of the dual-branch backbone network is , the feature map obtained by the lower branch network of the dual-branch backbone network is The channel feature interaction module first uses the coordinate attention mechanism to input feature maps and feature map Extract key features, and then perform global maximum pooling and global average pooling operations on the feature map along the channel to obtain four feature maps , , and Then, the four feature maps are concatenated to obtain the feature map , , represents a 3D feature map, 1×1×4C represents a feature map Finally, the multi-layer perceptron is used to operate on the feature map Perform deep interaction, then use the Sigmoid activation function to calculate the weight of the feature map after deep interaction and divide it into two weight feature maps of the same size and , , , and then according to the weight feature map And weight feature map Get the feature map after channel correction and ; The calculation formula of the feature map extracted by the channel feature interaction module is expressed as:

[0099] ,

[0100] ,

[0101] ,

[0102] In the formula, Represents the function of the coordinate attention mechanism to calculate the feature map, Indicates a global pooling operation on the feature map. Indicates the global maximum pooling operation on the feature map. Indicates the channel splicing operation on the feature map. represents the multi-layer perceptron operation, Represents the Sigmoid activation function of the channel feature interaction module, represents the channel-by-channel splitting function, Represents the feature map of the upper branch network after channel correction, Represents the feature map of the lower branch network after channel correction, Represents channel multiplication.

[0103] The spatial feature interaction module first performs the input feature map and feature map Perform channel splicing operation, and then use multi-layer perceptron operation to realize feature interaction and obtain feature map , , and then use the Sigmoid activation function to calculate the feature map The weights of , and divide it into two weight feature maps of the same size and , , , and then according to the weight feature map And weight feature map Get the feature map after spatial correction and The calculation formula for extracting feature maps from the spatial feature interaction module is:

[0104] ,

[0105] ,

[0106] ,

[0107] In the formula, Indicates the channel splicing operation on the feature map, MLP indicates the multi-layer perceptron operation, Represents the feature map obtained after the multi-layer perceptron operation. The Sigmoid activation function representing the spatial feature interaction module, represents the channel-by-channel splitting function, Represents the feature map of the upper branch network after spatial correction, Represents the feature map of the lower branch network after spatial correction, Represents channel multiplication. H×W represents the length and width of the feature map.

[0108] After the feature maps obtained by different branches of the dual-branch backbone network are interacted and corrected by the channel feature interaction module and the spatial feature interaction module, the formula for obtaining the corrected feature map is expressed as:

[0109] ,

[0110] In the formula, It represents the feature map obtained by the upper branch network after the multi-level feature interaction and aggregation module based on the coordinate attention mechanism outputs the feature map. It represents the feature map obtained by the lower branch network after the multi-level feature interaction and aggregation module based on the coordinate attention mechanism outputs the feature map.

[0111] The coordinate attention mechanism includes a parallel horizontal average pooling layer and a vertical average pooling layer, followed by a channel splicing operation, a 1×1 convolution layer, a normalization layer, a nonlinear layer, and then two parallel 1×1 convolution layers and two parallel Sigmoid activation functions. Among them, the input feature map first passes through the parallel horizontal average pooling layer and the vertical average pooling layer, and the two output feature maps obtained are subjected to channel splicing operations and then pass through a 1×1 convolution layer, a normalization layer and a nonlinear layer, and then pass through two parallel 1×1 convolution layers. The output feature map passes through the parallel Sigmoid activation function and then is element-wise added to the input feature map of the coordinate attention mechanism to obtain two output feature maps of the coordinate attention mechanism.

[0112] like Figure 5 As shown in the figure, the dense feature fusion module based on the coordinate attention mechanism includes a feature aggregation module and a feature fusion module. The feature aggregation module includes four 1×1 convolution layers, a 3×3 DW convolution layer (i.e., a depth-separable convolution layer), an activation function Relu, a batch normalization layer BN, a Softmax function, and a Split function. The feature fusion module includes a coordinate attention mechanism, a 7×7 convolution layer, a batch normalization layer BN, and Normalization.

[0113] The feature aggregation module first performs multi-level feature interaction based on the coordinate attention mechanism and the feature map output by the aggregation module. and feature map Perform channel splicing operations and use four 1×1 convolutional layers, a 3×3 DW convolutional layer, an activation function Relu, a batch normalization layer BN and a Softmax function to achieve deep aggregation of deep features. Secondly, use the Softmax function and Split function to obtain two weight maps and , and then use the weight map and weight graph Get the aggregated dense feature map and The calculation formula of the feature aggregation module aggregation feature map can be expressed as:

[0114] ,

[0115] ,

[0116] ,

[0117] In the formula, They represent the feature maps output by two 1×1 convolutional layers, represents the Softmax function, Represents a composite function consisting of a 1×1 convolution layer, a 3×3 DW convolution layer, a batch normalization layer BN, and an activation function Relu. Indicates the channel splicing operation on the feature map. represents the dense feature map of the upper branch network output by the feature aggregation module, Represents the dense feature map of the branch network output under the feature aggregation module output.

[0118] The feature fusion module is used to fuse the aggregated dense features to obtain a deep feature descriptor for matching. First, the dense feature map output by the feature aggregation module is and dense feature maps Perform channel splicing operations and use the coordinate attention mechanism to enhance the expression of key features and suppress irrelevant features. Secondly, use a 7×7 convolution layer and a batch normalization layer BN to obtain a 256-dimensional feature vector , and the eigenvector conduct Normalization operation (i.e. L2 Norm) obtains the deep feature descriptor The expression of the feature fusion module is:

[0119] ,

[0120] ,

[0121] In the formula, Represents a composite function consisting of a 7×7 convolutional layer and a batch normalization layer BN. represents the deep feature descriptor, express Normalization operation.

[0122] S3: Use the constructed multimodal remote sensing image dataset (optical-optical, optical-infrared, optical-SAR) to train the multi-level feature interaction and aggregation network, and use the loss function to calculate the minimum distance between the matching deep feature descriptor and the non-matching deep feature descriptor from the deep feature descriptors extracted by the multi-level feature interaction and aggregation network. Calculate the loss function value based on the distance and perform back propagation to optimize the network model, prevent gradient vanishing and overfitting, optimize the network model, and complete network training.

[0123] The multimodal remote sensing image dataset constructed in this embodiment includes three types of data: optical and optical image dataset, optical and infrared image dataset, and optical and SAR image dataset. The multimodal remote sensing image dataset has 60,000 pairs of image blocks, with 20,000 pairs of each type. (1) Optical-optical: The optical-optical image dataset consists of self-cropped Aowei interactive map satellite images and WUH datasets. 40 pairs of multi-temporal satellite images with a size of 5000×5000 pixels were obtained from the Aowei interactive map satellite images. The images were taken from 2020 to 2024, covering areas including cities, rural areas, farmlands, forests, rivers, coastlines, etc., with a resolution of 5-10m. After the images were geo-referenced and geometrically corrected, the large-size images were cropped into 14,000 pairs of small images of 256×256 pixels. The WUH dataset contains a set of multi-temporal aerial photographs, consisting of two pre-registered photos. The aerial photos were taken in the same area of ​​New Zealand in 2012 and 2016, respectively, with a resolution of 0.075m. This dataset is used for change detection, and the ground objects have changed. We manually selected areas with less ground object changes and cropped them to obtain 6,000 pairs of image blocks. (2) Optical-infrared: The optical-infrared image dataset mainly consists of Landsat-9 satellite images and Sentinel-2 satellite images. The Sentinel-2 satellite image data mainly comes from the SEN12MS-CR dataset. Several Landsat-9 satellite images were obtained from the Geospatial Data Cloud, covering plains, mountains and hilly areas with a resolution of 30m. The optical image is a true color image synthesized by band 2, band 3 and band 4 of the Landsat-9 satellite, taken in 2022, and the infrared image is a band 5 image of the Landsat-9 satellite, taken in 2023. After image preprocessing, the optical and infrared images were cropped into 10,000 pairs of small images with a size of 256×256 pixels. The SEN12MS-CR dataset is a large cloud-free dataset. We selected the cloud-free data in summer from this dataset. The dataset contains the full-band data of the Sentinel-2 satellite. We extracted the optical images synthesized from bands 2, 3, and 4 of the Sentinel-2 satellite. Band 8 of the Sentinel-2 satellite is an infrared image. We manually selected 10,000 pairs of images with good quality from the infrared images of band 8, with a size of 256×256. (3) Optical-SAR: The images of the optical-SAR image dataset come from the QXS-SAROPT dataset and the SAR2Opt dataset. The QXS-SAROPT dataset consists of 20,000 pairs of optical and SAR remote sensing image patches obtained from the Gaofen-3 SAR satellite images and the Aowei interactive map satellite images. The images cover three port cities: Santiago, Shanghai, and Qingdao.We manually selected 16,000 pairs of images with good image quality from the QXS-SAROPT dataset to enter the optical-SAR dataset, and the image size is 256×256. The SAR2Opt dataset mainly consists of TerraSAR-X satellite images and Aowei interactive map satellite images. TerraSAR-X satellite images collected SAR images with a spatial resolution of 1m for 10 cities in Asia, North America, Oceania and Europe from 2007 to 2013. The images obtained from the Aowei interactive map satellite images were aligned with the optical images by manually selecting control points, and the image size is 600×600. We only selected the data of the training set of the SAR2Opt dataset and cropped it into 6,000 pairs of 256×256 images. During training, the image blocks were rotated and translated, and finally randomly cropped into image blocks of size 224×224. The optical-optical image dataset, the optical and infrared image dataset, and the optical and SAR image dataset are combined into a multimodal remote sensing image dataset, and the multimodal remote sensing image dataset is randomly divided into a training set and a test set of the multi-level feature interaction and aggregation network in a ratio of 8:2. The multi-level feature interaction and aggregation network is trained using the image blocks in the training set and the test set, and the distance of the deep feature descriptors is calculated using the loss function, the distance between the matching deep feature descriptors is continuously reduced, and the distance between the non-matching deep feature descriptors is increased. The loss value is calculated by distance for back propagation, the model parameters are optimized, and the model training is completed.

[0124] A hybrid similarity metric triplet loss function is used to train the multi-level feature interaction and aggregation network. The loss function calculates the distance between the matching deep feature descriptors and the non-matching deep feature descriptors. The loss function continuously reduces the distance between the matching deep feature descriptors and increases the distance between the non-matching deep feature descriptors. This makes the deep feature descriptors extracted by the multi-level feature interaction and aggregation network robust to nonlinear radiation differences, thereby improving the matching performance of the deep feature descriptors.

[0125] The triplet loss function for the hybrid similarity metric is given by norm normalization and triplet loss, Norm normalization allows the network to generate better gradients during training and enhances the network's robustness to image intensity; The calculation formula for norm normalization is:

[0126] ,

[0127] In the formula, Descriptor pair and conduct The result after norm normalization is and yes A pair of matched deep feature descriptors before norm normalization, Indicates the number of deep feature descriptors.

[0128] In order to better generate gradients during network training, a hybrid similarity metric is used to improve the triplet loss. The calculation formula for the triplet loss is:

[0129] ,

[0130] ,

[0131] In the formula, represents the triplet loss value, represents the angle between the matching deep feature descriptor pairs, represents the angle between non-matching deep feature descriptor pairs, represents a hybrid similarity measure for matching pairs of deep feature descriptors, Representing a hybrid similarity metric for pairs of non-matching deep feature descriptors. represents a mixed similarity measure, represents the cosine value of the angle between the deep feature descriptors, , A function representing the angle of the deep feature descriptor, , is a scalar used to adjust the ratio, . is the maximum magnitude gradient normalization factor, is the angle between the deep feature descriptors, Represents the margin, which is a constant with a value of 1.

[0132] Finally, the hybrid similarity metric triplet loss function is defined as:

[0133] ,

[0134] In the formula, represents the hybrid similarity metric triplet loss function, is the regularization parameter.

[0135] S4: Use the trained multi-level feature interaction and aggregation network to construct a deep feature descriptor according to the key feature points detected in step S1, input the image block obtained with the key feature point as the center in step S1 into the trained multi-level feature interaction and aggregation network for feature extraction and feature fusion to obtain a deep feature descriptor, normalize the obtained deep feature descriptor, and use the nearest neighbor matching algorithm to match it to obtain an initial matching result corresponding to the deep feature descriptor.

[0136] S5: Utilize a multi-level mismatch elimination strategy to eliminate mismatches from the initial matching results to obtain the final matching results. The multi-level mismatch elimination strategy includes adaptive threshold constraints, rough screening based on a grid motion statistics algorithm, and fine matching based on a marginalized sample consensus algorithm (i.e., MAGSAC++ algorithm).

[0137] Adaptive threshold constraint: The initial matching result obtained by the nearest neighbor matching algorithm contains the nearest matching point pairs and the second nearest matching point pairs in Euclidean distance. The adaptive threshold constraint can filter out the initial correct matching point pairs from the initial matching result. If the distance between the nearest matching point pair and the deep feature descriptor is Much smaller than the distance between the next nearest matching point and the deep feature descriptor , then the corresponding initial matching point pair is considered to have good quality. The adaptive threshold constraint calculates the mean of the Euclidean distances between the nearest matching point pair and the second nearest matching point pair in all initial matching results as the basis for judging whether to retain the corresponding matching point pair. The calculation formula for the mean of the Euclidean distance is:

[0138] ;

[0139] in, is the number of key feature points on the reference image, is the mean of the Euclidean distance, is the distance between the next neighboring matching point and the depth feature descriptor, is the distance between the nearest matching point pair and the depth feature descriptor. For each set of nearest matching point pairs and second nearest matching point pairs, when , the corresponding nearest matching point pair is retained, otherwise the matching point pair is removed.

[0140] Rough screening based on grid motion statistics algorithm and precise matching based on marginalized sample consensus algorithm: After the initial matching results are subjected to adaptive threshold constraints, more mismatching points can be eliminated, but there are still mismatching points. First, the initial correct matching point pairs obtained after adaptive threshold constraint screening are subjected to motion smoothing constraint statistics based on the grid motion statistics algorithm (i.e., GMS algorithm) to obtain the matching point pairs after rough screening. The grid motion statistics algorithm can effectively solve the problem that the matching points are interfered by more feasible solutions. The marginalized sample consensus algorithm uses iterative weighted least squares to solve the estimated value to achieve robust estimation, and finally obtains the correct matching points, achieves precise matching, and improves the number and accuracy of matching points.

[0141] The rough screening based on the grid motion statistics algorithm is to find the matching points that meet the conditions from the initial correct matching point pairs obtained after the adaptive threshold constraint screening, find the nearest neighbor and the close neighbor with the coordinates of the input initial matching point as the center, count the number of matches in the neighborhood, and then count the probabilities of correct matches and wrong matches, and obtain the probabilities of correct matches and wrong matches in each neighborhood, and determine the correct and wrong matches according to the probability threshold. The more matching points there are in the image, the higher the degree of distinction between correct and wrong matches. In the process of implementing the grid motion statistics algorithm, the grid is used instead of the neighborhood to realize the fast and efficient operation of the algorithm.

[0142] The marginalized sample consensus algorithm is an iterative algorithm for correctly estimating the parameters of the mathematical model from a set of data containing incorrect matching points. The probability of a random point being a correct matching point is solved by marginalization, and the probability is used as a weight for weighted least squares optimization to optimize the model. The model with the largest number of correct matching points is found through iteration, and the final matching point can be obtained.

[0143] The embodiment of the present invention uses a trained multi-level feature interaction and aggregation network to obtain the deep feature descriptor of the image to be matched. The deep feature descriptor obtains the image matching result through feature matching and multi-level mismatch elimination strategy. The matching result is evaluated for accuracy. The evaluation index formula is as follows:

[0144] ;

[0145] ;

[0146] in, For difference determination, It indicates the position of the feature points matched on the image to be matched by the algorithm based on the key feature points on the reference image. It is the theoretical position of the key feature points on the image to be matched. Indicates the actual position of the key feature point on the image to be matched; if the theoretical position of the key feature point on the image to be matched The actual position of the key feature points on the image to be matched The pixel difference between the distances is less than a given threshold , then the feature points matched on the image to be matched are considered to be correct matching points. is the number of correct matching points, is the root mean square error, is the original coordinate of the key feature point on the image to be matched, is the coordinate of the matching point on the image to be matched after the change.

[0147] In addition, the matching method described in the embodiment of the present invention is compared with the existing more advanced image matching methods, including the RIFT method, MS-HLMO method and WSSF method of manually designed features, and the R2D2 algorithm and CMM-Net algorithm based on deep learning. 30 pairs of images are selected as test images, and the matching results are shown in Table 1 below.

[0148] Table 1 Statistics of matching results between the matching method according to the embodiment of the present invention and five existing image matching algorithms:

[0149]

[0150] Among them, the matching method described in the embodiment of the present invention achieves the best matching effect among 30 pairs of test images, with the highest mean of correct matching points and the highest root mean square error. All are less than 2 pixels, and the points are evenly distributed; the matching effect of the matching method described in the example of the present invention on the optical-optical image pair is comparable to that of the existing method, and the matching effect on the optical-infrared and optical-SAR image pairs is better than that of the existing method. On the optical-SAR image, the existing methods all show strong sensitivity, and the matching method proposed in the present invention can achieve good matching results in the matching of optical-SAR image pairs. The method described in the example of the present invention has good robustness in the face of complex geometric distortions and significant nonlinear radiation differences between multimodal remote sensing images, and can meet the needs of multimodal remote sensing image matching. The method proposed in the present invention has good applicability to multimodal remote sensing image matching tasks and can provide valuable reference for multimodal matching tasks. Figure 8 These are three groups of image pairs in the test images of the embodiment of the present invention, where: Figure 8 (a) is an optical image. Figure 8 (b) is an optical image. Figure 8 (c) in the figure is an optical image. Figure 8 (d) in the figure is an infrared image. Figure 8 (e) in the figure is an optical image. Figure 8 (f) in is the SAR image, that is, Figure 8 (a) and Figure 8 (b) in the figure constitutes an optical and optical image pair, Figure 8 (c) and Figure 8 (d) in the figure constitutes an optical and infrared image pair, Figure 8 (e) and Figure 8 (f) in which optical and SAR image pairs are composed; Fig. 9 , Fig.10 and Fig.11 for Figure 8 Schematic diagram of image matching results corresponding to the image pair shown in , where: Fig. 9is the matching result diagram of optical and optical images. Fig.10 is the matching result diagram of the optical and infrared image pairs, Fig.11 is the matching result diagram of the optical and SAR image pairs, Fig. 9 , Fig.10 and Fig.11 (a) is the image matching result of the RIFT method. Fig. 9 , Fig.10 and Fig.11 (b) is the image matching result of MS-HLMO method. Fig. 9 , Fig.10 and Fig.11 (c) in the figure is the image matching result of the WSSF method. Fig. 9 , Fig.10 and Fig.11 (d) in the figure is the image matching result of the R2D2 method. Fig. 9 , Fig.10 and Fig.11 (e) in the figure is the image matching result of the CMM-Net method. Fig. 9 , Fig.10 and Fig.11 (f) in FIG. 1 is the image matching result of the method of the present invention. It can be seen that the matching method described in the embodiment of the present invention can accurately and quickly achieve the matching of multi-modal remote sensing images, and performs best in matching effect.

[0151] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. For those skilled in the art, the present invention may have various modifications and variations. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention shall be included in the protection scope of the present invention.

Claims

1. A multimodal remote sensing image matching method coupling phase structure and depth features, characterized in that: The steps include: S1: Detect feature points of a multimodal remote sensing image using a phase structure feature detection algorithm, wherein the phase structure feature detection algorithm includes a point shape texture filtering algorithm, a phase consistency calculation and a Brisk algorithm, wherein the point shape texture filtering algorithm is used to obtain structural features of the multimodal remote sensing image, the phase consistency calculation is used to enhance edge information of the multimodal remote sensing image, and the Brisk algorithm is used to perform key point detection to obtain stable key feature points; S2: construct a multi-level feature interaction and aggregation network, wherein the multi-level feature interaction and aggregation network includes a dual-branch backbone network, a multi-level feature interaction and aggregation module based on a coordinate attention mechanism, and a dense feature fusion module based on a coordinate attention mechanism; The dual-branch backbone network includes an upper branch network and a lower branch network, the upper branch network is a convolutional neural network, which is used to obtain local features of the image, and the lower branch network is a Transformer network, which is used to extract global features of the image; The multi-level feature interaction and aggregation module based on the coordinate attention mechanism includes a channel feature interaction module and a spatial feature interaction module. The channel feature interaction module is used to interact and aggregate the features of different stages of the dual-branch backbone network in the channel dimension. The spatial feature interaction module is used to interact and aggregate the features of different stages of the dual-branch backbone network in the spatial dimension. The dense feature fusion module based on the coordinate attention mechanism includes a feature aggregation module and a feature fusion module. The feature aggregation module is used to aggregate the local and global features extracted by the dual-branch backbone network. The feature fusion module is used to fuse the aggregated features to obtain a deep feature descriptor and filter redundant features. S3: Construct a multimodal remote sensing image dataset and train a multi-level feature interaction and aggregation network. Use the hybrid similarity metric triplet loss function to calculate the minimum distance between the matching deep feature descriptor and the non-matching deep feature descriptor from the deep feature descriptors extracted by the multi-level feature interaction and aggregation network. Calculate the value of the hybrid similarity metric triplet loss function based on the distance and perform back propagation to optimize the network model, prevent gradient disappearance and overfitting, optimize the network model, and complete network training; S4: using the trained multi-level feature interaction and aggregation network to extract deep feature descriptors according to the key feature points detected in step S1, and using the nearest neighbor matching algorithm to match the deep feature descriptors to obtain an initial matching result; S5: Utilize a multi-level mismatch elimination strategy to eliminate mismatches from the initial matching result to obtain a final matching result; the multi-level mismatch elimination strategy includes adaptive threshold constraints, rough screening based on a grid motion statistics algorithm, and fine matching based on a marginalized sample consensus algorithm.

2. The multimodal remote sensing image matching method of coupling phase structure and depth features according to claim 1 is characterized in that: The specific method of step S1 is: S101: Image preprocessing of multimodal remote sensing images, including geometric correction, atmospheric correction and radiometric calibration; S102: calculating an edge confidence map of the preprocessed multimodal remote sensing image, and performing point shape texture filtering on the preprocessed multimodal remote sensing image according to the edge confidence map to obtain structural features of the multimodal remote sensing image; S103: performing phase consistency calculation on the filtered multimodal remote sensing image to obtain a maximum moment diagram of the multimodal remote sensing image and highlight edge information; S104: Use the Brisk algorithm to perform feature point detection on the maximum moment graph, obtain key feature point information, and obtain a stable feature graph; S105: With the key feature points detected in step S104 as the center, a number of image blocks of a size of 224×224 are cropped from the original multimodal remote sensing image to complete image block extraction.

3. The multimodal remote sensing image matching method of coupling phase structure and depth features according to claim 1, characterized in that: The convolutional neural network comprises five stages, the first stage is a 3×3 convolutional layer, the second stage, the third stage, the fourth stage and the fifth stage are formed by stacking multiple inverted residual structures, and the second stage, the third stage, the fourth stage and the fifth stage stack the inverted residual structure 3 times, 3 times, 7 times and 4 times respectively; wherein, the second stage is formed by stacking an inverted residual structure with a step length of 1, an inverted residual structure with a step length of 2 and an inverted residual structure with a step length of 1, the third stage is formed by stacking an inverted residual structure with a step length of 2 and two inverted residual structures with a step length of 1 in sequence, the fourth stage is formed by stacking an inverted residual structure with a step length of 2 and six inverted residual structures with a step length of 1 in sequence, and the fifth stage is formed by stacking an inverted residual structure with a step length of 2 and three inverted residual structures with a step length of 1; The Transformer network comprises five stages, the first stage is a 3×3 convolutional layer, the second stage is four inverted residual structures, and the third, fourth and fifth stages are each composed of an inverted residual structure and a MobileViT structure; wherein, the second stage is composed of three inverted residual structures with a step size of 1 and one inverted residual structure with a step size of 2 stacked together, and the step sizes of the inverted residual structures in the third, fourth and fifth stages are all 2.

4. The multimodal remote sensing image matching method of coupling phase structure and depth features according to claim 3 is characterized in that: The MobileViT structure includes three stages. The first stage is a 3×3 convolution layer and a 1×1 convolution layer. The second stage first flattens the feature map output by the first stage, and then extracts features through the Transformer module. The Transformer module is divided into two layers. The first layer normalizes the flattened features, and then calculates the normalized features by the self-attention mechanism. The initial input features of the first layer are added to the features output by the self-attention mechanism to obtain the output of the first layer. The second layer normalizes the features output by the first layer, and then performs multi-layer perceptron operations on the normalized features. Finally, the output features of the first layer are added to the output features of the multi-layer perceptron operation to obtain the output feature map of the final Transformer module. Then, the output feature map of the Transformer module is folded to the size of the input feature map of the second stage. The third stage first uses a 1×1 convolution layer to restore the feature map output by the second stage to the size of the input feature map of the first stage, and then splices the feature map input by the first stage with the feature map output by the 1×1 convolution on the channel, and finally uses a 3×3 convolution layer to fuse to obtain the output feature map.

5. The multimodal remote sensing image matching method of coupling phase structure and depth features according to claim 4 is characterized in that: The channel feature interaction module includes a coordinate attention mechanism, a global average pooling operation, a global maximum pooling operation, a channel splicing operation, a multi-layer perceptron operation and a Sigmoid activation function, and the spatial feature interaction module includes a channel splicing operation, a multi-layer perceptron operation and a Sigmoid activation function; The feature map obtained by the upper branch network of the dual-branch backbone network is , the feature map obtained by the lower branch network of the dual-branch backbone network is The channel feature interaction module first uses the coordinate attention mechanism to input feature maps and feature map Extract key features, and then perform global maximum pooling and global average pooling operations on the feature map along the channel to obtain four feature maps , , and ; Then, the four feature maps are concatenated to obtain the feature map , , Represents a 3D feature map; finally, the feature map is operated using a multi-layer perceptron Perform deep interaction, then use the Sigmoid activation function to calculate the weight of the feature map after deep interaction and divide it into two weight feature maps of the same size and , , , and then according to the weight feature map And weight feature map Get the feature map after channel correction and ; The calculation formula of the feature map extracted by the channel feature interaction module is expressed as: ， ， ， In the formula, Represents the function of the coordinate attention mechanism to calculate the feature map, Indicates a global pooling operation on the feature map. Indicates the global maximum pooling operation on the feature map. Indicates the channel splicing operation on the feature map. represents the multi-layer perceptron operation, Represents the Sigmoid activation function of the channel feature interaction module, represents the channel-by-channel splitting function, Represents the feature map of the upper branch network after channel correction, Represents the feature map of the lower branch network after channel correction, Indicates channel multiplication; The spatial feature interaction module first performs the input feature map and feature map Perform channel splicing operation, and then use multi-layer perceptron operation to realize feature interaction and obtain feature map , , and then use the Sigmoid activation function to calculate the feature map The weights of , and divide it into two weight feature maps of the same size and , , , and then according to the weight feature map And weight feature map Get the feature map after spatial correction and ; The calculation formula for extracting feature maps by the spatial feature interaction module is: ， ， ， In the formula, The Sigmoid activation function representing the spatial feature interaction module, Represents the feature map of the upper branch network after spatial correction, Represents the feature map of the lower branch network after spatial correction; After the feature maps obtained by different branches of the dual-branch backbone network are interacted and corrected by the channel feature interaction module and the spatial feature interaction module, the formula for obtaining the corrected feature map is expressed as: ， In the formula, It represents the feature map obtained by the upper branch network after the multi-level feature interaction and aggregation module based on the coordinate attention mechanism outputs the feature map. It represents the feature map obtained by the lower branch network after the multi-level feature interaction and aggregation module based on the coordinate attention mechanism outputs the feature map.

6. The multimodal remote sensing image matching method of coupling phase structure and depth features according to claim 5, characterized in that: The dense feature fusion module based on the coordinate attention mechanism includes a feature aggregation module and a feature fusion module. The feature aggregation module includes four 1×1 convolution layers, a 3×3 DW convolution layer, an activation function Relu, a batch normalization layer BN, a Softmax function and a Split function. The feature fusion module includes a coordinate attention mechanism, a 7×7 convolution layer, a batch normalization layer BN and Normalization; The feature aggregation module first performs multi-level feature interaction based on the coordinate attention mechanism and the feature map output by the aggregation module. and feature map Perform channel splicing operations and use four 1×1 convolutional layers, a 3×3 DW convolutional layer, an activation function Relu, a batch normalization layer BN and a Softmax function to achieve deep aggregation of deep features; Secondly, use the Softmax function and Split function to get two weight graphs and , and then use the weight map and weight graph Get the aggregated dense feature map and ; The calculation formula of the feature aggregation module aggregation feature map is expressed as: ， ， ， In the formula, They represent the feature maps output by two 1×1 convolutional layers, represents the Softmax function, Represents a composite function consisting of a 1×1 convolution layer, a 3×3 DW convolution layer, a batch normalization layer BN, and an activation function Relu. represents the dense feature map of the upper branch network output by the feature aggregation module, Represents the dense feature map of the branch network output under the feature aggregation module output; The feature fusion module is used to fuse the aggregated dense features to obtain a deep feature descriptor for matching. First, the dense feature map output by the feature aggregation module is and dense feature maps Perform channel splicing operations and use the coordinate attention mechanism to enhance the expression of key features and suppress irrelevant features; secondly, use a 7×7 convolution layer and a batch normalization layer BN to obtain a 256-dimensional feature vector , and the eigenvector conduct Normalization operation obtains deep feature descriptor ; The expression of the feature fusion module is: ， ， In the formula, Represents a composite function consisting of a 7×7 convolutional layer and a batch normalization layer BN. represents the deep feature descriptor, express Normalization operation.

7. The multimodal remote sensing image matching method of coupling phase structure and depth features according to claim 1, characterized in that: In step S3, the hybrid similarity metric triplet loss function is given by norm normalization and triplet loss, The calculation formula for norm normalization is: ， In the formula, Descriptor pair and conduct The result after norm normalization is and yes A pair of matched deep feature descriptors before norm normalization, Indicates the number of deep feature descriptors; The triplet loss is calculated as: ， ， In the formula, represents the triplet loss value, represents the angle between the matching deep feature descriptor pairs, represents the angle between non-matching deep feature descriptor pairs, represents a hybrid similarity measure for matching pairs of deep feature descriptors, Representing a hybrid similarity measure for pairs of non-matching deep feature descriptors; represents a mixed similarity measure, represents the cosine value of the angle between the deep feature descriptors, , A function representing the angle of the deep feature descriptor, , is a scalar used to adjust the ratio, ; is the maximum magnitude gradient normalization factor, is the angle between the deep feature descriptors, Represents the margin, which is a constant with a value of 1; Finally, the hybrid similarity metric triplet loss function is defined as: ， In the formula, represents the hybrid similarity metric triplet loss function, is the regularization parameter.

8. The multimodal remote sensing image matching method of coupling phase structure and depth features according to claim 1, characterized in that: The initial matching result obtained by the nearest neighbor matching algorithm contains the nearest matching point pair and the second nearest matching point pair in Euclidean distance. If the distance between the nearest matching point pair and the depth feature descriptor is Much smaller than the distance between the next nearest matching point and the deep feature descriptor , then the corresponding initial matching point pair is considered to have good quality. The adaptive threshold constraint calculates the mean of the Euclidean distances between the nearest matching point pair and the second nearest matching point pair in all initial matching results as the basis for judging whether to retain the corresponding matching point pair. The calculation formula for the mean of the Euclidean distance is: ； in, is the number of key feature points on the reference image, is the mean of the Euclidean distance, is the distance between the next neighboring matching point and the depth feature descriptor, is the distance between the nearest matching point pair and the depth feature descriptor. For each set of nearest matching point pairs and second nearest matching point pairs, when When , the corresponding nearest matching point pair is retained, otherwise the matching point pair is removed to filter out the initial correct matching point pair result; Then, the initial correct matching point pairs obtained after adaptive threshold constraint screening are subjected to motion smoothing constraint statistics based on the grid motion statistics algorithm to obtain the matching point pairs after rough screening; then, the marginalized sample consensus algorithm is used to solve the estimated value using iterative weighted least squares to achieve robust estimation, realize precise matching, and obtain the final correct matching points.

9. The multimodal remote sensing image matching method of coupling phase structure and depth features according to claim 8, characterized in that: The rough screening based on the grid motion statistics algorithm is to find the matching points that meet the conditions from the initial correct matching point pairs obtained after the adaptive threshold constraint screening, find the nearest neighborhood and the close neighborhood with the coordinates of the input initial matching point as the center, count the number of matches in the neighborhood, and then count the probabilities of correct matches and incorrect matches to obtain the probabilities of correct matches and incorrect matches in each neighborhood, and determine the correct and incorrect matches according to the probability threshold. The more the number of matching points in the image, the higher the degree of distinction between correct and incorrect matches. In the process of implementing the grid motion statistics algorithm, the grid is used instead of the neighborhood to realize the fast and efficient operation of the algorithm.

10. The multimodal remote sensing image matching method of coupling phase structure and depth features according to claim 9, characterized in that: The marginalized sample consensus algorithm is an iterative algorithm for correctly estimating the parameters of the mathematical model from a set of data containing erroneous matching points. The probability of a random point being a correct matching point is solved by marginalization, and the probability is used as a weight for weighted least squares optimization, thereby optimizing the model. The model with the largest number of correct matching points is found through iteration, and the final matching points are obtained.

Images