Fast yet accurate computation of radiances in shortwave infrared satellite remote sensing channels

Nan Chen; Wei Li; Tonomori Tanikawa; Masahiro Hori; Rigen Shimada; Teruo Aoki; Knut Stamnes

doi:10.1364/OE.25.00A649

1. Introduction

Radiative transfer (RT) simulations play an important role in the development of satellite remote sensing algorithms. Sensors like the Moderate-Resolution Imaging Spectroradiometer (MODIS), the Medium Resolution Imaging Spectrometer (MERIS), the Visible Infrared Imaging Radiometer Suite (VIIRS), and the Second Generation GLobal Imager (SGLI) are designed to collect multi-channel measurements in order to obtain information about the Earth in visible (VIS), near infrared (NIR), shortwave infrared (SWIR), and thermal infrared (TIR) parts of the spectrum. These measurements are used to infer the physical and chemical composition of the Earth’s atmosphere and surface [1–4]. Many inversion algorithms rely on simulations of the measured signals. In order to accurately simulate the signals received by a satellite sensor the radiative transfer of light in the Earth’s coupled atmosphere–surface system must be considered. Absorption by atmospheric gases represents a challenge due to rapid and erratic variations of gaseous absorption with wavelength. The SWIR wavelength range is of great importance to the remote sensing of cloud and snow propeties. Table 1 shows available SWIR channels for three different sensors. The lower panel of Fig. 1 shows the atmospheric transmittance in the shortwave infrared (SWIR) spectral range simulated with 1 cm⁻¹ spectral resolution. As shown in the top panel of Fig. 1, for MODIS, SGLI, and VIIRS (not shown, see Table 1), important SWIR channels for cloud/snow parameter retrieval have been selected to lie in atmospheric “window” regions, where the atmosphere is relatively transparent, to minimize attenuation by atmospheric gases. It can be seen that even in these window regions, there are many absorption lines due to different atmospheric gases as well as the continuum absorption by water vapor. A large number of monochromatic radiative transfer calculations may be needed to accurately account for the spectral variation in gaseous absorption. This circumstance makes accurate simulations very time consuming. To address this problem, much effort has been devoted to reduce the computation time while maintaining required accuracy. For this purpose, many methods have been developed, such as the Correlated k-distribution (Ckd), Exponential Sum Fitting of Transmittance (ESFT) [5], Optimal Spectral Sampling (OSS) [6], Principal Component Analysis (PCA) [7], and Neural Networks [8]. These methods are often “trained” by using a state-of-the-art line-by-line model (such as LBLRTM [9]) to deal with absorption by atmospheric gases and they can typically achieve better than 1% accuracy for a given satellite channel. The training process is usually quite time consuming, because in addition to the line-by-line computations, one must determine the required coefficients through numerical fitting. Also, the process may have to be repeated or adjusted for different atmospheric conditions.

Table 1. SWIR channel specifications of SGLI, MODIS and VIIRS

View Table | View all tables in this article

Fig. 1 The extraterrestrial spectral irradiance at the top of atmosphere (middle panel) and atmospheric transmittance (bottom panel) in the SWIR spectral region. MODIS and SGLI SWIR channel response functions are plotted in the top panel.

Download Full Size | PDF

On the other hand, for remote sensing of cloud or snow/ice properties one needs comprehensive radiative transfer simulations including not only atmospheric gases but also aerosol, cloud, and snow/ice materials. Due to the relatively “smooth” spectral variation of aerosol, cloud, and snow/ice IOPs compared to that of atmospheric gases, the treatment of aerosols, clouds, and snow/ice has received less attention in the past, and it is customary to average their IOPs to speed up radiative transfer calculations for narrow-band remote sensing applications (see, e.g. Baum et al. 2005 [10]). However, during the development of a cloud mask algorithm for the cryosphere mission of JAXA’s GCOM-C1/SGLI project [11], we noticed that there is significant spectral variation in the reflectivity of clouds and snow/ice in the SWIR wavelength range. This variation is related to the spectral change in refractive indices of water and ice. As shown in Fig. 2, the refractive indices of water and ice vary considerably within the SGLI SW03 channel (λ_ctr = 1.63 µm) and MODIS Ch07 (λ_ctr = 2.13 µm). Hence, there is a need to investigate how this spectral variation impacts our ability to create useful mean IOPs for these channels that can be used to reduce the number of quasi-monochromatic RT simulations. To address this problem, we performed detailed RT simulations to quantify how accurate top-of-the-atmosphere (TOA) radiances can be computed in a cost-effective manner.

Fig. 2 Refractive index of ice and liquid water in the SWIR spectral range. Solid lines: real parts; Broken lines: imaginary parts. The refractive indices of ice and liquid water are from Warren and Brandt 2008 [12] and Segelstein 1981 [13], respectively. Response functions of SGLI (blue) and MODIS (green) channels are also plotted.

Download Full Size | PDF

1.1. Radiative transfer simulations of satellite sensor signals

In order to simulate the light signal received by a satellite instrument, we need to solve the radiative transfer equation (RTE) pertinent for light propagation in the coupled atmosphere-surface system. The diffuse radiance I(τ, θ, ϕ) at wavelength λ is found by solving the following RTE:

