High Spatial Resolution Satellite Image Simulation Based on 3D Data and Existing Images

Abstract

This study proposes an approach for simulating high spatial resolution satellite images acquired under arbitrary sun-sensor geometry using existing images and 3D (three-dimensional) data. First, satellite images, having significant differences in spectral regions compared with those in the simulated image were transformed to the same spectral regions as those in simulated image by using the UPDM (Universal Pattern Decomposition Method). Simultaneously, shadows cast by buildings or high features under the new sun position were modeled. Then, pixels that changed from shadow into non-shadow areas and vice versa were simulated on the basis of existing images. Finally, buildings that were viewed under the new sensor position were modeled on the basis of open library-based 3D reconstruction program. An experiment was conducted to simulate WV-3 (WorldView-3) images acquired under two different sun-sensor geometries based on a Pleiades 1A image, an additional WV-3 image, a Landsat image, and 3D building models. The results show that the shapes of the buildings were modeled effectively, although some problems were noted in the simulation of pixels changing from shadows cast by buildings into non-shadow. Additionally, the mean reflectance of the simulated image was quite similar to that of actual images in vegetation and water areas. However, significant gaps between the mean reflectance of simulated and actual images in soil and road areas were noted, which could be attributed to differences in the moisture content.

1. Introduction

For more effective use of satellite measurements, it is necessary to understand the type of performance that can be expected from the sensors used. Therefore, generation of diverse images simulated under the same performance conditions is needed (Blaschke, 2010). Simulation of MS (MultiSpectral) high spatial resolution images requires detailed simulation of both spatial and spectral complexity of actual scenes (Schott et al., 1999), particularly 3D reconstruction of buildings under oblique viewing angles. Schott et al. (1999) developed a model known as DIRSIG (Digital Imaging and Remote Sensing Image Generation) for simulating multi- and hyperspectral images using a ray-tracing technique. This method requires the use of significant amounts of ground truth data to define the spatial and spectral properties of observed elements. An additional method was developed by Doz et al. (2010) to simulate optical images in urban areas. Similar to other high spatial resolution image simulators, this model uses a set of airborne data acquired at different viewing angles, which is not always available.

