Bias Compensation for Rational Polynomial Coefficients of High-Resolution Satellite Imagery by Local Polynomial Modeling

Xiang Shen 1,2,3, Qingquan Li 1,*, Guofeng Wu 1 and Jiasong Zhu 1 1 Key Laboratory for Geo-Environmental Monitoring of Coastal Zone of the National Administration of Surveying, Mapping and Geo-Information & Shenzhen Key Laboratory of Spatial-temporal Smart Sensing and Services, Shenzhen University, Nanhai Road 3688, Shenzhen 518060, China; shen@szu.edu.cn (X.S.); guofeng.wu@szu.edu.cn (G.W.); zhujiasong@gmail.com (J.Z.) 2 College of Information Engineering, Shenzhen University, Nanhai Road 3688, Shenzhen 518060, China 3 Beijing Key Laboratory of Urban Spatial Information Engineering, Beijing 100038, China * Correspondence: liqq@szu.edu.cn; Tel.: +86-755-2697-9741


Introduction
High-resolution satellite images (with a spatial resolution of several meters or finer) have been widely used in earth and environmental studies [1]. Limited by the unsatisfactory performance of onboard navigation sensors, the direct-georeferencing accuracy of the satellite imagery is commonly worse than several times of the Ground Sampling Distance (GSD) [2,3], which cannot meet the requirement of high-precision remote-sensing applications such as urban 3-D reconstruction [4] and orthophoto production [5]. To improve the geopositioning accuracy is one of the most critical issues in the photogrammetric processing chain of high-resolution satellite imagery, and it has been continuously studied for years [6][7][8].

The Rational Function Model
The rational function model used for satellite image georeferencing employs a ratio of two cubic polynomials to approximate rigorous imaging geometry [23,35].
with l n = (l − l 0 )/l S s n = (s − s 0 )/s S (2) and where (l, s) are the line and sample components of the image coordinate, respectively; (λ, ϕ, h) are the geodetic longitude, latitude, and height components of the ground coordinate, respectively; l n , s n , U, V, W are the normalized coordinates, l 0 , s 0 , λ 0 , ϕ 0 , h 0 and l S , s S , λ S , ϕ S , h S are the corresponding Remote Sens. 2017, 9,200 3 of 12 offset and scale terms, respectively [10]; and P i (i = 1, . . . , 4) is a third-order polynomial and has the following form.
P(U, V, W) = c 1 + c 2 U + c 3 V + c 4 W + c 5 UV +c 6 UW + c 7 VW + c 8 U 2 + c 9 V 2 + c 10 W 2 +c 11 UVW + c 12 U 3 + a 13 UV 2 + c 14 UW 2 + c 15 U 2 V +c 16 To eliminate unnecessary degrees of freedom, the zero-order terms of the denominator polynomials (i.e., P 2 and P 4 ) in Equation (1) are always fixed to be 1, and accordingly, there is a total of 78 independent coefficients in the RFM [35].
Because of the imperfect measurement of satellite ephemeris and attitude, the vendor-provided RPCs are commonly inaccurate and Equation (1) is no longer valid. The most common solution is to adjust the image coordinates to keep the equation balanced again [3,10,30]. Substituting Equation (2) into Equation (1) and adding the correction terms (∆l, ∆s) in the image space, it yields l = l + ∆l = P 1 (U,V,W) P 2 (U,V,W) · l S + l 0 s = s + ∆s = P 3 (U,V,W) P 4 (U,V,W) · s S + s 0 (5) where ∆l and ∆s are usually modeled as functions of the image coordinates (l, s).

Global Polynomial Models
In previous studies [3,23,26,28,30] the correction terms in Equation (6) are modeled as global polynomials of measured image coordinates and have the following general form. ∆l = a 0 + a 1 l + a 2 s + a 3 l 2 + a 4 ls + a 5 s 2 + . . .
More specifically, there are mainly four schemes with different independent variables being widely used [3,10].
(a 0 , a 1 , b 0 , b 1 ), which model shift and time-dependent errors (the global shift and drift model); 3.
Some researchers have reported that the affine model and the quadratic model usually produce better performance than the other two models [3]. Less sophisticated models are not routinely recommended for use unless the number of GCPs is very limited, primarily due to the lack of capacity to absorb radial distortions and some more complex errors.
with W the diagonal weight matrix Based on the recommendations of the literature [36,37], a tri-cube kernel is used in the weight function. 3 3 for with d i the Euclidean distance where (l i , s i ) refer to the image coordinates of the i-th nearest GCP around l p , s p , and h is the bandwidth of the kernel function and can be chosen by cross-validation [38]. As illustrated in Figure 1, the weights are zero for the GCPs far away from the image point p, and the local polynomial coefficients are only determined by the GCPs nearby. Benefiting from the introduction of the local weighting strategy, local polynomial models can compensate complex, non-rigid errors in RPCs, while global polynomial models treat all GCPs as equally weighted observations and can only describe rigid deformations. Given the previous studies on the global-polynomial-based models (cf. Section 3.1) for the bias correction of RPCs, the following two schemes of the local polynomial model are employed and tested in our experiments.
(α 0 , α 1 , . . . , α 5 , β 0 , β 1 , . . . , β 5 ), which describe a second-order polynomial transformation in a local area (the local quadratic model).  Given the previous studies on the global-polynomial-based models (cf. Section 3.1) for the bias correction of RPCs, the following two schemes of the local polynomial model are employed and tested in our experiments.

