A Multi-Source Heterogeneous Large-Scene Remote Sensing Image Matching Method for Global Optimal Solution

Abstract

The present invention discloses a multi-source heterogeneous large-scene remote sensing image matching method for global optimal solution, including: S1, preprocessing the input image information; S2, extracting the reference image point group features by constructing a sub-pixel point group feature detection algorithm; S3, constructing a multi-mode matching algorithm based on wavelet analysis for sub-pixel matching; S4, constructing a point group global registration algorithm for forward and inverse projection, hierarchical elimination, and obtaining the global optimal solution; S5, judging whether the sub-pixel requirement is met. If it is judged that the requirement is not met, return to step S4, otherwise output the result. According to the present invention, a high-precision sub-pixel is obtained, the consistency of image intensity and phase is ensured, an accurate matching result is obtained, and an algorithm system with good compatibility is provided.

Description

Technical Field

[0001] The present invention relates to the technical field of remote sensing image matching, and particularly relates to a multi-source heterogeneous large-scale scene remote sensing image matching method with a globally optimal solution. Background Art

[0002] In recent years, high-resolution satellite image resources have become increasingly abundant. In particular, radar and interferometric radar images are attracting wide attention due to their high resolution, all-weather conditions, and powerful penetration capabilities. In terms of optical satellites, the spatial resolutions of satellites such as GF-1, GF-2, and GF-7 in China's high-resolution series can reach the sub-meter level. These satellites will surely be widely applied in many fields such as military, geoscience, oceanography, agriculture and forestry, resource exploration, disaster monitoring, global change, and national economic construction. Multi-source heterogeneous remote sensing image matching is one of the key supporting technologies for these applications. In this system, multi-source images refer to images from different platforms, different sensors, and different acquisition methods, including satellite optical images, satellite radar images, satellite interferometric radar images, aerial images, etc. Heterogeneous images refer to two images for matching with different imaging mechanisms and image structure characteristics, such as satellite optical images and radar images, images with the same orbit but different imaging times, radar images and interferometric radar images with different orbits and different imaging periods.

[0003] Multi-source heterogeneous remote sensing image matching mainly includes region-based and feature-based matching methods. The common challenges faced by various matching methods include complex geometric deformations, low image quality, and high computational and storage costs. Therefore, further research is needed on the accuracy, speed, reliability, and adaptability of matching. Research shows that there are the following main problems to be solved in multi-source heterogeneous remote sensing image matching: Corner extraction methods. In image matching, traditional methods are to select matching points in the reference image according to regular figures, or use matrix detection methods, use the statistical graph characteristics of the image to determine matching points, and use edges and line moments to describe image features. These methods are difficult to meet the requirements of multi-source remote sensing image matching. Until Morevec first proposed the concept of "interest points" in 1977, many new algorithms have emerged for feature detection around corners. Among them, the famous Harris corner detector uses the second-order moment or autocorrelation matrix to detect corners. There are also methods such as the SUSAN operator, FAST operator, Shi Tomasi operator, and feature extraction methods based on pixel comparison, also known as binary features. These methods are greatly affected by image scale, affine deformation, and radiation distortion. It is difficult to accurately select thresholds, with great uncertainty, slow solution speed, and it is difficult to achieve sub-pixel accuracy in detection. Especially in the matching of heterogeneous images, due to the large differences in imaging mechanisms, noise levels, and image structures between images, it is difficult to accurately extract corners on the reference image using the above methods. Therefore, it is necessary to study efficient and reliable sub-pixel corner extraction methods to meet the requirements of multi-source heterogeneous image matching.

[0004] 2) High-precision matching methods. To improve the accuracy of image matching, many application fields have been exploring effective ways of sub-pixel matching and have made great progress in recent decades, such as "Automatic registration of low-altitude high-resolution remote sensing images", "Remote sensing image registration based on feature points and mutual information", "Remote sensing image target recognition based on matching technology", "Stereo matching and 3D reconstruction of urban remote sensing image pairs", "Inter-frame image registration of high-orbit high-resolution remote sensing satellites", "Automatic registration of low-altitude high-resolution remote sensing images", "Automatic correction method of remote sensing images based on matching", "Sub-pixel fine matching method for low signal-to-noise ratio images", "Matching method of fuzzy logic coupled histogram classification", "Sub-pixel matching combining genetic algorithm and least squares method", etc. These methods are mainly used for the matching of homologous images. Some algorithms, such as those in OPENCV, also use the least squares method, but only perform least squares solution for redundant observations in region matching and cannot effectively correct image radiation distortion and geometric deformation, which is a very important link for non-homologous images and even heterogeneous images. How to take into account the characteristics of multi-source heterogeneous images, how to suppress noise, correct deformation, and achieve sub-pixel matching quickly and accurately is an urgent problem to be solved.

[0005] 3) Similarity measure. In image matching, common similarity measures include the Normalized Cross-Correlation (NCC), Sum of Squared Differences (SSD), Mutual Information (MI), etc. Most of these similarity measures are constructed based on the measurement of the gray information of the image and are used to measure the similarity between points. However, factors such as terrain undulation and image noise will have a significant impact on the final matching result. For example, if the existing algorithm uses a correlation measure (such as the correlation coefficient) as a single criterion, the matching result will have great uncertainty. Practice shows that for heterogeneous and isomeric images, a large correlation coefficient does not necessarily mean a correct match.

