Fusion of Radarsat-2 polarimetric images for improved stereo-radargrammetric DEM

Abstract

The new full polarimetric modes of Radarsat-2 were evaluated for digital elevation model (DEM) generation using stereo-radargrammetry with Toutin's three-dimensional physical model. The fusion of the four polarimetric bands was thus addressed to improve the resulting DEMs. The stereo-radargrammetric DEMs were evaluated with 20 cm accuracy Lidar elevation data. Results on a test site north of Québec City, Canada showed a good accuracy: 5 m and 6–7 m (1σ) for the horizontal positioning and the elevation (on bare soils), respectively. However, the best results were obtained when the four polarimetric bands were fused with the restoration method and the second best results (around 5% worse) were obtained with the SPAN method (total power of the four intensity bands).

Keywords:

• Radarsat-2
• image fusion
• polarimetry
• radargrammetry
• stereo
• DEM

1. Introduction

Since the 1960s, stereo-radargrammetry methods have been used to extract elevation information from synthetic aperture radar (SAR) data (La Prade 1963). The research studies were, however, limited due to the lack of available satellite stereo-data. Radarsat-1 SAR with its various operating modes, beam positions and resolutions (Parashar et al. 1993) renewed the interest in stereo-radargrammetry for generating digital elevation models (DEMs; Sylvander et al. 1997, Toutin 1999, 2000). There was a general consensus on the achieved DEM extraction accuracy using Radarsat-1 data: a little more (12 m) than the sensor resolution with the fine mode and a little less (20 m) than the sensor resolution with the standard mode. These results were independent of the method used (digital stereo plotter or automatic image matching; Toutin and Gray 2000).

To obtain good geometry for stereo plotting, the intersection angle should be large in order to increase the stereo exaggeration factor or equivalently the observed parallax, which is used to determine the terrain elevation. Conversely, to have good stereo viewing, the interpreters (or the image matching software) prefer the images to be nearly identical as possible, implying a small intersection angle. Large geometric and radiometric disparities together hinder stereo viewing and precise stereo plotting and thus a compromise has to be reached between better stereo viewing (small radiometric differences) and more accurate elevation determination (large parallax), and the common compromise for any type of relief is to use a same-side stereo pair, thus fulfilling both the above conditions (Figure 1).

Figure 1. Same-side stereo-radargrammetry with RADARSAT-2 multi-polarisation data acquired with slant-range fine-quad mode (FQ5 in B and FQ18 in A). Note: ‘RADARSAT-2 Data © MacDONALD, DETTWILER AND ASSOCIATES LTD. (2008)–All Rights Reserved’ and courtesy of Canadian Space Agency (CSA).

With the launch of Radarsat-2 in December 2007, a new SAR stereo-data can be acquired in high-resolution and polarimetric modes (Figure 2). With the polarimetric modes, four polarimetric amplitude/intensity images (HH, VV, VH and HV) can be generated for each image of a stereo-pair (Figure 1), and can be used for generating four different DEMs over the same area. On the other hand, it is also possible to fuse the four polarimetric images into a single image, which preserves the ‘best information’ for the DEM generation process, e.g. mainly during the image matching for elevation parallax extraction. This article will thus investigate different known algorithms for polarimetric image fusion and evaluate their impact on DEM positioning and elevation accuracy generated from Radarsat-2 fine quad-polarisation (FQ) modes. The stereo-radargrammetry process uses Toutin's three-dimensional (3D) radargrammetric model, previously developed at the Canada Centre for Remote Sensing (CCRS) for Radarsat-1 (Toutin 1999, 2000) with a preliminary version for Radarsat-2 (Toutin and Chénier 2009)1.

Figure 2. Operating modes of Radarsat-2 C-band SAR. Courtesy of Canadian Space Agency, 2010.

2. Fusion of polarimetric images

