In this section, we will discuss the MAPOL retrieval algorithm based on multi-angle polarimetric measurements and the associated radiative transfer forward model. The retrieval algorithm has been validated using both synthetic data (Gao et al., 2018) and RSP field measurements (Gao et al., 2019, 2020). To apply the retrieval algorithm to AirHARP measurements, we will first discuss the AirHARP instrument characteristics.

AirHARP measures the total and linearly polarized radiance at 60 viewing angles at the 660 nm band and at 20 viewing angles at the 440, 550, and 870 nm bands. Different from AirHARP, HARP2 reduces the number of viewing angles to 10 at 440, 550, and 870 nm and maintains 20 viewing angles at 660 nm in order to fulfill the bandwidth requirement and preserve information content as much as possible. HARP instruments (AirHARP, HARP CubeSat, and HARP2) use a modified three-way Phillips prism located after the front lens to split the incident light into the three orthogonal linear polarization states (0^∘, 45^∘, and 90^∘), which can be recombined to obtain the Stokes parameters L_t, Q_t, and U_t at the observational altitude (Puthukkudy et al., 2020). Circular polarization (Stokes parameter V) is not measured by any of the polarimeters in ACEPOL as it is negligible for atmospheric studies (Kawata, 1978). We use the total measured reflectance (ρ_t(λ)) and degree of linear polarization (DoLP, P_t(λ)) at the height of the aircraft with spectral dependencies hereafter implied, which are defined as

where F₀ is the extraterrestrial solar irradiance, μ₀ is the cosine of the solar zenith angle, and r is the Sun–Earth distance correction factor in astronomical units.

Based on the MAP measurements, the MAPOL retrieval algorithm is developed to derive both the aerosol properties and the water-leaving signal simultaneously. The retrieval algorithm minimizes the difference between the MAP measurements and the forward model simulations computed from vector radiative transfer simulations (Zhai et al., 2009, 2010). By assuming the measurement and modeling uncertainties follow Gaussian statistical distributions, the retrieval parameters can be estimated through Bayesian theory using the cost function χ² to quantify the difference between the measurement and the forward model simulation (Rogers, 2000):

where ρ_t and P_t are the measured reflectance and DoLP as defined in Eqs. (1) and (2), and ${\mathit{\rho }}_{\mathrm{t}}^{\mathrm{f}}$ and $P_{\mathrm{t}}^{\mathrm{f}}$ are the corresponding quantities computed from the forward model. The state vector x contains all retrieval parameters, such as the aerosol size and refractive indices; the subscript i stands for the index of the measurements at different viewing angles and wavelengths; and N is the total number of the measurements used in the retrieval. For AirHARP measurements, the maximum value of N is 240, twice of the total number of viewing angles. The total uncertainties of the reflectance and DoLP used in the algorithm are denoted as σ_ρ and σ_P, which are contributed by both the measurement uncertainties σ_m and the forward model uncertainties σ_f (more details in Sect. 3.3):

One important component of σ_m is the calibration uncertainty. AirHARP was calibrated in the lab with an accuracy of 3 %–5 % for reflectance and 0.005 for DoLP (McBride et al., 2019). In-flight uncertainty for the AirHARP DoLP is conservatively estimated to be at most 0.01 without an onboard calibrator. In this study, we adopted the calibration uncertainty for reflectance as ${\mathit{\sigma }}_{\mathit{\rho },\text{cal}}=\mathrm{3}\mathrm{\%}{\mathit{\rho }}_{\mathrm{t}}$ and for DoLP as ${\mathit{\sigma }}_{P,\text{cal}}=\mathrm{0.01}$ for all four bands. The accuracy of the HARP2 measurements can be further improved through onboard calibration (McBride et al., 2020; Puthukkudy et al., 2020). In this study, we considered the total measurement uncertainties as the contributions only from the calibration (σ_cal):

However, other contributions such as spatial variability of the geophysical properties may also contribute to the measurement uncertainties, which will be discussed in Sect. 5. Furthermore, noise correlation is an import influence on the retrieval accuracy (Knobelspiesse et al., 2012) that is ignored in this study due to the lack of knowledge on this characteristic for AirHARP.

As observed by AirHARP (Puthukkudy et al., 2020) and RSP measurements (Gao et al., 2020), the sunglint angular pattern cannot be well modeled by an isotropic Cox–Munk model. Using these data will require characterization of the corresponding measurement and model uncertainties. To minimize the impact of sunglint in our discussions, we removed the signals within an angle range of 0^∘ to 40^∘ relative to the solar specular reflection direction.

The forward model uncertainties σ_f include the uncertainties of the radiative transfer calculation and uncertainties due to the incompleteness of the model to describe the system. However, the latter are difficult to quantify; we will discuss the possible sources for them in the next section. For convenience, we will only consider the uncertainties of the NN forward model (σ_NN) and the radiative transfer simulation used for generating the NN training data (σ_RT) as

σ_NN is evaluated by comparing with synthetic multi-angle AirHARP measurements discussed in Sect. 3.3.

To fully utilize the information contained in the AirHARP measurements, the forward model needs to achieve an accuracy level much better than the measurement uncertainties. This becomes the goal of the NN training in the next section. Detailed comparisons of the forward model uncertainties and the measurement uncertainties will be provided in the next section. To minimize the cost function defined in Eq. (3), we use an optimization method, called the subspace trust-region interior reflective (STIR) approach (Branch et al., 1999) as implemented in the Python SciPy package (Virtanen et al., 2020), to solve the state parameters x iteratively. The method is based upon the Levenberg–Marquardt method (Moré, 1978) and shows good stability for the boundary constraints.

2.1 Forward model

We used a vector radiative transfer model based on the successive order of scattering method for coupled atmosphere and ocean systems (Zhai et al., 2009, 2010) to model the measured reflectance and DoLP. The atmosphere is configured as three layers: a top molecular layer above the aircraft, a molecular layer below the aircraft in the middle, and an aerosol and molecular mixing layer on the bottom with a height of 2 km. Aerosols are assumed to be uniformly distributed in the mixing layer as shown in the left panel of Fig. 1. The same vertical structure of the atmosphere was successfully used in the inversion of RSP data (Gao et al., 2019, 2020).

Figure 1Panel (a) shows the coupled atmosphere and ocean system used in FastMAPOL including the atmosphere, ocean surface, and ocean body. Panel (b) represents a system used for atmospheric correction which only has atmosphere and ocean surface without scattering in the ocean body. The atmospheres in both systems are modeled as the same three layers. TOA indicates the top of the atmosphere. The bottom of the atmosphere (BOA) and the top of the ocean (TOO) indicate the locations just above and below the ocean surface, respectively. All quantities shown in the figures need to be computed from the forward model and represented by the NN for efficient calculations. Symbols are defined in Table 1.

Download

The atmospheric surface pressure is assumed to be 1 atm (standard atmosphere pressure), which is consistent with the value discussed in Sect. 5. Anisotropic molecular Rayleigh scatterings are accounted for in Hansen and Travis (1974). The molecular absorption properties are computed by the hyperspectral line-by-line atmospheric radiative transfer simulator (ARTS) (Buehler et al., 2005) with the molecular absorption parameters of oxygen, water vapor, methane, and carbon dioxide from the HITRAN database (Gordon et al., 2017). The gas absorption of ozone and nitrogen dioxide are from Gorshelev et al. (2014), Serdyuchenko et al. (2014), and Bogumil et al. (2003), respectively. The hyperspectral absorption coefficients are then averaged within the instrument spectral response function and used in the multiple scattering radiative transfer simulation (Zhai et al., 2009, 2010, 2018). The molecular profile used is the US standard atmospheric constituent profiles (Anderson et al., 1986). Ozone is the most important gas that influences the absorption transmittance at the AirHARP bands of 550 and 660 nm. For the application to AirHARP measurements in ACEPOL, we use the ozone column density as a free parameter with values from the Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA-2) developed by NASA's Global Modeling and Assimilation Office (Gelaro et al., 2017) to rescale the molecular absorption optical depth calculated under the abovementioned standard atmospheric profile.

Aerosols are diverse in size, composition, and morphology. To capture their variation in the atmosphere, we modeled the size and refractive index for both fine and coarse modes. The aerosol size is represented by the volume density distribution as a combination of five lognormal distributions:

