Validation of a Simplified Model to Generate Multispectral Synthetic Images

Abstract

A new procedure to assess the quality of topographic correction (TOC) algorithms applied to remote sensing imagery was previously proposed by the authors. This procedure was based on a model that simulated synthetic scenes, representing the radiance an optical sensor would receive from an area under some specific conditions. TOC algorithms were then applied to synthetic scenes and the resulting corrected scenes were compared with a horizontal synthetic scene free of topographic effect. This comparison enabled an objective and quantitative evaluation of TOC algorithms. This approach showed promising results but had some shortcomings that are addressed herein. First, the model, originally built to simulate only broadband panchromatic scenes, is extended to multispectral scenes in the visible, near infrared (NIR), and short wave infrared (SWIR) bands. Next, the model is validated by comparing synthetic scenes with four Satellite pour l'Observation de la Terre 5 (SPOT5) real scenes acquired on different dates and different test areas along the Pyrenees mountain range (Spain). The results obtained show a successful simulation of all the spectral bands. Therefore, the model is deemed accurate enough for its purpose of evaluating TOC algorithms.

1. Introduction

The use of remotely sensed data from mountainous regions generally requires additional pre-processing, including topographic correction (TOC). Specifically, variations in the solar incidence angle (γ_i) affect land cover discrimination, since the radiance observed for a given land cover varies depending on whether it is located on shadowed or non-shadowed areas [1]. This effect, normally referred to as the topographic effect, can adversely affect the usefulness of remote sensing data for different applications, such as Land-Use/Land-Cover mapping, vegetation cover monitoring, change detection or biophysical parameter estimation [2,3,4,5,6]. The objective of topographic correction algorithms is to compensate the differences in radiance between sunny and shaded areas caused by variations in the slope and aspect of terrain.

A number of TOC algorithms were proposed in the past (e.g., [1,4,7,8]), but their evaluation is not as simple as it might seem. In this sense, several strategies to evaluate TOC algorithms have been developed, i.e., visual assessment of the removal of the topographic effect in satellite imagery [5,9,10,11], quantification of the reduction of the dependence between cosγ_i and the radiance of each spectral band after the correction [12], analysis of the variations in the radiometry of the corrected scenes [5], measurement of the reduction of land cover class variability [2,13,14], and improvement on classification accuracy after topographic correction [4,11]. However, these procedures were not purely objective [15,16,17]. Therefore, Sola et al. [16] proposed a new methodology to quantitatively evaluate topographic correction algorithms based on synthetic imagery. In short, the approach proposed by Sola et al. was based on the generation of a pair of synthetic images a sensor would acquire for any given area, considering, on the one hand, its real topography and, on the other hand, a completely flat surface. Then, the latter could be used as a reference to compare against the TOC corrected synthetic scenes, using quantitative indexes, in a rigorous and consistent manner.

The approach proposed in [16] presents several advantages compared to traditional evaluation techniques, such as being objective, simple, and not requiring ancillary information on land-cover distribution to perform the TOC quality assessment. However, the model proposed for generating synthetic images was developed to simulate only broadband panchromatic scenes. This could imply an important limitation for applying this methodology to evaluate TOC algorithms over a range of sensor (e.g., spatial resolution, band frequency, etc.) and acquisition (e.g., solar geometry and acquisition time) configurations. In addition, the whole approach was based on the assumption of model validity, and this needs to be verified.

For that purpose, in this work, the model to generate synthetic images is extended to multispectral scenes, i.e., visible and infrared bands, to adapt the approach for a range of sensors. The model is then validated, since this is something necessary in any modeling attempt, by comparing simulated multispectral scenes with real SPOT5 imagery acquired in four different study areas and four different dates.

2. Study Area

Four different study areas of 15 × 15 km were selected, all of them located in Northern Spain, in regions of rough relief in different parts of the Pyrenees. The dominant land-covers in these areas were deciduous and mixed forest, pastures, and agricultural crops, with a sparse presence of bare soil and urban areas. For this work a Digital Elevation Model (DEM) of 5-m resolution was used, obtained from LIDAR data acquired by the Spanish National Geographic Institute (IGN).

As seen in Table 1, different dates and solar geometry were selected to validate the model over different acquisition configurations. Furthermore, these corresponded to real SPOT5 acquisitions. The effect of the topography was expected to be more severe in area 1 due to its lower solar elevation angle.

Table 1. Details of the study areas used for the simulation.