There are many ways to represent or decompose full polarimetric data. The goal of our polarimetric image fusion is to generate a composite image to preserve information useful for the DEM generation process, mainly during the image matching. It is thus important to preserve image intensity, contrast and texture while removing the speckle (Zakharov and Toutin 2009a). Different algorithms (total power of four polarimetric bands, the so-called SPAN, median/Frost/enhanced-Lee filtering with various window sizes, polarimetric whitening filtering (PWF), principal component (PC) analysis and restoration) were evaluated using a variety of image quality measures (equivalent number of look (ENL), texture, mutual information and auto- and cross-correlation coefficients). Based on these image quality measures, three fusion algorithms (SPAN, PC analysis and restoration) were considered to be the most appropriate for image matching (Zakharov and Toutin 2009a). Since the PC analysis is well known and the evaluations of these fused methods were already addressed in Zakharov and Toutin (2009a, b), only the SPAN and Restored fusions are described.

2.1 SPAN

The simplest fusion algorithm is SPAN, which corresponds to the sum of the four polarimetric intensity (power) images. SPAN image contains less speckle than each polarimetric image, because (1) the speckle is not correlated between the polarimetric images and (2) each polarimetric band contains different information and can be considered as an additional measurement such as multilook SAR imagery. The mutual information calculated between the different polarimetric channels and SPAN showed that HH and VV bands have a higher contribution than HV and VH bands in the SPAN image (Zakharov and Toutin 2009b).

2.2 Restoration

2.2.1 Theoretical aspect

The linear process of several images formation can be described by the following integral equation:

where (i, j) are the coordinates in image space, (x, y) are the coordinates in the object space, p is the number of registered intensity image, is the intensity image represented digitised signal after raw signal processing, Z(x, y) is the object (original image), which includes speckle noise, is the kernel which describes Point Spread Function (PSF), and is the additive noise which has the nature of the sensor noise and quantisation.

Theoretically, the decomposition of the object Z(x, y) is based on the regularised image restoration method (Dovnar and Predko 1986):

where are orthonormalised functions. The stabilised vector parameter minimises the mean square error (MSE) when the object Z(x, y) and the additive noise are not correlated. The decomposition coefficients are described by a recurrent relation (an equation, which defined a sequence recursively). This method requires a priori knowledge of the PSF, the power spectra of the object and the noise. They can be estimated similar to the Wiener method where the object and the additive noise are characterised by autocorrelation functions (Lim 1990). However, one advantage of the image restoration method is the capability to compute the MSE over a small image area using a reduced number of pixels (Dovnar and Predko 1986). Conversely to the Wiener method, the knowledge of eigenfunctions of the kernel of Equation (1) is not required and the method is still optimal for any discretisation, which is important when resolution reaches the Nyquist limit. Furthermore, the restoration method can also be used as an image filtering process with a filter Q:

Finally with the image restoration method, the information-theoretical assessment enables a potential improvement of the sensor resolution during the image fusion process (Fellgett and Linfoot 1955, Dovnar et al. 2002).

2.2.2 Applicability to SAR polarimetric data

The image restoration method can be applied for processing SAR data because it is a linear imaging system. However, the main problem for its direct application to coherent systems, such as SAR, is the multiplicative speckle noise. On the other hand, the additive noise in Equation (1) was estimated with an SAR image analysis and its dispersion was found to be equal to half the quantisation level for an 8-bit image (8 bits is a requirement in the matching process of DEM generation). The SAR PSF consists of the range and azimuth PSFs, which can be computed from the analysis of experimental data of a point source image (e.g. using a corner reflector). For different polarimetric bands of the same beam the PSF has identical shape. The computation of the main lobe of PSF's magnitude can be approximated by the sinc function (Cumming and Wong 2005). The information-theoretical assessment of the restoration method for processing Radarsat-2 SAR polarimetric images showed that the amount of information increased for higher spatial frequencies (Zakharov and Toutin 2009b).

