The impact and estimation of uncertainty correlation for multi-angle polarimetric remote sensing of aerosols and ocean color

Multi-angle polarimetric (MAP) measurements contain rich information for characterization of aerosol microphysical and optical properties that can be used to improve atmospheric correction in ocean color remote sensing. Advanced retrieval algorithms have been developed to obtain multiple geophysical parameters in the atmosphere–ocean system, although uncertainty correlation among measurements is generally ignored due to lack of knowledge on its strength and characterization. In this work, we provide a practical framework to evaluate the impact of the angular uncertainty correlation from retrieval results and a method to estimate correlation strength from retrieval fitting residuals. The Fast Multi-Angular Polarimetric Ocean coLor (FastMAPOL) retrieval algorithm, based on neural-network forward models, is used to conduct the retrievals and uncertainty quantification. In addition, we also discuss a flexible approach to include a correlated uncertainty model in the retrieval algorithm. The impact of angular correlation on retrieval uncertainties is discussed based on synthetic Airborne Hyper-Angular Rainbow Polarimeter (AirHARP) and Hyper-Angular Rainbow Polarimeter 2 (HARP2) measurements using a Monte Carlo uncertainty estimation method. Correlation properties are estimated using autocorrelation functions based on the fitting residuals from both synthetic AirHARP and HARP2 data and real AirHARP measurement, with the resulting angular correlation parameters found to be larger than 0.9 and 0.8 for reflectance and degree of linear polarization (DoLP), respectively, which correspond to correlation angles of 10 and 5^∘. Although this study focuses on angular correlation from HARP instruments, the methodology to study and quantify uncertainty correlation is also applicable to other instruments with angular, spectral, or spatial correlations and can help inform laboratory calibration and characterization of the instrument uncertainty structure.

1 Introduction

Satellite remote sensing is important for the study of the earth system at a global scale (National Academies of Sciences, Engineering, and Medicine, 2018). Remote sensing instruments are evolving rapidly, with increasing accuracy and spatial, spectral, and angular resolutions (Kokhanovsky et al., 2015; Dubovik et al., 2019). Multi-angle polarimeters (MAPs), measuring polarization states at multiple spectral bands and viewing angles, contain high information content for the study of aerosol and cloud optical and microphysical properties (Mishchenko and Travis, 1997; Chowdhary et al., 2001; Hasekamp and Landgraf, 2007; Knobelspiesse et al., 2012). The aerosol properties derived from MAP instruments can be used to assist atmospheric correction for ocean color remote sensing (Frouin et al., 2019; Gao et al., 2020; Hannadige et al., 2021).

Uncertainty quantification from MAP retrievals provides information on the quality of the data products and improves our understanding of retrieval sensitivities. These uncertainties depend on the retrieval algorithm as well as the instrument characterization, including the spectral bands, viewing angles, and polarization capability and the measurement accuracy. As shown in Fig. 1 and Table 1, MAP instruments collect a large number of high-quality measurements with differing numbers of spectral bands and viewing angles (Gao et al., 2021b). The number of spectral bands are mostly within 4–14 for MAP instruments. An exception is the SPEX sensors (SPEX airborne and SPEXone, Smit et al., 2019; Hasekamp et al., 2019), which acquire up to 400 spectral bands. The number of viewing angles for many instruments vary between 5 and 16 (e.g., 5 for SPEXone and 9 for SPEX airborne and MISR) but can be on the order of 100 for several hyper-angular instruments (90 for HARP2, 120 for AirHARP, 152 for RSP, and 250 for APS). More viewing angles are preferred for the observation of clouds (Waquet et al., 2009; McBride et al., 2020), and they can also be used to conduct multi-angle cloud masking and data screening to increase aerosol retrieval accuracy and coverage (Gao et al., 2021b).

Figure 1 Current MAPs in terms of the number of spectral bands and total number of viewing angles as summarized in Gao et al. (2021b). The bubble size for each instrument corresponds to the total number of measurements as indicated next to the instrument name. The full name and reference of the instruments are provided in Table 1. POLDER-1, 2, and 3 refer to the instruments on the ADEOS-I, ADEOS-II, and PARASOL missions. Note that MISR conducts multi-angle measurements without considering polarization.

Deschamps et al. (1994)Tanré et al. (2011)Fougnie et al. (2018)Diner et al. (1998)Diner et al. (2013)Diner et al. (2018)Cairns et al. (1999)Mishchenko et al. (2007)Martins et al. (2018)Martins et al. (2018)Smit et al. (2019)Hasekamp et al. (2019)

Table 1 Acronyms and their full name for the MAP instruments plotted in Fig. 1. The number of angles and spectral bands are summarized in Gao et al. (2021b). Details in the instrument characteristics are available in Dubovik et al. (2019).

To understand the retrieval uncertainties, an uncertainty model is required to describe the combined uncertainties from the MAP measurements, forward model, and a priori assumptions. These combined uncertainty sources are often assumed to be independent and without correlations; however, measurements with high angular or spectral resolution are likely to have correlated uncertainty, depending on instrument design. For example, a sensor may use the same detector to scan through all measurement view angles (e.g. the Research Scanning Polarimeter, RSP, Cairns et al., 1999), and thus the systematic errors due to calibration will be correlated for different measurement angles. The correlation property should be part of the MAP uncertainty model, but often it has not been sufficiently characterized. The characterization of measurement uncertainty correlation is affected by instrument calibration and data processing protocols, which is challenging to quantify. A sensitivity study considering the angular uncertainty correlation of the RSP data was conducted by Knobelspiesse et al. (2012), which showed that information content is affected by correlation strength. However, the actual impacts of uncertainty correlation in a retrieval algorithm have not been well explored, as this requires better understanding of the correlation characteristics and efficient implementation in a retrieval algorithm.

Retrieval algorithms that exploit correlation information in retrieval parameters and measurement uncertainties have shown benefits in improving remote sensing capabilities. The Generalized Retrieval of Aerosol and Surface Properties (GRASP) algorithm retrieves multiple pixels simultaneously while considering the spatial correlation of the retrieval parameters (Dubovik et al., 2014, 2021). Xu et al. (2019) developed a correlated multi-pixel inversion approach (CIMAP), which further considers the correlation between different retrieval parameters. Theys et al. (2021) developed a Covariance-Based Retrieval Algorithm (COBRA) based on an error covariance matrix estimated from measurements with spectral correlation, applied their approach to sulfur dioxide (SO2) retrievals from the TROPOspheric Monitoring Instrument (TROPOMI) data, and demonstrated improved retrieval performance. To accurately evaluate pixel-level uncertainty in ocean color retrievals, spectral correlations associated with the uncertainties in top-of-atmosphere reflectance are also accounted for OLCI (Lamquin et al., 2013) and MODIS (Zhang et al., 2022) in the uncertainty propagation.

In this study, we provide a practical framework to understand the measurement uncertainty structure, study the impact of correlation in MAP retrievals, and demonstrate the potential for improvement in geophysical retrieval performance when proper correlation information is incorporated into the retrieval algorithm. Angular uncertainty correlations in measurements from the AirHARP and HARP2 instruments are studied as examples. Both instruments measure 60 angles at 670 nm. AirHARP measures 20 angles for the 440, 550, and 870 nm bands, while HARP2 measures at 10 angles for these bands. Angular correlation within each band is considered and modeled separately. Two methods are used to evaluate the retrieval uncertainties under different correlation strengths: (1) the error propagation method is used to evaluate the optimal retrieval uncertainties, by mapping the input uncertainty model describing the total uncertainty of the measurement and forward model to the retrieval parameter domain, and (2) comparative analyses are performed between the retrieval results from synthetic MAP measurements and the “truth data” that were assumed in the generation of that synthetic MAP data. The Monte Carlo error propagation method (MCEP) is adopted to compare the retrieval uncertainties from these two methods. To efficiently conduct retrieval and uncertainty analysis, the FastMAPOL retrieval algorithm is employed in this study, which uses neural-network forward models for coupled atmosphere and ocean systems (Gao et al., 2021a). Analytical Jacobian matrices are derived based on the neural network and used to improve the efficiency of the retrieval (Gao et al., 2021b) and uncertainty quantification (Gao et al., 2022). To accurately evaluate the retrieval uncertainties of real measurements, an adaptive data screening approach is employed (Gao et al., 2021b). This ensures that only those measurements that can be sufficiently described by the forward model are used in this study, by avoiding uncharacterized uncertainty contributions due to contamination by cirrus clouds and other anomalies.