[0006] 4) Point-to-point matching method. As mentioned above, precise image matching is an important basis for a series of applications. Existing matching methods mainly perform point-to-point matching based on various measure functions and using local transformation relationships, without considering the global (overall) registration of the "point groups" between the two images. In the case of large terrain changes, especially when matching heterogeneous and isomeric images, due to the greatly different imaging mechanisms and noise levels, it leads to large deformations in the scale, distance, azimuth resolution, geometry, and radiation of the two images, resulting in "false matches". It is necessary to study a global matching method for the "point groups" of the two images to achieve the true correspondence of the "point groups" and images of the two images, conform to the state of the terrain surface, and meet the requirements of a series of subsequent applications. Summary of the Invention

[0007] Aiming at the deficiencies existing in the prior art, the object of the present invention is to provide a multi-source heterogeneous large-scale remote sensing image matching method with a global optimal solution, which obtains high-precision sub-pixels, ensures the consistency of image intensity and phase, and obtains accurate matching results, and has a very compatible algorithm system. To achieve the above object and other advantages according to the present invention, there is provided a multi-source heterogeneous large-scale remote sensing image matching method with a global optimal solution, including:

[0008] S1. Input image information for preprocessing;

[0009] S2. Extract the reference image point group features by constructing a sub-pixel point group feature detection algorithm;

[0010] S3. Construct a multi-mode matching algorithm based on wavelet analysis for sub-pixel matching;

[0011] S4. Construct a point group global registration algorithm for forward and inverse projection, hierarchical elimination, and obtain the global optimal solution;

[0012] S5. Judge whether the sub-pixel requirement is met. If the judgment does not meet the requirement, return to step S4; otherwise, output the result.

[0013] Preferably, in step S2, the sub-pixel point group feature detection algorithm uses the goodFeaturesToTrack function as the interface function of OpenCV, and performs sub-pixel matching by finding the perpendicularity of two vectors and then using the least squares method to determine the sub-pixel position of the intersection point.

[0014] Preferably, the sub-pixel point group feature detection algorithm includes the following steps:

[0015] S21. On the left and right images, perform automatic keypoint matching according to the image information to determine the relative position relationship between the two images and automatically determine the detection range;

[0016] S22. Extract the detection region of interest of the left image and take into account the deviation between the two images;

[0017] S23. Calculate the effective window size and determine the overall Manhattan distance and minimum point distance of the image;

[0018] S24. Extract the feature point groups by partition;

[0019] S25. Extract the sub-pixel feature point groups by least squares;

[0020] S26. Judge whether the extraction is completed. If it is judged that the extraction is not completed, return to step S25. Otherwise, take into account the relative position relationship between the two images, perform boundary detection, and output the feature point groups.

[0021] Preferably, in step S3, the multi-mode matching in the multi-mode matching algorithm based on wavelet analysis includes intensity matching, edge feature matching, and edge superposition intensity matching.

[0022] Preferably, the multi-mode matching algorithm based on wavelet analysis includes the following steps:

[0023] S31. Input the image feature point groups and perform matching mode selection;

[0024] S32. When selecting intensity and edge matching or edge feature matching, first perform edge feature extraction and then enter image decomposition; when selecting intensity matching, perform image decomposition;

[0025] S33. Perform wavelet decomposition at the top layer to establish a multi-scale pyramid image. The low-frequency image and high-frequency information of the pyramid image are stored in the same image memory;

[0026] S34. Use the low-frequency image at the top layer to perform initial matching according to the selected mode to obtain approximate positions for hierarchical matching. When performing hierarchical matching, automatically eliminate the worst matching points layer by layer;

[0027] S35. Determine whether the bottom layer has been reached. If it is determined that the bottom layer has not been reached, perform image reconstruction and return to step S34; otherwise, automatically adjust to the optimal matching window based on the result of hierarchical matching, perform least squares matching, use the left image as the reference image, perform radiometric correction and geometric correction on the right window image, and automatically converge.

[0028] Preferably, in step S4, the point cloud global registration algorithm obtains the global optimal solution through the projection model, taking into account the influence of the deformation of the two images, and performing inverse and forward projections; among them, the inverse projection has a positioning function, and the correct position of the right point of the image to be processed in the upper left of the reference image can be obtained; verification is performed through forward projection so that all points satisfy the characteristics of the ground surface corresponding to the two images.

[0029] Preferably, the point cloud global registration algorithm includes the following steps:

[0030] S41. Normalize the point cloud and correlation coefficient obtained by matching and form a weight matrix, which is used as a constraint condition for least squares global registration to form a projection model;

[0031] S42. Perform weighted LS global registration to obtain the inverse projection parameters;

[0032] S44. Calculate the coordinates of the point cloud using the parameters, compare them with the input point cloud coordinates, calculate the root mean square error mP, and determine a point with the largest error, record the point number Pn of this point and the maximum error value eMAX;

[0033] S45. Judge whether eMAX and mP are greater than sub-pixels. If so, discard this point, and use the new number of points N - 1 to perform weighted LS forward projection again. Iterate progressively in this way until mP is a sub-pixel, then the optimal registered point cloud can be obtained;

[0034] S46. Perform forward projection using the new point cloud to verify the reliability of the inverse projection.

