Adaptive Data Screening for Multi-Angle Polarimetric Aerosol and Ocean Color Remote Sensing Accelerated by Deep Learning

Remote sensing measurements from multi-angle polarimeters (MAPs) contain rich aerosol microphysical property information, and these sensors have been used to perform retrievals in optically complex atmosphere and ocean systems. Previous studies have concluded that, generally, five moderately separated viewing angles in each spectral band provide sufficient accuracy for aerosol property retrievals, with performance gradually saturating as angles are added above that threshold. The Hyper-Angular Rainbow Polarimeter (HARP) instruments provide high angular sampling with a total of 90–120 unique angles across four bands, a capability developed mainly for liquid cloud retrievals. In practice, not all view angles are optimal for aerosol retrievals due to impacts of clouds, sunglint, and other impediments. The many viewing angles of HARP can provide resilience to these effects, if the impacted views are screened from the dataset, as the remaining views may be sufficient for successful analysis. In this study, we discuss how the number of available viewing angles impacts aerosol and ocean color retrieval uncertainties, as applied to two versions of the HARP instrument. AirHARP is an airborne prototype that was deployed in the ACEPOL field campaign, while HARP2 is an instrument in development for the upcoming NASA Plankton, Aerosol, Cloud, ocean Ecosystem (PACE) mission. Based on synthetic data, we find that a total of 20–30 angles across all bands (i.e., five to eight viewing angles per band) are sufficient to achieve good retrieval performance. Following from this result, we develop an adaptive multi-angle polarimetric data screening (MAPDS) approach to evaluate data quality by comparing measurements with their best-fitted forward model. The FastMAPOL retrieval algorithm is used to retrieve scene geophysical values, by matching an efficient, deep learning-based, radiative transfer emulator to observations. The data screening method effectively identifies and removes viewing angles affected by thin cirrus clouds and other anomalies, improving retrieval performance. This was tested with AirHARP data, and we found agreement with the High Spectral Resolution Lidar-2 (HSRL-2) aerosol data. The data screening approach can be applied to modern satellite remote sensing missions, such as PACE, where a large amount of multi-angle, hyperspectral, polarimetric measurements will be collected.


