End-to-End Image Simulator for Optical Imaging Systems: Equations and Simulation Examples

Abstract

The theoretical description of a simplified end-to-end software tool for simulation of data produced by optical instruments, starting from either synthetic or airborne hyperspectral data, is described and some simulation examples of hyperspectral and panchromatic images for existing and future design instruments are also reported. High spatial/spectral resolution images with low intrinsic noise and the sensor/mission specifications are used as inputs for the simulations. The examples reported in this paper show the capabilities of the tool for simulating target detection scenarios, data quality assessment with respect to classification performance and class discrimination, impact of optical design on image quality, and 3D modelling of optical performances. The simulator is conceived as a tool (during phase 0/A) for the specification and early development of new Earth observation optical instruments, whose compliance to user’s requirements is achieved through a process of cost/performance trade-off. The Selex Galileo simulator, as compared with other existing image simulators for phase C/D projects of space-borne instruments, implements all modules necessary for a complete panchromatic and hyper spectral image simulation, and it allows excellent flexibility and expandability for new integrated functions because of the adopted IDL-ENVI software environment.

1. Introduction

Hyper-spectral imaging has dramatically changed the rationale of remote sensing of the Earth relying on spectral diversity.

Since the pioneering Hyperion mission launched in 2001 [1], hyper spectral imaging airborne and satellite sensors have shown their utility by obtaining calibrated data for determining a wide variety of bio- and geophysical products from the collected imagery.

However, all sensors have their own set of performance characteristics, response functions, noise statistics, and so on, which determine and can challenge the validity of the generated data products. Through simulation of the sensor response, the utility of a new sensor design can be ascertained prior to construction, by running algorithms on simulated remote sensing data sets. In the case of existing well-characterised sensors the generation of simulated data assists in debugging sensor problems and provides a better understanding of a particular sensor’s performance in new operational environments.

In this paper, an end-to-end Selex Galileo (SG) simulation tool developed in the ENVI-IDL [2] environment for the generation of simulated data from airborne/space-borne optical and infrared instruments, starting from high resolution imagery is presented.

High resolution hyper-spectral data from airborne campaigns can be typically used as input for space-borne sensors simulations. As an alternative, the input images can be completely synthesized by modelling the geometrical and spectral characteristics of the observed targets. The simulator is based on six different modules describing the reflectance scenario, the atmospheric conditions, the instrument models and the atmospheric inversion model.

The core modules aim to simulate instrument performances (spectral, spatial, and radiometric) from a variety of sensor parameters including optics, detector, scanning, and electronics characteristics. The Atmospheric module is based on the standard Modtran [3] model, whereas the scenario simulation module aims at associating a spectral signature to each pixel of a synthetic thematic map, whenever a high resolution image taken by an airborne instrument is not available.

Compared to a detailed instrument simulator, typically developed for the realization and commissioning phases (B/C phases) of a spaceborne/airborne payload, the proposed simplified end-to-end simulator is conceived as a tool (phase 0/A) to enable the rapid dimensioning of a new optical instrument and to trace the link between user and instrument requirements. SG simulator (SG_SIM) pursues a similar philosophy as other approaches useful for 0/A phases (e.g., SENSOR, MODO, CAMEO, and PICASSO), and it includes all main functions (implemented in the IDL-ENVI SW environment) necessary for a complete hyper spectral image simulation, which are not often simultaneously present in the others.

For instance, in comparison to SENSOR [4] the control of spectral mixing and the generation of synthetic scenes are also considered, whereas in comparison to US simulators, for example CAMEO [5] and PICASSO [6–8]), the extension to the MWIR/LWIR spectral bands, and a 3D reflectance rendering are missing.

After a detailed theoretical description of SG_SIM model equations and its key concepts (Section 2), some simulation examples for satellite and airborne hyper spectral and panchromatic data study cases are reported (Section 3).

2. Simulator Equations Description

The flow diagram of the software tool is shown in Figure 1. The input data can be either airborne reflectance images at high spatial, spectral, and radiometric resolution or synthetic reflectance maps, coming from a thematic map and a reflectance data base, and specifications for the instrument to be simulated (e.g., spatial and spectral response, sampling, transfer function, noise model, viewing geometry, and quantisation).

Figure 1 

Airborne/space-borne optical sensor data simulator (flow diagram).

The simulation procedure consists of four different processing steps. First the at-sensor radiance images are obtained by using the Atmospheric Modtran code, then the signal is spatially, spectrally, and radiometric degraded by applying the specific instrument response models to generate the instrument simulated radiance image.

2.1. Atmospheric Simulation

The Atmospheric Module ingests as input a reflectance image taken at high spatial and spectral resolution which is then transformed into sensor radiance images by using the atmospheric radiances and transmittances generated by the Modtran code.

A preliminary simplified atmospheric model has been used. It considers Lambertian surface scattering, near-nadir observation, no adjacency effects, and a flat Earth. The input spectral radiance 𝐿(𝜆,ℎ) for an observation sensor at altitude ℎ is obtained on the basis of the following relationship, derived from the radiative transfer model depicted in Figure 2. The radiance is described from the following: 𝐿(𝜆,ℎ)=𝐼_TOA-SUN(𝜆)⋅cos𝜃_sun𝑑²⋅𝜏_↓(𝜆)⋅𝜌(𝜆)𝜋⋅𝜏_↑(𝜆,ℎ)+𝐿_↑ATM(𝜆,ℎ), (1) with(i)𝐼_TOA-sun(𝜆) = Top of atmosphere sun irradiance (W/m²/μm);(ii)𝜌(𝜆) = Earth surface reflectance;(iii)𝜃_sun = Sun observation angle (function of latitude, longitude, day of the year, and time);(iv)𝑑² = Earth-Sun distance normalised to mean (depending from day of the year);(v)𝜏_↓(𝜆)= Total downwards atmosphere transmission;(vi)𝜏_↑(𝜆,ℎ) = Total upwards atmosphere transmission from ground to the observation altitude ℎ;(vii)𝐿_{↑𝐴𝑇𝑀}(𝜆,ℎ) = Scattered atmosphere radiance (W/m²/sr/μm) from ground to the observation altitude ℎ;(viii)𝐿(𝜆,ℎ) = Total Atmosphere radiance (W/m²/sr/μm) which represents the input to the instrument at altitude ℎ.

Figure 2 

General atmospheric model scheme: the simplified atmospheric simulator takes into account in first approximation the 1, 2, 4, and 6 paths. The adjacent effects described from 3 and 5 paths will be updated in the next version of simulator.

The downwards/upwards atmospheric transmittances 𝜏_↓(𝜆), 𝜏_↑(𝜆,ℎ) and the atmospheric radiance 𝐿_↑ATM(𝜆,ℎ) depend on the concentration of all atmospheric gases and the aerosols distribution. The simulator allows the control of the major variable atmospheric gases (i.e., the columnar water vapour and CO₂ contents), the aerosols visibility at a certain observation altitude ℎ, and the aerosols profile. These parameters can be controlled by means of Modtran code inputs, while the other parameters are considered constant. A dedicated graphical interface is used to create the Modtran input charts.

Generally the surface reflectance’s images 𝜌(𝜆), used as input to the simulator, come from a data base of experimental airborne or ground truth data acquired with other spectrometers, and they are affected by the spectral response of those instruments used for the database acquisition.