[0035] Compared with the prior art, the beneficial effects of the present invention are as follows: The present invention proposes the concept of core points, uses core point matching to automatically determine the deviation between the two images, obtains the information of the overlapping area of the two images, performs hierarchical extraction, and introduces the Manhattan distance to automatically determine the number of extracted points and density to ensure obtaining sufficient feature points. On this basis, an algorithm for sub-pixel point cloud feature detection, Subpixel Points feature detection (SPFD), is constructed, and least squares positioning is used to ensure sub-pixel accuracy. At the same time, it can automatically detect the boundary of the point cloud to ensure the reliability of subsequent matching.

[0036] A multi-mode matching based on wavelet analysis is constructed. The directional section detection method is proposed, and a multi-mode matching system of intensity matching, intensity + edge feature matching, and edge feature matching is constructed to meet various types of image matching. The wavelet transform is used to decompose and reconstruct the image, and multi-scale and coarse-to-fine hierarchical matching is performed, which is beneficial to suppressing noise and improving speed; instead of setting the threshold parameter of the correlation coefficient, the matching result is determined through hierarchical optimization, which improves the automation level of matching; at the same time, the matching window can be automatically adjusted, and the image radiation deformation and geometric distortion are simultaneously performed on the window image, and the signal-to-noise ratio measure test is performed step by step to ensure obtaining the sub-pixel matching accuracy.

[0037] The method of point group and global registration MTGRP (Match and Global Registration For Point) is proposed, and a projection model is constructed. An important breakthrough is that the correlation coefficient is not used as the only criterion for judging effective matching, because in complex terrains, especially in heterogeneous image matching, the magnitude of the correlation coefficient is not sufficient to guarantee the reliability of matching. Instead, in global registration, the normalized correlation coefficient is used as a weight and used as a constraint condition to perform progressive least squares global registration on the initial matching point group, gradually eliminating the points that do not conform to the group until the root mean square difference between the two groups of point groups is less than 1 pixel and automatically converges. This algorithm provides a basis for realizing accurate image registration and mosaic fusion, especially provides a reliable basis for the application of heterogeneous images and radar interferometry, and is also applicable to homologous images.

[0038] The present invention has good robustness; it is suitable for the matching and registration of standard and non-standard aerial images, satellite optical images, radar images and interferometric radar images, infrared images and visible light images, and other heterogeneous and heteromorphic images. BRIEF DESCRIPTION OF THE DRAWINGS

[0039] Figure 1 It is a flowchart of the multi-source heterogeneous large-scale remote sensing image matching method for the global optimal solution according to the present invention;

[0040] Figure 2 It is a feature point position diagram of the multi-source heterogeneous large-scale remote sensing image matching method for the global optimal solution according to the present invention;

[0041] Figure 3 It is a sub-pixel point group feature detection algorithm diagram of the multi-source heterogeneous large-scale remote sensing image matching method for the global optimal solution according to the present invention;

[0042] Figure 4 It is a wavelet decomposition and reconstruction diagram of the multi-source heterogeneous large-scale remote sensing image matching method for the global optimal solution according to the present invention;

[0043] Figure 5Flow chart of edge feature detection of gradient direction for the multi-source heterogeneous large-scale remote sensing image matching method of the global optimal solution according to the present invention;

[0044] Figure 6 Multi-mode matching algorithm diagram based on wavelet analysis for the multi-source heterogeneous large-scale remote sensing image matching method of the global optimal solution according to the present invention;

[0045] Figure 7 Initial matching and global registration diagram of three-channel images for the multi-source heterogeneous large-scale remote sensing image matching method of the global optimal solution according to the present invention;

[0046] Figure 8 Initial matching of single-channel image and global registration diagram of single-channel image for the multi-source heterogeneous large-scale remote sensing image matching method of the global optimal solution according to the present invention;

[0047] Figure 9 Initial matching and global registration diagram of TerraSAR in-orbit spaceborne images for the multi-source heterogeneous large-scale remote sensing image matching method of the global optimal solution according to the present invention;

[0048] Figure 10 Initial matching and global registration diagram of high-resolution satellite and spaceborne TerraSAR feature images, high-resolution satellite and spaceborne SAR images for the multi-source heterogeneous large-scale remote sensing image matching method of the global optimal solution according to the present invention;

[0049] Figure 11 Edge superposition intensity diagram of resource satellite image and spaceborne TerraSAR, initial matching and global registration diagram of image for the multi-source heterogeneous large-scale remote sensing image matching method of the global optimal solution according to the present invention;

[0050] Figure 12 Initial matching and global registration diagram of visible light and infrared feature images for the multi-source heterogeneous large-scale remote sensing image matching method of the global optimal solution according to the present invention;

[0051] Figure 13 Initial matching and global matching diagram of satellite interferometric cross-orbit TerraSAR images for the multi-source heterogeneous large-scale remote sensing image matching method of the global optimal solution according to the present invention;

[0052] Figure 14 Initial matching and global registration diagram of in-orbit SAR images for the multi-source heterogeneous large-scale remote sensing image matching method of the global optimal solution according to the present invention;

[0053] Figure 15 Initial matching and global registration diagram of Tianhui satellite interferometric SAR images for the multi-source heterogeneous large-scale remote sensing image matching method of the global optimal solution according to the present invention;

[0054] Figure 16 Registration point and sampling point verification diagram in the related algorithm of the multi-source heterogeneous large-scale remote sensing image matching method for the global optimal solution according to the present invention. Specific implementation mode

