A SINGLE-PASS AIRBORNE INTERFEROMETRIC CALIBRATION METHOD RESEARCH FOR DEM MAPPING

This paper presents a novel single-pass airborne interferometric calibration method with a limited Number GCPs for InSAR digital elevation model (DEM) mapping. The proposed method is based on a rigorous three-dimensional model for a single-pass airborne InSAR system. The corrected InSAR parameters of baseline length, baseline inclination, near slant range, and Doppler centroid frequency, as well as phase offset, can be jointly solved via a unified optimization procedure in terms of the constructed threedimensional geometric model using a limited number ground control points (GCPs). The proposed method is evaluated on real data of the CASMSAR system in X -band, the final DEMs generated by the calibration processing achieve a high accuracy level (1-3 m standard deviation), even in the presence of only 3-5 GCPs. .  Corresponding author Email: lulj@casm.ac.cn


INTRODUCTION
The airborne InSAR system calibration is a necessary work for DEM mapping. This means that all the InSAR parameters which could cause phase errors have to be corrected before DEM mapping. The significant error sources come from the radar system itself and the navigation system (Wimmer et al. 2000). Systematic phase errors and statistical phase errors constitute radar-induced phase errors. Systematic phase errors are induced by geometric misregistration due to the existence of a squint angle, and statistical phase errors are a high frequencyerror that can be regarded as noise. Navigation system phase errors play a critical role because these large errors need to be compensated for. Aircraft flight position and orientation, which are generally measured by a global positioning system (POS) instrument, have a direct influence on several important InSAR parameters, such as baseline length, baseline inclination, near slant range, etc. In certain circumstances (e.g., strong wind or airflow), the precision of the POS instrument may be reduced and affect the estimation of the system parameters; consequently, the height of the final terrain model cannot be retrieved correctly. The second prerequisite for high-precision DEM mapping is the fulfilment of the retrieval of the absolute phase from the wrapped interferometric phase. This procedure is named "phase unwrapping" in the InSAR processing chain (Bamler and Hartl 1998). However, there will be a residual phase-offset value in the interferogram after phase unwrapping, which is a constant value for the whole scene, and which must be estimated (Mura et al. 2012). In summary, from point of view of the high-precision DEM generation, airborne InSAR calibration involves the two aspects of InSAR parameter correction and phase-offset estimation.
For the errors induced by uncalibrated or poorly calibrated airborne InSAR systems, one effective solution is using a number of corner reflectors (CRs). The CRs are precisely measured as ground control points (GCPs) before being employed in the construction of the InSAR geometric model. A number of alternative optimization methods with respect to InSAR geometric models have been proposed to estimate the uncertainties in estimations of InSAR parameters (Sun et al. 2000), (Mallorqui, Rosado, and Bara 2001), (Sun et al. 2011). Furthermore, the residual phase offset is often estimated by use of GCPs within the scene, or through the use of an already calibrated area. The required phase offset can be retrieved by solving the optimization problem whose efficient implementation is provided by the least squares (LS) solution using GCPs (Sun et al. 2000). In order to save the cost and time involved with the deployment of CRs, phase-offset estimation methods that do not require the use of GCPs have been developed (Mura et al. 2012), (Ferraioli, Ferraiuolo, and Pascazio 2008), (Gatti et al. 2011) and (Perna et al. 2015). The condition of a well-calibrated interferometric airborne SAR system has to be fulfilled for the above methods to work well, as mentioned previously. This paper proposes a new single-pass airborne interferometric calibration (SAIC) method, which is suitable for single-pass airborne data pairs acquired by uncalibrated or poorly calibrated systems. The SAIC method involve integrating the InSAR parameter correction and phase-offset estimation into a unified optimization procedure, in terms of the reconstructed three-dimensional geometric model using a limited number of GCPs. In the SAIC method, a rigorous three-dimensional model for the single-pass airborne InSAR system, instead of the traditional two-dimensional InSAR geometry model constructed, and the least squares (LS) principle is introduced to solve the built equations using only a few GCPs. The performance of the proposed SAIC method is evaluated using data from the Chinese CASMSAR system in the X-band from a flat area of Pucheng in China.

