Designing Software Prototype for Digital Surface Model Extraction from Stereo Satellite Imagery Based on Rational Function Model

There are many commercial software that can be used to perform Digital Surface Model (DSM) extraction from stereo satellite imagery based on Rational Function Model (RFM), however the license price is quite expensive. Therefore, this research aims to build the new software prototype for DSM extraction from stereo satellite imagery based on RFM, by using free development tool and integrating various components or libraries that are noncommercial (free and/or open source). The method applied in the software prototype consists of 3 main steps, i.e. conjugate points collection manually and automatically using image matching with RFM based Cross Correlation method, Forward RFM based Virtual Intersection to perform 3D reconstruction of conjugate points, and interpolate data conjugate points with Ordinary Kriging method to establish DSM in the raster format. Stereo satellite imagery of WorldView-2 which is equipped with Rational Polynomial Coefficients (RPCs) from the vendor was applied to verify the DSM accuracy that obtained from the software prototype. Based on RFM formed from the RPCs can be obtained DSM with accuracy of 2.47639575534581 meters in planimetry (CE90) and 4.21487430387638 meters in


Introduction
Along development of science and technology, since the late 1950s application and utilization of the Digital Elevation Model (DEM) has continued to expand in various disciplines such as mapping, remote sensing, civil engineering, mining engineering, geology, geomorphology, military engineering, land planning and communications, and supported by the rapid development of GIS in 1990s which in this case also caused the application of DEM is increasingly widespread [1]. As DEM's increasingly widespread and important applications grow, DEM's data acquisition and processing methods are also being developed. Nowadays, there are a great number of satellite systems designed to capture stereo imagery which capable to establish the DEM, particularly in the form of Digital Surface Model (DSM). In this case, DigitalGlobe™ is among the vendors which currently operates a number of satellites in their constellation with the ability of capturing stereo satellite imagery with very high resolution namely Very High Resolution Satellite Stereoscopic Imagery (VHRS Satellite Imagery), such as QuickBird-2, GeoEye-1, WordView-1, WordView-2 and WordView-3.
For the 3D object reconstruction including DSM extraction using stereo satellite imagery, sensor orientation modeling is fundamental and absolutely necessary [2]. Sensor orientation modeling can be classified into two major groups, namely Rigorous Sensor Models (RSM) and Generalized Sensor Models (GSM) [3,4]. RSM is a parametric model which requires accurate physical sensor parameters, both for the Interior Orientation Parameters (IOP) of the scanner and the Exterior Orientation Parameters (EOP) of each scan line of the imagery's scene, so that the real imaging condition can be reconstructed accurately using the data. Hence, RSM so far has been recognized as the most precise sensor model as it can fully describe the actual imaging geometry between the image space and the object space [2].
These physical sensor parameters necessary for the RSM-based modeling usually can be obtained from the ephemeris, attitude and velocity data included in the imagery distribution. Nevertheless, because of a variety of causes and reasons, both relating to security and commercial purposes, most vendors currently do not include any longer the physical sensor parameters of an imagery they distribute directly. On the other hand, RSM forming in the linear array pushbroom scanner system based on the Ground Control Points (GCPs) in this case is also not easy to perform by users because it requires a great number of GCPs. Consequently, nowadays, most imagery vendors offer an alternative, which is including the Rational Function Coefficients (RFCs) in the imagery they distribute, so that their users can use them as a basis for establishing the Rational Function Model (RFM) to replace RSM of those imagery [5]. Basically, RFM is a type of GSM which in this case can be used as a substitute for RSM and the results generated are also fairly good. A great number of studies have been conducted in an attempt to develop this GSM, particularly RFM, until the resulting accuracy of the RFM modeling proves comparable to that of the RSM modeling [2,5,6,7,8,9]. Besides, RFM is also a general model deemed fairly simple and practical to apply in various types of sensors. Considering those reasons, now there has been lots of commercial digital photogrammetry or image processing software to facilitate the DSM extraction from stereo satellite imagery based on RFM model, such as PCI Geomatics, LPS, ENVI etc. However, the license of these commercial softwares are quite expensive that makes users encounter problems in accessing such commercial software. For those reasons, this research aims to develop a software prototype for DSM extraction from stereo satellite imagery based on RFM model, by using development tool, components and libraries which are non-commercial (free and/or open source).