\begin{array}{l} μ \frac{d I (τ, θ, ϕ)}{d τ} = I (τ, θ, ϕ) - \frac{ϖ (τ) F^{s} e^{- τ / μ_{0}}}{4 π} p (τ, θ^{'}, ϕ^{'}; θ_{0}, ϕ_{0}) \\ - \frac{ϖ (τ)}{4 π} \int_{0}^{2 π} d ϕ^{'} \int_{- 1}^{1} d μ^{'} p (τ, θ^{'}, ϕ^{'}; θ, ϕ) I (τ, θ^{'}, ϕ^{'}) . \end{array}

Here F^s is the TOA solar irradiance (normal to the beam), while the differential optical depth dτ = −(α + β)dz, the single scattering albedo ϖ = β/(α + β) = β/γ, and the scattering phase function p(τ, θ′, ϕ′; θ, ϕ) are the inherent optical properties (IOPs) of the scattering/absorbing medium. Note that we have used the Greek letters α, β, and γ = α + β to denote the absorption, scattering, and extinction coefficients, respectively. θ₀ and ϕ₀ represent solar zenith and azimuth angles, µ₀ = cos θ₀; θ′ and ϕ′ are sensor zenith and azimuth angles prior to a scattering event, and θ and ϕ the correspong angles after the scattering event, and µ = cos θ.

1.1.1. Atmospheric IOPs

We used the U.S. Standard atmospheric constituent profiles from Anderson et al. 1986 [14] divided into 14 layers to provide input to a band model based on MODTRAN [see for example, Thomas and Stamnes 1999 [15] for details] to generate absorption coefficients and optical depths due to atmospheric trace gases including H₂O, CO₂, O₃, CH₄, and NO₂. Layering is needed to resolve the vertical variation in the IOPs, and experience has shown that 14 layers is sufficient for this purpose. Molecular (Rayleigh) scattering optical depths are computed from the Rayleigh scattering cross section [15] multiplied by the bulk density of air available from Anderson et al. [14]. For aerosols and clouds (which can be assumed to consist of a polydispersion of spherical or nonspherical particles) the monochromatic IOPs can be calculated from Mie-Debye theory (for spherical particles) or T-matrix/ray-tracing techniques (for nonspherical particles) [16,17] given the size (and shape for nonspherical particles) distribution and the wavelength-dependent refractive index.

1.1.2. Surface IOPs

In this paper, we focus especially on snow surfaces and we consider a layer of snow sufficiently thick that it can be considered to be semi-infinite. Hence, we do not need to specify the surface below the snow. Snow particles were assumed to be ice spheres with the refractive index of ice obtained from Warren and Brandt 2008 [12]. The monochromatic IOPs can be calculated from Mie-Debye theory once the size distribution is specified or from a parameterization in terms of effective snow grain size [18]. For comparison with the snow surface, Lambertian surfaces with a fixed (spectrally independent) albedo will also be considered in the following.

1.1.3. Computation of the sensor signal

The TOA radiance simulation was carried out using the atmospheric and surface IOPs discussed above and DISORT3 [19] (the latest version of the DISORT code, Stamnes et al. 1988 [20]) was used to carry out the computations required in this paper based on input IOPs obtained as discussed in Sections 1.1.1 and 1.1.2. In order to simulate the signal received by an actual instrument channel with bandwidth Δλ, the radiance I(τ, θ, ϕ) must be multiplied by the response function R_s(λ) of that channel and integrated over the spectral range Δλ as follows:

\begin{array}{r} 〈 I (τ, θ, ϕ) 〉 \equiv I_{Δ λ} (τ, θ, ϕ) = \frac{\int_{Δ λ} R_{s} (λ) I (τ, θ, ϕ) d λ}{\int_{Δ λ} R_{s} (λ) d λ} \\ = \int_{Δ λ} {\tilde{R}}_{s} (λ) I (τ, λ, θ, ϕ) d λ \end{array}

where I_Δ_λ (τ, θ, ϕ) is the band-averaged radiance, R_s(λ) is the sensor response function (SRF) of the receiving instrument and

{\tilde{R}}_{s} (λ)

is the normalized SRF given by

{\tilde{R}}_{s} (λ) = \frac{R_{s} (λ)}{\int R_{s} (λ) d λ}

. In practice, this integration over wavelength is approximated as:

I_{Δ λ} (τ, θ, ϕ) \approx \frac{\sum_{i = 1}^{K} w_{i} R_{s} (λ_{i}) I_{i} (τ, θ, ϕ)}{\sum_{i = 1}^{K} w_{i} R_{s} (λ_{i})} = \sum_{i = 1}^{K} w_{i} {\tilde{R}}_{s} (λ_{i}) I_{i} (τ, θ, ϕ)

where the w_i’s are the weights needed in the numerical integration, and I_i (τ, θ, ϕ) ≡ I(λ_i, τ, θ, ϕ) is the TOA radiance at wavelength λ_i.

A large number of spectral points λ_i may be needed to resolve the rapidly varying optical properties of atmospheric gases (see Fig. 1). The more smoothly varying aerosol, cloud, and surface IOPs (see Fig. 2) must also be carefully considered as demonstrated below. These atmospheric and surface IOP variations make accurate evaluations of Eq. (2) or its approximate, discrete version [Eq. (3)] very time-consuming, because we need to solve the RTE [Eq. (1)] to obtain I_i (τ, θ, ϕ) at each spectral point i = 1, … K in Eq. (3). Hence, if a large value of K is required to get accurate results, the computations will be very expensive.

In satellite remote sensing, the quantity of interest is the light reflected at TOA (τ = 0) from the atmosphere-surface system. The TOA reflectance at wavelength λ, is usually defined as [I⁺(0, θ, ϕ) ≡ I⁺(τ = 0, θ, ϕ)]:

