THE EFFECT OF PANSHARPENING ALGORITHMS ON THE RESULTING ORTHOIMAGERY

This paper evaluates the geometric effects of pansharpening algorithms on automatically generated DSMs and thus on the resulting orthoimagery through a quantitative assessment of the accuracy on the end products. The main motivation was based on the fact that for automatically generated Digital Surface Models, an image correlation step is employed for extracting correspondences between the overlapping images. Thus their accuracy and reliability is strictly related to image quality, while pansharpening may result into lower image quality which may affect the DSM generation and the resulting orthoimage accuracy. To this direction, an iterative methodology was applied in order to combine the process described by Agrafiotis and Georgopoulos (2015) with different pansharpening algorithms and check the accuracy of orthoimagery resulting from pansharpened data. Results are thoroughly examined and statistically analysed. The overall evaluation indicated that the pansharpening process didn’t affect the geometric accuracy of the resulting DSM with a 10m interval, as well as the resulting orthoimagery. Although some residuals in the orthoimages were observed, their magnitude cannot adversely affect the accuracy of the final orthoimagery.


INTRODUCTION
Satellite optical sensors provide very valuable data about the earth surface through map-making and environmental monitoring.However, these applications are extremely expensive and the scientific community should use this data obtained from available sensors in the best way.A typical example of this is the pansharpening process, i.e. the fusion of multispectral satellite imagery of low spatial and high spectral resolution with the panchromatic imagery of high spatial and low spectral resolution.In recent years, pansharpening has become a task of great importance in the field of data fusion, as demonstrated by the increasing number of scientific contributions to this topic.
Pansharpening addresses the fusion of two optical remote sensing images characterized by different spectral and spatial features.Specifically, a Multi Spectral (MS) image with high spectral but low spatial resolution is considered along with a Panchromatic (PAN) image, which is obtained by sensing a single wide electromagnetic spectrum covering the visible and near infrared (VNIR) frequencies and has complementary characteristics with respect to MS: lower spectral but greater spatial resolution.Due to cost and complexity issues and limited on-board storage capabilities, the multispectral sensor has much smaller aperture than the panchromatic sensor thus reducing the spatial resolution of the sensed multispectral image (Shaw and Burke, 2003).For a typical modern multispectral satellite sensor, this ratio is 1 to 16, i.e., a single multispectral image pixel translates to 4 by 4 panchromatic pixels.
Pansharpening is the fusion of the images captured by the multispectral and panchromatic sensors.The objective of pansharpening algorithms is the generation of a fused product * Corresponding author.
characterized by the spectral content of the MS image and the spatial details of the PAN image.This means that the pansharpened image has the same number of pixels as the panchromatic image and also the same number of bands as the multispectral image.
In this paper are investigated the geometric effects of pansharpening algorithms on automatically generated DSMs and thus the resulting orthoimagery.This investigation is motivated by the fact that on automatically generated Digital Surface Models, image correlation is used to extract matching pixels in overlapping images.It is therefore inevitable that their accuracy is strictly related to image quality.Since the Pansharpening process sometimes results in low image quality, it may eventually affect the DSM generation and then the resulting orthoimage accuracy.
To this direction, an iterative methodology was applied in order to combine the methodology described in the prior work (Agrafiotis and Georgopoulos, 2015) with different pansharpening algorithms and check the accuracy of orthoimagery resulting from pansharpened data.

Antiparos Island
The study area, the small Cycladic island of Antiparos has an area of 35 km 2 and it is very close to Paros Island.It measures 12.5 km in length and 5.5 km in width and has a coastal perimeter of 54 km.Even though the island is almost flat, a few little hills in the centre reach a maximum height of 300 metres.The landscape is rather wild and varied including a main settlement and agricultural fields, while on the west coast there are steep cliffs.

The Pleiades 1B Imagery
Pleiades-1B satellite sensor was successfully launched on December 2, 2012.Built by AIRBUS Defence & Space, the satellite was launched from a Soyuz launcher at the European Space Centre in French Guiana.Pleiades-1A and 1B satellites are phased 180° apart in the same near-polar sun-synchronous orbit at an altitude of 694 km, enabling daily revisits to any location on the planet.The sensor can reach a ground resolution of 0.7m in panchromatic mode and 2.8m in multi-spectral mode in vertical direction.The images provided for this work were acquired in the tri-stereo mode for 3D information.They were provided by Astrium to the Laboratory of Photogrammetry of NTUA after a successful research proposal.According to this acquisition scheme, the satellite rotates around its axis and the HiRi camera scans a target area from three different viewing directions during one pass, thus resulting in an image triplet (Gleyzes et al., 2012).However, only the 2 external images of the tri-stereo were used.The images were acquired on 9 April 2013 in the morning within 22.5 seconds.The average viewing angles of the two selected images are, respectively, 6.30° and 9.47° in across-track direction with respect to the nadir and 1.63° and -12.39° in along-track direction (Figure 1).

