Abstract

In this paper, we have generated a 5m Digital Elevation Model (DEM) from GeoEye-1 along-track stero pair satellite imagery. The Rational Polynomial Coefficients (RPCs) provided by vendors are biased using the Rational Function Model (RFM) to improve the geo-positional accuracy. The texture features and edge features are extracted to efficiently identify the conjugate points. Epipolar resampling is performed and Normalized Correlation Coefficient (NCC) is used as a template matching technique. The patch transformation using the slope and aspect analysis is performed to match the conjugate points and whose success rate is more than without using patch transformation. The points which satisfy the defined correlation factor are accepted for DEM construction. The RMSE at five different random points of generated DEM are calculated using Google earth data and found better than LPS eATE algorithm.

1. Introduction

Digital Elevation Model (DEM) is a 3-dimensional digital representation of Elevation of Earth's surface. DEM is an essential data in many remote sensing applications like urban planning, hydrological modeling, military and security, mining and transportation, environmental studies, etc. It can also be used to produce orthophoto maps, perspective maps and contour maps.

The conventional methods of extracting DEM including digitizing topographic contour maps, manual outlining of stero modeling, performing autocorrelation for photogrammetric system and conversion of vector hypsographic data are discouraged due to expensive production procedure. Hence, the automatic extraction of DEM from stero satellite data using photogrammetry has been an active research topic from the past few years.

Stero-pair images are images of a point on the earth's surface acquired with different viewing angles, which is used to depict the elevation information. The stero images acquired in same earth's orbit are called as along-track stero satellite images. If the images are acquired with different orbit space are called as across track stero satellite images. In our proposed method, GeoEye-1 along-track stero satellite image is used to produce DEM.

Initially, researchers are focused on establishing a geometric relationship between the object space and image space by knowing the interior and exterior orientation of the physical sensor model. Though it produces high modeling accuracy, explicit understanding of each of the parameters of physical sensor model is difficult and it depends on the set of Ground Control Points (GCPs) to recover the exterior orientation parameter. Generalized senor model such as RPC model (Rational Polynomial Coefficient) has overcome this limitation. Vendors explicitly provide RPCs to end-users for photogrammetric processing. 0.5m resolution of GeoEye-1 stero data with RPC is used in the proposed work.

DEM extraction based on feature matching technique broadly involves feature extraction, feature matching and image to ground projection. The major limitation of feature-based DEM extraction are missing features due to illumination effects and occlusions leads to inappropriate DEM. Feature matching techniques are well sophisticated and computationally expensive. DEM extraction algorithms MATCH-T land INDUSURF leads poor result for texture variation regions and steep slope regions. Hence, we adopted both Feature-based technique (Textures) and patch transformation for Epipolar resampled images to improve DEM accuracy and time complexity.

The proposed method consists of three main modules: The first module consists of RPC biasing and feature extraction. The RPC biasing using generic sensor model such as Rational Function Model (RFM) pre-process the inappropriate altitude and satellite orbit measurements provided by image vendors. The earth's surface features such as edges and textures are extracted to accurately determine the conjugate points in both the epipolar images. The edges are extracted using canny edge detection algorithm and texture features are extracted using Gray Level Co-occurrence Matrix (GLCM).

The second module consists of Epipolar resampling: the Epipolar resampling is to establish the geometric relation between two images and an object (object space and image space). The construction of Epipolar geometry reduces search time in matching the conjugate points by aligning the conjugate points along the Epipolar line.

The third module is image matching using patch transformation. The Normalized Correlation Coefficient (NCC) is used as a template matching technique to extract the conjugate features along the Epipolar line of each layer in the pyramid image. The features that are not in the range of defined threshold are subjected to slope pattern analysis based on x-parallax. It involves repeated transformation from image to ground within the matching window.

The main contributions of the proposed method are:

1. Adopting RPC biasing improves geo-positional accuracy.

2. The texture features and edge features enable to accurately detect the conjugate points along Epipolar line.