Materials and Methods
The main purpose of this research aimed to build the new software prototype for DSM extraction from stereo satellite imagery based on RFM, by using free development tool and integrating various components or libraries that are non-commercial. In practice, this research is divided into 3 main steps, namely, development of software prototype, data processing, and accuracy test.

Data Sources
To test DSM extraction from the software prototype, in this research has been obtained stereopair of WorldView-2 from the DigitalGlobe Foundation with coverage area in Sadeng Beach, Special Region of Yogyakarta, Indonesia, which is complete with Rational Polynomial Coefficients (RPCs) and other Image Support Data files. From GPS surveying which was done on May 6, 2016, using RTK NTRIP (Real Time Kinematic Network Transport RTCM via Internet Protocol), also obtained 11 Independent Check Points (ICPs) for the accuracy test. The sampling technique chosen to determine the ICPs is selective sampling, in which the distribution (number and density) of ICPs will be determined by terrain variation. Locations with more complex terrain variations will be given ICPs with higher amounts and densities according to their complexity, as generally seen in Figure 1.

Data Processing
Data Processing consists of several processes performed in accordance with the concept applied in the prototype of the software which in this case consists of 3 main steps, namely.

Conjugate Points Collection
Conjugate points are corresponding points in the overlapping areas of stereo imagery in which case 3D reconstruction can be applied to it to obtain 3D ground coordinates respectively. In this case, a facility is provided to automatically collect conjugate points based on image matching and manually collect in order to select conjugate points at Very Important Points (VIPs) which necessary to represent the surface properly. The image matching technique applied in this software prototype is the cross correlation method based on RFM which in this case also adopt the implementation of feature-based keypoint and image pyramid. Feature-based keypoint is used as target point to be searched for on the paired imagery. RFM is used as the basis for finding matched point candidate as the center of search window formation on the paired imagery. Image pyramid is used to expand the search area as well as speed up the search process. The calculation of the RFM based Cross Correlation applied in the software prototype can be found detailed in the [10].   .

3D Reconstruction
Forward RFM based Virtual Intersection can be done by the least squares adjustment iteratively so that 3D coordinates of each conjugate point can be obtained. The least squares adjustment calculations for completing the RFM based 3D reconstruction can be found detailed in [10] and [11].

DSM Interpolation
In order to obtain a continuous model in raster format, any conjugate points that have such 3D ground coordinates must be interpolated. The interpolation method chosen in this case is the ordinary kriging method that can consider the existing spatial autocorrelation. The calculation of the ordinary kriging applied in the software prototype can be found detailed in [10].

Accuracy Test
To access the quality that can be achieved on the DSM generated through the prototype of the software will be done accuracy assessment in Accuracy Test, both absolute accuracy and relative accuracy. Absolute accuracy is obtained from accuracy test based on Independent Check Points (ICPs) that have been obtained from GPS surveying using RTK NTRIP (Real Time Kinematic Network Transport RTCM via Internet Protocol) method as mentioned before. Relative accuracy is obtained based on visual observations of the quality of DSM produced in representing the terrain surface it represents. Visual observations on DSM is aimed at the aspect of how well and realistically the DSM is in presenting terrain surface shapes.

RFM based DSM Extractor
The main result of this research is a software prototype which capable of performing DSM extraction from stereo satellite imagery (stereopair) based on Rational Function Model (RFM). The methods applied in the software prototype can be broadly divided into 3 main steps, ranging from conjugate points collection manually and automatically using image matching with RFM based Cross Correlation method, Forward RFM based Virtual Intersection to perform 3D reconstruction of conjugate points, and interpolate data conjugate points with Ordinary Kriging method to establish DSM in the raster format. By utilizing development tool and various components along with libraries that are free as mentioned in the research method, able to generate software prototype with Graphical User Interface (GUI) is quite interesting and user-friendly, hereinafter called "RFM based DSM Extractor".
Generally, the interface of "RFM based DSM Extractor" consists of a ribbon menu and two main workspaces, namely Imagery Workspace and Map Workspace. Imagery Workspace as can be seen in Figure 2 is a workspace that is used to perform stereo satellite imagery (stereopair) processing which has no spatial reference (raw data). Map Workspace is a special workspace designed to perform data processing and data analysis that already has a spatial reference (geospatial data), either raster data or vector data as can be seen in Figure 3.

