TECHNOLOGY ON HIGH-ACCURACY DEM EXTRACTION FROM AIRBORNE INTERFEROMETRIC SAR

Remote sensing image mapping (aviation, aerospace, UAV) is the main form of current mapping. Optical image mapping has been a certain research foundation and played an important role in the mapping task. However, cloud, ice, snow, rain, and haze conditions seriously affect the interpretation and processing of optical images, making it difficult to meet the needs of rapid emergency response. Synthetic Aperture Radar (SAR) is currently the only method that can achieve all-day, all-weather rapid imaging under various weather conditions such as day, night, cloud, snow, rain, etc. With unique advantages unmatched by traditional optical remote sensing technology, it is an important technical guarantee for improving emergency response capabilities and promoting remote sensing and geographic information industries, and an earth observation technology that has been invested and vigorously developed internationally. Digital Elevation Model (DEM) plays an increasingly important role in the application of terrain information acquisition such as rapid acquisition of sudden disasters and surface deformation information acquisition. Rapid acquisition of large area high-accuracy DEM is one of the important research contents of SAR mapping. Based on interferometric SAR data combined with sensor orbit data, external DEM and other auxiliary data, DEM is finally generated through complex image registration, interference image generation, removal of flat ground effects, interference phase filtering, phase unwrapping, interference calibration and geocoding. The registration accuracy of multiple images greatly affects the accuracy of DEM. This paper proposes an improved registration method based on the conventional rough registration, pixel-level registration and sub-pixel registration. When the motion compensation of the sensor platform is unknown during imaging, it can meet the requirements of high-accuracy and large-scale rapid registration of dual-antenna polarization interferometric SAR. The experimental data is X-band airborne SAR data in Weinan region of Shaanxi Province, and DEM products with the elevation accuracy better than 0.5m is generated, which verifies the effectiveness, accuracy and fast processing ability of the method.


INTRODUCTION
Synthetic Aperture radar (SAR) technology has the advantages of all-day, all-weather, high resolution, wide coverage and high flexibility, which has unique advantages and irreplaceable value in specific applications. With the emergence of high-resolution, multi-polarization, and multi-band SAR systems and the rapid development of spaceborne, airborne, and UAV platforms, especially UAV SAR systems, maneuverability, flexibility and low cost, have developed rapidly in the past two years. It has a wide range of applications in terrain mapping, natural resource survey and monitoring, disaster emergency and other fields, and has broad application prospects, which will also put forward higher requirements for SAR data processing technology. The main methods of DEM data acquisition include ground measurement, digitization of existing topographic maps, GPS measurement, photogrammetry, stereo remote sensing, and interferometric radar. Traditional DEM acquisition methods such as ground measurement and digitization of existing topographic  Corresponding author. E-mail address: kq0608@163.com (Kang Qi) maps cannot take into account the two aspects of accuracy and update speed. GPS measurement method is only suitable for small areas. Although stereo remote sensing and optical photogrammetry methods can quickly obtain data collection on a national and even global scale, the two methods are seriously affected by the weather and have a serious impact on topographic mapping in tropical regions and perennial rainy regions. InSAR uses phase information to obtain ground elevation information. By coherently combining signals from two antennas, according to the phase and intensity information recorded in the data, combined with the parameters of the radar sensor itself, platform attitude, and orbit information, the geometric relationship is used to calculate the position change or three-dimensional position information of the object point. Compared with traditional measurement methods, the use of InSAR technology to obtain DEM has irreplaceable advantages, making it the most studied DEM acquisition method at present. In 2000, the United States carried out the significant SRTM topographic mapping mission. It took 13 days to obtain DEM data of all land areas from 60°N to 58°S. The products obtained by SRTM can be obtained free of charge, which has great and far-reaching significance for scientific research and production practice in various countries around the world (Graham L C . ,2005). In 2007, the German Space Center launched the first high-resolution radar satellite TerraSAR-X. The TerraSAR-X satellite and the second TanDEM-X satellite launched in 2010 formed a formation mode for high accuracy ground elevation measurement (Jin Guowang, 2007). In 2012, the TanDEM-X system has completed the first global surface mapping of the formation flight (Huang Shiqi, 2015). Around 2010, China's State Bureau of Surveying and Mapping successfully applied InSAR technology to the 1:50,000 topographic blank mapping project in western China (Zhang, J. et al, 2010). The extraction of DEM by InSAR technology has made great progress under the joint research of scholars all over the world, and many high precision DEM results have been obtained, but some aspects are still worthy of further study. The data registration accuracy directly affects the subsequent data processing. As the flight speed, direction and attitude of each data acquisition are not exactly the same, translation, scaling, rotation and even irregular local deformation occur between image pairs. At present, the widely used registration methods are to select a certain number of points in the equal interval of the master image, oversampling the slave images and using a certain registration measure function as the evaluation standard to obtain the offset of the selected points, fitting the offset of the whole image with polynomial, so as to complete the registration of interference image pairs (Chureesampant, K. et al, 2011). There may be problems such as few effective matching points, many false matching, large residual value after polynomial fitting, and unable to meet the registration accuracy (Natsuaki, R. 2011). This paper proposes an improved registration method based on the conventional rough registration, pixel-level registration and subpixel registration. When the motion compensation of the sensor platform is unknown during imaging, it can meet the requirements of high-accuracy and large-scale rapid registration of dual-antenna polarization interferometric SAR. Through complex image registration, interference image generation, removal of flat ground effects, interference phase filtering, phase unwrapping, interference calibration and geocoding, the DEM is finally generated.