INTRODUCTION
Aerosols play a critical role in Earth's radiative balance by directly scattering and absorbing solar radiation, and indirectly interacting with clouds. Due to their complex micro-physical properties and spatio-temporal distributions, combined with measurement and modeling difficulties, aerosol-related processes are among the largest uncertainties in radiative forcing of the Earth climate (Boucher et al., 2013). Remote sensing measurements from multi-angle polarimeters (MAPs) contain rich information on aerosol microphysical properties compared to other sensor types and therefore can be used to improve aerosol retrieval accuracy (Mishchenko and Travis, 1997;Chowdhary et al., 2001;Hasekamp and Landgraf, 2007;Knobelspiesse et al., 2012). Meanwhile, ocean color remote sensing, which assesses water-leaving signals to infer chlorophyll concentration and other biogeochemical quantities, is important in the study of phytoplankton dynamics, global carbon cycle, and marine ecosystems (Frouin et al., 2019;Groom et al., 2019). Accurate estimation of the water-leaving signal requires the quantification and removal of the aerosol path radiance and the ocean surface reflectance from the remote sensing measurement (Mobley et al., 2016). To advance both aerosol and ocean color characterization based on MAP measurements, simultaneous multi-parameter retrieval algorithms have been developed over both open and coastal waters (Chowdhary et al., 2005;Hasekamp et al., 2011;Xu et al., 2016;Stamnes et al., 2018;Gao et al., 2018;Fan et al., 2019;Gao et al., 2019Gao et al., , 2021. To achieve various scientific objectives on the study of aerosols and clouds, several MAP instrument designs are available, with different spectral bands, viewing angles, and measurement accuracies (see Table 1). The Polarization and Directionality of the Earth's Reflectances (POLDER) instruments flew on Advanced Earth Observing Satellite missions (ADEOS-I;1996-1997 andADEOS-II;2002-2003, Deschamps et al. (1994) and the Polarization and Anisotropy of Reflectances for Atmospheric Sciences coupled with Observations from a Lidar (PARASOL;2004-2013 mission (Tanré et al., 2011). They measured up to 14 viewing angels in nine spectral bands from visible to near infrared (NIR) with three polarized bands. The 3MI mission, planned to launch in 2023, improves the POLDER design with a total of 12 spectral bands which include nine polarized bands and an extension to shortwave infrared (SWIR) (Fougnie et al., 2018). Similar to POLDER with three polarized bands, the Airborne Multiangle SpectroPolarimetric Imager (AirMSPI, Diner et al. (2018)) and its spaceborne version MAIA (planned to launch in 2022, Diner et al. (2018)) conduct measurements with eight spectral bands from visible to NIR and 14 bands from visible to SWIR, respectively, both with a selectable number of viewing angles typically ranging from five to nine. The airborne Research Scanning Polarimeter (RSP) measures a much higher angular resolution with 152 viewing angles in nine bands from visible to SWIR (Cairns et al., 1999); its space-borne analogue the Aerosol Polarimetry Sensor (APS) on NASA Glory mission (failed to reach orbit in 2011) would have had a total of 250 viewing angles (Mishchenko et al., 2007). The Airborne Hyper-Angular Rainbow Polarimeter (AirHARP) and its space-borne counterpart the HARP CubeSat (launched in 2020), and future HARP2 on NASA's Plankton, Aerosol, Cloud, ocean Ecosystem mission (planned for launch in 2024), also makes high angular resolution measurements (Martins et al., 2018;McBride et al., 2020). Unlike the other MAP designs, the viewing angles of the HARP polarimeters are slightly different for each band, with a total of 120 unique angles across the four bands for AirHARP and HARP CubeSat (60 angles for 670 nm band; 20 angles for each of the other three bands), and a total of 90 unique angles for HARP2 (60 angles for 670 nm band; 10 angles for each of the other three bands). The Spectro-Polarimeter for Planetary EXploration one (SPEXone) is the other MAP planned for PACE (Hasekamp et al., 2019). SPEXone and its airborne version (SPEX airborne, Smit et al. (2019)) measure at five and nine viewing angles in each band, respectively, with continuous spectral sampling (400 spectral bands, 50 polarized) mostly in visible bands. A thorough review of the MAP instruments and their corresponding retrieval algorithms can be found in Dubovik et al. (2019). Generally, SPEX, AirMSPI and MAIA, with a small number of viewing angles (5-9), are optimized for aerosol studies; POLDER and 3MI with up to 16 viewing angles are designed for monitoring aerosol and to some extent clouds; RSP, APS, and HARP instruments (AirHARP, HARP CubeSat and HARP2) with a much higher number of viewing angles (90-250) are capable of resolving cloud bow angular patterns and are optimized for cloud studies (Waquet et al., 2009;McBride et al., 2020) in addition to aerosols.
Extensive research has studied how aerosol retrieval uncertainties depend on the viewing angles available from MAP measurements. Hasekamp and Landgraf (2007) discussed the trade-off between spectral sampling and the number of viewing angles using synthetic data with spectral bands similar to POLDER, and found that the retrieval errors of aerosol optical depth (AOD), refractive index, and size are only slightly affected when the number of viewing angles is increased at the cost of the number of spectral sampling points, as long as at least three viewing angles are available. Based on synthetic high accuracy RSP measurements with either SWIR bands included or not, Wu et al. (2015) found that five angles equally spaced between ±60°provide sufficient accuracy for aerosol property retrievals, with performance gradually saturating as more angles are added beyond that. Similarly, Xu et al. (2017) concluded from AirMSPI measurements that the retrieval accuracy of AOD and single scattering albedo (SSA) do not improve significantly when more than around five viewing angles are used.
Although five viewing angles are sufficient for aerosol studies, in practice there will be ground pixels for which not all angles are optimal for aerosol retrievals due to impacts by cloud, sunglint, calibration artifacts, etc. Gao et al. (2019) used RSP measurements from the Ship-Aircraft Bio-Optical Research (SABOR) field campaign and found that viewing angles impacted by cloud needed to be removed in order to obtain high-quality retrieval results. Gao et al. (2021) found that AirHARP pixels impacted by cirrus clouds during the Aerosol Characterization from Polarimeter and Lidar (ACEPOL) field campaign produced large retrieval cost function values, indicative of larger deviation between the forward radiative transfer model and the measurements. Further, strong sunglint signals measured by both RSP and AirHARP instruments could not be modeled sufficiently well by the isotropic Cox and Munk (1954) model (Gao et al., 2019. Therefore, it is likely the number of suitable viewing angles for a given scene will be lower than those acquired by the sensor. The abundant angular measurements from the HARP instruments are useful to explore data quality screening by removing potentially problematic angles while maintaining sufficient useful angles for retrievals. Due to their frequent presence, clouds are one of the major factors reducing the number of suitable viewing angles from multi-angular measurements (whether polarized or not) for aerosol studies. The MAP often observes a maximum angular range of 110-120° (Table 1). Depending on the size and height of a cloud, it may obscure Earth from all view angles or only over a small angular range. To address the challenges in the aerosol retrieval with cloud present, three approaches are often taken: 1. Completely remove the pixels influenced by clouds, using a cloud mask before (e.g. Garay et al. (2020)) and/or additional filtering after (Stap et al., 2015;Chen et al., 2020) the retrieval. 2. Simultaneous cloud and aerosol retrievals. This assumes an optically thin cloud covers most of the region in the observation, or a mixing of cloudy and clear sky pixels (Hasekamp, 2010;Stap et al., 2016). 3. For partially cloud-covered pixels, the angles impacted by cloud can be removed, and then a retrieval assuming only the presence of aerosols can be conducted (Gao et al., 2019).
For a pixel with most or all of its angles covered by cloud the retrieval can be either discarded (approach 1), or have both cloud and aerosol properties included in the state vector (approach 2). However, for a partially cloudy scene, there are open questions as to whether cloud-influenced angles can be removed effectively, and how much removal is practical before retrieval performance is degraded. In this study we focus on approach 3) and aim to answer the following three questions: 1. How many angles are sufficient in order to retrieve aerosol properties and ocean water leaving reflectance from HARP instruments with sufficient accuracy? 2. How can view angles problematic for aerosol retrievals, such as those influenced by clouds, be efficiently identified and removed (screened)? 3. How does the data screening improve the performance of aerosol and ocean color retrievals?
To understand to what extent the aerosol retrievals are still valid, we will examine the impact of the number of viewing angles on retrieval uncertainties using synthetic data simulated for both AirHARP and HARP2, which have different polarimetric uncertainties. The impacts of the scattering angle range will also be considered. To identify the pixels impacted by water clouds, Stap et al. (2015) proposed using the retrieval cost function, a goodness-offit statistic based on how well the retrieval solution matches the observed measurements. Pixels whose cost function exceeded a threshold were deemed cloud-contaminated, and the whole pixel then discarded (approach 1 above). Building on this, we evaluate the fitting residuals (difference between the forward model and measurement) at each view angle. We implement a multi-angle polarimetric data screening (MAPDS) approach using these residuals as a metric to evaluate data quality, built upon the Fast Multi-Angular Polarimetric Ocean coLor (FastMAPOL) retrieval algorithm . FastMAPOL provides an efficient way to conduct aerosol and ocean color retrievals using efficient neural network (NN) forward models, while the data screening approach can automatically and adaptively identify and remove observations at view angles that cannot be wellrepresented by the forward model. Using AirHARP data, we will demonstrate the effectiveness of the data screening and how retrieval uncertainty changes after a portion of angles are removed.
MAPDS requires extra retrieval iterations and therefore slows down the retrieval. Retrieval speed can be improved using the neural networks for MAP measurements in both inverse model (Di Noia et al., 2015;Di Noia and Hasekamp, 2018), as well as forward models for ocean reflectance (Fan et al., 2019;Mukherjee et al., 2020) and coupled atmosphere and ocean model models . Another way to compensate for the increased computational cost, and also improve the overall FastMAPOL retrieval speed, is to improve the calculation of the Jacobian matrix that guides retrieval to a solution. This is challenging in MAP inversions due to the large number of retrieval parameters. To expedite the calculation, linearized radiative transfer models based on forward-adjoint perturbation theory have been developed (Spurr, 2008;Schepers et al., 2014). However, the linearization of the radiative transfer code requires sophisticated algorithm development. Instead, we employ automatic differentiation (AD) based on the NN forward model that acts as our radiative transfer (RT) emulator . AD Jacobians are based on the chain rule of derivatives (Baydin et al., 2018) and successfully reduce computational expense.
NASA's PACE mission includes a hyperspectral scanning radiometer named the Ocean Color Instrument (OCI) and the two above-mentioned MAPs: HARP2 and SPEXone . Together they will collect a large amount of multiangle, hyperspectral, polarimetric measurements essential for the characterization of atmosphere and ocean states (Remer et al., 2019a,b;Frouin et al., 2019). To facilitate cross-calibration and data synergy, the measurements from all three PACE instruments will be projected onto a common PACE Level-1C data format with a uniform spatial grid. Aerosol and ocean color retrievals based on multiple instruments, with larger data volume, require even higher computational efficiency. The combination of FastMAPOL and data screening provide an efficient approach to analyze and process the large volume of measurements from PACE, either from one instrument or in combinations with others.
The impacts of ocean color retrievals from data screening will be also discussed in this study. To derive ocean color signals from the MAP measurement, we conducted atmospheric correction by subtracting the contributions of the atmosphere and ocean surface from the spaceborne or airborne measurements taken at the top of atmosphere (TOA) (Mobley et al., 2016). Note that for nonpolarimetric ocean color sensors, neural networks trained on simulated TOA reflectance have been developed to conduct atmospheric correction. For example, the case 2 Regional Coast Colour (C2RCC) algorithm has been applied on ESA's MEdium Resolution Imaging Spectrometer (MERIS) measurements (Doerffer and Schiller, 2007), and the Ocean Color -Simultaneous Marine and Aerosol Retrieval Tool (OC-SMART) algorithm has been applied on NASA's Moderate Resolution Imaging Spectroradiometer (MODIS) measurements (Fan et al., 2021). In this study, neural networks are trained to represent the radiative transfer forward model over coupled atmosphere and ocean system which is used in the inversions of both aerosol and ocean parameters . The retrieved aerosol microphysical properties are then used to conduct the atmospheric correction on the multi-angle measurements. This procedure also allows the application of the MAP retrieved aerosol properties to atmospheric correction of co-located hyperspectral measurements from other sensors, such as PACE OCI Hannadige et al., 2021). This paper is organized as follow: Section 2 provides the methodology including the retrieval algorithm, NN forward model and AD Jacobian matrix; Section 3 describes synthetic HARP2 and AirHARP data retrievals to address question 1); Section 4 covers AirHARP filed data retrievals which address question 2) and 3); Section 5 discusses cloud mask and uncertainty analysis; Section 6 is the conclusions.

