Two-stage SAR image registration method based on SAR-SIFT and template matching

Abstract

The invention provides a two-stage SAR (Synthetic Aperture Radar) image registration method based on SAR-SIFT (Scale Invariant Feature Transform) and template matching, which comprises a pre-registration stage and a fine registration stage, and is characterized in that in the pre-registration stage, a down-sampling image is calculated through an SAR-SIFT algorithm to obtain an initial registration image; updating the initial registration image to obtain a final matching pair; after the final matching pair is mapped back to the original image size, a pre-registration image is obtained according to the final affine transformation matrix; in the fine registration stage, stable feature points are extracted from the pre-registration image blocks, and mismatching pairs are eliminated to obtain a fine registration image; and judging whether a region with relatively large topographic relief exists through the clustering mismatching pair, sampling a candidate topographic relief region, and then carrying out local correction to obtain a fine registration image again. The SAR-SIFT algorithm and iterative optimization are utilized to realize initial registration, a similarity measurement method is utilized to realize accurate registration of a flat region, adaptive stretching transformation is utilized to transform a topographic relief region to realize accurate registration of the topographic relief region, and finally accurate registration of the SAR image is realized.

Description

TECHNICAL FIELD

[0001] The present application belongs to the technical field of remote sensing image processing, and particularly relates to a two-stage SAR image registration method based on SAR-SIFT and template matching. BACKGROUND

[0002] SAR image registration is a key step in remote sensing image processing, and is widely used in change detection, image fusion, deformation monitoring and other tasks. Due to the inherent multiplicative noise, speckle interference and significant intensity difference of SAR images, there are often complex nonlinear intensity changes between SAR images acquired at different times or by different sensors. In addition, factors such as terrain undulations and sensor viewing angle differences will cause serious local geometric deformation, which makes the registration of SAR images, whether feature-based or template matching-based registration methods, problematic.

[0003] Region-based methods, also known as template matching, achieve geometric transformation estimation by optimizing the similarity measure between two images. Document one (A. Mahmood and S. Khan, "Correlation-Coefficient-Based Fast Template Matching Through Partial Elimination," IEEE Transactions on Image Processing, vol. 21, no. 4, pp. 2099-2108, April 2012.) points out that in the spatial domain, the traditional similarity measure is the normalized cross-correlation (NCC); document two (P. Kovesi, "Image features from phase congruency," Videre: Journal of Computer Vision Research, vol. 1, no. 3, pp. 1-26, 1999.) points out that in the frequency domain, the most commonly used similarity measure is phase congruency. In modern remote sensing sensors, using navigation devices such as global positioning systems, pre-registration can be achieved through direct georeferencing. Therefore, through this pre-registration, the obvious global geometric transformation (such as rotation and scale change) can be eliminated, leaving only a few tens of pixels of offset in the image. However, without pre-registration, these methods need to search the entire image, which usually requires high computational time. In addition, intensity-based methods perform poorly on images affected by strong noise and with significant geometric differences.

[0004] Feature-based methods find correspondences by matching local features extracted between images. The extracted features should be highly distinctive, stable and repeatable between images. These features include points, edges and regions, among which the study of points is the most, and the scale-invariant feature transform (SIFT) in document three (D. G. Lowe, "Distinctive image features from scale-invariant keypoints," International Journal of Computer Vision, vol. 60, no. 2, pp. 91-110, 2004.) is a classic feature point-based matching method. However, due to the existence of speckle noise, the SIFT method is unstable in feature extraction, and the Gaussian scale space of SIFT blurs a large amount of image information, which does not perform well in SAR images. Many scholars replace the Gaussian filter with several anisotropic filters, including bilateral filter SIFT (BF-SIFT) in document four (S. Wang, H. You, and K. Fu, "BF SIFT: A novel method to find feature matches for SAR image registration," IEEE Geoscience and Remote Sensing Letters, vol. 9, no. 4, pp. 649-653, 2011.) and the like to improve the matching performance. However, in these methods, the gradient is calculated by a simple difference equation, which is not suitable for SAR image matching with multiplicative speckle noise. In document five (F. Dellinger, J. Delon, Y. Gousseau, et al., "SAR-SIFT: A SIFT-like algorithm for SAR images," IEEE Transactions on Geoscience and Remote Sensing, vol. 53, no. 1, pp. 453-466, Jan. 2015.), Dellinger et al. proposed a SAR-SIFT algorithm, which uses a range of exponential weighted average (ROEWA) operator instead of a difference algorithm to calculate the gradient, and improves the gradient calculation. In finding feature points, the Harris corner and Gaussian filter are combined, and the scale space extreme point is used as the feature point. This method is more robust to multiplicative noise in SAR images. However, this image matching algorithm is a global image matching operation, and when the image content has a high degree of similarity, a large number of false matches will occur.

