An Image Matching Method for SAR Orthophotos from Adjacent Orbits in Large Area Based on SAR-Moravec

: In producing orthophoto mosaic in a large area from spaceborne synthetic aperture radar (SAR) images, SAR image matching from adjacent orbits is a technical di ﬃ culty due to the speckle noise and di ﬀ erent imaging mechanism between azimuth and range direction. In this paper, an area-based method, SAR-Moravec, is proposed for SAR orthophoto matching from adjacent orbits in a large area. Compared with the classical area-based Moravec, the template of SAR-Moravec is characterized by more directions for speckle noise restraint and a main direction consistent with the azimuth. In order to get evenly distributed matching points with high accuracy, the grid control mechanism and Gaussian pyramid from coarse to ﬁne are introduced in matching. The whole process contains three steps. First, the pyramid images are constructed by the down-sampling process. Second, under grid control, the feature points are evenly extracted by the modiﬁed template. Third, the transformation model is iteratively calculated from the ﬁrst to the last layer of the pyramid. After the matching process layer-by-layer, the ﬁnal matching points and transformation model can be obtained. In the experiments, we compare the SAR-Moravec with three widely used methods, including the Moravec, the SAR-scale invariant feature transform (SAR-SIFT), and the SAR-features from an accelerated segment test (SAR-FAST). The results indicate that the proposed method has the best global matching accuracy among these methods and the matching e ﬃ ciency is better than SAR-SIFT and SAR-FAST methods in large area.