In recent years, numerous MS and hyperspectral remote sensing satellites with diverse resolutions have been built, launched, and operated (Chen et al., 2013). In particular, a large number of satellite images can be freely obtained through open-access satellite data sources provided by the USGS (United States Geological Survey). Additionally, open data sources such as Google Street View (https://www.google.com/maps/streetview/) and Naver Street View (http://map.naver.com) provide a huge number of usable images of building sides that can be used in 3D model reconstruction of buildings. It is therefore necessary to conduct scientific research that takes advantage of datasets containing a substantial amount of image data.

This study develops a method for simulating high spatial resolution images acquired under arbitrary sun-sensor geometry based on existing images and 3D building models. Satellite images, which have significant differences in spectral regions compared with simulated images, are transformed to the same spectral regions as those of simulated images by using the UPDM. In addition, shadows cast by buildings or high features under the new sun position are modeled, and pixels changing from shadow into non-shadow areas and vice versa are simulated on the basis of existing images. Moreover, buildings viewed under the new sensor position are modeled on the basis of a 3D reconstruction program built from open libraries. An experiment is conducted to simulate a WV-3 image acquired under two different sun-sensor geometries based on a Pleiades 1A image, an additional WV-3 image, a Landsat image, and 3D building models.

 

2. Study Sites and Materials

Because this study focuses on high spatial resolution imagery, it is difficult to obtain DEMs (Digital Elevation Models) and building models covering a wide area; therefore, the study area was limited. In this case, the study area is Konkuk University campus located in Gwangjin-gu, Seoul, Republic of Korea (Fig. 1). Data used include a Pleiades 1A MS image, two Worldview-3 MS images, a Landsat 8 image downloaded from http://earthexplorer.usgs.gov/, building models with texture collected in September 2013, and a DEM established from LiDAR (Light Detection And Ranging) data acquired by the Optech ALTM 3070 system on September 1, 2013, at an altitude of 1300 m. Basic specifications of the images are shown in Table 1. The sub-Pleiades 1A image and two Sub-WV-3 images are shown in Fig. 2.

Fig. 1.Study area and WV-3 image in natural color

Table 1.Data specification in WV-3 simulation

Fig. 2.Pleiades 1A and WV-3 images, (a) Pleiades 1A acquired on October 16, 2013, (b) WV-3 image acquired on September 19, 2014, (c) WV-3 image acquired on August 7, 2015

The spectral patterns of vegetation and soil were obtained from the ASTER (Advanced Spaceborne Thermal Emission and Reflection Radiometer) spectral library (http://speclib.jpl.nasa.gov/), and the spectral pattern of sea water was obtained from the USGS spectral library (http://speclab.cr.usgs.gov/). Pleiades 1A, WV-3 images, and the Landsat image have different spatial resolutions and were acquired by different sensors on different dates. Therefore, they were geometrically registered to the UTM (Universal Transverse Mercator) zone 52N coordinate system and were resampled to 2 m. Furthermore, to prepare data for UPDM described below, atmospheric correction was performed using 6S (Second Simulation of the Satellite Signal in the Solar Spectrum) (Vermote et al., 1997) to obtain surface reflectance data. A sub-image covering Konkuk University–Seoul campus was then extracted.

Table 2.Sun-sensor geometry of simulation and satellite images used in the two experimental cases

 

3. Methodology

3.1 Universal pattern decomposition

UPDM is a sensor-independent method developed on the basis of the spectral decomposition method (Zhang et al., 2006), which is a type of spectral mixing analysis (Adams et al., 1995). This method decomposes each reflectance value of each pixel into a linear sum of three fixed standard spectral patterns: water, vegetation, and soil. If necessary for more detailed analysis, a supplementary standard spectral pattern can be used. The formulation is expressed as follows (Daigo et al., 2004):

where Ri : the reflectance of band i, Cw, Cv and Cs : the decomposition coefficients for water, vegetation, and soil, respectively, and Piw, Piv and Pis : the respective standard spectral patterns of water, vegetation, and soil, respectively. For simplicity, Eq. (1) can be expressed in a matrix form as follows:

where R = (R1,R2,…,Rn)T: the column vector of observations, n: the number of spectral bands, P = (Pw,Pv,Ps): the n × 3 matrix of which the row vector is the standard spectral pattern for band number n, C = (Cw,Cv,Cs)T: the column vector of coefficients, and r: the residual column vector for band i.

The pattern decomposition coefficients of Eq. (3) are evaluated using the least squares method:

In principle, in the sensor-independent UPDM, nearly identical values of Cw, Cv and Cs are obtained for data from the same object regardless of the sensor (Zhang et al., 2006). In practice, the precision of the coefficients improves as the number of bands increases (Daigo et al., 2004).

For each sensor band, the standard spectral patterns of each band Piw, Piv and Pis are calculated as follows (Liu et al., 2009):

where λs(i) and λe(i): the start and end wavelengths for band i, respectively, : the wavelength width of band i, Si: the SRF (Spectral Response Function) of a sensor of band i, and Pk(λ): the normalized standard pattern, which is fixed for use with all sensors and is defined as

where Rk(λ): the spectral reflectance patterns of standard objects, and ∫ dλ: the integration of the total wavelength range from 400 nm to 2500 nm.

According to Liu et al. (2009), Pk(λ) satisfies the following normalization Eq. (7):

The integral was calculated by interpolation between wavelengths using trapezoidal computation.

In recent studies, standard spectral patterns were either measured in the field with a radiospectrometer (Daigo et al., 2004; Zhang et al., 2006) or were collected from a historic dataset (Muramatsu et al., 2000). However, these types of data were not available for the present study.

3.2. Shadow-casting algorithm

In this study, shadows were modeled on the basis of the method reported by Richens (1997) and Ratti and Richens (2004). The approach simply computes shadow volumes as a DEM. First, three components of the vector pointing toward the sun are defined. Then, the components of an opposite vector are computed. Afterward, the DEM is translated by the x and y components; simultaneously, its height is reduced by subtracting the z component. The results derived from this step are considered as part of the shadow volume. This process is continued by multiples of this vector, the translated volume is compared with that previously calculated, and the maximum value is taken. The process can be stopped when all levels are zero or the translation has shifted the volume off the image. Subsequently, the original DEM is subtracted from the shadow volume. Finally, shadows are identified by pixel value of the shadow volume; pixels with negative or zero values are in light, and positive values are in shade.

3.3. Reconstruction of 3D building model and oblique image generation

One important phase in oblique image generation is to reconstruct and view 3D building models. For this purpose, a program was built on the various open libraries listed in Table 3.

Table 3.Open libraries used in 3D models reconstruction

The input data for this module includes the building models with texture information and a DTM (Digital Terrain Model) in 3DS format, which is one of the file formats used by the Autodesk 3DS Max 3D modeling, animation and rendering software. The Konkuk–Seoul Campus, houses 50 buildings; 50 3DS files consisting of the object information and texture image were used for buildings, and one 3DS file was used for terrain. The overall structure of the program is shown in Fig. 3.

Fig. 3.Structure of 3D building model reconstruction program

In this program, the Viewer module consists of two sub-modules for processing 2D texture image (2D Control) and the 3D model (3D Control). The 3D model viewer operates in multi-thread. The main thread and main shared data are shown in Tables 4 and 5.

Table 4.Main threads of 3DS processing

Table 5.Main shared data in 3DS processing

When shared data is being accessed by a thread, it is automatically locked, and the other thread cannot access this data. However, the access time is optimized to be as short as possible so that the accessing of all data can be taken place at almost the same time.

If the ground reflectance is Lambertian, the generation of simulated images under the off-nadir view angle requires simulation of pixels on the building side and changing them from shadow into non-shadow and vice versa. We labeled the simulated image as acquired under the sun-sensor geometry of Pleiades 1A image (Table 1) as image A and that acquired under the sun-sensor geometry of WV-3 image (Table 1) as image B. Shadows cast by buildings were generated by using the method described in section 3.2. For a given region, data acquired by different satellite sensors are comparable and correlated with each other after radiometric calibration, geometric rectification, and atmospheric correction (Zhu et al., 2010). Based on this statement, pixels turning from shadow into non-shadow and vice versa were simulated on the basis of the available high spatial resolution image (the other WV-3 image) or the Landsat 8 image acquired on September 19, 2013. If a pixel was in shadow in image A but not in shadow in image B, the pixel value in image B was replaced by the value in the remaining WV-3 image if possible, or the value in the Landsat 8 image. If a pixel was not in shadow in image A but in shadow in image B, the shadow is searched in pixels in a window of 31 × 31 pixels, or about 2 × 2 pixels of 30 × 30 m in Landsat 8. Then, linear regression was performed on these pixels to determine the relationship of reflectance in image A against that in the remaining WV-3 image or the Landsat 8 image resampled to 2 m. If the coefficient of determination is greater than 0.5, the linear regression equation will be used to compute the reflectance of the central pixel in image B. Otherwise, the central pixel reflectance will be replaced by the mean reflectance of the shadow pixels searched in the previous step. For buildings viewed from the off-nadir angle, pixel values will be replaced with building pixels generated by the 3DS model described above. An overall flowchart of oblique image generation is shown in Fig. 4.

Fig. 4.Oblique image generation

3.4 Overall experiment flow

In this experiment we simulated s WV-3 image acquired on October 16, 2013, under two different sun-sensor geometries: 1) a sensor-zenith angle of 41.6°, sensor-azimuth angle of 191.8°, sun-zenith angle of 39.40°, and sun-azimuth angle of 151°, which are the same values as those of the actual WV-3 image acquired on September 19, 2014; and 2) a sensor-zenith angle of 27°, sensor-azimuth angle of 90.1°, sun-zenith angle of 27.2°, and sun azimuth angle of 134.7°, which are the same as those in the actual WV-3 image acquired on August 07, 2015. These values are based on a Pleiades 1A image acquired on the same date but under different sun-sensor geometries (Table 1) and existing data including 3D building model and existing satellite images (Table 1). However, modeling the BRDF (Bidirectional Reflectance Distribution Function) requires many samples, which are not available in this research; therefore, the surface reflectance was assumed to be Lambertian. The overall flowchart used is illustrated in Fig. 5. First, several preprocessing steps were performed to obtain the Pleiades 1A-ground reflectance image. UPDM was then conducted on this ground reflectance image and on three typical spectral reflectance datasets derived from ASTER and USGS spectral library to transform the existing reflectance image (Pleiades 1A) into a WV-3 reflectance image. Because the building was simulated on the basis of the 3DS model, UPDM was not performed on the building areas. However, the shadows cast by buildings and pixels the changed from shadow into non-shadow were simulated. Subsequently, the WV-3 oblique image acquired at the new sun-sensor geometry was generated. Finally, the simulated image was compared with the actual image for quality assessment.

