CN107909018A - A kind of sane multi-modal Remote Sensing Images Matching Method and system - Google Patents

Abstract

The invention discloses a robust multimodal remote sensing image matching method and system, which integrates various local feature descriptors for automatic multimodal remote sensing image matching. First, dense grid sampling is performed on the image, and the local feature descriptor of each grid point is extracted to form a three-dimensional dense feature expression map, and then a fast matching is established in the frequency domain using three-dimensional Fourier transform based on the feature expression map Similarity measurement, and finally template matching is used to identify the same-name point. In addition, for the invented matching method, a new dense feature expression method Gradient Oriented Feature Channel (CFOG) is proposed, which is superior to existing methods based on HOG, LSS, SIFT, SURF, etc. in terms of matching performance and computational efficiency. Dense feature representation. The invention can effectively overcome the non-linear radiation difference among multimodal images such as visible light, infrared, laser radar, synthetic aperture radar, and map, quickly and accurately identify points with the same name among the images, and realize automatic image matching.

Description

A Robust Multimodal Remote Sensing Image Matching Method and System

technical field

The invention relates to the technical field of satellite image processing, in particular to a multimodal remote sensing image matching method and system for automatic matching of multimodal images such as visible light, infrared, laser radar, contract aperture radar, and maps.

Background technique

Image matching is the process of identifying points with the same name between two or more images. It is the basic preprocessing step of many remote sensing image analyzes such as image fusion, change detection, image mosaic, etc. The matching accuracy is important for subsequent analysis work. influences. At present, remote sensing satellite sensors are equipped with global positioning system and inertial navigation system, which can directly locate and roughly match remote sensing images, eliminate obvious rotation and scale changes between images, so that only a certain amount (such as within dozens of pixels) exists between images translation difference. However, due to the different imaging mechanisms of images, there are significant nonlinear radiation differences among multi-modal remote sensing images (visible light, infrared, lidar and contract aperture radar, etc.), resulting in automatic matching of points with the same name is still very challenging.

At present, the matching methods of multi-modal remote sensing images can be mainly divided into: feature matching and area matching. Feature matching is to achieve image matching through the similarity of salient features between images. Commonly used features include point feature methods, line feature methods and surface features. Recently, local invariant features such as Scale Invariant Feature Transform (SIFT), shape context, etc. have also been applied in remote sensing image matching. However, these methods usually need to extract features with a high repetition rate between images, and for multimodal remote sensing with significant radiation differences, the repetition rate of feature extraction is often low, so they are suitable for automatic multimodal remote sensing images. Matching also has certain limitations.

The region-based method mainly adopts the strategy of template matching, and uses a certain matching similarity measure as a criterion to identify points with the same name between images. In this process, the choice of similarity measure is very important, which directly affects the subsequent matching accuracy. Commonly used similarity measures include sum of squares of gray difference, normalized correlation coefficient and mutual information. However, these similarity measures use the gray similarity between images to identify points with the same name, and because there are large differences in gray information between multimodal remote sensing images, they cannot be well applied to multimodal remote sensing images. automatic matching. Compared with grayscale information, the structure and shape attributes of images have higher similarity, and related research uses local descriptions such as Histogram of Oriented Gradient (HOG) and Local Self-Similarity (LSS). The structure and shape features of the image are extracted by the symbol, and a similarity measure is established on this basis for image matching, which improves the matching performance. Nevertheless, HOG and LSS only use the local information of the image to construct features, or use the neighborhood information of the extracted feature points to calculate features, which is a sparse feature expression method, which is difficult to reflect the commonality between multi-modal images. attributes, and its computational efficiency is relatively low.

In view of this, the present invention constructs a robust multimodal remote sensing image matching method, which can integrate various local feature descriptors for automatic multimodal remote sensing image matching. This method first performs dense grid sampling on the image, and extracts local feature descriptors such as HOG, LSS, SIFT or Speeded-Up Robust Features (SURF) for each grid point to form a three-dimensional dense feature expression map to reflect the image. shared structural and shape properties. Then based on the feature expression graph, the similarity measure of image matching is established, and the strategy of template matching is used to identify the same name point. In addition, for the invented method, a dense feature descriptor based on the direction gradient feature is constructed, named direction. Gradient channel feature (Channel Feature of Orientated Gradient, CFOG). It is superior to feature expression techniques such as dense HOG, LSS, SIFT and SURF in terms of matching performance and computational efficiency.

Contents of the invention

The purpose of the present invention is to overcome the shortcomings of traditional matching methods and provide a robust multi-modal remote sensing image matching method. This method uses the dense feature expression technology to extract the common structure and shape features between images, and establishes a fast matching similarity measure based on it, which can quickly and accurately obtain a large number of uniformly distributed point of the same name. In addition, for the invented method, a novel dense feature expression technique named Channel Feature of Orientated Gradient (CFOG) is constructed.

On the one hand, the present invention provides a robust multimodal remote sensing image matching method, comprising the following steps:

A. Determine the resolution information of the reference image and the input image. If the two images have the same resolution, perform subsequent processing. If the resolutions are different, sample the two images to the same resolution;

B. Using the block strategy, detect a series of evenly distributed feature points on the reference image, which are recorded as P _1i (i=1,2,3,...,N), and select the template area centered on point P _1i AreaW _1i ;

C. According to the geographical coordinate information provided by the remote sensing image itself, predict the matching area AreaW _2i corresponding to the point set P _1i (i=1,2,3,...,N) on the input image;