3. The patch transformation using Normalized Correlation Coefficient (NCC) improves the matching speed and matching success rate.

4. Slope analysis model ensures pixel level accuracy.

The paper is organized as follows: Section 2 covers related work, which is a brief of existing works on extracting DEM using feature-based technique, Epipolar resampling and transformation techniques. Section 3 describes in detail about the system architecture adopted to extract DEM. Section 5 gives a detailed description of the data set used and experimental results. Section 6 is a discussion and Section 7 concludes the proposed method.

2. Related Work

In this section, we reviewed the papers related to generation of DEM from different across-track and along-track stereoscopic satellite images with GCPs and without GCPs.

Mojtaba Jannati et al. ¹ proposed a novel method to resample the Epipolar images of across-track linear pushbroom satellite images based on a rigorous sensor model. The z-parameter of the stero images are modified using polynomial function then the rotational operation is performed to make Conjugate Image Points (CIPs) and Baseline points (BPs) parallel. This, in turn, resamples hyperbola curve along scanline to straight lines and thus eliminates y-parallax. The main advantage of this method is, it needs a few CIPs and zero GCPs. The feasibility of the proposed method is validated on different terrain conditions of SPOT-1A and RapidEye stereo imagery. It fails to explain the model on along-track pushbroom stero images, normalization and generation of DEM and orthophoto for normalized images.

The y-parallax of high resolution across track satellite images are minimized using frame-based rectification method. The dense features are extracted and Fast Dense (FD) feature matching model was proposed by Wen-Liang Du ². It reduces the computational complexity to linear and compared against LLT, SOCBV and WGTM algorithms. The method is not tested for automatic Epipolar resampling of across/along-track pushbroom stero data.

The automatic extraction of DEM from across-track stero images was proposed in ³. The patch transformation for slope pattern was used to find the conjugate points. The overall matching rate was improved to 5% but elevation error was high when it was tested against LIDAR DEM.

The RPC bias compensation model for hierarchical image matching was developed to generate the DEM from worldview 1 stero satellite image ⁴. The NCC was adopted to image pyramids from coarsest to finest level to remove all outliers. The weighted adjustment technique was used to extract the relative bias terms which in turn improves the matching rate. The procedure was iteratively repeated for all the layers. This model decreases computational complexity and increases the success rate.

The accuracy of DEM from different stero satellite image (IRS-1D and IKONOS) for three different areas (Egypt, Toronto and Hong Kong) was studied using three different models (RFM, 3D-affine and polynomial model). Initially, sensor orientation process was performed followed by DEM extraction model, the 3D-affine achieved better accuracy than the polynomial model. 3D-polynomial achieves less than pixel-level accuracy but needs more GCPs whereas, higher accuracy with only one GCP is achieved by RFM. Overall, RFM was found to be more robust than 3D-affine and 3D polynomial ⁵.

Thierry Toutin et al ⁶ extracted Digital Surface Model (DSM) from High Resolution (HR) Radarsat-2 stereo-images of rural and hilly regions of Quebec, Canada. He used three dimensional empirical deterministic and hybrid approaches such as Toutin's Model (TM), Hybrid Toutin's Model (HTM) and empirical RFM model without GCPs. The obtained DSMs are compared with standard Lidar elevation point cloud data set. TM and HTM produce slight higher performance concerning elevation, systematic and random errors. HTM without GCP shows slightly less accuracy with respect to TM.

Paper ⁷ studied the changes in glaciers and ice caps of Nunavut, Canada by extracting Glaciers DEMs using RADARSAT-2 stereo data (Stero radargrammetry-2 Technology). The HTM and RFM geometric models with no GCPs were used to access the quality of DEM. The RFM shows better vertical bias error than HTM. The SRG-2 technology is suitable to monitor glacier changes in a cloudy and remote region.