Introduction
Due to the insensitivity to illumination and atmospheric conditions, synthetic aperture radar (SAR) provides an important data source for numerous research studies in all-weather types and all time points. The accessibility of spaceborne SAR images has been greatly improved for the growing number of SAR satellite missions launched in this century [1][2][3][4][5]. Rich data make it possible to mosaic a large and continuous SAR orthophoto for large scale applications [6], such as large area mapping [7], railway survey [8], and resource management [9]. However, forming such a large SAR orthophoto is not an easy task. The difficulties are mainly shown in three aspects. First, the incidence angles of adjacent orbits are quite different, which leads to different geometric distortion in the same area of two adjacent images. Second, the large area is usually covered by different types of landform. In some areas with less human activity, such as a forest, desert, and alpine regions, there are few ground control 2 of 19 points (GCPs) provided [10][11][12]. At the same time, the corresponding points for image matching are also not easy to acquire for the image coherence, which, in those areas, is relatively low. Third, to accurately geolocate the whole images in coordinate system, evenly distributed matching points with good accuracy are needed for GCPs densification in a follow-up block adjustment [13]. Moreover, a continuous matching process accumulate the error to a serious extent in areas where GCPs are widely lacked. For example, in a mountain area of Iran and Pakistan, shown in Figure 1, three images from adjacent orbits should be matching and mosaicking without GCPs. Therefore, to automatically and efficiently produce such a large SAR orthophoto, a SAR image matching method that is invulnerable to terrain types and capable of generating evenly distributed matching points is needed.
Remote Sens. 2020, 12, x FOR PEER REVIEW 2 of 19 control points (GCPs) provided [10][11][12]. At the same time, the corresponding points for image matching are also not easy to acquire for the image coherence, which, in those areas, is relatively low. Third, to accurately geolocate the whole images in coordinate system, evenly distributed matching points with good accuracy are needed for GCPs densification in a follow-up block adjustment [13]. Moreover, a continuous matching process accumulate the error to a serious extent in areas where GCPs are widely lacked. For example, in a mountain area of Iran and Pakistan, shown in Figure 1, three images from adjacent orbits should be matching and mosaicking without GCPs. Therefore, to automatically and efficiently produce such a large SAR orthophoto, a SAR image matching method that is invulnerable to terrain types and capable of generating evenly distributed matching points is needed. Unlike the optic image, SAR orthophoto has very different properties. The side-looking mechanism causes the difference of incidence angle from the near side to the far side, especially in large-width mode spaceborne SAR images [14][15][16]. Furthermore, after terrain correction, the azimuth and range directions are no longer consistent with the row and column directions in SAR orthophoto [15]. As a result, the quality of feature points and matching accuracy would be limited if the methods designed for optic image matching are directly applied to SAR [11,16]. Therefore, a necessary improvement should be made on those methods.
The methods of image matching can be generally classified into two types: feature-based methods and area-based methods [17]. Feature-based methods like scale invariant feature transform (SIFT) [18] and speeded up robust feature (SURF) [19] have been widely applied in image matching [20]. These methods can present convincing results in optic image matching but still have limits on the processing SAR image. One improved version for SAR image matching is proposed by Dellinger et al. [21], namely the SAR-scale invariant feature transform (SAR-SIFT). It relies on a new gradient computation method to decrease the influence of speckle noise and incidence angle [22]. In recent years, many improvements have been made on SAR-SIFT. Zhu et al. [23] proposed a novel SAR image registration method with SAR-SIFT for feature detection and the arborescence network matching (ANM) for feature matching. Dubois et al. applied SAR-SIFT on stereo-SAR images and combined it with the SAR-stereo model for extracting reliable tie points in TerraSAR-X HR SpotLight data. Wang et al. [22] proposed a uniform method for SAR image registration based on SAR-SIFT for generating enough robust, reliable, and evenly distributed features with the Voronoi diagram and feature scalespace proportional extraction. Paul et al. [24] modified SAR-SIFT and proposed the improved SAR-SIFT (I-SAR-SIFT) method with Delaunay-triangulation to enhance the matching performance in SAR images. However, since the SAR-SIFT lacks the controllability of feature number and spatial Unlike the optic image, SAR orthophoto has very different properties. The side-looking mechanism causes the difference of incidence angle from the near side to the far side, especially in large-width mode spaceborne SAR images [14][15][16]. Furthermore, after terrain correction, the azimuth and range directions are no longer consistent with the row and column directions in SAR orthophoto [15]. As a result, the quality of feature points and matching accuracy would be limited if the methods designed for optic image matching are directly applied to SAR [11,16]. Therefore, a necessary improvement should be made on those methods.
The methods of image matching can be generally classified into two types: feature-based methods and area-based methods [17]. Feature-based methods like scale invariant feature transform (SIFT) [18] and speeded up robust feature (SURF) [19] have been widely applied in image matching [20]. These methods can present convincing results in optic image matching but still have limits on the processing SAR image. One improved version for SAR image matching is proposed by Dellinger et al. [21], namely the SAR-scale invariant feature transform (SAR-SIFT). It relies on a new gradient computation method to decrease the influence of speckle noise and incidence angle [22]. In recent years, many improvements have been made on SAR-SIFT. Zhu et al. [23] proposed a novel SAR image registration method with SAR-SIFT for feature detection and the arborescence network matching (ANM) for feature matching. Dubois et al. applied SAR-SIFT on stereo-SAR images and combined it with the SAR-stereo model for extracting reliable tie points in TerraSAR-X HR SpotLight data. Wang et al. [22] proposed a uniform method for SAR image registration based on SAR-SIFT for generating enough robust, reliable, and evenly distributed features with the Voronoi diagram and feature scale-space proportional extraction. Paul et al. [24] modified SAR-SIFT and proposed the improved SAR-SIFT (I-SAR-SIFT) method with Delaunay-triangulation to enhance the matching performance in SAR images. However, since the SAR-SIFT lacks the controllability of feature Remote Sens. 2020, 12, 2892 3 of 19 number and spatial distribution, related experimental results of these research studies show that the matching points are concentrated in areas with rich texture information [19,21]. In a large-area matching process, the uneven distribution of matching points usually leads to an unsatisfactory global accuracy. Like other feature-based methods, the implementation of SAR-SIFT is very constrained by computational efficiency.
Compared with feature-based methods, the area-based methods do not need to construct the scale-space and feature descriptor, but simply use the pixel grayscale feature to optimize a similarity measure between two images. Acting as feature point detectors, some area-based methods can also participate in a feature-based matching process, such as Harris [25][26][27], Shi-Tomasi [28,29], and the features from an accelerated segment test (FAST) [30][31][32][33][34]. Usually, area-based methods have good computational efficiency in the optical image. However, the special imaging mechanism of SAR weakens their ability in the matching process. The SAR-features from accelerated segment test (SAR-FAST) is a modified area-based method, which utilizes the rolling guidance filter to reduce the impact of speckle noise. The dissimilarity between the detection windows and the center window is also introduced to locate the feature points [32,33]. As a new area-based method in the SAR image matching process, many research studies have focused on SAR-FAST. Jiao and Kang [35] reduced the impact of speckle noise by replacing one pixel with a Gaussian template and utilizing a ratio-based measure in GF-3 data. Yu and Yang [32] presented a coarse to the fine registration method for airborne SAR images by combining the SAR-FAST detector to the domain-size pooling of learned arrangements of a three patch (DSP-LATCH) descriptor. It makes the process only rely on the grayscale information of SAR data. Based on SAR-FAST, Xiong et al. [36] proposed the FAST-8 corner detector to suppress the points affected by speckle noise with the similarity between the candidates and surrounding pixels. Nevertheless, SAR-FAST is not appropriate for areas with weak texture information and the problem of direction inconsistency after terrain correcting still exists. Although the distribution performance and the restraint of speckle noise are improved, SAR-FAST still have limits on orthophoto matching from adjacent orbits.
In this paper, we tend to find an efficient matching method for SAR orthophotos from adjacent orbits. Meanwhile, the global matching accuracy should be promoted for follow-up tie points extraction in a large area. Widely used in detecting corner points, Moravec [5] is a classical area-based method, which takes corner points as the intensity-based feature points in image matching. Considering its well performance on efficiency, we propose a modified area-based method for SAR orthophoto matching, known as the SAR-Moravec. To overcome the disadvantages analyzed above, several improvements are made in the extraction template and matching procedure. In feature points extraction, more detecting directions are provided to reduce the influence of speckle noise. Furthermore, a rotated template with the main direction is applied in solving the problem of direction inconsistency. Due to the changes on the template, the quality of extracted feature points is improved. In a matching process, the whole procedure contains three steps. First, the Gaussian pyramid is constructed by the down-sampling process. Second, the feature points are extracted under grid control. Third, the transformation model is iteratively calculated from the first layer to the last and the traversal area of matching is iteratively predicted per layer. Due to the changes on the matching procedure, the even distribution of matching points in a large area can be guaranteed, and the transformation model is less vulnerable to terrain types.
In this paper, the improved extraction template and matching procedure are presented in Section 3 after the analysis and explanation of the study area and data selection in Section 2. The detailed experimental results and discussion are presented in Section 4. Lastly, Section 5 concludes the whole paper.

