Multi-Hypothesis Topological Isomorphism Matching Method for Synthetic Aperture Radar Images with Large Geometric Distortion

Abstract

The dramatic undulations of a mountainous terrain will introduce large geometric distortions in each Synthetic Aperture Radar (SAR) image with different look angles, resulting in a poor registration performance. To this end, this paper proposes a multi-hypothesis topological isomorphism matching method for SAR images with large geometric distortions. The method includes the Ridge-Line Keypoint Detection (RLKD) and Multi-Hypothesis Topological Isomorphism Matching (MHTIM). Firstly, based on the analysis of the ridge structure, a ridge keypoint detection module and a keypoint similarity description method are designed, which aim to quickly produce a small number of stable matching keypoint pairs under large look angle differences and large terrain undulations. The keypoint pairs are further fed into the MHTIM module. Subsequently, the MHTIM method is proposed, which uses the stability and isomorphism of the topological structure of the keypoint set under different perspectives to generate a variety of matching hypotheses, and iteratively achieves the keypoint matching. This method uses both local and global geometric relationships between two keypoints, hence it achieving better performance compared with traditional methods. We tested our approach on both simulated and real mountain SAR images with different look angles and different elevation ranges. The experimental results demonstrate the effectiveness and stable matching performance of our approach.

1. Introduction

About 24% of the earth’s land is covered by mountains [1]. Since NASA launched its first SAR satellite SEASAT in 1978, several countries have successively deployed multiple spaceborne SAR systems, accumulating massive amounts of SAR image data of mountain areas. In order to jointly exploit these data for elevation inversion, deformation detection, and biomass monitoring, an accurate matching performance becomes a prerequisite. However, the SAR imaging mechanism determines that a mountainous SAR image is a slope-distance mapping of the mountain from a three-dimensional space to a two-dimensional image. The difference in the viewing angles causes a relative geometric distortion between two images. In particular, the larger the difference in the angles, the larger the geometrical deformations. This poses great challenges to the registration of SAR images with large geometric distortion.

Increasing efforts have been made to improve the accuracy of registration. According to a measuring function, an acceptable classification [2] for existing SAR image matching methods is area-based [3,4,5,6,7,8,9] and feature-based [10,11,12,13,14,15,16,17] pipelines. The area-based methods either use image grayscale statistical information or transform domain statistical information as a measure, and register the image by searching for the maximum value of the measuring function. These types of methods are more suitable for the cases with a small difference in look angles, that is, the cases where the coherence between images is high, and widely used in interferometric SAR (InSAR) [18,19] and SAR tomography (TomoSAR) [20,21]. However, this type of SAR image accounts for only a small proportion. Since the topography of the imaging area increases sharply with increasing difference in SAR imaging angles, the relative deformations of the images increase and the correlation between them decreases. In this case, a feature-based method is more suitable. The feature-based methods first extract spatial features such as corners, contours, and line segments in an image, then produce homologous points by some matching methods, and finally calculate a transfer model between images based on the homologous points, so as to complete the matching. Among many feature-based methods, the Scale-Invariant Feature Transform (SIFT) [10] is the most widely used method. It is capable of coping with geometric transformations, such as scaling and rotation. However, its real-time performance is limited, and it cannot handle smooth edges well. Based on SIFT, Speeded-Up Robust Features (SURF) [16] use the Haar wavelet response and integral image to improve the computational efficiency. However, the performances of SIFT and SURF in SAR image matching are still not very satisfactory. In order to reduce the impact of coherent speckle on the SAR image processing, Dellinger et al. [12] proposed the SAR-SIFT method based on a new gradient definition. Subsequently, Ma et al. [22] proposed an improved SIFT and enhanced the feature matching method named PSO-SIFT to enhance the adaptability of features. Xiang et al. [23] proposed a high-level rotation invariant descriptor that does not specify the master direction to improve SIFT, and improved the discrimination of different feature descriptors. Apart from the SIFT-like methods, DASIY, BRIEF, ORB and other methods [24] also complete the optical image registration under small differences in viewing angles.

Most of the existing approaches work on the improvement of general feature descriptors. The matching algorithms for large geometric deformation of SAR images are still in their infancy. Although SAR-SIFT and PSO-SIFT [12,22] have achieved good results for mountain SAR image registration with small differences in viewing angles, they remain challenging when registering large geometric distortion of SAR images with only a small number of keypoints (refer to the experimental part in Section 3 for detail). For small-deformation SAR image registration, fewer keypoints can complete a global high-precision fitting of the matching transformation model. However, it is difficult for a matching transformation model of large distortion of SAR images to perform a global high-precision approximation through a linear affine transformation and a low-order nonlinear transformation. In this case, it is important to know how to improve the quality and quantity of keypoints, so as to increase the matching accuracy.

The ridge is one of the cardinal topographic features, which forms the skeleton line of an undulating topography and reflects the spatial topological structure of the topography. For solving the problem of large geometric distortion of SAR image matching in mountainous areas caused by large differences in look angles and severe terrain fluctuations, we propose a large geometric distortion SAR image multi-hypothesis topological isomorphism matching method. It is composed of keypoints of ridges Detection (Ridge-Line Keypoint Detection, RLKB) and Multi-hypothesis Topological Isomorphism Matching (MHTIM). The method can produce more matching keypoints with better stability, so as to achieve a higher matching precision.

This paper is arranged as follows. In Section 2, we introduce the motivation and ideas of the proposed method and discuss the two parts of the method in detail. In Section 3, we present the simulation results and measured data. Finally, in Section 5, we summarize the method and give directions for future improvements.

2. Methods