$$\begin{array}{}\text{(8)}& {\displaystyle }{\displaystyle \frac{\mathrm{d}V\left(r\right)}{\mathrm{d}\text{ln}r}}=\sum _{i=\mathrm{1}}^{\mathrm{5}}{\displaystyle \frac{V_i}{\sqrt{\mathrm{2}\mathit{\pi }}{\mathit{\sigma }}_{\mathrm{v},i}}}\exp \left[-{\displaystyle \frac{(\ln r-\ln r_{\mathrm{v},i}{)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{\mathrm{v},i}^{\mathrm{2}}}}\right],\end{array}$$

where V_i is the column volume density for each submode; the mean radius r_i is fixed with values of 0.1, 0.1732, 0.3, 1.0, and 2.9 µm; and the standard deviation σ_i is fixed with values of 0.35, 0.35, 0.35, 0.5, and 0.5 (Dubovik et al., 2006; Xu et al., 2016). The first three submodes are categorized as the fine-mode aerosol, while the last two submodes are the coarse mode. All aerosols are assumed to be spherical in the current forward model. The nonspherical particle shape is important in the aerosol model (Dubovik et al., 2006) and will be considered in future studies. The aerosol refractive index spectra for the fine and coarse modes are represented by the principal component analysis in MAPOL (Wu et al., 2015; Gao et al., 2018) as

$$\begin{array}{}\text{(9)}& {\displaystyle }{\displaystyle }m\left(\mathit{\lambda }\right)=m_{\mathrm{0}}+{\mathit{\alpha }}_{\mathrm{1}}{p}_{\mathrm{1}}\left(\mathit{\lambda }\right),\end{array}$$

where p₁(λ) is the first-order principal component computed from the aerosol refractive index dataset including water, sea salt, dust-like particles, biomass burning, soot, sulfate, water-soluble aerosols, and industrial aerosols (Shettle and Fenn, 1979; d'Almeida et al., 1991). m₀ and α₁ are two coefficients to determine the spectrum. For the application to AirHARP bands, p₁(λ) for the real part of the refractive index is approximately spectrally flat for both the fine- and coarse-mode aerosols within the AirHARP spectral range. We further assume the spectral shape for the imaginary refractive spectra is also flat. Two parameters can be combined into one to represent the refractive index. Hereafter, we only refer one independent parameter for each refractive index spectrum. Therefore, only four independent parameters are required to determine the real and imaginary refractive index spectra for the fine and coarse modes. With the aerosol size and refractive index, the polarimetric single scattering properties are modeled by the Lorenz–Mie theory and computed by the code developed by Mishchenko et al. (2002).

For the ocean layer shown in Fig. 1, two ocean bio-optical models are implemented in the forward model of MAPOL: one with chlorophyll a concentration (Chl a; mg m⁻³) as the single parameter applicable to open-ocean optical properties and the other with seven parameters more suitable to fully describe complex coastal waters (Gao et al., 2019). Since the waters are generally clear within the ocean scenes in this study (Gao et al., 2020), an open-ocean model is used for both NN training and retrievals. The optical properties of open-ocean waters include contributions from pure seawater, colored dissolved organic matter (CDOM), and phytoplankton, where the CDOM and phytoplankton absorption coefficients as well as the phytoplankton scattering coefficient and phase function are parameterized by Chl a (Gao et al., 2019). A complex costal water model for NN trainings will be investigated in separate studies. The ocean surface roughness is modeled by the isotropic Cox–Munk model with a scalar wind speed. Whitecap is not considered in the current study.

In summary, the parameters used to represent the forward model include five volume densities (one for each submode); four independent parameters for the refractive indices of fine and coarse modes; and one parameter for wind speed, ozone column density, and Chl a, respectively. Three additional geometric parameters are used to set up the system, including the solar zenith angle, viewing zenith angle, and relative viewing azimuth angle. Therefore, it requires a total of 15 parameters to conduct the radiative transfer calculation, with a total of 11 independent state parameters that can be retrieved from optimizing the cost function as defined in Eq. (3).

2.2 Remote sensing reflectance

An important task for the joint retrievals is to obtain the water-leaving signal, which is often represented in ocean color studies by the spectral remote sensing reflectance defined as $R_{\text{rs}}=L_{\mathrm{w}}^{+}/E_{\mathrm{d}}^{+}$, where $E_{\mathrm{d}}^{+}$ is the downwelling irradiance and $L_{\mathrm{w}}^{+}$ is the water-leaving radiance just above the ocean surface (Mobley et al., 2016). The remote sensing reflectance can be derived from the water-leaving reflectance reaching the sensor (ρ_w) via

$$\begin{array}{}\text{(10)}& {\displaystyle }{\displaystyle }{R}_{\text{rs}}=\left[{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{w}}({\mathit{\theta }}_{\mathrm{0}},{\mathit{\theta }}_{\mathrm{v}},{\mathit{\varphi }}_{\mathrm{v}})}{\mathit{\pi }{r}^{\mathrm{2}}}}\right]\times \left[{\displaystyle \frac{C_{\text{BRDF}}({\mathit{\theta }}_{\mathrm{0}},{\mathit{\theta }}_{\mathrm{v}},{\mathit{\varphi }}_{\mathrm{v}})}{T_{\mathrm{d}}^{\mathrm{f},+}\left({\mathit{\theta }}_{\mathrm{0}}\right)t_{\mathrm{u}}^{\mathrm{f},+}({\mathit{\theta }}_{\mathrm{0}},{\mathit{\theta }}_{\mathrm{v}},{\mathit{\varphi }}_{\mathrm{v}})}}\right],\end{array}$$

where θ₀ and θ_v are the solar and viewing zenith angles, and ϕ_v is the relative viewing azimuth angle. ρ_w represents the signals originating from scattering in the ocean that reached the sensor, which can be derived from the atmospheric correction process as

$$\begin{array}{}\text{(11)}& {\displaystyle }{\displaystyle }{\mathit{\rho }}_{\mathrm{w}}({\mathit{\theta }}_{\mathrm{0}},{\mathit{\theta }}_{\mathrm{v}},{\mathit{\varphi }}_{\mathrm{v}})={\mathit{\rho }}_{\mathrm{t}}({\mathit{\theta }}_{\mathrm{0}},{\mathit{\theta }}_{\mathrm{v}},{\mathit{\varphi }}_{\mathrm{v}})-{\mathit{\rho }}_{\mathrm{t},\text{atm}+\text{sfc}}^{\mathrm{f}}({\mathit{\theta }}_{\mathrm{0}},{\mathit{\theta }}_{\mathrm{v}},{\mathit{\varphi }}_{\mathrm{v}}),\end{array}$$

where ρ_t is the measured total reflectance as defined in Eq. (1), and ${\mathit{\rho }}_{\mathrm{t},\text{atm}+\text{sfc}}^{\mathrm{f}}$ is the reflectance from a system with only atmosphere and ocean surface (Mobley et al., 2016) as represented in the right panel of Fig. 1. The same formalism has been used to derive R_rs from RSP measurements Gao et al. (2019, 2020).

Table 1Definition of the symbols for the quantities computed from the forward model (indicated by the superscript f) as shown in Fig. 1.

Download Print Version | Download XLSX

The downwelling irradiance transmittance $T_{\mathrm{d}}^{\mathrm{f}}$ is for the solar irradiance from TOA to the surface, and the upwelling radiance transmittance $t_{\mathrm{u}}^{\mathrm{f},+}$ is for the water-leaving radiance from BOA to the sensor (Gao et al., 2019). Both $T_{\mathrm{d}}^{\mathrm{f}}$ and $t_{\mathrm{u}}^{\mathrm{f},+}$ are denoted in Fig. 1 and represented as follows:

$$\begin{array}{}\text{(12)}& {\displaystyle }{T}_{\mathrm{d}}^{\mathrm{f},+}\left({\mathit{\theta }}_{\mathrm{0}}\right)={\displaystyle \frac{E_{\mathrm{d}}^{\mathrm{f},+}\left({\mathit{\theta }}_{\mathrm{0}}\right)}{{\mathit{\mu }}_{\mathrm{0}}{F}_{\mathrm{0}}}},\text{(13)}& {\displaystyle }\begin{array}{rl}& t_{\mathrm{u}}^{\mathrm{f},+}({\mathit{\theta }}_{\mathrm{0}},{\mathit{\theta }}_{\mathrm{v}},{\mathit{\varphi }}_{\mathrm{v}})=\\ & \left({\displaystyle \frac{{\mathit{\rho }}_{\mathrm{t}}^{\mathrm{f}}({\mathit{\theta }}_{\mathrm{0}},{\mathit{\theta }}_{\mathrm{v}},{\mathit{\varphi }}_{\mathrm{v}})-{\mathit{\rho }}_{\mathrm{t},\text{atm}+\text{sfc}}^{\mathrm{f}}({\mathit{\theta }}_{\mathrm{0}},{\mathit{\theta }}_{\mathrm{v}},{\mathit{\varphi }}_{\mathrm{v}})}{{\mathit{\rho }}_{\mathrm{t}}^{\mathrm{f},+}({\mathit{\theta }}_{\mathrm{0}},{\mathit{\theta }}_{\mathrm{v}},{\mathit{\varphi }}_{\mathrm{v}})-{\mathit{\rho }}_{\mathrm{t},\text{atm}+\text{sfc}}^{\mathrm{f},+}({\mathit{\theta }}_{\mathrm{0}},{\mathit{\theta }}_{\mathrm{v}},{\mathit{\varphi }}_{\mathrm{v}})}}\right),\end{array}\end{array}$$

where ${\mathit{\rho }}_{\mathrm{t}}^{\mathrm{f},+}$ and ${\mathit{\rho }}_{\mathrm{t},\text{atm}+\text{sfc}}^{\mathrm{f},+}$ are reflectance just above ocean surface also denoted in Fig. 1 and Table 1.