[0005] In recent studies, many scholars combine region-based and feature-based methods and propose a new matching method including two stages of pre-matching and fine-matching. In Paper Six (Y. Xiang, F. Wang, and H. You, "An automatic and novel SAR image registration algorithm: A case study of the Chinese GF-3 satellite," Sensors, vol. 18, no. 2, p. 672, Feb. 2018.), Xiang et al. introduce an adaptive sampling method in the coarse matching stage. Combined with the SAR-SIFT algorithm, it realizes the effective coarse registration of two large SAR images and eliminates the obvious translation, rotation and scale differences. In the fine matching stage, the fast normalized cross-correlation of the improved phase consistency model is used as a new similarity measure to find the corresponding points in the local search area. In Paper Seven (S. Paul and U. C. Pati, "SAR image registration using an improved SAR-SIFT algorithm and Delaunay-triangulation-based local matching," IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, vol. 12, no. 8, pp. 2958-2966, Aug. 2019.), Paul and Pati propose an improved SAR-SIFT (I-SAR-SIFT) method that calculates the general weight of each feature by combining SAR-Harris response, entropy and distribution value, and constructs a multi-scale description factor. On the basis of coarse matching, the image is divided into blocks and the corresponding areas are matched respectively. In Paper Eight (M. Wang, J. Zhang, K. Deng, and F. Hua, "Combining Optimized SAR-SIFT Features and RD Model for Multisource SAR Image Registration," IEEE Transactions on Geoscience and Remote Sensing, vol. 60, pp. 1-16, 2022.), Wang et al. propose a local matching method assisted by the RD model, which converts the global matching mode to the local matching mode, avoiding rough errors and improving the matching efficiency.The method considers the influence of the terrain undulating area, combines the optimized SAR-SIFT feature with the imaging geometry of the SAR system, but needs a series of imaging parameters of the remote sensing satellite, including the elevation of the target, the slant range of the sensor and the target, and the like, and has great limitations.

[0006] SAR image registration is a key step in remote sensing data processing, and directly affects the accuracy of change detection, target recognition and multi-source information fusion. Due to the influence of sensor viewing angle, terrain undulation and time phase change, there are significant geometric changes and intensity differences between different SAR images. As a method depending on suitable feature points, the global transformation obtained by SAR-SIFT often results in low registration accuracy. Template matching can obtain relatively accurate registration results, but under the condition of no pre-registration, it needs to search from the global, and is time-consuming. However, the SAR-SIFT method is robust to initial alignment, and the template matching has accuracy. SUMMARY

[0007] The present application aims to solve the problems of the prior art that different SAR images are difficult to be registered or have low registration accuracy based on feature and template matching methods, and provides a two-stage SAR image registration method based on SAR-SIFT and template matching. The SAR-SIFT algorithm and iterative optimization are used to realize initial registration and avoid large registration errors. The phase consistency in the frequency domain and the gradient information similarity measurement method in the spatial domain are combined to realize template matching and achieve accurate registration of flat areas. The adaptive stretching transformation is used to transform the terrain undulating area to realize accurate registration of the terrain undulating area. Finally, the SAR images are accurately registered under the interference of different time phase speckles and terrain undulation.

[0008] To achieve the above purpose, the technical solution provided by the present application is:

[0009] A two-stage SAR image registration method based on SAR-SIFT and template matching includes a pre-registration stage and a fine registration stage, wherein:

[0010] In the pre-registration stage, the initial registration image is obtained by the SAR-SIFT algorithm for the down-sampled image. The affine transformation matrix of the initial registration image is iteratively updated by using the iterative optimization matching method, and the final matching pair is obtained. The final affine transformation matrix is obtained by mapping back to the original image size, and the pre-registration image is obtained according to the final affine transformation matrix.

[0011] In the fine registration stage, stable feature points are extracted from the pre-registered image block and combined with the similarity measure to eliminate the false matching pairs in the reference image template to obtain a fine registration image; whether there is a large terrain undulating area is judged by clustering the false matching pairs, and if there is a large terrain undulating area in the fine registration image, the candidate terrain undulating area in the fine registration image is then sampled and subjected to local correction of the undulating area to obtain a fine registration image again.

[0012] As a further limitation of the application, the method of the pre-registration stage specifically comprises: after the input SAR image is down-sampled, feature point sets are obtained based on Harris corner points in a spatial pyramid through Gaussian blur, feature description vectors are extracted by calculating gradient information in different directions for each feature point, and initial registration images are obtained through nearest neighbor distance ratio and fast sampling consistency; an iterative optimization matching method is used to iteratively optimize and update the affine transformation model of the initial registration image until the model tends to be stable, and corresponding matching pairs are obtained; the matching pairs are mapped back to the original image size to obtain a final affine transformation matrix; and the pre-registered image is obtained according to the final affine transformation matrix.

[0013] As a further limitation of the application, in the pre-registration stage, the initial registration image is obtained from the down-sampled image through the SAR-SIFT algorithm, comprising the following steps:

[0014] Feature point detection: the input SAR image is down-sampled to obtain a down-sampled image, multi-scale space is constructed, the down-sampled image is smoothed through Gaussian kernel functions of different scales, and a Harris scale space pyramid is established; the gradient amplitudes of pixel points in the multi-scale space are calculated based on local window regions, each pixel point is compared in each layer of the constructed Harris scale space pyramid, and candidate feature points are selected and screened out, and the candidate feature points are accurately positioned and false low contrast points are removed to obtain a stable feature point set;

[0015] Main direction assignment and descriptor extraction:

[0016] For each detected feature point, a gradient direction histogram is constructed and the peak maximum direction is selected as the main direction of the feature point; a circular neighborhood is constructed around the feature point, and the gradient information in each direction is counted and combined to form a feature description vector in the respective region;

[0017] Feature point matching:

[0018] For each feature point in the image to be matched, two feature points with the nearest neighbor distance and the second nearest neighbor distance in the reference image are found, and the ratio of the two is calculated, and if the ratio is allowed, the matching pair is accepted to improve the matching accuracy, an affine transformation model is calculated, and the optimal solution of the model with the largest number of corresponding inliers is selected by iterative optimization to perform feature point matching.

[0019] As a further limitation of the present application, in the pre-registration stage, an affine transformation model of the affine transformation, whose formula is:

[0020]