ρ (λ, θ, ϕ) = \frac{π I^{+} (0, θ, ϕ)}{μ_{0} F^{s} (λ)}

where I⁺(0, θ, ϕ) represents the upward radiance at the TOA. For an actual instrument, we may define the band-averaged reflectance similarly as:

\begin{matrix} ρ_{Δ λ} (θ, ϕ) \equiv 〈 ρ (θ, ϕ) 〉 = \frac{π I_{Δ λ}^{+} (0, θ, ϕ)}{μ_{0} F_{Δ λ}^{s}} = \frac{π \int_{Δ λ} {\tilde{R}}_{s} (λ) I^{+} (0, θ, ϕ) d λ}{μ_{0} \int_{Δ λ} {\tilde{R}}_{s} (λ) F^{s} (λ) d λ} \\ \approx \frac{π \sum_{i = 1}^{K} w_{i} {\tilde{R}}_{s} (λ_{i}) I_{i}^{+} (0, θ, ϕ)}{μ_{0} \sum_{i = 1}^{K} w_{i} {\tilde{R}}_{s} (λ_{i}) F^{s} (λ_{i})} \end{matrix}

2. Mean IOP methods

2.1. Chandrasekhar mean IOP

Since accurate radiative transfer computations are expensive, instead of evaluating Eq. (3) monochromatically one may average the IOPs of the atmosphere and surface media and perform radiative transfer calculation using the mean IOPs for a given satellite channel, which will greatly reduce the computational time. One commonly used approach is the Chandrasekhar mean IOP method [15] which was used by Baum et al. 2005 [10] to generate MODIS cloud products. For planetary media the Chandrasekhar mean IOPs can be obtained as follows:

γ_{Δ λ} = \frac{\int_{Δ λ} {\tilde{R}}_{s} (λ) F^{s} (λ) γ (λ) d λ}{\int_{Δ λ} {\tilde{R}}_{s} (λ) F^{s} (λ) d λ}

ϖ_{Δ λ} = \frac{β_{Δ λ}}{γ_{Δ λ}} = \frac{\int_{Δ λ} {\tilde{R}}_{s} (λ) F^{s} (λ) β (λ) d λ}{\int_{Δ λ} {\tilde{R}}_{s} (λ) F^{s} (λ) γ (λ) d λ} = \frac{\int_{Δ λ} {\tilde{R}}_{s} (λ) F^{s} (λ) γ (λ) ϖ (λ) d λ}{\int_{Δ λ} {\tilde{R}}_{s} (λ) F^{s} (λ) γ (λ) d λ}

g_{Δ λ} = \frac{\int_{Δ λ} {\tilde{R}}_{s} (λ) F^{s} (λ) β (λ) g (λ) d λ}{\int_{Δ λ} {\tilde{R}}_{s} (λ) F^{s} (λ) β (λ) d λ} = \frac{\int_{Δ λ} {\tilde{R}}_{s} (λ) F^{s} (λ) γ (λ) ϖ (λ) g (λ) d λ}{\int_{Δ λ} {\tilde{R}}_{s} (λ) F^{s} (λ) γ (λ) ϖ (λ) d λ}

where γ(λ), ϖ(λ) and g(λ), are the extinction coefficient, the single-scattering albedo, and the asymmetry factor at wavelength λ. For atmospheric gases the cumulative optical thickness due to gaseous absorption from TOA to atmospheric layer ℓ can be written as:

τ_{Δ λ, ℓ} = - \ln [T_{Δ λ, ℓ}]

where

T_{Δ λ, ℓ} = \frac{F_{Δ λ, ℓ}}{F_{Δ λ}^{s}} = \frac{\int_{Δ λ} {\tilde{R}}_{s} (λ) F^{s} (λ) T_{ℓ} (λ) d λ}{\int_{Δ λ} {\tilde{R}}_{s} (λ) F^{s} (λ) d λ}

is the mean transmission from the TOA to atmospheric layer ℓ. After calculating the mean molecular (Rayleigh) scattering optical depth τ_Ray,_ℓ in atmospheric layer ℓ (which is related to the pressure and the Rayleigh scattering cross section), the single-scattering albedo for atmospheric gases in atmospheric layer ℓ can be defined as:

ϖ_{Δ λ, ℓ} = \frac{τ_{Ray, ℓ}}{τ_{Δ λ, ℓ} + τ_{Ray, ℓ}} .

2.2. Refractive index mean IOP

For a narrowband satellite channel, one may define a band-averaged complex refractive index of a particle as:

n_{Δ λ, r} = \int_{Δ λ} {\tilde{R}}_{s} (λ) n_{r} (λ) d λ

and

n_{Δ λ, i} = \int_{Δ λ} {\tilde{R}}_{s} (λ) n_{i} (λ) d λ

where n_r(λ) and n_i(λ) are the real and imaginary parts of the particle’s refractive index at wavelength λ, respectively. Hence, only a single Mie/T-matrix/ray-tracing calculation is needed for one satellite channel by using the band-averaged refractive indices n_Δ_λ_,r and n_Δ_λ_,i in the IOP calculation. This approach can also reduce the computation time and we will refer to this method as the refractive index mean IOP method.

2.3. Central wavelength method