Furthermore, we study the angular uncertainty correlation in the measurements and demonstrate that the correlation property can be derived using the autocorrelation function from the retrieval fitting residuals. Studies on both synthetic data with various correlation strengths are conducted with results applied to the real measurement retrievals from AirHARP over multiple ocean scenes. Useful tools are provided to understand and analyze the angular correlated uncertainty structure and models. Note that autocorrelation analysis based on fitting residuals has been found useful in analyzing performance of machine learning algorithms such as using the Durbin–Watson test (Chatterjee and Simonoff, 2012).

In the following sections, we will discuss how to conveniently include an angular correlated uncertainty model in the retrieval algorithm (Sect. 2), evaluate the impact of correlated measurement uncertainty in retrieval uncertainties (Sect. 3), and estimate correlation strength by fitting residual analysis from both synthetic AirHARP and HARP2 data and real AirHARP measurements (Sect. 4). Discussion and conclusions are provided in Sect. 5. Although this study focuses on angular noise correlation, the conclusions on the impacts of correlations are also applicable to other instruments such as the hyperspectral measurements from both the SPEXone and the Ocean Color Instrument (OCI) that will be carried on NASA’s upcoming Plankton, Aerosol, Cloud, ocean Ecosystem (PACE) mission (Werdell et al., 2019).

2 Algorithm and methodology

2.1 FastMAPOL retrieval algorithm

In this study, the FastMAPOL algorithm is used to retrieve aerosol and ocean optical properties from HARP measurements. The algorithm includes three main components: (1) a set of neural-network-based radiative transfer forward models of the coupled atmosphere and ocean system (Gao et al., 2021a) and the corresponding analytical Jacobian matrix based on automatic differentiations derived from these neural networks (Gao et al., 2021b), (2) a multi-angle cloud masking and data screening module (Gao et al., 2021b), and (3) an efficient uncertainty quantification component (Gao et al., 2022). Water-leaving signals in terms of remote sensing reflectance are derived with an additional neural network trained for the atmospheric correction process (Gao et al., 2021a).

The neural-network forward models are trained for both reflectance and degree of linear polarization (DoLP) based on simulations from the successive orders of scattering radiative transfer model (RTSOS) developed by Zhai et al. (2009, 2010) and Zhai and Hu (2022). The atmosphere and ocean system is assumed to be a four-layer system, and the radiative transfer interactions among them are fully considered in the RTSOS model. The bottom layer is the ocean water body in which an open-ocean bio-optical model is used to parameterize scattering and absorption of ocean constituents based on the chlorophyll a concentration (Chl a; mg m⁻³) (Gao et al., 2019, 2021a). The second layer is the ocean surface with its roughness parameterized by wind speed through a scalar Cox–Munk model (Cox and Munk, 1954). The third layer is an aerosol layer mixed with Rayleigh scattering. This layer extends from ocean surface to a height of 2 km with a uniform aerosol vertical distribution. The last layer contains atmospheric molecules from 2 km to the top of the atmosphere. The US standard atmospheric constituent profile is used to describe the molecular distributions (Anderson et al., 1986).

A total of nine parameters are used to describe the aerosol microphysical properties. There are four parameters for the complex refractive index of fine and coarse mode. Aerosol size distributions are parameterized by five volume densities for five size submodes with fixed effective radius and variance (Dubovik et al., 2006; Xu et al., 2016; Gao et al., 2018). Absorption by atmospheric gases is considered in the RTSOS simulation, with ozone density as the only variable. The radiant path geometries are represented by the solar and viewing zenith angles and the viewing azimuth angle relative to the solar direction. Therefore, a total of 15 parameters are used as forward-model input, with 11 of them defined as retrievable parameters. Details of the parameter ranges are listed in Appendix A and discussed in Gao et al. (2021a, 2022). To represent the forward model accurately and efficiently, the neural-network (NN) architecture is optimized with an input layer of 15 parameters, followed by three hidden layers with 1024, 256, and 128 nodes and a final output layer with 4 nodes for each HARP band. The deep-learning Python library PyTorch is used for the training the NN (Paszke et al., 2019). The accuracy of the NN forward model is examined with an independent synthetic measurement dataset not used in training. An accuracy of better than 1 % for reflectance and better than 0.003 for DoLP has been achieved (Gao et al., 2021a). The uncertainties of the NN forward model are less than the instrument uncertainties of AirHARP and HARP2 (3 % in reflectance, 0.01 in DoLP for AirHARP, and 0.005 in DoLP for HARP2).

In this study, we will discuss the retrieval uncertainty and performance in aerosol properties, ocean surface wind speed, and Chl a in the ocean, as well as water-leaving signals based on the retrieval parameters. The water-leaving signal refers to the remote sensing reflectance (R_rs), which is the ratio of the upwelling water-leaving radiance and the downwelling solar irradiance just above the ocean surface (Mobley, 2022). R_rs can be estimated through the atmospheric correction process, which removes the contribution from the atmosphere and ocean surface from the total measurements at the sensor, and the additional BRDF (bidirectional reflectance distribution function) correction to reduce the dependency on the solar and viewing directions. Both atmospheric and BRDF corrections with their associated uncertainties are implemented using neural networks as discussed in Gao et al. (2021a, b) and followed by this study.

2.2 Retrieval cost function and uncertainty quantification

The maximum likelihood approach is used to retrieve the state parameters in FastMAPOL by minimizing a cost function that represents the difference between the measurements and the forward-model fitting (Rodgers, 2000)

$$\begin{array}{}\text{(1)}& {\mathit{\chi }}^{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{N}}{\mathit{y}}^T{\mathbf{S}}_{\mathit{ϵ}}^{-\mathrm{1}}\mathit{y},\end{array}$$

where $\mathit{y}=\mathit{m}-\mathit{f}\mathbf{\left(}\mathit{x}\mathbf{\right)}$ is the residual vector between measurement m and forward model f under retrieval parameters of x. Measurement vector m includes both reflectance (ρ_t) and DoLP (P_t), where the subscript t indicates the total signal measured by the instrument. The total number of measurements, N, at each pixel includes contributions from both reflectance and DoLP, which has been used in previous studies (Gao et al., 2021a, b).

The error covariance matrix S_ϵ in Eq. (1) specifies the uncertainties of each measurement and the correlation between different measurements at the same pixel, which is a symmetric matrix defined as

$$\begin{array}{}\text{(2)}& S_{\mathit{ϵ};i,j}=\mathbb{E}\left[\left(y_i-\mathbb{E}\left[y_i\right]\right)\left(y_j-\mathbb{E}\left[y_j\right]\right)\right],\end{array}$$

where i and j indicate the measurement at different angles and bands, and 𝔼 indicate the expectation values. To capture the angular uncertainty correlation, the autoregressive model of order of 1 (denoted as AR(1)) is used in the study of RSP data (Knobelspiesse et al., 2012) and adopted in this study for HARP data. AR(1) represents a linear Markov process with the error covariance matrix specified as

