Effective Determination of Optimal Regularization Parameter in Rational Polynomial Coefficients Derivation

Abstract

Recently, massive archives of ground information imagery from new sensors have become available. To establish a functional relationship between the image and the ground space, sensor models are required. The rational functional model (RFM), which is used as an alternative to the rigorous sensor model, is an attractive option owing to its generality and simplicity. To determine the rational polynomial coefficients (RPC) in RFM, however, we encounter the problem of obtaining a stable solution. The design matrix for solutions is usually ill-conditioned in the experiments. To solve this unstable solution problem, regularization techniques are generally used. In this paper, we describe the effective determination of the optimal regularization parameter in the regularization technique during RPC derivation. A brief mathematical background of RFM is presented, followed by numerical approaches for effective determination of the optimal regularization parameter using the Euler Method. Experiments are performed assuming that a tilted aerial image is taken with a known rigorous sensor. To show the effectiveness, calculation time and RMSE between L-curve method and proposed method is compared.

1. INTRODUCTION

Increasing demand for accurate spatial information, massive archives of ground information are provided by imagery from new sensors (e.g. IKONOS, GeoEye, and QuickBird satellite etc.). To establish a functional relationship between the image space and the ground space, sensor models are required. There are two categories for sensor models. One is a rigorous sensor model and the other is a generalized sensor model.

A rigorous sensor model produces high accuracy; however, it has several disadvantages. A rigorous sensor model, which is based on collinearity equations, presents the rigorous imaging geometric relationship between an image point and the homologous ground point, with parameters of physical meaning (Tong et al., 2010). The advantages of the rigorous sensor model are that it is very suitable for adjustment by analytical triangulation and normally yields high accuracy. However, it is complicated and the parameters used in the rigorous sensor model are usually kept confidential in the sensors. Some commercial satellite vendors have adopted a generalized sensor model, instead of the rigorous sensor model, to the end users (Di et al., 2003; Tao and Hu, 2002).

The rational functional model (RFM), which is one of the most popular generalized sensor models (Tong et al., 2010; Tao and Hu, 2001), has drawn special attention in the photogrammetry and remote sensing fields. One of the advantages of RFM is its sensor independence. The transformation between the image point and ground point is represented as a rational function without modeling the rigorous imaging process. Although rigorous sensor models produce more accurate results, the difference is negligible (Dial and Grodecki, 2005). Grodecki(2001) reported that the IKONOS rational model differs by no more than 0.04 pixel from the rigorous sensor model, with the RMS error below 0.01 pixel. Furthermore, the sensor information is effectively kept confidential by providing the RFM, due to the difficulty of transforming from the RFM to the rigorous sensor model.

The RFM expresses image coordinates as a ratio of two polynomials with variables of ground coordinate. The polynomial coefficients in the RFM are called rational polynomial coefficients (RPC). Comparing the rigorous sensor model and the RFM, the RFM would be overparameterized (Fraser et al., 2006). The ground to image transformation can be achieved with an 8-coefficient affine model. On the other hand, the usual number of RPC is 80. That could cause the design matrix to become almost rank deficient due to the complex correlation among RPC (Lin and Yuan, 2008) resulting in numerical instability in the least square adjustment (Lin and Yuan, 2008). To solve such an instability problem, a regularization technique with a regularization parameter (RP) is generally used.

The determination of an optimal RP has still been a challenging topic. There are several methods for the determination of an RP, including L-curve, U-curve, and ridge tracing methods. The L-curve (Hansen, 1992) is a loglog plot of the norm of a regularized solution versus the norm of the corresponding residual (fitting error) as the RP is varied (Choi et al., 2007). In Hansen(1992), the point on the L-curve that had the maximum curvature should be chosen as the corner of the curve and consequently became the optimal RP. This method can be automated without plotting the curve. An improvement to this approach using “U-curve” is reported in Krawczy-Stando and Rudnicki(2007). In this approach, the authors proposed a function that is expressed as a summation of the inverse number of the solution norm and the inverse number of the residual norm. They chose the optimal RP as the one producing the minimum value in the function. The limitation of L-curve method is that finding optimal RP is dependent on the number of k set in equation 19. With the small number of k set, finding optimal RP could be failure, because the curvature can be distorted. While the large number of k set, calculation time would be increased. The calculation time issue is similarly applied to U-curve method. In the ridge tracing method, root mean square errors are computed for a large number of different RPs. The best one is selected by suitable heuristics (e.q. trial method). Applying the ridge tracing methods to obtain an RP in RPC derivation is tried by Tao and Hu(2001) and Zhan et al.(2008). They chose the optimal RP within the range extracted by experiments, and proved that RMSE is not sensitive to a particular RP as long as the RP is within the specific range. The advantage of the ridge tracing method is its simplicity and generality. That means it is independent of condition number and noise. In this paper, we introduce a numerical approach for effective determination of an RP in the ridge tracing method using the Euler root-finding method. Our approach is independent of illcondition or well-condition, unlike the L-curve and U-curve methods. The RP can be determined at the target accuracy and can be calculated effectively. In this paper, we use the terminology “effective” as simple, general, and fast.