Study Area and Data Selection
In this paper, four pairs of SAR images are used in total, including Radarsat-2 Extra Fine single look complex (SLC) data and Sentinel-1 Interferometric Wide Swath SLC data. information of the data in detail. These image pairs are located in different places with 24 different acquisition types and different spatial resolutions. These data contain different terrains and complex texture information. Therefore, the data are effective for method comparison. To investigate the matching performance in different terrains, we choose 12 pairs of slices from the Sentinel-1 data, sized in 500 * 500 and grouped into four terrain types including desert, mountain, plain, and residential areas. Desert and plain areas represent terrains with weak texture information. Mountain and residential areas represent terrains with rich texture information. The Extra Fine acquisition mode (5 m spatial resolution and 125 km * 125 km swath size) is available with four product levels in Radarsat-2 [2,3]. Selective single polarization and incidence angle from 22 • to 49 • is available. The interferometric wide swath mode with high resolution (2.7 m * 22 m to 3.5 m * 22 m) is a wide swath product in Sentinel-1A and Sentinel-1B [4,37]. This product contains one image per sub-swath. Each image comprises a series of bursts and each burst was processed as a separate SLC image [38,39]. The incidence angle changes from 29.1 • to 46.0 • and its swath width is 250 km.
This paper focuses on SAR orthophoto matching from adjacent orbits in a large area, which covers various terrain types. Different complexity of the texture information in these types influences the extraction of feature points and distribution of matching points in area-based methods. Meanwhile, in areas without human activity, the lack of GCPs increases the difficulty of high-accuracy orthophoto production in a large area. Related characteristics in the data are presented in Figure 2. In this paper, the proposed method is designed for feature points extraction and SAR orthophoto matching in a large area. Section 3 presents the method in detail. and complex texture information. Therefore, the data are effective for method comparison. To investigate the matching performance in different terrains, we choose 12 pairs of slices from the Sentinel-1 data, sized in 500 * 500 and grouped into four terrain types including desert, mountain, plain, and residential areas. Desert and plain areas represent terrains with weak texture information. Mountain and residential areas represent terrains with rich texture information. The Extra Fine acquisition mode (5 m spatial resolution and 125 km * 125 km swath size) is available with four product levels in Radarsat-2 [2,3]. Selective single polarization and incidence angle from 22° to 49° is available. The interferometric wide swath mode with high resolution (2.7 m * 22 m to 3.5 m * 22 m) is a wide swath product in Sentinel-1A and Sentinel-1B [4,37]. This product contains one image per sub-swath. Each image comprises a series of bursts and each burst was processed as a separate SLC image [38,39]. The incidence angle changes from 29.1° to 46.0° and its swath width is 250 km.
This paper focuses on SAR orthophoto matching from adjacent orbits in a large area, which covers various terrain types. Different complexity of the texture information in these types influences the extraction of feature points and distribution of matching points in area-based methods. Meanwhile, in areas without human activity, the lack of GCPs increases the difficulty of highaccuracy orthophoto production in a large area. Related characteristics in the data are presented in Figure 2. In this paper, the proposed method is designed for feature points extraction and SAR orthophoto matching in a large area. Section 3 presents the method in detail.