The “central wavelength” method is another way to obtain the effective IOPs for a satellite sensor channel. It has been applied to construct MODIS aerosol retrieval algorithms [3,4]. In this method, one calculates a band-weighted “central wavelength” CWL_Δ_λ for a satellite channel and the IOPs for aerosols are calculated at this representative wavelength. The central wavelenght method is a fast approach that works similarly to the refractive index mean method since only one IOP calculation is needed for each channel. We used the following equation to determine the central wavelength of each channel:

{CWL}_{Δ λ} = \int_{Δ λ} {\tilde{R}}_{s} (λ) F^{s} (λ) λ d λ .

The complex refractive indices of ice and water are then interpolated to the central wavelength of each satellite channel to calculate the IOPs of clouds and snow needed for that channel. The central wavelength values for the SGLI and MODIS channels used in our calculations are listed in Table 1.

3. Tests of the mean IOP methods

We are now in a position to test the accuracy of the three mean IOP methods described in Section 2. In order to isolate possible sources of error we will test the clear atmosphere, cloud, and surface (snow) parameterizations separately. In all cases, the benchmarks for the evaluation, will be the TOA reflectance computed using Eq. (5) with 1 cm⁻¹ spectral resolution, which requires about 200 to 1000 monochromatic calculations for a SWIR channel depending on its spectral location and width. The computations were carried out as described in Section 1.1 with atmospheric and surface properties adopted as specified in Sections 1.1.1 and 1.1.2.

Figure 3 shows the simulated reflectance of snow using Chandrasekhar mean IOPs for gas and snow and the relative error of this method compared to the benchmark for MODIS channel Ch07 (2.13 µm). The grain size of snow was assumed to be 200 µm and the solar zenith angle (SZA) to be 50°. We note that although the reflectance has a strong anisotropic feature, the relative error resulting from the use of mean IOPs to reduce the number of monochromatic computations compared to the benchmark is a very “flat” function (Fig. 3, right panel) that changes very slowly with respect to the viewing zenith angle (VZA) for VZAs smaller than 75°, and it has very weak dependence on relative azimuth angle. The rapid increase of the relative error for VZAs larger than 75° may be related to the very long optical path of the light at large viewing angles, which tends to magnify the error. We repeated the tests for different settings (SZA from 0 – 70°, snow surface with different snow grain sizes, with and without atmosphere) and obtained similar results. For MODIS and SGLI the largest VZAs are 65° and 50°, respectively. Hence, it will be unnecessary to simulate the radiance at very large VZAs. Therefore, we will only consider VZAs smaller than 65° and the error will be plotted at 0° relative azimuth angle using the same scale to facilitate cross-comparison between different cases.

Fig. 3 Left: Polar plot of the snow reflectance when using Chandrasekhar mean IOPs for gas and snow. Right: Percentage relative error compared to the benchmark.

Download Full Size | PDF

3.1. Black-sky overlying a snow surface

We first did a black-sky (no atmosphere) test to investigate the performance of mean IOP methods for different satellite channels in the SWIR wavelength range. Since there is no influence from the atmosphere one can more easily quantify how use of mean surface IOPs impacts the accuracy. We averaged the IOPs for snow using the Chandrasekhar mean as well as the refractive index mean methods described in Sections 2.1 and 2.2. We define the relative error as

RE = \frac{ρ_{Δ λ}^{approx} - ρ_{Δ λ}^{accur}}{ρ_{Δ λ}^{accur}}

where the superscripts “accur” and “approx” refer to the benchmark and approximate results, respectively. Figure 4 shows the relative error in the reflectance incurred by using these mean IOP methods compared to the benchmark computations for the black-sky case for MODIS and SGLI SWIR channels at SZA = 50°. It can be seen that the mean IOP methods (especially the Chandrasekhar mean) yield quite accurate results for MODIS Ch05 (1.24 µm) and Ch06 (1.64 µm) channels as well as SGLI SW01 (1.05 µm) and SW04 (2.20 µm) channels. These results imply that one can save a significant amount of computing time since only a single (quasi-monochromatic) radiative transfer calculation can provide almost the same results as the benchmark. However, the error incurred by use of the mean IOP methods becomes quite large in SGLI SW03 (1.63 µm) and MODIS Ch07 (2.13 µm) channels. Since MODIS Ch06/Ch07 and SGLI SW03/SW04 have similar central wavelengths (λ_ctr), why would the results be so different? To explain these results we note the following differences between these SWIR channels:

The SGLI SW03 (λ_ctr = 1.63 µm, FWHM = 200 nm) channel is considerably wider than the MODIS Ch06 (λ_ctr = 1.64 µm, FWHM = 24 nm) channel as can be seen in Figs. 1 and 2.
The central wavelength of the SGLI SW04 (λ_ctr = 2.20 µm) channel is shifted to a longer wavelength by about 70 nm compared to the MODIS Ch07 (λ_ctr = 2.13 µm) channel.

Fig. 4 Relative error [see Eq. (15)] of surface IOP averaging under black-sky conditions. Upper panels: SGLI SWIR channels. Lower panels: MODIS SWIR channels.

Download Full Size | PDF

From Fig. 2 we note that compared to other SWIR channels the imaginary part of the refractive index of ice changes significantly within the wavelength range of the SGLI SW03 (λ_ctr = 1.63 µm, FWHM = 200 nm) channel and the MODIS Ch07 (λ_ctr = 2.13 µm, FWHM = 50 nm) channel, which would imply a significant variation in snow IOPs across these channels and consequently an increased error when applying mean IOP methods to these channels.

3.2. Atmosphere over a Lambertian surface