The radiances 𝐿(𝜆,ℎ) are generated from the Modtran code at the maximum spectral resolution (1 cm⁻¹) and are convolved with the spectral response (SR) of the instrument used to generate the data base. This spectral response (SR) is approximated with a Gaussian function with the centre wavelength 𝜆_𝑐 and the Full Width at Half Maximum (FWHM) equal to Δ𝜆/4, where the integral is performed in a Δ𝜆 spectral range centred in 𝜆_𝑐, and can be written as in the following: 𝐿_1𝑖(𝜆,ℎ)=^𝜆𝑐_𝜆^+Δ𝜆/2_{𝑐−Δ𝜆/2}𝐿_𝑖(𝜆,ℎ)⋅SR_𝑖𝜆_𝑐(𝑖)−𝜆⋅𝑑𝜆=Δ𝜆_{𝜆𝑐𝜆+Δ𝜆/2𝑐−Δ𝜆/2}𝐼_TOA-SUN(𝜆)⋅cos𝜃_sun𝑑²⋅𝜏_↓(𝜆)⋅𝜌_𝑖𝜆_𝑐(𝑖)𝜋⋅𝜏_↑(𝜆,ℎ)+𝐿_↑ATM(⋅𝜆,ℎ)SR_𝑖𝜆_𝑐(𝑖)−𝜆⋅𝑑𝜆≅Δ𝜆_{𝜆𝑐𝜆+Δ𝜆/2𝑐−Δ𝜆/2}𝐼_TOA-SUN(𝜆)⋅cos𝜃_sun𝑑²⋅𝜏_↓(𝜆)⋅𝜌_𝑖𝜆_𝑐,Δ𝜆𝜋⋅𝜏_↑(𝜆,ℎ)+𝐿_↑ATM(⋅𝜆,ℎ)SR_𝑖𝜆_𝑐(𝑖)−𝜆⋅𝑑𝜆=Δ𝜆_{𝜆𝑐𝜆+Δ𝜆/2𝑐−Δ𝜆/2}𝐿_𝑖(𝜆,ℎ)⋅SR_𝑖𝜆_𝑐(𝑖)−𝜆⋅𝑑𝜆⋅Δ𝜆 (2) With(i)𝐿_𝑖(𝜆,ℎ)=[𝐼_TOA-SUN(𝜆)⋅cos(𝜃_sun)/𝑑²⋅𝜏_↓(𝜆)⋅𝜌_𝑖(𝜆_𝑐(𝑖))/𝜋⋅𝜏_↑(𝜆,ℎ)+𝐿_↑ATM(𝜆,ℎ)] the output spectral radiance obtained from Modtran by using the real surface reflectivity 𝜌_𝑖(𝜆) for each spectral pixel 𝑖;(ii)SR_𝑖[𝜆,𝜆_𝑐(𝑖)] the normalised spectral response of instrument used to generate the data base for the 𝑖th spectral channel (with 𝜆_𝑐(𝑖) the central wavelength) as a function of wavelength 𝜆. Each SR_𝑖[𝜆,𝜆_𝑐(𝑖)] has been simulated with a Gaussian response centred at 𝜆_𝑐 and FWHM equal to Δ𝜆/4.(iii)𝐿_𝑖(𝜆,ℎ)=[𝐼_TOA-SUN(𝜆)⋅cos(𝜃_sun)/𝑑²⋅𝜏_↓(𝜆)⋅𝜌(𝜆_𝑐(𝑖),Δ𝜆)/𝜋⋅𝜏_↑(𝜆,ℎ)+𝐿_↑ATM(𝜆,ℎ)] the output spectral radiance from Modtran by using the data base reflectivity value.(iv)𝜌_𝑖(𝜆_𝑐(𝑖),Δ𝜆)=(∫^{𝜆𝑐+Δ𝜆/2}_{𝜆𝑐−Δ𝜆/2}𝜌_𝑖(𝜆_𝑐(𝑖))⋅SR_𝑖[𝜆_𝑐(𝑖)−𝜆]⋅𝑑𝜆)/Δ𝜆 the weighting mean of the Earth surface reflectivity within the SR(𝜆) spectral response, which represents the data-base reflectivity value.(v)𝐿_1𝑖(𝜆,ℎ)(W/m²/sr/μm) the mean spectral radiance within the SR(𝜆) spectral response.

The high resolution at sensor radiance 𝐿₁(𝜆,ℎ) is simulated with Modtran code for different values of the surface reflectivity and a 3D Look-Up-Table (reflectivity, radiance, and wavelength) is generated. Finally, for each wavelength, the simulation module determines the best linear fit between radiance and surface reflectivity, which is applied to all input reflectivity image pixels 𝜌_𝑖(Δ𝜆) to generate the at sensor radiances 𝐿_1𝑖(𝜆,ℎ).

2.2. Spectral Degradation

The second processing block applies a spectral degradation where the at-sensor radiance image is further spectrally degraded to the spectral channels and response of the airborne/satellite instrument to be simulated by means of a spectral interpolation and a convolution with the Instrument Chromatic Response (ICR). The ICR_𝑖[𝜆,𝜆_𝑐(𝑖)] represents the normalized to maximum instrument response for the 𝑖th spectral channel (with 𝜆_𝑐(𝑖) defined as the central wavelength) to a spatially uniform monochromatic source as a function of wavelength 𝜆. The at-sensor radiance 𝐿₂(W/m²/sr) is obtained from the following: 𝐿₂𝜆_𝑐(𝑖),ℎ=𝐿₁(𝜆,ℎ)⋅ICR_𝑖𝜆_𝑐(𝑖)−𝜆⋅𝑑𝜆. (3)

2.3. Spatial Degradation

The spatial degradation module ingests the at-sensor radiance image and degrades it to the required spatial sampling. This process is applied by means of a convolution between the input image and the Instrument Spatial Response (ISR) of the optical sensor to be simulated followed by a resampling process (decimation) (Figure 3). The ISR is defined as the response of the overall instrument in a given spatial pixel to a spectrally uniform point source as a function of its position in space. The spatially degraded radiance image (𝐿₃) is described by the following: 𝐿₃𝜆_𝑐(𝑖),ℎ,𝑥,𝑦=_image𝐿₂𝜆_𝑐(𝑖),ℎ,𝑥^′,𝑦^′⋅ISR𝑥−𝑥^′,𝑦−𝑦^′⋅𝑑𝑥^′⋅𝑑𝑦^′. (4)

Figure 3 

Spatial degradation module of the instrument simulator: convolution between the high spatial resolution image and the instrument spatial response (ISR) of the optical sensor that is to be simulated, followed by a resampling process (decimation).

The ISR is calculated as the Inverse Fourier Transform of the Modulation Transfer Function (MTF), which assumes the overall system is a linear shift invariant system. Then a “cascade model” for the system MTF is applied on the hypothesis of independent subsystems.

The hypothesis of independent subsystems is exact for many instruments, while the use of MTF, without taking into account the phase effects is valid only as a first approximation in incoherent imaging systems using well-corrected optics [10].

The “cascade model” (Figure 4) takes into account the hypothesis of separability of spatial frequency variables. Due to the properties of the Fourier Transform, the separability in the frequency domain corresponds to separability in the space domain. The along-track and across-track MTFs are calculated starting from a theoretical formulation and the Inverse Fourier Transform is calculated and normalized to a unit integral for both. In this way two unidimensional digital filters have been obtained and convolved with the high resolution image by means of the following: 𝐿₃𝜆_𝑐(𝑖)=𝐿,ℎ,𝑥,𝑦₂𝜆_𝑐(𝑖),ℎ,𝑥^′,𝑦^′⋅ISR_𝑥𝑥−𝑥^′⋅𝑑𝑥^′⋅ISR_𝑦𝑦−𝑦^′⋅𝑑𝑦^′, (5) where:(i)ISR_𝑥(𝑥)=IFFT(MTF_𝑥) the Instrument Spatial Response along 𝑥 (e.g., along-track),(ii)ISR_𝑦(𝑦)=IFFT(MTF_𝑦) the Instrument Spatial Response along 𝑦 (e.g., across-track).The along- and across-track MTF_𝑥,𝑦 are calculated taking into account the image degradation contributions reported in Table 1.

Table 1 

MTF subsystems contributions for a push broom system.

Figure 4 

Flux diagram to calculate the instrument spatial response of the cascade model.