[0021] In formula (1), a, b, c, d, e, and f are respectively affine transformation coefficients, (x1, y1) and (x2, y2) are respectively coordinates before and after transformation, and it is simply considered that the closer the coordinates (x1, y1) are to (x2, y2), the more stable the transformation is, otherwise the matching pair obtained is not stable, at this time the transformation matrix tends to be a unit matrix, and the iteration is controlled by limiting the affine transformation coefficients, and when the transformation coefficients are set as |a-1|<τ1&|b|<τ2&|c|<τ3 or greater than the maximum iteration number τ max , the iteration is stopped, and the coefficients d, e, and f have the same limiting conditions as b, a, and c respectively.

[0022] As a further limitation of the present application, in the fine registration stage, the pre-registration image is matched in a local area to obtain a fine registration image, comprising the following steps:

[0023] Feature point extraction: in the fine registration stage, the pre-registration image is divided into blocks, a Harris corner-based method is used to extract a set of feature points, and stable feature points are extracted by setting a threshold value;

[0024] Combination of phase consistency and gradient similarity measurement: based on phase consistency and gradient information, the stable feature points extracted by block are determined to have a similarity measurement, the stable feature points with high similarity value are removed as matching points in the reference image template to eliminate false matching pairs and obtain an initial fine registration image;

[0025] Local correction of terrain undulating area: the matching pairs that do not satisfy the final transformation model are clustered to obtain corresponding clustering matching pairs of different clusters, which are used as candidate terrain undulating areas, the boundary and the surrounding of the candidate terrain undulating areas are sampled and accurately calculated to further obtain more accurate local undulating areas, then each local undulating area is sampled and template matching is performed to obtain a final stretching transformation, and the local undulating areas are locally corrected to obtain a fine registration image.

[0026] As a further limitation of the present application, the similarity measurement combined with phase consistency and gradient information is calculated according to the following formula:

[0027] PC_G(x,y)=∑ o (PC 2o (x,y)+G o (c,y)) Formula (2)

[0028] In formula (2), PC 2o (x,y) and G o(x, y) represents the edge information obtained by using the phase congruency in frequency domain and the gradient in spatial domain respectively in the direction o; wherein:

[0029] The definition of the phase congruency model in frequency domain is:

[0030]

[0031] In formula (3), A n (x, y) represents the amplitude at the scale n, e n (x, y) represents the convolution of the image and the even-symmetry filter, n represents the number of scales, and I represents the image intensity, represents the even-symmetry filter at the scale n, o n (x, y) represents the convolution of the image and the odd-symmetry filter, represents the odd-symmetry filter at the scale n, represents the corresponding phase, W o represents the weighting function at the oth orientation, o represents the number of orientations, A no represents the amplitude at the scale n and the orientation o, T o represents the noise level at the oth orientation, and ε represents a small positive value to prevent the expression from becoming unstable when A no becomes very small, ΔΦ no (x, y) represents:

[0032] For PC2 (x, y), for a specific direction, it means the phase congruency in the direction, and the edge in the direction is obtained; according to this idea, the meaning of formula (4) is that the edge in the specific direction is taken after a certain value is taken; if 6 directions are taken together, the value of o is 0, 1, 2,..., 5, the orientation interval is 180° / 6 = 30°, and the direction is o*30°;

[0033]

[0034] The gradient information in the spatial domain is:

[0035]

[0036] In formula (5), I represents the image intensity, and respectively represent the local mean on both sides of the point (x, y), and the weight is distributed by using the Gaussian weighting method, wherein α represents the standard deviation of the Gaussian kernel. G o (x, y) represents the ratio gradient obtained in the direction o.

[0037] As a further limitation of the application, the matching pairs that do not satisfy the final transformation model are clustered to obtain different cluster corresponding matching pairs, which are used as candidate terrain relief areas; each relief area performs local adaptive stretching transformation to correct it locally, which specifically includes:

[0038] The rectangular boundary box of each terrain relief area drawn by the transformation range is calculated based on the distribution characteristics of the re-acquired matching pairs, and the boundary box is taken as the starting point to sample to both sides, if it is located in a flat area, it is sampled to the relief area, if it is located in a relief area, it is sampled to the flat area, and the accurate terrain relief area boundary is obtained through the spatial distribution characteristics of the sampling points; whether the adjacent sampling points belong to an independent area or the same continuous area is judged based on the distance measurement between the boundary sampling points, and it is further finely divided into multiple specific sub-areas;

[0039] The matching pairs obtained by uniform sampling are preliminarily roughly deleted using the affine transformation model, and the stretching transformation ratio of each finely divided sub-area is calculated from the correct matching pair with the maximum transformation distance in the area to ensure that the stretching transformation fully adapts to the local change characteristics of the terrain relief correction; the local adaptive stretching transformation is implemented on each specific terrain relief area according to the calculated stretching transformation ratio.

[0040] As a further limitation of the application, the transformation formula of the direction implemented in the above local adaptive stretching transformation is:

[0041]

[0042] In formula (6), x left represents the left boundary of the terrain relief area in the x direction, x right represents the right boundary of the terrain relief area in the x direction, x1 represents the x coordinate of the transformation pair determined by the pre-registration image, x2 represents the x coordinate of the transformation pair determined by the reference image, s left and s right represent the stretching ratio in the x direction, and x' represents the coordinate after stretching transformation at position x.

[0043] The advantages of the application are:

[0044] 1. The application uses SAR-SIFT algorithm and iterative optimization to realize initial registration and avoid large registration errors, uses similarity measurement obtained by combining phase consistency in frequency domain and gradient information in spatial domain to realize template matching for accurate registration of flat areas, and uses adaptive stretching transformation to transform terrain relief areas to realize accurate registration of terrain relief areas, thereby realizing accurate registration of SAR images under different temporal phase speckle interference and terrain relief.

[0045] 2、The application combines SAR-SIFT and iterative optimization to obtain a pre-registration image for a down-sampled image in a pre-registration stage, thereby reducing the amount of calculation while maintaining high registration accuracy and avoiding large initial registration errors; in the template matching method in the fine registration stage, a new similarity measurement method is studied by combining the phase consistency in the frequency domain and the gradient in the spatial domain, thereby improving the reliability of the similarity measurement; the image obtained by using the template matching method is further subjected to adaptive local stretching transformation in the fine registration stage, and the terrain undulating area is searched and corrected, thereby effectively implementing the registration of the area.

[0046] Additional aspects and advantages of the application will be set forth in part in the description which follows, and in part will become apparent to those skilled in the art upon examination of the following and / or can be learned by practice of the application. BRIEF DESCRIPTION OF DRAWINGS

[0047] The above and / or additional aspects and advantages of the application will become apparent and be readily appreciated from the following description, taken in conjunction with the accompanying drawings, in which:

[0048] Figure 1 The application provides a two-stage SAR image registration method based on SAR-SIFT and template matching;

[0049] Figure 2 The application provides a flowchart of the two-stage SAR image registration method in actual application;

[0050] Figure 3 The application provides a detection result diagram of the registration method shown in the SAR image without terrain undulating area; Figure 2

[0051] Figure 4 The application provides a detection result diagram of the registration method shown in the SAR image with terrain undulating area; Figure 2

[0052] Figure 5 The application provides a comparison table of SIFT, BF-SIFT, SAR-SIFT, I-SAR-SIFT and the method;

[0053] Figure 6 The application provides a terrain undulating stretching transformation diagram. DETAILED DESCRIPTION

[0054] The embodiments of the application are described in detail below, which are exemplary and intended to explain the application, and cannot be understood as a limitation of the application.

[0055] Reference is made to Figure 1 and Figure 2 ​​The embodiment of the present application provides a two-stage SAR image registration method based on SAR-SIFT and template matching, which comprises a pre-registration stage and a fine registration stage.

[0056] In the pre-registration stage, continuing to refer to Figure 1 and Figure 2 The embodiment of the present application processes the down-sampled image, obtains a pre-registration image through a SAR-SIFT algorithm, iteratively optimizes and updates an affine transformation model of the initial registration image until the model tends to be stable, and obtains a final matching pair; maps the final matching pair back to the original image size to obtain a final affine transformation matrix, and obtains a pre-registration image according to the final affine transformation matrix; specifically, the pre-registration image obtained by the down-sampled processing image through the SAR-SIFT algorithm replaces the current to-be-matched image to participate in the next round of matching, the iteration process is terminated when the pre-registration image gradually tends to be stable by setting a reasonable stopping rule, the registration error is reduced, the to-be-matched image and the reference image are roughly aligned, and the final matching pair of the down-sampled image is obtained.

[0057] Continuing to refer to Figure 2The feature point matching step in the embodiment of the application comprises the following steps: (1) feature point detection: firstly, a multi-scale space is constructed for an input SAR image, a Gaussian filter is combined with a Harris corner point detection to smooth the image by using a Gaussian kernel function of different scales, and a Harris scale space pyramid is established; an exponential weighted mean ratio operator is used to calculate the gradient amplitude and direction of a pixel point based on a local window region to suppress the multiplicative noise interference of the SAR image, then the pixel points in each layer of the established scale space pyramid are compared and candidate feature points are screened out, the response value of each pixel point and the adjacent points in the spatial neighborhood of the pixel point is calculated, the extreme points with significant are screened out as candidate feature points, accurate positioning and low-contrast point elimination are performed, and finally a stable feature point set is obtained; (2) main direction assignment and descriptor extraction: a gradient direction histogram is constructed for each feature point in the detected feature point set as the main direction of the feature point, a circular neighborhood is constructed around the feature point, and the gradient information of each direction is counted in the respective region to form a feature description vector; the gradient direction of the pixel within a certain radius range is calculated on the Gaussian smoothing image where each feature point in the detected feature point set is located, the gradient direction of 0-360 degrees is divided into 36 intervals by using a histogram counting method, a gradient direction histogram is constructed, and the direction with the maximum peak value in the histogram is selected as the main direction of the feature point; when there is a secondary peak value greater than 80% of the main peak value, the feature point is copied and a new main direction is assigned to ensure the rotation invariance, a circular neighborhood is constructed around the feature point in the main direction coordinate system, the neighborhood is divided into 9 sub-regions, the gradient information of 12 directions in each sub-region is counted, and finally a 108-dimensional feature description vector is formed to realize descriptor extraction; (3) feature point matching step: a nearest neighbor distance ratio matching strategy is adopted, the two feature points with the nearest and second nearest Euclidean distance are found in the target image for each feature point in the image to be matched, the ratio of the nearest neighbor distance to the second nearest neighbor distance is calculated, and when the ratio is less than a set threshold, the matching pair is accepted; a fast sampling consistency algorithm is used to screen the initial matching result to improve the matching accuracy, an affine transformation model is calculated by using a random sampling minimum point set, the number of inliers satisfying the model is counted, and after iterative optimization, the optimal solution of the transformation model with the largest number of inliers is selected as the feature point matching data.

[0058] The affine transformation model in the feature point matching step in the embodiment of the application has the formula:

[0059]