D. Perform dense grid sampling in the matching area and construct a three-dimensional dense feature expression map;

E. On the basis of the three-dimensional dense feature expression map, establish a similarity measure to match the same-name points;

F. For the obtained point with the same name, perform local extremum fitting on its similarity map, and solve the sub-pixel position of the matching point;

G. Repeat steps (C)—(F), traverse each point of P _1i (i=1,2,3,...,N), and obtain a point pair {PD _1i(x,y) with the same name with sub-pixel precision , PD _2i(x,y) }(i=1,2,3,...,N);

H. Eliminate {PD _1i(x,y) ，PD _2i(x,y) }(i=1,2,3,...,N) with the same name point pair with large error to obtain the final same name point pair {PID _1i(x,y) , PID _2i(x,y) }(i=1,2,3,...,S).

Wherein, the step D includes the following steps: performing dense grid sampling on the image data in the matching area, and calculating local feature descriptors such as HOG, LSS, SIFT or SURF of the grid points to form a three-dimensional dense feature expression map.

Further, the step D can also be to construct a channel feature of direction gradient (CFOG) in the matching area, which specifically includes the following steps:

D1. Perform dense grid sampling on the image data in the matching area, calculate the gradient information of each grid point in each direction, and form a three-dimensional direction gradient map;

D2. In the horizontal X direction and vertical Y direction, use a Gaussian filter or a triangular filter to perform convolution operations on the three-dimensional directional gradient map to obtain a feature map Then use the one-dimensional filter [1,2,1] to map the feature map in the Z direction Perform convolution operation to obtain feature map

D3. Pair feature map Perform a normalization operation to obtain the final directional gradient channel feature map.

Wherein said step D is to construct the direction gradient channel feature (CFOG) in the matching area, and the specific calculation process includes the following steps:

For all grid points in the area, use one-dimensional filters [-1,0,1] and [-1,0,1] ^T to calculate their gradients in the horizontal direction (X direction) and vertical direction (Y direction) g _x and g _y ;

Using g _x and g _y to calculate their gradient value g _θ in each direction, the calculation formula is as follows:

In the formula, θ represents the gradient direction of quantization, abs represents the absolute value, When the sign indicates that the value is positive, it takes itself, otherwise it takes 0;

Stack g _θ in all directions to form a three-dimensional directional gradient map g _o , and then use a two-dimensional Gaussian filter or a triangular filter with a standard σ in the X and Y directions to perform a convolution operation on g _o to obtain a feature map Then use the one-dimensional filter [1,2,1] in the Z direction to Perform convolution operation to obtain feature map

feature map Each pixel in Z corresponds to a feature vector f _i in the Z direction, traverse each pixel, and perform a normalization operation on its feature vector v _i to obtain the final directional gradient channel feature map, the normalized calculation formula as follows:

where ε is a number that avoids division by zero.

Wherein, the gradient direction is evenly divided into n equal parts in the range of 360°, and the angular interval of each equal part is 360/n°.

Further, the step E is to establish a similarity measure based on the three-dimensional dense feature expression map to obtain the matching position of the image.

Wherein said step E specifically comprises the following calculation steps:

After step D, the three-dimensional dense feature expression graphs D ₁ and D ₂ of AreaW _1i and AreaW _2i are respectively obtained, and D ₁ is used as a template to slide on D ₂ , and the similarity between D ₁ and D ₂ is used for matching, similar The sex measure can be difference sum of squares, Euclidean distance, normalized cross-correlation, phase correlation, etc.

The formula for calculating the sum of squared differences is:

The formula for calculating the Euclidean distance is:

The formula for calculating the normalized cross-correlation is:

The formula for phase correlation is:

In formulas (3)-( ₆ ), _c represents the coordinates in the _{three -} dimensional dense feature representation map, v represents the offset between D1 and D2, S _i represents the similarity measure between D1 and D2, F represents the fast Fourier forward transform, and F ^* represents the complex conjugate of F. For the difference sum of squares and Euclidean distance, when S _i takes the minimum value, the offset v _i between D ₁ and D ₂ will be obtained. While for normalized cross _- correlation and phase correlation, when S _i takes the maximum value _, the offset _v between D1 and D2 will be obtained

In addition, in order to perform fast matching, the present invention constructs a similarity measure based on fast Fourier transform, the invention of the measure originates from the difference sum of squares, as mentioned above, when the difference square sum is minimum, the matching position v _i will be obtained, the calculation formula as follows

Expand the formula (7) to get:

In formula (8), since the values of the first and second terms are close to constants, when the value of the third term is the largest, formula (7) will obtain the minimum value, therefore, the similarity function can be redefined as:

In the formula, is a convolution operation;

Considering that the point multiplication operation in the frequency domain is equivalent to the convolution operation in the space domain, fast Fourier transform is used in the frequency domain to improve its computational efficiency, and the similarity function based on Fourier transform is obtained as:

In the formula, F and F ^-1 represent fast Fourier forward transform and inverse transform respectively, and F ^* represents the complex conjugate of F, where F can be two-dimensional or three-dimensional Fourier transform.

On the other hand, the present invention also provides a robust multimodal remote sensing image matching system, which system includes the following units:

The preprocessing unit is used to judge the resolution information of the reference image and the input image. If the two images have the same resolution, execute the subsequent unit. If the resolutions are different, sample the two images to the same resolution;

The template area selection unit is used to detect a series of evenly distributed feature points on the reference image by adopting the block strategy, which is denoted as P _1i (i=1,2,3,...,N), and the point P _1i is Select the template area AreaW _1i in the center;