The image quality can be affected from many factors such as the size of detector (spatial aperture), the detector degradations (e.g., pixel cross talk or charge transfer & reading smearing in CCD), the integration time during image motion (temporal aperture) caused from satellite motion or the scanning system (resp. for a push broom or a whiskbroom system), the electronic filtering, the focal plane jitter (instrument micro-vibrations), the optics diffraction and aberrations [11].

These components can influence both the across-track and/or the along-track MTF depending on the direction of scanning and the disposition of detector. Some of these components are described in annex.

Examples of simulated MTF and SRF functions for airborne and spaceborne instruments that were generated with the simulator are reported in Section 3.1.

2.4. Radiometric Degradation

The fourth processing module accounts for radiometric degradation. A random noise term is added to the images to simulate the 𝑁𝑒Δ𝐿 (Noise Equivalent Difference Radiance in W/m²/sr) of the optical instrument. The radiance 𝐿₃(𝜆_𝑐(𝑖),𝑥,𝑦) of each pixel (𝑥,𝑦) and of the 𝑖th spectral band (with central wavelength 𝜆_𝑐(𝑖)) is substituted with a random value taken from a Gaussian distribution, where 𝐿₃(𝜆_𝑐(𝑖),𝑥,𝑦) represents the mean radiance value (W/m²/sr) and 𝑁𝑒Δ𝐿 (W/m²/sr) the noise equivalent radiance, which is the standard deviation of the instrument temporal noise. The relationship between 𝑁𝑒Δ𝐿 and the pixel radiance value 𝐿₃ is described in the following:

𝑁𝑒Δ𝐿𝜆_𝑐(𝑖)=,𝑥,𝑦𝐴𝜆_𝑐(𝑖)+𝐵𝜆_𝑐(𝑖)⋅𝐿₃𝜆_𝑐(𝑖),𝑥,𝑦, (6) where 𝐴(𝜆_𝑐(𝑖)) is the noise variance of the detector (dark current, read-out and Johnson noises) plus FEE/ADC (Front End Electronics/Analog to Digital Converter) for the 𝑖th spectral band and the product [𝐵(𝜆_𝑐(𝑖))⋅𝐿₃(𝜆_𝑐(𝑖),𝑥,𝑦)] is the photon noise variance, which is proportional to the input signal [𝐿₃(𝜆_𝑐(𝑖),𝑥,𝑦)]. A Gaussian distribution function for the noise is a good approximation also for the photon noise, because the Poisson distribution approximates a Gaussian function for a high number of generated photocarries.

𝑁𝑒Δ𝐿 is the minimum variation of the input radiance which can be measured and represents the radiometric resolution of the instrument. Another representation of the sensor noise can be derived from the signal-to-noise ratio (SNR) for each pixel and for each wavelength. The SNR can be obtained from the following: SNR_𝑖(𝜆)=𝐿₃𝜆_𝑐(𝑖),𝑥,𝑦𝑁𝑒Δ𝐿𝜆_𝑐(𝑖)=𝐿,𝑥,𝑦₃𝜆_𝑐(𝑖),𝑥,𝑦𝐴𝜆_𝑐(𝑖)+𝐵𝜆_𝑐(𝑖)⋅𝐿₃𝜆_𝑐(𝑖).,𝑥,𝑦 (7) The photon noise formulation reported in (6) ([𝐵(𝜆_𝑐(𝑖))⋅𝐿₃(𝜆_𝑐(𝑖),𝑥,𝑦)] is based on the relationship between the number of acquired electrons 𝑁𝑒 and the integrated input radiance 𝐿_IN (W/m²/sr) for each spectral channel as described in the following: 𝑁_𝑒(𝑒⁻)=𝐿_INWm²sr⋅𝐴_pupilm²⋅Ω(sr)⋅𝜏⋅𝜆Jℎ⋅𝑐⁻¹⋅𝑡_int(sec)⋅𝜂e⁻/phot=𝐾⋅𝐿_INW/m²sr, (8) with(i)𝜏 = Total mean in band (with ICR) instrument transmittance;(ii)𝐴_pupil = Input pupil area (m²);(iii)Ω= Scene pixel IFOV (sr);(iv)𝑡_int = Integration time (sec);(v)𝜂= Mean in band (with ICR) detector quantum efficiency (electrons/photon);(vi)ℎ𝑐/𝜆=(1,98E−19 Joule*μm)/𝜆 = Energy of a photon of wavelength 𝜆 (μm);(vii)𝐿_IN (W/m²/sr) = Spectrally integrated mean radiance in the ICR;(viii)𝐾(e⁻/(W/m²sr))=𝐴_pupil(m²)⋅Ω(sr)⋅𝜏⋅(𝜆/(ℎ⋅𝑐))⋅𝑡_int(sec)⋅𝜂(e⁻/phot)) the coefficient of proportionality between the number of acquired electrons 𝑁_𝑒 and the input radiance 𝐿_IN.

The photon noise equivalent difference radiance 𝑁𝑒Δ𝐿 is related to the photon noise equivalent difference electron 𝑁𝑒Δ𝑁_𝑒, which can be obtained from the standard deviation of the Poisson noise distribution. This standard deviation is equal to the square root of the number of electrons itself and is described from the following: 𝑁𝑒Δ𝑁_𝑒=𝑁_𝑒,𝑁𝑒Δ𝐿=𝑁_𝑒𝐾=𝐾⋅𝐿_IN𝐾=𝐿_IN𝐾≡𝐵⋅𝐿_IN,𝐵W/m²sre⁻=1𝐾=1𝐴_pupil⋅Ω⋅𝜏⋅𝜆/(ℎ⋅𝑐)⋅𝑡_int⋅𝜂. (9)

The 𝐴 and 𝐵 coefficients, which depend from the selected spectral channel and are fed as input to the simulator, can be derived from the radiometric model of the simulated optical sensor or they can be evaluated on the basis of acquired images of homogeneous targets acquired by the sensor [13–16].

Two additional procedures have been implemented to permit the analysis of the simulated images (Sections 2.5 and 2.6).

2.5. Atmospheric Correction

The first permits the retrieval of surface reflectance from airborne and spaceborne sensor radiances. Two standard methods can be used: one is based on Modtran code, by inverting (1), to obtain the surface reflectance from the instrument radiance and the second based on the standard ENVI-FLAASH [17] software, which allows aerosols to be estimated by means of the dark pixel method (water bodies, shadowed area, and dense vegetation) and the water vapour map to be estimated by means of the 820, 940, and 1135 nm absorption bands ratio method [18].

2.6. Synthetic Image Generation Module

The second procedure permits to quantitatively evaluate the impact of instrumental parameters on simulated image quality when a low noise airborne input image is not available.

In particular it allows the creation of black and white bar test images with different modulations (square or sinusoidal), periods and shading, to the scope to evaluate the impact on the image quality of instrument parameters such as MTF and noise as a function of the spatial sampling interval and the target reflectivity, and to analyse the minimum detectable albedo contrast as a function of spatial frequency and illumination conditions (Figure 5).

Figure 5 

Input mask to generate the synthetic black and white images.

It is also possible to generate synthetic hyper-spectral surface reflectance images at the desired spatial and spectral resolution by using as input a thematic map of the zone under investigation (derived synthetically or from a classification) and a spectral library of the surface materials of interest. A statistical mixing of spectral signatures for each zone with the Dirichlet method permits to control the percentage of the statistical variability [19].

The following further statistical variability, devoted to a better representation of a real scenario, can be introduced [20]: (i)a uniform or Gaussian variability for each spectral signature due to a possible spatial variation of the substance composition, such as contaminants, oxidation, and ageing, and so forth,(ii)a beta function distributed statistical variation of illumination, which takes into account possible image errors due uncompensated observation and surface slope angles,(iii)a Gaussian variability due to scenario noise, coming from uncompensated atmospheric and environment effects or uncompensated errors of sensors used to obtain the spectral library data.