This paper is organized as follows. Section 2 presents a mathematical background of RPC derivation in RFM with a terrain independent model. This is followed by a numerical approach for the automatic detection of the RP in the regularization process. Our algorithm is exercised with terrain independent scenarios.

 

2. GENERAL SOLUTION FOR RFM

The RFM relates ground coordinates (X, Y, Z) to image coordinates (r, c) in the form of rational functions that are ratios of polynomials. For the ground to image transformation, the defined ratios of polynomials have the following form for each section:

, where rn and cn are the normalized row and column indices of the pixels in the image space, and Xn , Yn, and Zn represent normalized coordinate values of the object points in the ground space. Here, aijk, bijk, cijk, and dijk are polynomial coefficients called RPC. The total number of aijk, bijk, cijk, and dijk is 80. b1 and d1 are usually set to 1. Therefore, the number of RPC is 78. The normalization of the coordinates is calculated as

, where r0 and c0 are row offset values and column offset values for the image coordinates. rs and cs are row scale values and column scale values for the image coordinates. Similarly, X0, Y0, and Z0 are offset values for the ground coordinates. In addition, Xs, Ys, and Zs are scale values for the ground coordinates. The offsets and scales normalize the coordinates to [-1. 1].

To obtain the aijk, bijk, cijk, and dijk, a least square solution is used. Eqs. (1) and (2) are rewritten as nonlinear condition equations,

, where l represents the observation vector (rn, cn, Xn , Yn, and Zn) and x represents the parameter vector (aijk, bijk, cijk, and dijk). For the least square adjustment, the linearization of Eq. (8) follows:

, in which f is given by:

,where l is the approximation vector for the observation and x0 is the approximation vector for the parameters. Δ is the vector of corrections to the approximations for the parameters. v is the vector of observational residuals and B is the design matrix of the partial derivation of F with respect to the parameters as follows:

The normal equation matrix and the vector of correction matrix are represented respectively as:

, in which fx is given by

, where lx is a vector of parameter observations, and W is weight matrix. x0 is updated x0 plus Δ . The iteration step is stopped when Δ is less than a threshold, which is insignificantly small. The final least square estimates of parameter is

 

3. REGULARIZATION

The RPC in the RFM are highly correlated between coefficents. As a result, the matrix B in Eq. (11) is usually illconditioned and matrix N in Eq. (12) can become singular. It happens often when high order (i.e., more than second-order) polynomials in the RFM are used (Tao and Hu, 2001). The negative impact of this is that the iterative solution cannot be converged (Tao and Hu, 2001).

In order to tackle the possible ill-conditioned problem during the least square adjustment, a regularization technique, in which a small multiplication (i.e., RP) of the identity matrix is added, is applied. With this modification, the normal equation matrix in Eq. (12), correction vector in Eq. (13), and final estimates of the parameter matrix in Eq. (15) are rewritten as:

, where λ is RP.

Determination of RP λ is non-trivial. The L-curve method is finding out the point with biggest curvature for L-curve which is a log-log plot of the norm of a regularized solution versus the norm of the corresponding residual as the RPs are varied (Lin, 2008). L-curve is presented as

L-curve is L-shaped, therefore it presents vertical for small k, and approximately horizontal for large k, with the corner providing the optimal RP (Lin, 2008). To find out the biggest curvature:

,where ξ' and ξ'' are the first and the second derivatives of ξ on k respectively, and η' and η'' are the first and second derivatives of η on k respectively. In Eq. (20), kopt R is the determined optimal regularized parameter.

Here, we propose the effective determination of the optimal RP using the Euler Method. First, a root mean square error function is presented as:

Here, NCK is the number of checkpoints. rλi and cλi are image coordinates obtained by using Eqs. (1) and (2) with the known ground coordinates and adjusted parameters that are from Eq. (18) at specific λ , rrigi and crigi are image coordinates obtained by using the rigorous sensor model with the known ground coordinates.

From the Euler method (Gerald and Wheatley, 1994), the first derivative of the root mean square error function in Eq. (21) can be numerically written as:

The Euler step for finding the optimum λ is derived as:

where h is the step size. This equation is iteratively solved and the iteration is stopped when the value of the first derivative of the λ function is insignificantly small. Theoretically, a small Δλ and h make a more accurate solution when the iteration time is increased(Gerald and Wheatley, 1994).

 

4. EXPERIMENT AND RESULTS

The RPC in the RFM can be solved for with or without knowing the rigorous sensor model (Tao and Hu, 2001). The rigorous sensor model being available, the terrain independent solution is able to be developed. In the terrain independent scenarios, the RFM can be solved using a ground grid with its grid-point coordinates determined using a rigorous sensor model. We modify the terrain independent scenarios proposed by Tao and Hu(2001).