Fig. 5.Experiment flowchart of WV-3 image simulation

3.5 Quality assessment

To assess the quality of the simulated WV-3 image, it was compared with the WV-3 images acquired on September 19, 2014, and August 07, 2015, for cases #1 and #2, respectively. However, the actual images were acquired on different dates, on September 19, 2014, and August 7, 2015, compared with the simulation date of October 16, 2013. After more than one year, several landcover changes could have occurred in the study area. For instants, according to the Konkuk University administration several new buildings have been built, this can be observed in Fig. 2(c) (red rectangle indicated the construction area). Moreover, the images were acquired under different sun-sensor geometries resulting in geometric and radiometric distortions. Therefore, it is not reasonable to compare the images by pixels. For this reason, the simulated image was compared with the actual image by visual visual inspection and comparison of the mean spectral profiles by four landcover classes: water, soil, vegetation, and road.

 

4. Experimental Results

After performing UPDM on the Pleiades 1A reflectance image by MATLAB, the WV-3 reflectance image simulated as that acquired on October 16, 2013, with the same sensor and sun geometry as those of the Pleiades 1A image (Table 1) was obtained as shown in Fig. 6.

Fig. 6.Simulated WV-3 image acquired on October 16, 2013, under the same sun-sensor geometry as that used for the Pleiades 1A image

