Generation of High Resolution Global Dsm from Alos Prism

Panchromatic Remote-sensing Instrument for Stereo Mapping (PRISM), one of onboard sensors carried on the Advanced Land Observing Satellite (ALOS), was designed to generate worldwide topographic data with its optical stereoscopic observation. The sensor consists of three independent panchromatic radiometers for viewing forward, nadir, and backward in 2.5m ground resolution producing a triplet stereoscopic image along its track. The sensor had observed huge amount of stereo images all over the world during the mission life of the satellite from 2006 through 2011. We have semi-automatically processed Digital Surface Model (DSM) data with the image archives in some limited areas. The height accuracy of the dataset was estimated at less than 5m (rms) from the evaluation with ground control points (GCPs) or reference DSMs derived from the Light Detection and Ranging (LiDAR). Then, we decided to process the global DSM datasets from all available archives of PRISM stereo images by the end of March 2016. This paper briefly reports on the latest processing algorithms for the global DSM datasets as well as their preliminary results on some test sites. The accuracies and error characteristics of datasets are analyzed and discussed on various fields by the comparison with existing global datasets such as Ice, Cloud, and land Elevation Satellite (ICESat) data and Shuttle Radar Topography Mission (SRTM) data, as well as the GCPs and the reference airborne LiDAR/DSM.


INTRODUCTION
The elevation data which represent the terrain on the ground is one of fundamental layers in the field of geographic information systems.It is essential not only for creating maps or ortho-photos but also for various applications e.g., infrastructure design, disaster monitoring, environmental monitoring, and natural resources survey.In recent years, global elevation datasets derived from spaceborne instruments are being widely used.The advantages of the satellite derived data are its wide range, continuity of quality, total cost, and so on whereas the aerial one has advantages in the spatial resolution and absolute accuracy.The methods to generate elevation data are roughly categorized as the traditional optical stereo matching, the radar interferometry, and the Light Detection and Ranging (LiDAR).The Shuttle Radar Topography Mission (SRTM) was the first global datasets made by spaceborne radar instruments and was released in 2003 first.The data has 3 arcsec (90 m) pixel spacing in which the absolute and relative height accuracies are ~9 m and ~10 m respectively (90% errors) (Rodriguez et al., 2006).The global elevation data derived from Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER/GDEM) is one of the datasets generated with the optical technique and was released in 2009 first.It has the height accuracy of 13 m (1) in 1 arcsec (30 m) pixel spacing (Tachikawa et al., 2011).Ice, Cloud, and land Elevation Satellite Geoscience Laser Altimeter System (ICESat/GLAS) measured land elevations from 2003 to 2009 by using the space-based LiDAR instrument.The footprint of the laser data is about 65 m in diameter at 172 m intervals along-track, namely, the data is point information rather than the raster data (Schutz et al., 2005).Its height accuracy is less than 1 m depending on characteristics of the ground (Duong et al., 2009).Panchromatic Remote-sensing Instrument for Stereo Mapping (PRISM) was an optical sensor onboard the Advanced Land Observing Satellite (ALOS) operated from 2006 to 2011, and was designed to generate worldwide elevation data.The sensor consists of three independent panchromatic radiometers for viewing forward, nadir, and backward in 2.5 m ground resolution producing a triplet stereoscopic image along its track (Tadono et al. 2009).Though the satellite already ended its mission life it left us huge amount of stereo images all over the world (Takaku et al. 2013).To utilize all these archive data we proceed to generate new global elevation datasets, named 'ALOS World 3D', which have finer ground resolution and higher accuracy than those existing ones.The data is Digital Surface Model (DSM) and its target height accuracy was set to 5 m (rmse) in the pixel spacing of 0.15 arcsec (approx.5 m).Table 1 shows the basic specification of the datasets.The process-chain of datasets should be full-automatic because the number of available archive data is approx.one million sets of stereo images which cover 35km square each on the ground.One of the most challenging tasks in the process-chain is the masking of outlier areas covered with water, snow, cloud, etc.This paper briefly reports on the latest algorithms for the fullautomatic processing of global DSM datasets with preliminary results on some test sites.The accuracies and error characteristics of datasets are analyzed on various fields by the comparison with existing global datasets such as ICESat point data and SRTM data, as well as the comparison with Ground Control Points (GCPs) and a reference airborne LiDAR/DSM.