[0060] In the formula, a, b, c, d, e, f are affine transformation coefficients respectively, (x1, y1) and (x2, y2) are coordinates before and after transformation respectively, it is considered that the closer the coordinates (x2, y1) are to (x2, y2), the more stable the transformation is, otherwise the matching pair obtained is not stable, the transformation matrix tends to be a unit matrix, the iteration is controlled by limiting the affine transformation coefficients, when the transformation coefficients are set as |a-1|<τ1&|b|<τ2&|c|<τ3, or greater than the maximum iteration number τ max , the iteration is stopped, the coefficients d, e, f have the same limiting conditions as b, a, c respectively.

[0061] In the fine registration stage, continuing to refer to Figure 1 and Figure 2 , the embodiment of the present application performs matching on the pre-registration image in a local area, firstly extracts stable feature points by blocking the pre-registration image, then completes template matching in the reference image by combining phase consistency and gradient similarity measurement, removes false matching pairs by fast sampling consistency, and generates a final transformation model; whether there is a large terrain fluctuation area is judged by clustering false matching pairs, if there is a large terrain fluctuation area in the image, then the matching pairs that do not satisfy the final transformation model are clustered by adaptive local stretching transformation to obtain a candidate terrain fluctuation area; the boundary and the surrounding thereof are sampled to further obtain a more accurate fluctuation area, then the final stretching transformation ratio is obtained by sampling and template matching for each fluctuation area, and local correction of the fluctuation area is completed.

[0062] The embodiment of the present application uses PCGCC as a similarity measurement method to perform template matching of the pre-registration image and the reference image, and removes false matching points by fast sampling consistency and affine transformation to obtain a transformation model and a registration result in this stage.

[0063] Continuing to refer to Figure 2The embodiment of the application matches the pre-registration image in a local area in a fine registration stage, and comprises the following steps: (1) feature point extraction: firstly, the image is divided into 5*5 non-overlapping blocks, and the Harris value in each block is calculated; then, the Harris value in each block is sorted from high to low, and the k points with the highest response value are selected as the feature points; (2) similarity measurement combining phase consistency and gradient: the stable feature points extracted by block are determined based on phase consistency and gradient information to determine the similarity measurement, and the stable feature points with high similarity value are removed as the matching points in the reference image template to eliminate the false matching pairs, so that the initial fine registration image is obtained by fine registration; (3) local correction of terrain undulating area: after global correction transformation, the registration image and the reference image have accurate registration in the flat area, and the correct matching pairs in the terrain undulating area do not satisfy the final global transformation model; firstly, the matching pairs that do not conform to the final transformation model are found, the DBSCAN algorithm based on density is used to cluster these matching pairs, different clusters are obtained, and are regarded as candidate terrain undulating areas, and then the local adaptive stretching transformation is performed.

[0064] The new similarity measurement obtained by combining the phase consistency and gradient information in the embodiment of the application is defined as follows:

[0065]

[0066] In the formula, A n (x,y) represents the amplitude at the scale n, e n (x,y) represents the convolution of the image and the even-symmetry filter, n represents the number of scales, and I represents the image intensity, represents the even-symmetry filter at the scale n, o n (x,y) represents the convolution of the image and the odd-symmetry filter, represents the odd-symmetry filter at the scale n, represents the corresponding phase, W o represents the weighting function at the oth orientation, o represents the number of orientations, A no represents the amplitude at the scale n and the orientation o, T o represents the noise level at the oth orientation, and epsilon represents a small positive value to prevent the expression from becoming unstable when A no becomes very small, and Delta Phi no (x,y) is represented as:

[0067] For PC2(x, y), for formula specific direction, its meaning is to represent the phase consistency in this direction, and get the edge in this direction. According to this idea, the meaning of the following formula is to get the edge in the specific direction after taking a certain value. For example, 6 directions are taken, and the value of o is 0, 1, 2,..., 5, and the orientation interval is 180° / 6=30°, and the direction is o*30°. It is represented as:

[0068] The above edge of different directions, the calculation of the spatial domain gradient, is calculated according to the formula of the embodiment of the application:

[0069]

[0070]

[0071] In the formula, And Indicates the local mean value on both sides of the point (x, y) in the direction o, and the weight is distributed by using the Gaussian weighting method, wherein α represents the standard deviation of the Gaussian kernel; G o (x, y) represents the ratio gradient obtained in the direction o;

[0072] Similar to the idea of the exponential weighted mean ratio, the embodiment of the application uses an exponential weighted mean ratio (ROEWA) operator instead of the difference algorithm to calculate the gradient, and the gradient calculation is improved. The embodiment of the application combines the frequency domain and the spatial domain to obtain new edge information PC_G , which is represented as:

[0073]

[0074] The NCC of the PC_G described above in the embodiment of the application is used as a similarity measure for image registration, and is represented as PCGCC, which is defined as:

[0075]

[0076] Wherein, E and T represent the PC_G response of the image and the template, Indicates the average value of T, Indicates the average value of E which is the same size as the template, and (x, y) represents the coordinates of the feature point in the pre-registration image I.

[0077] The embodiment of the application calculates the similarity of all pixel points in the local area of the corresponding point on the template image for the feature point in the pre-registration image, and forms a candidate matching pair with the point having the highest similarity.