Moravec and SAR-Moravec
Being one of the most common methods for feature points extraction, Moravec calculates the grayscale variance of adjacent pixels along four directions and compares the minimum with a specific Remote Sens. 2020, 12, 2892 5 of 19 threshold. Figure 3 presents the template of traditional Moravec in five pixels with four directions of detection. The variance of grayscale value can be calculated as follows.
where V denotes the grayscale variance of each direction, g presents the grayscale value of the image pixel, and m and n are a specific value decided in different directions.
Remote Sens. 2020, 12, x FOR PEER REVIEW 5 of 19 Being one of the most common methods for feature points extraction, Moravec calculates the grayscale variance of adjacent pixels along four directions and compares the minimum with a specific threshold. Figure 3 presents the template of traditional Moravec in five pixels with four directions of detection. The variance of grayscale value can be calculated as follows.
where denotes the grayscale variance of each direction, presents the grayscale value of the image pixel, and and are a specific value decided in different directions. There are two limits if the Moravec is directly applied to SAR orthophoto matching. First, Moravec only computes the grayscale information in four directions, which is insufficient to accurately determine the minimum grayscale variance of the pixel. Second, due to the side-looking imaging mechanism, the process of range direction is very different from the azimuth direction [40], and the Moravec is not qualified to deal with the difference. Therefore, in order to resolve these limits and boost the matching accuracy, based on Moravec, a modified SAR-Moravec template is proposed. First, four more directions are added to improve the ability of obtaining the minimum grayscale variance and decrease the influence of speckle noise. Second, since the characteristics in azimuth and range directions are different, we rotate the traditional template to make sure the main direction is consistent with the azimuth direction. The main direction is different from the average directions in its length and the size of the template is determined by specific data in experiments.
In the modified SAR-Moravec template, the direction lines after rotation are not necessarily across the center of other pixels. To calculate the grayscale variance, the simulated pixels in a red dotted frame (Figure 3c) are getting involved. Their grayscale values are the interpolation of overlapped pixels. The rotation angle can be calculated by location vectors of the satellite in image header files like XML file, and the formula is as follows.
where ( , ) denotes the coordinate of location vectors of the satellite, ( ̅ , ) is the mean value of these vectors, is the number of them, and presents the rotation angle in the template. As the angle is defined, the interpolation process can be performed. According to the positions of the pixel center and direction line, real or simulated pixels are confirmed to participate in the calculation of grayscale variance. The formula of the interpolation process is as follows. There are two limits if the Moravec is directly applied to SAR orthophoto matching. First, Moravec only computes the grayscale information in four directions, which is insufficient to accurately determine the minimum grayscale variance of the pixel. Second, due to the side-looking imaging mechanism, the process of range direction is very different from the azimuth direction [40], and the Moravec is not qualified to deal with the difference. Therefore, in order to resolve these limits and boost the matching accuracy, based on Moravec, a modified SAR-Moravec template is proposed. First, four more directions are added to improve the ability of obtaining the minimum grayscale variance and decrease the influence of speckle noise. Second, since the characteristics in azimuth and range directions are different, we rotate the traditional template to make sure the main direction is consistent with the azimuth direction. The main direction is different from the average directions in its length and the size of the template is determined by specific data in experiments.
In the modified SAR-Moravec template, the direction lines after rotation are not necessarily across the center of other pixels. To calculate the grayscale variance, the simulated pixels in a red dotted frame (Figure 3c) are getting involved. Their grayscale values are the interpolation of overlapped pixels. The rotation angle can be calculated by location vectors of the satellite in image header files like XML file, and the formula is as follows.
where (x i , y i ) denotes the coordinate of location vectors of the satellite, (x, y) is the mean value of these vectors, n is the number of them, and α presents the rotation angle in the template. As the angle is defined, the interpolation process can be performed. According to the positions of the pixel center and direction line, real or simulated pixels are confirmed to participate in the calculation of grayscale variance. The formula of the interpolation process is as follows.
where g denotes the grayscale value of real pixels overlapped, p is the weight of them, g is the grayscale value of simulated pixels, and n presents the number of real pixels overlapped.

Template Determining
The template size is one of the most critical parameters in this method. Generally, the main direction to average direction ratio of the template is inversely proportional to the ratio between the azimuth and range resolution. In the template determination of Radarsat-2 Extra Fine orthophoto, the average incidence angle is 36.5 • . After calculating by using the incidence angle, the ground range resolution is 4.650 m, the azimuth resolution is 2.9 m, and the ratio of them is 1.603. Our final choice (7 * 11) is the result of compromise between the closeness to resolution ratio and the computing cost.
To verify the effectiveness of the choice, three kinds of average directions are used (5, 7, and 9 pixels) in the following experiment, according to the resolution of Radarsat-2 [41]. Twenty-four templates with a different main direction size are tested, and the results are shown in Figure 4a,b. The 7 * 11 template size shows the best global matching accuracy and relatively good performance in computing time. The template size of Sentinel-1 Interferometric Wide Swath orthophoto is determined similarly. The average incidence angle is 38.2 • . Due to the process of multi-looking, the ground range resolution is 18.555 m, the azimuth resolution is 14.1 m, and the ratio of them is 1.316. Three kinds of average directions are used (3, 5, and 7 pixels), according to the resolution of Sentinel-1. Another 24 templates with a different main direction size are tested to process the slice of Data 3. The final choice is 5 * 7 and related experimental results are presented in Figure 4c,d.
Remote Sens. 2020, 12, x FOR PEER REVIEW 6 of 19 where denotes the grayscale value of real pixels overlapped, is the weight of them, is the grayscale value of simulated pixels, and presents the number of real pixels overlapped.