DATA PROCESSING
The software for DSM generation, which was exclusively developed for PRISM, is based on a traditional fashion of photogrammetry; it processes the data in a unit of 35km x 35km scene first, and then, they are mosaicked onto 1°x1° tiles of geographic latitude/longitude grids.Finally, the quality of data is inspected automatically by using reference information such as ICESat and SRTM as well as the visual check by human operators.

Scene data processing
For scene data processing we have developed software called DSM and Ortho-rectified image Generation Software for ALOS PRISM (DOGS-AP) (Takaku et al., 2009a) and semiautomatically processed DSM data of more than six thousand scenes in some limited areas as of 2013.The software briefly consists of the tie-point generation, the image orientation, the image matching for height-calculation, and the masking of outlier areas.Figure 1 shows the processing flow of scene data.It is equipped with a unique matching engine exclusively developed for triplet stereo images of PRISM.It also generates the Ortho-Rectified Image (ORI) data of nadir image on original 2.5m pixel spacing besides the DSM on 5m pixel spacing.The software does not need any GCPs to give the planimetric accuracy of 5m (rms) thanks to the direct geolocation accuracies of the PRISM physical sensor model.To implement the software into the processing-chain of global data generation it should be fully automated for mass processing.One of the most challenging tasks in the automation is the masking of outlier areas covered with water, snow, or cloud.In the process we apply a statistical classification algorithm analyzing the correlation derived from the image matching, the brightness value of the ORI data, and the roughness of DSM data itself.For open water areas the existing global water-bodydata in public domain such as SRTM Water Body Data (SWBD) (NASA/NGA, 2003) or Global Self-consistent, Hierarchical, High-resolution Shoreline Database (GSHHS) (Wessel et al., 1996) are utilized as initial masks.However, the initial masks may have some incorrect areas because their source dates were different from the ones of PRISM data.Moreover their ground resolution is larger than the one for PRISM.Therefore first the excessive masks in the initial waterbody masks are deleted by analysing the matching correlations on the areas.In the analysis a function R C representing a rate of pixels which have reliable correlations in local area is calculated as follows: where, (x, y) is the pixel address, N is the size of local area, Cr is the function of correlations, T C is the reliability-threshold of correlations.The matching algorithm for triplet images calculates two kinds of correlations per pixel: the one between nadir and forward and the other one between nadir and backward, so that the R C is calculated for each of them.The mask pixels which should be deleted are then selected where both of R C s (R Cf and R Cb ) are greater than threshold T R , namely, .
T C , T R , and N are set to 0.6, 0.7, and 33 respectively for the deletion based on preliminary experiments.The remaining mask segments are then trimmed with the morphological operation to delete small bits and to fill tiny holes.Figure 2 shows one of the results of the deletion for excessive masks in the SWBD.
The lack of masks in the initial existing masks is compensated in the next step which adds the cloud and snow masks as well as the water masks.In the step the masking areas are extracted first also by using Rcs regardless of whether they are cloud/snow masks or water masks.The extracting criteria is contrary to the one for the deletion namely, .
T C , T R , and N are set to 0.2, 0.8, and 33 respectively for the addition based on preliminary experiments.The extracted masks are then classified into the water or cloud/snow mask after morphological operations.The accuracy of the classification is important because the voids on cloud/snow masks in each scene can be possibly filled with other scenes, which were observed at different dates, in the following stacking process while the voids on water masks basically remain and are interpolated finally with their surrounding valid heights.The classification is performed based on statistical analyses of the image-brightness and of the DSM-height in each mask segment.The average, standard deviation, and maximum of image-brightness, and the standard deviation of DSM-height are calculated for each segment.Basically, high values of these parameters imply that the corresponding segment should be a cloud/snow mask, while low values of them imply that the corresponding segment should be a water mask.These thresholds of parameters were tuned and combined through preliminary experiments, and then, are used for the classification.Figure 3 shows the results of the added masks for water and cloud areas.