Control Dataset
In order to evaluate the produced orthoimages, accuracy specifications had to be set as a reference for the purpose of inferring about the actual metric efficiency of the product.As control dataset aerial imagery orthophotomosaic (Very Large Scale Orthophotos -VLSO) provided by the Hellenic National Cadastre & Mapping Agency S.A. (NCMA S.A) was used.It had 0.25m GSD and the imagery was taken on August 12th 2012.
According to NCMA S.A. the orthophotomosaic of the control dataset was created from scanned in 1700 dpi (1pixel = 15μm) aerial imagery of a scale of 1:15000.The initial aerial images were acquired using a ZEISS RMK TOP 30 film camera having a lens of 153.12mm.To produce this orthoimagery, a DEM was created having 5m grid interval.

PRIOR WORK
In a prior work (Agrafiotis and Georgopoulos, 2015) the accuracy and radiometric quality of orthorectified high resolution satellite imagery from Pleiades-1B and GeoEye-1 satellites were investigated through a comparative evaluation of their quantitative and qualitative properties with and LSO and a VLSO as the benchmark dataset.In addition, the advantages and limits of the Pleiades Imaging for producing Large Scale Orthophotos (LSO) were also investigated.The performed visual assessment of the orthoimagery revealed that Pleiades-1B orthoimagery is especially promising presenting much more information and clearer forms.On the contrary, abrupt changes of brightness and contrast and high radiometric saturation levels were observed on the GeoEye-1 orthophotomosaic.The geometric evaluation reveals that the used LSO and Geoeye-1 orthomosaic suffer from a geometric systematic error in Y axis.Both radiometric and accuracy test results show that Pleiades-1B orthoimage has almost the same absolute accuracy as the orthophotomosaic from aerial imagery (LSO from Hellenic National Cadastre & Mapping Agency S.A.).Hence, it could easily replace aerial imagery, when it comes to orthoimage production.Furthermore, all data are adequate for producing LSO for mapping and GIS, according to JRC (Kapnias et al., 2008) and NSSDA (FGDC, 1998) accuracy standards.The results also serve for a critical approach for the usability and cost efficiency of satellite imagery for the production of LSO.

METHODOLOGY
Having already exploited and measured a benchmark dataset, the orthomosaic of N.C.M.A. S.A. of 0.25m GSD (VLSO), was decided to use the same Control Points described in (Agrafiotis and Georgopoulos, 2015) and compare the measurements with measurements on the new Orthoimagery, produced from pansharpened images.To this end, the classical photogrammetric workflow available in ERDAS Photogrammetry Tool (former LPS) environment (ERDAS 2014) was followed for the processing of the stereo satellite imagery, DSM generation and orthorectification.In addition to the workflow applied in the prior work, pansharpening procedures were performed between the stages of Image Orientation and Image Matching.There, the pansharpened imagery created using HIS, PCA and Brovey transformations was replacing the PAN imagery, changing the file location.In this way, the initial Image Orientation Parameters remained the same for the whole block in order to produce more objective and reliable results.In this way, the only parameter differentiated, is the source image which may lead to different Image matching for DSM generation.
MRA methods are usually based on methods such as the Undecimated Discrete Wavelet Transform (UDWT) or other kinds of pyramidal or multi-scale representations.A more systematic overview of pansharpening methodologies belonging to these two categories can be found in Thomas et al., (2008) and Aiazzi et al., (2012).Additionally, during the last years, remarkable approaches which do not fit in this classifications have been described in the literature.Among those, Bayesian methods based on parameter estimation total variation penalization terms and sparse signal representation are listed (Fasbender et al., 2008;Palsson et al., 2014;Li and Yang, 2011;Zu and Bamler, 2013).

Implemented Algorithms and Results
Initially it was decided to check state of the art pansharpening algorithms such as HCS, IHS and HPF.However, this decision could not be realized since pansharpening procedures in ERDAS Imagine 2014 using the specific Pleiades 1B data produce images of wrong dimensions.To overcome this problem, pansharpening was performed using the Photomod photogrammetric software by Racurs.There, the three commonly used methods are offered, even not belonging to the state of the art: the Brovey transform, the IHS (Intensity, Hue, Saturation) method, and the PCA (Principal Component Analysis) method.
The Brovey transform is based on spectral modelling, while the IHS and PCA methods call upon projection techniques.These methods are available in most commercial software packages for satellite image processing, which may explain their large use by practitioners (Ranchin and Wald).Details on these methods can be found in e.g., Carper et al. (1990) or Pohl and Van Genderen (1998).

