Accuracy and Efficiency: The Comparison of Different RPC Parameters Solving Methods

Longhui Wang, Tao Wang, Yan Zhang, Kun Zhang

Information Engineering University, Zhengzhou, China.

DOI: 10.4236/gep.2020.810008                

Abstract

As a generalized sensor, the RPC model with its accuracy equally matches the physical sensor model. Moreover, the accurate positioning combining with the flexibility in application leads the RPC model to be the priority in photogrammetry processing. Generally, the RPC model is calculated through a control grid. Different RPC parameters solving methods and the operation efficiency all serve as variables in the accuracy of the model. In this paper, the ridge estimation iterative method, spectrum correction iteration, and conjugate gradient method are employed to solve RPC parameters; the accuracy and efficiency of three solving methods are analyzed and compared. The results show that ridge estimation iterative method and spectrum correction iteration have obvious advantages in accuracy. The ridge estimation iterative method has fewer iteration times and time con-sumption, and spectrum correction iteration has more stable precision.

Keywords

Rational Polynomial Coefficients (RPC), Ridge Estimation Iterative Method, Spectrum Correction Iteration, Conjugate Gradient Method

1. Introduction

Rational polynomial coefficients (RPC) model expresses the coordinates of image points via the ratio of ground coordinate polynomials a generalized sensor model and is a more generalized expression of the rigorous geometric model (Zhang & Liu, 2004; Grodecki & Dial, 2003). Tao and Hu (2001) introduce the solving principle and calculation process of the RPC model in his research and comes up with the least-squares of parameter solution. The positioning accuracy of the RPC model is evaluated by SPOT satellite image and aerial image, and the results show that its accuracy can approach the rigorous geometric model. Zhu, Li, and Li (2003) use the RPC model to carry out a geometric correction for the TH-1 panchromatic image. The experimental results show that the RPC model can achieve high correction accuracy for different landforms (flat, mountainous, and hilly areas). Yu and Tian (2017) managed to apply the RPC model to the geometric correction of the GF-1 panchromatic image with high accuracy. RPC model is suitable for all kinds of sensors and has a strong application prospect. The accuracy and efficiency of parameters solving methods become the focus of research.

In the process of RPC parameters solving methods, a large number of RPC parameters are related to each other. The coefficient matrix of the normal equation constructed is ill-conditioned and irreversible, which makes the least-squares adjustment unable to get a high-precision stable solution. Therefore, the following three frequently-used solving methods on the least-squares appear. Among them, ridge estimation (Yuan & Lin, 2008) modifies the state of the RPC model and ensures the stability of the solution by changing the equation equivalence relationship. For the experimental verification of SPOT image and Quick Bird image, ridge estimation based on L-curve method can effectively change the ill-condition of the normal matrix and the positioning accuracy of the RPC model. Spectrum correction iteration [Wang & Liu, 2002; Wang, Liu, & Huang, 2003] which adopts the idea of unbiased estimation estimates the parameters of the RPC model without changing the equation equivalence. Experiments indicate that the method can obtain stable results and is suitable for solving large-scale equations. The conjugate gradient method [Gu, Gui, Zhang, & Wei, 2014; Zhang, 2015] is commonly used in optimization and is suitable for solving large-scale equation optimization problems. It avoids the inversion of the normal matrix by solving iteratively. To overcome the ill condition of the normal matrix, the above algorithms intend to solve the ill condition of the model from different aspects, to obtain a stable solution. In this paper, ridge iterative estimation method, spectrum correction iteration and conjugate gradient method are studied in depth to form an objective and comprehensive comparative analysis, and the TH-1 three-line-array down view image and high-resolution segmented mosaic image are used to verify, and the accuracy and efficiency of different methods are comprehensively compared.

2. RPC Model

2.1. Definition

RPC model connects image coordinate ( $r,c$ ) with corresponding ground coordinate ( $P,L,H$ ) in the form of polynomial ratio [Liu, 2003; Meng, 2015]. To reduce the errors in the course of calculation and enhance the stability of the solution, the image coordinates and ground coordinates are normalized to the range of (−1, 1). RPC model is defined as follows:

$\begin{array}{l}r=N_S(P,L,H)/D_S(P,L,H)\\ c=N_L(P,L,H)/D_L(P,L,H)\end{array}$ (1)

where