Basic principles for DEM extraction in InSAR
The interferometry system has three antenna working modes: cross-track interferometry, along-track interferometry and repeat-track interferometry. At present, airborne interferometry mainly adopts cross-orbit interferometry. The working mode of cross-orbit interferometry is that two fixed antennas perpendicular to the flight direction are carried on the aircraft, which can accurately obtain three-dimensional information on the ground surface (Zebker, H. ,1986). The geometric relationship of airborne InSAR altimetry is shown in Figure1. A1 and A2 are the primary and secondary antennas on the flight platform, B is the baseline vector, αis the baseline Angle, is the perspective of A1 observing the ground target point P. The interference phase difference calculated from the radar complex image pair is: In Eq. (1), is radar wavelength, is the working mode of the antenna, when working in the "standard mode", = 1, when working in the "ping-pong mode", = 2, r1、r2 correspond to the slant distance from the antenna A1、A2 to point P.
The elevation value h of the ground target point P is: In 1 2 , calculated by the law of cosines: ( − ) = 1 2 + 2 − 2 2 2 1 (3) Based on the above three basic relations, the target point elevation information h in interferometric SAR altimetry is: The DEM of the observation area was finally obtained by solving the elevation of each pixel in the observation scene

Techniques for rapid registration
SAR image registration is the first step in InSAR processing, and its registration accuracy greatly affects the quality of the generated interferogram. The core idea is to restore points correspondence between SAR complex images, and achieve accurate registration between interferometric SAR complex images through geometric constraints between images, similarity criteria and sub-pixel registration methods. The specific process can be divided into: coarse registration, pixel registration and sub-pixel registration. Figure 2 shows the whole process of interference registration.

Figure 2.
Diagram illustrates the interferometric registration process.

Coarse registration:
Suppose that the row and column value corresponding to point P on the main image is ( , ), the row-column value corresponding to this point on the slave image is ( , ).
( , ) = ( , ) + ( , ) Coarse registration is to calculate the offset ( , ) of image pairs by using orbit data, external DEM and other auxiliary data.

Pixel level registration:
Using the window-based automatic registration technology, the two images are registered in the frequency domain or the spatial domain. Suppose M and S represent SAR master and slave images to be registered, respectively. Firstly, a pixel control point is selected in the master image M, and a matching window with × pixel points is selected with the this point as the center. Then select a search window with × pixels in the slave image S. The search window must be large enough to ensure that the window corresponding to the matching window can be found within it. In the search window, the corresponding window with the same size of the matching window is selected by different pixel offsets in the order of row and column, and the correlation coefficients Eq. (6) between the matching window and the corresponding window are calculated respectively.
( , ), ( , ) respectively represent the two complex data at the corresponding position ( , ) of the two matching windows; * represents a conjugate complex number; represents the size of the window. The maximum value is selected from the obtained ( − + 1) 2 correlation coefficient values through calculation. The window corresponding to the correlation coefficient value is the window that matches the matching window. The relative offset between the position of the corresponding window center pixel and the control point position of the matching window center pixel is the required registration offset.

