The Rational Polynomial Coefficients Modification Using Digital Elevation Models

The high-resolution satellite imageries (HRSI) are as primary dataset for different applications such as DEM generation, 3D city mapping, change detection, monitoring, and deformation detection. The geo-location information of HRSI are stored in metadata called Rational Polynomial Coefficients (RPCs). There are many methods to improve and modify the RPCs in order to have a precise mapping. In this paper, an automatic approach is presented for the RPC modification using global Digital Elevation Models. The main steps of this approach are: relative digital elevation model generation, shift parameters calculation, sparse point cloud generation and shift correction, and rational polynomial fitting. Using some ground control points, the accuracy of the proposed method is evaluated based on statistical descriptors in which the results show that the geo-location accuracy of HRSI can be improved without using Ground Control Points (GCPs). * Corresponding author


INTRODUCTION
Nowadays, high-resolution satellite imageries have been widely used in vast range of remote sensing and photogrammetric applications.There are many methods and models that are developed to extract the 3D geospatial information from HRSI such as Rational Function Models (RFM).The RFMs describe the relation between a 2D image and 3D coordinates of an object using rational polynomials which include 78 RPCs.There are basically two approaches to determine the RPCs including terrain-dependent and terrain-independent methods.In the former method, for the solution of the RPCs, a dense network of GCPs are required, whereas in the second method, known physical sensor models are used to generate the RPCs.Therefore, in the second method, the space resection problem is solved without using GCPs.However, the final 3D coordinates of the generated object points in the terrain-independent method, suffer from the shift and drift errors due to the inherent bias in the measured physical sensor parameters.This necessitates the use of a minimum number of GCPs to decrease the shift and drift errors and modify the RPCs.In this case, the number and accuracy of GCPs should be sufficient to have a precise mapping using RPCs.Since the gathering of GCPs is a costly, time consuming, and limited process in many application, significant efforts are dedicated for bias compensation of high-resolution satellite images using existing dataset such as reference maps, LiDAR points, ortho photos, and Digital Elevation Models (DEMs).The template-based edge matching technique can be used to determine the local shift between the image and large scale topographic map features and model the RPCs bias compensation parameters precisely.In this technique, the georeferencing accuracy is about one pixel (Oh, et. al., 2015).However, it seems that this method has some restriction to be employed in mountainous areas, if there is no sufficient edge information in matching window.Using 3D LiDAR point cloud, the accuracy of relative and absolute orientation of highresolution satellite images can be improved to sub pixel values through an automatic least squares 3D surface matching.For this, a 3D surface is generated using raw RPCs of stereo images, and after matching, precise ortho images can be generated using relative and absolute transformation (Teo, et. al., 2013).However, existing the high resolution LiDAR dataset is not always possible over the area of interest.Moreover, the accuracy and quality of matched points may affect the final results.In order to affine-based RPC correction, an image-toimage matching method can be used to reduce the geometric difference between raw image and elevation dataset.Based on the manually measurements, the horizontal accuracy can be improved to only a few meters (Oh, et. al., 2013).The automatic DEM matching technique can be applied between a global DEM and generated DEM using stereo satellite images to solve the absolute orientation of images and bias error without ground control points with the accuracy of about 2*GSD or better (Jeong, et. al., 2012).In another research, affine RPCs correction is performed in two steps based on the global DSM alignment and georefernced low resolution images.The results show a significant improvement on planimetry and vertical accuracy of 6.7 m and 5.1 m, respectively (d'Angelo, et. al., 2010).

PROPOSED METHOD
In this paper, a different approach is adopted for the RPC modification by incorporating the global Digital Elevation Models (DEMs) into the intersection solution without using any GCPs based on the considered idea by (Jeong, et. al., 2012).The main steps of this approach are as follows (cf. Figure 1  The summary of each steps and their main components are given in the next sub-sections.

Relative digital elevation model generation
In the first step, a so called relative digital elevation model (r-DEM) is generated using the raw RPCs and with the same spatial resolution of stereo images.For this, some corresponding points are extracted automatically in image space using Normalized Cross Correlation (NCC) or Fast Fourier Transformation (FFT) matching techniques.The accuracy of matched points should be in sub-pixel and existing mismatches should be removed by user or other efficient methods such as RANSAC.Next, using RPCs, the space intersection is solved for each corresponding point in image space and the 3D coordinates of points are calculated regarding to object space (WGS84 coordinate system) as Equation 1. (1) Where, x, y are the known image coordinates and {a0… d19} are raw RPCs.The X, Y, Z are unknown corresponding points in object space.Four equations are necessary to calculate the 3D coordinates of each matched image points.To r-DEM generation, the epipolar resampling technique is used to have a robust and simple dense image matching in 1D epipolar lines space.So, the shift errors in RPCs are transferred to r-DEM coordinates.

Shift parameters calculation
Before computing the shift parameters, the orthometric elevations of global DEMs with the EGM96 datum should be converted to ellipsoidal elevations regarding to World Geodetic System (WGS84) (Tachikawa, et. al., 2011b).To calculate the shift parameters between 3D coordinates of r-DEM and their corresponding absolute DEM (a-DEM), some matched points should be measured with sub-pixel accuracy.Because of the low information content of global DEMs, an efficient DEM matching approach is proposed here.First, the Scale Invariant Feature Transform (SIFT) (Lowe, 2004) is employed to detect corner points in each dataset.Then, some candidate of matched points between two 3D dataset are found based on applying a predefined threshold on sum of squared distance (SSD) of SIFT feature vectors.Finally, the correct matched points are selected as inliers by RANSAC.These correspondences are utilized to calculate the 3D coordinate differences.

Sparse point cloud generation and shift correction
In this step, a dense network of corresponding points are extracted from stereo images using the Fast Fourier Transformation (FFT) image matching method and with subpixel accuracy.Then, the 3D coordinates of these points is calculated using raw RPCs and the space intersection equations.These points are used to reproduce the modified RPCs.As shown in flowchart of Figure 1, the 3D coordinates of the matched points are corrected based on the calculated shift parameters of the second step.

Rational polynomial fitting
Finally, the rational polynomials are fitted to the image coordinates and their corresponding corrected object coordinates of the dense network points and new RPCs is calculated for each image in a bundle adjustment solution based on WGS84 coordinate system, as in Equation 1.

Case study
In this paper, the raw RPCs of IRS-P5 stereo images are used to assess the adopted strategy over the study area that is shown in Figure 2. The selected region is a mountainous area with ca.1500 meters elevation range and 20x24 km 2 area in west of Iran.Two different global DEM such as SRTM and ASTER GDEM are used to assess the "a-DEM's" spatial resolution effect on the accuracy of the modified RPCs.The ground sampling distance of SRTM is 90 meters (~ 3 arc-seconds) and the overall absolute vertical height accuracy is about 10 meters (Rodriguez et al. 2006).The ASTER GDEM ver2 was released by NASA and METI on 2011 to the public at a spatial resolution of 30 meters.The root mean square error of ASTER GDEM ver2 is 8.68 meters when compared against geodetic control points (Gesch, 2012).

Implementation
As shown in Figure 3, the r-DEM is generated using raw stereo images and about 50 matched points with a matching accuracy of 0.3 pixel.

Figure 3. The relative generated DEM using raw RPCs
The SIFT plus RANSAC method is used to calculate the shift parameters as shown in Figure 4. Table 2. Shift values between r-DEM and global DEMs The X and Y coordinate differences between SRTM and r-DEM is shown in Figure 5.As mentioned before, a dense network of corresponding points is generated in object space with a matching accuracy of 0.8 pixel.Then, 3D coordinates of this network are corrected regarding to the sift parameters.Consequently, the new RPCs can be calculated using image and object, corrected coordinates of dense network points.

Accuracy Assessments
To assess the accuracy of the modified RPCs, about 9 GCPs which are measured by GPS, are used to compare the calculated values of the 3D coordinates before and after RPCs modification with the measured coordinates of the corresponding GCPs.Based on the final results in Table 3 , the RMSE of the X coordinate differences between calculated values using raw RPCs and the original values of GCPs is about 169 m while the RMSE of the differences between the calculated X coordinates using modified RPCs and their corresponding original values of GCPs is about 15 m in both SRTM and ASTER GDEM cases.However, the mean error of the X coordinates is increased after RPCs modification.The RMSE and mean values of the Y coordinate differences are decreased from 36 m to 9 m and from 228 m to 8 m, respectively.Finally, the mean errors of the differences between the calculated Z coordinates using modified RPCs and their corresponding original values of GCPs is about -69 m while the mean errors are decreased after RPCs modification.There is no significant improvement in RMSE values of Z coordinates.

Table3. The results of RPC modification
As another result of assessment values in Table 2, the ASTER GDEM can improve the RPCs a little better than SRTM.In other assessment method, the RMSE of bundle adjustment of 9 GCPs using the raw and modified RPCs are compared.Using raw RPCs, the RMSE is about 9 m and 1.15 m for X and Y coordinates respectively, while these values are about 0.09 m and 0.31 m after RPCs modification using SRTM and ASTER GDEM, separately.

CONCLUSION
In order to precise mapping using the high-resolution satellite imageries, it is necessary to modify the RPCs.In this paper, an automatic and efficient approach is investigated based on global digital elevation models to modify the raw RPCs of IRS-P5 stereo images.The main steps of this approach are: relative digital elevation model generation, shift parameters calculation, sparse point cloud generation and shift correction, and rational polynomial fitting.As the final conclusion, the adopted method, as far as the planimetric coordinates are concerned, can be used to modify the RPCs of HRSI and improve the coordinates without using GCPs.Moreover, the test results also indicate that the role of spatial resolution of global DEM is not very significant on the RPC modification accuracy for generating X and Y components.For the data sets used in this investigation, no improvement is observed in the Z components of the generated 3D coordinates with the modified RPCs.This requires further assessments and scrutiny to find out the factors that inhibited the Z-component improvement.
): I. Relative digital elevation model generation, II.Shift parameters calculation, III.Sparse point cloud generation and shift correction, IV.Rational polynomial fitting.

Figure 1 .
Figure 1.The main steps of proposed method for RPCs modification

Figure 2 .
Figure 2. Study area in west of Iran

Figure 4 .
Figure 4. SIFT + RANSAC matching method The average values of the calculated shift parameters for SRTM and ASTER GDEM are presented inTable 2, separately.

Figure 5 .
Figure 5. Shift between SRTM and r-DEM