THREE-DIMENSIONAL RECONSTRUCTION MODEL
A moving coordinate system is first introduced for the threedimensional model construction, in which the InSAR parameters and phase offset can be clearly represented. In the following, a new three-dimensional model is constructed for the airborne InSAR system.

Moving Coordinate System
A moving coordinate system is a local coordinate system that is referenced to a Cartesian coordinate system. The master antenna phase center is positioned in the original point of the moving coordinate system, for which the three orthogonal basis vectors are defined as (1) where is the velocity vector at one moment, and is the baseline vector. The unit vector can be calculated by the ratio of vector and scalar . The unit vector is the cross product of two vectors of and , and the unit vectors is the cross product of three vectors of , and .
A conjunction of the three unit vectors of , , and constitutes the moving coordinate system, which is called the VNW system, as shown in Figure. 1. The transformation matrix from the Cartesian coordinate system to the VNW system is written as . The opposite is also addressed: (2) Given that is the orthogonal matrix, . In the Cartesian coordinate system, as shown in Figure. 1, the three-dimensional coordinates of the target is described as the sum of the master antenna phase center and the look vector : (3) The look vector can be expressed as: (4) where is the slant range corresponding to the master antenna and is the unit look vector. The parameters to be solved are implicated in the unit look vector , which is generally expressed in the moving coordinate system.

Three-Dimensional Model
The rigorous InSAR geometric model using the slant range, the Doppler and the interferometric equations can achieve a precise position for a pixel from the image coordinates to the object space coordinates (Crosetto 2002). However, the InSAR parameters and phase offset are implicated in the three nonlinear equations; in this case, an analytical solution cannot be acquired. The idea of orthogonal decomposition for look vector in the VNW system is an effective solution.
In the VNW system, the look vector can be decomposed into three orthogonal components: where is the radar wavelength, and is the Doppler centroid frequency. The baseline vector is decomposed into parallel component parallel to the velocity direction and perpendicular component perpendicular to the velocity direction. is the absolute interferometric phase, and is the antenna receiving and transmitting mode. When the antenna works in monostatic mode, ; when the antenna works in bistatic mode, .
For the purpose of the three-dimensional model construction for a given ground target position , the vectors of , , and have to be totally represented by the Cartesian coordinate system. In particular, the look vector represented by the VNW system in (5) needs to be converted to the Cartesian coordinate system. As a result, (3) can be written as: Equation (7) can be further transformed as: Therefore, the three-dimensional model can be established as shown in (8).

THE CALIBRATION SENSITIVITY EQUATIONS
The calibration sensitivity equations can be built in terms of the three dimensional model shown in (8): (9) where is the error value of and is the initial value of and are the threedimensional coordinates of the master antenna phase center and target , respectively. is the transformation matrix , and is the unknown vector including the InSAR parameters and phase offset needing to be solved.
We can assume that the variables related to three-dimensional coordinates of the ground target in terms of (6) and (7), i.e., radar wavelength, position of master antenna phase center, baseline vector, velocity vector, near slant range, Doppler centroid frequency, and absolute interferometric phase. Radar wavelength is deemed as a constant for a radar system working at the designated frequency. The position of the master antenna phase center and the ground target usually have a high accuracy because they are generally measured by the use of GPS ground stations and through carrying out differential GPS processing (DGPS), which can achieve an absolute position accuracy of around 5 cm. Consequently, baseline length, baseline inclination, near slant range, Doppler centroid frequency, and absolute interferometric phase together constitute the unknown vector (10) where and represent baseline length and baseline inclination (included angle between the Z-axis and Y-axis in the Cartesian coordinate system), which jointly denote the baseline vector . is the near slant range, is the Doppler centroid frequency, and is the absolute interferometric phase containing the phase-offset value.
The unknown vector can be solved by LS adjustment using GCPs. The number of GCPs determines the number of calibration sensitivity equations. Given that the number of GCPs is n, the n calibration sensitivity equations can be linearized as follows: (11) where is the coefficient matrix, is the increment of unknown vector X, and is the initial value of when the initial is referred.
The LS solution (Zhang et al. 2010) is employed for the parameter refinement of unknown vector X. The parameter refinement is an iterative process until the location errors are less than a constant.