Next we consider how the mean IOP methods perform when dealing with gas absorption which varies rapidly and erratically with wavelength in the SWIR channels compared to that of snow/ice and clouds. To investigate this issue we averaged the IOPs for atmospheric gases in each atmospheric layer ℓ using the Chandrasekhar mean [Eqs. (9) and (11)] and adopted a Lambertian surface with a fixed albedo of 0.1 to eliminate any spectral variation in the surface albedo. Figure 5 shows the relative error compared to the benchmark at SZA = 50°. Perhaps surprisingly, for atmospheric gases the Chandrasekhar mean IOP method performs quite well in the SWIR spectral region. For SGLI channel SW01 (1.05 µm) and MODIS channel Ch05 (1.24 µm), the error is well below 1% and for channels such as the SGLI SW04 (2.20 µm) and MODIS Ch07 (2.13 µm) with relatively more absorption, the error is about 1.5% at nadir (VZA = 0°). Generally, the Chandrasekhar mean IOP method tends to slightly overestimate the absorption by atmospheric gases since the radiance computed by this method is usually less than the benchmark value. This behavior is especially true for channels with a large number of absorption lines such as SGLI channel SW04 (2.20 µm) and MODIS channel Ch7 (2.13 µm). Although the treatment of gaseous absorption has traditionally been considered to be the most difficult and time-consuming part in radiative transfer computations due to the need for numerous “quasi-monochromatic” computations, our results show that the Chandrasekhar mean IOP method provides a computationally efficient alternative to CkD, ESFT, OSS, and PCA methods in the SWIR “window” channels.

Fig. 5 Relative error [see Eq. (15)] caused by gas IOP averaging for a clear atmosphere overlying a Lambertian surface with albedo = 0.1. Upper panels: SGLI SWIR channels. Lower panels: MODIS SWIR channels.

Download Full Size | PDF

4. New sub-band IOP method

From the tests in the previous section, we conclude that the mean IOP methods do not perform well enough when applied to a channel across which the surface IOPs change significantly, and that the main source of error in the simulation of satellite radiances may be due to the smooth, but considerable spectral variation in the surface IOPs. However, the Chandrasekhar mean IOP for atmospheric gases is shown to handle the gaseous absorption in the SWIR channels quite accurately. Hence, they can be used to speed up the traditionally most time-consuming gaseous absorption part in radiative transfer calculations. Therefore, we designed a new sub-band IOP method that is fast while being able to perform well under a variety of surface IOP conditions. This new sub-band IOP method may be described as follows:

Divide the wavelength range of a satellite channel into several sub-bands and perform one radiative transfer calculation for each sub-band.
Sample the surface IOP and sensor response function at the center of each sub-band.
Use the Chandrasekhar mean to compute gas IOPs in each sub-band.
Use the mean solar irradiance in each sub-band as the solar input.

The idea of dividing the wavelength range into several sub-bands, each weighted by the sub-band solar irradiance, was used by Aoki et al. 2011 [21] to calculate broadband albedo. In our new method the total number of “quasi-monochromatic” radiative transfer calculations is equal to the number of sub-bands, and the TOA reflectance is computed using Eq. (5). If this method can be used to achieve accurate results using a small number of sub-bands, the computational burden will be significantly reduced. In the following subsections we will test this new sub-band IOP method in order to quantify its accuracy and explore how many sub-bands are needed for accurate radiance simulations in commonly used SWIR channels.

4.1. New sub-band IOP method applied to clear-sky cases

In this section we will test the performance of the new method for a clear-sky atmosphere overlying a snow surface. For comparison, we added two other configurations using the mean IOP methods for the snow surface. The test configurations are described as follows:

Method 1: Chandrasekhar mean IOP for both atmospheric gases and snow;
Method 2: Chandrasekhar mean IOP for atmospheric gases and refractive index mean IOP for snow;
Method 3: Chandrasekhar mean IOP for atmospheric gases and central wavelength IOP for snow;
New sub-band IOP method with a few sub-bands per satellite channel.

Figure 6 and Table 2 show the relative error of all methods we have tested compared to the benchmark results. The new sub-band IOP method is superior to Methods 1, 2 and 3 in all channels [most notably for SGLI channel SW03 (λ_ctr = 1.63 µm, FWHM = 200 nm) and MODIS channel Ch07 (λ_ctr = 2.13 µm, FWHM = 24 nm)]. The mean absolute error is reduced to about 1% or less. Tests carried out with different numbers of sub-bands (4, 8, and 16) revealed that use of 8 sub-bands is sufficient to obtain satisfactory results (error of about 1%) for SWIR channels considering the trade-off between speed and accuracy. As an example, for MODIS Ch07 channel the error of the nadir reflectance incurred by using 4, 8, and 16 sub-bands were 3.8%, 0.80%, and 0.66%, respectively. Also, it is interesting to note that the errors of the new method in Fig. 6 seem to be very close to the errors shown in Fig. 5, which is based on a Lambertian surface rather than snow. d

Fig. 6 Relative error [see Eq. (15)] of the sub-band method (blue curve with 8 sub-bands) for SGLI and MODIS SWIR channels. Clear-sky atmosphere overlying a snow surface with grain size 200 µm.

Download Full Size | PDF

Table 2. Relative errors (%) incurred by use of mean IOP methods as well as the new method with 8 sub-bands.

View Table | View all tables in this article

4.2. New sub-band IOP method applied to cloudy cases