Then the surface reflectance 𝑅, represented from a column vector for each wavelength, is obtained by means of the following matrix mixing relationship: 𝑅=P⋅𝛾⋅Ψ⋅A⋅𝑡+M⋅(1−𝑡)+N, (10) with(i)P a matrix representing the end-members (pure elements) reflectance for each wavelength;(ii)A a column vector representing the statistical variability of abundances, according to a Dirichlet distribution;(iii)M a column vector representing the mean value of abundances for each wavelength;(iv)𝑡 a scalar parameter which describes the degree of statistical mixing (𝑡=0 no mixing, 𝑡=1 all random mixing);(v)𝛾 a statistical parameter obtained from a beta distribution, describing the illumination variation for each pixel;(vi)Ψ a diagonal matrix derived from a uniform or a Gaussian density, representing the spatial variation of end members;(vii)N a column vector representing a Gaussian scenario noise, that is uncompensated atmospheric retrieval and/or sensors errors used to obtain the library data set.

3. Simulations Examples

Several simulation tests were performed to assess the potential of the tool for the instrument image quality and applications evaluation in the framework of the study and the testing phases of the Selex Galileo SIMGA airborne hyper spectral camera and the HypSeo (ASI-PRISMA precursor [21]) spaceborne hyper spectral and panchromatic cameras phase A study.

Such activities allowed also the validation of the simulator by means of real SIMGA data acquired on clay soil targets during an airborne campaign of 2009 in Mugello (Tuscany, I) test site, where ground truth data were collected simultaneously at the same time of overflights [16].

Some examples of simulation of a 3D map representation of the ISR function have been produced for the purpose of evaluating the instrument image quality, which is generally given from the FWHM of the instrument spatial response or by the ratio between the integral of the spatial response within a delimitated spatial domain (e.g., 1 Spatial Sampling Interval) and the integral in all spatial domain, which is generally called integrated energy (in percentage unit).

A more detailed analysis based on another image quality parameter (SNR*MTF) has been done to trade-off the image quality of a panchromatic camera as a function of some instrument parameters (e.g., pupil diameter and spatial sampling) for different atmospheric conditions (summer/winter and rural/urban aerosol) aiming to a better definition of the instrument requirements.

Finally VIS and SWIR radiance and reflectance simulated images have been generated for some specific targets related to civil (land use) and dual use applications for terrestrial and marine environments to the scope of understanding the instrument capabilities for targets’ discrimination. These two dual use applications has been simulated during the testing phase of the airborne SG instrument (SIMGA) by means of targets of small green panels over vegetation cover and small grey panels under water, and then verified by means of an airborne campaign on a controlled area.

3.1. SRF 3D Maps for Integrated Energy Calculations

The simulator permits a 3D representation of the SRF map by using as input a delta function. As an example this representation has been done to evaluate the spatial resolution (defined in terms of percentage of integrated energy of SRF within a certain space domain) of the airborne SG SIMGA hyperspectral camera by taking into account both the laboratory measurements and the smearing effect introduced from the detector integration which occurs during platform motion. The along and the across-track MTF and SRF contributions are displayed in Figures 6(a) and 6(b), respectively, for the VIS and the SWIR channels. The instrument parameters used in the simulations are reported in Tables 2(a) and 2(b). From the tables it appears that the ratio between the FWHM of the SRF and the Spatial Sampling Distance (SSD) is much lower for SWIR channels (0.87 along scan*1.05 across scan) with respect to the VIS ones (2.70 along scan*1.49 across scan), showing that the SRF of VIS channels has a larger width (in ±2/3 pixels) with respect to that of the SWIR ones (±1 pixel) (see also Figures 6(a) and 6(b)).

Table 2 

(a)

(b)

Figure 6 

(a) Maps of Spatial Response Function (SRF) of SIMGA VIS channels ((a), (b)), plots of SRFs across and along-track (along and across-slit) for the central FOV pixel as a function of on-ground displacement (c). The simulation is obtained by using the static laboratory measurements after having introduced the smearing effect due to the along-track motion (Tables 2(a) and 2(b)). MTF motion and along/across-track components as a function of spatial frequency at detector level (d). (b) Maps of Spatial Response Function (SRF) of SIMGA SWIR channels (a), (b), plots of SRFs across and along-track (along and across-slit) for the central FOV pixel as a function of on-ground displacement (c). The simulation is obtained by using the static laboratory measurements after having introduced the smearing effect due to the along-track motion (Tables 2(a) and 2(b)). MTF motion and along/across-track components as a function of spatial frequency at detector level (d).

The integrated energy calculation performed within an area of 1 SSD*1 SSD of the VIS and SWIR 3D maps also confirms that the energy content within a pixel is much lower for VIS respect to SWIR channels (the same happens within the same ground size of 1.333 m*1.333 m):(i)Integrated Energy in 0.706 m*0.706 m for VIS (1 SSD*1 SSD) = 19%,(ii)Integrated Energy in 1.333 m*1.333 m for VIS (1.9 SSD*1.9 SSD) = 51%,(iii)Integrated Energy in 1.333 m*1.333 m for SWIR (1 SSD*1 SSD) = 61%.

In conclusion the spatial resolution of VIS channels is coarser with respect to that of SWIR channels, also if the spatial sampling is better (0.706 m respect to 1.333 m).

A further exercise was done to simulate the 3D Hypseo SRF [16] by using the instrument parameters reported in Table 3 and the Hypseo MTF model [21]. The FWHM of the spatial response is 24.4 m*20.6 m (along-scan*across-scan) while the integrated energy in 1 SSD*1 SSD (20 m*20 m) is 53%, which is a value substantially equal to that estimated for the airborne SIMGA instrument at 1.33 m*1.33 m of pixel size.

Table 3 

Characteristics of HYPSEO Hyper-spectral and PAN cameras at 620 Km satellite altitude.

3.2. Satellite Panchromatic Image Quality Requirements

The simulator permits to study the impact of system design parameters on the instrument image quality. To this scope a parametric analysis of the performance of the HypSEO-PAN Camera as a function of the pupil diameter dimension for different spatial sampling, atmospheric and illumination conditions was performed on the basis of simulated test images and instrument parameters (Table 4), to trade-off the instrument sizing with the image quality.

Table 4 

Illumination/atmospheric parameters used in the Hypseo PAN camera simulations.

We have adopted as a first approximation of image quality criterion the Minimum Resolvable Contrast (MRC) at a certain spatial frequency 𝑓, which is equal to the inverse of the product [MTF(𝑓)⋅SNR], where SNR is calculated for uniform scenes (spatial frequency 𝑓=0) [11]. As a rule of thumb, the adopted value of MRC = 10% gives the following threshold relationship for target identification with spatial frequency 𝑓: MTF𝑓⋅SNR>10. (11)

Two different kind of input images have been used for the simulations: a surface reflectance image, based on the IKONOS panchromatic camera at ~1 m spatial sampling and a bar synthetic image at a spatial sampling of 0.5 m.

In Figure 7 a comparison between simulations of the Hypseo panchromatic image, obtained from the IKONOS image, at different spatial sampling intervals and different Hypseo pupil diameters is shown for the low radiance case (case “B” in Table 4). The product MTF(𝑓)⋅SNR has been calculated and the results are displayed in Table 5. The simulation with a high pupil diameter of 300 mm (case (b) in Table 5) is better from image quality point of view (lower GSD, high SNR*MTF, and optimum targets discrimination in Figure 7(b)) but it has a large impact on instrument sizing. All other images are strongly affected from diffraction, due to the pupil size of 150 mm, but the case with 5 m of GSD (case (d) in Table 5) seems better for SNR*MTF parameter and targets discrimination (in case (a) and (c) of Figure 7 the instrument noise overlays all other possible image features).

Table 5 