Data Processing Simulation
To test the application of DSM extraction method that has been applied in the "RFM based DSM Extractor", in this research has been performed data processing simulation to perform DSM extraction using the software prototype of "RFM based DSM Extractor" based on the stereopair of WorldView-2 obtained from DigitalGlobe Foundation as mentioned before. Through the RFM based cross correlation applied in the software prototype has been obtained 684 conjugate points automatically in which the 110 mismatch points have been removed manually as seen in Figure 4. If observed in more detail, the distribution of conjugate points obtained from RFM based Cross Correlation tends to cluster in certain areas of high contrast as well as sharp borders, such as residential areas as can be seen in Figure 4. Therefore, manually collection of conjugate points is needed to spread the conjugate points, especially on the parts of VIPs (Very Important Points) which are important parts needed to properly represent the surface, such as ridge, valley bottom, etc, as can be seen in Figure 5.
With the conjugate points manually placed on VIPs such as these locations, it is expected that the resulting DSM can represent the surface of the area maximally. Through this manually collection has been obtained conjugate points as much as 5500 points. Thus the total number of conjugate points obtained are 6184 points with the distribution as can be seen in Figure 6.
After the conjugate points are collected, the next step is to perform 3D reconstruction to obtain 3D ground coordinates (longitude, latitude, heigh) of each conjugate point. 3D reconstruction applied in the software prototype is based on Forward RFM, otherwise known as Forward RFM based Virtual Intersection that can be done by the least squares adjustment iteratively. From the calculations that have been done in 10 times iteration is obtained a fairly satisfactory results, in which the RMSE graph shows a fairly good convergence as can be seen in Figure 7. In Figure 7 it can be seen that the RMSE value from iteration to iteration is getting smaller (convergent). The RMSE value has stabilized at the 2 nd iteration.  After each conjugate point has 3D ground coordinates (latitude, longitude, heigh), the next step is to build DSM (Digital Surface Model) based on the conjugate points through interpolation. The interpolation method applied in the software prototype is Ordinary Kriging. In the interpolation process, this Ordinary Kriging method will consider the existing spatial autocorrelation. In this case, spatial autocorrelation will be represented through a semivariogram which is a graph of semivariance value as a function of distance. This semivariogram will show the degree of dissimilarity of the value between points as a function of distance. The farther distance between points, the level of dissimilarity will be greater. Vice versa, the closer distance, the level of dissimilarity will be smaller. It is called spatial autocorrelation, which in this case is represented through the semivariogram. There are 2 types of semivariogram models, the first is empirical semivariogram model, while the second is the theoretical semivariogram model. Empirical semivariogram is built based on empirical data, while the theoretical semivariogram model is built based on common functions that exist to set the value range, nugget and sill to fit the empirical semivariogram model. General function of semivariogram model to choose from, among others circular, exponential, gaussian, linear, or spherical.
Based on these conjugate points, automatically the software prototype build an empirical semivariogram model which in this case divided into 12 lag bins as can be seen in Figure 8. Based on the empirical semivariogram model, the theorectical semivariogram model that is considered most appropriate is spherical model. The range, nugget and partial sill values applied can be seen in Figure 8. The theoretical semivariogram model will be used as the basis of weighting each neighborhood points that will determine the level of influence on interpolation point in ordinary kriging. In this case, the value at each interpolation point can be estimated with the average of the total value of neighborhood points multiplied by the weight of each. In this case the search area of neighbourhood points around the interpolation point will be divided into 4 sectors, where the neighborhood found in each sector will be used at least 2 and maximum 5. From the kriging interpolation process using this configuration obtained DSM in raster format with spatial resolution (cell size) of 5 meters as can be seen in Figure 10 and Figure 11.

Accuracy of the Generated DSM
Based on 11 Independent Check Points (ICPs) it can be seen that the generated DSM has great absolute accuracy of 2.47639575534581 meters in planimetry (CE90) and 4.21487430387638 meters in height (LE90), which the discrepancies (RMSE) of ICPs on the ground can be seen in Figure 9. When compared to the WorldView-2 spatial resolution of 0.5 meters, this absolute accuracy is still not satisfy. This is because the RPCs provided by the vendor still contain various biases and systematic errors. Therefore in the next research is strongly recommended to apply the RPCs bias compensation to improve the absolute accuracy.
Based on visual observations on the contour map, elevation profile and bird's eye view derived from generated DSM, it can be said that the DSM generated through the WorldView-2 stereopair has good relative accuracy which is visually good and representative in representing the terrain surface they represent as can be seen in Appendix 1.