$\begin{array}{c}{N}_S(P,L,H)=a_0+a_{1}L+a_{2}P+a_{3}H+a_{4}LP+a_{5}LH+a_{6}PH+a_7{L}^2\\ +a_8{P}^2+a_9{H}^2+a_{10}LPH+a_{11}{L}^3+a_{12}L{P}^2+a_{13}L{H}^3\\ +a_{14}{L}^{2}P+a_{15}{P}^3+a_{16}P{H}^2+a_{17}{L}^{2}H+a_{18}{P}^{2}H+a_{19}{H}^3\end{array}$ (2)

$\begin{array}{c}{D}_S(P,L,H)=b_0+b_{1}L+b_{2}P+b_{3}H+b_{4}LP+b_{5}LH+b_{6}PH+b_7{L}^2\\ +b_8{P}^2+b_9{H}^2+b_{10}LPH+b_{11}{L}^3+b_{12}L{P}^2+b_{13}L{H}^3\\ +b_{14}{L}^{2}P+b_{15}{P}^3+b_{16}P{H}^2+b_{17}{L}^{2}H+b_{18}{P}^{2}H+b_{19}{H}^3\end{array}$ (3)

$\begin{array}{c}{N}_L(P,L,H)=c_0+c_{1}L+c_{2}P+c_{3}H+c_{4}LP+c_{5}LH+c_{6}PH+c_7{L}^2\\ +c_8{P}^2+c_9{H}^2+c_{10}LPH+c_{11}{L}^3+c_{12}L{P}^2+c_{13}L{H}^3\\ +c_{14}{L}^{2}P+c_{15}{P}^3+c_{16}P{H}^2+c_{17}{L}^{2}H+c_{18}{P}^{2}H+c_{19}{H}^3\end{array}$ (4)

$\begin{array}{c}{D}_L(P,L,H)=d_0+d_{1}L+d_{2}P+d_{3}H+d_{4}LP+d_{5}LH+d_{6}PH+d_7{L}^2\\ +d_8{P}^2+d_9{H}^2+d_{10}LPH+d_{11}{L}^3+d_{12}L{P}^2+d_{13}L{H}^3\\ +d_{14}{L}^{2}P+d_{15}{P}^3+d_{16}P{H}^2+d_{17}{L}^{2}H+d_{18}{P}^{2}H+d_{19}{H}^3\end{array}$ (5)

where $a_i$, $b_i$, $c_i$ and $d_i$ ( $i=0,1,\cdots ,19$ ) are rational polynomial coefficients, $b_0$ and $d_0$ are usually equal to 1; $(r,c)$ are standardized image point coordinates, and $(P,L,H)$ are standardized geodetic coordinates: i.e.,

$\begin{array}{c}r=(r_u-r_0)/r_s,c=(c_u-c_0)/c_s\\ P=(P_u-P_0)/P_s,L=(L_u-L_0)/L_s,H=(H_u-H_0)/H_s\end{array}$ (6)

where $r_u$, $c_u$, $P_u$, $L_u$ and $H_u$ are original coordinates; $P_0$, $P_s$, $L_0$, $L_s$, $H_0$ and $H_s$ are standardized parameters of ground coordinate; $r_0$, $r_s$, $c_0$ and $c_s$ are standardized parameters of image coordinate.

2.2. RPC Model Parameter Solution

The steps to build a virtual control point grid are as follows:

1) According to the size of the image, establish the uniformly distributed image grid points.

2) Obtain the maximum and minimum elevation value of the image coverage area by using DEM data, and layer uniformly within the elevation range to get the elevation value of each layer.

3) According to the rigorous imaging geometry model, the coordinates of ground points corresponding to each image grid point on each elevation plane are solved (Zhang & Li, 2007).

After the coordinates are normalized by Equation (6), the RPC parameters can be solved. We can transform Equation (1) as

$\begin{array}{l}{F}_X=N_S(P,L,H)-r\times D_S(P,L,H)=0\\ F_Y=N_L(P,L,H)-c\times D_L(P,L,H)=0\end{array}$ (7)

Then the error equation is

$V=AX-L,P$. (8)

where

$A=\left[\begin{array}{l}\begin{array}{cccc}\frac{\partial F_X}{\partial a_i}& \frac{\partial F_X}{\partial b_j}& \frac{\partial F_X}{\partial c_i}& \frac{\partial F_X}{\partial d_j}\\ \frac{\partial F_Y}{\partial a_i}& \frac{\partial F_Y}{\partial b_j}& \frac{\partial F_Y}{\partial c_i}& \frac{\partial F_Y}{\partial d_j}\end{array}\end{array}\right],(i=0,\cdots ,19;j=1,\cdots ,19)$