Image quality parameters related to the simulations of Figure 7 with low illumination radiance (Case B in Table 4). The diffraction limit represents the Airy radius (𝜆=1 𝜇m, 𝐻=620 Km) of Rayleigh criterion. Case (b) is better respect to case (d) because SNR*MTF is related to a signal with lower Nyquist period.

(a)

(b)

(c)

(d)

Figure 7 

Simulated HypSEO Pan images (from IKONOS) for low TOA radiance (18 W*m−2*sr−1*μm−1 at albedo = 0.27, case “B” in Table 4), different spatial sampling (GSD) and pupil diameter (𝐷): (a) GSD=2.5 m, 𝐷=150 mm (SNR=11, MTF(𝑓Nyquist) = 0.09); (b) GSD=2.5 m, 𝐷=300 mm (SNR=32, MTF(𝑓Nyquist) = 0.18); (c) GSD=3.5 m, 𝐷=150 mm (SNR=27, MTF(𝑓Nyquist) = 0.17); (d) GSD = 5 m, 𝐷=150 mm (SNR=50, MTF(𝑓Nyquist) = 0.22).

Another simulation with synthetic bars has been done to verify the previous results, by changing the sampling and the illumination conditions avoiding any possible effect coming from the degraded characteristics of the IKONOS image quality. In Figure 8 an HypSEO PAN simulation, from a synthetic bar image, at different spatial resolution (GSD) and pupil diameter (𝐷) is shown for a high (case “A”) and a low (case “B”) radiance case, with parameters represented in Table 4. The input synthetic image in the horizontal direction is composed by 5 sequences of grey-black bars, each consisting of 10 cycles at fixed period (5, 7, 10, 15, and 20 m). In the vertical direction the albedo range of grey bars is between 10% and 20%, whereas the albedo of black bars is constant (10%). The up and down arrows in Figure 8 indicate the periods for which the criterion MTF(𝑓)⋅SNR>10 is satisfied or not.

Figure 8 

HypSEO PAN simulation (from synthetic bar image) at different spatial resolution (GSD), pupil diameter (𝐷), and illumination conditions (high radiance case “A” and low radiance case “B” in Table 4); the contrast along the vertical axis changes between 10% and 20%; the up and down arrows indicate the periods for which the criterion (SNR*MTF>10) is satisfied or not.

The radiance (𝐿), the Signal to Noise Ratio (SNR), and the Modulation Transfer Function (MTF) values corresponding to the extreme simulated albedo values are represented in Table 6 while the product MTF(𝑓)⋅SNR has been reported in Table 7.

Table 6 

Radiance, SNR, and across-track MTF values corresponding to the extreme simulated albedo values, calculated for different ground spatial sampling (GSD), pupil diameter (𝐷) and the two simulation conditions reported in Table 4.

Table 7 

Image quality parameter (SNR*MTF) calculated at low and high illumination radiance (Table 4), 20% of on ground albedo and different spatial periods of Figure 8. The diffraction limit represents the Airy radius of Rayleigh criterion. Values with SNR*MTF > 10 are bold.

The results confirm that image quality as defined from this kind of metric is improved by increasing pupil diameter (from 150 to 300 mm) at equal spatial sampling (GSD = 2.5 m), because of an increased SNR and a reduced diffraction effect on MTF. The image quality is also improved as spatial sampling decreases from 2.5 to 3.5 m (pupil diameter = 150 mm) because of an increased SNR.

The HypSEO PAN nominal case Ground Spatial Sampling (GSD) of 5 m and pupil diameter (𝐷) of 150 mm seems a good compromise in terms of image quality with respect to the others, because the simulation results are not so different with respect to the case with GSD=3.5 m (𝐷=150 mm) (SNR*MTF is higher for low radiance case) and from an instrument design point of view it appears more feasible with respect to the best case with GSD=2.5 m (𝐷=300 mm).

For PAN nominal case the above MTF(𝑓)⋅SNR criterion is satisfied only for periods larger than the Nyquist period of 10 m at high radiances, but some oscillations affected by aliasing can also be observed at low periods as 7 m.

Finally an example of simulated Hypseo PAN image (GSD=5 m, 𝐷=150 mm) obtained from airborne high resolution MIVIS data in a forest environment has been performed (Figure 9(c)) to the scope of testing image fusion methods based on the sharpening of hyperspectral image by means of panchromatic observations [22].

(a)

(b)

(c)

Figure 9 

Airborne and simulated space-borne image on S. Rossore Park (I) test site: (a) MIVIS airborne reflectance (2.5 m), used as input; (b) simulated hyper-spectral HYPSEO radiance (20 m); (c) simulated PAN HYPSEO radiance (5 m) image.

3.3. Satellite Hyperspectral Land Use Classification

Another important use of the simulator has regarded the demonstration of potential applications of the HYPSEO SG spaceborne hyperspectral camera.

A simulation of the HYPSEO SG space-borne hyperspectral camera was performed by using as input the airborne MIVIS reflectance images acquired on a Tuscany (I) test site (S. Rossore Park and Arno River mouth) at 2.5 m spatial resolution [23]. The instrumental parameters are reported in Table 3 [24].

A MIVIS reflectance image is transformed into the satellite HYPSEO radiance (Figure 9) (𝐻=620 Km) by using the atmospheric model parameters of Table 4. Then the HYPSEO radiance image is obtained by means of a spectral resampling of the MIVIS image to the 210 spectral bands of Hypseo, with a Gaussian Instrument Chromatic Response (with FWHM=10 nm) and a spatial resampling to the HYPSEO spatial sampling interval of 20 m by using the simulated spatial response, and adding the noise by means of parameters coming from the HYPSEO radiometric model. Moreover an HYPSEO reflectance image has been obtained after removal of atmospheric effects introduced by MODTRAN code (Table 4).

A land use classification map based on the Spectral Angle Mapper (SAM) algorithm [25] from the HYPSEO simulated reflectance image is shown in Figure 10. The confusion matrix shows a good correlation between classified and ground truth data, compared with multispectral sensors [23], confirming the instrument capabilities for this kind of application.

Figure 10 

Example of forest SAM classification map of S. Rossore test site (Italy) based on simulated Hypseo hyper-spectral image (20 m). Confusion matrix shows a good correlation between classified and ground truth data.

3.4. Target “Camouflage” in Rural Background

The dual use capability for targets discrimination with camouflage panels embedded in vegetation has been evaluated during the testing phase of the SIMGA airborne hyperspectral instrument. To this scope some simulations was done during the SIMGA project phase. The simulated instrument SIMGA reflectance images been obtained by using the MODTRAN code in a standard atmospheric condition (Table 8), the measured instrument spatial response (Figures 6(a) and 6(b) and Tables 2(a) and 2(b)) and the instrument noise [16]. In Figure 11 a SIMGA reflectance image of simulated green panels over vegetation after FLAASH inversion algorithm is shown. The result of simulation showed that green panels were clearly distinguished respect to vegetation, because of their higher reflectance in the SWIR bands (1.2 and 1.6 micron), so validating the utility of the hyperspectral sensor for this kind of application. Moreover a validation of the simulation was obtained during an airborne campaign performed in S. Rossore park (Tuscany, I), where different green panels were placed over green grass. In Figure 12 the green panels are clearly distinguished in the SWIR bands while the contrast in VIS bands is negligible.

Table 8 

Main parameters for Modtran simulations.

(a)

(b)

Figure 11 

Simulated false color SIMGA radiance image (SWIR bands) of green panels over vegetation (a). The scenario has a dimension of 740*740 m2 and it is composed of 4 classes of background scenario: class 1 = green grass, class 2 = bare soil, class 3 = linear mixing of green grass (1/3), dry grass (1/3), and bare soil (1/3); class 4 = pinewood [linear mixing of green (1/2) and dry (1/2) needles]. 18 targets are used grouping in three different sizes (T_1 = 1 × 1 m2, T_2 = 2 × 2 m2, and T_3 = 4 × 4 m2) and three different types (light green paint, dark green paint, and mimetic paint). Mimetic paint is composed of 1/3 light green paint, 1/3 dark green paint, and 1/3 light grey paint. The same type and size of targets are disposed in the upper and in the lower part of the figure at different distances (resp., 40 m and 4 m) respect to the class 4 scenario. Green panels are clearly distinguished as little dark points respect to vegetation, because of their lower radiance in the SWIR bands (at 1.25 and 1.6 micron). Comparison between green grass (upper black curve in the SWIR bands) and light green paint panel (lower red curve in the SWIR bands) (b).