In addition to the tests for clear-sky cases, we also tested the performance of the new method under cloudy-sky conditions. Since there could be many different combinations of clouds (liquid water/ice) and surface (lambertian/snow), we will only focus on some typical situations. Figure 7 shows the relative error of all methods we have tested compared to the benchmark results for a semi-transparent water cloud (optical depth τ = 0.5 at 550 nm) located at 2 km above a snow surface. It can be seen that the error for mean IOP methods become more dependent on the viewing zenith angle, which may due to the strong interaction between the cloud and the snow surface caused by multiple reflections between the cloud and the surface. Generally the error varies with cloud optical depth and solar/viewing geometries and the new sub-band IOP method performs consistently as shown in Figure 7.

Fig. 7 Relative error [see Eq. (15)] of the sub-band method (blue curve with 8 sub-bands) for SGLI and MODIS SWIR channels. Thin (τ = 0.5) liquid water cloud overlying a snow surface with grain size 200 µm.

Download Full Size | PDF

5. Further improvement

While the results in Fig. 6 and Table 2 show significant improvement of the new sub-band IOP method over the other methods, we note that some errors remain in MODIS Ch07 (λ_ctr = 2.13 µm, FWHM = 50 nm) and SGLI SW04 (λ_ctr = 2.20 µm, FWHM = 50 nm). Also, as mentioned before, the errors produced by the new method in Fig. 6 seem to be very close to the errors shown in Fig. 5, which is from use of the Chandrasekhar mean for gas IOP over a Lambertian surface with fixed albedo. This finding is very interesting because the two configurations (clear-sky atmosphere overlying a Lambertian surface with fixed albedo in Fig. 5, but a snow surface in Fig. 6) have very different TOA reflectances (snow yields a much more anisotropic reflectance than a Lambertian surface). To further explore this behavior, we did more tests for several solar illumination angles for these two configurations and the results are shown in Fig. 8. One can observe in the left panel that the error of our new sub-band IOP method increases with solar zenith angle. However, the error from the use of Chandrasekhar mean gas IOPs over a Lambertian surface with fixed albedo = 0.1 (right panel) has almost identical features. In fact, the error remains almost the same when we change the surface albedo from 0.1 to 0.5, which means that it is largely independent of the surface reflectance. This finding suggests that it may be possible to correct the TOA radiance for an arbitrary surface using the error resulting from the use of a Lambertian surface. We tested this possibility and our procedure is described as follows. First we did a benchmark calculation for a clear sky case with a Lambertian surface. Then we did a calculation with the same Lambertian surface but with Chandrasekhar mean gas IOPs, and we used Eq. (15) to quantify the error at each viewing zenith angle. Finally, since the error is nearly independent of the underlying surface condition (Lambertian versus snow) as shown in Fig. 8, we may invert Eq. (15) to obtain a “corrected” TOA reflectance for an arbitrary surface as follows:

ρ_{Δ λ}^{corr} = \frac{ρ_{Δ λ}^{approx}}{(1 + {RE}_{Lamb})}

where RE_Lamb is the error resulting from the use of Chandrasekhar mean gas IOP with a Lambertian surface. Figure 9 shows the relative error of the new method after correction for MODIS channel Ch07. It can be seen that the correction almost completely removes the error even at large solar zenith angles.

Fig. 8 Relative error [see Eq. (15)] for different solar zenith angles for MODIS channel Ch07. Left: Error from the new method (sub-band IOP method) for a snow surface (200 µm grain size). Right: Chandrasekhar mean gas IOP over a Lambertian surface (albedo = 0.1).

Download Full Size | PDF

Fig. 9 MODIS channel Ch07 relative errors [see Eq. (15)] for a clear-sky atmosphere overlying a snow surface using Method 3 before (black curve) and after (blue curve) correction for three different solar zenith angles.

Download Full Size | PDF

6. Conclusion

Detailed radiative transfer simulations were carried out to develop fast yet accurate ways to compute TOA radiances in satellite channels. While the approach is quite generic, and can, in principle, be applied to any spectral range in the solar spectrum, our focus in this paper is on the shortwave infrared (SWIR) range used for satellite remote sensing of atmospheric (cloud and aerosol) and surface (snow in particular) properties by existing sensors like MODIS and VIIRS, and future sensors like SGLI. By comparison with benchmark results computed at high spectral resolution, it was found that commonly used mean IOP methods (such as the Chandrasekhar mean) may yield unacceptably large errors when there is a significant spectral variation in the surface IOPs across a sensor channel. This finding suggests that one might need to consider absorption/scattering not only by atmospheric gases but also the underlying surface when selecting satellite channels and designing retrieval algorithms for specific targets such as ice/liquid water clouds, aerosols and snow/ice surfaces. On the other hand, for the SWIR channels of MODIS and SGLI use of Chandrasekhar mean IOPs for atmospheric gases yields good results. A new efficient approach to compute TOA radiances received by a satellite instrument has been established that employs the Chandrasekhar mean IOP for atmospheric gases, but samples the surface IOPs at a few wavelengths. It is found that 8 sub-bands may be sufficient to ensure accurate results for channels with significant spectral changes in surface IOPs. The error for solar zenith angles smaller than 70° is typically less than 1% for MODIS and SGLI SWIR channels. A special correction method was developed to remove the remaining error in MODIS Ch07 (λ_ctr = 2.13 µm) and SGLI SW04 (λ_ctr = 2.20 µm) channels. This correction, based on calculations using a Lambertian surface, leads to better than 0.1% accuracy for solar zenith angles smaller than 70°.

Funding

JAXA JX-PSPC-434863.

Acknowledgments