Data
A stereo triplet of Chinese ZY-3 satellite images, which was acquired on 23 February 2012 and was captured from different views (nadir, backward, and forward looking directions), was used in the experiments. The study area covers a region of about 2500 km 2 and is located in Beijing, China. As shown in Figure 2, the topography of the scene ranges from relatively flat urban areas with elevations near sea level in the eastern region, to undulating mountainous areas with the elevations up to around 1500 m in the western and northwestern regions. A total of 30 ground points was surveyed by differential GPS (Global Positioning System) with an average accuracy better than 5 cm, and the corresponding image coordinates were measured manually. A Leave-One-Out Cross-Validation (LOOCV) procedure [39] was employed to detect potential outliers. We adopted the following index.
where i e refers to the image-projection error when using the i-th ground point as the check point and all other ones as control points, N is the total number of ground points. Given that the values of e R were not larger than 3.0 in all test scenarios (cf. Table 1), no ground points were identified as

Data
A stereo triplet of Chinese ZY-3 satellite images, which was acquired on 23 February 2012 and was captured from different views (nadir, backward, and forward looking directions), was used in the experiments. The study area covers a region of about 2500 km 2 and is located in Beijing, China. As shown in Figure 2, the topography of the scene ranges from relatively flat urban areas with elevations near sea level in the eastern region, to undulating mountainous areas with the elevations up to around 1500 m in the western and northwestern regions. A total of 30 ground points was surveyed by differential GPS (Global Positioning System) with an average accuracy better than 5 cm, and the corresponding image coordinates were measured manually. A Leave-One-Out Cross-Validation (LOOCV) procedure [39] was employed to detect potential outliers. We adopted the following index.
where e i refers to the image-projection error when using the i-th ground point as the check point and all other ones as control points, N is the total number of ground points. Given that the values of R e were not larger than 3.0 in all test scenarios (cf. Table 1), no ground points were identified as outliers.
In most parts of the experiments, 15 points were used as Ground Control Points (GCPs), and the other 15 points were reserved as Independent Check Points (ICPs). Figure 3 shows the image-projection errors of the ground survey points by using the vendor-provided RPCs in the georeferencing. It can be seen that the errors in the three images are of the same magnitude and have a similar pattern, this is because the maximum difference between the acquisition times of the three images is only tens of seconds and the long-term variations of orbit errors are therefore almost the same. the other 15 points were reserved as Independent Check Points (ICPs). Figure 3 shows the image-projection errors of the ground survey points by using the vendor-provided RPCs in the georeferencing. It can be seen that the errors in the three images are of the same magnitude and have a similar pattern, this is because the maximum difference between the acquisition times of the three images is only tens of seconds and the long-term variations of orbit errors are therefore almost the same.    the other 15 points were reserved as Independent Check Points (ICPs). Figure 3 shows the image-projection errors of the ground survey points by using the vendor-provided RPCs in the georeferencing. It can be seen that the errors in the three images are of the same magnitude and have a similar pattern, this is because the maximum difference between the acquisition times of the three images is only tens of seconds and the long-term variations of orbit errors are therefore almost the same.

Schemes
Four image-based bias-correction models were tested and compared in the experiments, and they are the local affine model, the local quadratic model and their corresponding global polynomial versions. The four GCPs located at the corners of the scene (i.e., No. 1-4 in Figure 3) were used as the starting configuration, and then we gradually increased the number of GCPs in the test until all 15 points were added. In all test scenarios, the 15 fixed ICPs (i.e., the unnumbered points in Figure 3) were employed to evaluate the accuracies of different bias-correction models. Figure 4 shows the tendency of the residuals with an increasing number of GCPs. When the number of GCPs in use was less than the minimum requirement of a bias-correction model, the process results would be undefined and are not presented in the figure. In our algorithms, the minimum numbers of GCPs needed in the local polynomial models are somewhat larger than those in the global polynomial models (5 vs. 3 and 8 vs. 6 for the affine model and the quadratic model, respectively).

