TIE POINTS EXTRACTION FOR SAR IMAGES BASED ON DIFFERENTIAL CONSTRAINTS

Automatically extracting tie points (TPs) on large-size synthetic aperture radar (SAR) images is still challenging because the efficiency and correct ratio of the image matching need to be improved. This paper proposes an automatic TPs extraction method based on differential constraints for large-size SAR images obtained from approximately parallel tracks, between which the relative geometric distortions are small in azimuth direction and large in range direction. Image pyramids are built firstly, and then corresponding layers of pyramids are matched from the top to the bottom. In the process, the similarity is measured by the normalized cross correlation (NCC) algorithm, which is calculated from a rectangular window with the long side parallel to the azimuth direction. False matches are removed by the differential constrained random sample consensus (DC-RANSAC) algorithm, which appends strong constraints in azimuth direction and weak constraints in range direction. Matching points in the lower pyramid images are predicted with the local bilinear transformation model in range direction. Experiments performed on ENVISAT ASAR and Chinese airborne SAR images validated the efficiency, correct ratio and accuracy of the proposed method.


INTRODUCTION
Synthetic Aperture Radar (SAR) is widely used in rapid terrain mapping of cloudy, foggy and rainy areas because of its full-day, all-weather imaging capabilities (Jin et al., 2015).In commonly used terrain mapping techniques with SAR images, radargrammetry technology is not limited by image coherence and has a practical application value in some specific scenes.However, as a fundamental step of radargrammetry, tie points (TPs) are mainly extracted by the operator of human-computer interaction, which seriously reduces the efficiency of SAR image terrain mapping (Jin et al., 2014).
Automatic extraction of TPs is usually based on image matching technology.However, it is challenging for SAR image matching due to the multiplicative speckle noise and the nonlinear geometric deformations.Broadly, researches on SAR image matching can generally be classified into two types.The first type makes some improvements in the scale-invariant feature transform (SIFT) operator (Fan et al., 2015).SIFT operator is capable of overcome the effects of geometric distortions greatly.Thus, these improvements focus on diminishing speckle influence.Typical methods include the SIFT-OCT method, which skips the features detected at the first octave of the scale space pyramid (Schwind et al., 2010); the BFSIFT method, which builds scale space by a bilateral filter (Wang et al., 2012); the ISIFT method, which transforms the noise model using logarithmic transformation (Suri et al., 2010); and the SAR-SIFT method, which defines a new gradient obtained from the ratio of exponentially weighted averages (ROEWA) algorithm (Dellinger et al., 2015, Zhu et al., 2016).The second type uses the edge-feature-based matching method.Although the extracted edges tend to be discrete and dubious on SAR image, they are more robust compared with the feature points.Common methods of this type include the method based on crossroad and road junction (Dell'Acqua et al., 2004); the method based on the strength and direction of the edge points (Chen et al., 2014, Chen andChen, 2014); and the method based on edge features described by the distance and orientation to the center feature points (Zhang et al., 2015).
Although many of the aforementioned methods are effective to some extent, most of them are not validated by large-size SAR image experiments (Chen et al., 2007).As the image size increases, the efficiency of the two types of methods will be significantly reduced (Chen et al., 2008, Li et al., 2014).Moreover, the relative geometric distortions between large-size SAR images are not globally uniform.Therefore, the interest points may not be clearly distinguishable from one another, which results in mismatches.
To overcome these limitations, this paper mainly focuses on matching large-size SAR images obtained from approximately parallel tracks.This type of image is usually applied in radargrammetry and has obvious relative geometric distortion characteristics.The relative geometric distortions are small in azimuth direction and large in range direction.Taking into account the characteristics, this study proposes an automatic TPs extraction method based on differential constraints.Image pyramids are built firstly, and then corresponding layers of pyramids are matched from the top to the bottom.In the process, some improved technologies including normalized cross correlation (NCC) matching with rectangle window, false matches removal with differential constrained random sample consensus (DC-RANSAC) and matching point prediction with local modeling are presented to weaken the effects of geometry distortions.The efficiency, correct ratio and accuracy of the proposed method are confirmed by experiments on various types of SAR data obtained from platforms at different heights.

PROPOSED METHOD
The flow chart of the proposed method is shown in Fig. 1.The image pyramids are built for the reference and the sensed SAR image, respectively.In detail, the lower pyramid image is blurred by the Gaussian filter, and then 9 (3 × 3) pixels of the image are down sampled as 1 pixel of the upper pyramid image.Moreover, we use Moravec operator to extract interest points on each layer reference image since the proposed method only needs roughly uniformly distributed points on the reference image.Obviously, Moravec operator is simple and efficient, and can obtain the proper distribution and number of interest points with the local non-maximum suppressed.After the pyramids are built, the top layer is matched firstly, and then followed with the layer by layer matching.
Top layer matching: the interest points are extracted on the top layer reference SAR image.For each interest point, the corresponding candidate matching point is searched globally on the sensed image via the NCC with rectangle window.The principle of searching is to achieve the maximum correlation coefficient.After obtaining the candidate matches, the false matches are removed by the DC-RANSAC algorithm, and then obtaining the correct matches of the top layer.
Layer by layer matching: the interest points are extracted on the corresponding layer reference SAR image.The correct matches on the lower pyramid images are applied to establish the global bilinear transformation model in azimuth direction and the local bilinear transformation model in range direction, and then to predict the matching point of each interest point on the sensed image.The search window is centered at the predicted matching point, and the corresponding candidate matching point is searched locally in the search window via the NCC with rectangle window.After obtaining the candidate matches, the false matches are removed by the DC-RANSAC algorithm, and then obtaining the correct matches of the corresponding layer.Following the above steps, the pyramid images are matched layer by layer.Until the original images (pyramid bottom images) are matched, the TPs are obtained.
The key steps involved in the matching process include NCC with rectangle window, false matches removal with DC-RANSAC, matching point prediction with local modeling.

NCC with Rectangle Window
The normalized cross correlation measures the similarity between the interest points of the reference image and the image points of the sensed image by calculating the correlation coefficient of gray values in the neighborhood of two points (Ye et al., 2017).The correlation coefficient γ ranges from -1 to +1, and the closer γ is to +1, the more similar the two points will be.Thus, if the correlation coefficient γ is locally maximum and larger than a certain threshold λ, the image point is the matching point of the corresponding interest point.The neighborhood of the points is called the matching window, which is generally rectangular.The length (azimuth length) and the width (range length) of the window are m and n, respectively.As we known, the relative geometric distortions are small in azimuth direction and large in range direction for two SAR images.In the proposed method, the side of the matching window parallel to azimuth direction is larger than that parallel to range direction, that is m > n.
The correlation coefficient γ is calculated as: Where, M = INT (m/2) , N = INT (n/2) , INT(•) rounds the element of • to the nearest integer towards minus infinity.fi+x,j+y is the intensity value of the reference image at the image point (i + x, j + y) .g i+x ,j+y is the intensity value of the sensed image at the image point (i + x , j + y ).f and ḡ are the means of the intensity values in matching windows for the two image, respectively.

False Matches Removal with DC-RANSAC
RANSAC algorithm is a general parameter estimation approach, which is commonly used in the false matches removal of image matching.The algorithm uses the candidate matches set that contains a large number of false matches to estimate the transformation model of the two images.
The bilinear transformation model can be used in RANSAC: (2) (2) and ( 3) are the transformation model in range and azimuth direction, respectively.where, x1 and x2 are the range coordinates, y1 and y2 are the azimuth coordinates.ai and bi are the parameters of the transformation model.
If ρ and ε are the thresholds in range and azimuth direction, respectively, the criterions for determining the correct matches: Since the relative geometric distortions are small in azimuth direction and large in range direction for two images, the correct matches can accurately fit the bilinear transformation model in azimuth direction, but are difficult to fit the low order polynomial model (such as the bilinear transformation model) in range direction.Therefore, there is a coordinate offset x off between the calculated range coordinate and the correct range coordinate.The threshold ε in ( 5) is a small value, and the threshold ρ in ( 4) is related to the coordinate offset.
The RANSAC algorithm, with the bilinear transformation model (strong constraint model) in azimuth direction and the bilinear transformation model containing the coordinate offset (weak constraint model) in range direction, is called the DC-RANSAC algorithm.The algorithm is calculated as follows: 1. Randomly select a sample (4 matches) in the candidate matches set to initialize the bilinear transform model according to (2) and (3).
2. Determine whether each match in the candidate matches set is correct, and form the correct matches set.
3. If the size of the correct matches set exceeds the threshold τ , re-estimate the model parameters with the set data and end the algorithm.
4. Otherwise, select a new sample, and then repeat steps 1) to 3) K times.Update the correct matches set using the largest correct matches set, and re-estimate the model parameters with the set data.Then end the algorithm.
5. Calculate the coordinate offset of each correct matches in range direction, and get the maximum coordinate offset of all the correct matches x offmax .