This work was conducted as a part of the GCOM-C1/SGLI algorithm development effort and was supported by the Japan Aerospace Exploration Agency (JAXA) under contract number JX-PSPC-434863.

References and links

1. D. K. Hall, G. A. Riggs, and J. S. Barton, “Algorithm theoretical basis document (ATBD) for the MODIS snow and sea ice-mapping algorithms,” NASAGSFC Technical Report, September (2001).

2. M. D. King, S.-c. Tsay, S. E. Platinick, M. Wang, and K.-N. Liou, “Cloud retrieval algorithms for MODIS: optical thickness, effective particle radius, and thermodynamic phase,” MODIS Algorithm Theoretical Basis Document, December (1997).

3. R. C. Levy, L. A. Remer, S. Mattoo, E. F. Vermote, and Y. J. Kaufman, “Second-generation operational algorithm : Retrieval of aerosol properties over land from inversion of Moderate Resolution Imaging Spectroradiometer spectral reflectance,” J. Geophys. Res. Atmos. 112, 1–21 (2007 [CrossRef] ).

4. R. C. Levy, S. Mattoo, L. A. Munchak, L. A. Remer, A. M. Sayer, F. Patadia, and N. C. Hsu, “The collection 6 MODIS aerosol products over land and ocean,” Atmos. Meas. Tech. 6(11), 2989–3034 (2013 [CrossRef] ).

5. W. J. Wiscombe and J. W. Evans, “Exponential-sum fitting of radiative transmission functions,” J. Comput. Phys. 24, 416–444 (1977 [CrossRef] ).

6. J. L. Moncet and S. A. Clough, “Accelerated monochromatic radiative transfer for scattering atmospheres: Application of a new model to spectral radiance observations,” J. Geophys. Res. Atmos. 102, 21853–21866 (1997 [CrossRef] ).

7. V. Natraj, X. Jiang, R. lie Shia, X. Huang, J. S. Margolis, and Y. L. Yung, “Application of principal component analysis to high spectral resolution radiative transfer: A case study of the band,” J. Quant. Spectrosc. Radiat. Transf. 95, 539–556 (2005 [CrossRef] ).

8. R. Key and A. J. Schweiger, “Tools for atmospheric radiative transfer: streamer and FluxNet,” Comput. Geosci. 24, 443–451 (1998 [CrossRef] ).

9. S. A. Clough and M. J. Iacono, “Line-by-line calculation of atmospheric fluxes and cooling rates: 2. application to carbon dioxide, ozone, methane, nitrous oxide and the halocarbons,” J. Geophys. Res. Atmos. 100, 16519–16535 (1995 [CrossRef] ).

10. B. A. Baum, P. Yang, A. J. Heymsfield, S. Platnick, M. D. King, Y.-X. Hu, and S. T. Bedka, “Bulk Scattering Properties for the remote sensing of ice clouds. Part II: Narrowband models,” J. Appl. Meteorol. 44, 1896–1911 (2005 [CrossRef] ).

11. N. Chen, W. Li, T. Tanikawa, M. Hori, T. Aoki, and K. Stamnes, “Cloud Mask over snow/ice covered areas for the GCOM-C1/SGLI cryosphere mission: validations over Greenland,” J. Geophys. Res. Atmos. 19(21), 12287–12300 (2014 [CrossRef] ).

12. S. G. Warren and R. E. Brandt, “Optical constants of ice from the ultraviolet to the microwave: A revised compilation,” J. Geophys. Res. 113, D14220 (2008 [CrossRef] ).

13. D. J. Segelstein, “The complex refractive index of water,” Ph.D. thesis, University of Missouri–Kansas City (1981).

14. G. P. Anderson, S. A. Clough, F. X. Kneizys, J. H. Chetwynd, and E. P. Shettle, AFGL Atmospheric Constituent Profiles (0-120km), AFGL-TR-86-0110 (OPI)Optical Physics Division, Air Force Geophysics Laboratory, Hanscom AFB (Air Force Geophysics Laboratory, 1986).

15. G. Thomas and K. Stamnes, Radiative Transfer in the Atmosphere and Ocean, Atmospheric and Space Science (Cambridge University Press, 1999), 1 [CrossRef] st ed.

16. M. I. Mishchenko, L. D. Travis, and A. a. Lacis, “Scattering, absorption, and emission of Light by small particles,” NASA Goddard Institute for Space Studies, New York Institute for Space Studies, New York pace Studies, New York New York (2002), pp. 1–486.

17. P. Yang, L. Bi, B. a. Baum, K.-N. Liou, G. W. Kattawar, M. I. Mishchenko, and B. Cole, “Spectrally consistent scattering, absorption, and polarization properties of atmospheric ice crystals at wavelengths from 0.2 to 100 µ m,” J Atmos. Sci. 70, 330–347 (2013 [CrossRef] ).

18. K. Stamnes, B. Hamre, J. Stamnes, G. Ryzhikov, M. Biryulina, R. Mahoney, B. Hauss, and A. Sei, “Modeling of radiation transport in coupled atmosphere-snow-ice-ocean systems,” J. Quant. Spectrosc. Radiat. Transf. 112, 714–726 (2011 [CrossRef] ).

19. Z. Lin, S. Stamnes, Z. Jin, I. Laszlo, S.-C. Tsay, W. Wiscombe, and K. Stamnes, “Improved discrete ordinate solutions in the presence of an anisotropically reflecting lower boundary: Upgrades of the DISORT computational tool,” J. Quant. Spectrosc. Radiat. Transf. 157, 119–134 (2015 [CrossRef] ).