By using the 3DS building modeling program, buildings with textures were reconstructed, and a new viewing angle was configured. Shadows cast by buildings were generated as shown in Fig. 7(a). On the basis of this shadow image, we simulated pixels that changed from shadow to non-shadow and vice versa (section 3.3), as illustrated in Fig. 7(b). Subsequently, building models viewed with the sun-sensor geometry of WV-3 (Table 1) were generated. A subset of the final oblique image and the actual image, as indicated by red rectangle in Fig. 7(b), is shown in Figs. 7(c) and (d).

Fig. 7.Simulation of oblique WV-3 image with the aid of Landsat 8 only (a) shadow image, (b) simulated image under new sun position, (c) sub-oblique image, (d) sub-actual WV-3 image

As shown in Fig. 7(c) the area changing from shadow to non-shadow was not simulated effectively because the Landsat image has coarser spatial resolution. Therefore, the remaining WV-3 image was also used for this simulation. The images used for each simulation case are described in Table 3. The results of simulation with the aid of both the Landsat image and the remaining WV-3 image are shown in Figs. 8 and 9.

Fig. 8.Simulation of oblique WV-3 image under sun-sensor geometry #1 with aid of the second WV-3 image and the Landsat 8 image (a) shadow image, (b) simulated image under new sun position, (c) sub-oblique image, (d) sub-actual WV-3 image

Fig. 9.Simulation of oblique WV-3 image under sun-sensor geometry #2 with aid of the second WV-3 image and the Landsat 8 image (a) shadow image, (b) simulated image under new sun position, (c) sub-oblique image, (d) sub-actual WV-3 image

 

5. Quality Assessment and Discussion

A visual inspection of the simulated and the actual images in natural color composite are shown in Figs. 7–9(b)–(d). The bottom rows in the figures show sub-simulated and sub-actual WV-3 reflectance images; the sub area is indicated by a red rectangle in Figs. 7–9(b). The building shapes in the simulated image are very similar to those in the actual image. Some clear landcover changes can be seen in Fig. 9(c) and (d) compared with imaged acquired in 2013, for instants the edge of the playground has been paved. The length of shadow cast by building in the simulated image is very similar to that in the actual images. Nevertheless, some shadow pixels had different colors. Additionally, the color changes of the ground are quite similar in the natural color images. However, obvious differences can be observed in the areas that changed from shadow to not in shadow and vice versa (red circles in Fig. 7(c)) and in the surface textures of the buildings. According to the Konkuk University administration, the building roof was re-painted in 2014; therefore, the roof color of some buildings is different than that shown in Fig. 2(b) and (c).