Area	Center Lat/Long	Sensor Zenith	Date	Time	Sun Elevation	Sun Azimuth
1	43°08'35"N/1°42'54"W	12.6	15/10/2009	11:13	37.66	167.58
2	42°46'40"N/1°19'09"W	−24.3	15/08/2009	10:45	56.44	140.70
3	42°43'28"N/0°49'55"W	14.5	19/08/2009	11:08	57.97	152.54
4	43°06'06"N/2°06'33"W	15.0	30/08/2008	11:11	53.53	155.01

Table 1. Details of the study areas used for the simulation.

Area	Center Lat/Long	Sensor Zenith	Date	Time	Sun Elevation	Sun Azimuth
1	43°08'35"N/1°42'54"W	12.6	15/10/2009	11:13	37.66	167.58
2	42°46'40"N/1°19'09"W	−24.3	15/08/2009	10:45	56.44	140.70
3	42°43'28"N/0°49'55"W	14.5	19/08/2009	11:08	57.97	152.54
4	43°06'06"N/2°06'33"W	15.0	30/08/2008	11:11	53.53	155.01

3. Methods

3.1. Extension of the Model to Multispectral Images

The model to generate synthetic images was described in detail in [16]. In this work, the model was adapted to simulate multispectral SPOT5-like scenes with four bands, i.e., green, red, NIR, and SWIR bands. As already explained, a synthetic image represents the radiance the sensor would receive under certain conditions and within a spectral range. This spectral radiance (L_sen,λ) is obtained on the basis of the Lambertian reflectance law as follows [18]:

$$L_{sen,\lambda }=L_{p,\lambda }+\frac{{\rho }_{\lambda } T_{u,\lambda }{E}_{\beta ,g,\lambda }}{\pi }$$

(1)

where, L_p,λ is the path radiance of the corresponding spectral band, ρ_λ is land-cover spectral reflectance value, T_u,λ is the upward atmospheric transmittance, and E_β,g,λ is the global irradiance reaching each pixel. In [16], L_p,λ was calculated with Bird and Hulstrom’s model [19], and land-cover/land-use (LC/LU) cartography and spectral libraries were used to obtain ρ_λ. However, in this work L_p,λ and ρ_λ of each spectral band were obtained from real imagery to avoid introducing further uncertainties, since the focus is placed on validating the model, and in particular the simulation of the effect introduced by topography. In particular, SPOT5 images were converted from digital numbers (DN) to top of atmosphere radiance first, and to ground reflectance (ρ_λ) later, following the procedure explained in [20]. Moreover, the path radiance (L_p,λ) was extracted from the SPOT5 images themselves by using the Improved Dark Object Subtraction (IDOS) method [21].

The global irradiance of each band (E_β,g,λ) was obtained through the physical model of Sandmeier and Itten [22] as the sum of its three terms, i.e., direct, diffuse, and reflected irradiance:

$$E_{\beta ,g,\lambda }=\Theta \frac{E_{e,s,\lambda }\cos {\gamma }_i}{cos {\theta }_s}+E_{e,d,\lambda }\left[\Theta \frac{AI\cos {\gamma }_i}{cos{\theta }_s}+\left(1-\Theta AI\right)V_d\right]+E_{e,g,\lambda adj}{\rho }_{adj}{V}_t$$

(2)

where, Θ is the cast shadow’s binary factor, E_e,s,λ is the direct horizontal irradiance for each spectral band, γ_i is the solar incidence angle, θ_s is the solar zenith angle, E_e,d,λ is the diffuse horizontal irradiance, AI is Hay’s anisotropy index [23], V_d is the Sky View Factor [24], and V_t is the Terrain View Factor, that is, the portion of adjacent terrain seen from a certain location, ρ_adj is the average terrain reflectance reaching the adjacent slopes in a square box of 0.5 × 0.5 km assuming a Lambertian law for all the surrounding objects, E_{e,g,λ adj} is the average Global Horizontal Irradiance over a square box of the same size. The calculation of all these terms is explained in detail in Sola et al. [16].

To extend the original model to multispectral simulation some changes were introduced. First, both L_p,λ and ρ_λ were obtained for each spectral band. In addition, the direct and diffuse horizontal spectral irradiances were calculated for the considered bands, i.e., green, red, NIR, and SWIR. For that purpose, the fraction of direct and diffuse irradiance corresponding to each spectral band was calculated through the Simple Model of the Atmospheric Radiative Transfer of Sunshine (SMARTS2) spectral radiation model [25], considering mid-latitude summer/winter atmospheres, rural aerosol model and Thuillier solar spectrum. Finally, the obtained radiance values were converted from band-integrated values (W·m⁻²·sr⁻¹) to band-averaged values (W·m⁻²·sr⁻¹·μm⁻¹) dividing by the effective bandwidth of each band, obtained from the spectral response functions (SRF) of the High Resolution Geometric 2 (HRG2) sensor of SPOT5 [26].