The matching area selection unit is used to predict the matching area AreaW _2i corresponding to the point set P _1i (i=1,2,3,...,N) on the input image according to the geographical coordinate information provided by the remote sensing image itself;

A feature extraction unit is used to construct a three-dimensional dense feature expression map in the matching area;

The initial matching unit, on the basis of the three-dimensional dense feature expression map, establishes a similarity measure to match the points of the same name; for the obtained points of the same name, it performs local extremum fitting on the similarity map, and solves the sub-pixel position of the matching point; Repeat the operation of the above unit, traverse each point of P _1i (i=1,2,3,…,N), and get the same-name point pair {PD _1i(x,y) with sub-pixel precision, PD _{2i(x, y)} }(i=1,2,3,...,N);

The matching point screening unit is used to eliminate the same-name point pairs with large errors in {PD _1i(x,y) , PD _2i(x,y) }(i=1,2,3,...,N) to obtain the final Point pairs with the same name {PID _1i(x,y) , PID _2i(x,y) } (i=1, 2, 3, . . . , S).

Further, the feature extraction unit is used to perform dense grid sampling on the image data in the matching area, calculate the local feature descriptor of each grid point, and form a three-dimensional dense feature expression map.

Further, the initial matching unit establishes a similarity measure for image matching on the basis of a three-dimensional dense feature expression map, and performs calculations to obtain a similarity map, and takes the extreme value position of the similarity map as the matching position of the image.

In summary, owing to adopting above-mentioned technical scheme, the beneficial effect of the present invention is:

(1) The present invention constructs a robust multi-modal remote sensing image matching method, which extracts the local feature descriptor (such as HOG, LSS, SIFT or SURF, etc.) of each grid point by performing dense grid sampling on the image , to form a three-dimensional dense feature expression map, which can better reflect the common structure and shape attributes among multi-modal remote sensing images, and establish a similarity measure for image matching on the three-dimensional dense feature expression map. This method can quickly, accurately and automatically obtain a large number of uniformly distributed points of the same name among multi-modal images, which can effectively improve the actual production efficiency of matching and meet the needs of business operations. Moreover, it is a general method that can integrate various A local feature description is used for image matching (not limited to descriptors such as HOG, LSS, SIFT or SURF).

(2) In the construction method, the three-dimensional dense feature expression is formed by densely sampling the image and calculating the local feature description of each grid point such as HOG, LSS, SIFT or SURF, which is a kind of dense feature expression. This is different from traditional descriptors such as HOG, LSS, SIFT, or SURF. They only use the local area information of the image for feature construction, or use the neighborhood information of the extracted feature points to calculate features, which is a sparse feature expression method. In contrast, the dense 3D feature expression technology of the present invention can better reflect the common structure and shape attributes between multi-modal remote sensing images more accurately, and the matching performance is more robust. The similarity measure of expression technology can realize fast and accurate matching between multi-modal images.

(3) For the constructed method, the present invention constructs a novel dense feature descriptor-directional gradient channel feature (CFOG) using the directional gradient information of the image, which is superior to dense HOG, CFOG in terms of matching efficiency and accuracy. Feature expression methods such as LSS, SIFT, and SURF.

(4) On the basis of three-dimensional dense feature expression, the sum of squared differences, Euclidean distance, normalized cross-correlation, phase correlation, etc. can be used as similarity measures (not limited to these similarity measures) to carry out homonym points. In addition, the present invention also establishes a A fast similarity measure based on Fourier transform, which has higher calculation efficiency and better calculation effect.

(5) A large number of experimental results show that for images in flat areas, the overall matching accuracy can reach within 1 pixel, and for images in mountainous and urban areas, the overall matching accuracy can reach within 2.5 pixels. For large-scale remote sensing images (more than 20000×20000 pixels), image matching can be completed within 30 seconds. Moreover, compared with the currently popular commercial remote sensing image software (ENVI and ERDAS), the invented method has advantages in both matching accuracy and computational efficiency.

Description of drawings

Fig. 1 is the overall flowchart of the present invention.

Fig. 2 is a schematic diagram of three-dimensional dense feature expression in the present invention.

Fig. 3 is the construction process of the Channel Feature of Oriented Gradient (CFOG) in the present invention.

Detailed ways

In order to make those skilled in the art better understand the technical solution of the present invention, the technical solution of the present invention is clearly and completely described below in conjunction with the accompanying drawings of the present invention. Based on the embodiments in this application, those of ordinary skill in the art will Other similar embodiments obtained without creative work shall all fall within the scope of protection of this application.

As shown in Figure 1, it is a robust multimodal remote sensing image matching method, including the following steps:

Step A, according to the resolution information of the reference image and the input image, it is judged whether the resolutions of the two images are consistent, if they are consistent, follow-up processing is performed, and if they are inconsistent, the two images are sampled according to the same resolution.

Step B, using the block strategy, using the Harris or Forstner operator to extract a large number of uniformly distributed feature points on the reference image, specifically including:

Divide the reference image into n×n non-overlapping regions, calculate the Harris or Forstner eigenvalues of each pixel in each region, sort the eigenvalues from large to small, and select k with larger eigenvalues pixels as feature points. In this way, k feature points can be detected in each region, and there are n×n×k feature points on the entire image. The values of n and k can be set according to actual needs. Denote the feature point set detected on the reference image as P _1i (i=1, 2, 3, . . . , N).

In other embodiments, other operators may also be used for image feature extraction, which is not limited in the present invention.