[0078] Further, referring to Figure 2The execution of the above local adaptive stretch transformation of the embodiment of the application includes the following steps: (1) determination of the rectangular frame of the candidate terrain undulating area: the candidate terrain undulating area obtained through DBSCAN clustering is re-sampled and template matching is performed to re-obtain matching pairs; based on the distribution characteristics of the matching pairs, the transformation range is calculated in the x direction and the y direction respectively, and the rectangular boundary frame of each candidate area is drawn; (2) sampling and detection of the boundary of the undulating area: the center point of the rectangular boundary frame is taken as the starting point in the x direction and the y direction respectively, and sampling is performed to the two sides to detect the boundary position of the terrain undulating area; if each sampling point is located in a flat area, sampling is performed to the undulating area, and if it is located in the undulating area, sampling is performed to the flat area; the boundary area of the terrain undulation is preliminarily identified through the spatial distribution characteristics of the sampling points, and the boundary of the terrain undulating area is obtained more accurately; (3) fine division of the terrain undulating area: based on the distance measurement between the boundary sampling points, whether the adjacent sampling points belong to the same terrain undulating area is judged by using an adaptive threshold; if the distance exceeds the preset threshold, the terrain undulating area is divided into an independent area; otherwise, it is regarded as the same continuous area, and the candidate terrain undulating area is further refined into a plurality of specific sub-areas; (4) calculation of the local stretch transformation ratio: based on the local transformation characteristics of the matching pairs, the stretch transformation ratio of each subdivided terrain undulating area is calculated; the matching pairs are obtained by uniform sampling, and then the affine transformation model is used to preliminarily roughly remove the false matching pairs; the stretch transformation ratio is determined by the correct matching pairs with the maximum transformation distance in the area, which ensures that the stretch transformation can fully adapt to the local change characteristics of the terrain undulation and improve the correction accuracy; (5) application of the local adaptive stretch transformation: according to the calculated stretch transformation ratio, the local adaptive stretch transformation is implemented for each specific terrain undulating area.

[0079] In the implementation of the local adaptive stretch transformation, the transformation formula of the direction is represented as:

[0080]

[0081] wherein, x left represents the left boundary of the terrain undulating area in the x direction, x right represents the right boundary of the terrain undulating area in the x direction, x1 represents the x coordinate of the transformation pair determined by the pre-registration image, x2 represents the x coordinate of the transformation pair determined by the reference image, s left and s right represent the determined x direction transformation ratio, and x' represents the coordinate after the stretch transformation at position x.

[0082] In the pre-registration stage, the embodiment of the present invention combines the SAR-SIFT algorithm with iterative optimization to obtain a pre-registered image, eliminating significant translation, rotation, and scale differences between the pre-registered image and the reference image. In the fine registration stage, a template matching method is used to first segment the image into blocks to ensure a uniform distribution of feature points in the image. Then, a new similarity metric is designed, combining phase consistency in the frequency domain with gradient information in the spatial domain, to determine matching points within the local search area. Finally, the DBSCAN algorithm is used to cluster point pairs that do not satisfy the final transformation matrix to obtain candidate terrain relief areas. These areas are corrected using the proposed adaptive local stretching transformation to complete the registration.

[0083] See Figure 6 , an example diagram of terrain relief stretching transformation according to an embodiment of the present invention, Figure 6 The following is a simple local transformation plan view caused by terrain undulations, where points a1 and b1, and a5 and b5 are corresponding boundary areas, and a2, a3, and a4 are the corresponding correct matching pairs with b2, b3, and b4, respectively. At this time, the transformation distances of these three correct matching pairs are 1.5, 3, and 2, respectively. The stretching change ratio is determined by the correct matching pair with the largest transformation distance, that is, by a3 and b3. Specifically, Figure 6 (a) is the image to be corrected, Figure 6 (b) is the reference image, Figure 6 The triangle formed by a1-a2-a3-a4-a5 in (a) is Figure 6 The triangle formed by b1-b2-b3-b4-b5 in (b) undergoes stretching transformation, and the shape of the triangle changes. The embodiment of the present invention achieves precise registration of SAR images under different time phases, speckle interference, and terrain undulations, improving registration accuracy and reliability. The embodiment of the present invention also uses five pairs of polarimetric SAR image data for experiments. The sensors for image pairs 1, 2, and 4 are all GF3, with a resolution of 8m*8m. The size of the first image pair is 1000*1000, and the sizes of the second and fourth image pairs are 900*900. The sensors for the third image pair are TerraSAR-X and RADARSAT2, with resolutions of 2.06m*6.59m and 4.73m*4.80m, respectively, and sizes of 685*1351 and 821*631, respectively. The fifth image pair is from the ALOS2 and RADARSAT2 sensors, with resolutions of 6m*6m and 4.73m*4.82m, respectively, and sizes of 2876*2205 and 1749*1180, respectively. Image pairs 1, 2, and 3 depict flat terrain, while image pairs 4 and 5 depict undulating terrain.

[0084] In the pre-registration stage, the threshold used in the NNDR matching method is set to 0.9, τ1 is set to 0.05, τ2 is set to 0.05, τ3 is set to 1.5, and τ max is set to 20. When τ1, τ2 and τ3 are low, it is difficult to terminate the iteration, and when they are high, it is difficult to obtain a more accurate registration result. In the fine matching stage, the local search region is set to [-15, 15]. The pre-registered image obtained in the pre-registration stage is divided into 5*5 non-overlapping image blocks, k is set to 25, which refers to the number of feature points in each block. When constructing the similarity measure of PCGCC, an image block with a size of 71*71 is used.

[0085] Taking the third image pair as an example, the proposed method detects the results as shown in Figure 3 . Among them, (a) and (b) respectively represent the image to be registered and the reference image, (c) and (d) are the pre-registration and fine registration results displayed in the chessboard diagram. (e) and (f) are the first and second sub-images of the pre-registration, and it can be seen that the pre-registration stage has completed the rough alignment, and there is still an error in the connection of the edges. (g) and (h) are the first and second sub-images of the fine registration, and it can be seen that the edges are well connected, indicating that the overall accurate registration is achieved.