Aerosol and Ocean Color Retrieval Algorithm
This study is based upon the efficient aerosol and ocean color retrieval algorithm FastMAPOL , which uses a NN forward model to conduct vector radiative transfer calculations for a coupled atmosphere and ocean system (Zhai et al., 2009(Zhai et al., , 2010). The forward model assumes a three-layer atmosphere: the bottom layer from ocean surface to an altitude of 2 km containing aerosols and molecules, a molecular layer between 2 km and the aircraft (at 20 km), and an additional molecular layer above the aircraft.
A total of 15 parameters are used in the forward radiative transfer simulation, as listed in Table 4 (more details in Gao et al. (2021)). The aerosol optical properties are determined by the fine and coarse mode complex refractive indices (total 4 parameters), and the volume densities (μm 3 μm −2 ) for five aerosol size submodes with each submode following a lognormal distribution. The mean radii and variances of the submodes are prescribed and fixed in the retrieval. The three smaller submodes in combination represent the fine mode, and the two larger submodes in combination represent the coarse mode. In this study we only consider open ocean waters modeled as a uniform layer with its bio-optical properties parameterized in terms of the chlorophyll-a concentration (Chla, mg · m −3 (Gao et al., 2019). Complex water properties over coastal water require additional parameters in the biooptical model (Gao et al., 2018) which is subject to future investigations. The ocean surface roughness is represented by the isotropic model parameterized by wind speed (ms −1 ) (Cox and Munk, 1954). As mentioned in Section 1, the sunglint signal measured by AirHARP can not be modeled well by the isotropic Cox and Munk model. In addition, it is challenging to accurately represent sunglint in the NN training due to its large amplitude in comparison with non-glint angles . To minimize the impacts of the sunglint, we limited the wind speed within 10 ms −1 and excluded the view angles within 40°from the direction of specular reflection for the analysis in the next section. Larger wind speed and the associated white cap will be considered in future studies. Therefore, total 11 parameters are retrieved by FastMAPOL; the other four parameters needed by the forward model (viewing and solar zenith angle, relative azimuth angle, and ozone profile) are assumed to be known.
FastMAPOL determines optimal values of state parameters by minimizing the difference between the measurements and the NN forward model prediction. An iterative optimization approach is used as summarized in Figure 1A. The NN forward model computes the reflectance and degree of linear polarization (DoLP), defined as.
where L t , Q t , and U t are the Stokes parameters at the sensor altitude, F 0 is the extraterrestrial solar irradiance, μ 0 is the cosine of the solar zenith angle. The differences between measurements and model predictions are represented by a cost function χ 2 based on Bayesian theory (Rodgers, 2000): where ρ f t and P f t are the reflectance and DoLP computed from the NN forward model. The state vector x contains the forward-modeled parameters described above; the subscript i stands for the index of the measurements (where one "measurement" is defined as a pair of reflectance and DoLP) at different viewing angles and wavelengths; and N is the total number of measurements used in the retrieval. The maximum possible N assuming all angles are available is summarized for various MAPs in Table 1, which is 240 for AirHARP and HARP CubeSat, and 180 for HARP2.
The total uncertainties of the reflectance and DoLP used in the algorithm are denoted by σ ρ and σ P , respectively. These include contributions from the instrument uncertainties σ ins , the NN i − χ 2 i−1 | indicates the changes of the cost function between two iterations with η as threshold (Eq. 6). In panel (B), Δρ t and ΔP t indicate the difference between forward model and measurements for reflectance and DoLP with ξ as threshold (Eq. 8). The dashed box in (B) represents the same retrieval process as shown in the dashed box in (A). A maximum three passes (indicated by the loop in (B)) are used in the data screening process.
TABLE 2 | Uncertainties for reflectance (ρ) and DoLP (P) for both measurement and forward model used in this study including instrument uncertainty (σ ins ), the radiative transfer simulation uncertainty (σ RT ), and the NN uncertainty (σ NN ). Different σ P,ins is used for AirHARP and HARP2 but with the same σ ρ,ins .
Uncertainties 440 nm 550 nm 670 nm 870 nm forward model uncertainties σ NN , and the radiative transfer simulation uncertainties σ RT used to train the NN: Table 2 summarizes the uncertainties used in this study; the values of σ ρ,NN and σ ρ,RT were determined in Gao et al. (2021) with the same NN models as used in this study. Note that, at present, uncertainties are assumed spectrally and angularly uncorrelated.
The FastMAPOL retrieval follows the flow chart in Figure 1A. In each iteration, if χ 2 is larger than a threshold, the state parameter x is updated to compute both ρ f t and P f t , and a new value of cost will be computed. The Jacobian matrix (the gradient of the cost function with respect to x) is used to determine the magnitude and direction for the updated estimate of x. This process iterates until the convergence criterion is met: where i is the iteration index in optimization, and ϵ is the convergence tolerance, taken as 0.01 in this study. The convergence of the retrievals can be examined by the cost function histogram over an ensemble of retrievals (more discussion in the next section). If the retrievals converge well, the normalized cost function histogram can be represented by the theoretical probability density function (PDF) of the χ 2 distribution: where χ 2 is the cost function (Eq. 3), k is the degree of freedom (DOF), Γ(k/2) denotes the gamma function (James, 2006). The PDF is useful in understanding the behavior of the retrieval cost function distributions with different numbers of measurements used. After neglecting potential correlation between measurement uncertainties, we can approximate the DOF by the total number of measurements N, which shows good representation as demonstrated in the next section. Note that the DOF refers to the retrieval residuals in the cost function (Eq. 3), which are mostly contributed by the noise and uncertainties in the measurements and forward model. The DOF is likely between N and N-11, where 11 is the total number of retrieval parameters. The actual number of DOF can be determined from the Jacobian matrix and error covariance matrix (Rodgers, 2000).

Adaptive Data Screening
Based on the FastMAPOL retrieval algorithm, the MAPDS data screening approach, which conducts automatic data quality analysis and data screening, is summarized in the flowchart of Figure 1B. After each converged FastMAPOL retrieval, the residuals between each measurement and the forward model prediction are used to evaluate the data quality under the criteria: where the residuals are compared with the uncertainty model defined in the cost function of Eq. 3, and ξ is a threshold. When either the reflectance or DoLP does not satisfy the criteria, the corresponding measurement is excluded from the cost function calculation (i.e., the view angle which cannot be represented well by the forward model is removed), and a new FastMAPOL retrieval is performed. In practice, additional screening rules based on the criteria may be used depending on the data quality of the field measurements (illustrated with an example later in Section 4). This process is repeated until all angles remaining satisfy Eq. 8. Note that the whole data screening process will include several retrieval passes, each involving multiple iterations until convergence. At the end of each retrieval pass based on the new forward model fittings, all measurements used in the retrieval are evaluated through Eq. 8 and subsequently excluded from the next round if they failed to pass. The data screening approach is an adaptive process, since it depends on the fitting at each iteration for each pixel. A threshold value of ξ 3 is used in this study, which can be further tuned when more data is available. We found at most three passes of retrievals are sufficient to remove most of the problematic angles. Since the retrieved parameters from last retrieval can be used as the initial values to the next retrieval as shown in Figure 1, the total speed to conduct data screening are usually less than three times of the single round retrieval.

The NN Forward Model and Automatic Differentiation
The MAP retrievals are often computationally expensive due to their high dimensionality and iterative nature, with multiple forward model and Jacobian calculations. The data screening 3 | NN forward model and its Jacobian matrix computed by the forward mode and reverse mode of automatic differentiation (AD). The arrows indicate the order of steps for each process.
*For brevity, the summation of index l is implied following the Einstein notation.
Frontiers in Remote Sensing | www.frontiersin.org December 2021 | Volume 2 | Article 757832 approach developed here further increases the demand for CPU computations because the retrieval must be repeated several times. Therefore, fast forward model and Jacobian matrix calculations are advantageous for efficient processing, which was a motivation for the use of NN forward models for AirHARP in . In this work we further improve the efficiency through automatic differentiation to compute the Jacobian matrix analytically (as opposed to numerically through finite differencing), exploiting the differentiable properties of the NN models. The NN forward model developed in Gao et al. (2021) is a feed-forward neural network as defined in Table 3, where h 0 x is the input layer that contains all 15 forward model parameters. Two sets of weight matrix W p+1 and bias vector b p+1 have been determined from the NN training process for reflectance and DoLP respectively , both with three hidden layers (k 3) of 1,024, 256 and 128 nodes. Correspondingly, y is the output layer for either reflectance or DoLP at the four AirHARP bands.
For application to multi-angle measurements, the NN needs to be called to simulate y for each set of viewing and solar geometries for the state vector x. Elements of the Jacobian matrix are defined as follows: Here index m indicates the viewing and solar angles, index i indicates the wavelength, and index j indicates the state parameter. Table 3 illustrates the structure of the NN and how the Jacobian matrix is calculated based on the trained weights. In that table, the nonlinear activation function Φ is the LeakyReLU function, which is defined as where α 0.01, and Z is a matrix with m and i as its indices. The derivative with respect to each element in Φ is defined: The finite difference (FD) method was used to compute the Jacobian matrix in FastMAPOL in , where the NN forward model was called twice per input parameter under the central difference approximation of derivatives. To reduce the computational cost, the Jacobian matrix can be derived analytically from the NN forward model using AD based on the chain rule of differentiation (Baydin et al., 2018). Two recursive relations are obtained to compute the Jacobian matrix as shown in Table 3, where the forward mode indicates the evaluation sequence from the first layer to the last layer, and the reverse mode indicates the evaluation sequence from the last layer to the first layer. To represent the recursive relations, we define _ h p (tangent) and h p (adjoint) as follows: Note that h is defined in Table 3 (left column) as the output from each hidden layer of the NN. The Jacobian matrix can be represented by AD with either the tangent or the adjoint forms as the last step in Table 3: where Eqs. 14, 15 are computed from the forward and reverse mode AD, respectively. Forward and reverse AD produce identical results, but differ in computational efficiency due to the different sequence of matrix operations and NN architecture. For optimal efficiency, we implemented AD directly based on the formalism summarized in Table 3 using the Pytorch library (Paszke et al., 2019). The arrows in the table indicate the order of calculation steps. NN forward model can be computed layer by layer, with the output from the previous layer as the input to the next layer. The forward mode AD can be computed in the same sequence as the NN forward model as shown in Table 3. The reverse mode AD is computed from the last layer backward to the first layer. Note that AD in both modes requires the values of matrix D as defined in Eq. 11, which is determined by the output of the NN forward model at each layer. For the NN used in this study, AD in reverse mode provides the highest efficiency as investigated further in the next section. The AD methods are efficient and accurate in computing the Jacobian matrix, providing a convenient way to accelerate algorithms like FastMAPOL with a large number of retrieval parameters, making them more suitable for practical applications. 4 | Parameters used to represent the radiative transfer of the atmosphere and ocean system. θ 0 and θ v are the solar and viewing zenith angles. ϕ v is the relative viewing azimuth angle. V i denote the volume densities for the five aerosol submodes. m r and m i are the real and imaginary parts of the refractive index. Subscripts f and c refer to fine and coarse mode.

Parameters
Unit Min Max Initial Ozone column density (n O3 ) in the atmosphere, ocean surface wind speed(w), and chlorophyll a concentration (Chla) are also provided. The minimum (min) and maximum (max) values determine the parameter ranges used in the retrievals. The initial values are the ones used in the retrieval algorithm with θ 0 , θ v , ϕ v and n O3 are assumed to be known from inputs.

Synthetic Data
We performed radiative transfer simulations to generate 1,000 sets of synthetic polarized reflectances for a coupled atmosphere and ocean system as discussed in Section 2 and Gao et al. (2021).
Here, we use a fixed solar zenith angle of 50°as this is approximately the solar zenith angle in the AirHARP field campaign data used in Section 4. Input parameters and their range are listed in Table 4; Chla is randomly sampled from a loguniform distribution, AOD at 550 nm is sampled uniformly within [0,0.5], and the fine-mode volume fraction is sampled uniformly within [0, 1]. All the other parameters in Table 4 are sampled uniformly, except aerosol volume densities (which are determined by AOD and fine mode volume fraction). Simulations are performed with a view angle sampling of 1°. Random noise is added to the measurement according to the instrument uncertainties shown in Table 2.
In order to generate realistic viewing geometries we represent the orbit geometry as shown in Figure 2, where the viewing angles (θ v , ϕ v ) can be sampled according to the along-track and crosstrack angles indicated in the figure.
where θ a and θ c are the along-track and cross-track viewing angles. We use a set of predefined θ a values (120 angles for AirHARP, 90 for HARP2) and then randomly sampled θ c within ±47°according to HARP instrument characteristics. The resulting example geometries are shown in Figure 3A. Sunglint with view angles within a conservative angle of 40°relative to the direction of specular reflection are excluded as indicated in Figure 3. The simulated reflectance and DoLP are interpolated into this viewing geometry, each with a random value of θ c . Note that Figure 2 is a simplified PACE orbit geometry with the solar azimuth angle in the along-track direction; the curvature of Earth's surface is neglected.
To evaluate the retrieval uncertainties using different numbers of viewing angles, we reduce the number of available viewing angles by assuming the MAP measurements are blocked randomly by clouds. Two cloud configuration schemes are considered: 1. Cloud blocking the center part of the measurements, leaving available angles on both sides. Cloud location can move randomly from one side to another side as shown in Figure 3B. 2. Clouds blocking both sides of the measurements, leaving the middle part of cloud-free angles available which can also move randomly from one side to another side as shown in Figure 3C.
We generate viewing geometries with the same total number of viewing angles (N v ) for both schemes as shown in Figures 3B,C. We chose 11 different cases with N v of 5, 8,10,15,20,25,30,40,50,60, and 75 for AirHARP, and similarly for HARP2 with the two largest values as 56 and 58 instead. The largest averaged N v value for both AirHARP and HARP2 corresponds to the case shown in Figure 3A with only sunglint removed and no cloud present. For the same number of N v , Scheme 1 will result in a larger average range of scattering angles than Scheme 2. As shown in Table 1, the angular range of HARP2 and AirHARP measurements is 114°, therefore the cloud angular size (as viewed from ground) used to block the synthetic data can be estimated as 114°− N v × 114°/90 for HARP2, and 114°− N v × 114°/120 for AirHARP. Their physical size depends on both their angular size and height.
The retrievals are conducted for both AirHARP and HARP2 with their measurement uncertainties summarized in Table 2. This provides a total of 44,000 retrievals: 1,000 (RT simulation cases) × 2 (cloud schemes) × 2 (HARP2 and AirHARP) × 11 (total variations in angle number). Note that the total number of measurements N equals 10 for the case with N v 5, lower than the total 11 retrieval parameters, and therefore for that case the inverse problem is ill-posed with non-unique solutions.

Retrieval Efficiency Using FD and AD
The performance using FD with central difference is compared with the forward mode and reverse mode of AD as shown in Table 3 in Section 2. Figure 4 compares the retrieval χ 2 histogram and the retrieval time for the FD and AD methods. The results converge to similar χ 2 distributions; the χ 2 histograms change with respect to the available number of viewing angles, but can be well represented by the theoretical χ 2 distributions from Eq.7 with the same degree of freedom as the total number of measurements (N 2N v ).
Retrievals using conventional radiative transfer simulation with FD usually took 1 hour to converge on a CPU (AMD EPYC Processor). Using neural network forward models, the averaged retrieval time decreased to 3 s using the same FD method and CPU. With forward and reverse modes of AD, a further decrease to 0.6 and 0.3 s respectively was achieved, which is an increase of speed by a factor of 5-10. With GPU processing (GeForce GTX 1060), the retrieval times were 0.08 and 0.05 s, another factor of 6 faster than CPU. AD in reverse mode provides the highest efficiency and will be used in the following discussions as the default method.

Aerosol Retrievals
In this section, we compare the retrieval error as a function of the number of available viewing angles. The retrieval error (except Chla), shown in Figure 5, is defined as the root mean square error (RMSE) between the retrieval results and the simulated truth. The retrieval error of Chla is represented by the Mean Averaged Error (MAE) in log scale which is a better metric for Chla as recommended by Seegers et al. (2018), and indicates the averaged ratio between the retrieval and truth  values. The definition and performance of the retrieval error under different Chla range can be found in Gao et al. (2021). While it is well known that retrieval uncertainties depend on the aerosol loading (cf. Gao et al. (2021) for AirHARP data), for the convenience of discussion here, the retrieval errors are averaged over all cases with AOD within [0.01, 0.5]. The cloud blocking generally produces smaller errors than which has a smaller range of scattering angles. However, coarse mode v eff seems less sensitive to the number of viewing angles and does not have a clear advantage between S1 and S2. Errors for HARP2 are also smaller than AirHARP due to the lower expected DoLP calibration uncertainty in HARP2. Generally, the errors decrease rapidly as the total number of angles increases to 20 (corresponding to 2 × 20 unique measurements including both reflectance and DoLP across four HARP bands). For HARP2, improvement continues for AOD, SSA, m r , r eff , v eff up to 40-50 angles; Chla performance seems to plateau around 20 angles. AOD error decreases the most, by a factor of three from 0.06 to 0.02. The total SSA error also decreases by a factor of two (from 0.06 to 0.03), but for modal (either fine or coarse mode) SSA and m r , r eff and v eff , it only reduces by about 1/3 of its maximum error. Due to the removal of sunglint, the wind speed errors are relatively large with a value of 2-3 ms −1 . For AirHARP data, the decrease of the error is less pronounced due to the larger DoLP uncertainties of this sensor (Table 4), and 30 angles are generally needed, with marginal improvement after reaching 40 or more angles. Therefore a total 20-30 angles across all bands, i.e. five to eight viewing angles per band on average, seem sufficient to achieve good retrieval performance. As an example of having 30 continuous angles available for HARP2 measurements, the cloud angular size can be estimated to be 76°following the discussion in Section 3.1.

Remote Sensing Reflectance
The ocean signals are represented by the remote sensing reflectance, defined as the ratio between upwelling waterleaving radiance and the downwelling irradiance just above the FIGURE 5 | Retrieval errors in terms of root mean square error (RMSE) from the retrieval of 1,000 cases for AirHARP and HARP2, and with the two schemes of angle screening as shown in Figure 3. AOD and SSA RMSE are for the wavelength of 550 nm. Note that the retrieval errors for Chla is represented by the MAE error of Chla in log scale as used in Gao et al. (2021), which is different from other quantities. R rs RMSE at all the four wavelengths with AOD (550 nm) < 0.2 are shown in the last row with shaded area indicating the PACE mission R rs uncertainty goal. The horizontal axis, N v , is the total number of angles used in the retrievals, which is equivalent to the total number of unique measurements of 2N v including both reflectance and DoLP across four HARP bands as discussed in Section 3.1.
Frontiers in Remote Sensing | www.frontiersin.org December 2021 | Volume 2 | Article 757832 ocean surface (Mobley et al., 2016). Through atmospheric correction, the remote sensing reflectance can be derived as: where the reflectance contributed by the atmosphere and ocean surface (ρ f atm+sfc ) and the transmittance and BRDF corrections (C BRDF /T f,+ d t f,+ u ) can be evaluated using the NN model in Gao et al. (2021) based on the retrieved aerosol and ocean properties. The viewing angle closest to nadir is used to compute R rs .
Comparing the retrieved R rs with the truth data computed from radiative transfer simulation with Sun at zenith and viewing angle at nadir, the last row in Figure 5 shows R rs uncertainties in terms of RMSE for those simulations with AOD (550 nm) < 0.2. Note that the R rs uncertainties do not include the contribution from the additional random noise added in the synthetic data, therefore it is only used to represent the impacts from the aerosol and ocean surface property retrievals. For the application to PACE data, the calibration uncertainties for OCI is less than 0.5%, which can also be used to cross-calibrate PACE MAP measurements and therefore reduce the impacts of the MAP calibration uncertainties in R rs evaluation. To achieve the PACE mission R rs uncertainty goals (shown in Figure 5), the results suggest we need at least a total of 30 unique HARP angles for R rs at 440 nm, approximately 10-20 angles for all other bands. With larger AOD in the range of [0.3,0.5], the errors increase up to 0.001 even when all angles are used (not shown) as accurate atmospheric correction becomes more challenging due to the increased signal from aerosols.

AIRHARP FIELD DATA RETRIEVALS
The ACEPOL field campaign conducted aerosol and cloud measurements from the NASA high altitude ER-2 aircraft from October to November of 2017 over a variety of scenes including oceans, urban areas, intensive agriculture, forests and high desert around California, Nevada and Arizona . The field campaign involved four MAP instruments including AirHARP, AirMSPI, SPEX airborne, and RSP (Table 1), as well as two lidar sensors-HSRL-2 (Burton et al., 2015) and CPL (the Cloud Physics Lidar) (McGill et al., 2002). Aerosol retrievals have been conducted based on the MAP measurements from ACEPOL (Fu et al., 2020;Gao et al., 2020;Puthukkudy et al., 2020;Gao et al., 2021;Hannadige et al., 2021). Simultaneous aerosol and ocean color retrievals have been performed over three AirHARP scenes . In this study we focus on AirHARP Scene where cirrus clouds were reported to impact the aerosol retrievals. A portion of AirHARP viewing angles impacted by the water condensation on the instrument front lens has been removed from the dataset.
The RGB images of AirHARP scene at three different viewing directions with θ a − 20°, 0°, 20°are shown in Figure 7A1-C1. Two cloud patches are indicated in Figure 6 as A and B. Due to the parallax effect (the data are geolocated to the surface while the clouds are elevated), A and B appear at different locations when viewed from the three angles. Since our data screening criteria are applied to each angle and wavelength combination, the approach is not affected by the data being geolocated to an altitude other than cloud top. As discussed in the next section, a multi-angle cloud mask can be developed based on data geolocated at the surface.

Adaptive Data Screening
To identify and remove cloud-affected view angles, we apply the data screening approach (MAPDS) as explained in Figure 1B. Figure 7 provides an example comparing the measurement and fitted reflectance and DoLP from one pixel near cloud A after the first pass of retrievals. Here σ ρ and σ P indicate the uncertainties used in the cost function, Δρ t and ΔP t indicate the residuals (difference between the forward model and measurements), and Δρ t /σ ρ and ΔP t /σ P the normalized residuals as used in the cost function (Eq. 3).
The blue shaded region was identified as cloud contaminated due to the relatively large ρ t /σ ρ at all bands (with the largest at 870 nm), and the large ΔP t /σ P , at the 670 nm band. As a further buffer, we removed any angles within 4°of angles screened by either reflectance or DoLP criteria at 550 or 670 nm bands. In this way, both reflectance and DoLP measurements in the blue shaded region are excluded. We chose 550 and 670 nm bands as reference based on their higher data quality comparing with other bands  which can be adjusted for application to other datasets. For 440 and 870 nm bands outside of the blue and red shaded region, the screening criteria from Eq. 8 are still applied independently for reflectance and DoLP, which results in mostly DoLP angles removed. The final retrieval (at the third retrieval pass) is done using the remaining measurements after data screening and resulted in a value of χ 2 reduced from 4.98 to 0.6, and the retrieved total AOD reduced from 0.064 to 0.052.
For the AirHARP data used in Gao et al. (2021), the measurement uncertainties were larger than expected at 440 and 870 nm. In that study, DoLP measurements from the 440 nm band were excluded manually. In this study, the adaptive screening approach is able to automatically identify and remove the problematic angles at 440 nm band, providing an automatic data quality control mechanism. Figure 8 shows retrieval results with and without the adaptive data screening. As seen in Figure 8A1, the maximum number of total viewing angles is 120 around the center upper region. The elongated region with large χ 2 values in Figure 8C1 is due to the impacts of the cirrus cloud (Figure 7). Clouds can be viewed from a broader range of pixels in the along-track direction, which impacts the retrieval performance. The retrieved AOD at 550 nm in Figure 8D1 has a clear covariation with χ 2 with the highest AOD around clouds A and B. The R rs at 550 nm are mostly smooth with small variations that seems to also relate to cloud patterns. Therefore, the maximum number of viewing angles for DoLP is around 100 as shown in Figure 8B1 (recall as discussed above DoLP at 440 nm was excluded from retrievals in Figure 8). Figure 8A2-E2 applied the criteria of N v > 30 and χ 2 < 2 to ensure Frontiers in Remote Sensing | www.frontiersin.org December 2021 | Volume 2 | Article 757832 the retrieval quality, removing most of the pixels around the clouds due to the large χ 2 values and leaving an elongated holes in the retrieved AOD and R rs . After applying the MAPDS screening process, the number of available angles around clouds are reduced ( Figure 8A3, B3) and χ 2 ( Figure 8C3) becomes much more uniform. Similarly, as seen in Figure 8D3, the prominent large AOD regions around points A and B are reduced. Most of the cirrus cloud features for R rs from Figure 8E1 are removed as shown in Figure 8E3. There are still minor remaining artifacts for R rs in Figure 8E3 which can be further reduced by choosing a smaller χ 2 threshold in the analysis. Note that the data screening process also removed some of the original viewing angles closest to nadir and the new R rs can be computed using an angle further away which helps reduce the impacts of cirrus cloud on R rs . Future improvement may require using multiple angles to analyze R rs . After applying the criteria of N v > 30 and χ 2 < 2 in Figure 8A4-E4, the pixels removed are mostly around the edge of the image due to the small number of available viewing angles. The N v criterion is important to ensure a sufficient number of measurements used in the retrievals.

Retrieval Results
As shown in Figure 8D2, the criteria N v > 30 and χ 2 < 2 result in many retrievals being discarded when all angles are used in the processing. After applying the data screening, the χ 2 are greatly reduced; this reduction in χ 2 outweighs the decrease in N v , with the effect that most of these discarded retrievals are now retained ( Figure 8D4). The valid retrieval pixels almost doubled after MAPDS. These pixels mostly have N v larger than 50; synthetic retrievals in Figure 5 imply an AOD retrieval uncertainty of 0.03 for AOD up to 0.5. For smaller AOD around 0.1, the retrieval uncertainty can be as low as 0.01 .
The retrieved AirHARP AOD at 550 nm is compared with the HSRL AOD at 532 nm in Figure 9 for two lines (A1-A2) across FIGURE 7 | (A,B) Fitted reflectance and DoLP (circles) and AirHARP measurements (dots) for a pixel near cloud A after the first pass of retrievals; (C,D) uncertainty model as defined in Eqs 4 and 5; (E,F) fitting residuals; (G,H) normalized fitting residuals. The blue shaded region was identified by the data screening method as cloudcontaminated. Most DoLP measurements at 440 and 870 nm at viewing angles larger than 0°are also removed. Measurements within the red shaded region were removed due to water condensation on the instrument lens.
Frontiers in Remote Sensing | www.frontiersin.org December 2021 | Volume 2 | Article 757832 point A, and (B1-B2) across point B. The mean and standard deviation after averaging 4 × 4 pixels (2.2 km × 2.2 km) areas are used in the comparison, and only pixels satisfying N v > 30 and χ 2 < 2 are considered. To mitigate the influence from atmospheric turbulence in the field measurement, the HSRL AOD is estimated by multiplying an assumed lidar ratio of 40sr with the aerosol backscatter coefficient derived from the HSRL technique, and the corresponding AOD systematic uncertainty is estimated to be 50% (Fu et al., 2020). An additional cirrus cloud mask, which uses thresholds of backscatter ratio greater than 1 (backscatter ratio is the ratio of the particulate backscatter to molecular backscatter) and particulate depolarization greater than 0.2 above 8 km, has been applied on the HSRL-2 data. The AirHARP AOD fills in the data gaps due to the cirrus cloud mask between HSRL AOD results with a few overlapped pixels. The retrievals with all angles in Figure 9B1 show differences of 0.019 between the averaged FIGURE 8 | Retrievals on the AirHARP scene over ocean without data screening (all angles, first and second columns) and with MAPDS (third and fourth columns). Panels show the number of viewing angles (N v ) for the reflectance (A1-A4) and DoLP (B1-B4), retrieval cost function (χ 2 , C1-C4), retrieved AOD (D1-D4) and retrieved R rs (sr −1 , E1-D4). The first and third columns show all retrieved pixels; the second and the fourth columns show pixel satisfying criteria of N v >30 and χ 2 <2, respectively. The HSRL AOD at 532 nm is indicated as colored dots in the AOD plots (D1-D4). Same as Figure 6, two cloud patches are indicated by A and B.
Frontiers in Remote Sensing | www.frontiersin.org December 2021 | Volume 2 | Article 757832 AOD from AirHARP and HSRL, but after applying the data screening, these differences reduce to 0.005 in Figure 9B2. The averaged difference shown in (a1) and (a2) are similar to each other with a value of 0.01. However, a peak in Figure 9A1 with a magnitude of 0.17 is reduced to a value of 0.07 after the angles influenced by cloud are removed in Figure 9A2. Similarly, the peaks in Figure 9B1 with a value up to 0.15 is reduced to below 0.09 as shown in Figure 9B2. However, there are small variations of the retrieved AOD in the vicinity of cloud locations, which may be related to light scattering from the 3D structure of cloud and/or aerosol humidification effects. The small inhomogeneity of R rs near clouds may be due to the artifacts from the aerosol retrievals or the presence of cirrus clouds. The pixels which are influenced by the cirrus cloud can be identified from the remaining number of measurement angles (N v ).

Retrieval Uncertainties and Cost Function
Retrieval uncertainties can be evaluated in two ways: 1) the truth-in and truth-out studies as for the synthetic data in Section 3, where the retrieval errors are evaluated by comparing the retrieval results with the truth data in the synthetic data; 2) Using error propagation based on the Jacobian matrix. The second method is effective to evaluate the retrieval uncertainties and have been used to study MAP measurements  and more broadly in aerosol remote sensing (Sayer et al. (2020) for a review). However, it represents a best-case uncertainty, as it is based on various assumptions including that our knowledge on the uncertainty models in the cost function (Eq. 3) are accurate, that the retrievals converge to their global minimum, and that the forward model is close to linear near the solution (see Povey and Grainger (2015) for a review). Retrieval first guesses are important, as convergence to local minima when starting from different initial values can be a problem for high-dimensional retrievals  The histogram of χ 2 for retrievals with all angles and with the adaptive screening are shown in Figure 10, with most probable values of 1.85 and 1.35 respectively. Because these values are larger than 1, the retrieval uncertainty budget for σ ρ and σ P ( Table 2) may underestimate the real total uncertainties. The total number of reflectance and DoLP measurements used in the retrieval are 112 and 85 on average according to Figure 8A1, B1, and Figure 8A3, B3, respectively. The theoretical χ 2 distribution from Eq. 7 with degree of freedom of 85 can represent the histogram for MAPDS well after normalizing its most probable value to the same value of 1.35 as for the histogram. However, the width is slightly smaller, which may indicate overfitting of the data relating to the uncertainty correlation FIGURE 9 | Mean (colored dots) and standard deviation (bars) of the retrieved AirHARP AOD at 550 nm averaged over 4×4 pixels (2.2km ×2.2 km) for two along track lines across cloud patch A (line A, A1-A2) and cloud patch B (line B, B1-B2). HSRL AODs at 532 nm are indicated by circles. For the retrieval results with all angles (A1-B1), most pixels converge into χ 2 values larger than 2, but for retrievals with MAPDS (A2-B2) most pixels converged within χ 2 <2. The χ 2 value for the corresponding AirHARP AODs are indicated by its color label in the plots.
Frontiers in Remote Sensing | www.frontiersin.org December 2021 | Volume 2 | Article 757832 between angles. The retrieval results using all angles have a much wider distribution which cannot be approximated by theoretical χ 2 distribution. The fact that the histogram after adaptive screening matches the theoretical expectation more closely indicates that the rejected angles led to outliers as they did not fit the assumed forward model and (Gaussian) uncertainty assumption. To evaluate retrieval uncertainties using error propagation in future studies, the total uncertainties may need to be multiplied by a factor of 1.35 due to this cost mismatch. A further complication, not widely addressed in the remote sensing community, is that this work (and many others) assumes no spectral or angular correlation between measurement/forward model uncertainties, or between errors in nearby pixels. As discussed by Sayer et al. (2020), this has consequences in terms of uncertainty estimates and cost function when these correlations exist in the real world. However, these correlations are hard to incorporate into retrievals because their magnitudes are, in most cases, poorly known. The presence of correlations also complicates the theoretical χ 2 distribution, and a more complex generalized χ 2 distribution should ideally be used to determine optimal cost thresholds for data filtering, if correlations can be determined.