Sub-pixel level registration:
The steps and methods of sub-pixel level registration is basically same to pixel-level fine registration The sub-pixel level registration is divided into three steps: Step 1: Oversampling the master image and the slave image respectively. Before determining control points and searching for points with the same name, the interpolation method is used to oversample the image pairs, and the interpolation interval is selected according to specific requirements.
Step 3: Offset fitting and resampling of slave images. The registration method in this paper divides the image into several small blocks and takes the quadratic polynomial model as the local deformation model. The premise of the correct solution of the registration model is that the searched points with the same name are correct. The points with the same name are determined by counting the multiple correlation coefficients of the sliding window, and the point with the maximum coherence is the same name. However, affected by many factors, such as forest, grassland, water surface, relatively smooth road and so on, all belong to planar distribution targets, the point with the largest coherence coefficient may not be the corresponding point. Because the grid point statistical properties of these planar distributed targets are similar to those of their surroundings, there are multiple peaks in complex coherent images. The reliability of the master image and the corresponding slave image grid points can be judged by the following two conditions.
(1) Maximum correlation coefficient module | | and sublarge correlation coefficient module | _ | should have good separability. Define the Rate, The closer the Rate is to 1, the flatter the peak of the coherence coefficient, and the lower the reliability of the homonymous point. The closer the Rate is to 0, the closer the coherence coefficient graph is to the point response function, and the higher the reliability of the corresponding point is.
(2) Define the maximum and sublarge coherence coefficient distance: (x, y) is the image coordinate of maximum coherence coefficient; (j, k) is the image coordinate of sublarge coherence coefficient. If the coherence coefficient curve is similar to Sinc function, the closer Dis is to 0, the more reliable the corresponding point is.
For the surface distribution target, due to the similar statistical characteristics, the Dis is relatively large and the reliability of the corresponding point is low. If the grid point after precise registration satisfies the above two conditions, the grid point is regarded as the correct corresponding point and is used to solve the second registration model coefficient, otherwise, the registration accuracy of this point is poor, and the addition of the model coefficient will increase registration errors, grid points need to be processed by other methods. If the grid points with low reliability are directly eliminated from the block, the registration accuracy of the whole block can be effectively improved and the propagation of registration errors can be prevented. However, in this way, there is no grid point assistance in the vicinity of the culled point, and the geometric deformation can only be controlled by other grid points that are far away, which not only improves the overall registration accuracy, but also reduces the local registration accuracy. In this paper, a Delaunay triangulation network is constructed for the low reliability points in the block to control the error propagation as much as possible and improve the local registration accuracy. If the area of a triangle in the constructed triangular network is larger than the specified threshold, the triangle is discarded, and the registration relationship during resampling is based on the quadratic registration polynomial of the block. If the area of this triangle is smaller than the threshold, the registration relationship during resampling is based on the affine transformation formed by the triangle, which is called an effective triangle. The above mentioned method is improved on the basis of conventional three-step registration method. The image is divided into blocks so that different positions of the whole image are processed with different deformation models. The reliability of regular grid points is detected. Each high-reliability point is registered by a quadratic polynomial model, and the lowreliability point is registered by constructing a triangulation affine transformation method, which improves the adaptability of registration. The triangles in the triangulation that exceed the threshold area are culled, which limits the error propagation and improves the overall coherence well. The above improvements can meet the requirements of high-precision and large-scale fast registration of dual-antenna interferometric SAR when the motion compensation of the sensor platform is unknown during imaging.

Process for high-accuracy DEM extraction
The process of extraction DEM in InSAR mainly includes complex image registration, interference image generation, removal of flat ground effects, interference phase filtering, phase unwrapping, interference calibration and geocoding. The processing flow is shown in Figure 3.