To remove the dependency of R_rs on the solar and viewing geometries, a bidirectional reflectance distribution function (BRDF) correction C_BRDF is applied to adjust R_rs to the observation with a zenith sun and a nadir viewing direction as defined by Morel et al. (2002):

$$\begin{array}{}\text{(14)}& \begin{array}{rl}{C}_{\text{BRDF}}({\mathit{\theta }}_{\mathrm{0}},{\mathit{\theta }}_{\mathrm{v}},{\mathit{\varphi }}_{\mathrm{v}})& ={\displaystyle \frac{{\mathfrak{R}}_{\mathrm{o}}\left(W\right)}{\mathfrak{R}({\mathit{\theta }}_{\mathrm{v}}^{\prime },{\mathit{\varphi }}_{\mathrm{v}},W)}}\\ & \times {\displaystyle \frac{{\mathit{\rho }}_{\mathrm{t}}^{\mathrm{f},-}(\mathrm{0},\mathrm{0})}{T_{\mathrm{d}}^{\mathrm{f},-}\left(\mathrm{0}\right)}}{\left[{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{t}}^{\mathrm{f},-}({\mathit{\theta }}_{\mathrm{0}},{\mathit{\theta }}_{\mathrm{v}}^{\prime },{\mathit{\varphi }}_{\mathrm{v}})}{T_{\mathrm{d}}^{\mathrm{f},-}\left({\mathit{\theta }}_{\mathrm{0}}\right)}}\right]}^{-\mathrm{1}}.\end{array}\end{array}$$

${\mathfrak{R}}_{\mathrm{0}}/\mathfrak{R}$ accounts for reflection and refraction effects when light propagates through the ocean interface. (${\mathit{\theta }}_{\mathrm{v}}^{\prime },{\mathit{\varphi }}_{\mathrm{v}}$) is the direction of the upwelling radiance beneath the sea surface, where ${\mathit{\theta }}_{\mathrm{v}}^{\prime }$ is defined through Snell's law:

$$\begin{array}{}\text{(15)}& {\displaystyle }{\displaystyle }\sin {\mathit{\theta }}_{\mathrm{v}}^{\prime }=\sin {\mathit{\theta }}_{\mathrm{v}}/n_{\mathrm{w}},\end{array}$$

with water refractive index n_w. C_BRDF in its original form is defined using the radiance and irradiance just below the ocean surface (Morel et al., 2002); here we have converted all quantities into the radiance reflectance ${\mathit{\rho }}_{\mathrm{t}}^{\mathrm{f},-}$ and the irradiance transmittance $T_{\mathrm{d}}^{\mathrm{f},-}\left({\mathit{\theta }}_{\mathrm{0}}\right)$ similar to Eqs. (1) and (12).

To compute the remote sensing reflectance from the multi-angle AirHARP measurement, we only consider the reflectance at the minimum viewing zenith angle for each wavelength and apply the atmospheric correction and BRDF correction as discussed above. For ${\mathit{\theta }}_{\mathrm{v}}^{\prime }<\mathrm{15}{}^{\circ }$ (${\mathit{\theta }}_{\mathrm{v}}<\mathrm{20}{}^{\circ }$), the ${\mathfrak{R}}_{\mathrm{o}}/\mathfrak{R}$ factor is approximately a constant value of 1, but for larger θ_v angles, the ratio increases with both wind speed and θ_v (Morel and Gentili, 1996; Morel et al., 2002). In this study we ignored the ${\mathfrak{R}}_{\mathrm{o}}/\mathfrak{R}$ factor in Eq. (14), which will not impact R_rs calculation from synthetic data due to the small viewing zenith angle used but may cause underestimation of R_rs at the edge of the image, as will be discussed in Sects. 4 and 5. All quantities denoted in Fig. 1 and Table 1 need to be determined for the forward model and the calculation of remote sensing reflectance and will be represented by NN models.

Deep NN models are developing rapidly due to the advancement in machine-learning infrastructure and demands in broad applications (Goodfellow et al., 2016) and are demonstrated to be efficient in approximating physical functions (Lin et al., 2017). In this study, we employed the deep feed-forward NN (Goodfellow et al., 2016) to represent the MAP measurements. In this section, we will discuss the procedures to train the NN forward models for reflectance and DoLP respectively, with their performance evaluated.

3.1 Training data

To train a NN that can represent the forward model accurately for the AirHARP measurements from the ACEPOL field campaign, we generated the training data according to the average aircraft height of 20.1 km on the day of 23 October 2017 from ACEPOL. We simulated 21 000 cases according to the forward model as discussed in the previous section by considering general aerosol and ocean properties, as well as a large range of solar and viewing geometries with the minimum and maximum values of all parameters summarized in Table 2. The ranges of solar zenith angle θ₀, viewing zenith angle θ_v, and relative viewing azimuth angle ϕ_v are from 0^∘ to 70^∘, 60^∘ and 180^∘, respectively. The reflectance and DoLP with a viewing azimuth angle larger than 180^∘ can be evaluated by the corresponding value less than 180^∘ due to symmetry with respect to the principal plane (defined by ${\mathit{\varphi }}_{\mathrm{v}}=\mathrm{0}{}^{\circ }$ and ${\mathit{\varphi }}_{\mathrm{v}}=\mathrm{180}{}^{\circ }$). For each solar zenith angle, the polarized reflectance is calculated for all viewing angles within the aforementioned ranges with an angular resolution of 1^∘. The solar zenith angle, ozone column density, refractive index, and wind speed are randomly sampled from a uniform distribution. Chl a is randomly sampled from a log-uniform distribution. The fine-mode volume fraction is sampled uniformly within [0, 1], which is then randomly partitioned to each submode. To maintain a uniform distribution of the total AOD, we sampled the AOD at 550 nm within [0, 0.5] in a linear scale. The volume density V_i of each submode is determined by the total aerosol optical depth and volume fraction for each mode. Figure 2 shows one example simulation dataset for the angular distribution of reflectance and DoLP.

Table 2Parameters used to represent the atmosphere and ocean system for the radiative transfer simulation and NN training. θ₀ and θ_v are the solar and viewing zenith angles. ϕ_v is the relative viewing azimuth angle. V_i denotes the five volume densities defined in Eq. (8). m_r and m_i are the real and imaginary parts of the refractive index. Ozone column density ($n_{{\mathrm{O}}_{\mathrm{3}}}$) in the atmosphere, ocean surface wind speed, and Chl a are also provided. The minimum (min) and maximum (max) values determine the parameter ranges used to generate NN training data, which are also the constraints in the retrieval algorithm. The initial values are the ones used in the retrieval optimization algorithm, where θ₀, θ_v, ϕ_v, and $n_{{\mathrm{O}}_{\mathrm{3}}}$ are assumed to be known from inputs.

Download Print Version | Download XLSX

Figure 2The reflectance (a) and DoLP (b) from radiative transfer simulation with the wind speed of 4.13 m s⁻¹, the aerosol optical depth of 0.26, Chl a of 0.05 mg m⁻³, and ozone column density of 196 DU. The antisolar point is indicated by the red asterisk with a solar zenith angle ${\mathit{\theta }}_{\mathrm{0}}=\mathrm{46.41}{}^{\circ }$. θ_v and ϕ_v indicate the viewing zenith and relative azimuth angles. The principal plane is defined by the viewing azimuth angle of 0^∘ and 180^∘.

Download

We randomly selected 20 000 cases out of the total 21 000 simulated cases for the training and validation processes, and the remaining 1000 random cases will be used as test cases to evaluate the NN accuracy, which will be discussed in the next section. To enable the NNs to predict reflectance and DoLP at any given viewing geometry, for each case, we sampled 100 random pairs of viewing zenith and azimuth angles. If the sampled angles fall outside of the predefined angular grids, values from spline interpolation are used. The sunglint angles within an angle of 40^∘ to the solar specular reflection direction are removed. Approximately 1 million data points are obtained for each wavelength for training.

To maintain both flexibility and efficiency, we trained two NN models for reflectance and DoLP respectively in the next section. Reflectance and DoLP have different accuracy requirements as discussed in Sect. 2 and also differ in angular variations as shown in the Fig. 2; therefore, it is convenient to control their accuracy through separated training procedures.

3.2 Neural network training

A feed-forward NN can be defined recursively with one input layer, one output layer, and k hidden layers (Aggarwal, 2018):

$$\begin{array}{}\text{(16)}& {\displaystyle }& {\displaystyle }{\mathit{h}}_{\mathrm{1}}=\mathrm{\Phi }({\mathbf{W}}_{\mathrm{1}}^{\mathrm{T}}\mathit{x}+{\mathit{b}}_{\mathrm{1}}),\text{(17)}& {\displaystyle }& {\displaystyle }{\mathit{h}}_{p+\mathrm{1}}=\mathrm{\Phi }({\mathbf{W}}_{p+\mathrm{1}}^{\mathrm{T}}{\mathit{h}}_p+{\mathit{b}}_{p+\mathrm{1}}),p=\mathrm{1},\mathrm{\dots },k-\mathrm{1},\text{(18)}& {\displaystyle }& {\displaystyle }\mathit{y}={\mathbf{W}}_{k+\mathrm{1}}^{\mathrm{T}}{\mathit{h}}_k+{\mathit{b}}_{k+\mathrm{1}},\end{array}$$