3.2. SPOT5 Imagery

Four SPOT5 scenes were used for the validation. These were acquired under the same temporal and geometric conditions as the simulated scenes (Table 1). The sensor zenith angle ranges from −24° to 15°, (a negative incidence angle means the tilt direction is right of the flight direction), although these variations do not significantly affect the resulting synthetic scenes. The SPOT5 scenes, at a spatial resolution of 10 m, were orthorectified and converted from digital numbers (DN) to top of atmosphere radiance (W·m⁻²·sr⁻¹·μm⁻¹) by using the gain and offset provided in the metadata for each spectral band.

3.3. Validation

To validate the model, simulated scenes were compared with real SPOT5 scenes band per band. Three widely used statistical indexes were used to quantitatively evaluate the accuracy of the model: determination coefficient (R²) that measures the correlation between real and simulated spectral bands; mean structural similarity index (MSSIM) that measures their structural similarity [27]; and RMSE their root mean squared error. Additionally, scatterplots and histograms of both simulated and observed radiances were plotted to evaluate the quality of the simulation.

4. Results and Discussion

Four different synthetic scenes were generated (Figure 1, and supplementary results) and compared to their corresponding SPOT5 images. Only results of area 1 were shown here, while the rest were included as supplementary data.

Figure 1. Area 1 (a) RGB false color composite of the real scene; (b) RGB false color composite of the simulated scene; (c) RGB false color composite of terrain reflectance; (d) Cosine of solar incidence angle.

Figure 1. Area 1 (a) RGB false color composite of the real scene; (b) RGB false color composite of the simulated scene; (c) RGB false color composite of terrain reflectance; (d) Cosine of solar incidence angle.

Figure 1c showed the false color composite of ground reflectance used in the simulated scene. In Figure 1d the illumination, i.e., cosine of the solar incidence angle, was displayed. Both were used in the model to generate synthetic images. Visually, the simulated false color composites were very close to their corresponding real scenes (Figure 1a and 1b and supplementary results), but the former showed more spatial detail introduced by the 5-m DEM, while the real scenes looked slightly smoother. In this work the DEM was resampled to 10 m, without smoothing it. Although this issue is considered to be minor for the purpose of this work, further research is needed on the effect of spatial resolution and DEM smoothing on TOC performance, as discussed by several authors in the last years [13,28,29].

The shadowed areas introduced by the topography were well modeled, especially in areas 1 and 4. A more detailed analysis using scatterplots and histogram comparison (Figure 2 and supplementary results) confirmed this observation, with scatterplots following closely the 1:1 line and histograms of very similar shapes for the observed and simulated scenes.

Figure 2. Results for area 1. Scatterplots and histograms of observed and simulated radiance (in units of W·m⁻²·sr⁻¹·μm⁻¹) of (a) Band 1, (b) Band 2, (c) Band 3, (d) Band 4.

Figure 2. Results for area 1. Scatterplots and histograms of observed and simulated radiance (in units of W·m⁻²·sr⁻¹·μm⁻¹) of (a) Band 1, (b) Band 2, (c) Band 3, (d) Band 4.

Some limitations were observed though. In some areas, mainly area 2 and area 3, infrared bands, and, in particular, SWIR, seemed to introduce some more topographic effect than observed in real scenes, slightly underestimating low values of radiance and overestimating high radiances, i.e., slopes facing the sun. This effect had an influence in the results of statistical indexes, with lower values of MSSIM for the 4th band. This could be due to a higher influence of the direct irradiance term on the global irradiance impinging the surface for these bands, since this term is strongly influenced by cosγ_i, incrementing the variance in at-sensor radiances. However, this effect was not consistent in all the test sites, and, thus, in most cases, the topographic effect seemed to be well modeled.

In area 1, pixels with low radiance were slightly overestimated in all bands, so the false color composite in shadowed areas looked slightly darker in the real scene (Figure 1a) than in the simulated one (Figure 1b). This was also clearly visible in scatterplots (Figure 2) with values above the 1:1 line for low radiances and also in histograms, with simulated radiances slightly skewed to the right, but this effect was not apparent in areas 2, 3, and 4 (supplementary results).

In area 2 histogram shapes were well reproduced, showing a bimodal distribution due to the presence of agricultural crops on the one hand and forest areas on the other hand. Scatterplots showed a good coincidence between simulated and observed radiances. However, in band 3 and 4 the variance of the radiance was higher in the simulated scene. This effect was visible in Figure S2d, where the small peak at the lowest region of the histogram corresponding to a reservoir was underestimated, while the high radiances corresponding to agricultural crops in the left bottom corner of the area were overestimated.

