RESEARCH ON A DEM COREGISTRATION METHOD BASED ON THE SAR IMAGING GEOMETRY

Due to the systematic error, especially the horizontal deviation that exists in the multi-source, multi-temporal DEMs (Digital Elevation Models), a method for high precision coregistration is needed. This paper presents a new fast DEM coregistration method based on a given SAR (Synthetic Aperture Radar) imaging geometry to overcome the divergence and time-consuming problem of the conventional DEM coregistration method. First, intensity images are simulated for two DEMs under the given SAR imaging geometry. 2D (Two-dimensional) offsets are estimated in the frequency domain using the intensity cross-correlation operation in the FFT (Fast Fourier Transform) tool, which can greatly accelerate the calculation process. Next, the transformation function between two DEMs is achieved via the robust least-square fitting of 2D polynomial operation. Accordingly, two DEMs can be precisely coregistered. Last, two DEMs, i.e., one high-resolution LiDAR (Light Detection and Ranging) DEM and one low-resolution SRTM (Shutter Radar Topography Mission) DEM, covering the Yangjiao landslide region of Chongqing are taken as an example to test the new method. The results indicate that, in most cases, this new method can achieve not only a result as much as 80 times faster than the minimum elevation difference (Least Z-difference, LZD) DEM registration method, but also more accurate and more reliable results.


INTRODUCTION
Due to the globe coverage from various large-scale missions, such as Shutter Radar Topography Mission (SRTM)and the Advanced Spaceborne Thermal Emission and Reflection (ASTER) Global Digital Elevation Model (GDEM), and feasible access to technology, such as airborne light detection and ranging (LiDAR) and unmanned aerial vehicle (UAV) technology, recently, the DEM has been widely applied in many aspects of earth surface research, such as cartography, climate modelling, geophysics, geomorphology, geology and soil science.For deformation mapping using the Interferometric Synthetic Aperture Radar (InSAR) technique, DEM must provide the topography phase and geocode the SAR coordinate results into the geodetic coordinate system (Yang et al. 2011).However, before practical application, multi-source, multitemporal DEMs must be assessed due to their different height precisions and systematic horizontal deviations (Rodriguez et al. 2006).Furthermore, the deviation of the horizontal position between two DEMs will seriously affect the assessment of the DEM height precision and cause a height shift when multi-source, multi-temporal DEMs are merged.However, less attention has been paid on systematic horizontal deviation effect for the DEM precision assessment, and no measure has been taken to eliminate or reduce the systematic deviation effect.Therefore, in the applications of multi-source DEMs, the horizontal position of the DEM must be corrected systematically, namely, the DEMs must be well coregistered.
Currently, commonly applied DEM coregistration methods can be divided into two categories: one category is the coregistration-based ground control points (GCPs) approach, and the other category is the matching without control points approach.The former approach is actually the geodetic coordinate conversion problem.The prerequisite is that there existing a certain number of control points with a good distribution.The latter approach can also be divided into two classes: feature-information based methods and geometricrelation based methods (Akca et al. 2006).Due to the lack of requirement of feature points and high matching precision and efficiency, geometric-relation based methods are the primary DEM coregistration methods, which include the minimum height difference algorithm (Least Z-difference, LZD) (Rosenholm et al. 1988;Pilgrim 1996), Interactive Closest Points method (ICP) (Besl 1992), Least Squares 3D Surface algorithm (LS3D) and some improved algorithms (Karras et al. 1993;Pilgrim 1996;Zhllin et al. 2001).The ICP and LS3D algorithms are proposed for use on 3D terrain data with irregular distribution, such as 3D laser point cloud data; the LZD algorithm (which is the most commonly used DEM coregistration) is proposed for use on regular grid data, such as SRTM DEM data, because it provides an overall optimal performance (Yang 2012).However, these methods may lead to convergence errors, or even divergence, for the sake of local optimization, which in theory mainly depends on the choice of the initial value (Gruen et al. 2004).In practice, if the local difference between two DEMs is large, or the large DEM points are involved, then the computation is quite time-consuming when seeking the corresponding grid height value of the slave DEM through repeated two-dimensional interpolation during the iterative computation process.Therefore, this paper proposes a fast regular grid DEM coregistration method with the aid of the simulated DEM intensity image.Two different resolution DEMs are analysed as an example and compared with the performance of the LZD algorithm.