[0086] Taking the fourth image pair as an example, the proposed method detects the results as shown in Figure 4 . Among them, (a) and (b) respectively represent the image to be registered and the reference image, (c), (d) and (e) are the pre-registration results, the fine registration global transformation results and the fine registration local stretching transformation results displayed in the chessboard diagram. (f) and (g) are the first and second sub-images of the pre-registration, and it can be seen that the pre-registration stage has completed the rough alignment, and there is still an error in the connection of the edges. (h) and (i) are the first and second sub-images of the fine registration global transformation, and it can be seen that the edge connection of (h) is good, while the edge connection of (i) is poor, indicating that after the fine registration global transformation, the flat terrain area achieves accurate registration, while the terrain undulating area has poor registration results. (j) and (k) are the first and second sub-images of the fine registration local stretching transformation, and it can be seen that both (j) and (k) achieve accurate connection, indicating that the terrain undulating area is correctly corrected, and the overall accurate registration is achieved.

[0087] Referring to Figure 5 , the scale-invariant feature transform matching method (SIFT), the bilateral filter SIFT (BF-SIFT), the SAR-SIFT algorithm, the improved SAR-SIFT method (I-SAR-SIFT) and the registration method disclosed in the embodiments of the present application are compared, and the comparison results are shown in Figure 5It can be seen that the number of matching pairs of five pairs of images of the embodiment of the present application is 285, 103, 164, 235, and 228, respectively, which is much higher than that of other methods. The RMSE is 0.55, 0.74, 0.69, 0.66, and 0.85, respectively, all of which are less than 1, while the RMSE of other methods is all above 2, indicating that the method proposed in the embodiment of the present application can achieve more accurate registration. In terms of time, it is 27.77, 37.28, 32.27, 79.93, and 71.83, respectively, which is not necessarily optimal compared with other methods, because the proposed method is a two-stage method and is processed for the terrain undulating area, which requires more processing time. But overall, it can be seen that the method disclosed in the embodiment of the present application can achieve accurate registration.

[0088] The above is only a specific embodiment of the present application, but the protection scope of the present application is not limited thereto, and any person skilled in the art can easily think of various equivalent modifications or replacements within the technical range disclosed by the present application, and these modifications or replacements should be covered within the protection scope of the present application.

Claims

1. A two-stage SAR image registration method based on SAR-SIFT and template matching, characterized in that: It includes the pre-registration stage and the fine registration stage, where: In the pre-registration stage, the downsampled image is subjected to the SAR-SIFT algorithm to obtain an initial registration image; the affine transformation matrix of the initial registration image is iteratively updated using an iterative optimization matching method to obtain a final matching pair; the final affine transformation matrix is ​​mapped back to the original image size to obtain the final affine transformation matrix, and the pre-registration image is obtained based on the final affine transformation matrix; In the fine registration stage, the pre-registered image is divided into blocks to extract stable feature points and the similarity measurement is combined to eliminate the mismatched pairs in the reference image template to obtain a fine-registered image; by clustering the mismatched pairs, it is determined whether there is an area with large terrain undulations. If there is an area with large terrain undulations in the fine-registered image, the candidate terrain undulation areas in the fine-registered image are sampled and local correction of the undulation areas is performed to obtain the fine-registered image again.

2. The two-stage SAR image registration method based on SAR-SIFT and template matching according to claim 1, characterized in that: The method of the pre-registration stage specifically includes: after downsampling the input SAR image, a set of feature points is obtained based on Harris corner points in the spatial pyramid through Gaussian blurring, the gradient information in different directions is calculated for each feature point to extract the feature description vector, and the initial registration image is obtained through the nearest neighbor distance ratio and fast sampling consistency; the affine transformation model of the initial registration image is iteratively optimized and updated by an iterative optimization matching method until the model tends to be stable, and corresponding matching pairs are obtained; the matching pairs are mapped back to the original image size to obtain the final affine transformation matrix; and the pre-registration image is obtained according to the final affine transformation matrix.

3. The two-stage SAR image registration method based on SAR-SIFT and template matching according to claim 1, characterized in that: In the pre-registration stage, the downsampled image is subjected to the SAR-SIFT algorithm to obtain the initial registered image, which includes the following steps: Feature point detection: The input SAR image is downsampled to obtain a downsampled image, a multi-scale space is constructed, and the downsampled image is smoothed using Gaussian kernel functions of different scales, and a Harris scale space pyramid is established. The gradient amplitude of the pixel points in the multi-scale space is calculated based on the local window area. Each pixel in each layer of the constructed Harris scale space pyramid is compared and candidate feature points are screened. The candidate feature points are accurately located and low-contrast pixels are removed to obtain a stable feature point set. Main direction assignment and descriptor extraction: A gradient direction histogram is constructed for each detected feature point, and the direction with the largest peak value is selected as the main direction of the feature point. A circular neighborhood is constructed with the feature point as the center, and the gradient information of each direction in each area is statistically combined to form a feature description vector. Feature point matching: For each feature point in the image to be matched, the nearest neighbor distance and the next nearest neighbor distance of two feature points are found in the reference image and their ratio is calculated. If the ratio allows, the matching pair is accepted to improve the matching accuracy. The affine transformation model is calculated and iteratively optimized to select the optimal solution of the model with the largest number of corresponding internal points for feature point matching.