Template Determining
The template size is one of the most critical parameters in this method. Generally, the main direction to average direction ratio of the template is inversely proportional to the ratio between the azimuth and range resolution. In the template determination of Radarsat-2 Extra Fine orthophoto, the average incidence angle is 36.5°. After calculating by using the incidence angle, the ground range resolution is 4.650 m, the azimuth resolution is 2.9 m, and the ratio of them is 1.603. Our final choice (7 * 11) is the result of compromise between the closeness to resolution ratio and the computing cost.
To verify the effectiveness of the choice, three kinds of average directions are used (5, 7, and 9 pixels) in the following experiment, according to the resolution of Radarsat-2 [41]. Twenty-four templates with a different main direction size are tested, and the results are shown in Figure 4a,b. The 7 * 11 template size shows the best global matching accuracy and relatively good performance in computing time. The template size of Sentinel-1 Interferometric Wide Swath orthophoto is determined similarly. The average incidence angle is 38.2°. Due to the process of multi-looking, the ground range resolution is 18.555 m, the azimuth resolution is 14.1 m, and the ratio of them is 1.316. Three kinds of average directions are used (3, 5, and 7 pixels), according to the resolution of Sentinel-1. Another 24 templates with a different main direction size are tested to process the slice of Data 3. The final choice is 5 * 7 and related experimental results are presented in Figure 4c,d.

Matching Procedure
Usually, when the matching process is undertaken, a specific evaluation index is set up as a measurement in feature points extraction at first [42]. Then, with the similarity measure calculated, the matching points can be obtained from these feature points. Lastly, due to the random sample consensus (RANSAC) filtering the outliers, the correct matching points are reserved to construct a stable transformation model of image pairs [43].
However, SAR orthophoto has very different properties from optic images. The specific properties should be considered during matching. Therefore, after the optimization of the template for feature points extraction, we have made the following modifications in the matching procedure. It is a matching strategy from coarse to fine, which utilized the Gaussian pyramid [44] to improve the accuracy of the transformation model per layer. At the same time, the grid control is applied in feature points extraction for evenly distributed matching points in a large area. The matching procedure in detail is illustrated in Figure 5. The whole matching procedure contains three steps. First, the Gaussian pyramid images are constructed by the down-sampling process. Second, under grid control, the feature points are extracted by the modified template. Third, the transformation model is iteratively calculated from the first to the last layer of the pyramid and the traversal areas of the next matching layer are iteratively predicted by the model. After the matching process layer-by-layer, the final matching points and transformation model can be obtained.

Matching Procedure
Usually, when the matching process is undertaken, a specific evaluation index is set up as a measurement in feature points extraction at first [42]. Then, with the similarity measure calculated, the matching points can be obtained from these feature points. Lastly, due to the random sample consensus (RANSAC) filtering the outliers, the correct matching points are reserved to construct a stable transformation model of image pairs [43].
However, SAR orthophoto has very different properties from optic images. The specific properties should be considered during matching. Therefore, after the optimization of the template for feature points extraction, we have made the following modifications in the matching procedure. It is a matching strategy from coarse to fine, which utilized the Gaussian pyramid [44] to improve the accuracy of the transformation model per layer. At the same time, the grid control is applied in feature points extraction for evenly distributed matching points in a large area. The matching procedure in detail is illustrated in Figure 5. The whole matching procedure contains three steps. First, the Gaussian pyramid images are constructed by the down-sampling process. Second, under grid control, the feature points are extracted by the modified template. Third, the transformation model is iteratively calculated from the first to the last layer of the pyramid and the traversal areas of the next matching layer are iteratively predicted by the model. After the matching process layer-by-layer, the final matching points and transformation model can be obtained.
Based on the traditional matching process in area-based methods [42,45], in this procedure, three modifications are emphatically proposed. First, the original transformation model can be calculated by coarse matching, which is used for matching points prediction in a Gaussian pyramid. After down-sampling, the top layer of the pyramid has the lowest resolution and the simplest texture information. The global traversal is suitable to process the top layer image. Second, grid control for image matching per layer can relieve the influence of terrain types and improve the performance of matching points distribution [45,46]. Then, the global accuracy of matching in a large area can be guaranteed. Third, from a coarse type to a fine type, the maximum deviation of matching points and predicted points in the current layer is applied as the prediction range of fine matching in the next layer. This process improves the traversal efficiency with reduced and accurate location after prediction. Based on the traditional matching process in area-based methods [42,45], in this procedure, three modifications are emphatically proposed. First, the original transformation model can be calculated by coarse matching, which is used for matching points prediction in a Gaussian pyramid. After down-sampling, the top layer of the pyramid has the lowest resolution and the simplest texture information. The global traversal is suitable to process the top layer image. Second, grid control for image matching per layer can relieve the influence of terrain types and improve the performance of matching points distribution [45,46]. Then, the global accuracy of matching in a large area can be guaranteed. Third, from a coarse type to a fine type, the maximum deviation of matching points and predicted points in the current layer is applied as the prediction range of fine matching in the next layer. This process improves the traversal efficiency with reduced and accurate location after prediction.