Tile data processing
The DSM data processed on a scene basis are then automatically mosaicked on 1°x1° tiles in 0.15 arcsec (approx.5 m at equator) spacing as final dataset products.The number of 1°x1° tiles covering global land area is about twenty-three thousand.Figure 4 shows the basic processing flow of the tile data.The DSMs generated without GCPs include absolute vertical shift-errors to be corrected on a scene basis due to instability of sensor alignment angles for forward-and backward-viewing radiometers (Takaku et al. 2009a); they are corrected with existing global height control reference such as ICESat or SRTM, while maintaining the relative accuracy of details.The SRTM is used in the low and middle latitude areas while the ICESat is used in the high latitude areas.In the stacking/mosaicking process the grid data of scenes are stacked on each grid of the tile.The grid data of height are averaged after discarding outliers if the multiple data of different dates of scenes are available, while filling the void (cloud/snow) mask areas mutually.The detection of outliers is based on the majority decision and is performed if the number of data exceeds two.Next the height data in inland-water masks (i.e., rivers and lakes) are interpolated with surrounding valid height before the tile framing.The correlation data of scenes are also mosaicked for quality reference and are output as well as a map of stacking numbers.Finally the accuracy of every DSM tile is evaluated with globally available ICESat data and is annotated in the product metadata.To check the large matching errors, obvious missing masks, or planimetric systematic errors SRTM or ASTER/GDEM are used for the reference information as well.

VALIDATION
We preliminarily processed four 'ALOS World 3D' 1°x1° tiles which have various terrain-or texture-types for the validation.The tile-IDs which refer to the latitude/longitude of lower-left corners for these tiles are N27E089, N16W016, N11E105, and N36E140.Table 2 shows the basic specification of these tiles.
Figure 5 shows their DSM height data displayed in gray scale with the mask data displayed in colours.Figure 6 shows the

ICESat
ICESat data products GLA14 are used as main absolute reference to validate all tile data (Zwally et al. 2012).In the comparison the height data under the ICESat footprint in the tiles are averaged for each point.Since the size of footprint is much larger than the spacing in the tile data ICESat points which cover relatively flat terrain should be selected in the validation.To identify these points we use the standard deviation of the tile data under each footprint (Huber et al, 2009).The points in which the standard deviations are less than 5m are selected for the comparison.Furthermore, the ICESat data may include some outliers due to the cloud reflections or saturated waveforms.To identify the outliers we use a simple thresholding of height differences (Carabajal et al. 2006).The points that the heights of ICESat are higher in more than 100m than the ones in tiles are excluded from the validation results.Figure 7 shows the ICESat data distribution on tile N27E089 and its height difference histogram.Table 3 shows    while ones of others which include steep terrain are approx.4 m.All these results are enough satisfying the target accuracy of 5 m rms.

SRTM
Though the SRTM is not suitable to validate the detailed relative accuracy of the DSM tiles due to lack of ground resolution and of absolute accuracies, it is useful to validate global areas in perspective.Therefore we use the SRTM for the reference information to check systematic errors in large scale, large matching errors, or obvious missing masks.The SRTM-3 data in 3 arcsec spacing are up-sampled into 0.15 arcsec spacing of DSM tile frames by using the bilinear interpolation so that their difference image can be processed.Figure 8 shows the height difference images between the DSM tiles and SRTM-3.
In the visual check of them no obvious systematic errors or blunders were found except for the north east part of N11E105 in which the change of vegetation was the cause of errors.It also was confirmed that errors on ridges or ravines of steep mountains were due to the difference of the grid spacing between these two datasets.Table 4 shows the statistics of the height difference which reflect the trends confirmed in the visual check.