4. The two-stage SAR image registration method based on SAR-SIFT and template matching according to claim 1, characterized in that: In the pre-registration stage, the affine transformation model of the affine transformation is formulated as follows: In formula (1), a, b, c, d, e, and f are the affine transformation coefficients, and (x1, y1) and (x2, y2) are the coordinates before and after the transformation. It is simply assumed that the closer the coordinates (x1, y1) are to (x2, y2), the more stable the transformation is. Otherwise, the matching pair obtained is not stable. At this time, the transformation matrix tends to the unit matrix. The iteration is controlled by limiting the affine transformation coefficients. For the transformation coefficients, set |a-1|<τ1&|b|<τ2&|c|<τ3 or when it is greater than the maximum number of iterations τ max The iteration stops when , and the coefficients d, e, and f have the same restrictions as b, a, and c respectively.

5. The two-stage SAR image registration method based on SAR-SIFT and template matching according to claim 1, characterized in that: In the fine registration stage, the pre-registered image is matched in a local area to obtain a fine registered image, including the following steps: Feature point extraction: In the fine registration stage, the pre-registered image is divided into blocks, and a Harris corner point-based method is used to extract a set of feature points. A threshold is set to extract stable feature points. Combining phase consistency and gradient similarity measurement: The similarity measurement of the stable feature points extracted from the blocks is determined based on phase consistency and gradient information. The stable feature points with high similarity values ​​are used as matching points. Mismatched pairs are eliminated in the reference image template for fine registration to obtain the initial finely registered image. Local correction of terrain undulation areas: Cluster the matching pairs that do not meet the final transformation model to obtain cluster matching pairs corresponding to different clusters, which are used as candidate terrain undulation areas. The boundaries and surrounding areas of the candidate terrain undulation areas are sampled and accurately calculated to further obtain more accurate local undulation areas. Then, each local undulation area is sampled and template matched to obtain the final stretch transformation ratio. Local correction is performed on each local undulation area to obtain a finely registered image.

6. The two-stage SAR image registration method based on SAR-SIFT and template matching according to claim 5, characterized in that: The calculation formula of the similarity measure combining phase consistency and gradient information is as follows: PC_G(x,y) = ∑ o (PC 2o (x,y) + G o (x,y)) Formula (2) In formula (2), PC 2o (x,y) and G o (x,y) represents the edge information obtained in the direction o using phase consistency in the frequency domain and gradient in the spatial domain respectively; where: The phase consistency model in the frequency domain is defined as: In formula (3), A n (x,y) represents the amplitude at scale n, e n (x,y) represents the convolution of the image with an even-symmetric filter, n represents the number of scales, I represents the image intensity, represents an even symmetric filter at scale n, o n (x,y) represents the convolution of the image with an odd-symmetric filter, represents an odd-symmetric filter at scale n, Indicates the corresponding phase, W o represents the weighting function at the oth orientation, o represents the number of orientations, A no represents the magnitude at scale n and orientation o, T o represents the noise level at the oth orientation, and ε represents a small positive value to prevent no When it becomes very small, the expression becomes unstable, ΔΦ no (x,y) is represented as: For PC2(x,y), for a specific direction, it means the phase consistency in that direction, and the edge in that direction is obtained. The following formula is the edge in a specific direction after taking a certain value. If a total of 6 directions are taken, the value of o is 0, 1, 2, ..., 5, the orientation interval is 180° / 6=30°, and the direction is o*30°; The gradient information in the spatial domain is: In formula (5), I represents the image intensity, and They represent the local mean on both sides of the point (x, y) direction, and the weights are assigned using the Gaussian weighting method, where α represents the standard deviation of the Gaussian kernel; G o (x,y) represents the ratio gradient obtained in direction o.

7. The two-stage SAR image registration method based on SAR-SIFT and template matching according to claim 5, characterized in that: The matching pairs that do not meet the final transformation model are clustered to obtain the corresponding cluster matching pairs of different clusters, which are used as candidate terrain relief areas. Local adaptive stretching transformation is performed on each relief area to correct it locally, which specifically includes: Based on the distribution characteristics of the recovered matching pairs, a rectangular bounding box is calculated for each area of ​​undulating terrain drawn within the transformation range. Sampling is then done from the bounding box to both sides. If the area is in a flat region, the sample is sampled towards the undulating region; if the area is in an undulating region, the sample is sampled towards the flat region. The spatial distribution characteristics of the sampling points are used to identify the precise boundary of the undulating terrain area. The distance between the boundary sampling points is used to determine whether adjacent sampling points belong to independent regions or the same continuous region, and the region is further subdivided into multiple specific sub-regions. The affine transformation model is used to roughly delete the mismatched pairs obtained by uniform sampling. For each finely divided sub-region, the stretching transformation ratio is calculated based on the correct matching pair with the largest transformation distance in the region to ensure that the stretching transformation fully adapts to the local variation characteristics of terrain relief correction. Local adaptive stretching transformation is implemented for each specific terrain relief area according to the calculated stretching transformation ratio.

8. The two-stage SAR image registration method based on SAR-SIFT and template matching according to claim 7, characterized in that: In the local adaptive stretching transformation implemented above, the direction transformation formula is expressed as: In formula (6), x left Indicates the left boundary of the terrain relief area in the x direction, x right represents the right boundary of the terrain relief area in the x direction, x1 represents the x coordinate of the transformation pair determined by the pre-registered image, x2 represents the x coordinate of the transformation pair determined by the reference image, and s left and s right represents the determined stretching ratio in the x direction, and x′ represents the coordinate after stretching transformation at position x.

Images