In addition, in area 3 and 4 high radiances of band 1 and 2 were underestimated, probably corresponding to urban areas, with more complex reflective behavior and mainly located in the valleys, so therefore not affected by topography. Those were only a few pixels and they do not affect the quality of the simulation, as can be assessed both visually, and statistically (Table 2).

It was noticeable the presence of a clear high bound in the observed radiances in the scatterplots, which was not present in the synthetic images. This was caused by the format of the original imagery, stored at 8 bits. When the DN were transformed to radiance an upper limit was set, but this was not occurring in the synthetic images, where areas of high ground-reflectance in slopes facing the sun ended in higher values of at-sensor radiance. In addition, histograms of observed radiances seemed to be serrated, which is a typical effect when a smaller integer color space is expanded to a larger one.

Table 2. Statistical indexes to measure similarity between real and simulated scenes.

AREA	B1			B2			B3			B4
	R2	RMSE	MSSIM	R2	RMSE	MSSIM	R2	RMSE	MSSIM	R2	RMSE	MSSIM
1	0.92	1.38 (6.1%)	0.848	0.92	1.38 (10.6%)	0.802	0.87	5.93 (11.5%)	0.787	0.86	0.95 (14.1%)	0.633
2	0.96	2.85 (5.1%)	0.840	0.95	4.07 (9.1%)	0.806	0.86	4.10 (6.2%)	0.807	0.93	1.41 (9.9%)	0.681
3	0.92	2.51 (5.3%)	0.876	0.90	2.25 (9.9%)	0.803	0.97	2.15 (3.9%)	0.966	0.87	1.14 (14.3%)	0.701
4	0.99	1.09 (3.2%)	0.977	0.99	1.19 (5.8%)	0.966	0.98	3.75 (5.8%)	0.961	0.97	0.65 (7.5%)	0.911

Table 2. Statistical indexes to measure similarity between real and simulated scenes.

AREA	B1			B2			B3			B4
	R2	RMSE	MSSIM	R2	RMSE	MSSIM	R2	RMSE	MSSIM	R2	RMSE	MSSIM
1	0.92	1.38 (6.1%)	0.848	0.92	1.38 (10.6%)	0.802	0.87	5.93 (11.5%)	0.787	0.86	0.95 (14.1%)	0.633
2	0.96	2.85 (5.1%)	0.840	0.95	4.07 (9.1%)	0.806	0.86	4.10 (6.2%)	0.807	0.93	1.41 (9.9%)	0.681
3	0.92	2.51 (5.3%)	0.876	0.90	2.25 (9.9%)	0.803	0.97	2.15 (3.9%)	0.966	0.87	1.14 (14.3%)	0.701
4	0.99	1.09 (3.2%)	0.977	0.99	1.19 (5.8%)	0.966	0.98	3.75 (5.8%)	0.961	0.97	0.65 (7.5%)	0.911

As seen in Table 2, all the three indexes showed good results for the first three spectral bands, but a poorer performance of band 4, in particular in areas 1 and 3. This could be due to the lower influence of atmospheric scattering in this band. In all the cases the statistics were in line with results of other simulation models [30,31]. The coefficient of determination (R²) kept above 0.85, and MSSIM above 0.60. The RMSE between observed and simulated spectral radiances ranged from 1 to 6 W·m⁻²·sr⁻¹·μm⁻¹ (3% to 14% of the mean value), very close to the results obtained by Verhoef et al. [32] for similar spectral bands. Among the study areas, all of them performed good but the first seemed to perform slightly worse than others according to the statistical indexes, probably due to the more severe illumination conditions derived from the lower solar elevation angle.

5. Conclusions

The aim of this letter was to extend the synthetic image simulation model proposed in [16] to the multispectral case and to validate this model using real SPOT5 imagery. The results obtained using four test sites with different acquisition conditions illustrate a good behavior of the model. The comparison between simulated and real SPOT5 scenes yielded R² values above 0.90 for visible bands, and above 0.86 for the NIR and SWIR. Similarly, using MSSIM values could be ranked in order of accuracy as follows: green (>0.84), red (>0.80), NIR (>0.78), and SWIR (>0.63). These results were consistent for the four different test areas, although there were differences between them, with study area 1 achieving the lowest accuracies. This could be partly explained by the lower solar elevation angle of this area. All in all, relative RMSE values were normally below 10% of the observed radiance, which is considered accurate enough for the purpose this model was designed for. Thus, the TOC evaluation approach proposed in [16] could be subsequently used with multispectral data for evaluating TOC algorithms on different areas and acquisition conditions.