[0055] The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts shall fall within the protection scope of the present invention.

[0056] Refer to Figure 1-16 , a multi-source heterogeneous large-scale remote sensing image matching method for the global optimal solution, including: S1. Input image information for preprocessing;

[0057] S2. Extract the feature of the reference image point group by constructing a sub-pixel point group feature detection algorithm;

[0058] S3. Construct a multi-mode matching algorithm based on wavelet analysis for sub-pixel matching;

[0059] S4. Construct a point group global registration algorithm for forward and inverse projection, and eliminate layer by layer to obtain the global optimal solution;

[0060] S5. Determine whether the sub-pixel requirement is met. If the judgment does not meet the requirement, return to step S4; otherwise, output the result.

[0061] The sub-pixel point group feature detection algorithm uses the goodFeaturesToTrack function as the interface function of OpenCV, and an improvement is made for sub-pixel matching, that is, by finding the perpendicularity of two vectors and then using the least squares method to determine the sub-pixel position of the intersection point.

[0062] As Figure 2 shown, q is the sub-pixel point to be obtained, p i is the adjacent point, and the vector G i is the gradient at p i , and is perpendicular to the vector (p i -q)

[0063] Then, an error equation is listed from all points in the window

[0064] v = G i *(p i -q) (2-1)

[0065] That is

[0066] G i *p i = G i*q (2-2)

[0067] The sub-pixel solution of the q-point position can be obtained by the least squares method

[0068]

[0069]

[0070] The solution is an iterative process, and the iterative termination condition can be automatically set until the q position reaches 0.001 pixels. The specific implementation steps are as Figure 3 .

[0071] Another improvement is to use the Manhattan distance instead of the Euclidean distance. The Euclidean distance represents the straight-line distance between two points in two n-dimensional vectors

[0072]

[0073] In fact, in the real world, it is often impossible to reach the target point in a straight line from the origin. Therefore, the present invention introduces a Manhattan distance:

[0074]

[0075] The Manhattan distance is also called the City Block Distance, and it has the same mathematical properties as the Euclidean distance, such as non-negativity, identity, and symmetry. Its geometric meaning is more suitable for image matching because the image reading and image matching window traversal movement in the matching process are not straight lines but in the Manhattan way, and it is computationally simple, fast, and more suitable for large-scale image partitioning processing.

[0076] The most important features of the wavelet function in the multi-mode matching algorithm of wavelet analysis are orthogonality, scalability, translation invariance, symmetry, etc. In the present invention, it is mainly used to construct a multi-scale image pyramid, aiming to provide a basis for hierarchical matching to improve the matching efficiency, achieve low-pass filtering, improve the reliability of matching, and save memory for the processing of large-scale images.

[0077] The decomposition formula of the image f can be obtained from the wavelet function

[0078]

[0079] Reconstruction formula

[0080]

[0081] In the above formula, H is a low-pass filter and G is a high-pass filter, which are their dual forms respectively; dj，h d j，v d j，a respectively represent the characteristics of the image in the horizontal, vertical, and diagonal directions; H and G respectively represent low-pass and high-pass filters. The superscript j represents the layer number of the image pyramid; the subscripts r and c represent the filtering operations along the row and column directions. Figure 4 is an example.

[0082] For the directional cross-section detection method of edge features, in the matching of heterogeneous and isomeric images, the images are very different, even completely inverted, and reliable results cannot be obtained using intensity (gray level) matching. Edge feature matching is an effective way to solve this problem. The present invention proposes a directional cross-section detection algorithm to extract edge features and construct multi-mode matching using edge features to handle different image matching and improve reliability and robustness.

[0083] The directional cross-section edge detection method is described as follows.

[0084] The gradient direction of the centered image f(c, r) is

[0085]

[0086] According to the ε minimum criterion

[0087]

[0088] In the matching window s, all pixels are used to form an error equation, and the coefficient a is solved by the least squares method, and the gradient and direction (argument) are calculated

[0089]

[0090] Equation (2-8) shows that in each matching window, the edge feature points always appear on the cross-section with intensity mutation. Along the gradient direction cross-section, the point (c, r) with the maximum gradient is the edge point. The process is as Figure 5 shown. Edge features can better express the similarity of two images with completely inverted intensities and can be used as the basis for matching for edge feature matching. At the same time, in order to increase the information content of the matching window, the matching of edge superposed intensity is also a good choice. It will be proved in the subsequent examples.

[0091] Correlation measure based on signal-to-noise ratio and least squares matching: The sum of the squares of the left window image gray levels is defined as the signal power

[0092]

[0093] The sum of the squares of the gray level differences between the left and right images is defined as the noise power

[0094]

[0095] The signal-to-noise ratio is thus obtained as

[0096]

[0097] The relevant measure based on the signal-to-noise ratio is thus expressed as

[0098]

[0099] The goal of least-squares sub-pixel matching is to maximize the similarity between the left and right window images, which is achieved through geometric distortion correction and radiometric distortion. An image consists of two parts: signal and noise. In the ideal state, the following relationship holds.

[0100] g l (x l ,y l ) + n l (x l ,y l ) = h0 + h l g r (x r ,y r ) + n r (x r ,y r )

[0101] Considering the scale and rotation between the left and right window images

[0102] x r =a0 + a1x l + a2y l

[0103] y r =b0 + b1x l + b2y l

[0104] The error equation can be obtained from the above two equations as