Generally, due to speckle noise, the model (Equation (1)) for SAR data processing is simplified for despeckling and enhancement purposes. Thus the object Z(x, y) in Equation (1) is considered to be equal to the acquired intensity while the PSF and the additive noise are neglected. After simplification, the model for SAR image formation can be written as

where is the SAR backscatter and is the multiplicative speckle noise.

When the noise is considered as signal-independent, its parameter in Equation (4) can be logarithmically transformed into additive noise (Bovik 2000). The logarithmic transformation can be used for speckle reduction and enhancement processes of SAR imagery (Vidal-Pantaleoni and Martí 2004, Gupta and Gupta 2007).

The restoration was applied as a filter (Equation (3)) with an adaptation of the Wiener method to SAR image processing (Pascazio and Schirinzi 1996) using logarithmic scale: (1) to transform the multiplicative speckle noise into additive and (2) to compress the initial Radarsat-2 image from 16 bits to 8 bits. This scale reduction is required by some software limitation in the stereo-image matching.

The restoration algorithm (Equation (3)) reduces the speckle without degrading the original spatial resolution of each polarimetric image. The selection of the stabilised vector parameter in Equation (3) is used to enhance small features and details in the filtered image (Zakharov and Toutin 2009b). Finally, the fusion of the four polarimetric images improves the spatial resolution and enhances the fused image because the small details of each polarimetric image are preserved in the output.

3. Study site and data sets

The study site is located north of Québec City, Québec, Canada (47°N, 71°30′W) and spans different environments: urban and residential, semi-rural and forested (Figure 3). Québec City, St. Lawrence River and Ile d’Orléans are on the bottom right part of Figure 3. The elevation ranges almost from 0 m at the St. Lawrence River at the southeast, to around 1000 m in the Canadian Shield, located to the north. It is a Precambrian geomorphology with the oldest smoothed-glacial (eroded) topography. The northern part is mainly covered by forests (deciduous, conifer and mixed) with a hilly to mountainous topography (5–30° slopes) while the southern part is urban and residential areas with a semi-flat topography (0–5° slopes).

Figure 3. R(HH)/G(VH)/B(VV) composite with slant-range fine-quad image (FQ5): Québec City, St. Lawrence river with Île d’Orléans are on the bottom right. Note: Cities appear in reddish due to more double bounce scattering with HH and VV than VH polarisations. The Lidar elevation data coverage (5 × 13 km) is overlaid (yellow rectangle). ‘RADARSAT-2 Data © MacDONALD, DETTWILER AND ASSOCIATES LTD. (2008)–All Rights Reserved’ and courtesy of Canadian Space Agency.

The Radarsat-2 data set includes (Table 1) a stereo-pair formed with two images (25 by 25 km) FQ5 (Figure 3) and FQ18 (Figure 4) acquired with the C-band fine quad-polarisation mode (5.4 by 5–8 m resolution; 1 by 1 look). FQ5 and FQ18 were acquired between 23 June and 26 June 2008 from ascending orbits with view angles of 23.4–25.3° and 37.4–38.9° at the near-far edges, respectively. They only had 50% overlap with 13° intersection angle and a mean vertical parallax ratio (VPR) of 0.87 (Toutin 1999). The images were processed as slant-range single-look complex (SLC) product (with amplitude and phase data; 4.7 by 5.1 m pixel spacing). Cities (Figures 3 and 4) appear in red tone due to more double bounce scattering with HH than VV and VH polarisations.

Figure 4. R(HH)/G(VH)/B(VV) composite with slant-range fine-quad image (FQ18). Note: Cities appear in reddish due to more double bounce scattering with HH and VV than VH polarisations. ‘RADARSAT-2 Data © MacDONALD, DETTWILER AND ASSOCIATES LTD. (2008)–All Rights Reserved’ and courtesy of Canadian Space Agency.

Table 1. Description of Radarsat-2 single look complex fine-quad stereo-images.