(a)

(b)

(c)

(d)

(e)

Figure 12 

SIMGA false colour radiance images on panels over vegetation (S. Rossore park): VNIR false colour (a), SWIR false colour (b). Comparison between reflectance of green panel (lower black curve at 1.6 micron) and green grass (upper green curve at 1.6 micron) (c), photo of green panel over grass (d), and comparison between SIMGA radiance acquired on green grass (upper green curve at 1.6 micron) and green panel (lower black curve at 1.6 micron) (e). Panels are clearly distinguished in the SWIR bands, while the contrast respect to vegetation in VIS bands is negligible.

3.5. Underwater Submerged Targets

Another dual use capability regarding the discrimination of underwater submerged targets was tested by means of SIMGA image simulations and verified with overflights in a controlled zone.

In order to test the detection capabilities of small grey panels under water a direct bathymetric model has been developed to simulate the total reflectance of shallow waters on the basis of chlorophyll, sediment, and yellow substance content, the bottom and panel reflectance and the water depth height.

The total reflectance (Figure 13) 𝑅(𝜆) has been calculated by means of the surface 𝑅_surface(𝜆) and subsurface 𝑅_sub-surface(𝜆) reflectance with the following relationship [26]: 𝑅(𝜆)=𝑅_surface(𝜆)+𝑅_sub-surface(𝜆)∗1−𝜌_𝐹𝜗^′𝑛², (12) With (i)𝜌(𝜃^′) the Fresnel reflectivity at the interface air-water which takes into account the reflection of the subsurface radiance into the water (~0.021),(ii)𝜃^′ is the incident angle of the radiation coming from below the water which generates a refraction in air at the angle 𝜃 in the observation direction [𝑛∗sin(𝜃^′)=sin(𝜃) with 𝑛=1.8368 the water refraction index],(iii)𝑅_surface reflectance depends from surface roughness and foam, but in this analysis has been taken constant and equal in first approximation to 0.021.

Figure 13 