Matching Point Prediction with Local modeling
The purpose of predicting the position of the matching point on the sensed image is to improve the matching efficiency.In this step, firstly map the correct matches of the upper pyramid to the lower pyramid.Then, establish the transformation model between the two images to complete the prediction.
The transformation model has been established when remove false matches, namely (2) and (3).However, there is a coordinate offset when predict matching points with (3).In this case, to search for the matching point in the neighborhood of the predicted point, the local search window for the NCC: Where, k is the multiple of the down sampling.ε is the azimuth threshold of DC-RANSAC for the upper pyramid.x offmax is the maximum coordinate offset of all the candidate matches for the upper pyramid, x offmax ≥ x offmax .
Since the range coordinates predicted by (2) are not accurate enough, the local search window represented by ( 6) is oversize, which reduces the matching efficiency.To overcome this problem, the local transformation model in range direction is needed.For each feature point i , 4 correct matches closest to it is searched on the lower reference image.These matches can estimate the parameters of a local bilinear transformation model, that is: The global bilinear transformation model ( 2) and the local bilinear transformation model ( 7) can predict the matching point more accurately.Therefore, we set the local search window:

Test Datasets
To evaluate applicability, three different types of SAR datasets obtained from platforms at different heights are tested.Dataset 1 is the spaceborne SAR data of the Henan Dengfeng mountain The parameters of the three datasets (the reference and sensed images) including the image type, the wave band, the azimuth pixel spacing (Aps), the range pixel spacing (Rps), the azimuth size (As) and the range size (Rs) are shown in table 1.