The DEMs of high latitude and Polar Regions with varying textures due to shadow cover and ice melt was generated by Myoung-Jong Noh et al ⁸. The Triangulated Irregular Network-based Search-space Minimization (SETSM) algorithm was developed with RPCs of worldview 1 and 2 stero pair images. The elevation accuracy of the generated DEM was validated with the Lidar point cloud data set. Approximately 20cm RMS error was reduced from 5m actual DEM.

et al ⁹ improved the geo-positional accuracy of ZY-3 Satellite Imagery using statistical learning theory without the use of Ground Control Points (GCPs). The RFM and affine transformation functions are used for bundle block adjustment. Fast Iterative Shrinkage-Thresholding (FIST) algorithm is used as statistical theory approach to define and solve the equation matrix. The residual errors were minimized using DBSCAN, the proposed method improves the geo-positional accuracy of ZY-3 Imagery.

The baseline of Interferometric Synthetic Aperture Radar (InSAR) was calibrated to generate GF-3 real-time DEM. The elevation error of generated DEM is directly proportional to baseline estimation hence, the orbit fitting was executed to eliminate nonlinear error effects. The height errors were generated to estimate real-time baseline errors. Finally, DEM was generated based on the modified baseline. The accuracy of the generated DEM was validated by selecting random GCPs of google earth ¹⁰.

Tianen Chen et al ¹¹ derived DEM from ALOS PRISM stero pair images. Initially, the image pyramids are constructed and point features are extracted using intersect operator. Relaxation matching technique for grid-based was used to extract conjugate point. Generalized bundle block adjustment was performed to construct DEM.

As consult with the published literature, most of the works adopted RFM for bias compensation and Epipolar resampling for geometric calibration. Pyramid layer construction reduces complexity in matching technique. Feature-based DEM extraction alone found a poor result. Hence, we adopted a pyramid layer strategy and feature-based template matching using NCC for Epipolar resampled images to accurately extract DEM.

3. Methodology

The system architecture of the proposed method is as shown in Figure 2. It consists of feature extraction and Epipolar resampling in first phase. Template matching using NCC is second phase. To improve the matching rate, the patch transformation based on slope and aspect analysis is proposed in the third phase.

The high resolution along-track stereoscopic satellite imagery consists of RPCs which are biased to accurately align the three-dimensional information (latitude, longitude and height) from object space to two-dimensional information (line and sample) in image space. Biasing is an essential pre-processing step in generating accurate Digital Elevation Model (DEM). RPCs provided by image vendors are biased to improve the geo-positioning accuracy. The Rational Function Model (RFM) is used to bias the RPCs using Ground Control Points (GCPs). It is a ratio of two cubic polynomial functions as described in equation 1.

(1)

(2)

(3)

Where, ‘L’ and ‘S’ are line and sample coordinates in image space. F₁, F₂, F₃ and F₄ are polynomial equations of degree 3. ‘Ln’, ‘Sn’, ‘X’, ‘Y’, ‘Z’ are normalized coordinates. (L₀, S₀, a₀, b₀, c₀) are corresponding offset and (L_s, S_s, a_s, b_s, c_s) are corresponding scale terms respectively.