Automatic DSM Generation
Even if the same procedure is applied in the processing stage of the prior work (Agrafiotis and Georgopoulos, 2015) the automatic DSM generation stage or Automatic Terrain Extraction (ATE) is described here in order to better understand how the Pansharpening process may affect the extracted DSM and thus the resulting Orthoimagery.
In the literature several techniques are presented suitable for DSM extraction.Some of them are digital aerial and terrestrial photogrammetry, airborne and terrestrial laser scanning, GPS methodology with its different measurement approaches and active and passive remote sensing, with optical satellite imagery systems (Fraser et al., 2002).In particular, aerial and satellite photogrammetry is a powerful tool for surface model generation extracting high resolution DSMs by means of automated image matching procedures.For these tools, DSM can be either in vector or raster format and its density ranges from 1/3 to 1/10 of density of the original image pixels.
The standard procedures to generate DSMs are based on fundamental steps which consist in internal orientation, external orientation (registration into a defined reference system) and point extraction.For the tests performed in the context of this project, the standard procedures described above are the same for all blocks.The only parameter which is differentiated is the used imagery and thus the extracted points.The seven different DSMs were generated in the ERDAS LPS module Classic ATE (Classic Automatic Terrain Extraction) with stereo image matching, which aims at finding dense and robust correspondences between stereo images and estimate their 3D coordinates in object space.Since the elevation interval was chosen to be 10 meters and the external accuracy of the buildings was not to be checked, there was no need of using the eATE tool of the software which is more time consuming.
The matching procedure in Classic ATE is pyramid based, that is, it starts from high-pyramid levels and continues to low pyramid levels.At each level the search range used stays the same on image space but is actually reduced to half on object space, which means search range goes smaller and smaller by brute-force.This method is proven to be a reliable solution (Xu et al., 2008) since the approach is based on multi-band image matching for increasing the precision of the results.Since this automatic DSM generation approach of the used software is mainly characterised by the feature-based matching technique being hierarchically applied in image pyramids and a robust surface reconstruction with finite elements, it is expected the resulting DSMs from different pansharpened data to be differentiated between them.Remarkable is that all the DSMs resulting from the 4-Band Pansharpened data are displayed by almost the same grey values.However, there values are differentiated by the values consisting the DSMs resulting from the 3-Band data.This is probably occurs due to the existence of the 4 th Band, the NIR band.

SPATIAL ACCURACY ASSESSMENT OF PANSHARPENING ALGORITHMS
In order to assess and compare the pansharpening algorithms in terms of accuracy, the same 25 checkpoints used in the previous work (Agrafiotis and Georgopoulos, 2015) were measured on the resulting orthoimagery.In more detail, the 25 checkpoints were measured on the orthoimagery created using the IHS, PCA and BROVEY pansharpening algorithms for three (RGB) and four spectral bands (RGB-NIR).The coordinates of these checkpoints were determined (i) on the reference/control data i.e., the aerial orthophoto (VLSO) provided by the Hellenic National Cadastre & Mapping Agency S.A. (NCMA S.A), with 25cm GSD, acquired during August 2012 and also (ii) on the orthoimagery created using the Panchromatic data.