Step C, according to the geographical coordinate information provided by the remote sensing image itself, predict the matching search area corresponding to the point set P _1i (i=1,2,3,...,N) on the input image, specifically including the following steps:

(1) extract a point P _1i (x, y) of P _1i (i=1,2,3,...,N), x and y represent the image coordinates of point P _1i , and point P _1i (x, y ) as the center to select a template area AreaW _1i with a size of M×M, and obtain the geographic coordinates Geo _i corresponding to the point;

(2) According to the geographic coordinate Geo _i and combined with the geographic coordinate information of the input image, calculate its corresponding image coordinate P _2i (x,y) on the input image, and determine a size centered on P _2i (x,y) A square window of M ₁ ×M ₁ is used as the matching search area AreaW _2i .

Step D, constructing a three-dimensional dense feature expression map in AreaW _1i and AreaW _2i , see FIG. 2 .

In one embodiment, step D comprises the following steps:

Sampling the pixels in the area with a dense equidistant grid, and calculating local feature descriptors such as HOG, LSS, SIFT or SURF for each grid point, and arranging the feature vectors corresponding to each grid point in the Z direction, A three-dimensional dense feature representation map is formed. For the grid sampling interval η∈[0.5,4] pixels, the smaller the sampling interval, the higher the matching accuracy, but the greater the computational complexity. When the sampling interval is 1 pixel, that is, sampling pixel by pixel, the dense feature expression map formed is the pixel-by-pixel feature expression map. If the sampling interval is a decimal, it is also necessary to interpolate the pixel values in the neighborhood of the grid point to facilitate the calculation of local feature descriptors. In addition, the specific local feature descriptor used is not limited in this embodiment.

In another embodiment, step D is to construct a directional gradient channel feature (Channel Feature of Orientated Gradient, CFOG) in the areas AreaW _1i and AreaW _2i , referring to Figure 3, specifically comprising the following steps:

(1) For all pixels in the area, use one-dimensional filters [-1,0,1] and [-1,0,1] ^T to calculate their horizontal direction (X direction) and vertical direction (Y direction) Gradients g _x and g _y of .

(2) Using g _x and g _y to calculate their gradient value g _θ in each direction, the calculation formula is as follows:

In the formula, θ represents the quantized gradient direction. Here, the gradient direction is evenly divided into n equal parts within the range of 360°, and the angular interval of each equal part is (360/n)°. abs means to take the absolute value, When the sign indicates that the value is positive, it takes itself, otherwise it takes 0.

(3) Stack g _θ in all directions together to form a three-dimensional direction gradient map g _o . Then in the X and Y directions, use a two-dimensional Gaussian filter or a triangular filter with a standard of σ to perform a convolution operation on g _o to obtain a feature map Then use the one-dimensional filter [1,2,1] in the Z direction to Perform convolution operation to obtain feature map

(4) Feature map Each pixel in corresponds to a feature vector f _i in the Z direction. Each pixel is traversed, and its feature vector v _i is normalized to further eliminate the influence of illumination changes, and the final CFOG feature map is obtained. The normalized calculation formula is as follows:

where ε is a number that avoids division by zero.

Step E, on the basis of the three-dimensional dense feature expression graph, establish a similarity measure to match the points with the same name, specifically including the following steps:

(1) After step D, the dense feature expression graphs D ₁ and D ₂ corresponding to areas AreaW _1i and AreaW _2i can be obtained respectively, and D ₁ is used as a template to slide on D ₂ , using the similarity between D ₁ and D ₂ The similarity measure can be sum of squared difference, Euclidean distance, normalized cross-correlation, phase correlation, etc. In addition, the specific similarity measure used is not limited in this embodiment.

The formula for calculating the sum of squared differences is:

The formula for calculating the Euclidean distance is:

The formula for calculating the normalized cross-correlation is:

The formula for phase correlation is:

In formulas (3)-( ₆ ), _c represents the coordinates in the _{three -} dimensional dense feature representation map, v represents the offset between D1 and D2, S _i represents the similarity measure between D1 and D2, F represents the fast Fourier forward transform, and F ^* represents the complex conjugate of F. For the difference sum of squares and Euclidean distance, when S _i takes the minimum value, the offset v _i between D ₁ and D ₂ will be obtained. For normalized cross-correlation and phase correlation, when S _i reaches the maximum value, the offset v _i between D ₁ and D ₂ will be obtained

In addition, the present invention constructs a fast similarity measure based on Fast Fourier Transform. The invention of this measure originates from the sum of squared differences. As mentioned above, when the sum of squared differences is minimized, the matching position v _i will be obtained. The calculation formula is as follows

Expand the formula (7) to get:

In formula (8), since the values of the first and second terms are close to constants, when the value of the third term is the largest, formula (7) will obtain the minimum value, therefore, the similarity function can be redefined as:

In the formula, is a convolution operation;

Considering that the point multiplication operation in the frequency domain is equivalent to the convolution operation in the space domain, fast Fourier transform is used in the frequency domain to improve its computational efficiency, and the similarity measure based on Fourier transform is obtained as:

In the formula, F and F ^-1 represent fast Fourier forward transform and inverse transform respectively, and F ^* represents the complex conjugate of F, where F can be two-dimensional or three-dimensional Fourier transform. The matching process using formula (10) is as follows: firstly perform fast Fourier transform on D ₁ to obtain DF(D ₁ ), then perform Fourier transform on D ₂ and take complex conjugate to obtain DF ^* (D ₂ ), Then perform dot multiplication of DF(D ₁ ) and DF ^* (D ₂ ), and inverse Fourier transform the result to obtain the similarity graph between D ₁ and D ₂ , and the position of the maximum value of the similarity graph Then it corresponds to the offset v _i between D ₁ and D ₂ , that is, the offset between point P _1i (x,y) and point P _2i (x,y).