Results
Four image-based bias-correction models were tested and compared in the experiments, and they are the local affine model, the local quadratic model and their corresponding global polynomial versions. The four GCPs located at the corners of the scene (i.e., No. 1-4 in Figure 3) were used as the starting configuration, and then we gradually increased the number of GCPs in the test until all 15 points were added. In all test scenarios, the 15 fixed ICPs (i.e., the unnumbered points in Figure 3) were employed to evaluate the accuracies of different bias-correction models. Figure 4 shows the tendency of the residuals with an increasing number of GCPs. When the number of GCPs in use was less than the minimum requirement of a bias-correction model, the process results would be undefined and are not presented in the figure. In our algorithms, the minimum numbers of GCPs needed in the local polynomial models are somewhat larger than those in the global polynomial models (5 vs. 3 and 8 vs. 6 for the affine model and the quadratic model, respectively).  It can be found from Figure 4a-c that there was not a global quadratic trend in the image-projection errors, because the performance of the global second-order polynomial model was commonly worse or only marginally better than that of the global affine model. The correction residuals of the local models were always noticeably smaller than those of the corresponding global models, which indicates that local polynomials have greater ability to model complex RPC errors, such as the distortions caused by orientation variations. When sufficient GCPs were available (the number of GCPs exceeded 10), the correction accuracy of the local quadratic model was significantly better than those of all other test models. Figure 5 further illustrates an example of the spatial distribution of correction residuals in the image space. The test scenario employed all 15 GCPs and the nadir-view image. Large errors in Figure 5a mainly arise at the edges of the scene, which indicates that some complex deviation existed in the results of the global affine model. The correction residuals of the global quadratic model (Figure 5b) were somewhat smaller than those of the global affine model, mainly observed at the GCPs (e.g., point No. 1, 3, and 6) but not at the ICPs (cf. Figure 4a). The local polynomial models removed almost all errors of the GCPs (Figure 5c-d), and their residuals at the ICPs were also noticeably smaller than those of the global-polynomial-based models. The local quadratic model performed the best among the four test methods. At most checkpoints, correction residuals of the local quadratic model were significantly smaller than those of all competitors.

Results
significantly better than those of all other test models. Figure 5 further illustrates an example of the spatial distribution of correction residuals in the image space. The test scenario employed all 15 GCPs and the nadir-view image. Large errors in Figure  5a mainly arise at the edges of the scene, which indicates that some complex deviation existed in the results of the global affine model. The correction residuals of the global quadratic model (Figure 5b) were somewhat smaller than those of the global affine model, mainly observed at the GCPs (e.g., point No. 1, 3, and 6) but not at the ICPs (cf. Figure 4a). The local polynomial models removed almost all errors of the GCPs (Figure 5c-d), and their residuals at the ICPs were also noticeably smaller than those of the global-polynomial-based models. The local quadratic model performed the best among the four test methods. At most checkpoints, correction residuals of the local quadratic model were significantly smaller than those of all competitors.  The performance of different bias-correction models was also compared in the object space. After correcting RPC errors in the image space by using different methods, the geographic coordinates of the object points were then calculated by forward intersection, and they were finally transformed to a Cartesian coordinate system. The Euclidean distances between the calculated object coordinates and their corresponding true values (i.e., GPS survey results) were computed, and the statistical results of the ICPs are shown in Figure 4d. It can be seen that the accuracies of the local polynomial models were always better than those of the global polynomial models when the same polynomial form was used. Correction residuals of the local affine model were overall 5%~10% less than those of the global affine model. The performance of the local quadratic model improved more rapidly than that of the other Remote Sens. 2017, 9, 200 9 of 12 three models with an increase in the number of redundant observations. When all 15 GCPs were used, the correction residuals of the local quadratic model were about 20% and 35% smaller than those of the local affine model and the global quadratic model, respectively. More experiments were conducted using different training sets to reduce the effect of the spatial distribution of GCPs on the performance assessment. Figure 6 shows the box plots of 100 tests with using 15 GCPs. In each test, the four points at the corners of the scene were fixed as GCPs, and the other 11 GCPs were randomly chosen from the remainder of ground points. It can be seen that the standard deviations of affine models were always smaller than those of the corresponding quadratic models, which shows that the performance of lower-order polynomial models is less influenced by the spatial distribution of GCPs. The local polynomial models achieved better accuracy than the global models for all three images and object coordinates. The local quadratic model performed the best in the four methods. Its mean residual in the object space was about 9% and 27% smaller than those of the local affine model and the global quadratic model, respectively. the object points were then calculated by forward intersection, and they were finally transformed to a Cartesian coordinate system. The Euclidean distances between the calculated object coordinates and their corresponding true values (i.e., GPS survey results) were computed, and the statistical results of the ICPs are shown in Figure 4d. It can be seen that the accuracies of the local polynomial models were always better than those of the global polynomial models when the same polynomial form was used. Correction residuals of the local affine model were overall 5%~10% less than those of the global affine model. The performance of the local quadratic model improved more rapidly than that of the other three models with an increase in the number of redundant observations. When all 15 GCPs were used, the correction residuals of the local quadratic model were about 20% and 35% smaller than those of the local affine model and the global quadratic model, respectively. More experiments were conducted using different training sets to reduce the effect of the spatial distribution of GCPs on the performance assessment. Figure 6 shows the box plots of 100 tests with using 15 GCPs. In each test, the four points at the corners of the scene were fixed as GCPs, and the other 11 GCPs were randomly chosen from the remainder of ground points. It can be seen that the standard deviations of affine models were always smaller than those of the corresponding quadratic models, which shows that the performance of lower-order polynomial models is less influenced by the spatial distribution of GCPs. The local polynomial models achieved better accuracy than the global models for all three images and object coordinates. The local quadratic model performed the best in the four methods. Its mean residual in the object space was about 9% and 27% smaller than those of the local affine model and the global quadratic model, respectively.