First, a 3D object grid in the ground space is established. The maximum and minimum easting-ground UTM coordinates of the grid are 506,600 and 507,460 meters, respectively. The maximum and minimum northing-ground UTM coordinates of the grid are 4,474,600 and 4,475,630 meters, respectively. The relief range is from zero to 400 meters. The grid intervals of easting, northing, and height are 100, 100, and 50 meters, respectively. Next, an image grid in the image space is determined. Using the available rigorous sensor model, corresponding image coordinates of the ground coordinates are calculated. Fig. 1 shows the used fitting points and checking points of the object grid. The location of the exposure station in the ground space is also presented in Figs 1, 2, and 3. In Fig. 4, we present fitting points and checking points in the image space. The number of fitting points and checking points are 71 and 284, respectively. The parameters used for the rigorous sensor model are shown in Table 1. For the details of the notation used in Table 1, see Wolf and Dewitt(2000).

Fig. 1.Fitting points, checking points, and exposure station in the ground space(3D)

Fig. 2.Fitting points, checking points, and exposure station in the ground space (X-Y plane)

Fig. 3.Fitting points, checking points, and exposure station in the ground space (X-Z plane)

Fig. 4.Fitting points and checking points in the image space

Table 1.Parameters used for rigorous sensor model

The unknown RPC can be calculated using the corresponding image and object grid points with Eqs. (1-18) presented in the previous section. To show the RP effects on the RMSE, we calculate the RMSE of checkpoints with increasing RPs. Figures 5 and 6 show the results of RP versus RMSE. In Figures 5 and 6, the x-axis denotes the RP and the y-axis denotes the RMSE calculated by Eq. (21). Fig. 5 shows the RMSE of the checkpoints with an RP range from 4x10-7 to 3.2x10-5. The step size of the RP is 1x10-7. The RMSE of the check points with an RP range from 4x10-5 to 1x10-1 are presented in Fig. 6. Figs. 5 and 6 show that the RMSE curve decreases when the RP increases from near zero RP, and continues nearly flat with increasing RPs. The RMSE curve is increased from the specific RP. Therefore, we conclude that the optimal RP can be determined using the Euler root finding method.

Fig. 5.RP vs. RMSE, a range of RPs from 4x10-7 to 3.2x10-5

Fig. 6.RP vs. RMSE, a range of RPs from 4x10-5 to 1x10-1

The optimal RP is calculated using Eqs. (21), (22), and (23). The initial RP is set to 1x10-1, and the step size h is set to 1x10-3 in Eq. (23). Δλ is 1x10-5 in Eq. (22). In this experiment the determined RP with the proposed method is 1.368x10-4 , and the RMSE for the RP is 0.1014 pixel. The calculation time for the proposed method is 0.0862 second. To show the effectiveness of the proposed method, optimal RP is obtained by using L-curve method. L-curve method is finding maximum curvature in the plot of the norm of a regularized solution versus the norm of the corresponding residual as the RP is varied. Candidates of optimal RP, k set in Eq. (19), are varied from 1x10-6 to 1x10-1. The interval of k set is 5x10-5. With the k set, solutions and residuals for each k is calculated by using Eqs. (1-18). From the Eq. (19), we obtain the L-curve presented in Fig. 7. Determined optimal RP being calculated by using Eq. (20) is 2.510x10-4, and the RMSE for the RP is 0.0949 pixel. The calculation time for the L-curve method is 1.0631 second. All algorithms are programmed by MATLAB software. The results of proposed method and L-curve method are shown in Table 2. In Table 2, RMSE differences for two methods are not a crucial, because those are less than 1 pixel. We are focusing the calculation times. In this experiment, the number of k set for L-curve is 2000, and the calculation time is 1.0631 second. The calculation time for the proposed method is only 0.0862 second. If we want to calculate the optimal RP with the L-curve method as fast as with the proposed method, the number of k set should be decreased about 170. With the small size of k set, calculated RP cannot be secured as optimal one, because points are not densely distributed in the curve. Therefore, we can conclude that the propose method is more time efficient than L-curve method.

Fig. 7.Residual norm vs. solution norm, a range of ks from 1x10-6 to 1x10-1

Table 2.Results of proposed method and L-curve method

 

5. CONCLUSIONS

This paper presents effective determination of optimal RP in rational polynomial coefficients derivation. The curve of RMSE versus RP has a tendency with the variation of RP. The curve decrease when RP increase from near zero RP, continues nearly flat with increasing RPs, and is increased from the specific RP. Therefore, optimal RP can be determined by using the Euler root finding method. First, we define the RMSE function for RP. Next, optimal RP, which produces minimum RMSE, is numerically calculated by Euler method. The calculation time for the proposed method is 0.0862 second, and RMSE is 0.1014 pixel in the experiment. Experiment shows that the proposed method is time efficient comparing with L-curve method. Our future plan include the how calculation time and RMSE have a tendency upon various rigorous sensor models and initial value for RP and step size.