20. K. Stamnes, S. C. Tsay, W. Wiscombe, and K. Jayaweera, “Numerically stable algorithm for discrete-ordinate-method radiative transfer in multiple scattering and emitting layered media,” Appl. Opt. 27, 2502–2509 (1988 [CrossRef] [PubMed] ).

21. T. Aoki, K. Kuchiki, M. Niwano, Y. Kodama, M. Hosaka, and T. Tanaka, “Physically based snow albedo model for calculating broadband albedos and the solar heating profile in snowpack for general circulation models,” J. Geophys. Res. Atmos. 116, 2156–2202 (2011 [CrossRef] ).

Sensor Name	SWIR Channels
Sensor Name	Band No.	CWL^a (µm)	FWHM^b (µm)	Res.^c
SGLI	SW01	1.052	0.02	1km
	SW02	1.385	0.02	1km
	SW03	1.625	0.20	250m
	SW04	2.208	0.05	1km

MODIS	5	1.242	0.020	500m
	26	1.382	0.030	1km
	6	1.629	0.024	500m
	7	2.126	0.050	500m

VIIRS	M8	1.238	0.020	750m
	M9	1.375	0.015	750m
	I3/M10	1.600	0.060	375/750m
	M11	2.257	0.050	750m

MODIS ch05	mean error (VZA 0 – 65°)	error (nadir)	error (VZA = 65°)
Method 1	−0.40	−0.32	−0.59
Method 2	−0.38	−0.29	−0.58
Method 3	−0.75	−0.80	−0.84
New method	0.23	0.23	0.22

MODIS ch06

Method 1	−0.421	−0.36	−0.53
Method 2	−0.37	−0.31	−0.49
Method 3	−0.24	−0.19	−0.39
New method	−0.14	−0.018	−0.21

MODIS ch07

Method 1	−8.5	−7.9	−9.4
Method 2	−14	−13	−15
Method 3	−11	−10	−12
New method	−1.1	−0.80	−1.7

SGLI SW01

Method 1	−0.049	−0.073	−0.033
Method 2	−0.026	−0.041	−0.018
Method 3	−0.15	−0.22	−0.10
New method	0.0051	0.0040	0.0066

SGLI SW03

Method 1	−8.6	−9.5	−7.6
Method 2	−13	−13	−11
Method 3	−1.2	−1.9	−0.64
New method	−0.94	−0.52	−1.7

SGLI SW04

Method 1	−3.2	−3.0	−3.9
Method 2	−3.5	−3.3	−4.1
Method 3	−1.8	−1.5	−2.6
New method	−1.6	−1.1	−2.4

Sensor Name	SWIR Channels
Sensor Name	Band No.	CWL^a (µm)	FWHM^b (µm)	Res.^c
SGLI	SW01	1.052	0.02	1km
	SW02	1.385	0.02	1km
	SW03	1.625	0.20	250m
	SW04	2.208	0.05	1km

MODIS	5	1.242	0.020	500m
	26	1.382	0.030	1km
	6	1.629	0.024	500m
	7	2.126	0.050	500m

VIIRS	M8	1.238	0.020	750m
	M9	1.375	0.015	750m
	I3/M10	1.600	0.060	375/750m
	M11	2.257	0.050	750m

MODIS ch05	mean error (VZA 0 – 65°)	error (nadir)	error (VZA = 65°)
Method 1	−0.40	−0.32	−0.59
Method 2	−0.38	−0.29	−0.58
Method 3	−0.75	−0.80	−0.84
New method	0.23	0.23	0.22

MODIS ch06

Method 1	−0.421	−0.36	−0.53
Method 2	−0.37	−0.31	−0.49
Method 3	−0.24	−0.19	−0.39
New method	−0.14	−0.018	−0.21

MODIS ch07

Method 1	−8.5	−7.9	−9.4
Method 2	−14	−13	−15
Method 3	−11	−10	−12
New method	−1.1	−0.80	−1.7

SGLI SW01

Method 1	−0.049	−0.073	−0.033
Method 2	−0.026	−0.041	−0.018
Method 3	−0.15	−0.22	−0.10
New method	0.0051	0.0040	0.0066

SGLI SW03

Method 1	−8.6	−9.5	−7.6
Method 2	−13	−13	−11
Method 3	−1.2	−1.9	−0.64
New method	−0.94	−0.52	−1.7

SGLI SW04

Method 1	−3.2	−3.0	−3.9
Method 2	−3.5	−3.3	−4.1
Method 3	−1.8	−1.5	−2.6
New method	−1.6	−1.1	−2.4

Fast yet accurate computation of radiances in shortwave infrared satellite remote sensing channels

Abstract

1. Introduction

1.1. Radiative transfer simulations of satellite sensor signals

1.1.1. Atmospheric IOPs

1.1.2. Surface IOPs

1.1.3. Computation of the sensor signal

2. Mean IOP methods

2.1. Chandrasekhar mean IOP

2.2. Refractive index mean IOP

2.3. Central wavelength method

3. Tests of the mean IOP methods

3.1. Black-sky overlying a snow surface

3.2. Atmosphere over a Lambertian surface

4. New sub-band IOP method

4.1. New sub-band IOP method applied to clear-sky cases

4.2. New sub-band IOP method applied to cloudy cases

5. Further improvement

6. Conclusion

Funding

Acknowledgments

References and links

Cited By

Figures (9)

Tables (2)

Equations (16)

Optics Express