Experiments and Results
Experiments on three datasets were designed.In the experiments, 4-layer pyramid was built for Dataset 1 and 2, and 5-layer pyramid was built for Dataset 3. In addition, all the following experiments were carried out on a laptop with Intel Core i5 2.40GHz processor and 1 GB RAM.
The threshold of NCC matching λ, the range and azimuth thresholds of DC-RANSAC (ρ and ε) as well as the size of matching window were shown in table 2. The range threshold of DC-RANSAC was obtained by extending the maximum coordinate offset x offmaxj by k times (j was the layer number, j = 1, 2, 3, 4).
The length and width of the each layer matching window were 1.5 times of the upper layer matching window.The size and shape of the top layer matching window was given directly here and would be discussed in section 3.4.
The matching results of the proposed method were analyzed quantitatively and shown in table 3. The correctness of each match was identified by artificially comparing the corresponding points on the reference and the sensed SAR images.The evaluation criterion includes the total time (Tt), the total time by parallel computation (Ttpc), the number of matches (Nm), the number of correct matches by artificial comparison (Ncmac) and the correct ratio by artificial comparison (Crac).The positions of the  In table 3, the correct ratio of Dataset 1 and 2 was 100.0%, and of Dataset 3 was 94.4%, where a false match appeared.Therefore, the correct ratio of the proposed method was high.However, the method could not ensure that all matches were correct, especially for Dataset 3, because the local texture of Dataset 3 was similar in some areas.After multi-thread parallel computing, the total time was reduced from over 5 minutes to about 1 minute indicating that the proposed method could achieve high computational efficiency through parallel computing.
In Fig. 2, the correct matches or TPs were uniformly distributed in the overlapping area because the interest points were extracted with the local non-maximum suppressed.Fig. 3 showed the details of 6 matches.The TP "Tm121" in Fig. 3 (c) was the false match.

Accuracy Analysis
To objectively evaluate the accuracy of the extracted TPs, the stereoscopic orientation experiments were curried out for Dataset 3. In the experiment, we constructed the error equations using the range-Doppler (r-D) model.r-D model is a rigorous geometric model widely used in processing of SAR images (Schmitt et al., 2013).There were 5 ground control points (GCPs) on the reference image and 6 GCPs on the sensed image. 2 of them were conjugate GCPs, which were in overlapping areas between the image pair.5 independent check points (ICPs), named ICP1, ICP2, ICP3, ICP4 and ICP5, were also collected.
The ground coordinates of GCPs and ICPs were measured by differential GPS.The distribution of them was shown in Fig. 4. Subsequently, GCPs, ICPs and extracted TPs were used as input data for the stereoscopic orientation.The root mean square error (RMSE) of ICPs was computed to evaluate the accuracy of the stereoscopic orientation and it could also reflect the accuracy of TPs.
In table 4, the RMSE of ICPs was less than 0.2 m, which was a smaller value.It indicated that the accuracy of the extracted TPs was high, which can meet the needs of SAR image geometric processing.