Record the offset of v _i in the X and Y directions obtained according to the above measurement as (Δx, Δy), then the point with the same name corresponding to the point P _1i (x, y) is P _2i (x-Δx,y- Δy), denoted as P ^* _2i (x, y), and the obtained point pair with the same name is {P _1i (x, y), P ^* _2i (x, y)}.

Step F, for the above point pairs with the same name {P _1i (x, y), P* _2i (x, y)}, use binary quadratic polynomials to perform local interpolation to obtain sub-pixel accuracy, specifically including the following steps:

(1) Select a local small window of 3×3 pixels centered on the point P* _2i (x,y), and count the similarity values of all pixels in the window, and then use the binary quadratic polynomial according to the principle of least squares Establish the corresponding relationship between the similarity value and the pixel position;

(2) Calculate the partial derivative of the binary quadratic polynomial, and find the position where the partial derivative is 0, that is, obtain the sub-pixel precision point pair {PD _1i (x, y), PD _2i (x, y)};

Repeat step CF to traverse each point in P _1i (i=1,2,3,...,N), and obtain the same-named point pair {PD _1i (x,y), PD _2i (x, y)} (i=1,2,3,...,N).

Step G: Eliminate {PD _1i (x, y), PD _2i (x, y)} (i=1,2,3,….,N) with the same name point pair with large error to obtain the final same name point pair , including the following steps:

(1) Based on the principle of least squares, use the coordinates of {PD _1i (x, y), PD _2i (x, y)} (i=1,2,3,...,N) to establish a projection transformation model.

(2) Calculate the root mean square error and residual error of the point pairs with the same name, and remove the point pair with the same name with the largest residual error.

(3) Repeat the above two steps until the root mean square error is less than 1.5 pixels, and obtain the final point pair {PID _1i (x, y), PID _2i (x, y)} (i= 1,2,3,...,S).

Step H, remove {PD _1i(x,y) ，PD _2i(x,y) }(i=1,2,3,….,N) with the same name point pair with large error to obtain the final same name point pair {PID _1i(x,y) , PID _2i(x,y) }(i=1,2,3,...,S), and used for matching.

In another embodiment, the present invention also provides a robust multimodal remote sensing image matching system, which includes the following units:

The preprocessing unit is used to judge the resolution information of the reference image and the input image. If the two images have the same resolution, execute the subsequent unit. If the resolutions are different, sample the two images to the same resolution;

The template area selection unit is used to detect a series of evenly distributed feature points on the reference image by adopting the block strategy, which is denoted as P _1i (i=1,2,3,...,N), and the point P _1i is Select the template area AreaW _1i in the center;

The matching area selection unit is used to predict the matching area AreaW _2i corresponding to the point set P _1i (i=1,2,3,...,N) on the input image according to the geographical coordinate information provided by the remote sensing image itself;

The feature extraction unit is used for three-dimensional dense feature expression map in the matching area;

The initial matching unit is used to establish a similarity measure on the basis of the three-dimensional dense feature expression map to match the points with the same name; for the obtained points with the same name, perform local extremum fitting on the similarity map to solve the sub-pixels of the matching points Position; repeat the operation of the above unit, traverse each point of P _1i (i=1,2,3,...,N), and obtain the same-named point pair {PD _1i(x,y) with sub-pixel precision, PD _{2i( x,y)} }(i=1,2,3,...,N);

The matching point screening unit is used to eliminate the same-name point pairs with large errors in {PD _1i(x,y) , PD _2i(x,y) }(i=1,2,3,...,N) to obtain the final Point pairs with the same name {PID _1i(x,y) , PID _2i(x,y) } (i=1, 2, 3, . . . , S).

Further, the feature extraction unit is used to perform dense grid sampling on the image data in the matching area, calculate the local feature descriptor of each grid point, and form a three-dimensional dense feature expression map.

Further, the initial matching unit uses the basis of three-dimensional dense feature expression to establish a similarity measure to perform correlation operations to obtain a similarity map, and take the position of the extreme value of the similarity map as the matching position of the image.

The above is the description of the specific implementation of the present invention. Through the matching method constructed in the present invention, various local feature descriptors (such as HOG, LSS, SIFT, or SURF, etc.) can be used for effective three-dimensional dense feature expression, thereby effectively reflecting Compared with similarity measures based on gray information, such as normalized correlation coefficient and mutual information, the calculation efficiency is higher by sharing structure and shape features between images, and establishing a similarity measure for image matching on this basis. In addition, the Oriented Gradient Channel Feature (CFOG) constructed by the present invention is a novel three-dimensional dense feature expression, which is superior to dense feature expressions based on descriptors such as HOG, LSS, SIFT, and SURF in terms of matching efficiency and accuracy. technology. The technical solution of the present invention can make up for the deficiency that the traditional matching method is sensitive to nonlinear radiation differences between multimodal images, and can effectively solve the matching of multimodal remote sensing data such as visible light, infrared, laser radar, synthetic aperture radar, and maps. problem.

The technical solution proposed by the present invention is a general technical method, which can integrate various local feature descriptors (not limited to CFOG, HOG, LSS, SIFT and SURF, etc.) to express three-dimensional dense features, and can use Various similarity measures (not limited to difference sum of squares, Euclidean distance, normalized cross-correlation, phase correlation, and similarity measures based on fast Fourier transform constructed by the present invention, etc.) for image matching

