Effective Determination of Optimal Regularization Parameter in Rational Polynomial Coefficients Derivation*

Journal of the Korean Society of Surveying, Geodesy, Photogrammetry and Cartography.
2013.
Nov,
31(6_2):
577-583

This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/3.0) which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.

- Received : November 25, 2013
- Accepted : December 31, 2013
- Published : November 28, 2013

Download

PDF

e-PUB

PubReader

PPT

Export by style

Share

Article

Metrics

Cited by

TagCloud

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.
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.
, where
r_{n}
and
c_{n}
are the normalized row and column indices of the pixels in the image space, and
X_{n} , Y_{n}
, and
Z_{n}
represent normalized coordinate values of the object points in the ground space. Here,
a_{ijk}, b_{ijk}, c_{ijk}
, and
d_{ijk}
are polynomial coefficients called RPC. The total number of
a_{ijk}, b_{ijk}, c_{ijk}
, and
d_{ijk}
is 80.
b_{1}
and
d_{1}
are usually set to 1. Therefore, the number of RPC is 78. The normalization of the coordinates is calculated as
, where
r_{0}
and
c_{0}
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,
X_{0}, Y_{0}
, and
Z_{0}
are offset values for the ground coordinates. In addition,
X_{s}, Y_{s}
, and
Z_{s}
are scale values for the ground coordinates. The offsets and scales normalize the coordinates to [-1. 1].
To obtain the
a_{ijk}, b_{ijk}, c_{ijk}
, and
d_{ijk}
, a least square solution is used. Eqs. (1) and (2) are rewritten as nonlinear condition equations,
, where
l
represents the observation vector (
r_{n}, c_{n}, X_{n} , Y_{n},
and
Z_{n}
) and
x
represents the parameter vector (
a_{ijk}, b_{ijk}, c_{ijk}
, and
d_{ijk}
). 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
x^{0}
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
f_{x}
is given by
, where
l_{x}
is a vector of parameter observations, and
W
is weight matrix.
x^{0}
is updated
x^{0}
plus Δ . The iteration step is stopped when Δ is less than a threshold, which is insignificantly small. The final least square estimates of parameter
is
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),
k_{opt 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,
N_{CK}
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
λ
,
r_{rigi}
and
c_{rigi}
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).
Fitting points, checking points, and exposure station in the ground space(3D)
Fitting points, checking points, and exposure station in the ground space (X-Y plane)
Fitting points, checking points, and exposure station in the ground space (X-Z plane)
Fitting points and checking points in the image space
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.
RP vs. RMSE, a range of RPs from 4x10^{-7} to 3.2x10^{-5}
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.
Residual norm vs. solution norm, a range of ks from 1x10^{-6} to 1x10^{-1}
Results of proposed method and L-curve method

Sensor model
;
Rational functional model
;
Rational polynomial coefficients
;
Regularization Parameter

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
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:
PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

3. REGULARIZATION

The RPC in the RFM are highly correlated between coefficents. As a result, the matrix
PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

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).
PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

Parameters used for rigorous sensor model

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

PPT Slide

Lager Image

Results of proposed method and L-curve method

PPT Slide

Lager Image

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.
Acknowledgements

This work was supported by a grant from a Strategic Research Project (Development of Site Investigation and Monitoring System for Extreme Cold Region) funded by the Korea Institute of Construction Technology.

Choi H. G.
,
Thite A. N.
,
Thompson D. J.
2007
Comparison of methods for parameter selection in Tikhonov regularization with application to inverse force determination
Journal of Sound and Vibration
304
(35)
894 -
917
** DOI : 10.1016/j.jsv.2007.03.040**

Di K.
,
Ma R.
,
Li R. X.
2003
Rational functions and potential for rigorous sensor model recovery
PE&RS
69
(1)
33 -
41

Dial G.
,
Grodecki J.
2005
RPC replacement camera models
ASPRS
Proceedings of the 2005 ASPRS Annual Conference
Baltimore, Maryland
7-11 March

Fraser C. S.
,
Dial G.
,
Grodecki J.
2006
Sensor orientation via RPCs
ISPRS Journal of Photogrammetry and Remote Sensing
60
(3)
182 -
194
** DOI : 10.1016/j.isprsjprs.2005.11.001**

Gerald C. F.
,
Wheatley P. O.
1994
Applied Numerical Analysis
Addison-Wesley Publishing Company
624 -

Grodecki J.
2001
IKONOS stereo feature extraction - RPC approach
ASPRS
Proceedings of the 2001 ASPRS Annual Conference
St. Louis, Missouri
23-27 April

Hansen P. C.
1992
Analysis of discrete ill-posed problem by means of L-curve
SIAM Review
34
(4)
561 -
580
** DOI : 10.1137/1034115**

Krawczy-Stando D.
,
Rudnicki M.
2007
Regularization parameter selection in discrete ill-posed problems - the use of the U-curve
International Journal of Applied Mathematics and Computer Science
17
(2)
157 -
164

Lin X.
,
Yuan X.
2008
Improvement of the stability solving rational polynomial coefficients
ISPRS
Beijing, China
International Archives of the Photogrammetry, Remote Sensing and Spatial Science
3-11 July
XXXVII
(B1)
711 -
716

Tao C. V.
,
Hu Y.
2001
A comprehensive study of the rational function model for photogrammetric processing
PE&RS
67
(12)
1347 -
1357

Tao C. V.
,
Hu Y.
2002
3D reconstruction methods based on the rational function model
PE&RS
68
(7)
705 -
714

Tong X.
,
Liu S.
,
Weng Q.
2010
Bias-corrected rational polynomial coefficients for high accuracy geopositioning of QuickBird stereo imagery
ISPRS Journal of Photogrammetry and Remote Sensing
65
(2)
218 -
226
** DOI : 10.1016/j.isprsjprs.2009.12.004**

Wolf P. R.
,
Dewitt B. A.
2000
Elements of Photogrammetry with Application in GIS
McGraw-Hill Companies
624 -

Zhan Y.
,
Liu C.
,
Qiao G.
2008
Accuracy evaluation of rational polynomial coefficients solution for QuickBird imagery based on auxiliary ground control points
International Archives of the Photogrammetry, Remote Sensing and Spatial Science
Beijing, China
3-11 July
XXXVII
(B7)
1287 -
1294

Citing 'Effective Determination of Optimal Regularization Parameter in Rational Polynomial Coefficients Derivation*
'

@article{ GCRHBD_2013_v31n6_2_577}
,title={Effective Determination of Optimal Regularization Parameter in Rational Polynomial Coefficients Derivation*}
,volume={6_2}
, url={http://dx.doi.org/10.7848/ksgpc.2013.31.6-2.577}, DOI={10.7848/ksgpc.2013.31.6-2.577}
, number= {6_2}
, journal={Journal of the Korean Society of Surveying, Geodesy, Photogrammetry and Cartography}
, publisher={Korean Society of Surveying, Geodesy, Photogrammetry and Cartography}
, author={Youn, Junhee
and
Hong, Changhee
and
Kim, TaeHoon
and
Kim, Gihong}
, year={2013}
, month={Nov}