CALIBRATION RESULT
The proposed SAIC method was tested on real data acquired by the airborne CASMSAR system (J. X. Zhang 2012), which operates in the X-band and P-band. The experimental area, Pucheng area located in the north-west of Shanxi province, China, is a plain area covered by an amount of planted crops, paddy fields, and buildings. The terrain variation is less than 30 m. Nine CRs were uniformly deployed over the range direction within the illuminated area, as shown in Figure 2, and DGPS measurements of their positions achieved an accuracy of 3 cm. Moreover, 70 Check Points (CPs) were also measured by DGPS in the test area for precise validation of the generated InSAR DEM after calibration. The accuracy of the CPs measurements is 0.1m. In the calibration result, the corrected values for all the five variables converge through eight iterations. In the Figure 3, we can see the most of corrected values tend to be stable after three iterations. For the variables of and , their converge speeds are basically the same and their corrected values vary with a larger amplitude when the number of CRs is less than seven. It is worth stressing that the corrected value of not only contains the error itself of , but also a constant offset for the phase-offset estimation in the interferometric phase unwrapping. For the variables of , and , their corrected values vary with a smaller amplitude when they converge by iteration computation.
With regard to different numbers of CRs, we chose the final calibration result with nine CRs because the result has the more stable corrected values and the smallest 3D positioning errors on CRs, i.e., x-, y-and z-direction error. Seventy CPs were used to directly evaluate the accuracy of the obtained InSAR DEMs after completion of the calibration. Table 1  numbers of CRs. As expected, the greater number of CRs, the higher the accuracy of the DEM. However, it is found the mean and maximum of the elevation error are less than 0.1 m and 2 m, respectively, and the standard deviation of the elevation error remains in the same order of magnitude of decimeters (0.5-0.6 m). That reflects that the increase of the number of CRs has only a slight impact on the promotion of DEM accuracy.  Figure 4 shows the distribution of the 70 CPs in test area. It can be observed that the elevation of the test area gradually rises from the south to the north. The variation of the 70 CPs is within a range of 4 m, exclusive of some objects, e.g., paddy fields and buildings. Some of the CPs cover the route of the CR deployment in the south-north direction, in accordance with the near-far range direction, while some of the points are distributed in the test area along the west-east direction. The final DEM result shown in Figure 5 was obtained by the use of the calibrated variables with the smallest elevation error obtained with the CPs. We can analyze the DEM result from two perspectives. The trend of the terrain between the DEM and CPs is consistent, excluding the low coherence regions, i.e., paddy fields, for which the elevation value appears to be zero. From the point of view of the distribution of the elevation the majority of the elevation errors are within the range of 0-1 m. Exception in the elevation errors (more than 1 m) arise on two CPs; one is located on the verge of the test area, and the other is close to the buildings. However, the larger error values are reasonable since they are either far away from the position of the CRs or they are affected by ground objects. The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B1-2020, 2020 XXIV ISPRS Congress (2020 edition) We could offer two aspects of recommendations with regard to testing result in Pucheng area. In plain area, at least three GCPs are necessarily available for obtaining solvable calibration results; On the other hand, if the number of GCPs is enough, etc., more than three, one set of calibration result, with the smallest 3D positioning error, i.e., x-, y-and z-direction error decided by the CPs measured by field surveying or extracted on flat area of external DEM without objects covering, could be selected as the optimum calibration result.

CONCLUSION
A new single-pass airborne interferometric calibration (SAIC) method to correct InSAR parameters and estimate the phase offset has been presented in this paper. In the method, differing from the traditional two-dimensional InSAR geometric model, a rigorous three-dimensional model for the single-pass airborne InSAR system is first constructed. Based on the constructed model, the corrected InSAR parameters of baseline length, baseline inclination, near slant range, and Doppler centroid frequency, as well as phase offset, can be jointly solved in calibration sensitivity equations by virtue of the LS principle with a limited number of GCPs. The proposed SAIC method has been validated in two ways: the efficiency from computation converge speed and the effectiveness from the derived DEMs accuracy. It demonstrates that the proposed method has the capability of high-accuracy DEM generation using only a few GCPs in plain area, more calibration experiments in the various terrains will be tested in the future.