MAPDS as a Multi-Angle Cirrus Cloud Mask
The RGB images in Figure 6 indicated two cloud patches, which can be identified and removed by MAPDS. In this sense the screening can serve as an angle-by-angle cloud mask for multi-angle data. Figure 11A1-C1 shows the reflectance at 670 nm for the same three viewing angles of Figure 6, and indicates the angles identified by MAPDS at 670 nm, which correspond to clouds A and B very well. The lengths for cloud A and B in the along-track direction are estimated to be around 5-8 km. The displacement between the cloud A locations in Figure 11A2, C2 is estimated through the cross correlation of the masked reflectance from Figure 11A2, C2 in the along-track lines, which shows a change of distance (Δd) of 7-9 km. With the angular span of 2 × 20°, the height of the cloud can be estimated as ctan (20°) ×Δd/2 which is 10-12 km. This assumes no significant cloud motion in the along-track direction during the time interval between the measurements at the two viewing directions. Note that we have assumed the cloud is not moving during the time interval when conducting measurements at the two viewing directions. A more precise treatment considering cloud motion may be conducted by involving more viewing angles at different observing times. The screening could be used to generate a cloud fraction by dividing the number of viewing angles impacted by clouds by the total number of angles. For PACE, OCI's viewing angle will be 20°, and SPEXone viewing angles include both ±20°viewing angles (and other three angles), therefore a HARP2-derived MAPDS cloud mask could be used to complement cloud masks generated by OCI and/or SPEXone. This possibility will be assessed in future studies when PACE data are available.