The overall flowchart of our proposed method is shown in Figure 1. In the rest of this section, we first analyze the deficiency of existing methods and present the motivation of the proposed method in Section 2.1, and then explain RLKB and MHTIM in Section 2.2 and Section 2.3, respectively.

2.1. Problem Description

Despite its long success in optical image matching, the feature-based method still has vulnerability in detecting and matching the ridge features of most of the parts of SAR images with large geometric distortion. In the case of a SIFT-like method, there are two reasons: (1) The method constructs the Difference of Gaussian (DoG) pyramid of the image when considering the image scale. The position of the keypoints in the image has an offset relative to that of the ridge. So, the keypoints cannot represent the position of the ridge. (2) The method generally uses information such as the gray gradient direction of the small image blocks around the keypoint as the descriptor, and calculates the similarity of the descriptor (distance between two vectors) to determine whether the keypoints represent a homologous object. In fact, in mountain areas, when the look angle of the SAR image changes significantly, except for the ridge line with a larger scale, other areas of the image have significant changes in brightness, shape and even phase, which make the similarity of the descriptor invalid. Similar to the intuitive experience, the topological structure of the ridge features in the SAR image at different look angles is isomorphic.

Analyzing Figure 2, it is observed that the distributions of the extreme points of the image intensity formed by SAR images with different look angles on ridges are isomorphic. Figure 3 shows the SAR image of the mountainous area of the Sichuan-Tibet Plateau in China, where the DEM data of the DEM map, ascending stripe mode SAR image from Sentinel 1 and descending stripe mode SAR image from Sentinel 1 are shown in a–c, respectively. Sub-figures a–c in Figure 3 have undergone rough geometric registration. It is worth mentioning that geometric registration can roughly overlap the regions to be registered to enhance the efficiency of subsequent algorithms, and when the images overlap (like the images produced by TanDEM-X), geometric registration is not necessary. In Figure 3, two images are taken from the opposite-side of the terrain, and the angle between the line of sight of (b,c) is greater than 90°. We can find even more intuitively that when SAR images are taken from opposite-side, the topological structures composed of yellow circles and yellow lines in the three figures are still isomorphic. Therefore, this paper attempts to match the ridge features through the topological relationship between the ridge features. It divides the method into two parts: ridge keypoint detection and matching.

In geometry, the feature information that characterizes the mountain morphology includes extreme points, saddle points, ridge lines, valley lines, etc. [25]. The most apparent feature in mountain SAR images is the ridge line and the intersection of two or more ridge lines. Since the SAR image is a slant distance mapping from a three-dimensional space to a two-dimensional image, it is necessary to find a stable ridge line detection method. The ridge line is a kind of edge information. Commonly used edge detection operators include algorithms such as Sobel, LoG, Canny, etc. [26], followed by deep learning methods such as Holistically-nested Edge Detection (HED) [27] and Richer Convolutional Features (RCF) [28]. The location of the intersection of lines is generally obtained by calculating its surrounding curvature, gradient, etc. [13]. However, the time complexity of the above methods is relatively high. This paper proposes a ridge feature keypoint detector based on the LoG operator, which can produce stable ridge keypoints while detecting ridge lines at a lower time complexity.

The keypoint matching is a key step to complete image matching. Existing feature-based pipelines generally use some algorithms to solve the assignment problem [29] after calculating the similarity matrix of the master and slave image descriptor sets, and then use the RANdom SAmpling Consensus (RANSAC) [30] algorithm for outlier culling. This type of pipelines generally pays attention only to the similarity between keypoint descriptors (the distance between vectors) and ignores the inherent geometric topological relationship between the keypoints. Although RANSAC is simple and clear algorithm, the calculation cost in its iterative process is relatively large and the optimal solution cannot always be found [31].

Analyzing Figure 3, it is found that the distributions of the extreme points of the image intensity formed by SAR images with different look angles on ridges are still isomorphic. Therefore, this paper proposes a Multi-Hypothesis Topological Isomorphism Matching (MHTIM) method. This method converts the stable keypoint matching pairs generated by RLKD into an initial topological structure graph hypothesis according to its topology. Based on this, the method iteratively introduces the remaining unmatched keypoints to form a hypothesis tree. When the hypothesis tree reaches a certain depth, the hypothesis score is calculated, and the hypothesis tree is pruned to gradually complete the matching process.

2.2. Ridge Line Keypoint Detection Method

The RLKD method is divided into three parts: (1) Quick detection of the ridge line intersection point, which ridge detection is performed in the distance and azimuth direction, respectively, to quickly obtain the ridge intersection point; (2) keypoint generation and description, which cluster the intersection point pixels to produce the keypoint, and a keypoint descriptors are designed to measure their similarity; and (3) fast matching, which calculates the distance matrix of ridge keypoints through the descriptor, and uses the simulated annealing algorithm to solve the two-allocation problem for obtaining a small number of stable keypoint matching pairs. As there exist many mathematical operators in the following passage, for convenience, we define all the notations in

Table 1. Notation.