Checkpoints Selection and Distribution
According to the NSSDA and JRC Guidelines, accuracy testing should be performed using an independent source of higher positional accuracy.The accuracy of the independent test points should fall within one-third of the intended accuracy (95% confidence level) of the examined dataset.A minimum of 20 well defined test points should be used to evaluate the accuracy of the dataset.The check points will be ideally evenly distributed and located across the image (Figure 4).The selected check point positions may be located with reference to the positions of the GCPs used to correct the imagery in order to ensure that the two sets of points are independent (CPs should not be close to the GCPs).The location or the distribution of the checkpoints is also specified in NSSDA and JRC guidelines.These Standards assume that the area to be evaluated is a rectangle and is divided into four quads and a diagonal is to be established across the area.At least 20% of the points should lie in each quarter whereas the optimum distance between points (is related to the diagonal distance of the area (1/10th of the diagonal length).
The results are presented with a statistical analysis and they are evaluated in order to present the effect of the pansharpened imagery on the orthoimage production process.To this direction the Standard Deviation (σ) or sigma are computed as an indicator of how well the measurements fit each other and a measure of precision.In addition, the Root Mean Squares Error (RMSE) is computed for Northing (X) and Easting (Y) coordinates as well as the combined RMSE (XY).It is assumed that errors in the spatial data have random behavior and that systematic errors have been eliminated as best as possible. (2) (3) where d = the deviation n = the number of check points X,Ycheck = the check points coordinates measured on control dataset and X,Ydata = the points coordinates measured on test dataset

Overall Evaluation
Tables 1 and 2 present the resulting residuals for all performed tests for the Pleiades 1B dataset.More specifically, the residuals of the comparison of the control dataset (VLSO) with (a) the Panchromatic orthoimage (Table 1), (b) the orthoimages resulting from four band pansharpened data (Table 1) and (c) three band pansharpened data (Table 2) are presented.In addition, results of the comparison of the panchromatic orthoimage with the orthoimages from pansharpened data of three and four bands are included.

RESIDUALS (m)
However, the reported residuals are within the accuracy limits that are specified by the JRC Guidelines and NSSDA Accuracy Standards.Moreover, it is observed that orthoimageries resulting from four band data have almost the same residual errors as the orthoimageries produced from three band data.
For the comparison between the Panchromatic orthoimage and the orthoimages produced using Pansharpened data (Table 2), it is obvious that the three band orthoimagery presents better accuracy 25%.This was initially not anticipated.It was expected that the 4 bands would provide more data, thus enhancing the matching results, which sees that it is not the case.However, the residual errors of both datasets are about the size of the GSD and within the accuracy limits for producing LSO.Among the pansharpening algorithms, it appears that the Brovey fusion presents the largest residuals at 75% of the cases (3 out of 4) while the PCA is covering the rest 25%.The IHS pansharpening algorithm residual errors seem to be the most times similar with the average of the residual errors in most of the cases.The above results suggest that the Pansharpening process doesn't seem to affect the geometric accuracy of the resulting DSM with a 10m interval and the resulting orthoimagery.Although some residuals in the orthoimages produced by the methods examined are observed, their magnitude cannot adversely affect the overall accuracy of the orthoimagery.For accurate orthoimages, DSMs could be generated either by panchromatic data or by pansharpened data as regards the tested algorithms without compromising the geometric accuracy of the produced orthoimage.This leads the users to choose the appropriate fusion algorithm based only on its radiometric properties and capabilities.The above observations may lead users to acquire the pansharpened products directly and thus saving processing time.

CONCLUDING REMARKS AND FUTURE WORK
In this paper the effects of the IHS, PCA and Brovey pansharpening algorithms on resulting orthoimagery were quantitatively evaluated.In particular, an iterative methodology was applied in order to combine the methodology described in prior work (Agrafiotis and Georgopoulos, 2015).DSMs and Orthoimages were created using: (a) only the Panchromatic and (b) the Pansharpened data for each block.In this way, for each different block, the same control points, check points and tie points were used.
Future work includes testing more state of the art pansharpening algorithms on more Very High Resolution Satellite imagery such as GeoEye-1 and WorldView-2.In addition, these tests will be applied using various software and freeware in order to acquire more reliable and objective results for this investigation.

Figure 1 .
Figure 1.Stereo acquisition for 3D applications (Google Earth preview of the footprints and the satellite's position)

Figure 2 :
Figure 2: Workflow chart for checking the effects of pansharpening process recent years, pansharpening has become a task of great importance in the field of data fusion, as demonstrated by the increasing number of scientific contributions to this topic.It is possible to group them into two main families (Aiazzi et al., 2012): i.The methods based on the projection of the MS image into a new space and the substitution of a component with a histogram matched version of the PAN image (the socalled Component Substitution (CS) class) and ii.The approaches based on the extraction of spatial details from the PAN image and their injection into the MS one (this class is called Multi-Resolution Analysis (MRA) because details are usually obtained through a multiscale decomposition of the original image).There are also hybrid methods based on both MRA and CS.Then there are methods based on varied techniques, such as the P+XS method.Widely known examples of the CS methods are Intensity Hue Saturation (Haydn et al., 1982), Principal Component Figure 3: A zoom of the original MS image (a) and the Pansharpening results from IHS (b), Brovey (c) and PCA (d) algorithms for 4 Bands.

Figure 4 .
Figure 4. Checkpoints distribution on the reference/control dataset

Figure 5 :
Figure 5: Graphic representation of the resulting residuals in meters

Table 1 :
Resulting residuals for the performed experiments on Pleiades 1B when compared with the control dataset As it can be observed in Table

Table 2 :
Resulting residuals for all experiments on Pleiades 1B panchromatic data compared with the pansharpened data.