[0105] v = n l (x l ,y l ) - n r (x r ,y r ) = h0 + h1g r (x r ,y r ) - g1(x l ,y l ) (2 - 13)

[0106] The corrected error equation after linearizing the above equation is

[0107] v = c0dh0 + c1dh1 + c2da0 + c3da1 + c4da2 + c5db0 + c6db1 + c7db2 - ⊿g (2-14)

[0108] Among the eight unknowns in the above formula, dh is the radiation correction coefficient; da is the translation, scale, and rotation correction coefficient in the x direction; db is the translation, scale, and rotation correction coefficient in the y direction.

[0109] Matrix forms of the error equation, normal equation, and solution

[0110] V n*n,1 = A n*n,1 X 8,1 - L n*n,1

[0111]

[0112] By iteratively solving (2-14), image radiation and geometric deformation can be achieved, and precise sub-pixel matching can be accomplished.

[0113] The point cloud global registration algorithm uses the projection model, takes into account the influence of the deformation of the two images, performs forward and inverse projections, and obtains the global optimal solution. The inverse projection has a positioning function: obtaining the correct position of the right point in the image to be processed in the upper left of the reference image. And it is verified by the forward projection to make all points satisfy the characteristics of the ground surface corresponding to the two images. The processing process does not resample the images, keeping the coordinates of the matching point cloud of the two images and the geometric and physical characteristics of the images unchanged. This is particularly important for the subsequent processing of interferometric radar images, so that the information such as the phase, imaging angle, distance, and direction resolution of the original image based on the overall matching point cloud position remains unchanged for a series of applications such as DEM construction and interferogram generation. Moreover, only with point cloud registration can precise image correction be carried out instead of just "mosaic", and then image mosaicing in the same coordinate system can be realized.

[0114] To describe the distribution form of the point cloud in the image, a bi-binary 4th-order homogeneous equation is introduced, which describes a series of morphological changes such as the displacement, scale difference, rotation, hyperboloid, and paraboloid of the two images.

[0115] The relationship between the number of coefficients (unknowns) n of the homogeneous polynomial and the order k of the polynomial is as follows:

[0116]

[0117] For the 4th-order homogeneous, n = 15, and there are 30 unknowns for x and y in total.

[0118] When the number of matching point clouds is m, the point cloud error equation is

[0119]

[0120] In the formula, the coefficient matrix is

[0121]

[0122]