where x is the input parameter vector including all 15 parameters needed to define the forward model as listed in Table 2. Here x not only contains the retrieval parameters in the state vector defined in Eq. (3) but also includes additional non-retrieval parameters such as the solar zenith angle, viewing zenith and azimuth angles, and the ozone column density. y is a four-dimensional output vector for reflectance or DoLP at the four AirHARP bands. The weight matrix W_p+1 connects the pth and (p+1)th NN layers. The bias vector for the (p+1) layer is defined as b_p+1. The output of each layer h_p+1 becomes the input of the next layer as shown in Eq. (17). k is the number of hidden layers, and k+1 refers to the output layer. In this study, we tested several NN architectures and eventually chose three hidden layers with the number of nodes of 1024, 256, and 128 as shown in Table 3. The nonlinear activation function Φ used in this model is the Leaky ReLU function, which is defined as

$$\begin{array}{}\text{(19)}& {\displaystyle }{\displaystyle }\mathrm{\Phi }\left(z\right)=max(\mathrm{0},z)+\mathrm{0.01}\times min(\mathrm{0},z).\end{array}$$

Both Leaky ReLU and ReLU (defined as $max(\mathrm{0},z)$) activation functions are simple in their mathematical forms and are tested in our NN trainings. Leaky ReLU is eventually chosen due to its slightly better accuracy achieved than the NN with ReLU.

Table 3The accuracy of the NN for the corresponding quantities in terms of the RMSE (σ_NN) of the difference between the NN-predicted values and the truth values from radiative transfer simulation. The NN architecture denotes the number of the nodes in each layer. The corresponding NN data sizes are indicated. Remote sensing reflectance is computed by Eq. (10) using the NNs for ${\mathit{\rho }}_{\mathrm{t},\text{atm}+\text{sfc}}^{\mathrm{f}}$ and $[C_{\text{BRDF}}/T_{\mathrm{d}}{t}_{\mathrm{u}}]$ as discussed in Sect. 3.4. (The percentage values listed below in the parenthesis are the percentage uncertainties defined as the RMSE of the percentage difference between the RT simulation and NN predictions.). n/a means not applicable.

Download Print Version | Download XLSX

The training process is to minimize the cost function defined as the mean square error between the training data generated from radiative transfer simulations and the NN-predicted values (Aggarwal, 2018). All parameters in the neural network weight matrices and bias vectors, over 670 000 numbers, need to be trained. With this large number of parameters, it is a challenging task to avoid overfitting where the model works well for the training dataset but poorly for the dataset not used in the training process. Several training procedures are performed for reflectance and DoLP data to avoid overfitting and improve NN performance.

1. Both input and output data are normalized before training. We normalize the input data into the range of [0, 1] using the minimum and maximum values from the datasets as listed in Table 2. The reflectance and DoLP in the output layers are normalized by dividing their standard deviation of the training data at each wavelength.

2. The Adam (short for Adaptive Moment Estimation) optimization algorithm (Kingma and Ba, 2015) with weight decay regularization (Loshchilov and Hutter, 2019) is used to update the weights and bias of the NN. The training dataset is divided into multiple mini-batches, each with 1024 random samples. The training iterations loop through all mini-batches in the training data (each loop is called an epoch). Convergence requires training through multiple epochs, where mini-batches are resampled in each epoch.

3. The learning rate determines the step size in the parameter update. We use an exponential decay schedule to reduce the learning rate: we start with a learning rate of 0.005 and reduce the learning rate by a factor of 10 every 200 epochs.

4. To monitor overfitting in the training process, we split the data into 70 % for training and 30 % for validation. We conduct the optimization based on the training dataset, and in the meantime we monitor the performance of training by applying the NN model to the validation dataset. To avoid overfitting, the early-stopping approach is employed where the training is stopped when the cost function on the validation dataset stops to reduce for a threshold of 50 epochs.

The machine-learning Python library PyTorch is used for the training (Paszke et al., 2019). The trained NN model is used to replace the radiative transfer model to compute the reflectance and DoLP in the retrieval algorithm. The Jacobian matrix used in the optimization is computed by the finite difference approximation of the partial derivatives of reflectance and DoLP with respect to the retrieval parameters. Here central difference method is used. Note that the Jacobian matrix can also be computed analytically from the NN model using the automatic differentiation techniques based on the chain rule of differentiation (Baydin et al., 2018). This will be a topic in our future studies.

3.3 Neural network accuracy

After training the NN model, we evaluated its accuracy using synthetic AirHARP measurements generated from the 1000 simulation cases which have not been used in the training and validation process. Each simulation dataset includes polarized reflectance on regular viewing angle grids, which are interpolated to the viewing geometry of AirHARP to create synthetic measurement data and compare with the NN predictions. Glint angles are excluded from the comparison because the NNs are not trained over these angles. As the example shown in Fig. 3, both the reflectance and DoLP are in good agreement between the synthetic data and the NN results, where the maximum absolute differences for reflectance and DoLP are within 0.001 and 0.0025. This translates to a difference for both reflectance and DoLP mostly less than 1 % for bands 440, 550, and 670 nm. The maximum percentage difference can be as large as 3 % for 870 nm bands due to the small reflectance magnitude.

Figure 3The synthetic HARP reflectance (a) and DoLP (b) sampled from the radiative transfer data shown in Fig. 2. Panels (c) and (d) indicate the difference between the NN predictions and RT simulations. Panels (e) and (f) indicate their percentage differences. The positive and negative signs of the viewing zenith angles indicate the azimuth angles of ${\mathit{\varphi }}_{\mathrm{v}}=\mathrm{116.2}{}^{\circ }$ and $\mathrm{180}{}^{\circ }+{\mathit{\varphi }}_{\mathrm{v}}$.

Download

Figure 4Comparison between the radiative transfer simulation and NN prediction: (a, c, e) reflectance (ρ) and (b, d, f) DoLP (P). The scatter plots are shown in panels (a) and (b), the difference in panels (c) and (d), and the percentage difference in panels (e) and (f). For each plot, the data points for the 550, 660, and 870 nm bands are shifted upward by constant offsets consecutively as indicated by the solid cyan lines.

Download

The comparison with all 1000 synthetic datasets and their NN predictions are shown in Fig. 4. The mean absolute error (MAE) and the root mean square error (RMSE) between the simulation data (T_i) and the NN-predicted data (R_i) shown in Fig. 4 are defined as