Discussions
The proposed method was also tested using some IKONOS, QuickBird, and other ZY-3 images. No obvious non-rigid RPC deformations were detected in these data sets, and the correction accuracies of the local and global polynomial models were at the same level. Figure 7 shows an IKONOS example. The study area is located in Hobart, Australia, and a total of 113 ground points were surveyed. As illustrated in Figure 7a, 77 points were reserved as ICPs, and part or all of the other 36 points were used as GCPs in the experiments. The results in Figure 7b show that, when more than 20 GCPs were employed, the image residuals of all four test models barely reduced with an increase in the number of redundant observations. The bias-correction accuracies of the local models were only slightly better than those of the global models, which implies that no serious non-rigid bias exists in the RPC data of the test IKONOS image. Given that strong non-rigid RPC errors were detected only in all three images of one ZY-3 stereo triplet, but not in other test images, it can be inferred that the errors are most likely caused by complex navigation errors induced by temporary short-period orbital perturbations.
The proposed method was also tested using some IKONOS, QuickBird, and other ZY-3 images. No obvious non-rigid RPC deformations were detected in these data sets, and the correction accuracies of the local and global polynomial models were at the same level. Figure 7 shows an IKONOS example. The study area is located in Hobart, Australia, and a total of 113 ground points were surveyed. As illustrated in Figure 7a, 77 points were reserved as ICPs, and part or all of the other 36 points were used as GCPs in the experiments. The results in Figure 7b show that, when more than 20 GCPs were employed, the image residuals of all four test models barely reduced with an increase in the number of redundant observations. The bias-correction accuracies of the local models were only slightly better than those of the global models, which implies that no serious non-rigid bias exists in the RPC data of the test IKONOS image. Given that strong non-rigid RPC errors were detected only in all three images of one ZY-3 stereo triplet, but not in other test images, it can be inferred that the errors are most likely caused by complex navigation errors induced by temporary short-period orbital perturbations.

Conclusions
In this paper, we have presented a local-polynomial-based method for bias compensation of rational polynomial coefficients, which can be used in the photogrammetric processing of high-resolution satellite imagery. Benefiting from the flexibility of local polynomials on data modeling, the new method has great potential for removing complex errors in vendor-provided RPCs, such as the distortions caused by perturbations in satellite motion. Our preliminary results on a Chinese ZY-3 satellite data set show that, when the same polynomial form was used, accuracy of the local polynomial model can be noticeably better than that of the corresponding global polynomial model. More specifically, correction residuals of the local affine model were about 15% smaller than those of the global affine model. The performance of the local quadratic model was significantly higher than that of the other test models when sufficient GCPs were available. In the test scenario using 15 GCPs, the accuracy of the local quadratic model was about 9% and 27% better than those of the local affine model and the global quadratic model, respectively. Our future research on this topic will focus on testing the performance of the local polynomial models on more data from different sensors.

Conclusions
In this paper, we have presented a local-polynomial-based method for bias compensation of rational polynomial coefficients, which can be used in the photogrammetric processing of high-resolution satellite imagery. Benefiting from the flexibility of local polynomials on data modeling, the new method has great potential for removing complex errors in vendor-provided RPCs, such as the distortions caused by perturbations in satellite motion. Our preliminary results on a Chinese ZY-3 satellite data set show that, when the same polynomial form was used, accuracy of the local polynomial model can be noticeably better than that of the corresponding global polynomial model. More specifically, correction residuals of the local affine model were about 15% smaller than those of the global affine model. The performance of the local quadratic model was significantly higher than that of the other test models when sufficient GCPs were available. In the test scenario using 15 GCPs, the accuracy of the local quadratic model was about 9% and 27% better than those of the local affine model and the global quadratic model, respectively. Our future research on this topic will focus on testing the performance of the local polynomial models on more data from different sensors.