By the least-squares algorithm, the normal equation is displayed in Equation (9):

$(A^{T}PA)X=A^{T}PL$ (9)

Simplify Equation (9) to:

$NX=Y$ (10)

where $N=A^{T}PA$, $Y=A^{T}PL$ ; LS estimation of unknown parameter X is calculated as:

$X_{LS}=N^{-1}Y$ (11)

When the least-squares algorithm is employed, the normal equation is ill-conditioned because of the huge condition number of the normal matrix. The inversion of the normal matrix is not stable, so it is difficult to get a better least squares estimation. The following describes three RPC solving methods used in this paper.

3. Three RPC Solving Methods

3.1. Ridge Estimation Iterative Method Based on the L-Curve Method

Ridge estimation is an improved biased estimation of the least-squares algorithm. The narrow ridge estimation of the parameter is:

$X_0={(N+kI)}^{-1}Y$. (12)

where k is the ridge parameter, usually positive decimal; I is the unit matrix; $X_0$ is the ridge estimator of the ridge parameter.

It can be proved that there is always a value that makes $MSE(X_{(k)})<MSE(X)$ (MSE is the mean square error). That is, for a linear model, there is always a positive decimal k, which makes $X_0$ better than the least square estimation. Therefore, the key problem of ridge estimation is to find the optimal value k. The L-Curve method is to utilize $\lg \parallel AX-L\parallel$ as abscissa $\lambda$ and $\lg \parallel X\parallel$ as ordinate $\varphi$ to draw a graph. Abscissa and ordinate are all functions of k, which fit a curve (Wang & Ou, 2004). i.e.,

$f(\lambda ,\varphi )=f(\lg \parallel AX-L\parallel ,\lg \parallel X\parallel )$ (13)

Equation (13), $\parallel ·\parallel$ is the Euclidean two norm, which is the same below; The optimal solution k is at the maximum curvature of the curve. The calculation formula is as follows

$k=\arg \max \frac{{\lambda }^{\text{'}}{\varphi }^{\text{'}\text{'}}-{\lambda }^{\text{'}\text{'}}{\varphi }^{\text{'}}}{{({({\lambda }^{\text{'}})}^2+{({\varphi }^{\text{'}})}^2)}^{3/2}}$ (14)

where ${\lambda }^{\text{'}}$, ${\lambda }^{\text{'}\text{'}}$, ${\varphi }^{\text{'}}$ and ${\varphi }^{\text{'}\text{'}}$ represent the first and second-order forms of the derivatives of $\lambda$ and $\varphi$ respectively.

The L-curve method is used to obtain the approximate optimal ridge parameter k, which emphasizes the balance between the data fitting part $\parallel AX-L{\parallel }^2$ and the part $\parallel X{\parallel }^2$.

Generally, the ridge estimation method cannot obtain the last levels square error through a single solution, so the iterative method is used to substitute the parameter estimation into the model. Then the ridge estimation method is used to solve the parameter (Wu & Lu, 2019). By substituting $X_0$ into Equation (8), the observation vector is updated too. By updating the observation vector, the iterative equation becomes

$X_{m+1}=(N+k_{m+1}I)N{X}_m$. (15)

where $X_{m+1}$ and $X_m$ are the iteration values of unknown parameters. The purpose of ridge estimation is to diminish the mean square error and make the coefficient test better, but its disadvantage is that the estimator is biased. Equation (12) adds kI to the normal matrix N and the introduction of kI changes the unbiasedness of the least-squares estimator into a biased estimator. The larger k is, the larger the deviation is. The mean square error can be effectively reduced and the stability of parameter estimation can be improved by an iterative solution.

3.2. Spectrum Correction Iteration

Now add the parameter solution to both sides of Equation (10) simultaneously,

$(N+I)X=Y+X$. (16)

Since the symmetric positive characterization of N, there is no matter in a good or ill condition, $Rank(N+I)=n$ is always true, i.e., $N+I$ is a full rank square matrix.

It is noted that both sides of Equation (16) have a parameter X, so it can be solved iteratively as follows:

$X_{(n)}={(N+I)}^{-1}(Y+X_{(n-1)}),(n=1,2,\cdots )$. (17)

where $X_{(n)}$ and $X_{(n-1)}$ are the iteration. Equation (17) is the spectrum correction iteration.

Adding X to both ends of Equation (10) at the same time does not change the equivalence relation of the equation, and weaken the ill-condition of the normal equation. It can get a good approximate solution and is suitable for solving large-scale equations.

3.3. Conjugate Gradient Method