THE PRINCIPLE AND METHOD OF COREGISTRATION
The principle of the newly proposed DEM coregistration method is as follows.First, without obvious terrain features on DEM images, intensity images are simulated based on a given SAR imaging geometry and the ground backscattering characteristics, which can strengthen the DEM terrain feature information.Second, two-dimensional cross-correlation is calculated based on the simulated intensity images, which can be divided into two steps: coarse registration and fine registration.Coarse registration is used to obtain the initial DEM offset with FFT operation of the whole image; the fine registration is used to obtain 2D offsets for several distributed windows, within which the FFT operation is performed to obtain the 2D offsets by searching for the maximum SNR.Next, the 2D polynomial robust least squares fitting method is applied to the high-quality 2D offsets, which enables 2D high precision coordinate mapping relationship between two DEMs to be obtained through iterative calculation.Finally, the slave DEM is interpolated and resampled to the reference DEM.The flowchart of DEM coregistration method based on the SAR imaging geometry is shown in Figure 1.
Figure 1.Flowchart of the DEM coregistration method based on the SAR imaging geometry

DEM intensity image simulation
The terrain has a strong influence on the geometric and radiometric characteristics of a SAR image (Wegmuller 1999;Loew et al. 2007;Small et al. 2009;Frey et al. 2013).Given a SAR imaging geometry, SAR backscatter images (i.e., DEM intensity images) from two DEMs in geographic coordinates can be simulated.If we compare figure 5 with figure 6, there is no significant features can be seen from both high and low resolution DEMs, while obvious features, such as ridge, valley line, and slope changing points can be well recognized from the DEM simulated intensity images.Subsequently, two DEMs are coregistered based on their intensity images.
Because the ground objects scattering properties from DEM are unknown, the simulated DEM intensity image mainly depends on two factors: the actual ground area corresponding to the unit pixel and the scattering caused by a local incidence angle (Loew et al. 2007).Wegmuller (1999) noted that SAR image simulation can be regarded as the product of the actual surface area represented by the unit pixel and the local incidence angle related experience function, and Zhang et al. (2002) reported that the SAR radiation distortion resulted from area effect and a local incidence angle effect.To account for these two factors, the SAR intensity simulation includes two steps: the first step is to obtain the relative size of the pixel area, i.e., the establishment of the pixel normalization factor; the second step is to construct the local incidence angle related experience function.
The satellite position vector and velocity vector at a given epoch can be abstracted from the SAR parameter file.Next, according to the rigid SAR Range-Doppler (RD) positioning model and the available DEM, the SAR ( ) , ( SAR j i ) satellite position vector and velocity vector, which correspond to a ground point ) , ( P j i , can be accurately calculated.As shown in Figure 2, the local incidence angle  is determined by the surface normal vector n  at point ) , ( P j i and the look vector .In Figure 2, i represents the north-south direction, and j represents the eastwest direction.As the SAR satellite flight direction is approximately north-south, the number of adjacent points at the direction of east-west's corresponding relative satellite position moving is much less than the number of adjacent points of north-south direction.Therefore, we assume that the satellite positions to the adjacent point of east-west direction are negligible and that the surface normal vector is calculated as follows: ) Figure 2. Schematic plan of Figure 3. SAR imaging normal surface vector.geometry (Wegmuller 1999).
Figure 3 shows the SAR imaging geometry, where an arbitrary point on the ground is taken as the origin of coordinates, the satellite flight direction is the X axis, the spherical surface normal vector is the Z axis, r  is the radar line of sight vector, n  is the surface normal vector, m  is the horizontal normal vector of the image, R represents the distance between ground point and the earth centre, S is the distance from the satellite to the centre of the earth ,  is the local incident angle, and  is the angle between the ground surface normal vector and the horizontal normal vector of the image.
The projection angle  is determined by the image surface normal vector and the ground surface normal vector, as shown in the following formula.
Where, m  is the image plane normal vector which is determined by the line of sight vector r  and the satellite velocity vector (i.e., the azimuth direction).
For a generic ground grid, the intensity value of the grid can be regarded as the corresponding pixel intensity value, which is mainly determined by the pixel size.The flat surface area corresponding to one pixel can be expressed as where, r is the SAR satellite slant range pixel size, and az is the azimuth pixel size, and  is the incidence angle.
The actual ground area corresponding to the unit pixel is (Wegmuller 1999) The pixel size normalized factor is defined as

The intensity image registration
The intensity image registration is mainly based on the intensity FFT cross-correlation method.Compared with other image registration methods, FFT is easy to implement, sensitive to noise, and has high calculation efficiency (Brown 1992).
The basic principle of the FFT method applied to image matching is described in several references (Brown 1992;Reddy et al. 1996).When a signal shifts the time axis or space axis, the FFT amplitude spectrum remains unchanged, whereas an additional phase shift is produced.Similarly, when a 2D signal moves in two directions, the FFT amplitude spectrum also remains unchanged, and the phase spectrum shifts along the two directions.Assuming that the relative shift of image In most cases, the shifts for different image sections are not constant between two images, so distributed small windows will be set to calculate different shifts.We will discuss it in next section.
The DEM simulated intensity image registration includes the following five steps.
(i) One large window centred at the area with a clear topographic feature is chosen during the coarse registration step.By implementing the FFT operation, the initial offset is obtained between the slave DEM (DEM2) intensity image and the reference DEM (DEM1) intensity image, which will be set as the initial value during the fine coregistration step.

(ii)
During the fine coregistration step, m n  uniformly distributed grid points in the DEM1 intensity image are set as the reference points, and centred at each point, a window with N 2 pixels ( N is between 6 to 9) is determined.
On the basis of the initial offsets, a template of identical size to the above-mentioned window is extracted in DEM2.These two windows are subjected to the following operations: two or more times oversampled, FFT, conjugate multiplication, and inverse FFT; subsequently, 2D offsets and their SNR value can be readily obtained based on the peak value of the real part.
(iii) The offsets ) , ( with SNR values larger than the given threshold, such as 6-7 dB, will be preliminary obtained.Next, the robust least squares method is applied to fit the second-order polynomial by using M 2D offset measurements in the range and azimuth directions, as given by equations ( 11) and ( 12), respectively.
The error equations in matrix form become  13) can be obtained as follows,