Two fluxes direct bathymetry reflectance model (𝑅surf = surface reflectance due to roughness and foam, 𝑅bulk = Sub-surface water reflectance due to particulate concentration 𝐶 = Chlorophyll, 𝑆 = sediment, 𝑌 = yellow substance, and 𝑅bottom = reflectance of bottom.

The subsurface reflectance is obtained by means of a two-flux algorithm (Figure 13) which yields the following analytical relationship [27] for a water layer of uniform optical properties and thickness 𝐻 (m), above a reflecting bottom with reflectance 𝜌_bottom:𝑅(𝜆)=𝑅⁰(𝜆)+𝜌_bottom(𝜆)−𝑅⁰(𝜆)⋅1−𝑅⁰(𝜆)²𝑅⁰(𝜆)⋅𝜌_bottom(𝜆)−𝑅⁰(𝜆)+1−𝑅⁰(𝜆)𝜌_bottom(𝜆)⋅𝑒^{2𝐻⋅(𝐴(𝜆)2+2𝐴(𝜆)𝐵(𝜆))}, (13)where(i)𝑅⁰(𝜆)=𝐵(𝜆)/(𝐴(𝜆)+𝐵(𝜆)+𝐴(𝜆)²+2𝐴(𝜆)𝐵(𝜆)),(ii)𝐴(𝜆) and 𝐵(𝜆) represent, respectively, the total absorption and scattering coefficients, including that of water, chlorophyll, and yellow substance,(iii)the range of validity is 𝐵/(𝐴+𝐵)≤0.6.

The total absorption and backscattering coefficients are calculated from a three component water colour model ([28, 29]), which has been adapted for class 1 and class 2 waters [30].

In this model the total absorption and backscattering coefficients (m⁻¹) are obtained as a linear combination of that of water, chlorophyll, sediment, and yellow substance with the following relationships: 𝐵(𝜆)=0.5⋅𝐵_𝑤(𝜆)+0.005⋅𝐵_𝑐(𝜆)+0.015⋅𝐵_𝑠(𝜆),𝐴(𝜆)=𝐴_𝑤(𝜆)+𝑎_𝑐(𝜆)⋅𝐶+𝑎_𝑠(𝜆)⋅𝑆+𝐴_𝑌(𝜆), (14) where(i)the suffix 𝑤, 𝑐, 𝑠, 𝑌 means, respectively, water, chlorophyll, sediment, and yellow substance,(ii)𝐶, 𝑆, and 𝑌 represent the chlorophyll, sediment, and yellow substance content (𝐶 in mg/m³, 𝑆, and 𝑌 in g/m³),(iii)𝑎_𝑐,𝑎_𝑠,𝑎_𝑌, respectively, the chlorophyll, sediment and yellow substance specific absorption coefficients (m²/mg and m²/g, resp.) shown in Figure 14 [30](iv)𝐵_𝑤(𝜆)=0.002m⁻¹ (𝜆(nm)/550)^−4.3,(v)𝐵_𝑠(𝜆)=1(m²/mg)∗𝑆(g/m³)∗[(𝜆(nm)/550)⁻¹], (vi)𝐵_𝑐(𝜆)=0.12𝐶^0.63 [𝑎_𝑐(550)/𝑎_𝑐(𝜆)](m⁻¹),(vii)𝑎_𝑐(550)=0.0189m²/mg if 𝐶<1 mg/m³,(viii)𝑎_𝑐(550)=0.00486 m²/mg if 𝐶>1 mg/m³,(ix)𝐴_𝑌(𝜆)=𝐴_𝑌(443nm)∗exp[−0.014∗(𝜆(nm)−443)](m⁻¹),(x)𝐴_𝑌(443nm)=(0.12m⁻¹)∗𝑌, (xi)𝑎_𝑠(443nm)=0.034 m²/g,(xii)𝑎_𝑐(443nm)=0.07m²/mg if 𝐶<1 mg/m³,(xiii)𝑎_𝑐(443nm)=0.018 m²/mg if 𝐶>1 mg/m³.

Figure 14 

Specific absorption spectra of chlorophyll a, suspended sediment, and yellow substance normalized to the value at the 440 nm wavelength (left scale). Absorption spectrum of water is in m−1 (right scale) [12].

The concentration of the three water components can be divided in(i)completely correlated type 1 waters characterised by a rather stable correlation between optically active substances, with phytoplankton concentration as dominant, in this case the yellow substance backscattering and sediment absorption coefficients has been considered as related to the chlorophyll with the following relationships [12, 29]: 𝐵𝑠(550nm)=0.2⋅𝐶^0.63⋅m⁻¹,𝐴_𝑌()=443nm0.2𝐴0.8^𝑤(443nm)+𝑎_𝑐(443nm)⋅𝐶+𝑎_𝑠(),𝑆443nm⋅𝑆g/m³=10^−0.252⋅𝐶^0.57mg/m³, (15)(ii)completely uncorrelated coastal type 2 waters, with no correlation between the three water components, when high concentration of sediments and yellow substances exist,(iii)partially correlated coastal type 2 waters for which it is possible to retrieve a partial correlation between the three water components [30]. Examples are given by the following relationships:(1)Gulf of Naples [31] (𝑆)(𝐶)𝐴Log=−0.25+0.57⋅LogLog_𝑌(440)=−1.20+0.47⋅Log(𝐶), (16)(2)Northern basin of the Adriatic Sea [32] (𝑆)(𝐶)𝐴Log=−0.026+0.59⋅LogLog_𝑌(440)=−1.28+0.38⋅Log(𝐶), (17)(3)Tirrenian Sea near Migliarino-S. Rossore (Tuscany) [33] Log(𝐶)=Log(0.0206)+2.0615⋅Log(𝑆), (18)A reflectance simulation of S. Rossore waters at 1 m of bottom depth obtained by means of the two flux model is displayed in Figure 15. The Total Reflectance represents simulations of correlated waters type 2 model (𝐶=3 mg/m³, 𝑆=1.8 g/m³, 𝑌=0.6 g/m³), with the S. Rossore (I) bottom sand and grey panels reflectance measurements performed with a Field Spec portable spectrometer.

(a)

(b)

Figure 15 

Total Reflectance curves (green curves in the middle part of Figures 15(a) and 15(b)) from the two fluxes direct bathymetry reflectance model obtained from simulations with different type of S. Rossore bottom targets: grey panel (upper maroon part of Figure 15(a) (a) sand (upper blue part of Figure 15(b)) (b). The simulation is performed at a bottom depth (𝐻) of 1 m with partially correlated waters type 2 model (𝐶=3 mg/m3, 𝑆=1.8 g/m3, 𝑌=0.6 g/m3) (lower pink reflectance curves in Figures 15(a) and 15(b)).

Finally simulated SIMGA reflectance and radiance images of a marine environment (sand, waters with sand bottom at 2 m and 8 m of depths, panels of 1 m*1 m and 2 m*2 m under water at 1 m and 0.2 m of depth) have been performed (Figure 16). SIMGA radiance has been simulated at 1.5 km airborne altitude with MODTRAN code and at the SIMGA spatial (1 m for VIS and 2 m for SWIR bands) and spectral resolution (2.4 nm for VIS and 10.8 nm for SWIR). The simulated SIMGA reflectance image has been obtained through the inversion of MODTRAN parameters used for the direct simulation and results show that all grey panels (both at 0.2 m or 1 m depth) can be clearly distinguished both in low (2 m) and high (8 m) depth waters (Figure 16).

Figure 16 

Simulated radiance/reflectance image of marine environment (sand, waters with sand bottom at 2 m & 8 m depths, panels of 1 m * 1 m & 2 m ∗ 2 m under water at 2 m & 0.2 m depth, chlorophyll of 3 mg/m3, sediment of 1.8 g/m3, and yellow substance of 0.6 g/m3) obtained with SIMGA model at 1.5 km altitude (spatial resolution 1 m VIS, 2 m SWIR, spectral resolution 2.4 nm VIS, 10.8 nm SWIR). Simulated synthetic reflectance (waters with sand bottom at 2 m or 8 m of depth, panels under water at 0.2 m and 1 m of depth, beach surface sand) from two-fluxes model (a), radiances at 1.5 Km altitude (b), surface reflectance (also beach sand and concrete pier are represented in the upper part of Figures 16(b) and 16(d)) after atmospheric inversion algorithm based on Modtran code with the same parameters as the direct simulation (d), and RGB synthetic radiance image representing the simulated scenario (c).

This result was validated (Figure 17) by means of SIMGA overflights on the Morto mouth river (S. Rossore park in Tuscany, I), where two different grey panels were submerged. The two panels are clearly detectable in the visible part of the spectrum, so demonstrating the capability of the SIMGA hyperspectral instrument for this kind of application.

Figure 17 

SIMGA RGB radiance images ((a), (b), (c)) on panels under water (S. Rossore park, Fiume Morto river), comparison between radiance measurements of grey panel under water and water (d).

4. Conclusions

An end-to-end software tool (SG_SIM) for the simulation of airborne/satellite optical sensors images has been implemented in ENVI-IDL environment. Input images can be either high resolution airborne or synthetic data. The simulator features three separate modules: the reflectance scenario, which generates a desired reflectance image with spectral mixtures, the atmosphere module, which converts the input reflectance map into the at-sensor radiance image, and the instrument module, which simulates the main degradations introduced by the instrument (ISR, MTF, ICR and noise). As other end-to-end simulators the SG_SIM Simulator integrates a complete atmospheric radiative transfer modelling which could easily refined through the implementation of most MODTRAN options and it includes all main functions and features necessary for a complete hyperspectral image simulation such as ISR&MTF, ICR and noise sources. Compared to the other simulators (e.g., SENSOR, [4]), SG_SIM allows also the control of spectral mixing and the generation of synthetic scenario, but is lacking of MWIR/LWIR spectral bands, 3D reflectance simulation, and DEM ray-tracing functions as included in CAMEO [5]. The implementation and further development of the SG_SIM approach was boosted significantly by the Selex Galileo S.p.A. airborne imaging system SIMGA and by other phase 0/A studies carried out for preliminary evaluations of image quality and product accuracy from new classes of space-borne optical sensors. The validation of the simulator is reported in [16], whereas in this paper the simulator’s theoretical basis and some simulation examples have been described. For the simulated cases the following results can be outlined:(i)the 3D representation of the SRF allows the visual inspection of the spatial pixel response for image quality analysis,(ii)the potentials of the simulator for the HYPSEO Panchromatic camera trade-off analysis between project parameters (pupil diameter, optics degradations, detector noise, etc.) and system performances (SNR, spatial resolution, etc.) have been demonstrated by using as inputs synthetic bars and an IKONOS images with different radiance and surface albedo levels,(iii)the potentials of the HYPSEO hyperspectral camera for vegetation mapping has been demonstrated on the basis of a MIVIS airborne scene rescaled at satellite level and ground truth data,(iv)the potentials for detection of camuffled targets in a rural background has been demonstrated in the SWIR bands by means of a simulation of a synthetic scenario with green panels at different size,(v)the potentials for the identification of submerged targets in the visible spectral range at airborne level (1 m of spatial resolution) have been demonstrated by means of the simulation of a synthetic scenario with submerged grey panels and the implementation of a direct bathymetric-water color model to generate surface reflectance as input to the scene simulator,(vi)real airborne data on submerged and camuffled targets have confirmed the results from simulations performed before the flight campaign,

These results demonstrate the potentials of the proposed simplified end-to-end simulator as a preliminary aid tool (during phase 0/A) for the dimensioning of new optical instruments to trace the link between user and instrument requirements.

5. Annex

The formulation of the following MTF components implemented in the SG_SIM model are described in the following paragraphs:(i)Detector pixel size(ii)Detector cross talk(iii) CCD detector charge transfer(iv) Image motion during integration time(v)Electronic filtering(vi) Focal plane random jitter during integration time(vii)Optics diffraction and aberrations.

5.1. Detector Pixel Size

The finite size of detector pixel permits the spatial integration of the signal coming from a finite spatial region on ground and this introduces a sort of degradation of the original high resolution image. This effect is analogous to a spatial filter windowing with a rect function which is 1 within a certain spatial rectangular domain and 0 outside: PSF_pixel𝑥,𝑦=rectΔ𝑥,Δ𝑦=rect_𝑥(Δ𝑥)⋅rect_𝑦Δ𝑦. (19) The transfer function of this function obtained by its Fourier Transform is represented with the following relationship: MTF_pixel𝑘_𝑥,𝑘_𝑦=sin𝑐_𝑥𝜋⋅𝑘_𝑥⋅Δ𝑥⋅sin𝑐_𝑦𝜋⋅𝑘_𝑦⋅Δ𝑦, (20) where(i)𝑘_𝑥,𝑘_𝑦 are the spatial frequencies along the 𝑥 and 𝑦 direction,(ii)Δ𝑥, Δ𝑦 are the detector size along 𝑥 and 𝑦 directions,(iii)the sinc functions are expressed from the following relationships: sin𝑐_𝑥𝜋⋅𝑘_𝑥⋅Δ𝑥=sin𝜋⋅𝑘_𝑥⋅Δ𝑥𝜋⋅𝑘_𝑥,⋅Δ𝑥sin𝑐_𝑦𝜋⋅𝑘_𝑦⋅Δ𝑦=sin𝜋⋅𝑘_𝑦⋅Δ𝑦𝜋⋅𝑘_𝑦.⋅Δ𝑦 (21)

5.2. Detector Cross Talk

The detector cross talk between two successive pixels is taken into account as a first approximation by assuming a trapezoidal spatial windowing filter, instead of a rectangular one, which can be obtained by means of a convolution between two rect functions, one representing the detector size (Δ𝑥,Δ𝑦) and the other representing the cross talk size (𝛿𝑥,𝛿𝑦) between two successive pixels (𝛿𝑥<Δ𝑥,𝛿𝑦<Δ𝑦): PSF_pixel−2=𝑥,𝑦=rectΔ𝑥,Δ𝑦⊗rect𝛿𝑥,𝛿𝑦rect_𝑥(Δ𝑥)⊗rect_𝑥(⋅𝛿𝑥)rect_𝑦Δ𝑦⊗rect_𝑦𝛿𝑦. (22) The Transfer function is obtained from the following relationship: MTF_pixel−2𝑘_𝑥,𝑘_𝑦=sin𝑐_𝑥𝜋⋅𝑘_𝑥⋅Δ𝑥⋅sin𝑐_𝑥𝜋⋅𝑘_𝑥⋅⋅𝛿𝑥sin𝑐_𝑦𝜋⋅𝑘_𝑦⋅Δ𝑦⋅sin𝑐_𝑦𝜋⋅𝑘_𝑦⋅𝛿𝑦. (23)

5.3. CCD Detector Charge Transfer

For a CCD (Charge Capacitance Device) detector the reading of electrons acquired in each pixel of the matrix is performed by mean of a charge transfer from a pixel to the other. In this way the total transfer efficiency is related to the pixel-to-pixel Charge Transfer Efficiency (CTE) and the total number of transfers 𝑁_trans. As the CTE is not 1 some losses are present at the end of the reading time which implies a reduction of image contrast and then an MTF less than 1 in the direction of output register: MTF_CCD-CTE𝑘,𝑘_Nyq=𝑒^{−𝑁trans⋅(1−CTE)⋅[1−cos(𝜋⋅𝑘/𝑘Nyq)]}. (24) With(i)CTE the pixel Charge Transfer Efficiency (greater then 99,99%),(ii)𝑁_trans the total number of charge transfers,(iii)𝑘 the spatial frequency and 𝑘_Nyq the Nyquist frequency of the system (half of the sampling frequency 𝑘_Nyq=1/(2𝑝), with 𝑝 the detector pitch).

5.4. Image Motion During Integration Time

The effect of temporal acquisition (integration time 𝑡_int greater than zero) during the image motion (with velocity 𝑣), which happens along the scan direction (the satellite velocity in a push broom system) introduces an image blur, which can be taken into account with a PSF similar to a rect function, which represents the temporal aperture along the motion: PSF_temporal(𝑥)=rect_𝑥𝑣⋅𝑡_int. (25) With the following MTF: MTF_along-motion𝑘_𝑥,𝑘_𝑦=sin𝑐_𝑥𝜋⋅𝑘_𝑥⋅𝑣⋅𝑡_int. (26) The worst case happens when the integration time is equal or larger than the dwell time, that is, the spatial displacement equivalent to a pixel size, while the best case (MTF≅1) is for a short integration time.

5.5. Electronic Filtering

An electronic system can introduce a temporal smoothing due to a finite frequency bandwidth and thus a reduced space/temporal response. This effect has been simulated by using a general formulation based on the following Butterworth filter response: MTF_filter𝑓=11+𝑓/𝑁_𝑓⋅𝑓_Nyq^2⋅𝑁𝑏. (27) With(i)𝑓 the frequency,(ii)𝑓_Nyq the Nyquist frequency,(iii)𝑁_𝑏 the order of Butterworth filter,(iv)𝑁_𝑓 the ratio between the 3 dB filter frequency and the Nyquist frequency 𝑁_𝑓=𝑓_3dB/𝑓_Nyq; this ratio should be between 2.2 and 3 for a good reproduction of a square wave.

The above equation correctly reproduces the behaviour of the classical low-pass filter for 𝑁_𝑏=1.

5.6. Focal Plane Random Jitter during Integration Time

For high frequencies random vibration of the focal plane a Gaussian spatial response (PSF) can be taken into account. The Fourier Transform of PSF is still a Gaussian function, representing the MTF, with the following relationship: MTF_{jitter-random}𝑓=𝑒^{−2⋅𝜋2⋅[jitter⋅𝑝⋅𝑓]2}. (28) With(i)𝑝 the detector pitch,(ii)jitter∗𝑝 the fraction of pixel representing the rms values of random fluctuations,(iii)𝑓 the spatial frequency at detector level.

5.7. Optics Diffraction and Aberrations

The MTF related to diffraction from optics has been evaluated by using the O’Neill formulas, valid for diffraction in presence of a telescope with central obscuration.

The following formulation for MTF diffraction term is used [34]: MTF_diffraction𝑓=(𝐴+𝐵+𝐶)1−𝜂², (29a) with 𝜂,𝑓, and other parameters defined as follows:(i)𝜂 = obscuration factor = ratio between the obscuration diameter and the pupil diameter,(ii)𝑓 = spatial frequency at detector level (cm⁻¹),(iii)𝜔 = 𝑓/𝑓cut-off, (iv)𝑓 cut-off = optics cut-off (cm⁻¹) at detector level = 1/(𝜆*𝑓-number),(v)𝜆 = wavelength;(vi)𝑓-number = 𝐹/𝐷, ratio between the focal length 𝐹 and the pupil diameter 𝐷.

The 𝐴, 𝐵, and 𝐶 parameters are defined from the following relationships: 2𝜔≤1,𝐴=𝜋⋅arccos(𝜔)−𝜔√1−𝜔²𝜔𝜔>1,𝐴=0,𝜂=0,𝑡=0,𝜙=0,𝐵=0;𝜂≠0,𝑡=𝜂,𝜙=arccos1+𝜂²−4𝜔²2𝜂𝑡>1,𝐵=0𝑡≤1,𝐵=2𝜂²arccos(𝑡)−𝑡√1−𝑡²𝜋−2𝜂²𝜔≤1−𝜂2,𝐶=−2𝜂²𝜔<1+𝜂2,𝐶=2𝜂sin𝜙+1+𝜂²𝜙−21−𝜂²arctan1+𝜂𝜙1−𝜂tan2⋅1𝜋−2𝜂²,𝜔≥1+𝜂2,𝐶=0. (29b)(29c)(29d) The above MTF formulation for optics diffraction can be simplified to the following well known diffraction relationship in absence of central obscuration, which is zero for 𝑓≥𝑓_cut-off∶𝑓≤𝑓_cut-off,=2MTF𝜋⋅𝑓arccos𝑓_cut-off−𝑓𝑓_cut-off𝑓1−𝑓_cut-off²,𝑓>𝑓_cut-off𝑓,MTF=0,_cut-off=1𝜆⋅𝑓_number. (29e)Regarding possible optics aberrations the model takes into account, as a first approximation, the following exponential fitting function: MTF_aberration𝑓=𝑒^{−𝑘⋅(𝑓/𝑓cut-off)𝑥}, (30) with 𝑘 and 𝑥 representing two empirical parameters used to approximate all optics degradation effects.

Acknowledgments

The authors wish to thank L. Tommasi of Selex Galileo, R. Bonsignori of Eumetsat (formerly with Selex Galileo) and F. Pecchioni of University of Florence for the useful technical discussions and contributions. The authors also wish to acknowledge the ASI HYPSEO Phase A/B study under ASI/CSM/vdc/299/00 Contract for instrument specification and modelling.