When the linear equations are large-scale equations, the conjugate gradient method is a commonly used effective solution. The steps are as follows:

Step 1: Select the initial value $X_0$.

Step 2: Calculate residual $r_0=Y-N{X}_0$, so $p_0=r_0,{\alpha }_0=\frac{r_0^T{r}_0}{p_0^{T}N{p}_0}$. With the first iteration result $X_1=X_0+{\alpha }_0{p}_0$, compute $r_1=Y-N{X}_1$ for standby.

Step 3: For $k=1,2,\cdots ,n$, repeat the following iterative process: Calculate $p_k=r_k+{\beta }_{k-1}{p}_{k-1}$, where ${\beta }_{k-1}=-\frac{p_{k-1}^{T}N{r}_k}{p_{k-1}^{T}N{p}_{k-1}}$ ; Calculate $a=\frac{r_k^T{p}_k}{p_{k-1}^{T}N{p}_{k-1}}$ to get the iterative result $X_{k+1}=X_k+{\alpha }_k{p}_k$ ; Estimate $r_{k+1}=Y-N{X}_{k+1}$. If $\parallel r_{k+1}{\parallel }^2<\epsilon$ ( $\parallel \cdot {\parallel }^2$ is the second norm of the European formula), stop iteration and go to step 4. Otherwise, continue iteration.

Step 4: Put the last $X_{k+1}$ as the approximate solution ${\hat{X}}_{CG}=X_{k+1}$.

The conjugate gradient iterates along the gradient direction to find the local optimal solution and obtains the optimal solution of the equations after finite iterations [Li, 2018]. In the process of solving, the inverse of the matrix is avoided.

4. Experiment and Result Analysis

The experimental hardware environment is Dell G7, CPU main frequency 2.60 G, and 16 G memory. The software platform is Windows 10 and VS 2015.

4.1. Experimental Data

To compare the accuracy and efficiency of the three RPC parameters solving methods mentioned in this paper, two TH-1 remote sensing images are selected as experimental data. The first image is a three-line-array nadir image covering Henan Province obtained in 2014, with a size of 12,000 pixels × 12,000 pixels and a ground resolution of 5 m. 8 ground points are accurately measured as control points of image. The accuracy of these points is better than 2 cm, and their distribution is shown in Figure 1(a). The second image is the regional high-resolution segmented mosaic image of Henan Province taken in 2016. There are 8 segments in total, each of which is 35,000 pixels × 4096 pixels in size. The ground resolution is 2 m, including 23 ground control points with measurement accuracy better than 2 cm (Figure 1(b)).

(a) (b)

Figure 1. Distribution of ground control points of TH-1 image: (a) three-line-array nadir image; (b) high-resolution mosaic image.

4.2. Accuracy Verification

In the experiment, the third-order RPC model with different denominators is used, and the RPC parameters solving experiment is carried out under the conditions of establishing the three-line-array image and high-resolution rigorous image respectively. In this paper, three methods are used for RPC parameters solving to test their efficiency. For 8 segments of the high-resolution image, each segment uses a rigorous imaging model to establish virtual control points and then conduct RPC solving respectively to test the accuracy of the RPC model of each segment. The test method is to take the control point (Figure 1) as the standard and use the coefficient calculated by the RPC parameters solving methods to back projection the control point to the image plane coordinate. The accuracy of the calculated image point coordinate value and the measured value is analyzed, and the maximum value, minimum value, mean value, and root mean square error are counted as the standard to measure the RPC model. The generated virtual grid has a plane size of 50 × 50 and an elevation of 5 layers.

It can be seen from Table 1 and Figure 2 that: 1) The accuracy of the three methods in the column direction is higher than that in the row direction, which may be caused by the small errors such as platform flutter that affect the positioning quality in the row direction. 2) The positioning accuracy of the ridge estimation iteration method and the spectrum correction iteration is better, and the accuracy of spectrum correction iteration reaches 10⁻² pixel in the row direction and 10⁻³ pixel in the column direction. 3) The positioning accuracy of the control points is observed, and there are different degrees of fluctuation between ridge estimation iteration and conjugate gradient method, and the performance of spectrum correction iteration is relatively stable, which may be related to the unbiased estimation characteristics, which is closer to the true value after iteration.

The solution time and the number of iterations is shown in Table 2, and the

(a) (b) (c) (d)

Table 1. Accuracy of different RPC parameters solving methods (Units: pixels).

(a) (b) (c) (d)