$$\begin{array}{}\text{(20)}& {\displaystyle }& {\displaystyle }\text{MAE}={\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N|R_i-T_i|,\text{(21)}& {\displaystyle }& {\displaystyle }\text{RMSE}=\sqrt{{\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N(R_i-T_i{)}^{\mathrm{2}}}.\end{array}$$

Both MAE and RMSE are useful metrics, where MAE is less dependent on outliers compared to RMSE.

Analysis shows that the statistics of the differences between the NN prediction and the RT simulations as shown in Fig. 4 can be well modeled by Gaussian distributions and characterized by RMSE. Therefore the RMSE is used to represent the NN uncertainties for both reflectance (σ_ρ,NN) and DoLP (σ_ρ,NN) and will be incorporated into the total uncertainties in the cost function. Table 3 summarizes the uncertainties of the NN models. The σ_ρ,NN at 440 nm is 0.0006, which decreases to 0.0004 at 870 nm. However, due to the smaller reflectance magnitude at 870 nm, the corresponding RMSE for the percentage reflectance difference as shown in Fig. 4 is increased from 0.4 % at 440 nm to 1.0 % at 870 nm. For DoLP, the maximum σ_P,NN is 0.003 at 870 nm, which decreases to 0.0016 at 440 nm. The uncertainties can be further improved with more training data points. For the readers' information, RMSE of the NN model trained with 20 000 cases (1 million data points) decreases by a factor of $\sqrt{\mathrm{2}}$ in comparison with the one using 10 000 cases (0.5 million data points). It takes 0.01 s in a single-core CPU (AMD EPYC processor) or 1 ms in a GPU (GeForce GTX 1060) to predict all 120 angles for both reflectance and DoLP in the NN forward model. Furthermore, the data sizes for the NNs are minimal, which are 1.2 MB for the reflectance and DoLP and less than 100 KB for the factor $[C_{\text{BRDF}}/T_{\mathrm{d}}{t}_{\mathrm{u}}]$ (details in Sect. 3.4 as shown in Table 3).

The assessment of the NN accuracy is relative to the synthetic measurements simulated by the vector radiative transfer simulations. To account for the modeling uncertainties of the forward model σ_f, we consider both the NN accuracy σ_NN and the numerical accuracy of the radiative transfer simulations σ_RT for reflectance and DoLP, respectively. Uncertainties due to incomplete assumptions in the forward model are not considered. Several internal parameters determine the numerical accuracy of the radiative transfer simulations. In the framework of the successive order of scattering (Zhai et al., 2008, 2009), these parameters include the number of scattering orders (N_s), the number of Gaussian quadratures for discretizing the viewing zenith angle in the atmosphere (P_a) and ocean (P_o), the order of Fourier decomposition (M) for the viewing azimuth angle, and the order of Legendre expansion (L) of the single scattering phase function. In this study, we chose N_s=20, P_a=32, P_o=64, M=32, and L=32, which has a higher accuracy than the radiative transfer forward model directly used in our previous retrieval studies (Gao et al., 2020).

To quantify the accuracy of the radiative transfer calculation used for generating training data (σ_RT), we simulated an additional 1000 synthetic AirHARP datasets with all internal parameters doubled as the most rigorous calculations, and the viewing angular resolution was reduced from 1^∘ to 0.5^∘ in order to reduce interpolation errors. The resultant reflectance and DoLP values are compared between these two sets of radiative transfer calculations. The RMSE for each band can be used as a measure of the accuracy for the radiative transfer calculation used to generate the training data (σ_RT). The uncertainties σ_RT for reflectance and DoLP are summarized in Table 4, with reflectance uncertainties less than 0.00015 and DoLP uncertainties less than 0.0007 for all AirHARP bands. σ_ρ,RT is about 4 times smaller than the NN uncertainties, and σ_P,RT is about 4 to 10 times smaller. The measurement uncertainties from calibration (σ_cal) are also summarized in Table 4. The total forward model uncertainties ${\mathit{\sigma }}_{\mathrm{f}}=\sqrt{{\mathit{\sigma }}_{\text{RT}}^{\mathrm{2}}+{\mathit{\sigma }}_{\text{NN}}^{\mathrm{2}}}$ as in Eq. (7) are much smaller than σ_cal. The overall uncertainties used in the retrieval cost function in Eq. (3) are dominated by the measurement contributions.

Table 4Comparisons of the uncertainties for reflectance (ρ) and DoLP (P) for both measurement and forward model including calibration uncertainty (σ_cal), the radiative transfer simulation uncertainty (σ_RT), and the NN uncertainty (σ_NN). The same σ_ρ,NN and σ_P,NN have been shown in Table 3 and are repeated here for comparisons. (As in Table 3, the percentage values listed in the table indicate the percentage uncertainties.)

Download Print Version | Download XLSX

Furthermore, in this study higher accuracies from the radiative transfer simulations are used for the NN training for comparison with the accuracies from the radiative transfer model directly used in our previous retrieval algorithm. Since the simulations of the training data can be conducted independent of the retrieval algorithm, higher computational costs can be accommodated to improve NN forward model accuracy. After the NN model is trained, the model can be applied to the retrieval algorithm through efficient matrix operations.

3.4 Neural network model for remote sensing reflectance

As discussed in Sect. 2.2, the water-leaving signals are represented by the remote sensing reflectance as defined in Eq. (10) (Mobley et al., 2016). To conduct the atmospheric correction in Eq. (11), we need to determine the reflectance ${\mathit{\rho }}_{\mathrm{t},\text{atm}+\text{sfc}}^{\mathrm{f}}$ at the aircraft level, transmittance $t_{\mathrm{u}}^{\mathrm{f}}$ and $T_{\mathrm{d}}^{\mathrm{f}}$, and the BRDF correction coefficient C_BRDF. Based on Eq. (10), we combined $T_{\mathrm{d}}^{\mathrm{f},+}$, $t_{\mathrm{u}}^{\mathrm{f},+}$, and C_BRDF into a single term denoted as $[C_{\text{BRDF}}/T_{\mathrm{d}}{t}_{\mathrm{u}}]$. To efficiently determine R_rs, two NNs need to be trained to represent ${\mathit{\rho }}_{\mathrm{t},\text{atm}+\text{sfc}}^{\mathrm{f}}$, and $[C_{\text{BRDF}}/T_{\mathrm{d}}{t}_{\mathrm{u}}]$, respectively.

Following similar NN training schemes as discussed previously, we conducted 10 000 simulations to determine the reflectance at aircraft altitude ${\mathit{\rho }}_{\mathrm{t},\text{atm}+\text{sfc}}$ from a system with only atmosphere and ocean surface (right panel of Fig. 1) and trained the NN for ${\mathit{\rho }}_{\mathrm{t},\text{atm}+\text{sfc}}$ in the same way as ${\mathit{\rho }}_{\mathrm{t}}^{\mathrm{f}}$. Since this system only includes the atmosphere and ocean surface but without the ocean body, there are a total of 14 input parameters (without Chl a). To train a NN for $[C_{\text{BRDF}}/T_{\mathrm{d}}{t}_{\mathrm{u}}]$ with $T_{\mathrm{d}}^{\mathrm{f},+}$, $t_{\mathrm{u}}^{\mathrm{f},+}$, and C_BRDF defined in Eqs. (12)–(14), we obtained five additional quantities corresponding to the abovementioned 10 000 cases with and without ocean body: for the fully coupled system with atmosphere, ocean surface, and ocean body (left panel of Fig. 1), we computed the reflectance just above and below the ocean surface (${\mathit{\rho }}_{\mathrm{t}}^{\mathrm{f},+}$ and ${\mathit{\rho }}_{\mathrm{t}}^{\mathrm{f},-}$) and irradiance transmittance just above and below the ocean surface ($T_{\mathrm{d}}^{\mathrm{f},+}$ and $T_{\mathrm{d}}^{\mathrm{f},\mathrm{0}}$); for the system without ocean body but with ocean surface (right panel of Fig. 1), we computed the reflectance just above the ocean surface (${\mathit{\rho }}_{\mathrm{t},\text{atm}+\text{sfc}}^{\mathrm{f},+}$). The accuracies of the NNs for ${\mathit{\rho }}_{\mathrm{t},\text{atm}+\text{sfc}}^{\mathrm{f}}$ and $[C_{\text{BRDF}}/t_{\mathrm{u}}{T}_{\mathrm{d}}]$ are evaluated and shown in Table 3, which has an accuracy around 1 % similar to other quantities.

To evaluate the overall accuracy for the R_rs after the BRDF correction, we conducted radiative transfer simulations with a zenith sun and a nadir viewing direction and obtained the truth remote sensing reflectance using the upwelling radiance and downwelling irradiance just above the ocean surface as the examples shown in Fig. 5. The predicted R_rs values were computed from Eq. (10) after the application of two NNs. The RMSEs of the difference between the simulated and NN-predicted R_rs are shown in Table 3, with a maximum value of 0.0004 at wavelength 440 nm and smaller than 0.0002 in other bands.

Figure 5Comparison of the truth R_rs (RT) and the neural network (NN) computed R_rs. The truth R_rs is computed from radiative transfer simulations with a zenith sun and nadir viewing direction. The NN computed R_rs follows Eq. (10).

Download

The NN forward models for reflectance (${\mathit{\rho }}_{\mathrm{t}}^{\mathrm{f}}$) and DoLP ($P_{\mathrm{t}}^{\mathrm{f}}$) are used in the FastMAPOL retrieval algorithm as discussed in Sect. 2. To evaluate the performance of the retrieval algorithm, we conducted retrievals on the synthetic AirHARP data. The creation of the synthetic data is discussed in Sect. 3.3. To account for the measurement uncertainties, random noise is added to the simulated data according to the calibration uncertainties as listed in Table 4. The total uncertainties in the cost function include contributions from calibration (σ_cal), radiative transfer simulation (σ_RT), and NN model (σ_NN).

Using the initial values as listed in Table 2, a total of 1000 synthetic AirHARP cases are retrieved with the cost function values (χ²) summarized in Fig. 6. Retrievals with χ²<1.5 are chosen in our following discussion, which includes 96 % of all retrieval cases. Gao et al. (2020) showed that the retrieval results depend on the initial values. Testing with several random sets of initial values shows that the statistics of the retrieval results from the 1000 synthetic cases are robust. As demonstrated by Di Noia et al. (2015) and Di Noia et al. (2017), a better choice of initial values for each pixel in the optimization may further improve the overall retrieval accuracy.

Figure 6Histogram of the cost function values (χ²) with initial values as specified in Table 3 with a total of 1000 cases. The most probable χ² is 0.82. A threshold of χ²<1.5 is used in the discussion.

Download

Figure 7The comparisons of the retrieved and truth values for total AOD (550 nm), SSA (550 nm), wind speed, and Chl a are shown in the top panels. The dashed line indicates the linear regression fitting with $y=\mathit{\beta }x+\mathit{\alpha }$, where β is the slope and α is the intercept. The lower panels show the difference between the retrieved and truth values of the corresponding upper panel parameters as a function of the total AOD at 550 nm.

Download

Figure 8The comparisons of the retrieved and truth values for the fine-mode aerosol parameters including AOD, SSA, refractive index (m_r), and effective radius (r_eff) and variance (v_eff).

Download

With the directly retrieved aerosol refractive index and volume densities (see Table 2) as inputs, the aerosol optical depth (AOD) and single scattering albedo (SSA) for both the fine and coarse modes were computed using additional NNs to represent the Lorenz–Mie calculations in Appendix A. The retrieved total AOD, SSA, wind speed, and Chl a are compared with the truth values as shown in Fig. 7. Total AOD indicates the summation of the fine- and coarse-mode AODs, and total SSA is the ratio of the total scattering and extinction cross sections; both are specified in Appendix A. For fine aerosol, the AOD, SSA, refractive index (m_r), and effective radius (r_eff) and variance (v_eff) are shown in Fig. 8. The color plots indicate the data point density (normalized by its maximum value) approximated by a kernel density estimation method (Silverman, 1986).

Figure 9The retrieval uncertainties at various aerosol loadings for AOD, SSA, refractive index (m_r), effective radius (r_eff) and variance (v_eff), wind speed, and Chl a. AOD values at the x axis from 0.1 to 0.5 indicate the five ranges of total AOD including [0.01, 0.1], [0.1, 0.2], [0.2, 0.3], [0.3, 0.4], and [0.4, 0.5], which are used to compute the corresponding uncertainties. Chl a uncertainties are evaluated in terms of MAE in log scale (see Eq. 22), and all other parameters are evaluated in terms of RMSE. AOD (%) indicates the percentage AOD uncertainties comparing to the truth AOD.

Download

In order to quantify the variation of the retrieval uncertainties with respect to different aerosol loadings, we computed the RMSE between the retrieved and truth values at five AOD ranges including [0.01, 0.1], [0.1, 0.2], [0.2, 0.3], [0.3, 0.4], and [0.4, 0.5]. Each AOD range includes approximately 200 cases. Note that as discussed in Sect. 3.3 the total AOD and the fine-mode volume fraction are uniformly sampled for the simulated data; therefore, there is an equal mixing fraction of fine- and coarse-mode aerosol for each AOD range. The retrieval uncertainties for aerosols are shown in Fig. 9, with the corresponding ranges indicated by AOD values from 0.1 to 0.5. All discussions regarding the AOD and SSA are for a wavelength of 550 nm in this section.

As shown in Figs. 7 and 9, the errors of the retrieved total AOD increase with aerosol loadings: the uncertainty (evaluated using RMSE) is 0.008 and 0.015 for the AOD range [0.01, 0.1] and [0.1, 0.2] and increases to 0.035 for the AOD range [0.4, 0.5]. Similar absolute uncertainties are found for both the fine- and coarse-mode AODs, with a slightly smaller value. In percentage, the total AOD uncertainties is 28.3 % at the AOD range [0.01, 0.1], where the large uncertainties are due to the cases with small AODs. For the AOD range from [0.1, 0.2] to [0.4, 0.5], the AOD uncertainties further decrease from 14.4 % to 5.6 %.

Similar to the total AOD uncertainties, the total SSA uncertainties decreases with AOD from 0.05 to 0.02. The fine-mode SSA uncertainties reduce similarly from 0.05 to 0.03. The uncertainties for coarse-mode SSA reduces slightly from 0.1 to 0.08, which is more than twice as large as the fine-mode SSA uncertainties. The uncertainties for the fine-mode m_r, r_eff, and v_eff show a larger value in the AOD bin of [0.01, 0.1] of 0.06, 0.024 µm, and 0.08 and then remain close to a constant at larger AOD ranges with a smaller value of 0.03, 0.01 µm, and 0.03 respectively. The averaged uncertainties for coarse-mode m_r, r_eff, and v_eff are approximately 0.08, 0.5 µm, and 0.15 respectively with less AOD dependency at all AOD ranges. The coarse-mode m_r uncertainty is more than twice the fine-mode uncertainty. The larger uncertainty values for coarse-mode r_eff and v_eff are also related to their large particle size.

For wind speed retrievals as shown in Fig. 7, the agreements between the truth and retrievals depend strongly on the wind speed value itself: when the wind speed is small, there is less retrieval sensitivity due to the removal of glint; for a higher wind speed, the agreements are improved, likely due to the larger range of angles influenced by wind speed. The retrieval uncertainties are shown in Fig. 9; for a wind speed (WS) lower than 3 m s⁻¹, the uncertainty increases from 1.5 to 2.1 m s⁻¹ for AOD ranges from [0.1, 0.2] to [0.4, 0.5]. For wind speed higher than 3 m s⁻¹, the averaged retrieval uncertainty is 1.2 m s⁻¹ with a small variation of less than 0.1 m s⁻¹.

The retrieved and truth Chl a is compared in Fig. 7, where the MAE in log scale is used with the definition

where

where R_i and T_i denote the retrieval and truth values. MAE(log) is recommended by Seegers et al. (2018) as a better metric for Chl a, which indicates the averaged ratio between the retrieval and truth values. The dependency of the MAE(log) for Chl a with the aerosol loadings is shown in Fig. 9. The Chl a retrieval performance depends on the magnitude of the Chl a. In this work, we chose four ranges of Chl a according to the trophic regions discussed in Seegers et al. (2018). Note that Chl a from in situ measurements is typically larger than 0.01 mg m⁻³, and we chose a broader range of Chl a with its minimum value of 0.001 mg m⁻³ as listed in Table 2 for sensitivity studies. For $\mathrm{0.01}\mathrm{mg}{\mathrm{m}}^{-\mathrm{3}}<\text{Chl}a<\mathrm{0.1}\mathrm{mg}{\mathrm{m}}^{-\mathrm{3}}$ and $\mathrm{0.1}\mathrm{mg}{\mathrm{m}}^{-\mathrm{3}}<\text{Chl }a<\mathrm{1}\mathrm{mg}{\mathrm{m}}^{-\mathrm{3}}$, Chl a retrieval uncertainties vary within 1.3 to 1.6 when AOD < 0.3 and then increase to 2.3 at AOD range [0.4, 0.5]. For $\text{Chl }a>\mathrm{1}\mathrm{mg}{\mathrm{m}}^{-\mathrm{3}}$ and $\text{Chl }a<\mathrm{0.01}\mathrm{mg}{\mathrm{m}}^{-\mathrm{3}}$, the uncertainties are generally larger, with a value around 2 to 3.

With the retrieved aerosol and ocean properties, the atmospheric correction procedures can be applied to compute the remote sensing reflectance as discussed in Sect. 3.4. The comparison of the retrieved R_rs with the truth data is shown in Fig. 5. To account for the various solar geometries, the BRDF correction has been applied to the retrieved R_rs as discussed in Sect. 3.4. Note that the ${\mathfrak{R}}_{\mathrm{o}}/\mathfrak{R}$ factor will not impact the BRDF correction in computing R_rs for the synthetic data, because of the small viewing zenith angles used at the four AirHARP bands, which are 1.22^∘, 1.17^∘, 0.03^∘, and 3.52^∘, respectively. The truth R_rs was computed with a zenith sun and a nadir viewing direction, emphasizing the need for the latter correction to the MAP observations. Overall R_rs uncertainties for the four bands are 0.007, 0.0004, 0.0002, and 0.0002 as shown by the RMSE in Fig. 10. MAE showed values of 0.0006, 0.0003, 0.0002, and 0.0001, which are less sensitive to outliers. Note that the atmospheric correction is applied to the synthetic measurements without adding additional random noise in order to evaluate the impacts on R_rs uncertainties from only aerosol and ocean surface property retrievals. The retrieval uncertainties for R_rs for each AirHARP bands are shown in Fig. 10 depending on the aerosol loadings: larger uncertainties are found with larger aerosol optical depth.

Figure 10The difference between the retrieval and truth R_rs with respect to AOD. The truth R_rs is computed with a zenith sun and a nadir viewing direction. The retrieved R_rs follows Eq. (10) with the BRDF correction considered. RMSE and MAE are for all retrievals cases at each wavelength.

Download

The PACE accuracy requirements on ocean color are specified in terms of the water-leaving reflectance, which can be converted to those of R_rs by dividing them by a factor of π. The resultant requirements in terms of R_rs are 0.0006 sr⁻¹ or 5 % from 400 to 600 nm and 0.0002 sr⁻¹ or 10 % from 600 to 710 nm (Werdell et al., 2019). As shown in Fig. 11, R_rs values at 550 nm are within the requirement of 0.0006 sr⁻¹ for all AOD ranges. For the 440 nm band, when AOD is less than 0.3, the R_rs retrieval uncertainties are less than 0.0005 sr⁻¹, but the uncertainties become as high as 0.001 sr⁻¹ at a larger AOD of 0.5. R_rs at 670 and 870 nm varies in a very small dynamical range and is less impacted by the aerosol retrievals. R_rs uncertainties at 670 and 870 nm are slightly larger than the requirement of 0.0002 sr⁻¹ when AOD (550 nm) is larger than 0.4 and 0.3 respectively. Further work is needed to understand how the uncertainties of the retrieved aerosol properties influence the retrievals. Note that from Table 3, the uncertainties of the R_rs computed using NNs have an uncertainty of 0.0004 to 0.0001 from 440 to 870 nm, which may be further minimized with better training and help to reduce the R_rs retrieval uncertainties.

Figure 11Retrieval uncertainties for R_rs at the four AirHARP bands. The uncertainties are computed in the same way as for Fig. 9 in terms of RMSE.

Download

The ACEPOL field campaign, conducted from October to November of 2017, included a total of six passive and active instruments on the NASA ER-2 high-altitude aircraft (Knobelspiesse et al., 2020) with four MAPs – AirHARP (McBride et al., 2020), AirMSPI (Diner et al., 2013), SPEX airborne (Smit et al., 2019), and the RSP (Cairns et al., 1999) – and two lidars – HSRL-2 (Burton et al., 2015) and CPL (the Cloud Physics Lidar) (McGill et al., 2002). Aerosol retrieval algorithms have been applied for all four MAPs (Fu et al., 2020; Puthukkudy et al., 2020; Gao et al., 2020). The measurement datasets are available from the ACEPOL data portal (Knobelspiesse et al., 2020). In this work, we focus on the study of the AirHARP measurements over ocean scenes as shown in Fig. 12 on 23 October 2017. The viewing angles are within $\pm \mathrm{57}{}^{\circ }$ along the track and $\pm \mathrm{47}{}^{\circ }$ across the track as shown in the polar plots in Fig. 12. Figure 13 shows the RGB images (670, 550, and 440 nm) for the three scenes at near-nadir viewing direction. AirHARP conducted high-spatial-resolution measurements with a grid size of 55 m and swath width of 42 km at nadir (up to 60 km at far angles). We averaged the reflectance and DoLP respectively within a bin box of 10×10 pixels (550 m×550 m).

Figure 12The location of the three ocean scenes from AirHARP from ACEPOL on 23 October 2017. The flight track color shows the UTC time along the flight path. The aircraft flew at an altitude of 20.1 km. The viewing zenith and relative azimuth angle (relative to the solar azimuth angle) for the 440 nm band from all pixels in the corresponding scene are shown in the bottom polar plots. The central portion of viewing angle plot is removed due to water condensation on the lens. The averaged solar zenith angles for the three scenes are 47.0^∘, 45.6^∘, and 52.9^∘, respectively, as indicated in the polar plots by the red asterisks.

Figure 13The RGB images (670, 550, and 440 nm bands) for the three scenes at near-nadir viewing directions. Scene 1 and Scene 2 observe only ocean, while Scene 3 observes both ocean and land. Sparse thin clouds are visible from Scene 1 and Scene 2. Sunglint can be observed at the lower portion of the Scene 2 image.

Download

The HSRL-2 instrument from ACEPOL provided useful aerosol optical depth ground truth at 355 and 532 nm (Hair et al., 2008; Burton et al., 2016), which was used for the validation of the aerosol retrieval algorithm using the AirHARP data. The HSRL-2 measures the pixels along the track as shown in Fig. 12, where an assumed lidar ratio of 40 sr is multiplied by the aerosol backscatter coefficient derived from the HSRL technique to produce aerosol extinction and AOD at 532 nm. For the low-aerosol-loading cases considered in this study, the assumed lidar ratio approach produces a systematic uncertainty of ±50 % (Fu et al., 2020). In Scene 3, the aircraft also flew over an AERONET USC_SEAPRISM site (33.564^∘ N, 118.118^∘ W) which is equipped with a CIMEL-based system called the Sea-viewing Wide Field-of-view Sensor (SeaWiFS) Photometer Revision for Incident Surface Measurements (SeaPRISM) that collects radiances at eight wavelengths of 412, 443, 490, 532, 550, 667, 870, and 1020 nm (Zibordi et al., 2009). AOD from the AERONET data product version 3 level 2.0 data was used in this study, which is also consistent with the HSRL-2 AOD at 532 nm as shown latter in Fig. 20. The estimated AERONET AOD uncertainty is from 0.01 to 0.02, with the maximum uncertainty in the UV channels (Giles et al., 2019). We compared AOD from AirHARP retrievals with those from both HSRL and AERONET. Furthermore, to validate the atmospheric correction procedure in the retrieval algorithm, we compared the retrieved remote sensing reflectance with the AERONET ocean color (OC) products as reported by the SeaPRISM measurements at the USC_SEAPRISM site.

Table 5The standard deviations of the reflectance (σ_ρ,avg) and DoLP (σ_P,avg) are calculated within a 10×10 pixel box and averaged over all angle and pixels from the three scenes (over ocean pixels only). The percentage values listed in the table indicate the percentage uncertainties.

Download Print Version | Download XLSX

To assess the spatial variability of field measurement, we computed the standard deviations for the reflectance (σ_ρ,avg) and DoLP (σ_P,avg) within the bin box. Representative values are provided in Table 5. The values of σ_ρ,avg and σ_P,avg at the 870 nm band are 4.5 % and 0.05, respectively, which suggest larger measurement uncertainties at 870 nm than other bands probably due to small radiometric magnitudes. Meanwhile, our retrieval tests showed larger coarse-mode retrieval uncertainties than synthetic data results. To better constraint retrievals, we assume the coarse-mode aerosol as sea salt by setting its imaginary refractive index to zero. All other retrieval parameter ranges are kept the same as in Table 2. Furthermore, we found our forward model cannot predict the angular variation of DoLP in the 440 nm band well (with an estimated MAE of 0.04), which contributes a major portion to the cost function and increases both fine- and coarse-mode retrieval uncertainties. Therefore, we exclude DoLP in the 440 nm band from our retrievals in this study.

Figure 14The histograms of the cost function values over the three scenes as shown in Fig. 12 with total pixel numbers of 13 491, 13 226, and 9159. Only pixels over ocean are considered.

Download

We conducted retrievals with similar procedures as discussed for synthetic data. The solar and viewing geometries as shown in Fig. 12 and the ozone column density ($n_{{\mathrm{O}}_{\mathrm{3}}}$) are assumed to be known inputs to the retrieval algorithm. The averaged values of $n_{{\mathrm{O}}_{\mathrm{3}}}$ from MERRA2 over each of the three scenes are obtained, which are 277.5, 278.6, and 281.3 DU respectively. The averaged surface pressures from MERRA2 over the three scenes are 1.008, 1.006 and 1.003 atm, which is consistent with our assumption in the atmosphere model as discussed in Sect. 3. The histograms of χ² for all pixels retrieved in each scene are shown in Fig. 14. The most probable χ² values are 1.1, 1.7, and 1.1 respectively.

Figure 15The number of viewing angles used in the retrieval (N_v), the cost function value (χ²), the retrieved AOD (550 nm), and R_rs (550 nm) for all pixels in Scene 1. The HSRL AODs at 532 nm are indicated by the colored dots in the AOD plot.

Download

Figure 16Same as Fig. 15 but for Scene 2. For R_rs, viewing angles at least 40^∘ away from the solar specular reflection direction are used to avoid sunglint as shown in Fig. 13.

Download

Figure 17Same as Fig. 15 but for Scene 3. The pixels with large χ² are mostly influenced by the land (upper region) and island (lower left). The retrieved AOD and R_rs over land pixels are not shown. The location of the AERONET USC_SEAPRISM site is indicated by a red star symbol.

To evaluate the retrieval performance, we plotted the map of the total number of viewing angles used in the retrieval (N_v), the cost function χ², the retrieved AOD (550 nm), and R_rs for each scene in Figs. 15–17. As discussed in Sect. 2, the maximum number of viewing angles is 120 for AirHARP measurements. Figures 15–17 show the number of available viewing angles vary from 0 to 120 due to the removal of glint and other data quality control measures. Discontinuity in the number of angles can be seen as a stripe, due to the removal of angles influenced by water condensation on the lens, which can also be observed in the polar plots in Fig. 12 with the nadir region removed. All three figures show that the number of viewing angles are smaller at the edges parallel to the flight track, where small χ² can be achieved but may be less reliable.

For pixels with large χ² as shown in Fig. 14, the forward model cannot fit the measured reflectance or DoLP well, which may be due to cloud contamination (Stap et al., 2015), land, or residuals of glint. In Scene 1, the top region with large χ² values is impacted by the thin cloud which is visible from Fig. 13. Larger residuals in the 870 nm band between measurement and forward model are also observed. The retrieved AOD are overestimated in this region. In Scene 2, the region with χ²>3 correlates closely with the thin clouds (Fig. 13), which influence nearby AOD and R_rs retrievals. For Scene 3, χ² become larger than 2 when close to the coast. This may be due to complex water properties which are not well represented by the open-water bio-optical model used in the simulation (Gao et al., 2019). The pixels near the coast are also potentially impacted by the bottom effect and adjacency effect of land pixels.

Figure 18Comparison of the retrieved AOD (550 nm) from AirHARP measurement with the AOD (532 nm) from HSRL for Scene 1 to 3. The AOD (550 nm) from the AERONET USC_SeaPRISM site is shown in Scene 3. The AirHARP retrieved AOD is averaged with 4×4 pixels (2.2 km×2.2 km). The averaged AODs and their standard deviations for both AirHARP retrievals and HSRL products are provided in the text. Pixels with N_v>30 and χ²<10 are considered.

Download

To compare with the HSRL AOD in the along-track direction, the retrieved AOD (550 nm) is averaged within a box of 4×4 pixels (2.2 km ×2.2 km). The averaged AOD (550 nm) values and the corresponding standard deviations are shown in Fig. 18. Pixels with N_v>30 and χ²<10 are considered. The overall averaged AODs and their standard deviation are also computed and indicated in the plots. The averaged HSRL AODs are 0.079, 0.071, and 0.037 for Scene 1 to 3. The averaged retrieved AOD (550 nm) are 0.096, 0.078, and 0.049, with relatively larger retrieval variation of 0.02 to 0.03. For Scene 1, most χ² values are larger than 2, while for the other two scenes, more pixels are less than 2 except for those pixels very close to cloud and coast. In Scene 1, the retrieved AOD is larger than that of the HSRL AOD by 0.03 in the overlapped region, which may be influenced by the remaining effect of water condensation. In Scene 2, the peaks of the retrieved AOD values correspond to the χ² values larger than 2, which are influenced by the nearby thin cloud. There are no overlapped pixels except the ones associated with high AOD peaks, but the general trend of the retrieved AOD agrees with the HSRL results. For Scene 3, the retrieved AOD values agree well with the HSRL AOD with an average difference of less than 0.01 and χ² mostly less than 2. However when the pixels are close to the coast, both χ² and AOD increased significantly as discussed previously.

Figure 19Similar to Fig. 18, the retrieved R_rs values are computed for the AirHARP band of 440, 550, and 670 nm bands. The averaged R_rs and its standard deviation are shown in the legends. For Scene 3, R_rs values from the AERONET USC_SeaPRISM site at wavelengths corresponding to AirHARP bands are indicated by the star symbols.

Download

Figure 19 shows the mean value and standard deviation of R_rs averaged in the same way as AOD discussed above. There is similar spatial variation between the retrieved R_rs and AOD. Pixels with large R_rs uncertainties are mostly associated with the large AOD uncertainties shown in Fig. 18. The R_rs values at 440 nm for the three scenes are 0.0055, 0.0072, and 0.0030 sr⁻¹, where the decrease in R_rs from Scene 2 to Scene 3 may be due to the increase in CDOM close to the coast as demonstrated in Fig. 5. Moreover, R_rs values at Scene 1 are likely to be underestimated due to the large χ² and retrieved AOD over the center of Scene 1. The averaged R_rs values at 550 nm remain approximately constant with a value of 0.0003 sr⁻¹ over all three scenes. R_rs values from the AERONET USC_SeaPRISM site are indicated in Scene 3 of Fig. 19 and also compared in Fig. 20. As discussed in Sect. 2.2, we chose the minimum viewing zenith angle available from the measurements after removing the sunglint. The removal of sunglint improves the R_rs calculation for Scene 2 as shown in Fig. 16. Moreover, we have ignored the ${\mathfrak{R}}_{\mathrm{0}}/\mathfrak{R}$ factor in Eq. (14), which may cause underestimation of the R_rs at the edge of the image where θ_v can reach as large as 60^∘. However it is challenging to analyze its impact at large θ_v angles. ${\mathfrak{R}}_{\mathrm{0}}/\mathfrak{R}$ has a strong dependency on wind speed, but the retrieved wind speeds from current retrievals show large uncertainties. Further work may require a better treatment of sunglint and improved accuracy in wind speed.

Figure 20Comparisons of the AOD and R_rs from AirHARP retrievals with AERONET products. The retrieval results are averaged with 4×4 pixels (2.2 km×2.2 km) around the AERONET USC_SeaPRISM site. This is similar to Figs. 18 and 19, with error bars indicating the standard deviations, but only pixels with χ²<2 are considered. The AERONET AOD and R_rs spectra are taken from 23 October 2017, with the error bars indicating daily variations. HSRL AOD at 532 nm is also shown.

Download

To better compare with AERONET results, we only considered the pixels with χ²<2 and conducted the same averaging (4×4 pixels) around the USC_SeaPRISM site for the retrieved AOD and R_rs. The averaged values and their standard deviations are plotted in Fig. 20. The overall retrieved AOD spectrum is similar to AERONET results with a difference smaller than 0.01. The results are similar to the retrieval results from the Research Scanning Polarimeter as reported by Gao et al. (2020). The retrieved R_rs agrees well with the AERONET R_rs, with a difference of less than 0.001 sr⁻¹. Note that this study is done with a possible AirHARP measurement uncertainty of 3 % in reflectance, which may impact the atmospheric correction accuracy.

The complete retrieval results, including the aerosol microphysical properties, wind speed, Chl a, and atmospheric-correction-related datasets, are provided in the “Data availability” section. The retrieval uncertainties for aerosol microphysical properties are relatively large due to the small aerosol optical depths. Chl a retrievals are sensitive to the aerosol retrievals and are more challenging to retrieve accurately at small values as discussed in Sect. 4.

The NN model greatly improved the speed of the forward model used in the iterative optimization approach and boosted the efficiency of the FastMAPOL retrievals. The average retrieval speed for one pixel with FastMAPOL is approximately 3 s with a single CPU core, or approximately 0.3 s with a GPU using the same hardware as mentioned in Sect. 3.3. Comparing to the retrieval speed of approximately 1 h per pixel using conventional radiative transfer forward model, the computational acceleration is 10³ times faster with CPU or 10⁴ times with GPU processors. Meanwhile, the NN model maintains a high accuracy as shown in Tables 3 and 4. For retrieval algorithms running the radiative transfer simulation, the accuracy is often tuned down to reduce the simulation time. By training a NN, however, the high accuracy of the radiative transfer model simulation can be achieved, as has been demonstrated in this work. Thus, despite the one-time high computational costs in generating the training datasets and conducting the training, the trained NN can be applied efficiently to large observational datasets in the retrieval algorithm.

In the retrieval of the AirHARP field measurement, each pixel has multiple viewing angles that are aggregated from measurements at different times with slightly different solar zenith angles. The NN used in FastMAPOL computes the polarimetric measurement for specific solar zenith angles for each viewing direction, and, therefore, the variation of the solar zenith angle can be captured. These effects may be small for AirHARP measurement in ACEPOL, with a maximum solar zenith angle difference within 0.6^∘. However, this capability can help to minimize the impacts of the solar angle for HARP2 in spaceborne measurements, which can reach to a maximum difference of 1.5^∘ for HARP2 observations.

With the efficient retrievals from FastMAPOL, we have discussed the retrieval performance and uncertainties for the aerosol properties, including AOD, SSA, refractive index, and particle sizes. Since the AirHARP measurements share many similar characteristics with HARP2 as planned for the PACE mission, the knowledge from the retrieval analysis can help to understand the retrieval performance for the HARP2 instrument in spaceborne measurements. Note that HARP2 is likely to have higher accuracy due to the onboard calibration capability and the potential to conduct cross calibration with the OCI instrument. For the development of the NN forward model for spaceborne measurements, similar training procedures can be applied with the sensor altitude at the top of the atmosphere instead of the aircraft altitude used in this study. Due to the impact of retrieval capability by geometry (Fougnie et al., 2020), solar and viewing geometries according to the PACE orbits need to be considered.

The water-leaving reflectance is obtained from the atmospheric correction process using the aerosol and ocean properties retrieved from the AirHARP measurements, and a similar procedure can be applied to future HARP2 retrievals. Since the hyperspectral OCI in PACE will provide high-accuracy measurements, the retrieved information can be applied to OCI and therefore assist hyperspectral atmosphere corrections as demonstrated by Gao et al. (2020) and Hannadige et al. (2021). However, aerosol retrieval and atmospheric correction are challenging in the UV spectral range (Remer et al., 2019a). For the ocean bio-optical model in this study, the water properties are modeled as open-ocean waters parameterized by a single Chl a value. For complex coastal water, complex bio-optical models are preferred in the retrieval of both accurate aerosol properties and water-leaving signals as demonstrated by Gao et al. (2019).

The data files for AirHARP and HSRL-2 used in this study are available from the ACEPOL website (https://doi.org/10.5067/SUBORBITAL/ACEPOL2017/DATA001, ACEPOL Science Team, 2017). Complete retrieval results as well as data related to atmospheric correction from the three AirHARP scenes discussed in Sect. 5 are available from NASA Open Data Portal: https://data.nasa.gov/Earth-Science/FastMAPOL_ACEPOL_AIRHARP_L2/8b9y-7rgh (Gao, 2021).

MG, BAF, KK, and PWZ formulated the original concept for this study. MG developed the NN model and generated the scientific data. PWZ developed the radiative transfer code. KK, PWZ, BC, OH, and YH advised on the uncertainty model and retrieval algorithm. BAF, AI, and PJW advised on the atmospheric correction. KK, AI, and YH advised on the NN model. JG supported the parallel computing of the training data. VM, XX, BM, and AP provided and advised on the AirHARP data. RF and SB provided and advised on the HSRL-2 data. All authors participated in reviewing and editing this paper.

The authors declare that they have no conflict of interest.

This research has been supported by the NASA PACE project, NASA (grant nos. 80NSSC18K0345 and 80NSSC20M0227), NASA (ACE and CALIPSO missions), SRON, NWO/NSO project ACEPOL (project no. ALW-GO/16-09), and the Radiation Sciences Program.

This paper was edited by Alexander Kokhanovsky and reviewed by Yingxi Shi and Reed Espinosa.