Data Quality Control
In this study, we discussed the removal of angles for the HARP measurements by filtering out view angles that are problematic for aerosol retrievals to improve data quality. Analogous concepts may be applied to other instruments, such as in the hyperspectral rather than angular domains for OCI and SPEXone. The removed measurements can also be used to identify instrument artifacts or cases involving excess amounts of sunglint due to insufficient information on wind direction, surface slopes, and other factors (Gao et al., 2019). The data with better-modeled glint characteristics can then be used for further analysis in simultaneous aerosol and glint retrievals .

CONCLUSION
In this study we developed and analyzed a data thinning technique to improve performance of aerosol and ocean color retrievals from hyper-angular MAPs by screening and removing problematic measurements affected by cloud and other anomalies. We investigated the impact of the number of viewing angles on retrieval uncertainty for the AirHARP and HARP2 instruments based on synthetic data, finding that a total of 20-30 unique angles across all bands (five to eight viewing angles per band in average) were sufficient to achieve good retrieval performance. Therefore, the total number of viewing angles of HARP2 (90) and AirHARP (120) allows for some screening while retaining sufficient observations for high-quality retrievals.
We further developed and applied an automatic adaptive data screening approach called MAPDS to the AirHARP field measurements. We found that it effectively identified and removed angles influenced by thin cirrus clouds, increasing the level of agreement with independent AOD data from airborne HSRL-2, and providing an additional angle-by-angle cloud mask. We also found that this screening resulted in a better match of retrieval cost to theoretical χ 2 distributions. To improve processing efficiency and accuracy, deep learning techniques including neural network and automatic differentiation are explored. The FastMAPOL algorithm using neural network has demonstrated a factor of 1,000 speed improvement. After we replaced numerical calculation of Jacobians with automatic differentiation of the neural network forward model, the retrieval speed is improved by another factor of 10, more than compensating for the additional overhead of the screening process. The FastMAPOL algorithm and the adaptive screening approach provide efficient ways to process multiangle polarimetric measurements for NASA's Plankton, Aerosol, Cloud, ocean Ecosystem (PACE) mission, and can be applied to other modern remote sensing studies where a large amount of data collected.

AUTHOR CONTRIBUTIONS
MG, KK, BF, and P-WZ formulated the original concept for this study. MG developed the algorithm and generated the scientific data. P-WZ developed the radiative transfer code used to train the NN model. KK, AS, and OH advised on the uncertainty analysis and cloud mask. BF, AI and P-WZ advised on the ocean products. KK, P-WZ, AS, BC and YH advised on the aerosol products. VM and XX provided and advised on the AirHARP data. RF, SB, and MF provided and advised on the HSRL-2 data. All authors provided critical feedback and contributed to the final manuscript.

FUNDING
MG, KK, BF, BC, AI, AS, and P-WZ have been supported by the NASA PACE project. P-W has been supported by NASA (grants 80NSSC18K0345 and 80NSSC20M0227). Funding for the ACEPOL field campaign came from NASA (ACE and CALIPSO missions) and SRON. Part of this work has been funded by the NWO/NSO project ACEPOL (project no. ALWGO/16-09). The ACEPOL campaign has been supported by the NASA Radiation Sciences Program.