GPS-Track/LiDAR-DSM
For ensuring final detailed accuracies of the products the global GCPs, which were distributed in several areas on every continent and were used in the sensor calibration of PRISM, are used for corresponding tiles as well as the airborne LiDAR/DSM datasets used in the past validations (Takaku et al. 2009b).Only tile N36E140 includes the GCPs and the airborne LiDAR/DSM site among the four tiles.
The absolute accuracy of the GCPs derived from GPS tracks is better than 1 m in both of planimetry and height.The height of the DSM tile at the planimetric position of each GCP is interpolated from its surrounding grid data for their comparison.
The number of available GCPs is 122 in the tile N36E140.
Table 5 shows the statistics of the height difference between the DSM tile data and GCPs.The rmse is approx.4 m and is almost consistent with the result from the ICESat reference.The airborne LiDAR/DSM dataset has the height accuracy of 0.25 m in the ground spacing of 1 m.The height range of the data is approx.0 to 900 m in an area of 8 km x 8 km including Mt. Tsukuba in Japan.The area includes flat terrain e.g., paddy or truck-farm and some manmade structures at the base of the mountain as well.The data are down-sampled onto the DSM tile frame of 0.15 arcsec spacing by using the bilinear interpolation for their comparison.Figure 9 shows the ORI, DSM tile data, and the height difference image between the DSM tile data and the airborne LiDAR/DSM on the 8 km x 8 km site area.The errors depending on different slope angles are evaluated as well as the one in the whole site area.The slope angle in each pixel on the down-sampled LiDAR/DSM is calculated with the third-order finite difference weighed by reciprocal of squared distance (Horn, 1981).Table 6 shows the statistics of the height difference between the DSM tile data and the airborne LiDAR/DSM depending on the slope angles.The rmse in whole area is approx.4 m and enough consistent with results from ICESat and GCP reference.The errors increase as the slope angles become larger, nonetheless, the rmse is approx.5 m in which the slope angles are over 30 degrees.The rmse in flat areas in which the slope angles are less than 10 degrees is approx.2 m and is almost consistent with the results from ICESat reference in flat areas such as the desert or basin.

CONCLUSIONS
The automatic processing algorithm to generate high resolution global elevation datasets, named 'ALOS world DEM', is briefly presented.The height accuracies and their characteristics of the datasets processed with the algorithm are analyzed and discussed on various fields by the comparison with existing global datasets such as ICESat and SRTM, as well as the GCPs and the airborne LiDAR/DSM.The results were enough consistent among these different sources of reference and are almost satisfying the target accuracy of 5 m rms.We plan to process all global archive data of PRISM stereo images by the end of March 2016 and to continue the improvement of the algorithms during the operation of the mass processing.The datasets which have 5 m spacing will be released basically in commercial base, while the low resolution data (e.g., 30m spacing, TBD) will be generated by averaging the original ones and will be provided free of charge.

Figure 1 .
Figure 1.Processing flow of PRISM/DSM scene data

Figure 2 .
Figure 2. The result of deletion of excessive masks in SWBD.Left to right: DSM gray image, nadir ORI, correlation map between FWD and NDR, correlation map between BWD and NDR.Blue lines depict outline of original SWBD while green lines depict outlines of corrected mask segments.

Figure 3 .
Figure 3. Added masks on water and cloud areas.Left to right: water masks on DSM gray image, water masks on nadir ORI.cloud mask on DSM gray image, cloud mask on nadir ORI. Green and red lines depict outlines of segments for the water and cloud masks respectively.

Figure 4 .
Figure 4. Processing flow of PRISM/DSM tile product

Figure 5 .
Figure 5. DSM height data in gray scale with mask data in colours: (a)N27E089, (b)N16W016,(c) N11E105, (d)N36E140.The red, blue, and green in the images indicate cloud/snow, inland-water, and sea masks respectively.Figure7.ICESat data distribution of 3709 points on tile N27E089 and its height difference histogram.

Figure 7 .
Figure 5. DSM height data in gray scale with mask data in colours: (a)N27E089, (b)N16W016,(c) N11E105, (d)N36E140.The red, blue, and green in the images indicate cloud/snow, inland-water, and sea masks respectively.Figure7.ICESat data distribution of 3709 points on tile N27E089 and its height difference histogram.

Figure 6 .
Figure 6.Maps of stacking numbers (left) and correlation (right) in the tile of N36E140.

Table 4 .
Stats of height difference between the DSM tiles and the SRTM-3 (DSM minus SRTM-3)

Table 5 .
Stats of height difference between the DSM tile and GCPs (DSM minus GCP)