Evaluation Criteria of Matching
In this paper, the evaluation criteria of matching results are classified as matching efficiency, matching accuracy, the correct ratio of matching points, and the distribution of matching points. With these evaluation criteria calculated and analyzed, we can make an intuitive comparison on SAR-Moravec and other methods (Moravec, SAR-FAST, and SAR-SIFT) of their matching results. The criteria are defined as follows.

Evaluation Criteria of Matching
In this paper, the evaluation criteria of matching results are classified as matching efficiency, matching accuracy, the correct ratio of matching points, and the distribution of matching points. With these evaluation criteria calculated and analyzed, we can make an intuitive comparison on SAR-Moravec and other methods (Moravec, SAR-FAST, and SAR-SIFT) of their matching results. The criteria are defined as follows.

1.
Matching efficiency: Time consumption of the whole procedure is used to indicate the matching efficiency of each method.

2.
Matching accuracy: To quantitatively evaluate the accuracy of matching points, the standard deviation (SD) and root mean square error (RMSE) are formulated. SD can be calculated for evaluating the average pixel error. RMSE can be calculated for measuring the alignment of all matching points. These two criteria are formulated as follows.
In Equation (4), (x m,i , y m,i ) represents the location of matching point i and x p,i , y p,i represents the predicted location of the matching point i. The predicted location is calculated by the transformation model. d i is the matching error value of point i. In Equations (5) and (6), d denotes the average error value of all matching points. N represents the number of matching points.

3.
Correct ratio of matching points: The extraction of feature points is the first step in area-based matching and feature-based matching. The quality of them is one important factor of the correct ratio of matching points. To utilize this criterion, we calculate the correct matching ratio (CMR) [22] as follows. CMR = N corr/N orig Remote Sens. 2020, 12, 2892 9 of 19 where N corr denotes the number of correct matching points after the process of RANSAC, which we called inliers. N orig represents the number of all matching points before RANSAC, which includes false matching points (outliers).

4.
Distribution of matching points: Due to the existence of different terrain types in wide swath SAR orthophotos, the distribution of matching points is extremely impacted by the texture information. In the matching procedure of SAR-Moravec, we design the regular grid to optimize the distribution of matching points. The grid control in the SAR orthophoto can be utilized in a large area. The global distribution density (GDD) can evaluate the performance of distribution, which is formulated as follows.
where G i represents the number of matching points in the regular grid i, and G denotes the mean number of matching points in all grids. N is the number of regular grids in experimental data.

Feature Points Extraction of Moravec and SAR-Moravec
The quality of feature points can seriously influence the matching results [47,48]. Ordinarily, the qualified feature points with high grayscale variance are closer to the edge of objects. In feature-based methods, the strong feature descriptor is often located in a sharp gradient area [48]. Calculated by area-based feature points extraction of Moravec and SAR-Moravec, the grayscale variance of feature points shows a higher value than other points. To compare the quality of features points extracted by Moravec and SAR-Moravec, Data 3 are processed and the detailed extraction results are zoomed in for clear presentation. Figure 6  In Equation (4), ( , , , ) represents the location of matching point and ( , , , ) represents the predicted location of the matching point . The predicted location is calculated by the transformation model. is the matching error value of point . In Equations (5) and (6), ̅ denotes the average error value of all matching points. represents the number of matching points.
3. Correct ratio of matching points: The extraction of feature points is the first step in areabased matching and feature-based matching. The quality of them is one important factor of the correct ratio of matching points. To utilize this criterion, we calculate the correct matching ratio (CMR) [22] as follows.
where denotes the number of correct matching points after the process of RANSAC, which we called inliers.
represents the number of all matching points before RANSAC, which includes false matching points (outliers).
4. Distribution of matching points: Due to the existence of different terrain types in wide swath SAR orthophotos, the distribution of matching points is extremely impacted by the texture information. In the matching procedure of SAR-Moravec, we design the regular grid to optimize the distribution of matching points. The grid control in the SAR orthophoto can be utilized in a large area. The global distribution density (GDD) can evaluate the performance of distribution, which is formulated as follows.
where represents the number of matching points in the regular grid , and ̅ denotes the mean number of matching points in all grids.
is the number of regular grids in experimental data.

Feature Points Extraction of Moravec and SAR-Moravec
The quality of feature points can seriously influence the matching results [47,48]. Ordinarily, the qualified feature points with high grayscale variance are closer to the edge of objects. In feature-based methods, the strong feature descriptor is often located in a sharp gradient area [48]. Calculated by area-based feature points extraction of Moravec and SAR-Moravec, the grayscale variance of feature points shows a higher value than other points. To compare the quality of features points extracted by Moravec and SAR-Moravec, Data 3 are processed and the detailed extraction results are zoomed in for clear presentation. Figure 6 shows the extraction results of two slices in Data 3 by Moravec and SAR-Moravec. The detailed presentation should be analyzed to confirm the improvement of SAR-Moravec feature points extraction.