The present invention is not limited to the foregoing specific embodiments. The present invention extends to any new feature or any new combination disclosed in this specification, and any new method or process step or any new combination disclosed.

Claims (10)

1. A robust multimodal remote sensing image matching method, characterized in that it comprises the following steps: A. Determine the resolution information of the reference image and the input image. If the two images have the same resolution, proceed to step B. If the resolutions are different, sample the two images to the same resolution; B. Using the block strategy, detect a series of evenly distributed feature points on the reference image, which are recorded as P _1i (i=1,2,3,...,N), and select the template area centered on point P _1i AreaW _1i ; C. According to the geographical coordinate information provided by the remote sensing image itself, predict the matching area AreaW _2i corresponding to the point set P _1i (i=1,2,3,...,N) on the input image; D. Construct a three-dimensional dense feature expression map in the matching area; E. On the basis of the pixel-by-pixel feature expression map, establish a similarity measure to match the points with the same name; F. For the obtained point with the same name, perform local extremum fitting on its similarity map, and solve the sub-pixel position of the matching point; G.. Repeat step CF, traversing each point of P _1i (i=1,2,3,...,N), and obtain the same-named point pair {PD _1i(x,y) with sub-pixel precision, PD _{2i( x,y)} }(i=1,2,3,...,N); H. Eliminate {PD _1i(x,y) ，PD _2i(x,y) }(i=1,2,3,….,N) with the same name point pair with large error to obtain the final same name point pair{ PID _{1i(x, y)} , PID _{2i(x, y)} } (i=1, 2, 3, . . . , S). 2. multimodal remote sensing image matching method as claimed in claim 1, is characterized in that, described step D comprises the steps: Dense grid sampling is performed on the image data in the matching area, the local feature descriptor of each grid point is calculated, and the corresponding feature vectors are arranged in the Z direction to form a three-dimensional dense feature expression map. When the grid sampling interval is 1 pixel, pixel-by-pixel sampling is performed, and the formed three-dimensional dense feature expression map is a three-dimensional pixel-by-pixel feature expression map. 3. The multimodal remote sensing image matching method according to claim 2, wherein the local feature descriptor is one of HOG, LSS, SURF or Scale Invariant Feature Transform (SIFT). 4. The multimodal remote sensing image matching method according to claim 1, wherein the step D is to construct a channel feature of oriented gradients (CFOG) in the matching area, specifically comprising the following steps: D1. Perform dense grid sampling on the image data in the matching area, calculate the gradient information of each grid point in each direction, and form a three-dimensional direction gradient map; D2. In the horizontal X direction and vertical Y direction, use a Gaussian filter or a triangular filter to perform convolution operations on the three-dimensional directional gradient map to obtain a feature map Then use the one-dimensional filter [1,2,1] to map the feature map in the Z direction Perform convolution operation to obtain feature map D3. Pair feature map Perform a normalization operation to obtain the final directional gradient channel feature map. 5. multimodal remote sensing image matching method as claimed in claim 4, is characterized in that, the construction direction gradient channel feature (CFOG) in the described step D, specifically comprises the following computing steps: For all pixels in the area, use the one-dimensional filter [-1,0,1] and [-1,0,1] ^T to calculate their gradient g in the horizontal direction (X direction) and vertical direction (Y direction) _x and g _y ; Using g _x and g _y to calculate their gradient value g _θ in each direction, the calculation formula is as follows: In the formula, θ represents the gradient direction of quantization, abs represents the absolute value, When the sign indicates that the value is positive, it takes itself, otherwise it takes 0; Stack g _θ in all directions to form a three-dimensional direction gradient map g _o , and then use a two-dimensional Gaussian filter with standard σ to perform convolution operation on g _o in the X and Y directions to obtain the feature map Then use the one-dimensional filter [1,2,1] in the Z direction to Perform convolution operation to obtain feature map feature map Each pixel in Z corresponds to a feature vector f _i in the Z direction, traverse each pixel, and perform a normalization operation on its feature vector v _i to obtain the final directional gradient channel feature map (CFOG), and normalize The calculation formula is as follows: <mrow><msub><mi>f</mi><mi>i</mi></msub><mo>=</mo><mfrac><msub><mi>f</mi><mi>i</mi></msub><msqrt><mrow><mo>|</mo><mo>|</mo><msub><mi>f</mi><mi>i</mi></msub><mo>|</mo><msup><mo>|</mo><mn>2</mn></msup><mo>+</mo><mi>&amp;epsiv;</mi></mrow></msqrt></mfrac><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>2</mn><mo>)</mo></mrow></mrow> where ε is a number that avoids division by zero. 6. The multimodal remote sensing image matching method according to claim 1, wherein said step E is based on the three-dimensional dense feature expression map, using the difference sum of squares, Euclidean distance, normalized cross-correlation, phase One of the correlations is used as a similarity measure for homonym matching, or a fast Fourier transform-based similarity measure is used for homonym matching. 7. The multimodal remote sensing image matching method according to claim 6, wherein said step E specifically comprises the following calculation steps: After step D, the dense three-dimensional feature expression graphs D ₁ and D ₂ of AreaW _1i and AreaW _2i are respectively obtained, and D ₁ is used as a template to slide on D ₂ , and the similarity between D ₁ and D ₂ is used for matching. ; The calculation formula of the difference sum of squares is: <mrow><msubsup><mi>S</mi><mi>i</mi><mn>1</mn></msubsup><mrow><mo>(</mo><mi>v</mi><mo>)</mo></mrow><mo>=</mo><munder><mo>&amp;Sigma;</mo><mi>c</mi></munder><msup><mrow><mo>&amp;lsqb;</mo><msub><mi>D</mi><mn>1</mn></msub><mrow><mo>(</mo><mi>c</mi><mo>)</mo></mrow><mo>-</mo><msub><mi>D</mi><mn>2</mn></msub><mrow><mo>(</mo><mi>c</mi><mo>-</mo><mi>v</mi><mo>)</mo></mrow><mo>&amp;rsqb;</mo></mrow><mn>2</mn></msup><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>3</mn><mo>)</mo></mrow></mrow> The formula for calculating the Euclidean distance is: <mrow><msubsup><mi>S</mi><mi>i</mi><mn>2</mn></msubsup><mrow><mo>(</mo><mi>v</mi><mo>)</mo></mrow><mo>=</mo><munder><mo>&amp;Sigma;</mo><mi>c</mi></munder><msup><mrow><mo>&amp;lsqb;</mo><msub><mi>D</mi><mn>1</mn></msub><mrow><mo>(</mo><mi>c</mi><mo>)</mo></mrow><mo>-</mo><msub><mi>D</mi><mn>2</mn></msub><mrow><mo>(</mo><mi>c</mi><mo>-</mo><mi>v</mi><mo>)</mo></mrow><mo>&amp;rsqb;</mo></mrow><mn>2</mn></msup><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>4</mn><mo>)</mo></mrow></mrow> The formula for calculating the normalized cross-correlation is: <mrow><msubsup><mi>S</mi><mi>i</mi><mn>3</mn></msubsup><mrow><mo>(</mo><mi>v</mi><mo>)</mo></mrow><mo>=</mo><mfrac><mrow><munder><mo>&amp;Sigma;</mo><mi>c</mi></munder><mrow><mo>(</mo><msub><mi>D</mi><mn>1</mn></msub><mo>(</mo><mi>c</mi><mo>)</mo><mo>-</mo><msub><mover><mi>D</mi><mo>&amp;OverBar;</mo></mover><mn>1</mn></msub><mo>)</mo></mrow><mrow><mo>(</mo><msub><mi>D</mi><mn>2</mn></msub><mo>(</mo><mrow><mi>c</mi><mo>-</mo><mi>v</mi></mrow><mo>)</mo><mo>-</mo><msub><mover><mi>D</mi><mo>&amp;OverBar;</mo></mover><mn>2</mn></msub><mo>(</mo><mrow><mi>c</mi><mo>-</mo><mi>v</mi></mrow><mo>)</mo><mo>)</mo></mrow></mrow><msqrt><mrow><munder><mo>&amp;Sigma;</mo><mi>c</mi></munder><msup><mrow><mo>(</mo><msub><mi>D</mi><mn>1</mn></msub><mo>(</mo><mi>c</mi><mo>)</mo><mo>-</mo><msub><mover><mi>D</mi><mo>&amp;OverBar;</mo></mover><mn>1</mn></msub><mo>)</mo></mrow><mn>2</mn></msup><munder><mo>&amp;Sigma;</mo><mi>c</mi></munder><msup><mrow><mo>(</mo><msub><mi>D</mi><mn>2</mn></msub><mo>(</mo><mrow><mi>c</mi><mo>-</mo><mi>v</mi></mrow><mo>)</mo><mo>-</mo><msub><mover><mi>D</mi><mo>&amp;OverBar;</mo></mover><mn>2</mn></msub><mo>(</mo><mrow><mi>c</mi><mo>-</mo><mi>v</mi></mrow><mo>)</mo><mo>)</mo></mrow><mn>2</mn></msup></mrow></msqrt></mfrac><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>5</mn><mo>)</mo></mrow></mrow> The calculation formula of the phase correlation is: In formulas (3)-( ₆ ), _c represents the coordinates in the _{three -} dimensional dense feature representation map, v represents the offset between D1 and D2, S _i represents the similarity measure between D1 and D2, F represents the fast Fourier forward transform, and F ^* represents the complex conjugate of F; when the similarity measure is the difference sum of squares or Euclidean distance, when S _i takes the minimum value, the distance between D ₁ and D ₂ will be obtained Offset v _i ; when the similarity measure is normalized cross-correlation or phase correlation, when S _i reaches the maximum value, the offset v _i between D ₁ and D ₂ will be obtained; The fast Fourier transform-based similarity measure is to obtain the matching position v _i when the sum of squared differences is the smallest, and the calculation formula is as follows: <mrow><msub><mi>v</mi><mi>i</mi></msub><mo>=</mo><munder><mi>argmin</mi><mi>v</mi></munder><mo>{</mo><munder><mo>&amp;Sigma;</mo><mi>c</mi></munder><msup><mrow><mo>&amp;lsqb;</mo><msub><mi>D</mi><mn>1</mn></msub><mrow><mo>(</mo><mi>c</mi><mo>)</mo></mrow><mo>-</mo><msub><mi>D</mi><mn>2</mn></msub><mrow><mo>(</mo><mi>c</mi><mo>-</mo><mi>v</mi><mo>)</mo></mrow><mo>&amp;rsqb;</mo></mrow><mn>2</mn></msup><mo>}</mo><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>7</mn><mo>)</mo></mrow></mrow> Expand the formula (7) to get: <mrow><msub><mi>v</mi><mi>i</mi></msub><mo>=</mo><munder><mrow><mi>arg</mi><mi>min</mi></mrow><mi>v</mi></munder><mo>{</mo><munder><mo>&amp;Sigma;</mo><mi>c</mi></munder><msubsup><mi>D</mi><mn>1</mn><mn>2</mn></msubsup><mrow><mo>(</mo><mi>c</mi><mo>)</mo></mrow><mo>+</mo><munder><mo>&amp;Sigma;</mo><mi>c</mi></munder><msubsup><mi>D</mi><mn>2</mn><mn>2</mn></msubsup><mrow><mo>(</mo><mi>c</mi><mo>-</mo><mi>v</mi><mo>)</mo></mrow><mo>-</mo><mn>2</mn><munder><mo>&amp;Sigma;</mo><mi>c</mi></munder><msub><mi>D</mi><mn>1</mn></msub><mrow><mo>(</mo><mi>c</mi><mo>)</mo></mrow><mo>&amp;CenterDot;</mo><msub><mi>D</mi><mn>2</mn></msub><mrow><mo>(</mo><mi>c</mi><mo>-</mo><mi>v</mi><mo>)</mo></mrow><mo>}</mo><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>8</mn><mo>)</mo></mrow></mrow> In formula (8), since the values of the first and second terms are close to constants, when the value of the third term is the largest, formula (7) will obtain the minimum value, therefore, the similarity function can be redefined as: <mrow><msub><mi>v</mi><mi>i</mi></msub><mo>=</mo><munder><mrow><mi>arg</mi><mi>max</mi></mrow><mi>v</mi></munder><mo>&amp;lsqb;</mo><munder><mo>&amp;Sigma;</mo><mi>x</mi></munder><msub><mi>D</mi><mn>1</mn></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mo>&amp;CenterDot;</mo><msub><mi>D</mi><mn>2</mn></msub><mrow><mo>(</mo><mi>x</mi><mo>-</mo><mi>v</mi><mo>)</mo></mrow><mo>&amp;rsqb;</mo><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>9</mn><mo>)</mo></mrow></mrow> In the formula, is a convolution operation; Considering that the point multiplication operation in the frequency domain is equivalent to the convolution operation in the space domain, fast Fourier transform is used in the frequency domain to improve its computational efficiency, and the similarity measure based on Fourier transform is obtained as: <mrow><msub><mi>v</mi><mi>i</mi></msub><mo>=</mo><munder><mrow><mi>arg</mi><mi>max</mi></mrow><mi>v</mi></munder><mo>{</mo><msup><mi>F</mi><mrow><mo>-</mo><mn>1</mn></mrow></msup><mo>&amp;lsqb;</mo><mi>F</mi><mrow><mo>(</mo><msub><mi>D</mi><mn>1</mn></msub><mo>(</mo><mi>c</mi><mo>)</mo><mo>)</mo></mrow><mo>&amp;CenterDot;</mo><msup><mi>F</mi><mo>*</mo></msup><mrow><mo>(</mo><msub><mi>D</mi><mn>2</mn></msub><mo>(</mo><mrow><mi>c</mi><mo>-</mo><mi>v</mi></mrow><mo>)</mo><mo>)</mo></mrow><mo>&amp;rsqb;</mo><mo>}</mo><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>10</mn><mo>)</mo></mrow></mrow> In the formula, F and F ^-1 represent fast Fourier forward transform and inverse transform respectively, and F ^* represents the complex conjugate of F, where F is a two-dimensional or three-dimensional Fourier transform. 8. A robust multimodal remote sensing image matching system, characterized in that it comprises the following units: The preprocessing unit is used to judge the resolution information of the reference image and the input image. If the two images have the same resolution, execute the subsequent unit. If the resolutions are different, sample the two images to the same resolution; The template area selection unit is used to detect a series of evenly distributed feature points on the reference image by adopting the block strategy, which is denoted as P _1i (i=1,2,3,...,N), and the point P _1i is Select the template area AreaW _1i in the center; The matching area selection unit is used to predict the matching area AreaW _2i corresponding to the point set P _1i (i=1,2,3,...,N) on the input image according to the geographical coordinate information provided by the remote sensing image itself; A feature extraction unit is used to construct a three-dimensional dense feature expression map in the matching area; The initial matching unit is used to establish a similarity measure to match points with the same name on the basis of the three-dimensional dense feature expression map; for the obtained points with the same name, perform local extremum fitting on the similarity map to solve the sub-pixel position of the matching point ;Repeat the operation of the above unit, traverse each point of P _1i (i=1,2,3,...,N), and obtain the same-name point pair {PD _1i(x,y) with sub-pixel precision, PD _{2i(x ,y)} }(i=1,2,3,...,N); The matching point screening unit is used to eliminate the same-name point pairs with large errors in {PD _1i(x,y) , PD _2i(x,y) }(i=1,2,3,...,N) to obtain the final Point pairs with the same name {PID _1i(x,y) , PID _2i(x,y) } (i=1, 2, 3, . . . , S). 9. The multimodal remote sensing image matching system as claimed in claim 8, wherein the feature extraction unit is used to perform dense grid sampling on the image data in the matching area, and calculate the local features of each grid point Descriptors form a three-dimensional dense feature representation map. 10. The multimodal remote sensing image matching method according to claim 8, wherein the initial matching unit is based on a three-dimensional dense feature expression map, and establishes a similarity measure for image matching to obtain a similarity map, Take the extreme value position of the similarity map as the matching position of the image.