Notation	Definition
$rDec$,$aDec$	LoG operator in the range and azimuth direction
$G\left(r,a,\sigma \right)$	Two-dimensional Gaussian filter
$\sigma$	Standard deviation
$I_{rDec}$,$I_{aDec}$	Responses of image function I through $aDec$ and $rDec$
$I_{re}(r,a)$, $I_{ae}(r,a)$	Ridge detection in range and azimuth direction
$I_E$	The final edge detection result
$k_i=\left(r_i,a_i\right)$	The coordinate of edge intersection in $I_E$
v	Keypoint produced by RLKD
s	Keypoint descriptor produced by RLKD
$d\left(s_1,s_2\right)$	The similarity between two descriptors
${\left(\right)}^{\ast }$	The conjugate operation
$\mathcal{C}$$\left(s_1,s_2\right)$	The real correlation coefficient matrix of $s_1$ and $s_2$
D	The similarity matrix
$V^m=\left\{v_i^m\right\},$$i=1,...,N_m$	Keypoints detected in the master image
$V^s=\left\{v_j^s\right\},$$j=1,...,N_s$	Keypoints detected in the slave image
$V_{matched}^m=\left\{v_k^m\right\}$, $k=1,...,r$	The matching keypoint set of master image
$V_{unmatched}^m=\left\{v_p^m\right\},$$p=r+1,...,N_m$	The candidate keypoint set of master image
$V_{matched}^s=\left\{v_k^s\right\}$, $k=1,...,r$	The matching keypoint set of slave image
$V_{unmatched}^s=\left\{v_q^s\right\},$$q=r+1,...,N_s$	The candidate keypoint set of slave image
$D^m={\left(d_{i_1,i_2}\right)}_{r\times r}$	The putative self-distance matrices of $V_{matched}^m$
$D^s={\left(d_{j_1,j_2}\right)}_{r\times r}$	The putative self-distance matrices of $V_{matched}^s$
$G_m$	The master undirected weighted graph
$G_s$	The slave undirected weighted graph
$\mathcal{S}\left(p,q\right)$	The closeness between $v_p^m$ and $v_q^s$ in their respective graphs
$\alpha (p,q)$	The node angle similarity
$\mathfrak{C}=\left\{v_k^m,v_k^s\right\}$, $k=1,2,...,r$	The final matching pair set

.

2.2.1. Quick Detection of Intersection of Ridge Lines

Our method is based on the LoG to rapidly detect the intersection of ridge lines by using two detectors $rDec$(Detector in range) and $aDec$(Detector in azimuth), which are defined as follows:

$$\left.\begin{array}{c}rDec=\frac{{\partial }^2}{\partial r^2}G\left(r,a,\sigma \right)=\frac{r^2-{\sigma }^2}{{\sigma }^4}{e}^{\frac{-\left(r^2+a^2\right)}{2{\sigma }^2}}\\ aDec=\frac{{\partial }^2}{\partial a^2}G\left(r,a,\sigma \right)=\frac{a^2-{\sigma }^2}{{\sigma }^4}{e}^{\frac{-\left(r^2+a^2\right)}{2{\sigma }^2}}\end{array}\right\}$$

(1)

Among them, $G\left(r,a,\sigma \right)$ is a two-dimensional Gaussian filter:

$$G\left(r,a,\sigma \right)=\frac{1}{2\pi {\sigma }^2}{e}^{-\left(r^2+a^2\right)/2{\sigma }^2}$$

(2)

In the above formula, $\sigma$ is the standard deviation. Due to the low-pass characteristics of the Gaussian filter, fine textures can be eliminated and large-scale ridge features can be retained while suppressing the influence of coherent speckles. In addition, if not otherwise stated, r and a represent the distance pixel index and azimuth pixel index of the image, respectively. Next, the intersection point is obtained based on the detected ridge line. Assuming that the image gray function is I, $I_{rDec}$ and $I_{aDec}$ as the responses of I can be obtained through $aDec$ and $rDec$ as follows:

$$\begin{array}{c}{I}_{rDec}=I\ast rDec,\hfill \\ I_{aDec}=I\ast aDec\hfill \end{array}$$

(3)

where ∗ is the convolution operation. The edges of the image along the distance and the azimuth directions exist at the position where the two adjacent pixel values have different signs. The results of edge detection along the two directions can be produced by using the following two equations:

$$\begin{array}{c}{I}_{re}(r,a)=\left\{\begin{array}{cc}0,\hfill & I_{rDec}(r,a)·I_{rDec}(r,a+1)\ge 0\hfill \\ 1,\hfill & I_{rDec}(r,a)·I_{rDec}(r,a+1)<0\hfill \end{array}\right.\hfill \\ I_{ra}(r,a)=\left\{\begin{array}{cc}0,\hfill & I_{aDec}(r,a)·I_{aDec}(r+1,a)\ge 0\hfill \\ 1,\hfill & I_{aDec}(r,a)·I_{aDec}(r+1,a)<0\hfill \end{array}\right.\hfill \end{array}$$

(4)

The pixel with value of 1 in $I_{re}$ and $I_{ra}$ is the edge. Now, $I_{re}$ and $I_{ra}$ are to be superimposed, that is:

$$I_E=I_{re}+I_{ra}$$

(5)

$I_E$ represents the final edge detection result. In $I_E$, the pixel with the value greater than 0 is the ridge, and the pixel with the value of 2 is the intersection of the ridge lines.

2.2.2. Keypoint Generation and Descriptor

A real ridge intersection generates a concentrated area in $I_E$ with a value greater than 1, which contains a lot of information about the ridge intersection. This paper uses the block of this area as a keypoint descriptor, clusters these pixels, and uses the cluster center as the keypoint. Suppose that in $I_E$, the number of intersections generated is T, and $k_i=\left(r_i,a_i\right)$ are the coordinates of the intersection in the image, where $i=1,2,\dots ,T$. For the first intersection, a traversal search is performed on the image block centered on this point with a fixed radius, any other intersection in the image block will center on itself, continue to search for other intersections within the same radius, and cluster to C. Additionally, $C_j$ represents the set of subscripts of the intersection points contained in the jth cluster. Equation (6) can be used to calculate the jth cluster center, which is the keypoint required by the algorithm in this paper:

$$v_j=\frac{1}{\left|C_j\right|}{\sum }_{i\in C_j}{k}_i$$

(6)

This paper uses a $R\times R$9×9 $I_E$ block s with the keypoint position as the center of the keypoint descriptor, where R is the block size, which we set to 9. This descriptor retains the shape information of the ridge line around the keypoint. For two image blocks, the measures of similarity include correlation, mutual information and cross entropy. The correlation and mutual information are two basic and effective methods of similarity measurement. The cross entropy [32] has high robustness and it is widely used in the loss function in the framework of deep learning algorithms. However, it is usually used to calculate the similarity of multi-spectral images with complex textures. As shown in

Table 2. Results of different descriptor measuring functions in SAR image registration on the synthetic data set with a difference of 5${}^{°}$ in look angles.

Method	MAE (pixel)	NKM	TIME
Correlation	0.59	21	0.90 s
Mutual information	0.67	18	0.84 s
Cross entropy	0.63	20	0.88 s

, for SAR images with large geometric distortion, we get the best result by using cross-correlation as the descriptor similarity measurement method. The reason is explained in more detail in Table 2 in Section 3. In particular, the larger the maximum value of the correlation coefficient, the more similar the descriptors. Therefore, our method uses the maximum value of the correlation coefficient to represent the similarity between two descriptors, which is defined as follows:

$$d\left(s_1,s_2\right)=max\left(\mathcal{C}\left(s_1,s_2\right)\right)$$

(7)

In the above formula, $\mathcal{C}$ is the real correlation coefficient matrix of descriptors $s_1$ and $s_1$, obtained by using formula (8), where ${\left(\right)}^{\ast }$ represents the conjugate operation and R is the block size of s:

$$\mathcal{C}\left(s_1,s_2\right)=\left(\sum _{u=1}^R\sum _{v=1}^R\left|\sum _{r=1}^R\sum _{a=1}^R{s}_1\left(r,a\right)e^{-2j\pi \left(\frac{ur+va}{R}\right)}·{\left(\sum _{r=1}^R\sum _{a=1}^R{s}_2\left(r,a\right)e^{-2j\pi \left(\frac{ur+va}{R}\right)}\right)}^{\ast }\right|e^{2j\pi \left(\frac{ur+va}{R}\right)}\right)$$

(8)

In Figure 4, we show the SAR-SIFT keypoint detection results of the SAR image pair in the mountain area under the typical parameters reported by Dellinger et al. [12]. We present the keypoint detection and description results of our method in Figure 5. We can find in Figure 4 that the detection results of keypoints in the mountain area, obtained by the SAR-SIFT method, have two characteristics. Firstly, for images with different look angles, the number of characteristic points is obviously different, and second, the keypoints mainly exist on the back slope and are not even.

Figure 5a,b show the ridge line and keypoint detection results of the SAR image shown in Figure 4a,b, respectively. In (a,b), the green line is the detection result of the ridge line, and the yellow dot is the keypoint before clustering; (c,d) are corresponding results after clustering. The red boxes marked in (a,b) and (c,d) are the positions of the manually selected keypoints to reveal the descriptor and enlarged patches in the 3rd row. Keypoints 1 and 2 of the master image correspond, respectively, to Keypoints 1 and 2 of the same ground object in the slave image. The 4th row shows the calculated maximum value of the two-dimensional correlation coefficient of the descriptor. It can be seen that the maximum values of the correlation coefficient between keypoint 1 of the master image and keypoints 1 and 2 of the slave image are 0.86 and 0.74, respectively; while the maximum correlation coefficients between keypoint 2 in the master image and keypoints 1 and 2 in the slave image are 0.64 and 0.84, respectively. That is, keypoint 1 of the master image matched more with keypoint 1 of the slave image, and keypoint 2 of the master image matched more with keypoint 2 of the slave image, which is in line with the actual situation. Figure 5 shows that the descriptor designed in this paper describes the similarity of keypoints effectively.

2.2.3. Quick Matching

Assuming that M and N keypoints are detected in the master and slave images, respectively, an $M\times N$-dimensional similarity matrix can be generated as $D=\left\{d\left(s_i,s_j\right)|i=1,2,...,M,j=1,2,\dots ,N\right\}$, where s is the keypoint descriptor. Finding the best match between the keypoints of the master and slave images in the similarity matrix is a typical optimization problem. Our method uses the simulated annealing algorithm [33] to solve the optimization problem. Finally, the method proposed in this section still uses the RANSAC [30] algorithm to delete outliers of matched pairs, and quickly finds stable keypoints. It is worth mentioning that the application of the transformation model to the image distortion mode in the RANSAC algorithm is very important for the registration result. The two most basic transformation models, used in the task of large-distortion SAR image registration, are projection transformation and affine transformation. For more complex scenes, polynomials can be used to fit a transformation model. Since the geometric distortion of a SAR image is more in the slant range direction, the Local Weighted Mean (LWM) [34] transformation model is used. The LWM model introduced by Goshtasby [34] is capable of the condition when parts of the image appear distorted differently, and the distortion varies locally and piecewise in a nonlinear manner.

In Section 3, we compare the pros and cons of several models, and verify that the LWM model is indeed the most suitable model for the scenario targeted in this article.

2.3. Multi-Hypothesis Topological Isomorphism Matching Method

The main aim of the MHTIM method is to transform the stable keypoints matching from RLKD into two initial graph hypotheses according to their topological structures. Based on the ridge feature topological isomorphism, it iteratively adds the remaining unmatched keypoints to form a hypothesis tree. When the hypothesis tree reaches a certain depth, the hypothesis score is calculated. Then, the hypothesis tree is pruned according to the score, and subsequently the matching process is completed.

In this section, we first introduce the initialization of the hypothesis tree and the sorting of candidate points. The subsequent iteration process is introduced one by one in three steps as follows:

• Multi-hypothesis generation: According to the ranking results, select the top candidate keypoints and add them to the graph to generate new hypotheses.
• Hypothesis score calculation: We use five graph indicators and node angle similarity indicators to rank the new hypotheses. The hypothesis that ranks first in a single indicator gets a certain score. The final score of a new hypothesis is the sum of the scores assumed under each indicator.
• Pruning: New hypotheses are sorted in terms of their final scores, pruning low-scoring hypothesis branches, retaining high-scoring hypothesis branches, and updating the root node, matching, and candidate point sets.

A schematic diagram of the above steps is shown in Figure 6. Through the above steps, the keypoints that are not matched but have potential matches are gradually added to the matching set.

2.3.1. Hypothesis Initialization and Candidate Keypoints Sorting

Before starting iterations, the MHTIM method first initializes the root node of the hypothetical tree. The specific process is shown in Figure 7.

We use $v=\left(r,a\right)$ to represent the coordinates of the keypoint in the image. Suppose that after RLKD, keypoints detected in the master and slave images are defined as point sets $V^m=\left\{v_i^m\right\};i=1,\cdots ,$ $N_m$ and $V^s=\left\{v_j^s\right\};j=1,\cdots ,$ $N_s$, respectively. After the RLKD algorithm quick matching, $V^m$ and $V^s$ are divided into two parts: matching point set $V_{\mathrm{matched}}^m=\left\{v_k^m\right\}$ and $V_{\mathrm{matched}}^s=\left\{v_k^s\right\};k=1,\cdots ,r$, and candidate point set $V_{\mathrm{unmatched}}^m=\left\{v_p^m\right\};$$p=r+1,\cdots ,N_m$ and $V_{\mathrm{unmatched}}^s=\left\{v_q^s\right\};$$q=r+1,\cdots ,N_s$. $D^m={\left(d_{i_1,i_2}\right)}_{r\times r}$ and $D^s={\left(d_{j_1,j_2}\right)}_{r\times r}$ are the putative self-distance matrices of $V_{\mathrm{matched}}^m$ and $V_{\mathrm{matched}}^s$, respectively, where $d_{i_1,i_2}={\Vert v_{i_1}^m-v_{i_2}^m\Vert }_2$ and $d_{j_1,j_2}={\Vert v_{j_1}^s-v_{j_2}^s\Vert }_2$.

First of all, this article quickly matches the keypoint pair initialization, generates the master and slave undirected weighted graphs $G_m$ and $G_s$, respectively, and $v^m$ and $v^s$ are regarded as nodes in $G_m$ and $G_s$, respectively. After that, if $d_{i_1,i_2}$ or $d_{j_1,j_2}$ is greater than $\epsilon$, an edge $e_{i_1,i_2}^m=(v_{i_1}^m,v_{i_2}^m)$ or $e_{j_1,j_2}^s=(v_{j_1}^s,v_{j_2}^s)$ is generated, where $\epsilon$ is the threshold value set according to the number of pixels in the image. Generally, $\epsilon$ is set to one half of the diagonal length of the image as follows:

$$\epsilon =0.5\sqrt{h^2+w^2}$$

(9)

In the above formula, h and w are the number of pixels of the image in the range and azimuth directions, respectively. As shown in Figure 7, the master and slave initial graphs show strong topological isomorphism. Please note that the distance in Figure 7 is the unitless distance after a normalization operation, and $\epsilon$ is scaled in the same proportion.

We sort the candidate points by $\mathcal{S}$ to ensure the quality of the hypothesis generated by the selected candidate points in each iteration of the algorithm. The definition of $\mathcal{S}$ is as follows:

$$\mathcal{S}(p,q)=\left|\frac{1}{{\zeta }_m}\times {∥\begin{array}{c}{v}_p^m-{\displaystyle \sum _{v\in V_{\mathrm{matched}}^m}}\frac{v}{\left|V_{\mathrm{matched}}^m\right|}\end{array}∥}_2-{∥\begin{array}{c}{v}_q^s-{\displaystyle \sum _{v\in V_{\mathrm{matched}}^s}}\frac{v}{\left|V_{\mathrm{matched}}^s\right|}\end{array}∥}_2\times \frac{1}{{\zeta }_s}\right|$$

(10)

The right-hand side of Equation (10) contains the first and second normalized value terms, in which the distance from the candidate point $v_p^m$ of the master image to the geometric center of the matching keypoint set $V_{\mathrm{matched}}^m$, and the distance from the candidate point $v_q^s$ of the slave image to the geometric center of the matching keypoint set $V_{\mathrm{matched}}^s$. ${\zeta }^m$ and ${\zeta }^s$ are the normalization factors, which are set according to the area covered by the image on the ground and size of the image:

$$\zeta =\sqrt{\frac{h^2+w^2}{L_r^2+L_a^2}}$$

(11)

where, $L_r$ and $L_a$ are the length with a unit of meter of the area covered by the image on the ground in the range and azimuth direction, respectively. It is worth mentioning that when the master and slave images are geometrically registered and their scales are the same, ${\zeta }^m$ and ${\zeta }^s$ can be set to 1 at the same time.

Equation (10) can be understood in terms of the similarity of the distance from the candidate keypoints in the master and slave images to the respective geometric centers. The higher the similarity, the more likely the two candidate keypoints represent the same ridge feature.

2.3.2. Multi-Hypothesis Generation

Referring to Figure 6, we assume that the maximum depth of the tree is $H=3$, and the leaf nodes of the tree generate at most $W=2$ new hypothetical nodes at each iteration to illustrate the iteration process. After initialization, suppose that at the beginning of the $(k-3)$th iteration, the root graph of each of master and slave trees has 4 nodes (as shown in the $(k-3)$th layer in Figure 6). After the first $(k-2)$ iterations, the first node $v_{1st}^m$ in the sequence of the remaining candidate keypoints after sorting in the master graph is added to $G_m$. For the slave tree, the two points in $V_{unmatched}^s$ with the highest similarity to $v_{1st}^m$ are added to $G_s$ to form two hypotheses. At this point, the depth of the target hypothesis tree from the root node is 2. The above steps are reproduced sequentially in the $(k-1)$th and kth iterations. At k, the target hypothesis tree has a depth of 4, and there are at most 8 leaf nodes in the fourth layer. So far, in this example, the hypothesis tree has been generated.

We can find that the hypothesis tree retains multiple matching combinations. The following steps are to calculate the scores of the hypotheses for evaluating their qualities, and for pruning the hypothesis tree so as to remove the low-quality hypotheses and retain the correct ones.

2.3.3. Hypothesis Score Calculation

The score of a hypothesis comes from the similarity of the newly added vertices of the master and slave hypotheses. We use 5 common graph indicators and a custom indicator of the newly added nodes in the graph to measure the similarity of hypotheses. The five graph indicators are node centrality, betweenness centrality, proximity centrality, K kernel number, and eigenvector centrality. In addition to the above general graph indicators, the use of geometric constraints can enhance the matching accuracy of graph nodes in a realistic scene [35]. We define a node angle similarity index. First, an angle vector is used to describe the position of a node in the graph relative to the rest of the nodes. The vector is defined as:

$$\mathsf{\Phi }\left(v_i,G\right)=\left\{{\varphi }_l\right\},{\varphi }_l=arctan\left(v_i-v_l\right),v_l\in G,v_l\ne v_i$$

(12)

After that, the node angle similarity $\alpha (p,q)$ is defined in terms of the mean absolute value of the difference between the angle vectors of the newly added vertices of the master and slave graph, namely:

$$\alpha (p,q)=\frac{\left|\right|\mathsf{\Phi }(v_p^m,G_m)-\mathsf{\Phi }(v_q^s,G_s){\left|\right|}_2}{|G_m|-1}$$

(13)

where $|G_m|$ represents the number of nodes in the master graph at the kth iteration.

After the calculation of the 6 indicators of the newly added nodes is completed, we rank each hypothesis according to the similarity of each indicator.

After we get the 6 index values of the newly added vertices of the new hypothesis, we sort each new hypothesis according to the similarity of each index. The hypothesis that ranks high in an indicator gets a certain score, and the final score of the hypothesis is the sum of the scores on each indicator.

2.3.4. Pruning

The purpose of pruning is to remove hypotheses with lower scores in the hypothesis tree, and to update the root node to output a new pair of matched keypoints. After obtaining the hypothesis score at the kth iteration, the branch that does not contain the highest scoring hypothesis is deleted. After that, the k-H layer node of the reserved branch is used as the new root node, and the matching point which added by this node is output to $V_{\mathrm{matched}}^m$ and $V_{\mathrm{matched}}^s$. Finally, the point is deleted from $V_{\mathrm{unmatched}}^s$ and $V_{\mathrm{unmatched}}^s$, and the next iteration is started.

After all unmatched keypoints participate in the iterative process, the MHTIM terminates the iteration process and outputs the final hypothesis with the highest score. At this point, all the vertices in the master and slave graphs correspond one-to-one, and the final matching pair set $\mathfrak{C}=\left\{v_k^m,v_k^s\right\}$, $k=1,2,...,r$ is formed.

3. Experiment

In this section, we use both simulated and measured data to verify the performance of the two steps of our proposed method. The simulated data simulates a pair of SAR images in mountain areas and different look angles to verify the matching effect of our method under different look angles. The measured data is used to verify the matching effect of our method for different types of terrain SAR images under the same difference in look angles. We compare the performance of our method with those of SAR-SFIT and PSO-SIFT to show that our method is more suitable for SAR images with large geometric distortion matching than these two methods. More specifically, we compare the Mean-Absolute Error (MAE) of matching results of these algorithms, Number of Keypoints Matched (NKM), and Proportion of Keypoints Matched (PKM) to demonstrate the pros and cons of these algorithms.

3.1. Data Set

The simulated SAR data is generated by the Space-borne Radar Advanced Simulator (SRAS) system [36,37]. This batch of data is shown in Figure 8. It is a simulation of four types of mountain terrains. The size is 512 × 512 pixels, and the range and azimuth resolution is about 1 m. From area 1 to area 4, their elevation ranges are 350 m–470 m, 320 m–470 m, 390 m–550 m, and 395 m–570 m, respectively. Their look angles are ${15}^{°}$, ${20}^{°}$, ${25}^{°}$, ${30}^{°}$, and ${40}^{°}$, respectively. The measured data are taken from the TerraSAR-X system, which are L1A-level SAR images collected in the Alps and its vicinity. Master images were collected on 12 January 2009, with a look angle of 35.8153${}^{°}$, and slave images were collected on 9 December 2008, with a look angle of 20.7765${}^{°}$. As shown in Figure 9, we use four terrain image blocks with a size of 512 × 512 pixels.

3.2. Implementation Details

Refer to Dellinger et al. [12] and Ma et al. [22] for SAR-SIFT and PSO-SIFT, respectively. When constructing the scale space, use the initial scale $\delta =2$, ratio coefficient $k=1.26$, and number of scale space layers $N_{max}=8$. The arbitrary parameter d of the SAR-Harris function is set to 0.04, and the threshold is set to 0.8. For RLKD, we set the radius of the search space to 5. For the SAR image after geometric registration, the feature scale and direction in the image are almost the same. Therefore, the standard deviation of the Gaussian function of the algorithm in this paper is set to $\sigma =k^{N_{max}-1}$ for producing large-scale features. In addition, for SAR-SIFT, PSO-SIFT and the method proposed in this paper, the LWM model is set as the default transformation model between the reference and the image.

We tested all the programs on an Ubuntu 18.04 system computer with 128 GB RAM, which is equipped with an Intel i9-9700X CPU and two Nvidia RTX3090 graphics cards.

3.3. Evaluation Index

• Mean-Absolute Error (MAE):
  MAE is capable to measure the alignment error of keypoints, which is defined as follows:

  $$\mathrm{MAE}={\displaystyle \sum _{\left(v_i^m,v_j^s\right)\in \mathfrak{C}}}\frac{{∥v_i^m\mathsf{\Gamma }-v_j^s∥}_2}{\left|\mathfrak{C}\right|}$$

  (14)

  where, $\Gamma$ is the transfer model, and $\left|\mathfrak{C}\right|$ is the number of keypoint pairs that are correctly matched, that is, NKM.
• Number of Keypoints Matched (NKM):
  We use the final number of matching keypoints generated by each method as the number of keypoints matched to measure the effectiveness of the transfer model fitting.
• Proportion of Keypoints Matched (PKM):
  In order to evaluate whether the keypoints detected by the method are efficient, we also use PKM as one of the evaluation indicators. PKM is defined as follows:

  $$\gamma =\frac{\left|V_{\mathrm{matched}}^s\right|}{\left|V^s\right|}$$

  (15)
  In the equation, $\left|V_{\mathrm{matched}}^s\right|$ represents the number of matching keypoints in the master image, and $\left|V^s\right|$ represents the number of all keypoints detected in the master image.

3.4. Result Analysis

In order to verify the performance of the algorithm in this paper, we designed the following experiments. First, in order to confirm the correctness of our choice of measurement function and transformation model in the algorithm, we designed the experiments and presented the results in Table 2 and

Table 3. Differences in the RLKD matching results of our method obtained by using different transformation models.

Method	RLKD		RLKD + MHTIM
	MAE (pixel)	NKM	MAE (pixel)	NKM
Similarity	1.44	11	2.61	20
Polynomial (order 2)	0.63	18	1.19	26
Affine	0.78	16	0.84	18
LWM	0.59	21	0.55	31

. Second, in order to verify the pros and cons of the algorithm compared with other methods, we compared the MAE, NKM and PKM values of the registration results of the four methods on SAR images with different incident angle differences and different terrain undulations in Figure 8, Figure 9, Figure 10, Figure 11, Figure 12 and Figure 13. Then fusion result of our method on real data was showed in Figure 14. The rest of this section will provide a detailed discussion of the results of these experiments.

For the descriptors proposed in this paper, Table 2 shows the effect of using different similarity measurement functions on the fast matching results of RLKD on the simulation data set under a difference of 5${}^{°}$ look angles. In Table 2, MCC is the abbreviation for the maximum value of the correlation function. We can observe that using the correlation function, the proposed method produces the smallest MAE and the largest NKM, but the average time is not increased significantly. Therefore, we choose the correlation function as the measuring function of the descriptor similarity in the RLKD step.

Table 3 shows the influence of different transformation models on the final matching results of our method for the simulated data set with a difference of 5${}^{°}$ in look angles. It can be seen that after RLKD, the LWM model has the smallest MAE and the largest NKM, which are significantly better than those of other models. The similarity model is the simplest fitting model, but its result is the worst. The polynomial and affine models are worse than LWM, and the simplest similarity model is the worst. After MHTIM, increased NKM was produced by all the models. When the LWM model was used alone, the MAE of the matching results decreased. This is because, in the RLKD stage, the LWM model creates more stable matching result. The results presented in Table 3 also confirm our analysis in Section 3, that is, when registering mountain SAR images, the LWM model may achieve better results than other models.

Figure 8 shows the simulated data. In order to give an intuitive experience of the matching effects of different algorithms, and to take into account the clarity of the figure, this article shows only the matching keypoint results of the RLKD method and the SAR-SIFT in the figure. We can observe that compared to SAR-SIFT, the RLKD method achieves more matching points under large difference in look angles. It shows that the descriptors proposed in this paper are more suitable for the matching problem of mountain SAR images having large geometric distortion, which confirms the Section 2.1 of this paper. It also confirms that the ridge structure remains relatively unchanged in SAR images with different look angles.

Figure 10 shows the MAE values of the three registration algorithms for the simulated data under different differences in look angles. Among them, the MAE of SAR-SIFT is the largest, and the results are unstable in areas with different elevation differences. Overall, except for the case where DIA is equal to ${25}^{°}$ in district 2, our method obtains the smallest MAE, and when the DIA increases, the MAE of the RKLD and MHTIM methods remain stable. This result shows that the method proposed in this paper can effectively overcome the geometric distortion caused by the large look angle difference and ensure the stability of the registration process. From the comprehensive analysis of NKM in Figure 11, it can be seen that for PSO-SIFT, the result of MAE is smaller than SAR-SIFT and its NKM is higher. But from the perspective of PKM in Figure 12, only about $10\%$ of the keypoints detected by PSO-SIFT are matched. This value is not higher than SAR-SIFT. This suggests that the keypoints detected by PSO-SIFT may not be more applicable than SAR-SIFT in SAR images with large geometric distortion. The MAE of RLKD in this paper is smaller than those of SAR-SIFT and PSO-SIFT. The RLKD+MHTIM method further reduces the MAE. This shows that the use of the isomorphism of the distribution structure of keypoints increases the stability of our method.

Figure 11 shows the NKM of the matching results of each method for the simulated data. Intuitively, PSO-SIFT is close to RLKD+MHTIM in district 1 and district 4, that is because these two areas have richer textures. The NKM of the RLKD method is between those of SAR-SIFT and PSO-SIFT, but the RLKD+MHTIM method greatly increases the NKM. Comprehensive analysis of Figure 10 and Figure 11 reveals that the RLKD+MHTIM method has the largest number of matching keypoint and the smallest error.

Figure 12 shows the PKM value of the matching results of each method for the simulated data. Except for district 3, the RLKD method produced slightly higher PKM values than SAR-SIFT and PSO-SIFT. However, the PKM value of the RLKD+MHTIM method is significantly higher than those of the other two methods.

In addition, the MAE of the matching result of the correlation-based method is shown in Figure 13. Since this method does not look for keypoints, we present only the MAE results. Compared with the feature-based method, the MAE of this method is at a higher level.

Figure 9 shows the measured TerraSAR-X image data. As in Figure 8, in order to give an intuitive experience of the matching effects of different algorithms and take into account the clarity of the figure, this paper shows only the matching keypoint results of the RLKD and SAR-SIFT methods. Figure 15 shows the histogram of NKM after matching the TerraSAR image data set with three methods. We can see that in the Mountain (Big) area with large terrain undulations, the number of matching keypoints produced by SAR-SIFT is less than that produced by our RLKD method and the distribution of keypoints is uneven. The performance of PSO-SIFT on this type of terrain is the most unstable. In the Mountains (Big) 1, 2 and 4 areas, the number of matching point pairs it obtains is less than one half of those obtained by the SAR-SIFT and RLKD methods. However, on the Mountains (Big) 3 area, the number of matching keypoint reached 39, surpassing SAR-SIFT’s 7 and RLKD + MHTIM’s 33. In the Mountains (Small) area with slightly smaller terrain undulations, NKM obtained by the RLKD method is slightly more than that of SAR-SIFT, and, more than that of PSO-SIFT except in area 3. The MHTIM method further matches the keypoints of the ridge detected by the RLKD method and produces at least twice the number of matching keypoints produced by the RLKD method.

In the areas of towns and others where the terrain is less undulating and more common, the NKM obtained by the RLKD+MHTIM method still has an absolute advantage over those of PSO-SIFT and SAR-SIFT, which proves that the RLKD+MHTIM method proposed in this paper is not only effective and efficient in matching SAR images with large geometric distortion, but also has advantages in registering general types of SAR images.

Finally, Figure 14 shows the matching and fusion results of the TerraSAR-X data based on the RLKD + MHTIM method. The shape of a grid can be seen in every plot because we use “checkerboard” display mode. We can find that in the mountain area, the shape of the ridge in the center of the image fits well. In areas with small terrain undulations, roads, rivers and other terrains on the ground are well matched, which proves that our method has good effectiveness to different types of terrains.

4. Discussion

In this paper, a Ridge Line Keypoint Detection (RLKD) method and a Multi-Hypothesis Topological Isomorphism Matching method (MHTIM) are proposed to match SAR images with large geometric distortion. We designed adequate experiments to verify the effectiveness of the intermediate links and final results of this method. Judging from the experimental results, the LWM and MCC methods are the best choices in the middle part of the method in this paper. Considering the registration results, the method proposed in this paper is more stable than traditional methods when the relative geometric distortion between SAR images increases. Thanks to the inherent isomorphism of the distribution of ridge structures under different viewing angles, the MHTIM method outputs more keypoints and obtains a smaller MAE. Compared with other methods, MHTIM also uses the keypoints detected in RLKD more efficiently.

As shown by experimental results on simulated and real SAR images, the merits of using RLKD and MHTIM are pretty well demonstrated. However, the key points obtained by several types of methods are still unable to obtain a high-precision transformation model to completely correct the SAR images. This shows that the isomorphism between the distribution of ridge features from different look angles has not been completely explored utilized.

5. Conclusions

Aiming at the problem of SAR image registration with large differences in look angles and large terrain undulations, a Ridge Line Keypoint Detection (RLKD) method and a Multi-Hypothesis Topological Isomorphism Matching method (MHTIM) are proposed. First of all, this paper designs a method for detection and similarity description of ridge keypoints. This method can quickly generate a small number of stable matching keypoint pairs under large differences in look angles and large terrain undulations. This method uses the local and global geometric relationship information between keypoints at the same time, therefore, it is more effective than traditional methods when registering SAR images with large geometric distortion. Experiments show that the proposed method can match SAR images with large geometric distortions and can process different types of terrains. It is worth mentioning that before using our method, it is best to perform rough geometric registration first, so that the image regions to be registered roughly overlap. In this way, it will be more efficient when using our method for higher precision registration. In future work, we will conduct in-depth research on the features of ridge images from the electromagnetic model and imaging geometry of the ridge terrain, in order to find a feature description that can cope with the distortion caused by changes in scale, rotation, and imaging look angle.