$$\begin{array}{}\text{(3)}& S_{\mathit{ϵ},i,j}=\left\{\begin{array}{ll}{\mathit{\sigma }}_{t,i}^{\mathrm{2}}& \text{if}i=j\\ {\mathit{\sigma }}_{c,i}{\mathit{\sigma }}_{c,j}{r}^{{\mathrm{\Delta }}_{\mathit{\theta }}|i-j|}& \text{if}i\ne j\text{but at the same band and}\\ & \text{polarization state}\\ \mathrm{0}& \text{otherwise},\end{array}\right.\end{array}$$

where σ_t is the total uncertainty, which includes both random noise and calibration uncertainty (σ_c). Only σ_c is assumed to be correlated between measurements at different viewing angles.

The ratios between random and calibration uncertainties may be different for reflectance and polarized signals (Knobelspiesse et al., 2019). The synthetic data are generated directly using the forward model; therefore, the contribution of forward-modeling uncertainty is not considered for the synthetic data study. Δ_θ is the average angular grid size, which depends on the spectral bands. We model the correlation properties using the Δ_θ estimated from the viewing angles in the along-track direction to better represent the stripe filter characteristics used to conduct HARP angular measurement. The averaged Δ_θ is approximately 6.0^∘ for AirHARP and 12^∘ for HARP2 at 440, 550, and 870 nm bands and 2.0^∘ for the 670 nm bands for both HARP instruments. r in Eq. (3) is the correlation parameter with a value between 0 and 1. For uncertainties with more complex structures, a general autoregression and moving average (ARMA) model can be used (Priestley, 1983; Brockwell and Davis, 1991). However, from our analysis based on the retrieval results from real AirHARP measurements, AR(1) works well for most cases. Detailed analysis can be found in Sect. 2.5 for theoretical basis, Sect. 4.3 for real data applications, and Sect. 5 for general discussions.

To better represent stronger correlations when it is close to 1, we define the correlation angles θ_c based on the correlation parameter r as

$$\begin{array}{}\text{(4)}& r^{{\mathrm{\Delta }}_{\mathit{\theta }}|i-j|}=e^{-{\mathrm{\Delta }}_{\mathit{\theta }}|i-j|/{\mathit{\theta }}_{\mathrm{c}}}.\end{array}$$

Therefore, θ_c indicates the angular range where magnitude in the correlation between angles is reduced by a factor of e. Similarly, correlation angles can be derived from r as

$$\begin{array}{}\text{(5)}& {\mathit{\theta }}_{\mathrm{c}}=-\mathrm{1}/\ln r.\end{array}$$

2.3 Uncertainty quantification

The pixel-wise retrieval uncertainty can be quantified by mapping the measurement and forward-model uncertainties into retrieval parameter space (Rodgers, 2000):

$$\begin{array}{}\text{(6)}& {\mathbf{S}}^{-\mathrm{1}}={\mathbf{K}}^T{\mathbf{S}}_{\mathit{ϵ}}^{-\mathrm{1}}\mathbf{K}+{{\mathbf{S}}_{\mathbf{a}}}^{-\mathrm{1}},\end{array}$$

where the Jacobian matrix K represents the partial derivatives of the measurements with respect to all the retrieval parameters. In this study, each retrieval parameter can only vary in a limited range, which imposes an implicit a priori constraint on the retrieval parameters. To capture its influence on retrieval uncertainties, we assume the a priori error matrix S_a in Eq. (6) to be diagonal with the a prior uncertainty for each state parameter approximated by its permitted range in retrievals (Gao et al., 2022). Both the parameter ranges and a priori values are listed in Table A1. The uncertainties are defined as the standard deviation (1σ) around the retrieval solution, which is estimated by the square roots of the diagonal elements of S. The uncertainties of variables which are a function of the retrieval parameters can also be derived from S and their derivatives. Due to the large number of retrieval parameters used in the retrieval, the evaluation of the retrieval uncertainties can be time consuming. The speed to compute uncertainties is improved using automatic differentiations based on neural-network forward models (Gao et al., 2022). For example, the uncertainty of remote sensing reflectance (R_rs) can be derived using the automatic differentiation applied on the neural networks for BRDF correction and atmospheric correction components as discussed in Appendix A from Gao et al. (2021b).

The retrieval uncertainties estimated by error propagation (hereafter called theoretical retrieval uncertainty) as shown in Eq. (6) represent the optimal scenarios, with limitations such as the assumption that the retrieval parameters successfully converged to the global minima (more discussions in Sayer et al., 2020; Gao et al., 2022). When the retrievals converge to a local minimum, both the retrieval results and associated Jacobians can be less representative of the truth values and therefore lead to inaccurate error propagation and uncertainty estimation. To quantify the retrieval uncertainties based on actual retrieval results, the retrieval errors are defined as the difference between the retrieval results and the truth from synthetic data, which are then used to compute the retrieval uncertainty (hereafter called real retrieval uncertainty). The comparison between the theoretical and real uncertainties is useful to assess the optimal and actual performance of a retrieval algorithm. The Monte Carlo error propagation (MCEP) method is used in this study to conduct such comparison (Gao et al., 2022). MCEP samples the retrieval errors from theoretical retrieval uncertainties and then directly compares the error distributions between theoretical estimation and real retrievals. This method provides additional flexibility in analyzing their statistics. Multiple sets of random samples are generated from the theoretical uncertainties with their variations analyzed, which therefore provides a way to evaluate the impact of sample size in estimating uncertainties (Gao et al., 2022). This method is used to quantify the retrieval uncertainties with various correlation strengths in the next section. A summary of the terminology used for the error and uncertainties for measurement and retrieval results is provided in Table 2.

Table 2 Error and uncertainty definitions.

2.4 Eigenvector decomposition on error covariance matrix

The error covariance matrix with non-diagonal terms is challenging to implement efficiently in optimization algorithms, which typically operate in diagonal space with no correlation between measurements. The error covariance matrix also creates barriers to understand the retrieval uncertainties, as the input uncertainties are not for a single measurement but rather related to multiple measurements. To overcome these issues, we convert the measurements into a new space where the error covariance matrix is diagonalized. Therefore, conventional optimization and error analysis techniques can be readily used.

To achieve this goal, eigenvector decomposition is applied on the error covariance matrix (Rodgers, 2000) as

$$\begin{array}{}\text{(7)}& {\mathbf{S}}_{\mathit{ϵ}}={\mathbf{U}}^T{\mathbf{D}}_{\mathit{ϵ}}\mathbf{U},\end{array}$$

where D_ϵ is a positive diagonal matrix defined by the eigenvalues of S_ϵ and U is a unitary matrix. Examples of how the eigenvalues vary with correlation strength and the associated impact of those variations on Shannon information content and retrieval uncertainties are discussed in Appendix B. Based on Eq. (7), the cost function in Eq. (1) and the error propagation in Eq. (6) can be written as

$$\begin{array}{}\text{(8)}& {\displaystyle }{\mathit{\chi }}^{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{N}}{{\mathit{y}}^{\mathbf{\prime }}}^T{\mathbf{D}}_{\mathit{ϵ}}^{-\mathrm{1}}{\mathit{y}}^{\mathbf{\prime }},\text{(9)}& {\displaystyle }{\mathbf{S}}^{-\mathrm{1}}={{\mathbf{K}}^{\prime }}^T{\mathbf{D}}_{\mathit{ϵ}}^{-\mathrm{1}}{\mathbf{K}}^{\prime }+{\mathbf{S}}_a^{-\mathrm{1}},\end{array}$$

where the original set of measurements, y, are converted to a new set of measurements, y^′, without any correlation, along with the respective Jacobian matrices:

$$\begin{array}{}\text{(10)}& {\displaystyle }{\mathit{y}}^{\prime }=\mathbf{U}\mathit{y},\text{(11)}& {\displaystyle }{\mathbf{K}}^{\prime }=\mathbf{UK}.\end{array}$$

Equations (8) and (9) are mathematical transformations that conveniently allow for working with diagonal matrices, with advantages and applications summarized below.

• A clear expression of measurement uncertainty. The diagonal terms in matrix D_ϵ represent the uncertainties for the new measurement y^′, therefore providing insights on the accuracy of the measurements impacted by correlation. An example of this is shown in Appendix B.

• Conducting minimization on the retrieval cost function. Equation (8) represents the cost function for non-correlated measurement y^′ and therefore can be used for conventional optimization algorithms such as the subspace trust-region interior reflective (STIR) algorithm (Branch et al., 1999) as used in current FastMAPOL algorithm (Gao et al., 2021a, b, 2022).

• Generating correlated errors. To study and visualize angular uncertainty correlations, correlated errors need to be generated and then added to the synthetic data. To achieve this goal, we generated the errors in the space for y^′ with random parameters sampled from a normal distribution assuming the eigenvalues in D_ϵ as its variance. These errors in y^′ are then transformed back to the original space y through $\mathit{y}={\mathbf{U}}^T{\mathit{y}}^{\mathbf{\prime }}$. To demonstrate how angular correlations impact the errors, the correlated error samples with a correlation angle of ${\mathit{\theta }}_{\mathrm{c}}=\mathrm{10}{}^{\circ }$ (r=0.9) and correlation angle of ${\mathit{\theta }}_{\mathrm{c}}=\mathrm{60}{}^{\circ }$ (r=0.98) are shown in Fig. 2a and c. A value of r=0.9 has been assumed in the study of RSP angular correlation by Knobelspiesse et al. (2012). With a larger θ_c the errors start to form a longer range of correlation with smoother variations. Note that the overall magnitude of the errors can vary within the full range as described by the calibration uncertainties. These errors are then added to the synthetic data and used to study retrieval. The fitting residuals from retrievals on the synthetic data with the added correlated errors are shown in Fig. 2c and d and discussed in the next section.

Figure 2 Examples of simulated measurement errors generated for reflectance at the 660 nm band with a correlation angle of 10^∘ (a) and 60^∘ (b), respectively. The fitting residuals are shown in panels (c) and (d). A total of 1000 sets of errors are generated and added to the simulation data. Three sets of error examples are shown in different colors. The right side of viewing angle ends around 25^∘ due to the removal of sunglint as shown in Fig. 4.

2.5 Correlation strength estimation using autocorrelation

Autocorrelation is a useful function to quantify correlation in a discrete data sequence and is defined as (Priestley, 1983)

$$\begin{array}{}\text{(12)}& R_{i,j}=\mathbb{E}\left[y_i{y}_j\right],\end{array}$$

where i and j are two indices of the datasets. Comparing with Eq. (2), the autocorrelation is equivalent to the autocovariance when 𝔼[y_i]=0. This method can be applied to the simulated noise generated in Sect. 2.4 and used to analyze the fitting residuals. However, the mean values and variance in the fitting residuals often vary with respect to the angular grids. This type of signal is classified as non-stationary and difficult to study by the AR models (Priestley, 1983). To overcome this issue, the original residual data y are processed by removing their mean and normalizing by their standard deviation. These normalized data are denoted as $\stackrel{\mathrm{̃}}{\mathit{y}}$. For the data within the same band and polarization state, the autocorrelation function on the normalized data is equal to their covariance as defined in Eq. (3),

$$\begin{array}{}\text{(13)}& {\stackrel{\mathrm{̃}}{R}}_k=\mathbb{E}\left[{\stackrel{\mathrm{̃}}{y}}_i{\stackrel{\mathrm{̃}}{y}}_{i+k}\right]=r^{{\mathrm{\Delta }}_{\mathit{\theta }}k}.\end{array}$$

We can estimate the correlation by analyzing the residuals between the measurement and forward model. The autocorrelation function is averaged over multiple pixels to reduce uncertainties for the analysis in both synthetic data and real retrieval residuals. The correlation parameter can then be derived as

$$\begin{array}{}\text{(14)}& r=({\stackrel{\mathrm{̃}}{R}}_{\mathrm{1}}{)}^{\mathrm{1}/{\mathrm{\Delta }}_{\mathit{\theta }}}.\end{array}$$

Correlation angles, θ_c, are then computed based on r following the formula in Eq. (5). Furthermore, a partial autocorrelation function from a sequence of data can be computed, which removes the correlation due to lags higher than 1 (Priestley, 1983). If the AR(1) model is sufficient to describe the noise structure, only one additional term would be left besides the zero-order term in the partial autocorrelation results. Therefore, partial autocorrelation can be used to validate our assumption in the uncertainty model. The Python packages, statsmodels (Seabold and Perktold, 2010) and SciPy (Virtanen et al., 2018), are used to conduct the autocorrelation analysis.

An example is shown in Fig. 3; the autocorrelation function and partial autocorrelation function are applied on the simulated errors and the retrieval residuals with examples from Fig. 2b and d. The autocorrelations are shown in Fig. 3a and c for the simulated errors and retrieval residuals in reflectance data. The partial autocorrelations for Fig. 3a and c are shown in Fig. 3b and d, respectively. For both cases in Fig. 3b and d, only the first order of data are prominent, which confirms that the data can be represented by the AR(1) process. If higher orders in the AR process are presented, more prominent data points will appear in Fig. 3b and d. The estimated correlation angles for the errors in Fig. 3a and b and residuals in Fig. 3c and d are approximately 30 and 15^∘, respectively, after converting correlation parameters to correlation angles but less than the actual correlation angle of 60^∘. The results show that autocorrelation can be a useful way to estimate correlation strength but with a tendency to underestimate due to the finite length of the data and overfitting of the retrievals (more discussion in the Sect. 4). Note that the Bartlett confidence interval that corresponds the autocorrelation of uncorrelated white noise can be calculated but is found to be much smaller than the results in Fig. 3 and is therefore not shown in the figure (Brockwell and Davis, 1991).

Figure 3 The autocorrelation and partial autocorrelation on the 1000 sets of simulated errors with a correlation angle of 60^∘ and its corresponding fitting residuals as shown in Fig. 2b and d for the 670 nm band. The maximum value in the curve is normalized to 1. The horizontal axis indicates the angular step k as defined in Eq. (13). The Bartlett confidence interval on white noise is not shown due to its small value.

3 Retrieval uncertainties for AirHARP and HARP2

3.1 Synthetic data generated using the NN forward model

The neural-network forward model discussed in Sect. 2.1 is used to generate 1000 sets of synthetic AirHARP data, and then the number of viewing angles at 440, 550 and 870 nm are downsampled to 10 to represent HARP2 data. A fixed solar zenith angle of 50^∘ is used to represent the solar geometries of the AirHARP scenes over ocean from the ACEPOL field campaign (more information in Sect. 4.3). The aerosol properties, wind speed, and Chl a values are randomly sampled based on their allowed range, as discussed in Sect. 2 and Appendix A. The same sampling approach discussed in Gao et al. (2022) is conducted assuming that the aerosol optical depth (AOD) and fine-mode volume fraction are uniformly distributed within [0.01,0.5] and [0,1], respectively. A larger range of AOD values will be needed for applying this algorithm to cases of smoke and plume events. Realistic HARP-like viewing geometries are constructed by sampling the along-track and cross-track viewing angles randomly and then converting to the actual viewing zenith and azimuth angles following the formulas provided in Gao et al. (2021b). Example viewing geometries for AirHARP and HARP2 for bands at 550 nm are provided in Fig. 4, with geometries at other bands constructed similarly.

Figure 4 Example viewing zenith (θ_v) and relative azimuth (ϕ_v) in a polar plot at 550 nm for AirHARP (20 angles) and HARP2 (10 angles). Five sets of examples are provided in different colored dots. The 440 and 870 nm bands are similar. At 670 nm, there is a total of 60 angles for both instruments. The anti-solar point is indicated by the red asterisk. The viewing angles within the sunglint region indicated by the blue shaded area are removed. Note that the viewing angles from HARP2 are downsampled from AirHARP.

To generate realistic measurements, correlated uncertainties with a correlation angle (θ_c) from 0, 1, 2, 5, 10, 20, 30, 60, and 120^∘ are considered in this study, with corresponding correlation parameters of 0, 0.368, 0.607, 0.819, 0.905, 0.951, 0.967, 0.983, and 0.992, respectively. The corresponding correlated error samples are generated based on the error covariance matrix using the method discussed in Sect. 2.4. Examples of the correlated errors are already shown in Fig. 2 for AirHARP at 660 nm, with correlation angles of 10 and 60^∘, respectively.

Correlated errors for both the AirHARP and HARP2 instruments are generated according to the same 3 % uncertainty for reflectance but 0.01 in DoLP for AirHARP and 0.005 in DoLP for HARP2. These errors are added to the corresponding simulated reflectance and DoLP ($\sqrt{Q^{\mathrm{2}}+U^{\mathrm{2}}}/I$) for further studies. The reflectance is more likely to be dominated by systematic uncertainty like calibration (possibly correlated), while DoLP defined as the ratio between two measurements is more likely dominated by randomly generated uncertainty like shot noise (probably less correlated) (Knobelspiesse et al., 2012). Therefore, we assume two scenarios: (1) angular correlation only existed in reflectance measurement, not in DoLP measurement, and (2) both reflectance and DoLP have angular correlations with the same strength. Since the actual amount of correlation is not known, we designed our studies with the assumed correlation in the synthetic measurement but with either no information or full information on the correlation angle in the retrieval cost functions. Therefore, four scenarios are discussed in this study as summarized in Table 3 denoted by C1 to C4. We will discuss whether better retrieval results can be obtained if accurate correlation angles are considered in the retrieval cost function and whether we can estimate correlation from retrieval residuals in Sect. 4.

Table 3 Four scenarios of simulated uncertainties are considered in the synthetic data and retrievals. C1 and C2 indicate simulated errors with a correlation angle of θ_c added to synthetic reflectance data. C1 assumes the correlation property is unknown and no correlation is considered in the cost function model, but C2 assumes the correlation as in the simulated errors is known. This is also the case for C3 and C4, which considered correlated uncertainty in both reflectance and DoLP.

3.2 Retrieval uncertainties impacted by uncertainty correlation

Using the Monte Carlo error propagation (MCEP) method discussed in Gao et al. (2022), we compared both real and theoretical retrieval uncertainties with different correlation angle and testing scenarios as summarized in Table 3 for both synthetic AirHARP and HARP2 measurements. Figure 5 demonstrates the basic approach in MCEP by comparing the AirHARP retrieval results with ${\mathit{\theta }}_{\mathrm{c}}=\mathrm{60}{}^{\circ }$ for scenarios C3 and C4, where the real retrieval errors are sampled based on their retrieval uncertainties. The real uncertainties in both the root mean square error (RMSE) and the mean average error (MAE) are mostly larger when uncertainty is correlated (comparing b and a), with exceptions possibly due to statistical fluctuations in the Monte Carlo sampling. The theoretical uncertainties are similar because in both cases the correlation angles are assumed to be zero. After considering the same correlation angle in the retrieval cost function model as shown in Fig. 5c, both theoretical and real uncertainties are reduced. The real and theoretical uncertainties are similar to each other as shown in Fig. 5a, but the agreement degrades when correlation is considered (Fig. 5b, c). There are more real retrieval errors in the negative side as shown in Fig. 5c. This may be related to the convergence for cases with small AOD values as discussed in Gao et al. (2022). Although real uncertainties are generally larger than theoretical uncertainties, the differences are mostly associated with cases where AODs are small and underestimated relative to the truth values.

Figure 5 The histogram of the AirHARP retrieval errors for the fine-mode AOD from theoretical and real uncertainty estimations based on the MCEP method. (a) For the case without any correlation in the uncertainty; (b) for the cases with correlation only in the synthetic data (C3) and (c) also in the retrieval cost function (C4) with a correlation angle of ${\mathit{\theta }}_{\mathrm{c}}=\mathrm{60}{}^{\circ }$. Both mean average error (MAE) and root mean square error (RMSE) are computed for the theoretical and real errors indicated in red and blue text, respectively.

Both real errors and theoretical uncertainties have occasional outliers with large values possibly due to convergence to local minima instead of global minima, and this has large impacts on the RMSE values. For example, in Fig. 5b there are a few large estimated uncertainty values outside of the plot range that result in large RMSE for the theoretical error, almost comparable to the real uncertainty. Comparing MAE, however, shows that the theoretical value is much smaller than the real error, consistent with the histogram shape where the theoretical curve is narrower than the real curve. Given its robustness to outliers, MAE is used as the primary metric in this study.

Based on the MCEP method, we analyzed the retrieval uncertainties for synthetic AirHARP measurements from real errors and theoretical errors for various properties in Fig. 6, including the AOD, single scattering albedo (SSA), the real part of the refractive index (m_r), effective radius (r_eff), and effective variance (v_eff) for both fine and coarse modes and their combinations. fvf is the fine-mode volume fraction. Ocean surface wind speed and Chl a concentration are also retrieved. As recommended by Seegers et al. (2018), the uncertainty of Chl a is represented by a log-transformed metric:

$$\begin{array}{}\text{(15)}& {\displaystyle }\text{MAE}\left(\log \right)={\mathrm{10}}^Y\text{(16)}& {\displaystyle }Y={\displaystyle \frac{\mathrm{1}}{M}}\sum _i^M\left[{\log }_{\mathrm{10}}\left(\text{Chl}{a}_{i,\text{retrieval}}\right)-{\log }_{\mathrm{10}}\left(\text{Chl}{a}_{i,\text{truth}}\right)\right],\end{array}$$

where M is the total number of samples, which equals to the total number of synthetic measurement cases. The values of ${\log }_{\mathrm{10}}\left(\text{Chl}{a}_{i,\text{retrieval}}\right)-{\log }_{\mathrm{10}}\left(\text{Chl}{a}_{i,\text{truth}}\right)$ are sampled for both the real and theoretical uncertainties similar to Fig. 5, with detailed discussions provided in Gao et al. (2022). Furthermore, remote sensing reflectance at the four HARP wavelengths is derived after conducting the atmospheric correction using the retrieved aerosol and ocean properties (Gao et al., 2021a). Multiple sets of random errors are sampled and averaged to estimate the impact of sample size in estimating retrieval uncertainties and are shown in Fig. 6.

Figure 6 Retrieval uncertainties averaged for the AOD within [0.01,0.5] based on the MCEP method for synthetic AirHARP measurements. Ten sets of random samples are conducted, with the mean values shown in the plot and standard deviation indicated by the width of the shaded lines. The blue lines indicate the real uncertainty for Scenario C3 and Scenario C4 (Table 3), and the red line indicates the theoretical uncertainty for Scenario C4.

The real retrieval uncertainties for Scenario C3, in which correlation is considered in the simulated errors but not in the retrieval cost function, are found to be always increasing with the correlation angle. When the correct correlation angle is considered in the retrieval cost function (C4), the real retrieval uncertainty increases until θ_c, reaching 10 to 20^∘, and then slightly decreases at higher correlation angles, for most retrieval parameters. Similar behavior of the information content has been reported by Knobelspiesse et al. (2012) on the study of error correlation in RSP measurements. To understand how the correlation influences retrieval accuracy, we further analyze its impacts on the eigenvalues of the error covariance matrix with details discussed in Appendix B. The theoretical uncertainties for Scenario C4 are similar to the real uncertainties for ${\mathit{\theta }}_{\mathrm{c}}<\mathrm{10}{}^{\circ }$ except for refractive index, for which the theoretical values are almost half that of the real uncertainties. When ${\mathit{\theta }}_{\mathrm{c}}>\mathrm{10}{}^{\circ }$, the theoretical uncertainties then decrease much faster than the real uncertainties as θ_c increases. The difference mostly comes from the retrieval cases which underestimate the truth values as shown in Fig. 5. Note that the real retrieval uncertainties are mostly larger than the theoretical uncertainties without any correlation within the range of ${\mathit{\theta }}_{\mathrm{c}}<\mathrm{120}{}^{\circ }$. However, theoretical uncertainties predict that the retrieval uncertainties increase slightly with θ_c around 10^∘ (or r=0.9) and then decrease. Its values can be smaller than those with zero correlation when θ_c is larger than around 20^∘.

The results for HARP2 are similar to that of AirHARP as shown in Fig. 7. The overall HARP2 retrieval uncertainties are slightly smaller than AirHARP retrieval uncertainties, but the difference is mostly within 20 % for the theoretical uncertainties and mostly within 10 % for the real uncertainties as shown in Fig. 7 for both ${\mathit{\theta }}_{\mathrm{c}}=\mathrm{60}{}^{\circ }$ and ${\mathit{\theta }}_{\mathrm{c}}=\mathrm{60}{}^{\circ }$. Although HARP2 measures fewer viewing angles at 440, 550, and 870 nm bands, its better DoLP accuracy still results in slightly smaller uncertainties in most cases. Note that the retrieval accuracies also depend on the total number of viewing angles used and the range of scattering angles as discussed in Gao et al. (2021b).

Figure 7 The ratio of the real (R) and theoretical (T) retrieval uncertainties between the AirHARP and HARP2 measurements for correlation angles equal to 0 and 60^∘. In the legend, $R\left(\text{HARP2}\right)/R\left(\text{AirHARP}\right)$ denotes the ratio of real retrieval uncertainties between HARP2 and AirHARP, and $T\left(\text{HARP2}\right)/T\left(\text{AirHARP}\right)$ denotes the ratio of theoretical uncertainties.

To understand more quantitatively how the correlation angle impacts both AirHARP and HARP2 retrievals, the ratios of the real uncertainties between scenarios C3 and C4 are represented as $R\left(\text{C4}\right)/R\left(\text{C3}\right)$ for each retrieval quantity as shown in Fig. 8. Note that both C3 and C4 considered uncertainty correlations in both reflectance and DoLP, but only C4 considered the same correlation in its retrieval cost function, and C3 assumes no correlation in its retrievals. Therefore, the ratio $R\left(\text{C4}\right)/R\left(\text{C3}\right)$ represents how much retrieval uncertainties can be reduced if the correct amount of correlations are known in the retrieval process. As shown in Fig. 8, the ratio is close to 1 for most parameters for ${\mathit{\theta }}_{\mathrm{c}}=\mathrm{10}{}^{\circ }$, with slightly larger impacts for effective variance for both fine and coarse mode as well as for wind speed uncertainty. Such ratios almost double for ${\mathit{\theta }}_{\mathrm{c}}=\mathrm{60}{}^{\circ }$ for both AirHARP and HARP2 measurements. Note that sunglint has been removed in the retrieval as shown in Fig. 4, which may contribute to a larger wind speed retrieval uncertainty. When comparing the theoretical uncertainty T(C4) with R(C3), their ratio reduces to a value of 0.5 to 0.7 for ${\mathit{\theta }}_{\mathrm{c}}=\mathrm{10}{}^{\circ }$ for aerosol properties and further decreases to a value of 0.3 to 0.5 with ${\mathit{\theta }}_{\mathrm{c}}=\mathrm{60}{}^{\circ }$. The impacts on the remote sensing reflectance are generally smaller for ${\mathit{\theta }}_{\mathrm{c}}=\mathrm{10}{}^{\circ }$ but more significant for ${\mathit{\theta }}_{\mathrm{c}}=\mathrm{60}{}^{\circ }$. This is because the theoretical uncertainty with a correlation angle decreases faster than the real uncertainties (see Fig. 6). These results demonstrate that there is potentially more room to reduce the retrieval uncertainty as predicted by the error propagation theory represented by Eq. (6), which warrants future investigations. The largest gap between the real and theoretical uncertainties is still for refractive index in both fine mode and coarse mode, which is the same as observed in Gao et al. (2022). The gaps between the real and theoretical uncertainties increase with increasing θ_c, which indicates degrading retrieval performance in the presence of correlated uncertainties.

Figure 8 The ratio of real retrieval uncertainties between scenarios C4 and C3 as denoted in the legend by $R\left(\text{C4}\right)/R\left(\text{C3}\right)$ for AirHARP and HARP2 measurements with a correlation angle equal to 10 and 60^∘; similarly, the real and theoretical uncertainties for C4 are compared with the real uncertainties for C3 denoted as $T\left(\text{C4}\right)/R\left(\text{C3}\right)$.

To understand retrieval performance when the uncertainty correlation is only in reflectance for scenarios C1 and C2, similar ratios to Fig. 8 are plotted in Fig. 9. The ratios between the real uncertainties are almost always equal to a value of 1, with slightly larger impacts (less than 10 %) for ${\mathit{\theta }}_{\mathrm{c}}=\mathrm{60}{}^{\circ }$. However, the theoretical errors are much smaller, with values ranging from approximately 0.75 for ${\mathit{\theta }}_{\mathrm{c}}=\mathrm{10}{}^{\circ }$ to 0.5 for ${\mathit{\theta }}_{\mathrm{c}}=\mathrm{60}{}^{\circ }$ and with even smaller values observed for refractive index. This suggest that correlation in reflectance alone has the potential to be harvested to improve retrieval performance but is even harder to realize in real retrievals from current algorithms.

Figure 9 Similar to Fig. 8 but for the scenarios of C1 and C2 when correlations are only in reflectance uncertainties.

4 Estimating correlation from residual analysis

4.1 Cost function and fitting performance

As discussed in Gao et al. (2021b), the retrieval cost function can be characterized well by the χ² distribution for synthetic data and for real data after removing anomalies such as cirrus cloud. Note that the cost function in Eq. (1) is called a χ² function, but its histogram may not always follow a χ² statistical distribution, which depends on how well the fitting residuals can represent the real uncertainty and the degree of freedom of the χ² distribution (Gao et al., 2021a). Therefore, the histogram of successful retrievals is a useful indicator on how well the retrieval residuals compare with the assumed input uncertainty model (Rodgers, 2000). In this study, we found that the retrieval residuals can be represented well by the χ² distribution with a degree of freedom equal to the total number of measurements used (N), which is twice that of the total number of viewing angles (N_v) when there is no correlation, as shown in Fig. 10. Note that the degree of freedom is defined for the χ² distribution (Gao et al., 2021b).

Figure 10 Cost function histogram with a correlation angle of 0, 10, and 60^∘ for scenarios C3 (a) and C4 (b). The average number of viewing angles from the 1000 simulation cases, after removing sunglint, is N_v=75 for the AirHARP measurement. The red line indicates the χ² distribution with a degree of freedom of 2N_v=150. The green and blue lines indicate the χ² distribution with a reduced degree of freedom of 40 and 20 fitted to the corresponding histogram.

However, the cost function histogram shifts to a smaller value for Scenario C3 (Fig. 10a), where correlated uncertainty is included in the synthetic data but no correlation is considered in the retrieval cost function. Smaller cost function values indicate smaller retrieval residual, which may be caused by overfitting of the data and also possibly lead to the larger real retrieval uncertainties as shown in Fig. 6. When the correct value of correlation is considered in the retrieval cost function as shown in Fig. 10b for Scenario C4, the histogram better approaches a χ² distribution. For example, in the right panel in Fig. 10, the χ² distributions with degrees of freedom of 40 and 20 is found to better fit the cost function histogram with ${\mathit{\theta }}_{\mathrm{c}}=\mathrm{10}{}^{\circ }$ and ${\mathit{\theta }}_{\mathrm{c}}=\mathrm{60}{}^{\circ }$, which compares to results using all the measurement degrees of freedom (150). This suggests that the correlation in the uncertainty reduces its degree of freedom.

To understand how much overfitting impacts retrievals under a different strength of correlations, we compare the standard deviation of the retrieval residuals with the original simulated uncertainty for all 1000 cases. For reflectance, we consider the ratio between the simulated uncertainty and the reflectance as a convenient way to compare with the 3 % (or 0.03) uncertainty model for reflectance (Fig. 11). For the simulated uncertainty, the value of 0.03 is confirmed; however, the retrieval residuals become smaller with increasing θ_c for C3. The 670 nm band shows the largest reduction, where a ratio of 0.025 is observed for ${\mathit{\theta }}_{\mathrm{c}}=\mathrm{10}{}^{\circ }$, reaching 0.015 for ${\mathit{\theta }}_{\mathrm{c}}=\mathrm{120}{}^{\circ }$, which may be due to the large number of angles (60) in the 670 nm bands versus the other bands (10 or 20). This behavior indicates overfitting, where the uncertainties are partially removed as real signals and lead to reduced residuals. When the correct correlation is considered in the model, the result is much like the assumption, with a very slight indication of overfitting (within 0.005 in most cases). Similarly for DoLP, the standard deviations of the fitting residuals at the 670 nm band are smaller than other bands for C3, which indicate larger overfitting at the 670 nm band (Fig. 11c). The overfitting is much reduced for C4 with a standard deviation close to 0.01 (Fig. 11d). The results for HARP2 are very similar, but with the assumed uncertainty for DoLP of 0.005 (not shown).

Figure 11 The standard deviation of simulated measurement errors and fitting residual for cases C3 and C4. The model uncertainty of 3 % for reflectance and 0.01 for DoLP are indicated by dashed lines. Results of scenarios C1 and C2 for reflectance are similar, but with estimated DoLP uncertainties closer to the 0.01 line.

4.2 Correlation estimation results with synthetic data

Different amount of over fitting also partially removed the correlation in the fitting residuals. As shown in Fig. 12a and b, the estimated correlation is smaller than the truth. The estimated θ_c for 670 nm bands are smaller than other bands, which is consistent with previous studies where these bands overfit the data the most. Green bands seem to show the best results, which estimate ${\mathit{\theta }}_{\mathrm{c}}<\mathrm{10}{}^{\circ }$ (r=0.90) well but underestimate true ${\mathit{\theta }}_{\mathrm{c}}=\mathrm{60}{}^{\circ }$ (r=0.983) as 20^∘ (r=0.951) and ${\mathit{\theta }}_{\mathrm{c}}=\mathrm{120}{}^{\circ }$ (r=0.992) as 30^∘ (r=0.967). Note that the value of correlation parameter r asymptotically approaches 1 and becomes harder to distinguish with a finite length of data. Meanwhile, we also computed the correlation angles from the added simulated errors, which also resulted in a smaller correlation angle compared to the truth but much closer than using the residual data. This may due to the finite length of the measurement, which is 90^∘ after the glint is removed (total 120^∘). When the true correlation angle is considered in the model (C4), the estimated correlation angle is improved. The results of DoLP are slightly better where the retrieval residuals seem to estimate the results from simulated errors in a similar way for both C3 and C4 as shown in Fig. 12c and d. In real data, the correlation strength may be different for reflectance and DoLP; it would be useful to estimate the correlation angles for all the four bands and both reflectance and DoLP and analyze the difference by comparing it with the synthetic data study.

Figure 12 The estimated correlation angle θ_c for scenarios C3 and C4 from the simulated measurement errors and fitting residuals for reflectance and DoLP are compared with the truth values. The dashed line indicates the one-to-one line. Results of scenarios C1 and C2 for reflectance are similar, but with the estimated correlation angle close to zero since no correlation is considered in DoLP (not shown).

4.3 Correlation estimation results with real AirHARP data

The Aerosol Characterization from Polarimeter and Lidar (ACEPOL) field campaign was conducted from October to November of 2017 with the NASA's ER-2 aircraft at a high altitude of approximately 20 km (Knobelspiesse et al., 2020). Measurements over a variety of scenes are conducted from four MAPs – AirHARP, AirMSPI, SPEX airborne, and RSP – and two lidar sensors – HSRL-2 (Burton et al., 2015) and CPL (the Cloud Physics Lidar) (McGill et al., 2002). There is a total of five AirHARP ocean scenes available in the ACEPOL measurements with retrieval uncertainties studied by Gao et al. (2022) without considering angular uncertainty correlation. Gao et al. (2021a) reported that the retrieval results of both aerosol and ocean color signals have been found to be in good agreement with the AERONET Ocean Color site (Zibordi et al., 2009). However, the cost function histogram was much wider than expected due to the impacts of cirrus clouds. After removing the cirrus cloud impacts from the multiple-angle measurement using an adaptive data screening method, the cost function histogram improved significantly, with much higher similarity with a χ² distribution (Gao et al., 2021b). In this study, the same adaptive data screening methods are applied on all the five AirHARP ocean scenes, which removes cirrus clouds and other anomalies that could not be represented adequately by the current forward model. The resulting total number of measurements, including both reflectance and DoLP, are shown in Fig. 13, where the spatial distribution of pixels with many valid viewing angles is not uniform.

Figure 13 The number of total viewing angles (N) considering both reflectance and DoLP for the five AirHARP ocean scenes from ACEPOL.

Correlation properties from real measurements are difficult to quantity. However, we show in Sect. 2.5 that, as also demonstrated by Sect. 4.2, the autocorrelation analysis on the retrieval residual can be used as a good estimator when correlation is not too strong with consistent behaviors to derive the correlation parameters. Therefore, using the ACEPOL data, we can estimate angular uncertainty correlation from retrieval residuals. From each scene we selected 200 pixels with the most available number of angles, which are clustered together around the region with the maximum number of measurements as shown in Fig. 14. The retrieval residual data were then normalized as discussed in Sect. 2.5 to remove impacts by the non-uniform mean and variance, which are often observed from the real data residuals. The autocorrelation and partial autocorrelation are calculated to assess whether the AR(1) model is sufficient, with examples shown in Fig. 14, for the fitting residuals of reflectance and DoLP from Scene 3 at the 670 nm band. Partial autocorrelation for reflectance showed similar results for the synthetic data in Fig. 3b, with only the first-order term prominent, which suggests that the AR(1) model is sufficient to describe the fitting residual for reflectance, with higher-order contributions negligible. Figure 3d for DoLP shows that higher-order terms may also contribute to the uncertainty model; however, the overall correlation strength is small. From these plots, we can estimate the correlation parameters r following Eq. (14), where ${\stackrel{\mathrm{̃}}{R}}_{\mathrm{1}}$ corresponds to the first-order point in the autocorrelation plots. The values of r are approximately 0.9 and 0.7 for reflectance and DoLP, respectively, corresponding to correlation angles of approximately 10 and 15^∘ following Eq. (5).

Figure 14 Autocorrelation (a, c) and partial autocorrelation (b, d) for the fitting residuals of reflectance and DoLP from Scene 3 at 670 nm bands.

We analyzed the fitting residuals for all the five AirHARP scenes at the four bands with results summarized in Fig. 15. The minimal number of available angles used in the correlation estimation at each band are shown in Fig. 15a. Note that different scenes have various numbers of angles removed due to the impacts of cirrus cloud or other anomalies. Those values are filled with zeros, which may reduce the strength of angular correlation. Therefore, the real correlation angle is likely larger than these values because only a subset of the total measurement data are used to estimate these angles. The estimated correlation angles for reflectance and DoLP vary mostly between 5 and 20^∘. The correlation angles are smaller for the 670 nm band, probably due to overfitting that partially removes the correlated errors as real retrieval signals, consistent with our observations based on synthetic data (Fig. 12). The correlation angles for the 440 and 550 nm bands are largest, which may be mostly close to the truth. The correlation angles are generally larger for reflectance, with a value between 10 and 20^∘ compared to DoLP with a value between 5 and 15^∘, which suggests different amounts of correlation in the reflectance and DoLP data. Comparing the retrieval uncertainties for θ_c around 10 to 20^∘ as shown in Fig. 6, the impact of the correlation to the real retrieval is small, but there is potential for large reduction of theoretical uncertainties by 25 % to 50 % if the correct correlation is considered in the retrieval cost function. This requires, however, that the retrieval algorithm is capable of achieving its optimal performance as described by its theoretical uncertainties. Furthermore, the information on the correlation properties is useful to parameterize realistic measurement uncertainties into synthetic data. The correlation angles for reflectance and DoLP are likely larger than 10 and 5^∘, respectively, which correspond to an estimated correlation parameter of 0.9 and 0.8.

Figure 15 (a) The minimal number of angles used to estimate correlation parameters from each scene and wavelength; (b, c) the estimated correlation angles.

5 Discussions and conclusions

In this study, we evaluate the impacts of angular correlation on the retrieval uncertainties for various aerosol microphysical and optical properties, ocean surface properties, and water-leaving signals. Theoretical uncertainties are derived based on error propagation, and the real uncertainties are obtained through the comparison of retrieved and true values. The theoretical and real uncertainties are compared and discussed. Only small angular correlation impacts are found on the real retrieval uncertainties unless the correlation strength is large (such as with a correlation angle larger than ${\mathit{\theta }}_{\mathrm{c}}>\mathrm{10}{}^{\circ }$). The impacts vary with different retrieval parameters. Theoretical uncertainties are more impacted by angular correlation, which suggests the retrieval algorithm may not always converge to the global minima and that there is potential room for algorithm improvement.

Studies on the fitting residuals from both synthetic data and real AirHARP measurements were conducted. Autocorrelation is useful to estimate the angular correlation, though it tends to underestimate when the correlation strength is strong, and thus overfitting of the measurements is likely. Analysis on the real data showed that the angular correlation is stronger in the reflectance data than DoLP, which makes sense because we expect that DoLP is less sensitive to systematic uncertainties that are more likely to be correlated. Partial autocorrelation analysis suggests that the uncertainty model considering a linear Markov process (AR(1)) is sufficient for reflectance but may need to be further studied for DoLP. From AirHARP retrieval residual analysis, the correlation angles for reflectance and DoLP are estimated to be larger than 10 and 5^∘, corresponding to correlation parameters larger than 0.9 and 0.8, respectively.

This work intends to provide basic methodology to analyze the measurement uncertainties with angular correlations, but the methods can also be applied in the spatial and spectral domains that may be more appropriate for other instruments. There are several remaining issues that need discussion in future works.

• Application to real data. It is complex to analyze real data as we discussed in previous sections. The major challenges and possible issues that may impact uncertainty correlation estimation are summarized below.

  • 1.

    The retrieval is based on a forward model which also has uncertainties, a portion of which may be correlated. This uncertainty will contribute to the fitting residuals and may impact correlation analysis, but it is difficult to quantify.

  • 2.

    The fitting residuals are often not stationary with uniform mean and variance. To reduce this issue, the residuals are normalized, but it would be valuable to analyze how the mean value and variance depend on the angle, as this may provide insight into the modeling uncertainties.

  • 3.

    Some residuals are not continuous with angle due to removed cirrus clouds, which may reduce the correlation.

  • 4.

    Synthetic data analysis has demonstrated that the retrieval is likely to overfit the data when the correlation is strong.

  • 5.

    The angular grids for HARP measurements are slightly non-uniform, which is likely to further reduce the correlation strength from autocorrelation analysis. To evaluate impacts of this feature, an uncertainty model considering the impact of the real angular grids needs to be built. But the variation of the angular grids is less than 1^∘ (670 nm band) or 2^∘ (other bands), which may impact more the cases with small correlation angles.

• Lab calibration. Although the correlation strength is estimated from fitting residuals from real AirHARP measurements, it may be only used for qualitative discussions due to various issues discussed above. To obtain the actual correlation properties, lab characterizations are desired to separate measurement characteristics. Lab measurement signals may also be evaluated through autocorrelation and partial correlation functions, or more general ARMA models. It is also interesting to discuss the possible impacts in the uncertainty model due to the binning and collocation that happen in later processing steps.

• Correlation strength as a fitting metric. Due to the limitation discussed above, our analysis on the fitting residuals may only provide the lower boundary for the correlation strength in terms of correlation parameter (r) or correlation angle (θ_c). Since the correlation strength represents properties of the fitting residuals, it can be also used as a metric to represent retrieval fitting performance, together with the cost function χ² and variance of the fitting residuals as discussed in Sect. 4.1.

• Signal correlation vs. uncertainty correlation. Nature as measured by the instrument and expressed in the forward model has inherent correlation, which becomes part of the retrieval process. The phenomena we observe tend to be only slowly changing with respect to view angle, and thus measurements at different view angles do not necessarily express retrieved parameters independently. This correlation is related to the actual signal in the measurement rather than its uncertainties. This type of correlation is captured by the Jacobian matrix. The overall information content of the measurements with respect to the set of retrieval parameters is determined by both the Jacobian matrix and the correlated uncertainty model as further discussed in Appendix B.

• Future retrieval algorithm development. The pixel-wise theoretical uncertainties achieve a reasonably good performance to represent real retrievals when no correlation is presented. Their performances on various retrieved geophysical properties are quantified by comparing with the real retrieval errors. The difference grows bigger when the angular correlation is stronger, which suggests convergence to local minima and indicates that more development is needed to improve the retrieval optimization. For the development of future algorithms with more retrieval parameters, such as aerosols with more complex shape and absorption properties and coastal waters with more complex bio-optical properties, a better characterized error model, such as the one considering angular or spectral correlations, will be helpful to identify information useful for the retrievals and therefore improve retrieval performance and uncertainty assessment.

Appendix A: Input parameters of the neural-network forward model

A total of 15 parameters are used as input of the forward model as discussed in Sect. 2.1 as listed in Table A1. The solar zenith (θ₀), viewing zenith angle (θ_v), viewing azimuth angle relative to solar direction (ϕ_v), and ozone density ($n_{{\mathrm{O}}_{\mathrm{3}}}$) are assumed to be known input. All the other 11 parameters are retrieval parameters in the FastMAPOL algorithm, including the aerosol volume density for each sub-mode (V_i), the real (m_r) and imaginary parts (m_i) of the refractive index for fine- and coarse-mode aerosols, ocean surface wind speed (w), and chlorophyll a concentration (Chl a).

Table A1 Parameters used to train the FastMAPOL forward model. The minimum (min) and maximum (max) values of each parameter are also shown. The a priori uncertainties (σ_a) are estimated as the difference between the maximum and minimum values for the study, except the four parameters as indicated, which are assumed to be known input.

Appendix B: Eigenvector decomposition in the error covariance matrix

Section 2.4 discussed that the error covariance matrix with non-diagonal terms can be diagonalized through eigen-decomposition as shown in Eq. (7). The original measurements can be transformed into a new space without correlations, with uncertainty variance described by the eigenvalues in the diagonal matrix D_ϵ denoted as $d_i^{\mathrm{2}}$. To understand how the uncertainties vary with different correlation strength, the square root of the diagonal term in D_ϵ is used to represent the new measurement uncertainties. Results for different correlation angles are shown in Fig. B1.

Figure B1 The uncertainties for reflectance and DoLP of the error covariance matrix at the 670 nm band after eigenvector decomposition as discussed in Sect. 2.4. These values correspond to the square root of the diagonal values in matrix D_ϵ from Eq. (7).

The uncertainties in the new measurement space (y^′) show components both larger and smaller than the original uncertainties from different combinations of measurements. The corresponding Shannon information content (SIC) (Rodgers, 2000) is defined as

$$\begin{array}{}\text{(B1)}& {\displaystyle }\text{SIC}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\ln \left|{\mathbf{S}}^{-\mathrm{1}}{\mathbf{S}}_{\mathbf{a}}\right|,\text{(B2)}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\ln |{{\mathbf{S}}_{\mathbf{a}}}^{\mathrm{1}/\mathrm{2}}{{\mathbf{K}}^{\prime }}^T{\mathbf{D}}_{\mathit{ϵ}}^{-\mathrm{1}}{\mathbf{K}}^{\prime }{{\mathbf{S}}_{\mathbf{a}}}^{\mathrm{1}/\mathrm{2}}+\mathbf{I}|,\end{array}$$

where the error covariance matrix S and the a priori matrix S_a in Fig. B2 are from Eq. (6). The Jacobian matrix and the error covariance matrix are converted into the diagonal space as shown in Sect. 2.4 and are used to represent the SIC in Eq. (B1). Different correlation strength will lead to a different unitary matrix, which is used to transform the Jacobian and error covariance matrix (Eq. 7). We analyzed the SIC and retrieval uncertainties for a sequence of correlation angles from 0 to 120^∘. As shown in Fig. B2, when the correlation strength is strong, the SIC is increasing with respect to correlation, and the theoretical retrieval uncertainties estimated from error propagation are decreasing, such as for fine-mode refractive index, wind speed, and Chl a. However, when the correlation is relatively weak, SIC decreases and uncertainties increase with the correlations. The correlation angles with maximum uncertainties are different with different retrieval parameters but generally fall within 30^∘. This behavior may relate how the measurement uncertainties are mapped to the retrieval parameters space through the Jacobian matrix. Similar behavior of the SIC has been reported by Knobelspiesse et al. (2012) for the RSP measurement. Understanding of the uncertainty properties such as correlation strength is useful to further exploit the information in the measurements.

Figure B2 The Shannon information content (SIC) and corresponding theoretical retrieval uncertainties for refractive index (m_r,f), wind speed (w) and Chl a with respect to various correlation angles. The results are for one example case with Chl a =0.1 mg m⁻³, w=3.0 m s⁻¹, $m_{\mathrm{r},\mathrm{f}}=\mathrm{1.55}$ and AOD (550 nm) =0.2 – the same as the case discussed in Gao et al. (2022) (ID: 201).