1. Introduction
The Gaofen-3 (GF-3) satellite is the first civil C-band multipolarization synthetic aperture radar (SAR) satellite in China. This satellite has 12 imaging modes, including strip, spotlight and scanning modes. The image resolution ranges from 1 to 500 m, while the observation width ranges from 10 to 650 km. GF-3 can monitor global ocean and land information during all types of weather at all times and can expand the range of Earth observations through left–right attitude maneuvering. The satellite data can be used in many fields, such asoceanography, disaster reduction, water conservancy and meteorological and mapping research, to improve rapid response capabilities. This satellite has also filled the gap in civil autonomous high-resolution multipolarization SAR remote sensing data in China [
1,
2]. Topographic mapping is an important application direction in the field of remote sensing. Digital elevation model (DEM) data play important roles in national economic construction and scientific research [
3,
4]. Precise surveying and mapping (SAM) calibration technology has become one of the most effective means of global mapping; by incorporating SAR interferometry (InSAR) to acquire interferometric phase information of the land surface, high-precision elevation information is retrieved [
5,
6]. In 2000, the National Aeronautics and Space Administration (NASA), the National Imagery and Mapping Agency (NIMA) and the German Aerospace Center (DLR) jointly conducted the Shuttle Radar Topography Mission (SRTM), successfully completing approximately 80% of the proposed terrestrial topographic mapping mission on Earth; in addition, the absolute elevation accuracy of the DEM obtained by the SRTM was 16 m [
7,
8]. In 2016, the DLR acquired global TanDEM-X DEM data with an absolute elevation accuracy of 4 m; these data were improved over other global DEM data in their accuracy, thus introducing new opportunities for global topographic mapping [
9]. In 2018, the DLR released global TanDEM-X DEM data at a 90 m grid size [
10,
11]. Therefore, obtaining global DEM data constitutes major engineering projects in the field of InSAR topographic mapping [
11]. The abilities of high-precision satellite measurements in this application are important and depend not only on the design of the satellite platform index but also on the calibration of the geometric and interferometric parameters [
6,
12]. Geometric calibration steps improve the plane positioning accuracy of SAR images by calibrating the geometric parameters, while interferometric calibration technologies improve the vertical accuracy of the interferometric results by calibrating the acquired interferometric parameters [
13,
14]; together, these calibration steps can be used to determine the accuracy level of a global DEM. At present, many research teams in China and globally have published geometric and interferometric calibration methods that are applicable for mainstream SAR satellites [
9,
13,
15].
After applying the geometric calibration method to the data obtained by the Japanese Advanced Land-Observing Satellite-Phased Array L-band SAR (ALSO-PALSAR), the geometric positioning accuracy of the strip mode reached 9.7 m [
16]. A 5.74 m plane positioning accuracy was obtained for the Canada Radarsat-2 satellite after a geometric calibration method using corner reflectors (CRs) was applied [
17]. German TerraSAR-X satellite data were subjected to high-precision geometric calibration with CRs and achieved a plane positioning result with an accuracy of 0.3 m [
18]. In 2017, the geometric calibration precision of the GF-3 SAR satellite was studied, and the positioning accuracy was found to be better than 3 m [
14]. In 2017, a multimode hybrid geometric calibration method was proposed for GF-3 images, and the geometric positioning accuracy was found to be better than 3 m [
19]. In 2018, a sparse-control multiplatform geometric calibration method was proposed for TerraSAR-X/TanDEM-X and GF-3 images. The geometric positioning accuracy of the TerraSAR-X/TanDEM-X images found using this method was better than 3 m, while that of the GF-3 images was better than 7.5 m following calibration [
15]. In 2020, a wide-area joint geometric calibration method was proposed for GF-3 SAR images, and the positioning error of the calibrated SAR images was better than 9 m [
20]. Thus, this geometric calibration method has been successfully applied in engineering and scientific research in the past.
Although interferometric calibration technologies are complex, extensive research work has been conducted on this topic both in China and across the globe, and many interferometric calibration methods have been proposed. In 1999, an airborne interferometric calibration method based on a sensitivity equation was deduced using flat-ground geometric relationships [
21]. In 2001, based on this research, an airborne interferometric calibration method was formulated based on the interferometric sensitivity equation under squint angle conditions [
22]. In 2003, a three-dimensional reconstruction (TDR) model was introduced to derive the interferometric sensitivity equation, and a CR placement algorithm was proposed to minimize the condition numbers in sensitivity matrices [
23]. In 2010, an airborne interferometric calibration model was proposed to accommodate the phase offset, baseline length and baseline angle [
24]. Then, an airborne dual-antenna SAR interferometric calibration method based on the three-dimensional information of control points was proposed [
25]. In 2011, a distributed satellite SAR interferometric altimetry model was proposed based on the TDR model [
26]. An interferometric calibration optimization method for spaceborne distributed SAR systems was proposed to solve the problem of coupling the contribution of interferometric parameters to the elevation error [
27]. In 2016, the interferometric calibration method derived based on the InSAR TDR model and sensitivity equation was verified using TanDEM-X data, thus proving the universality of this spaceborne SAR interferometric calibration method. However, in this technology verification, each error was not solved independently [
28].
Table 1 further describes the current status of research on geometric and interferometric calibration methods.
To date, geometric calibration technologies have gradually matured, and interferometric calibration technologies have developed from airborne to spaceborne systems. Although extensive research has been conducted on this topic, most analyzed parameter errors exhibit strong coupling effects during the interference process, and no joint solution has been accurately distinguished, leading to inaccurate interferometric parameter corrections [
28]. To solve this problem, a precise SAM calibration method is proposed herein based on independent parameter decomposition (IPD). We established a geometric calibration model and resolved the initial imaging time error as well as the initial slant range error to improve the location accuracy. On this basis, the interferometric phase error and baseline vector error were precisely estimated to ensure the vertical accuracy of the interferometric results by establishing the interferometric calibration model. The accuracy and reliability of the proposed method were verified by processing GF-3 SAR data obtained on the same orbit in six scenes, and high-accuracy DEM products were obtained, thus providing a technical reference for the global mapping missions of domestic SAR satellites.
2. Geometric Calibration Model
Due to the measurement errors of satellite orbit and imaging parameters, the plane positioning accuracy of SAR images is always restricted. Geometric calibration technologies are the key to achieving high-precision SAR image positioning information and are also the main method by which the generation of high-precision DEMs can be ensured. The Range–Doppler (R-D) model, the mainstream positioning model used for SAR geometry processing, is a rigorous model that conforms to the SAR imaging mechanism and describes the geometric relationships between satellite sensors and object points in the geocentric coordinate system. The R-D model is a nonlinear system of equations composed of the slant range equation, Doppler equation and Earth ellipsoid equation; these three equations can be expressed as follows [
15]:
where
is the distance between the main satellite sensor and the ground object,
is the position vector of the satellite sensor,
is the position vector of the ground object,
is the Doppler frequency,
is the radar wavelength,
is the velocity vector of the satellite sensor and
is the velocity vector of the ground object. Because the ground object point is located on the Earth ellipsoid, the velocity vector of the ground object point in the geocentric coordinate system is 0; that is,
.
and
are the major semiaxes and minor semiaxes of the Earth reference ellipsoid, respectively, and
is the elevation of the ground object.
The position and velocity vectors of the satellite sensor can be described as functions of time using cubic polynomials, as shown in Formula (4).
According to the R-D positioning model, the variables that are necessary for geometric positioning include the slant range (R), satellite sensor position vector (), satellite sensor velocity vector (), Doppler frequency (), radar wavelength () and ellipsoid parameters (). Among them, the radar wavelength is a known value that can be obtained according to the formula , where c is the speed of light and f is the frequency. The ellipsoid parameters can also be obtained after the ellipsoid model is selected.
The slant range can be obtained by measuring the temporal difference between the transmission and reception of the radar pulse; the slant range error causes the position of the field of view of the target to move along the Doppler line, thus resulting in a positioning error. In the geometric positioning model, the slant range can be expressed as a linear function of the image column number as follows:
where
is the initial slant range,
is the range pixel spacing and
is the range pixel coordinate. In Equation (5), only the effect of the initial slant range error is considered when assessing the positioning accuracy in the range direction.
The linear function between the azimuth direction time and image line number can be expressed as follows:
where
is the initial imaging time,
is the azimuth time spacing and
is the line number of the image along the azimuth. Equation (6) considers only the effect of the initial imaging time error on the positioning accuracy in the azimuth direction.
In summary, we consider only the impacts of the initial slant range and imaging time on the plane positioning accuracy, so the geometric calibration parameters of interest in this paper are . Then, the geometric calibration model is established according to the R-D model, and the geometric calibration parameters are estimated.
First, a conditional equation based on the slant range equation is established as follows:
Next, a conditional equation based on the Doppler frequency is established as follows:
Then, the geometric calibration parameters are estimated using the nonlinear least squares method by considering several calibration points derived based on the conditional equation. The calibration parameter
is used as an unknown observation value, and the error equation is obtained by linearizing the observation value derivation; these equations can be expressed as follows:
where the true observation value is
. The observation errors
and
can be calculated by substituting approximate values of each undetermined geometric calibration parameter
into Equations (7) and (8), respectively. The error equation can be rewritten as a matrix and expressed as follows:
where
,
,
and
.
The coefficients in matrix
B are derived from the partial derivatives of Equations (7) and (8) to obtain the geometric calibration parameter
, which can be expressed as follows:
where
,
, and
Finally, the normal equation can be established based on the initial value (
) of the geometric parameter.
The geometric correction parameter can then be corrected based on the normal equation; the result can be written as follows:
The first geometric calibration parameter estimation can be obtained as follows:
The geometric correction parameter solution is obtained in a gradual iterative approach process. First, the initial values of calibration parameters and are obtained; then, the error equation is established according to these initial values, and the geometric calibration parameter is corrected using the normal equation. Next, the first approximate value of the calibration parameter is obtained, and this value is taken as the initial value in the next iteration. When the calibration parameter gradually stabilizes after many iterations and the parameter correction is less than the specified threshold (), the iteration is said to have converged. At this time, the cumulative parameter correction result is considered to be the final error in the geometric calibration parameters.
3. Three-Dimensional Reconstruction Model
The TDR is the basic model applied to spaceborne InSAR. Its basic principle involves obtaining the three-dimensional coordinates of GCPs by combining the slant range equation and Doppler equation in the R-D positioning model through the interferometric equation [
29].
Herein, we describe the geometric relationship of the InSAR data in geocentric space rectangular coordinates (O-XYZ), as shown in
Figure 1, where
and
are satellites;
and
are the position vectors of the satellite to the GCP;
is the position vector of the GCP;
and
are the slant ranges of
and
, respectively;
is the velocity vector;
is the Doppler frequency;
is the radar wavelength;
is the absolute phase;
is the wrapped phase;
is the baseline vector; and
is the antenna receiving mode. The Q values of 2 and 1 represent the repeated-orbit interferometric mode and the double-antenna interferometric mode, respectively.
We mainly focus on the cross-track interferometry (XTI) mode of the spaceborne InSAR. The slant range equation, the Doppler equation in the R-D model and the interferometric phase equation are combined to establish the TDR model, and the geometric relationships can be expressed as follows:
According to the geometric relationships of the interferometric SAR, the GCP can be expressed in O-XYZ coordinates as follows:
where
is the look vector and
is the unit look vector. In Equation (18), the look vector decomposition (LVD) method is generally used to obtain the GCP coordinates. In this work, we convert the look vector to the unit look vector during the calculation [
29], and the unit look vector needs to be represented in the Madsen moving coordinate (MMC) system [
30]. The MMC system is based on applying the phase center of the main antenna as the origin of the coordinate system; the velocity vector direction is set to be equal to the V axis; the velocity and baseline cross product direction are set to be equal to the W axis; and V, N and W are set to meet the requirements of the right-hand coordinate system. To more conveniently describe the system, we refer to the MMC system as VNW in this work. The terms
,
and
are considered components of the unit look vector on the V, N and W axes, respectively, and
can be written as follows:
In the LVD method, the look vector is converted from VNW coordinates to O-XYZ coordinates to obtain the three-dimensional coordinates of the GCPs, as shown in Equation (20) below:
where
is the transformation matrix,
is the baseline component in the velocity direction and
is the baseline component in the direction perpendicular to the velocity.
According to the TDR expression, the main factors that affect the precision of interferometric SAR data are the interferometric phase error and the baseline vector error. Due to coupling among interferometric parameter error sources, to improve the accuracy of the interferometric parameter calculations during the calibration process, we adopted the IPD method, thus allowing us to precisely estimate the interferometric phase error and baseline vector error, improve the vertical accuracy of interferometric results and ultimately obtain high-precision DEM products.
4. Interferometric Calibration Model
The method utilized in this paper is based on the transmission characteristics of the interferometric phase error and baseline vector error to the elevation error; these characteristics are used herein to establish the interferometric calibration model. In spaceborne InSAR topographic mapping, interferometric phase errors generally consist of random and systematic errors. The random errors are caused by incoherent sources in the InSAR system, and the corresponding phases are randomly distributed; these errors are not involved in the calibration process. The systematic errors consist of phase drift, phase synchronization, temperature and absolute phase deviation errors [
12,
27]. Because the phase drift, phase synchronization and temperature errors are not the main sources of interferometric phase errors, they are omitted from the calibration process in this work. In a spaceborne InSAR system, when the interferometric phase is filtered and unwrapped, fixed deviations still exist among the flat phase, the unwrapped phase and the absolute interferometric phase; this fixed deviation can be considered the absolute phase deviation [
31]. Therefore, the main component of the interferometric phase error is the absolute phase deviation, calculated as the interferometric parameter in this paper.
First, the absolute interferometric phase of each point was obtained from the GCPs. Then, the unwrapped phase and flat phase were subtracted from the absolute phase. Finally, the fixed deviation of the phase data of the SAR image was compensated based on the mean estimation method. In this way, the interferometric phase error of each image was obtained as follows:
where
is the absolute interferometric phase,
is the unwrapped phase,
is the flat phase and
is the interferometric phase error.
In spaceborne InSAR topographic mapping, the baseline vector error directly affects many factors, such as the elevation ambiguity and the transmission coefficient of the elevation error. The baseline is a three-dimensional vector, and to more conveniently express the baseline parameters, the track–cross–normal (TCN) coordinate system is generally adopted. In this system, the master antenna phase center is set as the origin of the coordinates, and the N axis is derived from the main antenna phase center in the direction pointing toward the geocenter; the cross-product direction between the N axis and the velocity vector is considered the positive direction of the C axis, and T, C and N are set to meet the requirements of the right-hand coordinate system [
28]. The baseline vector is then decomposed into the initial baseline and time-related baseline rate in the TCN system. Considering that all calculations are projected to the zero-Doppler plane of the main image during the XTI process conducted in this paper, the baseline along-track direction is always zero (
) [
32]. The baseline vectors are calculated as follows:
where
,
and
are the baseline three-dimensional (3D) vectors of the TCN coordinate system;
and
are the initial baseline and baseline rate of baseline
in the C direction; and
and
are the initial baseline and baseline rate of baseline
in the N direction.
In this paper, when the baseline vector parameters are calibrated, we first start from the TDR model expression of the InSAR data; then, the sensitivity equation is established. Finally, we use the adjustment model to correct the interferometric parameters based on the least-squares criterion. To improve the accuracy of the baseline calculation, we adopt the iteration method. The root mean square error (RMSE) of the coordinate difference derived before and after the iteration is defined as the objective function. When the objective function is less than a specified threshold (
in this paper), iterative convergence occurs. Finally, high-precision baseline 3D vectors are obtained to ensure the generation of high-precision DEM products. The process by which the baseline vectors are precisely estimated is shown in
Figure 2 [
13].
First, TDR equations representing the three coordinate axis directions are derived from the spaceborne InSAR model; these equations can be expressed as follows:
where
are the coordinates of the master antenna phase center,
are the coordinates of the GCP,
is the slant range between the phase center of the master antenna and the GCP and
is the transformation matrix.
Second, the partial derivative of the baseline vector is calculated based on the TDR equation, and the interferometric sensitivity equation is established; this equation can be expressed as follows:
The error equation is then established as follows:
where
and
denote the correction of the initial baseline and the baseline rate in the C direction, respectively, and
and
denote the correction of the initial baseline and the baseline rate in the N direction, respectively. The observation error can be calculated as follows:
where
are the three-dimensional reconstruction coordinates of the targets,
are the coordinates of the known GCP, and
is the true observed value.
In the interferometric calibration process, the observation accuracy usually varies due to different observation conditions; thus, the introduced errors are often different. The coherence coefficient is an important criterion used to evaluate the quality of SAR images. Therefore, in this work, the coherence coefficients of the GCPs are weighted, and the objective function is expressed as follows:
where
is the diagonal matrix. The diagonal matrix can be written as follows:
where
is the coherence coefficient of each GCP.
Finally, given the initial values of the baseline vector parameters, the normal equation can be established as follows:
The baseline vector was corrected based on the normal equation, which could be denoted as follows:
The first estimation of the baseline vector was obtained as follows:
6. Conclusions
In this paper, we proposed an IPD method to support the application of SAM approaches using GF-3 SAR data. In the proposed IPD method, geometric and interferometric calibration models are separately utilized to prevent parameter coupling phenomena from arising. We decomposed the sources of geometric and interferometric parameter errors in the GF-3 SAR system by establishing IPD models, and the derived errors could also be considered the main factors affecting the resulting DEM products. Herein, the plane positioning accuracies of the SAR images, the vertical accuracies of the interferometric results and the spatiotemporal distribution characteristics and elevation accuracies of DEM products were analyzed. The correctness, reliability and practicability of the proposed method were verified. The results show that the average plane positioning accuracy of the SAR images was 5.09 m following the geometric calibration process, and the average post-calibration vertical accuracy of the interferometric results was 4.18 m. The corresponding average absolute elevation accuracy of the DEM products was better than 3.09 m. These results confirm that the method proposed herein can improve the elevation accuracies of Gaofen-3 DEM products, which met the 1:50000-scale topographic map SAM precision, or even higher precisions, in the flat areas of China. Moreover, there are still some problems to be further studied in this paper. Due to the lack of high-precision CRs, the method proposed cannot fully verify the on-orbit calibration of satellites. Follow-up research will mainly focus on this aspect.
The follow-up satellites to GF-3, i.e., the 1 m C-SAR 01 and 02, were launched on 23 November 2021 and on 7 April 2022. This ocean-monitoring SAR satellite constellation consists of three satellites. The 1 m C-SAR satellites were designed in consideration of their interferometric capabilities. An increasing number of interferograms can be provided by this new satellite constellation, and these novel interferometric data can also be obtained more easily. We hope that the proposed method will be useful in supporting the usage of these new satellites.