Removal of flat ground effects:
The phase of interferogram obtained by conjugate multiplication of complex image includes flat ground phase, topographic phase, displacement phase, atmospheric phase and noise phase. Flatland effect refers to the phenomenon that the interference phase changes periodically in the direction of distance and azimuth caused by the uneven flatland on the reference earth surface. Flat ground phase makes the interferogram appear dense fringe, which to a certain extent covers the change of interference fringe caused by terrain change. The removal of flat phase is convenient for subsequent interference filtering and phase unwrapping. The flattening method based on orbital parameters is used to eliminate the ground effect. The specific algorithm is as follows: (1) Calculate the coordinates of the image center point corresponding to the point on the ellipsoid. According to Doppler equation, distance equation and ellipsoid equation, the geodetic coordinates (X,Y,Z) of the master image pixel ( 1 , 1 ) corresponding to the point P on the reference ellipsoid are calculated.
(2) Calculate the effect of flat ground. According to the coordinates of point P and satellite, the oblique distance is calculated, then the flatland effect can be expressed as the phase difference between the point corresponding to the center point on the ellipsoid and two satellites, the phase represented by the parallel baselines: (3) Polynomial fitting. According to the resolution of oblique distance, the phase values of a series of reference points covering the image are derived from the center point of the imaging region, and the two dimensional polynomial is formed by the coordinate of the image points, and then the polynomial coefficients are solved by the least square method.
(4) Removal of the phase of the flat earth.
Using the polynomial coefficients obtained by (3) and the coordinates of each image point, the flat ground effect of each point was calculated and removed.

Interference phase filtering:
Due to the system noise generated in the SAR imaging process, the interference phase has very serious interference. The interference phase filtering adopts the goldstein filtering method, which not only filters out the noise interference, but also maintains the continuity of the interference fringes.

Phase unwrapping:
The phase value in the interferogram is the main value of the interference phase, and its value range is [− , ), and the interference phase value will have a step jump at the edge of the interference fringe. In order to calculate the elevation value of the ground point from the interference phase, the relative integer number of cycles between each interference phase in the whole interferogram must be determined. That is the phase unwrapping process. Phase unwinding mainly adopts the minimum cost flow algorithm, which is a network planning algorithm for optimizing objective function and solving optimal objective, and has good unwinding effect.

Interference calibration:
The purpose of interferometric SAR calibration is to correct the phase offset and interference parameter deviation of system equipment so as to improve the accuracy of DEM generation. The correction of the phase offset of system equipment is called the equipment calibration of the interference system, and the correction of the interference parameter deviation is called the interference parameter calibration. There are two methods of interference calibration: interference calibration based on sensitivity equation and interference calibration based on phase correction. Interference calibration based on sensitivity equation is to establish the relationship between the interference parameter deviation and the interference elevation measurement deviation through sensitivity equation, and then estimate the interference parameter deviation by solving equations, and then use the parameter deviation estimation value to correct the interference parameters. The basic idea of interferometric calibration based on phase correction is to establish the relationship between the interference parameter deviation and the the interference phase and the reference phase's difference through sensitivity equation. The estimation value of the interference parameter deviation is obtained by solving the equations, and the interference parameter deviation is corrected by using the estimated value, and the interference height is reconstructed by using the corrected interference parameter, so as to improve the accuracy of the interference height.

Geocoding:
Elevation calculated with InSAR is only a numerical set of ground elevation corresponding to each point of the reference image, the height value is based on the coordinates of the SAR image. SAR image coordinates must be converted into a more general reference coordinate system (such as geographical coordinate system or UTM plane coordinate system) before practical application of DEM and SAR image can be carried out. Geocoding mainly goes through location and resampling. Positioning is divided into direct positioning and indirect positioning, both of which adopt R-D model.

EXPERIMENTS
Type text single-spaced, with one blank line between paragraphs and following headings. Start paragraphs flush with left margin.

Data
In this paper, we used X-band airborne interferometric SAR data taken in Weinan City, Shaanxi Province, China. The data resolution is 0.2m, the image size is 24,576 pixels in azimuth, 20,000 pixels in range, 0.101491m in azimuth sampling, 0.136236m in range sampling, the azimuth direction of the image bandwidth is 2.5km, with range direction is 3.5km, and the central viewing angle of the image is 53°. The detailed parameters are as follows in Table 1 Parameter The data has a total of 10 strips, which are divided into two sideview directions, of which 5 strips fly from east to west and look to the north side, and the other 5 strips fly from west to east and look to the south. The data coverage is 18km from east to west, 10km from north to south, 20% within the overlapping strips and 30% between the strips.