Figure 2. Comparison of positioning accuracy of control points. (a) Three-line-array image row direction; (b) Three-line-array image column direction; (c) High-resolution image row direction; (d) High-resolution image row direction.

Table 2. Solution time and iteration numbers of different RPC sloving method.

solution time is the shortest when the number of iterations of the ridge estimation iteration method is less; the time of the spectrum correction iteration is longer than that of the ridge estimation iteration method. Compared with the other two methods, the conjugate gradient method has lower efficiency in solving time and iterations.

5. Summary and Prospect

In this paper, three parameters solving methods are used to calculate the parameters of the RPC model, which are ridge estimation iterative method, spectrum correction iteration, and conjugate gradient method. The verification of TH-1 satellite three-line-array nadir image and high-resolution mosaic image data shows that the ridge estimation iterative method and spectrum correction iteration have obvious advantages in the accuracy of the solution and can better improve the positioning accuracy of RPC model. The advantage of the ridge estimation iterative method is that the number of iterations is less and the time is shorter. The advantage of the spectrum correction iteration is an unbiased estimation and the accuracy is more stable. The next step is to study the error propagation law between parameters, further improve the accuracy of the model solution, and improve the efficiency of large-scale ill-conditioned equation solutions.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Grodecki, J., & Dial, G. (2003). Block Adjustment of High-Resolution Satellite Images Describedby Rational Polynomials. Photogrammetric Engineering & Remote Sensing, 69, 59-68. https://doi.org/10.14358/PERS.69.1.59
[2]	Gu, Y., Gui, Q., Zhang, X., & Wei, M. (2014). Iterative Solution of Regulation to Ill-Conditioned Problems in Geodesy and Geophysics. Acta Geodaetica et Cartographica Sinica, 43, 331-336.
[3]	Li, M. (2018). Research on Two Types of Conjugate Gradient Methods. Ph.D. Thesis, Xi’an: XIDIAN University.
[4]	Liu, J. (2003). Research on High Resolution Satellite CCD Stereo Image Positioning Technology. Ph.D. Thesis, Zhengzhou: Information Engineering University.
[5]	Meng, W. (2015). Geometric Model for Push-Broom Satellite Imagery of Butted TDI-CCDs. Ph.D. Thesis, Zhengzhou: Information Engineering University.
[6]	Tao, C. V., & Hu, Y. (2001). A Comprehensive Study of the Rational Function Model for Photogrammetric Processing. Photogrammetric Engineering & Remote Sensing, 67, 1347-1357.
[7]	Wang, X., & Liu, D. (2002). Iteration Method by Correcting Characteristic Value for Ill-Conditioned Equations in the Least Square Estimations. Journal of Hubei Institute for Nationalities (Natural Science Edition), 20, 1-4.
[8]	Wang, X., Liu, D., & Huang, H. (2003). The Co-Factor Matrix of the Iteration Method by Correcting Characteristic Value. Geomatics and Information Science of Wuhan University, 28, 429-431.
[9]	Wang, Z., & Ou, J. (2004). Determining the Ridge Parameter in a Ridge Estimation Using L-Curve Method. Geomatics and Information Science of Wuhan University, 29, 235-238.
[10]	Wu, G., & Lu, T. (2019). Ridge Estimation Iterative Method for Illness Data Processing. Journal of Geodesy and Geodynamics, 39, 178-183.
[11]	Yu, B., & Tian, S. (2017). Geometric Correction Research on Satellite GF-1 Data. Remote Sensing Technology and Application, 32, 133-139.
[12]	Yuan, X., & Lin, X. (2008). A Method for Solving Rational Polynomial Coefficients Based on Ridge Estimation. Geomatics and Information Science of Wuhan University, 33, 1130-1133.
[13]	Zhang, G., & Li, D. (2007). The Algorithm of Computation RPC Model’s Parameters for Satellite Imagery. Journal of Image and Graphics, 12, 2080-2088.
[14]	Zhang, L. (2015). Nonlinear Conjugate Gradient Methods for Optimization Problems. Ph.D. Thesis, Changsha: Hunan University.
[15]	Zhang, Y., & Liu, J. (2004). The Positioning Algorithm Based on RPC Models and Its Optimizing of Stereo Images from High Resolution Remote Sensing Satellites. Engineering of Surveying and Mapping, 13, 1-4.
[16]	Zhu, Q., Li, X., & Li, S. (2003). Accuracy Assessment of TH-1 Satellite Image Geometric Correction Using Rational Function Model and Polynomial Model. Journal of University of Science and Technology of China, 43, 110-114.