Experiment of Adaptability Compared with SAR-SIFT
To compare the performance of SAR-Moravec with feature-based methods in different terrain types, we process four pairs of slices in Data 1 (200 * 200) by SAR-Moravec and SAR-SIFT. Table 2 presents the matching results with the evaluation criteria above (Matches, SD, RMSE, CMR, and Time consuming). The matching results of these two methods are illustrated in Figure 7.

Experiment of Adaptability Compared with SAR-SIFT
To compare the performance of SAR-Moravec with feature-based methods in different terrain types, we process four pairs of slices in Data 1 (200 * 200) by SAR-Moravec and SAR-SIFT. Table 2 presents the matching results with the evaluation criteria above (Matches, SD, RMSE, CMR, and Time consuming). The matching results of these two methods are illustrated in Figure 7.
Seen from the results of Figure 7, the matching points extracted by SAR-Moravec distribute more evenly than SAR-SIFT. The grid control process contributes good matching results in each terrain type of the slices. According to the evaluation criteria in Table 2, SAR-Moravec can extract more matching points than SAR-SIFT regardless of terrain types. Furthermore, SAR-Moravec shows higher SD and RMSE, even in areas with weak texture information. In the residential area, the improvement is not clear enough. Both methods can reach high-accuracy results in the rich texture information area. However, the SAR-Moravec has higher CMR than SAR-SIFT, which means the SAR-Moravec has better ability to obtain matching point. In time consuming, SAR-Moravec performs significantly better than SAR-SIFT. Therefore, SAR-Moravec is more effective and efficient than feature-based SAR-SIFT in SAR orthophoto matching.

Experiment of Adaptability Compared with SAR-SIFT
To compare the performance of SAR-Moravec with feature-based methods in different terrain types, we process four pairs of slices in Data 1 (200 * 200) by SAR-Moravec and SAR-SIFT. Table 2 presents the matching results with the evaluation criteria above (Matches, SD, RMSE, CMR, and Time consuming). The matching results of these two methods are illustrated in Figure 7.  Seen from the results of Figure 7, the matching points extracted by SAR-Moravec distribute more evenly than SAR-SIFT. The grid control process contributes good matching results in each terrain type of the slices. According to the evaluation criteria in Table 2, SAR-Moravec can extract more matching points than SAR-SIFT regardless of terrain types. Furthermore, SAR-Moravec shows higher SD and RMSE, even in areas with weak texture information. In the residential area, the improvement is not clear enough. Both methods can reach high-accuracy results in the rich texture information area. However, the SAR-Moravec has higher CMR than SAR-SIFT, which means the SAR-Moravec has better ability to obtain matching point. In time consuming, SAR-Moravec performs significantly better than SAR-SIFT. Therefore, SAR-Moravec is more effective and efficient than feature-based SAR-SIFT in SAR orthophoto matching.

Experiment of Adaptability in Different Terrain Types
To verify the SAR-Moravec method applicable to SAR orthophoto matching in different terrain types, slices of representative terrains from Sentienl-1 are prepared for this experiment in Section 2. Twelve pairs of slices classified into four terrain types are processed in SAR-Moravec matching. Four

Experiment of Adaptability in Different Terrain Types
To verify the SAR-Moravec method applicable to SAR orthophoto matching in different terrain types, slices of representative terrains from Sentienl-1 are prepared for this experiment in Section 2. Twelve pairs of slices classified into four terrain types are processed in SAR-Moravec matching. Four evaluation criteria above, including Matches, SD, RMSE, and CMR, can be utilized in analyzing the matching results of different terrain types. The matching results are illustrated in Figure 8 and the detailed indicators are presented in Table 3.
Illustrated in Figure 8, the number of matching points (Matches) in desert areas and plain areas are relatively low for the weak texture information, especially in desert areas. In mountain areas, the change of elevation makes the texture information complicated, which increases the number of matching points. In resident areas, the existence of numerous artificial features with rich texture information also increases the number of matching points. Presented in Table 3 and Figure 8, the number and distribution of matching points can impact the matching criteria of SD, RMSE, and CMR. With good feature points extracted in SAR-Moravec, the matching accuracy and CMR can be promoted. Furthermore, because of the shadow and layover phenomena in mountain areas [49], the matching accuracy (SD, RMSE) and CMR are both lower than resident areas. Due to the weak texture information, less feature points are extracted in desert and plain areas. It leads to less matching points in these slices. However, the matching accuracy of different slices all reach a subpixel level. Therefore, SAR-Moravec can deal well with different terrain types and is suitable to the matching procedure of complex-terrain SAR orthophotos in a large area.
Remote Sens. 2020, 12, x FOR PEER REVIEW 12 of 19 evaluation criteria above, including Matches, SD, RMSE, and CMR, can be utilized in analyzing the matching results of different terrain types. The matching results are illustrated in Figure 8 and the detailed indicators are presented in Table 3.