The corresponding offset residual equation and standard deviation equation are
After setting the offset residual threshold (e.g., 2.5 ˆ), the iterative robust least squares solution is performed.
(iv) By using equation ( 14), the offset corresponding to each grid of the DEM is calculated as follows.
Where A and B are the polynomial coefficients, N  and E  represent offsets corresponding to   E N, at DEM2. ( 5) After the coordinate transforming, re-sampling and gridding on DEM2, the two DEMs are well coregistered.

LZD coregistration method
To evaluate the newly proposed DEM coregistration method, we compare the results of the new method to the LZD method.The basic idea of the LZD method is briefly introduced as follows.Taking the points from two DEMs with the same horizontal coordinates to be the corresponding points, the objective function that minimizes the sum of height difference square between corresponding points is established.Finally, 7 transformation parameters, including 3 rotation parameters ) , , ( and 1 scale parameter S , are calculated using least-squares criteria.This procedure is iteratively processed until the changes of the parameters satisfy the given thresholds.The details of the LZD method can be found in the related literature (Rosenholm et al. 1988;Karras et al. 1993;Pilgrim 1996;Pilgrim 1996).

TEST AREA AND DATA
The Yangjiao landslide area of Chongqing municipality, China is selected as the test area; the area is located downstream of the Wujiang River of south-western China.The area has dense vegetation and steep terrain, where maximum elevation is approximately 1370 m and the minimum height is approximately 160 m.This region is prone to rock collapse and landslide because of the cliff terrain and the water during the rainy season (Figure 4) .The airborne LiDAR DEM is set as the reference DEM and the SRTM DEM is tested for coregistration using the proposed method and the LZD method.    .At each iteration times, two-dimensional cubical polynomial interpolation is implemented to obtain the height value according to the solutions of the last iteration at the given grid point.In this test, the thresholds of are set to be less than 0.001 of a grid point, is less than 1″, and is less than 0.01.In total, 150 iterations are processed.The runtime is approximately 47900 seconds using the Matlab program under the same computational environment as described above.Theoretical speaking, this method may result in local optimization solution, once the data volume is too big and the height changes are complex (Akca and Gruen 2005 ). is less than 0.01.In total, 150 iterations are processed.The runtime is approximately 47900 seconds using the Matlab program under the same computational environment as described above.Theoretical speaking, this method may result in local optimization solution, once the data volume is too big and the height changes are complex (Akca and Gruen 2005 ).

Analysis of the test results
Figure 7 shows the height difference image between the SRTM DEM and the LiDAR DEM before and after the coregistration process.After coregistration, both approaches can well eliminate the systematic errors, that is, the two DEMs are well coregistered.Because the difference of two standard deviations is only 0.54 m, the precisions of the two methods are quantitatively similar.However, the mean value of our method is approximately 1.4 m, whereas the mean value of the LZD method is approximately -3.3 m, which indicates that our results are more reliable.after coregistration using our method, and (c) after coregistration using LZD method, and (d) is the comparison of consuming time between Our method and LZD method.
For close comparison of the coregistration differences between the two methods, two profiles along AA' and BB' in Figure 5(a) are analysed.Figure 8 shows the height differences along line AA' over a small terrain change area.Figure 9 shows the height differences along line BB' over a steep topography area, where the slope angle reaches almost 90 degrees.From Figure 8(a) and Figure 9(a), significant systematic deviations are observed between the two DEMs before coregistration, whereas after registration, both methods are found to eliminate the systematic offsets significantly.Furthermore, if we compare Figures 8(c) and 8(d), the mean height difference of our method is determined to be approximately 0.5 m, which is 10 times better than the LZD method, and the variance of our method is also smaller than that of LZD method.Similar conclusions can also be drawn from the steep terrain profile shown in Figure 9(b).

CONCLUSIONS
The SAR imaging geometry based DEM coregistration method was proposed in this paper, which takes advantages of twodimensional image patching.First, the procedure to simulate DEM intensity image based on the given SAR imaging geometry was discussed.The DEM intensity image more clearly shows the surface information compared with the DEM height map.Next, the FFT cross-correlation method is applied to calculate the 2D offsets, which accelerate the DEM coregistration process.Comparing with the traditional LZD method, the proposed method is found to operate as much as 80 times faster on real DEM data of 3435× 2348 pixels.Moreover, the new method can achieve more precise and more reliable results, especially for a complicated terrain area.In addition, we can imagine our method operating more efficiently when the DEM data is large enough.
the surface normal vector n  is determined by the look vectors and the satellite position vectors, which correspond to adjacent this stage, we have obtained the first part of the simulated intensity image.The second part of the method is determining the empirical function associated with the local incidence angle, that is, the change of the scattering mechanisms caused by the local incidence angle change, which can be represented by the cosine function model, the cotangent function model, or the polynomial function model.We use the cotangent function model to highlight the terrain features.According to equation (6), the intensity value of a DEM grid point is b As a result, we can obtain the following: parameters in equation ( of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLII-3, 2018 ISPRS TC III Mid-term Symposium "Developments, Technologies and Applications in Remote Sensing", 7-10 May, Beijing, China This contribution has been peer-reviewed.https://doi.org/10.5194/isprs-archives-XLII-3-1333-2018| © Authors 2018.CC BY 4.0 License.

Figure 4 .
Figure 4. Geographic image of the Yangjiao landslide region.

Figure 5
Figure 5 (a) shows a 2-meter resolution airborne LiDAR DEM, which is generated from the LiDAR point cloud data by interpolation.The height accuracy of airborne LiDAR DEM is 0.7 meter (China Institute of Geo-Environment Monitoring 2014).The size of the test region is 3435 × 2348 pixels.The horizontal and height data of the test region are both World Geodetic System 1984 (WGS-84).Figure 5(b) shows the SRTM DEM with 90-m resolution.The height accuracy of SRTM is approximate 16 meter globally (CGIAR-CSI 2009).The horizontal datum is WGS-84 and the height datum is earth gravitational model 1996 (EGM96) geoid.

Figure 5 .
Figure 5. Two DEM images in Yangjiao region, Chongqing: (a) LiDAR DEM and (b) SRTM DEM. 4. RESULTS AND ANALYSIS 4.1 Test procedure 4.1.1Our method: (i) The height data are unified by taking into consideration the EGM96 height anomaly, which is calculated through an online service provided by the non-profit organization UNAVCO.Next, the SRTM DEM elevation is converted to WGS-84 geodetic height.(ii) Basically, regarding two different resolution DEMs comparison, two ways can be usually adopted to unify the DEM resolution.The one is downsampling, while the other is oversampling.As the accuracy of image coregistration is highly correlated to pixel size, in this test, SRTM DEM is oversampled to the same resolution as that of the LiDAR DEM with the bicubic spline method to hold the high resolution of LiDAR DEM.(iii) The two DEMs are simulated as intensity images based on a given SAR imaging geometry, as shown in Figure 6.(iv) The SRTM DEM is coregistered to the LiDAR DEM.The process takes approximately 600 seconds using the Matlab program under the computational operating environment of a Dell E5420 server, Ubuntu 12.04 OS, CPU Intel(R) XEON 2.5 GHZ, 24 GB RAM.

Figure 6 .
Figure 6.Simulated intensity images of two DEMs: (a) LiDAR DEM and (b) SRTM DEM.4.1.2LZD method: Steps (i) and (ii) are the same as the ones discussed above, that is, the same nominal coordinates are used as the initial coordinates.(iii) The registration equation is established.Without loss of generality, we use the four parameters model for LZD coregistration because here we only consider the horizontal position derivation of the two DEMs.The initial values of four parameters are set to ) 1 , 0 , 0 , 0 ( .At each iteration times, two-dimensional cubical polynomial interpolation is implemented to obtain the height value method: Steps (i) and (ii) are the same as the ones discussed above, that is, the same nominal coordinates are used as the initial coordinates.(iii) The registration equation is established.Without loss of generality, we use the four parameters model for LZD coregistration because here we only consider the horizontal position derivation of the two DEMs.The initial values of four parameters are set to ) to obtain the height value according to the solutions of the last iteration at the given grid point.In this test, the thresholds of y x T T   , are set to be less than 0.001 of a grid point, h R  is less than 1″, and S

Figure 7 .
Figure 7. (a), (b) and (c) are height difference images between SRTM DEM and LiDAR DEM: (a) before coregistration, (b)after coregistration using our method, and (c) after coregistration using LZD method, and (d) is the comparison of consuming time between Our method and LZD method.