Feasibility of Replacing the Range Doppler Equation of Spaceborne Synthetic Aperture Radar Considering Atmospheric Propagation Delay with a Rational Polynomial Coefficient Model

Usually, the rational polynomial coefficient (RPC) model of spaceborne synthetic aperture radar (SAR) is fitted by the original range Doppler (RD) model. However, the radar signal is affected by two-way atmospheric delay, which causes measurement error in the slant range term of the RD model. In this paper, two atmospheric delay correction methods are proposed for use in terrain-independent RPC fitting: single-scene SAR imaging with a unique atmospheric delay correction parameter (plan 1) and single-scene SAR imaging with spatially varying atmospheric delay correction parameters (plan 2). The feasibility of the two methods was verified by conducting fitting experiments and geometric positioning accuracy verification of the RPC model. The experiments for the GF-3 satellite were performed by using global meteorological data, a global digital elevation model, and ground control data from several regions in China. The experimental results show that it is feasible to use plan 1 or plan 2 to correct the atmospheric delay error, no matter whether in plain, mountainous, or plateau areas. Moreover, the geometric positioning accuracy of the RPC model after correcting the atmospheric delay was improved to better than 3 m. This is of great significance for the efficient and high-precision geometric processing of spaceborne SAR images.


Introduction
Since 1999, for the IKONOS satellite, considering technical confidentiality and other factors, American space imaging companies have started to provide rational polynomial coefficient (RPC) models instead of rigorous geometry models to end users as basic imaging products [1]. The RPC model is a universal geometric model of remote sensing satellite sensors, which is usually fitted by a range Doppler (RD) model in spaceborne synthetic aperture radar (SAR) imaging systems. The difference between the RPC and RD models is that the former is a purely mathematical model, is independent of the sensor, and involves simple calculations. Furthermore, with the emergence of various imaging sensors, it is difficult for end users to add new sensor models to existing software systems to process new sensor data, while the RPC model solves this difficulty well. Moreover, from the perspective of the integral of the sum of the dry, wet, and ionosphere components of the refractivity N, from h t to the top of the atmosphere h ∞ [13,14]: where P is the total pressure (mbar), T is the temperature (K), e is the partial pressure of water vapor (mbar), and W cloud is the cloud water content (g/m 3 ). These meteorological data can be downloaded from the U.S. National Center for Environmental Prediction (NCEP) website (https: //rda.ucar.edu/datasets/ds083.2/). In addition, n e is the electronic density of the ionosphere, which can be obtained by bilinear interpolation in time and space based on the global ionospheric map provided by the European Center for Orbit Determination (CODE) (ftp://cddis.gsfc.nasa.gov/gnss/products/ionex/). Here, f and θ are the frequency and incidence angle of the radar, respectively, which can be obtained from the auxiliary file of the SAR satellite image. The constant coefficients have the following values: k 1 = 77.6 K/mbar, k 2 = 26.0 K/mbar, k 3 = 3.75310 5 K 2 /mbar, k 4 = 1.45 m 3 /g, and k 5 = 24.03310 7 m 3 /s 2 [15][16][17].

Spaceborne SAR RD Equation
The RD equation establishes the correspondence between the image and object point coordinates from the perspective of the spaceborne SAR imaging geometry. It consists of an earth model equation, SAR doppler equation, and SAR range equation [18]: where R s and → R t = [X t Y t Z t ] are the position vectors of the SAR satellite and target, respectively; → V s and → V t are the satellite and target velocity vectors, respectively; R is the slant range; λ is the radar wavelength; A = 6378.137 km is the mean equatorial radius; and B = (1 − 1/ f )A is the polar radius with a flattening factor f = 298.255.
Considering the slant range measurement error caused by the atmospheric propagation delay, the slant range can be expressed as

RPC Model
The RPC model uses the ratio polynomial to represent the correspondence between the ground coordinate D(Latitude, Longitude, Height) (the polar stereographic coordinate system is more applicable at 70-75 degrees latitude north or south) and the pixel coordinates d(line, sample) [19]: where N l (P, L, H), N s (P, L, H), D l (P, L, H), and D s (P, L, H) are all third-order polynomials. These functions have the form F(P, L, H) = a 1 + a 2 L + a 3 P + a 4 H + a 5 LP + a 6 LH + a 7 PH + a 8 L 2 + a 9 P 2 + a 10 H 2 + a 11 PLH +a 12 L 3 + a 13 LP 2 + a 14 LH 2 + a 15 L 2 P + a 16 P 3 + a 17 PH 2 + a 18 L 2 H + a 19 P 2 H + a 20 H 3 (6) where a i is the RPC model coefficient, i = 1, 2, . . . , 20; Generally, the first parameter coefficient of the denominator term is set to 1, so there are 78 parameters for RPC model. Here, (P, L, H) and (X, Y) are the normalized ground and image coordinates, respectively. These coordinates can be expressed as where is the normalized parameter of the ground coordinates and s j , l j ( j = o f f , scale) are the normalized parameters of the image coordinates. The RPC fitting results obtained based on a variety of spaceborne SAR sensors have shown that RPC models with third-order and unequal denominators have the highest replacement accuracies [7].
The following experiments were based on this model.

RPC Model Fitting and Accuracy Verification Methods
The RPC model is essentially a mathematical model. In this study, we used a terrain-independent solution algorithm to calculate the RPC parameters. The core of the algorithm involves establishing a virtual control grid based on the SAR image, DEM data of the coverage area, and forward transformation of the RD model. Then, the RPC parameters are fitted according to the least squares fitting method and the correspondence between the ground coordinates and image pixel coordinates of the virtual control points. However, every certain virtual grid point in space is affected by atmospheric propagation delay, which causes slant range measurement error. Therefore, we used two atmospheric delay correction schemes to correct the slant range values at all virtual grid points, then re-fitted the RPC model and evaluated the fitting and positioning accuracies.
The steps of the algorithm are listed in Figure 1. Step 3: Correct Atmospheric Propagation Delay Plan 1: The atmospheric delay correction value at the point in the center of the scene and the average elevation of the coverage area is used to correct the atmospheric propagation delay at all virtual points.

Δs=s'-s'' Δl=l'-l''
Output : maximum error, root mean square error of range, azimuth and plane, respectively   The steps of the algorithm can be described as follows.
(1) Firstly, the input data sets, including the SAR image auxiliary file and DEM of the study area, were prepared; (2) Then, a virtual plane grid was built according to the image size, the elevation range of the coverage area was interpolated based on the DEM data and RD model, and layering processing was performed in the elevation direction to obtain a virtual space control grid. The number of layers was required to be greater than three to prevent the design matrix from becoming ill-conditioned [1]. Next, a virtual space check grid was constructed by interpolating between the centers of four adjacent virtual control points and the two adjacent elevation layers; (3) The atmospheric propagation delay values were corrected at all virtual grid points using two plans: one in which the atmospheric delay correction value at the point at the center of the scene and the average elevation of the coverage area were employed to correct the atmospheric propagation delay at all virtual points, and one in which the atmospheric delay correction value at each virtual grid point was used to correct its own atmospheric delay; (4) The RPC model was fitted separately using the original RD model, RD model modified using atmospheric delay correction plan 1, and RD model modified using atmospheric delay correction plan 2, and then the fitting accuracy was evaluated. A detailed flowchart of the RPC model fitting accuracy evaluation is shown in Figure 2. The steps of the algorithm can be described as follows.
(1) Firstly, the input data sets, including the SAR image auxiliary file and DEM of the study area, were prepared; (2) Then, a virtual plane grid was built according to the image size, the elevation range of the coverage area was interpolated based on the DEM data and RD model, and layering processing was performed in the elevation direction to obtain a virtual space control grid. The number of layers was required to be greater than three to prevent the design matrix from becoming illconditioned [1]. Next, a virtual space check grid was constructed by interpolating between the centers of four adjacent virtual control points and the two adjacent elevation layers; (3) The atmospheric propagation delay values were corrected at all virtual grid points using two plans: one in which the atmospheric delay correction value at the point at the center of the scene and the average elevation of the coverage area were employed to correct the atmospheric propagation delay at all virtual points, and one in which the atmospheric delay correction value at each virtual grid point was used to correct its own atmospheric delay; (4) The RPC model was fitted separately using the original RD model, RD model modified using atmospheric delay correction plan 1, and RD model modified using atmospheric delay correction plan 2, and then the fitting accuracy was evaluated. A detailed flowchart of the RPC model fitting accuracy evaluation is shown in Figure 2.

Experimental Data
In this study, GF-3 images of four imaging modes, namely Spot-Light (1 m resolution), Ultra-Fine-Strip (3 m resolution), Fine-Strip-I (5 m resolution), and Fine-Strip-II (10 m resolution), were used as experimental data. To verify the positioning accuracy of the RPC model considering the atmospheric propagation delay, GF-3 satellite images and GCP data from Anping, Zhanjiang, Xianning, Zhangye, and Inner Mongolia in China were used. The GCPs were derived from global navigation satellite system (GNSS) receivers, which were selected in locations with obvious features, such as the intersections of pathways. Centimeter-level positioning accuracy could be achieved using real-time kinematic static observation technology. The details of the GF-3 satellite images are provided in Table 1, and the distribution of the experimental data and GCPs is depicted in Figure 3.

Experimental Data
In this study, GF-3 images of four imaging modes, namely Spot-Light (1 m resolution), Ultra-Fine-Strip (3 m resolution), Fine-Strip-I (5 m resolution), and Fine-Strip-II (10 m resolution), were used as experimental data. To verify the positioning accuracy of the RPC model considering the atmospheric propagation delay, GF-3 satellite images and GCP data from Anping, Zhanjiang, Xianning, Zhangye, and Inner Mongolia in China were used. The GCPs were derived from global navigation satellite system (GNSS) receivers, which were selected in locations with obvious features, such as the intersections of pathways. Centimeter-level positioning accuracy could be achieved using real-time kinematic static observation technology. The details of the GF-3 satellite images are provided in Table 1, and the distribution of the experimental data and GCPs is depicted in Figure 3.

Experimental Results and Analysis
To evaluate the feasibility of replacing the range Doppler equation of spaceborne SAR by considering the atmospheric delay with the RPC model, the following three experiments were conducted. (1) For the terrain-independent calculation algorithm of the RPC model of spaceborne SAR, the effects of the virtual grid pattern, number of elevation layers, and orbit fitting order on the fitting accuracy of the RPC model were analyzed. (2) Based on the appropriate grid pattern, number of elevation layers, and orbital fitting order, the fitting accuracy of the RPC model was assessed using the two proposed atmospheric delay correction plans. (3) Based on the fitting results of the RPC models, with and without considering the atmospheric delay of the radar signal, combined with the GCP data in the image coverage area, the positioning accuracy of the RPC model considering the atmospheric delay was verified.

Impact of RPC Model Fitting Parameter
(1) Grid Size Assuming there to be five elevation layers, six different grid styles were set (style 1: 4000 pixels × 4000 pixels; style 2: 2000 pixels × 2000 pixels; style 3: 1000 pixels × 1000 pixels; style 4: 500 pixels × 500 pixels; style 5: 200 pixels × 200 pixels; style 6: 100 pixels × 100 pixels), here "N pixels × N pixels" refers to a grid spacing of N pixels, or a grid size. We take N pixels as the sampling interval in the range direction and azimuth direction, and take the sampling points as the virtual control points. Based on the fine-strip product data of the GF-3 satellite, DEM data with 30 m resolution, and the RPC fitting method described in Section 3, the influence of the grid size on the fitting accuracy of the RPC model was analyzed. The experimental results are shown in Figure 4a. (

2) Number of Elevation Layers
Using the appropriate grid pattern identified in the previous step, the number of different elevation layers (two, three, four, five, or six) was set to evaluate the influence of the number of elevation layers on the fitting accuracy of the RPC model. The experimental results are presented in Figure 4b.

(3) Orbit Fitting Order
With the appropriate grid pattern and number of elevation layers, different orders of satellite orbital fitting (third, fourth, or fifth order) were set to analyze the influence of the orbit fitting order on the fitting accuracy of the RPC model. The experimental results are depicted in Figure 4c. Figure 4 shows that as the control grid becomes smaller and the number of elevation layers increases, the plane RMSE of the checkpoint gradually decreases. Eventually, the accuracy of the RPC model becomes close to that of the rigorous imaging geometry model. As the orbital parameter fitting order increases, the plane RMSE of the checkpoint gradually increases. To balance the fitting accuracy and computational efficiency of the RPC model, we set the grid size to 500 pixels × 500 pixels, the number of elevation layers to five, and the orbit fitting order to third-order for the subsequent study.

Experimental Results and Analysis
To evaluate the feasibility of replacing the range Doppler equation of spaceborne SAR by considering the atmospheric delay with the RPC model, the following three experiments were conducted. (1) For the terrain-independent calculation algorithm of the RPC model of spaceborne SAR, the effects of the virtual grid pattern, number of elevation layers, and orbit fitting order on the fitting accuracy of the RPC model were analyzed. (2) Based on the appropriate grid pattern, number of elevation layers, and orbital fitting order, the fitting accuracy of the RPC model was assessed using the two proposed atmospheric delay correction plans. (3) Based on the fitting results of the RPC models, with and without considering the atmospheric delay of the radar signal, combined with the GCP data in the image coverage area, the positioning accuracy of the RPC model considering the atmospheric delay was verified.

Impact of RPC Model Fitting Parameter
(1) Grid Size Assuming there to be five elevation layers, six different grid styles were set (style 1: 4000 pixels × 4000 pixels; style 2: 2000 pixels × 2000 pixels; style 3: 1000 pixels × 1000 pixels; style 4: 500 pixels × 500 pixels; style 5: 200 pixels × 200 pixels; style 6: 100 pixels × 100 pixels), here "N pixels × N pixels " refers to a grid spacing of N pixels, or a grid size. We take N pixels as the sampling interval in the range direction and azimuth direction, and take the sampling points as the virtual control points. Based on the fine-strip product data of the GF-3 satellite, DEM data with 30 m resolution, and the RPC fitting method described in Section 3, the influence of the grid size on the fitting accuracy of the RPC model was analyzed. The experimental results are shown in Figure 4a.

(2) Number of Elevation Layers
Using the appropriate grid pattern identified in the previous step, the number of different elevation layers (two, three, four, five, or six) was set to evaluate the influence of the number of elevation layers on the fitting accuracy of the RPC model. The experimental results are presented in Figure 4b.   Figure 4 shows that as the control grid becomes smaller and the number of elevation layers increases, the plane RMSE of the checkpoint gradually decreases. Eventually, the accuracy of the RPC model becomes close to that of the rigorous imaging geometry model. As the orbital parameter fitting order increases, the plane RMSE of the checkpoint gradually increases. To balance the fitting accuracy

RPC Model Fitting Accuracy Evaluation
Based on the atmospheric propagation delay correction model described in this paper combined with the meteorological data stored by the NCEP every 6 h, the global ionospheric map provided by CODE every 1 h, and the average elevation and size of the image coverage area provided by the GF-3 image auxiliary file, the atmospheric propagation delay correction of the radar signal was calculated to analyze the influence of atmospheric propagation delay correction on the RPC fitting accuracy.
In this study, atmospheric propagation delay correction of the radar signal was performed using two methods: one with a GF-3 image of a single scene and unique atmospheric delay correction parameters (plan 1) and one with an image of a single scene and spatially varying atmospheric delay correction parameters (plan 2). The two plans can be described as follows. In plan 1, the atmospheric delay correction at the point at the center of the scene and the average elevation of the coverage area are calculated to correct the atmospheric propagation delay at all virtual grid points. In plan 2, the atmospheric propagation delay correction value is calculated at each virtual grid point to correct its own atmospheric delay. The delay results are presented in Table 2. In Table 2, the calculation results for plan 1 show that the maximum slant range correction value at the center point among the different images is −3.795 m, and the delay error in the ground range direction was calculated based on the incident angle of the image as approximately −5 m, which indicates that the atmospheric delay error has a significant influence on the positioning accuracy. The calculation results for plan 2 show that the atmospheric delay correction values of the virtual grid points with different spatial distributions on a scene image are different, and the slant range correction values of the virtual grid points in all of the image scenes are between 1.521 m and −3.919 m. Thus, the atmospheric delay correction values differ between imaging angles and imaging regions. It was, therefore, concluded that the atmospheric propagation delay is mainly affected by the radar incidence angle and the topography of the study area.
To analyze the main influencing factors of the atmospheric propagation delay in the SAR image scenes further, we visualized the atmospheric propagation delay correction values of the virtual grid points calculated using plan 2, assuming this to be 0 pixels in the azimuth direction and sampling every 500 pixels in the range direction, with five elevation layers. The trend of the atmospheric delay correction value in the range direction is shown in Figure 5. To analyze the influence of the elevation on the atmospheric propagation delay, we increased the number of elevation layers of the virtual control grid to 300 separately in the plains and mountains area (ZJ-5851) and the plateau area (ZY-7668). The maximum variation ∆ϕ of the corresponding incident angle at the same sampling point in each elevation layer in the experimental region is shown in Figure 6, where the orbital height of the satellite is known to be approximately 755 km. Here, ∆ϕ was approximated using Equation (9), and was found to be less than 0.03 • in the ZJ-5851 image coverage area and less than 0.15 • in the ZY-7668 image coverage area. Thus, the variation of the incident angle in the elevation direction in the study area is negligible, and the same sampling point can be selected in each elevation layer to analyze the atmospheric propagation delay law as a function of elevation, as shown in Figure 7.
where H s is the orbital height of the GF-3 satellite; H max and H min are the maximum and minimum elevations of the image coverage area, respectively.   Figure 6, where the orbital height of the satellite is known to be approximately 755 km. Here, Δφ was approximated using Equation (9), and was found to be less than 0.03° in the ZJ-5851 image coverage area and less than 0.15° in the ZY-7668 image coverage area. Thus, the variation of the incident angle in the elevation direction in the study area is negligible, and the same sampling point can be selected in each elevation layer to analyze the atmospheric propagation delay law as a function of elevation, as shown in Figure 7.  Some researchers have demonstrated that the elevation-atmospheric delay relationship is usually a linear or exponential function [20,21]. Figure 7a shows that in the plains and mountains area, the elevation-atmospheric delay relationship is approximately linear. Meanwhile, Figure 7b reveals that in the plateau region, the elevation-atmospheric delay relationship is also approximately linear.  Some researchers have demonstrated that the elevation-atmospheric delay relationship is usually a linear or exponential function [20,21]. Figure 7a shows that in the plains and mountains area, the elevation-atmospheric delay relationship is approximately linear. Meanwhile, Figure 7b reveals that in the plateau region, the elevation-atmospheric delay relationship is also approximately linear.  According to the RPC model fitting parameter impact analysis results in Section 4.2.1, we established a control grid by setting the grid size to 500 pixels × 500 pixels and the number of elevation layers to five. The checkpoint was determined based on the centers of the four adjacent control points and two adjacent elevation layers, and the fitting accuracy of the RPC model was evaluated after correcting the atmospheric delay grid by grid. The results are presented in Table 3.
In general, as the plane error of the checkpoint does not exceed 5% of the pixels, we believe that the RPC model can be used instead of the RD model for photogrammetric processing of spaceborne SAR satellites. Table 3 shows that in the same scene, both the uncorrected atmospheric delay and atmospheric delay correction with plan 1 yield exactly the same RPC model fitting accuracy, and that the maximum plane RMSE of the checkpoint is 0.00356 pixels. For the atmospheric delay correction with plan 2, the fitting accuracy of the RPC model is very approximate, so both approaches achieve precision better than 1% of a pixel. This indicates that the RPC model still provides high accuracy for the RD model after correcting the atmospheric delay using two plans, in plain, mountainous, and plateau areas.
The reason for obtaining these experimental results can be confirmed by research that has shown that the third-order RPC model can better fit a function whose highest power is no more than 5 [22]. For a single-scene SAR image, plan 1 uses a unique atmospheric delay correction parameter to correct the atmospheric propagation delay at all virtual grid points, which is equivalent to adding the atmospheric propagation delay as a constant function to the slant range (R) term of the RD model. Therefore, the fitting accuracy of the RPC model is unchanged. Meanwhile, plan 2 corrects the atmospheric propagation delay at all virtual grid points with the correction parameters of the corresponding virtual grid points, which is equivalent to adding the atmospheric delay correction model to the R term with the elevation as the dependent variable. Figure 7a shows that the elevation-atmospheric delay function is approximately linear in the plains and mountains area, so it can be correctly fitted with a third-order RPC model. Meanwhile, Figure 7b shows that in the plateau area, the elevation-atmospheric delay function also shows an approximate linear relationship, which indicates that RPC model can not only fit complex terrain, but also fit the terrain-related atmospheric propagation delay.

RPC Model Positioning Accuracy Verification
From the perspective of RPC model application, we carried out RPC model positioning experiments to verify the effectiveness of the RPC model considering atmospheric propagation delay. Based on the experimental data of the GF-3 satellite and GCPs, geometric positioning verification experiments were performed for the four RPC models, and the results are shown in Table 4.
(1) Firstly, for the RPC model fitted using the original RD model without atmospheric delay correction, positioning accuracy verification was conducted; (2) After correcting the systematic error in the slant range direction of the GF-3 SAR satellite, geometric positioning accuracy verification of the RPC model was performed. This systematic error was obtained by geometric calibration of the spaceborne SAR [23]; (3) After correcting the systematic error in the slant range direction of the GF-3 SAR satellite and correcting the atmospheric propagation delay error using plan 1, geometric positioning accuracy verification of the RPC model was performed; (4) After correcting the systematic error in the slant range direction of the GF-3 SAR satellite and correcting the atmospheric propagation delay error using plan 2, geometric positioning accuracy verification of the RPC model was performed. Table 4 shows that the positioning accuracy of the RPC model is between 21.241 m and 26.004 m without correction of system error and atmospheric delay error. After correcting the system error, the geometric positioning error of the RPC model was between 3.064 m and 5.529 m. After correcting the atmospheric delay error using plans 1 and 2, respectively, the maximum geometric positioning error of the RPC model was found to be 2.948 m and 2.957 m, and the positioning accuracy of a single scene image could be increased by 2.7 m, such as for ZJ-5492. Thus, it is effective to improve the geometric positioning accuracy of the RPC model by correcting the atmospheric propagation delay in the process of fitting the RD model to the RPC model. Furthermore, the positioning accuracy of the four RPC models in the azimuth direction is almost unchanged, which indicates that the slant range error caused by the atmospheric propagation delay mainly affects the positioning accuracy in the range direction. Figure 8 shows that the positioning error of all GCPs in the SAR image of a single scene is approximately 2.5 m after correcting the atmospheric propagation delay using plans 1 and 2 separately, and this remains relatively stable, which indicates that the geometric positioning results of the RPC model are reliable.
It is worth noting that although both plans are feasible, it is still difficult to produce a RPC model considering atmospheric delay directly and in real time in the SAR system, because sufficiently precise atmospheric models typically have latency. We may be able to choose a general model in which the atmospheric delay is a function of the SAR incidence angle and the elevation of a given point, in order to replace the precise atmospheric model [24]. This will be our next research direction.

Conclusions
The atmospheric delays of SAR satellite radar signals can reach several meters, which is nonnegligible in high-precision geometric positioning. Thus, two methods of atmospheric delay correction during RPC fitting of an RD model were designed and tested in this study. Based on the GF-3 SAR satellite data, a global digital elevation model, global meteorological data, and ground control data from several regions in China, the feasibility of the two methods was verified by conducting RPC model fitting and positioning experiments. We found that the RPC model exhibited high substitution and positioning accuracies after using plan 1 or plan 2 to correct the atmospheric propagation delay, in both the plains and mountains area where the relative height difference is less

Conclusions
The atmospheric delays of SAR satellite radar signals can reach several meters, which is non-negligible in high-precision geometric positioning. Thus, two methods of atmospheric delay correction during RPC fitting of an RD model were designed and tested in this study. Based on the GF-3 SAR satellite data, a global digital elevation model, global meteorological data, and ground control data from several regions in China, the feasibility of the two methods was verified by conducting RPC model fitting and positioning experiments. We found that the RPC model exhibited high substitution and positioning accuracies after using plan 1 or plan 2 to correct the atmospheric propagation delay, in both the plains and mountains area where the relative height difference is less than 1000 m, and in the plateau area where the relative height difference is more than 1000 m. Therefore, it is feasible to use plan 1 or plan 2 to correct the atmospheric delay error under any terrain conditions. Moreover, the geometric positioning accuracy of the RPC model after correcting the atmospheric delay was improved to better than 3 m.