[0123] The parameter X to be solved = [a1 a2 … a n T (2-18)

[0124] The constant term L is the coordinates of the matching point pairs and serves as the observed value

[0125] L = [(x r -x l )1 (x r -x l )2 … (x r -x l ) m T (2-19)

[0126] As mentioned above, in images, especially in heterogeneous images, the image noise varies greatly. The correlation coefficient obtained in hierarchical matching cannot be used as a reliable measure to judge whether the matching is correct, and can only be used as a constraint condition in the projection model to construct the weight matrix and participate in the adjustment.

[0127] The normalized correlation coefficient of the matching point pairs is

[0128]

[0129] The weight matrix is

[0130]

[0131] The solution of the parameter to be solved can be written as

[0132]

[0133] According to Equation (2-21), weighted least squares global registration is performed. Usually, a system of equations with thousands of equations needs to be solved. Through progressive iteration, unqualified points are eliminated one by one to obtain the global optimal registration of the left and right image point groups. It should be noted that the convergence condition of MTGRP is that the root mean square error is less than the sub-pixel to ensure high-precision matching.

[0134] Matching of the first type of remote sensing image in Embodiment 1

[0135] Matching of three-channel aerial remote sensing images:

[0136] Table 4.1 Main parameters of three-channel aerial remote sensing images

[0137] ​​ Image Overlap degree Ground resolution m Image size (height * width) Type Left / Right 60％ 0.5 13824*7680 1

[0138] As can be seen from the above table, the three-channel aerial remote sensing image is a typical first-class homologous image and operates automatically throughout the process.

[0139] The program running process includes: number of wavelet layers = 4

[0140] Approximate deviation of the core point: offsety = 176 offsetx = -2908

[0141] Manhattan distance = 18596 point distance = 90

[0142] Total number of sub-pixel strong feature points: 3628 Sub-pixel strong feature points within the boundary: 3432

[0143] **********Lowest layer - Number of valid points matched in the first layer**********: 3432

[0144] Number of points for global registration (back-projection R-L) = 1122 Root mean square variance of point positions = 0.73974

[0145] Number of points for forward projection L-R (testing back-projection) = 1122 Difference between the root mean square of forward and back projections = 0.00087

[0146] Global deviation of the image offy = -0.0 offx = -2959.6

[0147] Running time of the program segment: 63.69s

[0148] It can be seen from the process that the aerial image has a standard overlap degree, and both the scale and rotation are within the specified requirements. Therefore, the number of points eliminated in the hierarchical matching is limited. After global matching (forward projection), the error is within sub-pixels. After back and forward projections, the root mean square difference is less than 0.74 pixels, fully demonstrating the reliability of global registration.

[0149] After core point matching, the approximate positions of the core points of the left and right images are obtained, which determines the basic relationship between the left and right images and provides a basis for the extraction of the effective point group, the estimation of the initial positions of the matching points, and even the hierarchical matching. After global registration of the point group, the overall deviation state of the image is determined.

[0150] From Figure 7 It can be seen that the number of initially matched points is 3432. From the initial matching results, the aerial image has a standard overlap degree and a small rotation. After hierarchical matching and step-by-step optimization, a good matching effect is obtained. However, under the conditions of buildings and undulating terrain, the influence of projection differences is relatively large. After global registration by forward and back projections, the number of selected points is 1122. These points conform to the terrain change state and are sufficient to form an effective control network, providing a basis for subsequent applications.

[0151] Matching of Channel Aerial Remote Sensing Images:

[0152] Table 4.2 Main Parameters of Single-Channel Aerial Remote Sensing Images

[0153] Image Overlap degree Ground resolution m Image size (height * width) Type Left / Right 60％ 0.5 5515*5554 / 5515*5518 1

[0154] Program Running Process

[0155]

[0156]

[0157] The running process shows that in the hierarchical matching process, almost no points are eliminated, which is the characteristic of homologous images. It should be noted that this is only point-to-point matching. In order to make the point clusters fit the terrain of the two images, global matching is carried out, and a large number of points will also be eliminated. And the root mean square differences of back-projection and forward-projection are both less than 0.71 pixels.

[0158] Finally, the number of globally matched points reaches 2673 pairs, which is better than the result of the three-channel. Figure 8 and Figure 9 , which are the results of initial matching and global matching respectively.

[0159] Matching of TerraSAR Constellation Spaceborne Images:

[0160] Table 4.3 Main Parameters of TerraSAR Constellation Spaceborne Images

[0161] Image Overlap degree Acquisition date Image size (height * width) Type Left / Right Basic overlap 2020.12.22 / 2020.11.30 2590*4980 1

[0162] The revisit period of TerraSAR is 11 days. The example images are images of the same orbit but different times. Although the images deviate, they basically coincide; the noise level is high, but at the same level. They are the first type of images and can complete the matching automatically through large-scale search without collecting any points.

[0163] Program Running Process

[0164]

[0165] It can be seen from the process that compared with aerial images, the imaging mechanism and noise level of TerraSAR images are completely different, and there are deviations and scale differences. Therefore, in hierarchical matching, more points are eliminated than in aerial images, and in global registration, the elimination rate is smaller. This is because SAR is range imaging and the perspective contraction and other deformations of constellation images have little influence. This characteristic is also manifested as that the root mean square differences of forward and back projections are very small. Such as Figure 8 .

[0166] Example 2 Matching of the Second Type of Remote Sensing Images

[0167] Matching of High-Resolution Satellite Images and Spaceborne TerraSAR Images:

[0168] Table 4.4 Main Parameters of High-Resolution Satellite Images and Spaceborne TerraSAR Images

[0169] Image Overlap degree Channel (left / right) Image size (height * width) Type Left / Right Basic overlap 3 / 1 2200*1800 (height * width) 2

[0170] In this example, the high-resolution satellite image has three channels and the SAR image has a single channel. Although there are characteristics of heterogeneous images, the images basically overlap, with small rotation and deviation, so they are classified into the second type of images. Different from the first type of images, only a smaller search area is needed. Due to the small size of the images, two-layer matching is sufficient.

[0171] Program Running Process

[0172]

[0173] For heterogeneous and isomeric images with different imaging mechanisms and significant differences in contour details, the noise of spaceborne optical images is much smaller than that of spaceborne SAR images, and the intensity maps of the two types of images are in an inverted state. Therefore, an edge feature + intensity superposition matching mode is adopted. Figure 9 Extract the edge superposition intensity map from the first image in

[0174] The global registration is very different from that of aerial images. Although the image format is not large, sufficient registration points are still obtained, and the root mean square errors of forward and inverse projections are both less than 0.94, proving the reliability of global registration.

[0175] Figure 9 The second and third images are the initial matching result and the global registration result respectively.

[0176] Matching of Resource Satellite Images and Spaceborne TerraSAR Images:

[0177] Table 4.5 Main Parameters of Resource Satellite Images and Spaceborne SAR Images

[0178] Image Overlap degree Channel (left / right) Image size (height * width) Type Left / Right Basic overlap 3 / 1 4400*2200 2

[0179] Program Running Process

[0180]

[0181] In Table 4.5, the left image is a resource satellite three-channel image, and the right image is a spaceborne single-channel SAR image. It is the same as the previous pair and belongs to the second type of images. The difference is that the noise levels of the two images are a little closer. For heterogeneous images with isomeric characteristics, an edge feature + intensity matching mode is adopted. Figure 10 The first one is the edge superposition intensity map.

[0182] From the process and Figure 10 As can be seen from the second and third images, since the noise between the resource satellite image and the SAR image is small; while the noise difference between the high-resolution image and the SAR image is relatively small, it is natural to obtain a larger group of registration points.

[0183] Visible and infrared image matching:

[0184] Table 4.6 Main parameters of resource satellite images and spaceborne SAR images

[0185] Image Overlap degree Channel Image size (height * width) Type Left / Right Basic overlap Visible light 1 / Infrared 1 512*512 2

[0186] Visible and infrared images. In this example, it has the smallest format, is a heterologous and heterogeneous image, the imaging is completely inverted, is divided into two layers, and edge feature matching is used.

[0187] Program running process

[0188]

[0189]

[0190] For two images in an inverted state, it is difficult to obtain good results using general intensity matching. The present invention uses a matching algorithm that combines edge features and intensity information, effectively overcoming the problem of image inversion, and the root mean square error of the forward and inverse projections of global registration is less than 0.68. Global precise registration is achieved.

[0191] Example 3 Matching of the third type of remote sensing image

[0192] The third type of remote sensing image is the most complex type in this example, which is basically a radar image. The main feature is that the acquisition time interval is long. Affected by factors such as the earth's rotation, the distance, azimuth resolution, scale factor, incident angle, etc. of the two images are very different, resulting in large perspective contractions and other deformations in the distance direction of the right image. In order to ensure image and phase consistency, the image cannot be resampled, and only a few corresponding points are roughly read on the image pair to achieve hierarchical matching and global registration from approximate to precise.

[0193] Satellite interferometric cross-track TerraSAR image matching:

[0194] Table 4.7 Main parameters of satellite interferometric cross-track TerraSAR images

[0195]

[0196] The so-called different orbits refer to ascending orbits and descending orbits. The descending orbit scans and forms images from north to south, and the ascending orbit scans and forms images from south to north. From the parameter analysis in Table 4.7, the distance and azimuth resolutions, scales, resolution, and incident angles of the two images are very different, which are typical heterogeneous images.

[0197] From Figure 12 It can be seen that there are significant differences in the scales and rotations of the two images, and the incident angle difference is nearly 23 degrees, causing great deformation. At the same time, in order to maintain the consistency between the image and the ground coordinates, the ascending orbit image is rotated by 180 degrees (i.e., changing the image reading order), without changing the correspondence between the image points and the interference phase, to ensure the requirements of user interference processing.

[0198] Only approximately 5 homologous points are collected as the basis for determining the initial position, and hierarchical matching from approximation to precision is realized.

[0199] Program running process

[0200]

[0201] Process and Figure 13 It shows that global registration is carried out on the basis of hierarchical matching. In order to adapt to terrain changes, direct back-projection is carried out, and a large number of eliminations are made. This is an inevitable result of the large differences in the deformations such as scales, rotations, and perspective contractions of the two images. Through gradual optimization, 312 pairs of matching point groups are still obtained, forming a dense control network. The root mean square errors of both forward and backward projections are less than 0.89 pixels, and the difference between the root mean square errors of the forward and backward projection point positions is also less than 0.11 pixels, ensuring the reliability of global registration and providing a basis for subsequent applications (such as image registration, mosaicking, DEM extraction, and new target point recognition).

[0202] Satellite interferometric along-track SAR image matching:

[0203] Table 4.8 Main parameters of satellite interferometric along-track SAR images

[0204]

[0205] Program running process

[0206]

[0207] The biggest feature of the example image is that the imaging times are 110 days apart, experiencing 10 imaging cycles. The distance and direction resolutions, distance and direction scale factors, especially the distance resolution and distance scale factor, are very different, forming an obvious contraction. In order to ensure that the results of matching and global registration correspond to the original image phase and other parameters unchanged, no image resampling is performed. Approximately 5 homologous points are collected to complete matching and registration, and 1650 sufficiently dense point pairs are obtained. The result parameters can be seen in the process.

[0208] Tianhui satellite interferometric SAR image matching:

[0209] Table 4.9 Main parameters of Tianhui satellite interferometric SAR images

[0210]

[0211] Program running process

[0212]

[0213] The Tianhui-2 satellite system launched by my country is my country's first microwave mapping satellite system based on interferometric synthetic aperture radar technology, and is also the second microwave interferometric mapping satellite system after the German TanDEM-X system. This example uses the Tianhui-2 01 satellite (launched on April 30, 2019) to study the matching of domestic interferometric SAR satellite images, which is of great significance in many future application fields. As can be seen from Table 4.9, the two images have great differences in range and azimuth resolution and incident angle, and are obviously heterogeneous images.

[0214] From the running process and Figure 15 It can be seen that the two images are greatly offset, with an overall deviation of nearly 10,000 pixels in the y direction. The effects of scale, resolution, rotation, etc. are clearly revealed. In order to correct these deformations, an accurate point group (738) is obtained in the overlapping area through global alignment, forming an accurate control point network. Both the forward and reverse projections meet the sub-pixel requirements, and the difference in the root mean square difference of the point positions is less than 0.00478.

[0215] In Table 4.10, the calculation time is no longer an important indicator as computer performance continues to improve, but it can reflect the influence of comprehensive factors, such as image size, complexity, number of points extracted, etc. In the present invention, the number of points extracted is automatically determined without restriction, and is strictly automatically screened based on the projection model. The maximum error and root mean square error must meet the sub-pixel requirements.

[0216] The above 9 images can all be checked from the global registration result map and can be verified by third-party software. In order to compare and ensure its authenticity, especially for the 6th image, Photoshop is used to collect 82 uniform points on the left projection image and compare the coordinates of the valid matching points to calculate the root mean square difference.

[0217] Table 4.10 Comparison of three algorithms

[0218]

[0219]

[0220] As can be seen from Table 4.10, the effective number of matching points of the present invention is related not only to the image size but also to the degree of image distortion; the root mean square error is related to factors such as image distortion and image deviation.

[0221] In the examples, the global registration results of the 9 images can all be manually verified in an image tool (such as Photoshop). To compare with the other two algorithms, 86 points were collected from the left image (back-projection result) in Figure 4 using Photoshop, as Figure 16 shown. The first figure is the left image of global registration, with green for the registration points and red for the 86 points collected; Figure 16 the second figure is the distribution map of the 86 points. By comparing the coordinates of the collected points with the coordinates of the global registration points, the root mean square error is 0.306.

[0222] All comparison results verify that:

[0223] 1. The present invention has good robustness for compatible multi-source heterogeneous images; and the function of multi-mode matching.

[0224] 2. The accuracy of sub-pixel matching and registration, with the maximum error not exceeding 0.9 pixels among the global registration points, and can obtain a control network with sufficient density.

[0225] 3. Global forward and back-projection verification, independent of the reliability of relevant measures.

[0226] For heterogeneous images such as radar, no resampling is performed, and the original image structure is maintained, which is beneficial for subsequent applications.

[0227] The number of devices and the processing scale described here are used to simplify the description of the present invention, and it is obvious to those skilled in the art for the application, modification and variation of the present invention.

[0228] Although the embodiments of the present invention have been disclosed above, it is not limited to the applications listed in the specification and the embodiments. It can be fully applied to various fields suitable for the present invention. For those familiar with the field, additional modifications can be easily achieved. Therefore, without departing from the general concept defined by the claims and the equivalent scope, the present invention is not limited to the specific details and the illustrated examples described here.

Claims

1. A multi-source heterogeneous large-scale remote sensing image matching method for global optimal solution, characterized in that It includes the following steps: S1. Input image information and perform preprocessing; S2. Extract the reference image point group features by constructing a sub-pixel point group feature detection algorithm; S3. Construct a multi-mode matching algorithm based on wavelet analysis for sub-pixel matching; The multi-mode matching in the multi-mode matching algorithm based on wavelet analysis includes intensity matching, edge feature matching, and edge superposition intensity matching; The multi-mode matching algorithm based on wavelet analysis includes the following steps: S31. Input the image feature point group and select the matching mode; S32. When selecting intensity and edge matching or edge feature matching, first extract the edge features and then enter the image decomposition; When selecting intensity matching, perform image decomposition; S33. Perform wavelet decomposition at the top layer to establish a multi-scale pyramid image. The low-frequency image and high-frequency information of the pyramid image are stored in the same image memory; S34. Use the low-frequency image at the top layer to perform an initial match according to the selected mode to obtain approximate point positions for hierarchical matching. When performing hierarchical matching, automatically eliminate the worst matching points layer by layer; S35. Determine whether the bottom layer has been reached. If it is determined that the bottom layer has not been reached, perform image reconstruction and return to step S34; Otherwise, automatically adjust to the best matching window through the results of hierarchical matching, perform least squares matching, use the left image as the reference image, perform radiometric correction and geometric correction on the right window image, and automatically converge; S4. Construct a point group global registration algorithm for forward and inverse projection, hierarchical elimination, and obtain the global optimal solution; S5. Determine whether the sub-pixel requirements are met. If it is determined that the requirements are not met, return to step S4; Otherwise, output the result.

2. The multi-source heterogeneous large-scale remote sensing image matching method for global optimal solution according to claim 1, wherein In step S2, the sub-pixel point group feature detection algorithm uses the goodFeaturesToTrack function as the interface function of OpenCV, and performs sub-pixel matching by finding the perpendicularity of two vectors and then using the least squares method to determine the sub-pixel position of the intersection point.

3. The multi-source heterogeneous large-scale remote sensing image matching method for global optimal solution according to claim 2, wherein, The sub-pixel point group feature detection algorithm includes the following steps: S21. On the left and right images, automatically match the core points according to the image information to determine the relative position relationship between the two images and automatically determine the detection range; S22. Extract the detection region of interest of the left image and consider the deviation between the two images; S23. Calculate the effective window size and determine the overall Manhattan distance and minimum point distance of the image; S24. Extract the feature point group by partition; S25. Extract the sub-pixel feature point group by least squares; S26. Determine whether the extraction is completed. If it is determined that the extraction is not completed, return to step S25; Otherwise, consider the relative position relationship between the two images, perform boundary detection, and output the feature point group.

4. The multi-source heterogeneous large-scale remote sensing image matching method for global optimal solution according to claim 1, characterized in that, In step S4, the point group global registration algorithm uses the projection model, considers the influence of the deformation of the two images, performs forward and inverse projection, and obtains the global optimal solution; Among them, the inverse projection has a positioning function and will obtain the correct position of the right point of the image to be processed at the upper left of the reference image; Verify through forward projection to make all points meet the characteristics of the ground surface corresponding to the two images.

5. The multi-source heterogeneous large-scale remote sensing image matching method for global optimal solution according to claim 4, characterized in that, The point group global registration algorithm includes the following steps: S41. Normalize the matched point cloud and correlation coefficient to form a weight matrix, which is used as a constraint condition for least squares global registration to form a projection model. S42. Perform weighted LS global registration to obtain back-projection parameters. S44. Calculate the coordinates of the point cloud using the parameters, compare them with the input point cloud coordinates, calculate the root mean square error mP, determine a point with the largest error, record the point number Pn of this point and the maximum error value eMAX. S45. Judge if eMAX and mP are greater than sub-pixels. If so, discard this point, and use the new number of points N-1 to perform weighted LS forward projection again. Iterate progressively like this until mP is a sub-pixel, then the optimal registered point cloud is obtained. S46. Use the new point cloud to perform forward projection to verify the reliability of the back-projection.

Images