The RPCs provided by image vendors are erroneous. Hence, to avoid complexity, shift terms (correction terms) ( set zero orders of denominator polynomials to unity.

(4)

(5)

The ‘n’ Ground Control Points are used to bias RPCs. The local polynomial model is used to model the correction terms which is given by (6) and (7).

(6)

(7)

For ‘n’ GCPs, M=QN.

Where, are the coordinates of the image point being processed. ‘Q’ is () the coefficients of polynomial equation. ‘N’ is the image-coordinate vector of GCP. ‘M’ is the vector of polynomial coefficients which are calculated using weighted least squares approach ^{12, 13}.

In the proposed method, the Epipolar resampling has three main advantages: it reduces time in searching conjugate feature in the right image, ensures zero- y-parallax and maintains proportionality between x-parallax and height, hence x-parallax itself is used for slope estimation.

Typical Epipolar geometry is explained in Figure 1. It is the intersection of two image planes with a common baseline. The baseline is the line joining the camera centres. The object ‘C’ on earth surface is viewed from two different cameras. ‘A’ and ‘B’ are the projections of camera centres. Epipolar plane is the plane consisting of these projection points and object point ‘C’. ‘C_L’ and ‘C_R ‘ are the projection points of ‘C’ on image planes. The image points B_L and B_R are called Epipolar points. The projections ‘C1’, ‘C2’, ‘C3’ passes through a common line, which is an intersection of the image plane and Epipolar plane and the line is called an Epipolar line. All Epipolar lines intersect at an epipole. This Epipolar geometry is motivated by contributing to searching the conjugate point in the matching process.

Existing studies proved that the Epipolar line of pushbroom satellite sensor images is hyperbola like curve because of non-linearity of Epipolar planes ¹⁴. Hence, the Epipolar curve is to be approximated over an entire image to find the conjugate features efficiently. The piecewise linear Epipolar resampling ¹⁵ is adopted in the proposed method followed by iterative transformation of linear segments of Epipolar curves from object space to image space. Select point ‘p’ in the left Epipolar image, the corresponding control point’s ‘q₁’ and ‘q₂’ for right Epipolar image is constructed for maximum height ‘h₁’ and minimum height ‘h₂’ respectively ¹⁶.

Curve between ‘q₁’ and ‘q₂’ is approximated as line segments. Second order polynomial equation is used to transform approximated points from object space to image space. Similarly, relation is formed for ‘q₁’ in right image space. Procedure is iteratively performed for all the points until Epipolar images are entirely delineated.

The High-resolution stero images of Geo-Eye 1 is several hundred bytes of data and hence the images are decomposed into several levels based on coarse-to-fine strategy ¹⁷. This pyramid layer construction is best suitable for image matching algorithm. This also reduces repeated initializations in the matching algorithm.

The texture features help to identify conjugate points in Epipolar images. The texture features for stero images are extracted using pixel coordinates. Gray Level Co-occurrence Matrix (GLCM) is used to extract different texture features (mean, variance, energy, and entropy) shown in Table 1. It differentiates surface objects based on density, shape, size, and illumination. The Canny edge detection technique is applied to the stero images to extract the edges of an object.

The conventional feature matching techniques such as point feature matching technique, Speed Up Robust Feature (SURF) and Scale Invariant Feature Transform (SIFT) have large computational complexity and not suitable for real-time applications. Therefore, in the proposed method Normalized Correlation Coefficient (NCC) is used to find the conjugate points based on features. The efficient points are extracted based on the correlation factor and are given in equation 8. The correlation factors greater than a defined threshold are accepted for DEM construction.

(8)

Where, ‘M’ is the window of the left Epipolar image patch, ‘N’ is the window of right Epipolar image patch. ‘M_a,b’ and ‘N_a,b’ are the components of image patches ‘M’ and ‘N’ respectively, at line ‘a’ and sample ‘b’. ‘M’ and ‘N’ are the mean of kernel window of image patches ‘M^|’ and ‘N^|’ respectively. The ‘p’ is the window size of image patches.

The correlation factors less than the defined threshold are subjected to main slope analysis. The y-parallax is reduced to zero in Epipolar resampling and x-parallax is directly related to height information. The x-parallax are used for slope analysis.

Step 1: Resampling of x-parallax using bilinear interpolation technique:

Without the resampling of x-parallax, the slope gradient yields poor result, so x-parallax is refined using a bilinear interpolation technique. Bilinear Interpolation Technique (BIT) results in sharp DEM and it delineates steeper downslope.

The BIT produces the output raster by calculating the weighted average of the cell values of four nearest neighbors. It results in a smoother-looking surface ¹⁸.

Step 2: Slope estimation is performed in this step. It is the maximum rate of change in elevation information. For each cell, the rate of change of elevation between the current cell and its neighbors is calculated. The steepest downhill drop is extracted by calculating the maximum rate of change of elevation from the current cell to its eight neighbor cells over a distance cover.

The average maximum technique is used to calculate the slope value for the selected 3×3 window. The aspect of the processing cell is the direction that the plane faces. If the slope value is low, then the terrain is flat and vice versa.

The z-value for the selected cell will be assigned to the cell location if the neighborhood cell has NoData or low-value. The above method is applied only when, at least three cells (outside the raster's extent) at the edge of the raster contains NoData. This results in fall in the slope due to flattening of the 3 x 3 plane fitted to these edge cells.

Step 3: based on the drop in slope and/or slope value (degree of slope) the direction of the slope is estimated using aspect parameter. Aspect is the downslope direction of the slope value. It can also be called as the slope direction. The aspect value indicates the compass direction that the surface faces at that location. It is measured clockwise in degrees from 0-360 (due north to again due north). -1 is assigned to the flat areas having NoData in z-values. The rate of change of elevation in a horizontal and vertical direction at point ‘α' is given by (equation 9 and 10)

(9)

(10)

Where, l,m,n…t are the cell values of 3×3 matrix. dZ/dX and dZ/dY are the rate of change of elevation value in horizontal and vertical directions respectively.

The rate of change in elevation in both the direction is given by. (equation 11)

(11)

The aspect value is converted into compass direction values based on the following set of rules.

if aspect<0

Compass_value=90.0-aspect

else if

Aspect>90

Compass_value =360.0-aspect+90

else

Compass_value=90.0-aspect

end if

end if

Slope value is found to be constant along the line direction and varying along with the sample. The aspect value (Compass Value) is grouped into four bins based on the degree of direction from 0-360.

The aspect values belonging to Bin-1 is the slope facing up or down in the left image, its consequent slope direction is the same for its conjugate point in the right image with distance ‘d’. Hence the point ‘q’ is scaled by considering ‘d’ factor. The distance ‘d’ is determined using Euclidian’s distance method.

The coordinates of left image for point ‘p’ is (x,y). Let the scaling factor be (S_x,S_y). The resulted coordinates are (x^|, y^|). Mathematically it is given by,

(12)

Or

The aspect value belonging to bin 2 is the slope facing towards the right or left, in the left resampled image. The consequent slope direction in the right image is found to be in the same direction with a displacement of x-coordinate system. Hence, the shear transformation is performed by considering the displacement factor ¹⁹.

The shear transformation for x coordinate with displacement ‘S_h ‘ by preserving y coordinate is given by,

(13)

The aspect value belongs to bin-3, is the slope facing diagonal direction (top right to bottom left) and aspect belonging to bin-4 is the diagonal facing top left to bottom right in the left image, its consecutive point is found to be in the same direction with a displacement of angle theta. The rotational operation is performed to the left image to extract its conjugate points in the right image.

In rotational transformation, the rotation of a point is performed concerning its coordinates at a particular angle ‘θ' (theta) from its origin ²⁰.

The translated points from the above-discussed process are subjected to template matching technique (NCC). The points whose correlation factor greater than the threshold are accepted for DEM construction. The unsatisfied points are again subjected to the main slope analysis. This procedure is iteratively repeated until all the pixels are done.

4. Experimental Results

The Geo Eye-1 along track stero satellite imagery with 0.5m resolution was acquired from Digital Globe. The 8-band ortho-ready panchromatic stereoscopic imagery of typical urban scenario with RPCs file is used for the study. The study area covers 100sqkm with stero overlap of 19 with elevation and azimuth angle. The stero images are shown in Figure 6.

Initially, the RPCs of Stero satellite imagery are biased concerning ground control points by adding shift terms. The RPC bias correction drops the geo-positional error from seven pixels to 0.3 pixels.

The pyramid layer was constructed using a coarse-to-fine strategy, staring with 50m resolution (lower pyramid layer) to the highest pyramid of 3cm resolution. Lower pyramid does not contain too much information and unable to perform transformation during a matching technique. Next higher layers have large information and perform transformation by dividing the layers into several patches of size 3×3. The resampled Epipolar lines and edges using canny edge detection algorithm for 7×7 window is shown in Figure 3.

The four different texture features and edges are extracted for each window patch of size 3×3 for all the pyramid layers. Totally 1400968850 features are extracted. Resampled Epipolar images and features enable us to extract the conjugate points efficiently. NCC is used as a matching technique, the resulted correlation factors were in the range of -1 to +1. The threshold of 0.8 is fixed, values greater than 0.8 was found to be accurate.

Without slope analysis and patch transformation, 1000470011 conjugate points were extracted. With main slope analysis and patch transformation, 1400088510 conjugate feature points were extracted with an improved success rate of 97.5%. The observed 26 tie points (Figure 4) are extracted using LPS eATE algorithm to validate our extracted conjugate points.

The patch transformation in slope analysis based on distance, displacement and angular rotation is explained in Figure 5 and the 5m generated DEM is shown in Figure 7.

The Root Mean Square Error for the constructed DEM at five different random points are calculated using google earth data. The contour map for the generated DEM is generated which is shown in Figure 8. We extracted DEM for the same stero data using LPS eATE algorithm and checked the RMSE for the same points, we found that elevation error for the proposed method at four points out of five points is producing lower error than LPS eATE algorithm, which is tabulated in Table 2.

5. Discussions

Stero-pair images of 0.5m with RPCs are acquired from DigitalGlobe. RPC biasing is performed using shift terms of RFM model, which improves the pixel positional accuracy concerning GCPs. As the image consists of thousands of bytes, it has a high risk of computational complexity. The pyramid layer image is constructed using a coarse- to- fine strategy. Lower layer of pyramid image bears the least information whereas higher layer contains more information.

The texture features and edge features are extracted for all the layers using GLCM and Canny edge detection algorithm respectively. These features make us extract the conjugate points efficiently. The Epipolar resampling is carried out in an Epipolar geometry to align the conjugate features on to the Epipolar line. The NCC is used as a template matching technique. The points whose correlation value greater than 0.8 are accepted for DEM construction. Residual points are subjected to the main slope and aspect analysis module.

Initially, the degree of slope value is calculated, then the direction of the slope is calculated using aspect parameter. Based on the aspect values (0.0-360.0) four bins are constructed. The four bins possess eight different directions for each point. The value belongs to bin-1 has a direction facing up or down, bin-2 has a direction facing either left or right, bin-3 and bin 4 has direction facing diagonally, top right to bottom left and top left to bottom right.

Based on the direction and displacement (distance, angle, and shift), the transformation (scaling, shearing or rotational) is performed with the selected point has a center cell of a 3×3 matrix of the image patch.

The NCC is applied to all the transformed points of each image patches in the pyramid image layers. The correlation greater than 0.8 are accepted and the procedure is iteratively performed until all the points are done. Finally, the RMSE of the generated 5m DEM with the google data is calculated and found a better value than LPS eATE algorithm.

6. Conclusions

We extracted 5m DEM from GeoEye-1 stereoscopic satellite imagery. The texture feature extraction using GLCM and Epipolar resampling allows us to extract efficient conjugate points. The RPC biasing reduces the geo-positional error, thereby it improves the matching success rate during the template matching process. The patch transformation in slope and aspect analysis explicitly improves the matching rate than without patch transformation. The result shows that the success rate was improved by up to 97%. We proved that NCC can also be used for feature matching in Epipolar geometry. The accuracy of the generated DEM is calculated using RMSE parameter. Google earth data is used for RMSE calculation. The contour map is generated for the constructed DEM. The proposed method performs better than LPS eATE algorithm.

In the future, we try to use edge feature effectively in the matching process to find the conjugate feature for the shaded region.

References