Experiment of SAR Orthophoto Matching
In this part, the whole swathes of SAR orthophotos from adjacent orbits are used. Figure 9 illustrates parts of the matching results by SAR-Moravec. Table 4 presents the matching results of four data by methods including Moravec, SAR-FAST, SAR-SIFT, and SAR-Moravec. With evaluation criteria of SD, RMSE, CMR, GDD, and time consuming, the different matching results can be analyzed in detail. GDD in feature-based matching methods is obtained by the same regular grids utilized in area-based methods.
Remote Sens. 2020, 12, x FOR PEER REVIEW 13 of 19 Illustrated in Figure 8, the number of matching points (Matches) in desert areas and plain areas are relatively low for the weak texture information, especially in desert areas. In mountain areas, the change of elevation makes the texture information complicated, which increases the number of matching points. In resident areas, the existence of numerous artificial features with rich texture information also increases the number of matching points. Presented in Table 3 and Figure 8, the number and distribution of matching points can impact the matching criteria of SD, RMSE, and CMR. With good feature points extracted in SAR-Moravec, the matching accuracy and CMR can be promoted. Furthermore, because of the shadow and layover phenomena in mountain areas [49], the matching accuracy (SD, RMSE) and CMR are both lower than resident areas. Due to the weak texture information, less feature points are extracted in desert and plain areas. It leads to less matching points in these slices. However, the matching accuracy of different slices all reach a subpixel level. Therefore, SAR-Moravec can deal well with different terrain types and is suitable to the matching procedure of complex-terrain SAR orthophotos in a large area.

Experiment of SAR Orthophoto Matching
In this part, the whole swathes of SAR orthophotos from adjacent orbits are used. Figure 9 illustrates parts of the matching results by SAR-Moravec. Table 4 presents the matching results of four data by methods including Moravec, SAR-FAST, SAR-SIFT, and SAR-Moravec. With evaluation criteria of SD, RMSE, CMR, GDD, and time consuming, the different matching results can be analyzed in detail. GDD in feature-based matching methods is obtained by the same regular grids utilized in area-based methods.
(a)  Shown in Figure 9, the SAR-Moravec has a good performance on distribution of matching points in these data. Compared with the other three methods, SD and RMSE criteria both show that the global matching accuracy of SAR-Moravec is the best. SAR-FAST and SAR-SIFT show moderate matching accuracy. The SAR-SIFT, which does not contain a grid control process, shows the lowest matching accuracy. In terms of GDD, the area-based methods with a grid control process perform  Shown in Figure 9, the SAR-Moravec has a good performance on distribution of matching points in these data. Compared with the other three methods, SD and RMSE criteria both show that the global matching accuracy of SAR-Moravec is the best. SAR-FAST and SAR-SIFT show moderate matching accuracy. The SAR-SIFT, which does not contain a grid control process, shows the lowest matching accuracy. In terms of GDD, the area-based methods with a grid control process perform well in matching points distribution. Especially in SAR-Moravec matching, the quality of feature points is better than other methods and the distribution performance is clearly promoted. In terms of time consuming, though slightly slower than Moravec, area-based SAR-Moravec is still more efficient than the other two modified methods. Considering all the evaluation criteria in the experiment, it can be concluded that SAR-Moravec shows good matching accuracy and matching efficiency in SAR orthophoto matching from adjacent orbits in a large area.

Conclusions
This paper proposed a modified area-based Moravec method for SAR orthophoto matching from adjacent orbits in a large area. In the feature points extraction, detecting more directions are added into the template to relieve the impact of speckle noise. Meanwhile, the main direction of the template is set up to adapt the rotation angle after terrain correction in SAR orthophoto. This template can improve the quality of feature points in SAR orthophoto matching. In the matching procedure, a coarse-to-fine matching process with Gaussian pyramid and grid control is designed to improve its adaptability in a large area and the homogeneity of feature points distribution. The coarse-to-fine matching in the Gaussian pyramid can effectively decrease the time consumption in a traversal process. The even distribution of matching points can ensure the stability of the transformation model and further improve the adaptability to different terrain types.
In the experiments, Radarsat-2 and Sentienl-1 data are used to analyze the effectiveness of the proposed method in SAR orthophoto in a large area. First, according to the experiment results of the comparison with Moravec, it can be demonstrated that the modified template improves the performance on feature points extraction in SAR orthophoto. Then, according to the experiment results of the adaptability to different terrain types, the proposed method is better than SAR-SIFT and, thus, more suitable for SAR image matching in complex terrain areas. At the end, according to the experiment results of the SAR orthophoto matching, it can be concluded that, though the time consuming is slightly longer than the original Moravec, it still has better global matching accuracy and better feature points distribution than other methods (SAR-SIFT and SAR-FAST). In conclusion, it is invulnerable to terrain types and capable of generating evenly distributed and high accuracy matching points. The proposed method can be widely applied in orthophoto matching from adjacent orbits in a large area.