Mode and beam	Acquisition date	Orbit	Look angles	Ground coverage (km)	Sensor resolution (m)	Pixel spacing (m)
FQ5	23/06/2008	Ascending	23.4–25.3°	25 × 25	5.4 × 8.0	4.7 × 5.1
FQ18	26/06/2008	Ascending	37.4–38.9°	25 × 25	5.4 × 8.0	4.7 × 5.1

The cartographic data, obtained from the Ministère des Ressources naturelles du Québec, Canada were the 1 : 20,000 raster maps with 10 m contour lines and 1 m ortho-photos. Ground control points (GCPs) were collected from the raster maps with 10 m contour lines and ortho-photos: their planimetric and elevation accuracy are 3–5 m. They are generally road intersection in the low lands and urban areas and lake-river borders in the mountains. Theoretically seven sub-pixel accurate GCPs are enough to compute the stereo-radargrammetric model, 2but more GCPs were acquired either to reduce the propagation of errors (planimetric/elevation and image pointing) in the least-squares adjustment (due to an overestimation) or to perform accuracy tests with independent check points (ICPs).

To evaluate the accuracy of the stereo extracted DEM, accurate spot elevation data was obtained from a Lidar survey by GPR Consultants, Quebec, Canada (http://www.lasermap.com) on 6 September 2001 (yellow square in Figure 3). The Optech ALTM-1020 system is comprised of a high-frequency optical laser coupled with a Global Positioning System (GPS) and an Inertial Navigation System (INS; http://www.optech.on.ca). A 3D GPS solution is used to position the laser scanner at each second or half second, while the INS data are used to determine the systems’ orientation. The GPS solution is computed from differential kinematic processing, using data collected simultaneously at the aircraft, and at base stations near the project site (Fowler 2001). From a fixed-wing airborne platform, the laser emits pulses at frequencies of up to 5000 Hz and the first echoed pulses are reflected off vegetation or human-made structures and recorded. Consequently, in most of the forest areas, the elevation approximately corresponds to the top of the canopy or a little lower. With 700–850 m flying height, 70 m s⁻¹ velocity, 5000 Hz pulse rate, 12 Hz scanning frequency and ±20° scan angle (510–630 m wide swath), the ground point density was about 300,000 3D points per minute and the accuracy was 0.30 m in planimetry and 0.15 m in elevation. Ten swaths were acquired covering an area of 5 km by 13 km, which nearly corresponds to the western part of FQ5. The results are then an irregular-spacing grid (around 3 m), due also to no echo return in some conditions (buildings with black roof, some roads and lakes). Since the objectives were to evaluate the stereo FQ5–18 DEM, the Lidar elevation data were not interpolated into a regular spacing grid to avoid the propagation of interpolation error into the checked elevation and the final accuracy evaluation.

4. Generation of stereo-radargrammetric DEMs

The Toutin's 3D radargrammetric model used for the DEM generation is based on the integrated and unified 3D physical model for multi-sensors developed at CCRS (Toutin 1995), adapted for Radarsat-1 (Toutin 1999, 2000) and Radarsat-2 (Toutin and Chénier 2009). The precision of the mathematical modelling for SAR is better than 25 cm (Toutin and Chénier 2009). Since the processing steps of DEM generation from Radarsat-2 stereo images are roughly the same as for other stereo images (data collection and pre-processing, stereo bundle adjustment with GCPs, elevation parallax measurements, DEM generation and editing), the six processing steps (Figure 5) are summarised below (Toutin 1995):

Figure 5. Flow chart of the DEM generation process.

Before the DEM evaluation, there are six main processing steps.

1. Acquisition and pre-processing of the SLC SAR data (images and metadata) to determine an approximate value for each parameter of 3D physical model for the two images. The SLC data were pre-processed into amplitude images, filtered with the enhanced Lee filter. The non-filtered FQ images are fused using the three algorithms previously mentioned (SPAN being the total power of four bands, principal components and Restoration).

2. Collection on HH images of 90 ground points with their 3D cartographic coordinates and two-dimensional (2D) image coordinates. These ground points were used for all images (single pol. and fused) and spanned the total surface with points at the lowest and highest elevation to avoid extrapolations, both in planimetry and elevation. However, only 28 GCPs are stereo points (yellow and pink in Figures 3 and 4) belonging to both images due to only 50% stereo-coverage. The remaining 62 points are monoscopic points (blue in Figures 3 and 4) belonging to only one image. The image pointing accuracy was generally one pixel (5 m) in the residential areas but sometimes two pixels in the rural areas (10 m). Most of the points will be used as Independent Check Points (ICPs) to perform accuracy tests.

3. Computation of the 3D stereo model, which is valid for all polarimetric and fused images, initialised with the approximate parameter values and refined by an iterative least-squares stereo bundle adjustment with the GCPs (step 2) and orbital constraints. Stereo-radargrammetric equations are used as observation equations and weighted as a function of input errors. The ground points are used either as GCPs to have an overestimation in the adjustment and to reduce the impact of errors or ICPs to perform accuracy evaluations.

4. Extraction of elevation parallaxes using multi-scale mean normalised cross-correlation method with sub-pixel computation of the maximum of the correlation coefficient.

5. Computation of 3D cartographic coordinates from elevation parallaxes (step 4) using the previously computed stereo model parameters (step 3) with 3D least-squares stereo intersection.

6. Generation of regular grid spacing with 3D automatic and visual editing tools: automatic for blunders removal and for interpolating the small mismatched areas, but visual for filling the large mismatched areas and for the lakes in operational environments (not performed in this research study).

Error propagation can be tracked along the processing steps with stereo bundle adjustment results (step 3) and during the DEM generation (steps 4 and 6). The final DEMs were thus compared to Lidar elevation data and linear errors with 68% level of confidence (LE68), minimum and maximum errors of elevation differences were computed over the Lidar-stereo overlap (around 5 by 13 km; yellow rectangle in Figure 3).

5. Analysis of results

5.1 Evaluation of stereo bundle adjustments

Because the polarimetric images and their different fusions have the same geometry, the stereo bundle adjustment results are valid for all. These first results enable the applicability of the 3D CCRS radargrammetric model to Radarsat-2 FQ stereo images to be confirmed (Toutin and Chénier 2009). As a function of numbers, location and accuracies of GCPs used in the stereo bundle adjustment, two tests were performed: Test 1 was conducted with all 90 GCPs: some are monoscopic points (blue points in Figures 3 and 4) and only 28 are stereo. Test 2 was performed with 10–8 stereo GCPs (yellow points in Figures 3 and 4) to assess the robustness and to maintain a sub-meter internal precision of the stereo model (Toutin and Chénier 2009) while the stereo-model was verified with the remaining 18–20 stereo ICPs (pink points in Figures 3 and 4) not used in the 3D stereo model calculation and which enabled an unbiased validation of the modelling accuracy to be performed. Table 2 summarises these results obtained from the residuals on GCPs (root mean square (rms) and minimum/maximum for Tests 1 and 2) and in addition the errors on ICPs (bias, rms and minimum/maximum for Test 2).

Table 2. Root mean square and absolute maximum residuals/errors on GCPs/ICPs, respectively and bias on ICPs.

Test number	GCPs ICPs	GCP rms residuals (m) X, Y, Z	Absolute GCP maximum residuals (m) X, Y, Z	ICP bias (m) X, Y, Z	ICP rms errors (m) X, Y, Z	Absolute ICPs maximum errors (m) X, Y, Z
Test 1	90/0	7.3, 8.1, 6.8	20.2, 17.8, 15.0	N/A	N/A	N/A
Test 2a	10/18	6.7, 2.2, 4.4	15.6, 3.3, 10.3	7.3, 6.5, 3.8	3.7, 5.5, 7.1	12.5, 21.0, 15.9
Test 2b	9/19	2.4, 1.6, 1.6	4.1, 3.3, 2.5	7.5, 6.2, 2.1	3.4, 5.4, 8.7	12.3, 20.9, 16.7
Test 2c	8/20	0.7, 1.8, 1.4	1.5, 3.4, 2.6	7.3, 6.3, 1.6	3.4, 5.5, 9.0	12.6, 22.0, 16.8

The rms residuals in Test 1 (7.3, 8.1 and 6.8 m in X, Y and Z, respectively corresponding to a little more than one pixel) were generally in the same order of magnitude as the input data errors (5–10 m). When much more GCPs than the minimum (7) required are used for computing the stereo-model with a least-squares adjustment, GCP residuals mainly reflect the GCP error budget (Toutin 1995, 2000): the largest being 1–2 pixel image-pointing error (5–10 m) and the smallest the cartographic error (3–5 m). Consequently, it is thus normal and ‘safe’ with a precise radargrammetric model to obtain rms residuals converging to the same order of magnitude as the GCP errors, here the image-pointing error, and more accurate cartographic coordinates (e.g. with GPS) would not change the results. However, the radargrammetric modelling is thus more accurate than the rms residuals because the GCP errors did not propagate into the radargrammetric model but into the residuals.

Test 2, using only 10–8 GCPs, confirmed these statements because the rms errors on ICPs (being an unbiased evaluation) are around one pixel (4–5 m) in planimetry and less than two ‘equivalent’ pixels (using the 0.87 VPR) in elevation. Because these rms errors include the image pointing or extraction error of ground features, they are thus a good estimation of the future 3D restitution accuracy for topographic mapping. Finally, it is normal that GCP rms residuals of Test 2c were very small (1 m) due to the use of only one extra GCP (8 vs. 7, the theoretical minimum) in the least-squares adjustment. Results with few GCPs, however, should never be considered as an accuracy indication of the modelling. On the other hand, because (1) there was a large-enough overestimation in the least-squares adjustment with 9–10 GCPs and (2) the rms residuals and errors on GCPs and ICPs, respectively, were about the same order of magnitude, the GCP rms residuals in Tests 2a and 2b can be used as a priori 3D mapping error in operational environments. The largest maximum error (around 20 m) on GCPs/ICPs was related to the same ground point (a less accurate feature), showing consistency and robustness of the 3D radargrammetric model, even with a reduced number of GCPs. Consequently, these overall results well confirmed the previous results (applicability and robustness) with ultra-fine mode Radarsat-2 data (Toutin and Chénier 2009).

5.2 Evaluation of DEMs

The second result is the qualitative and visual evaluation of the full DEM (10 m pixel spacing) and the quantitative evaluation over the Lidar coverage. Figure 6 shows only the four best over six stereo-radargrammetric DEMs (15 by 30 km; 10 by 5 m pixel spacing) extracted from HH, SPAN, PC1 and Restored FQ5–FQ18 images without any editing. The large black area at the bottom-right corner corresponds to the lowest elevation near the water level. All DEMs have the same global appearance and well reproduced the general terrain relief and the different cartographic and topographic features, which can be seen in FQ5–FQ18 images (Figures 3 and 4): the north–south largest and smallest valleys and the different surrounding mounts. Even small relief features in and between the valleys were captured. However, the following observations can be made:

• The HH DEM displays small non-continuous and/or sharp-angle topographic features (not representative of Precambrian topography), which do not appear on the three other ones and look like artefacts.

• The HH and SPAN DEMs appear more fuzzy than the PC1 and Restored DEMs.

• The SPAN and PC1 DEMs appear blockier than the Restored DEM.

• The SPAN and Restored DEMs look equivalent with more topographic details, which better correspond to the ground topography.

• The Restored DEM is smoother with more regular topographic feature shapes.

Figure 6. The four stereo-radargrammetric DEMs (15 by 30 km; 10 by 5 m pixel spacing) extracted from HH, SPAN, PC1 and Restored FQ5–FQ18 polarimetric images.

Figure 6 showed the good matching performance, regardless the stereo-pair; there was only few percents of mismatched areas where no elevation parallax was found during the correlation process and which were automatically interpolated. Because the mismatched areas were reduced for all stereo-pair, it is difficult to conclude on the image matching performance impact of the fused stereo-pair versus single-polarisation stereo-pair. More tests with different study sites and data set should be performed to look at this impact.

The quantitative evaluation of six DEMs was conducted with the comparison of the Lidar elevation data and around 1 million of DEM points were used in statistical computations. Table 3 gives these results computed from elevation comparison (Lidar minus DEM): the biases, the linear errors with 68% level of confidence (LE68) and the minimum and maximum errors (in metres). The first three lines are for the single polarisation stereo-pairs (HH, VV and VH) and the last three are for the fused stereo-pairs (SPAN, PC1 and Restored). There are no large error differences (bias, LE68 minimum/maximum) between all stereo-pair DEMs: it thus demonstrated that the statistical distribution of errors is quite equivalent. All biases (1–2 m) are positive but with no significant differences between the four methods (1.2 m being the largest). These positive biases are normal due to less penetration of the Lidar (1st echoed pulse reflected off vegetation) than the SAR C-band (double-bounced and volume scattering increasing the range) in the boreal forests, which represent around 70% of the study site. The minimum/maximum errors are consistent between the different methods but the differences are not significant. On the other hand, the 2 m maximum difference for LE68 is more significant to draw conclusions related to the fusion method for improved stereo DEMs.

Table 3. Errors results of DEMs: biases, errors with 68% level of confidence (LE68), minimum and maximum errors computed over around 1 million of DEM points.

DEM FQ5–FQ18	Bias (m)	LE68 (m)	Minimum errors (m)	Maximum errors (m)
HH	0.9	11.7	−77	102
VV	1.0	13.1	−90	113
VH	0.7	12.6	−96	108
SPAN	1.3	11.4	−83	110
PC1	1.9	12.1	−82	104
Restored	1.3	11.1	−81	113

According to higher textural, mutual information and correlation parameters (Zakharov and Toutin 2009b) combined with a higher signal-to-noise ratio of HH, it was expected a priori that HH achieved the best matching performance and elevation accuracy results between the three single polarisation stereo-pairs and these a posteriori results with 11.7 m LE68 versus 13.1 and 12.6 m for VV and VH, respectively, confirmed the a priori statements. On the other hand, VH with lower textural, mutual information and correlation parameters as well as signal-to-noise ratio should have achieved the worse results. However, less double-bounced scattering with vertical structures in VH stereo-images maybe offsets the impact of the previous parameters on image matching. It is also normal that PC1 method gave LE68 errors comprised between HH and VV LE68 results because PC1 is equally composed of these two components.

The image matching software used in this study is based on image correlation. When applying to fused images as a combination of image components and with weighting coefficients the correlation C can be written with the distributive property of convolution * as:

Because HH, HV, VH and VV are correlated (Boerner et al. 1998), the correlation C between all image components A_i and B_j in Equation (5) is always positive for SPAN (sum of power), except for specific conditions (e.g. water bodies). In addition, all weighting coefficients are equal to one (Equation (1)) for SPAN. Consequently, the SPAN correlation C is higher than the PC1 correlation C because weighting coefficients are applied in Equation (5) for PC1. In fact computing the two correlations C, the SPAN one was twice higher than the PC1 one: it thus explains the better performance and results with SPAN versus PC1 methods during the image matching for elevation extraction.

As expected (Zakharov and Toutin 2009a), the Restored fusion achieved the best results of all fusion methods due to its three advantages: resolution improvement, image enhancement and speckle filtering. In addition, because the equivalent number of look, the textural parameters as well as the mutual information between the two images were the highest parameters in the different fusion methods (Zakharov and Toutin 2009a), it was a priori indicator for obtaining the best DEM results.

However, as mentioned, these results do not reflect the true stereo DEMs accuracy since one source of errors comes from the fact that both sensors (SAR and Lidar) do not have the same penetration in the vegetated cover (few metre variations). Consequently, the same error evaluation was performed only on bare soils, where the stereo SAR and Lidar elevation points are the same. Only the three DEMs (HH, SPAN and Restored), which achieved the best results, were evaluated (Table 4). According to the biases, the stereo DEMs are consistently too high of 4–5 m. The minimum and maximum errors were drastically reduced (−11 and 2 m, respectively), confirming the largest errors in the total DEMs were due to some problems and errors during the matching process over forested areas. In an absolute point of view, LE68 over bare soils (around 6–7 m) are almost twice better than over the total DEMs, which are mainly offset by the land cover heights. These LE68 on bare soils thus correspond to the true elevation accuracy, which can be achieved with Radarsat-2 fine quad stereo data. On the other hand, the relative comparison of the three DEMs confirms the previous results and explanations of Table 3, e.g. the Restoration method, which fused/enhanced/filtered the four (4) polarimetric bands still achieved the best results.

Table 4. Errors results of DEMs on bare soils for the best methods (HH, SPAN, Restored): biases, errors with 68% level of confidence (LE68), minimum and maximum errors computed over around 1 million of DEM points.

DEM FQ5–FQ18	Biases (m)	LE68 (m)	Minimum errors (m)	Maximum errors (m)
HH	−4.9	6.7	−56	75
SPAN	−4.4	6.6	−52	69
Restored	−4.4	6.2	−57	74

While few-percent land-cover changes (mainly forests to residential but almost non-existent on bare soils) occurred between the two data sets (2001 for Lidar and 2008 for Radarsat-2), they were not statistically significant to impact the final results bias, LE68 and even minimum/maximum errors, especially when computed over bare soils. Only the tails of histograms could have been a little different. In addition, these changes did not impact at all on the relative comparisons of the six DEMs and the final conclusions.

6. Conclusions

Toutin's 3D radargrammetric model was adapted to Radarsat-2 and tested on fine-quad polarimetric images, FQ5 and 18, acquired over a hilly relief study site north of Quebec City. The general results confirm Toutin's 3D radargrammetric model accuracy and robustness to Radarsat-2 data, previously proven to be precise at sub-metre level (Toutin and Chénier 2009). Different fusion methods (SPAN, PCA and Restoration) of the polarimetric images were thus introduced into the stereo-radargrammetric process for DEM generation. The qualitative and quantitative evaluations of six DEMs from HH, VV, VH, SPAN, PC1 and Restoration stereo-pairs demonstrated a slight improvement with the SPAN and Restoration methods when compared to the single polarisation (HH). However, the best results (6.2 m LE68 on bare soils) were consistently obtained with the Restoration method, in accordance with the results of previous evaluations using different image quality measures (Zakharov and Toutin 2009a,b). More tests will be performed with other data sets (study site, SAR polarimetric FQ data, looks and stereo angles) to confirm these results.

Acknowledgements

This research work was financed by the Science and Operational Applications Research program between the CSA and MDA, by the Government Related Initiative Program of CSA, and supported by the Natural Sciences and Engineering Research Council of Canada (NSERC). The ESS contribution number is: 20090339.

Notes

Notes

1. An advanced version of 3D Toutin's radargrammetric model has been recently developed to take advantage of the more precise orbit determination of Radarsat-2 since September 2008.

2. With the advanced version of 3D Toutin's radargrammetric model, no GCP is required now.