To evaluate the stability of spectral features, a comparison of spectral profiles was performed. However, the actual images and simulated image have different acquisition dates and viewing geometries; therefore, the comparison was performed only on mean spectral profiles by main landcover classes including water, vegetation, soil, and road. The mean spectral profiles are shown in Figs. 10 and 11. As shown in the figures, the same pattern was revealed in both simulation cases. The mean water-spectral profile of the simulated image is very similar to that of the actual image. In the cases of vegetation and soil class, the mean spectral patterns are also similar; however, gaps are present between them. This result is attributed to the following explanation: According to the KMA (Korean Meteorological Administration), in the early morning of October 16, 2013, a light rain event of 3.5 mm occurred, whereas on September 19, 2014, August 7, 2015, and on some days before these dates, no rain occurred. Therefore, higher moisture content in soil and vegetation could have been present on October 16, 2013, compared with that on September 19, 2014 and August 7, 2015. This resulted in a decrease in reflectance throughout the VIS (VISible) to NIR (Near-InfRared) region (Hoffer, 1978). Similarly, in the case of vegetation, an increase in the moisture content led to a decrease in reflectance throughout the VIS to NIR region (Carter, 1991). However, the gap between the mean simulated reflectance profile and the mean actual reflectance profile in the vegetation class is very narrow. In the case of a road, its surface after a rain could be cleaner, which could have led to the slight increase in reflectance. In addition, BRDF was not considered. Moreover, the atmospheric correction result of multi-temporal high spatial resolution images from different satellite sensors was inconsistent owing to unstable radiometric calibration coefficients and different sun-sensor geometries (Lee and Lee, 2015).

Fig. 10.Mean spectral reflectance of simulated and actual images by landcover class and simulation under sun-sensor geometry #1 (September 14, 2014)

Fig. 11.Mean spectral reflectance of simulated and actual images by landcover class and simulation under sun-sensor geometry #2 (August 7, 2015)

The procedure for simulating high spatial resolution imaged is very complicated and is not a trivial task, particularly in urban areas with numerous high buildings. The replacement of a Landsat pixel value when a pixel turns from shadow to non-shadow in the simulated image is acceptable if the Landsat pixel is a pure (homogenous) pixel covered by only one land type. However, if it is a mixed (heterogeneous) pixel, this replacement could not simulate the actual landcover type of high-resolution pixel. This resulted in significant errors in the simulated image. However, if a high spatial resolution image was used in this replacement, the result was improved. It is expected that a substantial number of high spatial resolution images can be used to determine a relationship between the reflectance of a pixel that changed from shadow to non-shadow as well as that of a normal pixel and changes in sun-sensor geometry. This can significantly improve the quality of high spatial resolution image simulation. One important phase in high spatial image simulation is building simulation viewed from an off-nadir angle. Although the building shape can be effectively modeled, the texture information is not that simple. Because the texture image of a building is captured by a common camera or is extracted from some open data sources, this method is simple and provides high resolution texture images. However, the spectral information is poor.

On August 7, 2015, the level of suspended minerals or chlorophyll in water could be increased owing to environmental changes over a long period of time; as a result, the reflectance of water slightly increased in the VIS region (Jensen, 2007). However its mean reflectance was almost the same. Moreover, the mean reflectance of the road was significantly reduced in the NIR region, which could be attributed to changes in quality of the asphalt road over a one-year period.

 

6. Conclusion

The objective of this study was to develop a model for satellite high spatial resolution image simulation based on existing images and 3D data. The model was used to simulate a WV-3 image acquired under two different sun-sensor geometries based on a Pleiades 1A image, an additional WV-3 image, a Landsat image, and 3D building models. The results showed that the building shape was effectively modeled, although problems were noted in the simulation of pixels changing from shadows cast by buildings into non-shadow. Additionally, the mean reflectance of the simulated image was quite similar to that of actual images in vegetation and water areas. However, significant gaps were noted between the mean reflectance of simulated and actual image in soil and road areas, which may be attributed to differences in moisture content, BRDF, or landcover changes over long periods of time.

Given its flexibility, the proposed approach considering 3D geometry of the scene can be used to simulate images acquired under arbitrary sun-sensor geometry; therefore, it can be used in flight path simulation. Moreover, when additional super-high-resolution satellite sensors are launched in the future, satellite images with finer resolution will become easily accessible and will therefore advance the capabilities and quality of the proposed approach. However, developing an optimal method for simulation of pixels changing from shadow into non-shadow, or a simple method for 3D reconstruction, is likewise important for future studies.