Experimental processing results
The master and slave SLC data are registered by the improved registration method using the parameter file. After accurate registration, the phase difference can be obtained at each pixel by conjugate multiplication of the corresponding pixel values to form an interference fringe pattern. Figure 4 below shows the original SAR image, the interferogram, the interferogram obtained after filtering and flattening, and the unwrapped phase map, and finally the generated DEM. According to the relevant parameters provided by the parameter file, the method based on the precise track is used to remove the phase difference caused by the flat ground effect, and the Fourier transform is used for the interferogram to convert the spatial domain image into the frequency domain image, and then filter the image by the filter function, and then through Fourier inverse transformation to the spatial domain, the filtered flattened interferogram is obtained. It can be seen in Figure 4. (c) that in the interferogram after removing flat earth effect and filtering, the speckle noise is significantly reduced. For the filtered interferogram, set the coherence threshold and other parameters, and use the interferometric quality map as the weight. Because the experimental area is located in a hilly area, the minimum cost network flow phase unwrapping method is used to list the unwrapping equations and solve the real phase difference. Then we obtain the unwrapped phase map. Figure 4 shows that the unwrapping effect is good, and the true value of the phase is restored. According to the relationship between phase value and ground elevation, the mathematical model is used to calculate the elevation value on each pixel, and according to the relationship between the slope range measured by radar and ground range, the slope file is converted into ground range file, and the elevation map and ground range intensity image are obtained. Based on the earth ellipsoid model and range Doppler model, the image coordinates of elevation map are located, and the corresponding relationship between the image coordinates of elevation map and DEM is obtained. Finally, the final DEM is obtained by resampling. It can be seen that there is a reversal in the direction of the original image. This is because the original data is the range direction from left to right, and the azimuth direction is from top to bottom, which does not conform to the coordinate expression of the topographic map. However, the longitude of the geocoded DEM increases from left to right and latitude from bottom to top, which is consistent with the actual topographic map expression. DEM results of 5 strips from east to west and 5 strips from west to east in the survey area were spliced and rendered respectively to produce DEM products. DEM products with a standard size of 10,000 map were obtained by cutting according to the size of map. Figure 5. and Figure 6. show two DEM products.  In order to evaluate the accuracy of DEM obtained by experiment, 70 elevation checkpoints measured in the field in this region were selected for comparison, and the absolute accuracy of two 1:10,000 SARDEM results was tested and evaluated, as shown in the following As can be seen from Table 2., the max root mean square error (RMSE) of DX_I49G039026-DEM is -2.183 meters and RMSE is 0.931 meters; the max RMSE of XD_I49G039026-DEM is -2.298 meters and RMSE is 1.024 meters. The DEM products extracted in this paper meet the industry standard "Basic Geographic Information Digital Achievements 1:5000, 1:10000, 1:25000, 1:50000, 1:100000 Digital 1:10000 digital elevation model in "Elevation Model" requires a first-class accuracy of 1.2m for hilly land.

CONCLUSIONS
Using InSAR to extract DEM, due to its high degree of automation and high efficiency, in the case of good coherence of most radar images, the acquired DEM has high precision and can meet the needs of users. At present, it is more and more widely used. Aiming at the problem of DEM extraction with high precision, this paper uses interferometric SAR data combined with sensor orbit data, external DEM and other auxiliary data to obtain DEM products through image registration, removal of flatland effect, interference filtering, phase unwrapping, geocoding and other processes. The high-precision registration technology meets the needs of large-scale rapid registration, and the high-precision DEM extraction achieves optimization and breakthrough in the efficiency and accuracy of key processing processes.
There are many factors that can affect the InSAR processing results, such as the coherence of data, the influence of the noise and the effect of phase unwrapping. Different regions may also have different interference results due to terrain fluctuations and differences in surface coverage. Further research on the above aspects can be carried out in the future to improve the DEM extraction accuracy.