Discussion on Optimal Matching Window
To determine the size and shape of the NCC matching window, the top layer matching experiments were performed under different size and shape of matching window.In the process, the DC-RANSAC correct ratio of matches P , which was the ratio of the correct matches to the candidate matches, was obtained after the false matches were removed.
Then, we could give the principle of selecting the optimal matching window (OMW).If P was higher than 50%, the matching window with the minimum matching time was OMW.If P was lower than 50% and greater than 25%, the matching window with the maximum correct ratio was OMW.If P was lower than 25%, there was no OMW.
The matching time included the time of NCC matching and DC-RANSAC calculation.Compared to the NCC matching time T , the DC-RANSAC calculation time was expected to be very small, so the matching time here only included T .
From Fig. 5, we could find the optimal square matching windows for the three datasets.For Dataset 1, P could exceed 50%.The optimal square matching window was 15×15 where the minimum matching time happened.For Dataset 2, P did not exceed 50% but greater than 25%.The optimal square matching window was 31×31 where P reaches the maximum value of 40.2%.For Dataset 2, P was less than 25%, so there was no optimal square matching window.From Fig. 6, we could also get the optimal rectangular matching windows for the three datasets, which was 7×23, 23×43 and 7×39, respectively, according to the principle of selecting OMW.
We could compare the optimal square matching window and the optimal rectangular matching window for the three datasets according to Fig. 5 and 6.For Dataset 1, when the matching window was 15×15 or 7×23, P exceeded 50%, and T was 335 seconds and 262 seconds, respectively.Hence, the rectangular window had higher matching efficiency.For Dataset 2, when the matching window was 31×31 or 23×43, P was 40.2% and 43.8%, respectively, and T was 249 seconds and 202 seconds, respectively.Hence, the rectangular window had higher matching accuracy and efficiency.For Dataset 3, there was no optimal square matching window, but there was an optimal rectangular matching window of 7×39.The above comparison showed that the rectangular window was more suitable for matching, and the improvements in the matching accuracy and efficiency were more obvious for Dataset 2 and 3.The reason was that the airborne SAR images used in the experiment were slant images.The relative geometric distortions in azimuth direction were larger than those in range direction.
Experiments on different size and shape of matching windows showed that the top layer matching window for the three datasets should be 7×23, 23×43 and 7×39, respectively.

CONCLUSION
This paper focuses on matching large-size SAR images obtained from approximately parallel tracks to extract TPs.The characteristics of the relative geometric distortion of such images are emphasized.In view of these characteristics, an automatic TPs extraction method based on normalized cross correlation matching with rectangle window, false matches removal with DC-RANSAC and matching point prediction with local modeling is proposed.
The related experiments are performed on spaceborne and airborne, respectively.The proposed approach weakens the effects of geometry distortions.Hence, it can achieve high matching correct ratio and efficiency and accuracy under the optimal matching windows.In addition, it can extract TPs with high accuracy which can meet the needs of SAR image geometric processing.Also, it is suitable for various types of SAR data obtained from platforms at different heights.However, the DC-RANSAC algorithm in this paper does not completely remove false matches, and poor results will be obtained especially for local texture-like SAR images because there is only a strong constraint model in azimuth direction.This will be the subject of further work.

Figure 1 .
Figure 1.Flow chart of the proposed method.

Figure 2 .
Figure 2. Matching results for three datasets.

Figure 3 .
Figure 3. Six matching points on enlarged pictures for three datasets.

Figure 4 .
Figure 4.The distribution of GCPs and ICPs for Dataset 3.

Figure 5 .
Figure 5. Top layer matching with square matching window.

Figure 6 .
Figure 6.Top layer matching with rectangular matching window.

Table 1 .
Parameters of three datasets.
area obtained by European Space Agency ENVISAT ASAR system.Dataset 2 is the airborne SAR data of the Shanxi Weinan hilly area obtained by the airborne SAR System of Chinese Academy of Surveying and Mapping (CASM).Dataset 3 is the airborne SAR data of the Shanxi Yanliang hilly area obtained by the airborne SAR system of 23rd Institute, China Aerospace Science & Industry Corp.

Table 3 .
Matching results for three datasets.

Table 4 .
Errors and RMSE of ICPs for Dataset 3.