Feng Xu
email: feng.xu@jpl.caltech.edu

David J Diner
email: david.j.diner@jpl.nasa.gov

Oleg Dubovik
email: oleg.dubovik@univ-lille.fr

Yoav Schechner

remote sensing

Keywords: correlated aerosol inversion, radiative transfer, multiangle radiometry, polarimetry

Aerosol retrieval algorithms used in conjunction with remote sensing are subject to ill-posedness. To mitigate non-uniqueness, extra constraints (in addition to observations) are valuable for stabilizing the inversion process. This paper focuses on the imposition of an empirical correlation constraint on the retrieved aerosol parameters. This constraint reflects the empirical dependency between different aerosol parameters, thereby reducing the number of degrees of freedom and enabling accelerated computation of the radiation fields associated with neighboring pixels. A cross-pixel constraint that capitalizes on the smooth spatial variations of aerosol properties was built into the original multi-pixel inversion approach. Here, the spatial smoothness condition is imposed on principal components (PCs) of the aerosol model, and on the corresponding PC weights, where the PCs are used to characterize departures from the mean. Mutual orthogonality and unit length of the PC vectors, as well as zero sum of the PC weights also impose stabilizing constraints on the retrieval. Capitalizing on the dependencies among aerosol parameters and the mutual orthogonality of PCs, a perturbation-based radiative transfer computation scheme is developed. It uses a few dominant PCs to capture the difference in the radiation fields across an imaged area. The approach is tested using 27 observations acquired by the Airborne Multiangle SpectroPolarimetric Imager (AirMSPI) during multiple NASA field campaigns and validated using collocated AERONET observations. In particular, aerosol optical depth, single scattering albedo, aerosol size, and refractive index are compared with AERONET aerosol reference data. Retrieval uncertainty is formulated by accounting for both instrumental errors and the effects of multiple types of constraints.

Introduction

Aerosol retrievals performed on remotely-sensed radiometric and polarimetric imagery are subject to non-uniqueness of the solutions and a large computational burden. The former is caused by insufficient information in the observations and the latter is associated with the high dimensionality of the parameter space. To stabilize the inversion process, optimization-based inversions are often informed by extra constraints that supplement the observations. Exploration and application of these constraints are primary tasks in designing a reliable inversion algorithm. A particularly effective physical constraint for remote sensing applications is the so-called smoothness condition (see [START_REF] Phillips | A technique for numerical solution of certain integral equation of first kind[END_REF][START_REF] Twomey | On the numerical solution of Fredholm integral equations of the first kind by the inversion of the linear system produced by quadrature[END_REF][START_REF] Twomey | Introduction to the Mathematics of Inversion in Remote Sensing and Indirect Measurements[END_REF][START_REF] Tikhonov | On the solution of incorrectly stated problems and a method of regularization[END_REF][START_REF] Tikhonov | Solution of Ill-Posed Problems[END_REF]), which takes advantage of the general observation that some quantities (such as ozone vertical distribution and atmospheric temperature profile) vary smoothly in certain dimensions. For aerosols, their loading and properties (e.g., size distribution or spectral real and imaginary refractive index) tend to vary smoothly in the spatial (horizontal and vertical) dimensions, while the underlying surface reflectance tends to vary smoothly in the temporal dimension. Imposition of this type of constraint leads to the "multi-pixel" inversion approach [START_REF] Dubovik | Statistically optimized inversion algorithm for enhanced retrieval of aerosol properties from spectral multi-angle polarimetric satellite observations[END_REF]. In this case, many pixels are stably inverted at one time and the dimensionality of the retrieval space is proportional to the number of inverted pixels.

Additional constraints can further reduce the large dimensionality of the parameter space, with the expected benefits of higher stability and accelerated retrievals. In this study, we introduce a "correlation constraint" on the aerosol inversions. This constraint takes advantage of empirical dependencies between different sets of aerosol parameters. The cross-correlation coefficient between a pair of aerosol fields quantifies the degree of correlation-as long as this coefficient is not zero, the aerosol variables are not independent of each other, i.e., they are correlated to some extent. However, development of a sophisticated physical model to parameterize these correlations is intractable due to the complexity of the real physical world. On the other hand, empirical methods, such as principal component analysis (PCA) or empirical orthogonal function (EOF) analysis, can be used to capture phenomenological correlations between aerosol fields. By implementing PCA over a training dataset, the original fields can be transformed to superpositions of a set of mutually orthogonal vectors, or principal components (PCs). Then, the original aerosol fields can be approximately reconstructed using a linear combination of a certain number of PCs, with accuracy dependent on the number of PCs employed.

Figures 1 and2 show a PC analysis of AERONET Level 2 aerosol inversions within a circular domain of diameter 2000 km around the Fresno, California, and Namibe, Angola, sites, respectively. The aerosol fields include total volume concentration, fraction of fine mode particles, effective radii of fine and coarse mode particles, spectral real and imaginary parts of refractive index, and volume concentration of the spherical particles. To generate Figure 1, a total of 691 retrievals from 60 AERONET stations within a circular domain of 2000 km around the AERONET site in Fresno, California, at longitude -119.787 • and latitude 36.738 • for the years 1994-2018 were used in the PC analysis. For Figure 2, a total of 780 retrievals from 18 AERONET stations within a circular domain of 2000 km around the AERONET site in Namibe, Angola, at longitude 12.178 • and latitude -15.159 • for the years 2000-2017 were used for the PC analysis. Although there are 13 parameters and the aerosol fields vary both temporally and spatially, a few PCs dominate the variations about the mean. For these two cases, 85-90% of the variance in the aerosol fields is captured by the first four PCs. Corresponding to the AERONET analysis in Figure 1 (Fresno domain), Table 1 lists the mutual correlation coefficients among all 13 parameters. Many different types of parameters show correlation coefficients exceeding ~0.3, and among related parameters, such as real and imaginary part of the refractive index, even stronger correlation is observed. While non-absorbing and weakly absorbing particles are dominant in the Fresno domain (Figure 1), absorbing (smoke) aerosols prevail in the Namibe domain (Figure 2). Nevertheless, mutual correlation among aerosol parameters is found in the latter domain as well, as indicated in Table 2. Moreover, comparison of Tables 1 and2 shows dependence of correlation on specific sets of parameters and regions. For example, the correlation between total aerosol concentrations and spherical particle concentration in the 2000 km Fresno region is weaker than that in the Namibe region. The correlation between fine mode aerosol fraction and real part of refractive index, however, is stronger in Fresno region.

intractable due to the complexity of the real physical world. On the other hand, empirical methods, such as principal component analysis (PCA) or empirical orthogonal function (EOF) analysis, can be used to capture phenomenological correlations between aerosol fields. By implementing PCA over a training dataset, the original fields can be transformed to superpositions of a set of mutually orthogonal vectors, or principal components (PCs). Then, the original aerosol fields can be approximately reconstructed using a linear combination of a certain number of PCs, with accuracy dependent on the number of PCs employed. effective radii of fine ( r eff, fine ) and coarse ( r eff, coarse ) mode aerosols, column effective real part of refractive indices ( n r, 1-4 ) at 0.439, 0.675, 0.870, and 1.018 m, respectively, imaginary part of refractive indices ( n i, 1-4 ), and spherical particle volume concentration ( C v,sphere ). Though not specified by use of legend, each color is associated with an independent set of retrieval parameters. Top right panel: spatial and temporal mean of all retrievals and the first four principal components (PCs). Bottom left panel: percentage variance of the aerosol fields captured by PCs, indicating that 41%, 62%, 77%, and 85% of variance is captured by the first one, two, three, and four PCs, respectively. Bottom right panel: Aerosol optical depth (AOD) and single scattering albedo (SSA) at 440 nm that are reported in the AERONET retrievals analyzed here. (b) Regression of the derived aerosol properties from four PC against input AERONET values. The scatter plot is colored by the density of the points. Each panel corresponds to a specific aerosol property indicated in the title. AERONET inversion of aerosol fields and properties (in natural logarithm space) including total volume concentration (C v,tot ), fraction of fine mode aerosols ( f fine ), effective radii of fine (r eff, fine ) and coarse (r eff, coarse ) mode aerosols, column effective real part of refractive indices (n r, 1-4 ) at 0.439, 0.675, 0.870, and 1.018 µm, respectively, imaginary part of refractive indices (n i, 1-4 ), and spherical particle volume concentration (C v,sphere ). Though not specified by use of legend, each color is associated with an independent set of retrieval parameters. Top right panel: spatial and temporal mean of all retrievals and the first four principal components (PCs). Bottom left panel: percentage variance of the aerosol fields captured by PCs, indicating that 41%, 62%, 77%, and 85% of variance is captured by the first one, two, three, and four PCs, respectively. Bottom right panel: Aerosol optical depth (AOD) and single scattering albedo (SSA) at 440 nm that are reported in the AERONET retrievals analyzed here. (b) Regression of the derived aerosol properties from four PC against input AERONET values. The scatter plot is colored by the density of the points. Each panel corresponds to a specific aerosol property indicated in the title. This empirical correlation among aerosol parameters motivates reduction of the dimensionality of the retrieval parameter space and acceleration of the radiative transfer (RT) computations. Benefiting from the empirical aerosol correlations and orthogonality of PCs, increased computational efficiency is achieved by developing a PC-based fast multi-pixel RT computation scheme, as described in Section 4. A straightforward way to reduce dimensionality of the retrieval space is to first derive PCs from a training dataset, such that the retrieval problem reduces to determination of the pixel-dependent PC weights. To account for possible errors in the precomputed PCs, the correlation-based inversion should allow adjustment of the elements of the PC vectors. This is not straightforward, as the vector elements vary significantly in magnitude. A methodology for handling this issue is discussed in Section 3.

PCA has been used to significantly accelerate hyperspectral RT computations [START_REF] Liu | Principal component-based radiative transfer model for hyperspectral sensors: Theoretical concept[END_REF][START_REF] Natraj | A fast linearized pseudo-spherical two orders of scattering model to account for polarization in vertically inhomogeneous scattering-absorbing media[END_REF][START_REF] Spurr | Linearization of the principal component analysis method for radiative transfer acceleration: Application to retrieval algorithms and sensitivity studies[END_REF]. Following a determination of representative components from a training dataset of trace-gas and solar-induced chlorophyll fluorescence, the hyperspectral inversions of these quantities are sped up by retrieving PC weights [START_REF] Liu | Case-study of a principal-component-based radiative transferforward model and retrieval algorithm using EAQUATE data[END_REF][START_REF] Liu | Retrieval of atmospheric profiles and cloud properties from IASI spectra using super-channels[END_REF][START_REF] Liu | New spectral fitting method for full-spectrum solar-induced chlorophyll fluorescence retrieval based on principal components analysis[END_REF]. In another retrieval application, Multi-angle Imaging SpectroRadiometer (MISR) operational aerosol retrieval over land uses spatial contrasts to derive PCs of the surface related contribution to the top-of-the-atmosphere (TOA) radiances [START_REF] Martonchik | Techniques for the retrieval of aerosol properties over land and ocean using multiangle data[END_REF][START_REF] Martonchik | Retrieval of aerosol properties over land using MISR observations[END_REF]. The similarity in the angular shapes of surface bidirectional reflectance factors among MISR's four spectral bands in the visible and near-infrared [START_REF] Diner | Using angular and spectral shape similarity constraints to improve MISR aerosol and surface retrievals overland[END_REF] is an example of a prior empirical constraint used to stabilize the aerosol retrievals. In contrast to MISR's application of PCA to the surface boundary, this paper applies PCA to the aerosol parameters. An optimization approach is developed, which allows both PC vectors and PC weights to be retrieved.

Our retrieval approach is referred to as the "correlated multi-pixel inversion", where "multi-pixel" reflects that our approach simultaneously uses multiple pixels within image-based measurements (as in the original multi-pixel inversion [START_REF] Dubovik | Statistically optimized inversion algorithm for enhanced retrieval of aerosol properties from spectral multi-angle polarimetric satellite observations[END_REF]), while "correlated" captures the extra correlation constraint. In contrast to most optimization approaches that retrieve the individual aerosol fields (see reviews in [START_REF] Dubovik | Polarimetric remote sensing of atmospheric aerosols: Instruments, methodologies, results, and perspectives[END_REF][START_REF] Kokhanovsky | The modern aerosol retrieval algorithms based on the simultaneous measurements of the intensity and polarization of reflected solar light: A review[END_REF]) from multi-spectral radiometric or polarimetric remote sensing observations, we retrieve PC weights and PC vectors of those fields that have spectral, spatial, or temporal correlation with each other, with the purpose of reducing parameter space and improving algorithm efficiency. The individual fields are then constructed from the PC weights and PC vectors. The smoothness constraints on those fields, as implemented in the original "multi-pixel algorithm" used with airborne data at JPL [START_REF] Xu | Joint retrieval of aerosol and water-leaving radiance from multispectral, multiangular and polarimetric measurements over ocean[END_REF][START_REF] Xu | Coupled retrieval of aerosol properties and land surface reflection using the Airborne Multiangle SpectroPolarimetric Imager[END_REF], are adapted to the current correlated multi-pixel inversion approach. Both "multi-pixel" and "correlated multi-pixel" retrieval approaches utilize a one-dimensional (1D) code based on the independent pixel approximation [START_REF] Cahalan | The albedo of fractal stratocumulus clouds[END_REF] to ensure forward modeling efficiency.

This paper is organized as follows. Following an algorithm overview in Section 2, the correlated multi-pixel retrieval algorithm and error analysis are formulated in Section 3. A faster forward vector (polarized) RT model for a coupled atmosphere-surface system is presented in Section 4. In Section 5, retrieval tests are performed using 27 datasets acquired by the Airborne Multiangle SpectroPolarimetric Imager (AirMSPI) during multiple NASA field campaigns. A summary is presented in Section 6.

General Structure of the Algorithm

The expected advantages of correlated multi-pixel inversion in PC space are two-fold. First, it reduces the number of parameters to be retrieved and utilizes observed correlations as empirical constraints to improve the retrieval efficiency and accuracy. Specifically, the correlated parameter space reduces from N corr × N pixel to (N PC + 1) × N corr + N PC × N pixel , where N pixel , N corr , and N PC are the number of image pixels in the scene, the number of correlated aerosol parameters per pixel, and the number of PC vectors, respectively. Reduction of parameter in conjunction with the imposition of proper constraints is expected to mitigate ill-posedness. Morever, establishment of the PC vectors using ground-based measurements, climatology, or other data sources enables characterization of aerosol properties for which the remote-sensing observations may have insufficient sensitivity (such as aerosol chemical composition or vertical profile).

An overview of the correlated multi-pixel approach and its mapping to different sections of this paper is shown in Figure 3. Symbols and abbreviations used in this work are found in Appendix A (Table A1). A training dataset of correlated parameters is used to derive the spatial and temporal mean of correlated parameters over a targeted area, along with PCs of the correlated departures from the mean. Together with the uncorrelated retrieval parameters, a state vector is initialized with elements specified in Section 3.1. In addition to observational constraints provided by the remote-sensing observations used in the retrievals (Section 3.2.1), convergence and robustness retrieval are optimized by imposing additional constraints. For correlated fields, the additional constraints include: (a) the spectral variation of certain correlated fields (e.g., aerosol optical properties); (b) the spatial variation of certain correlated fields across neighboring image pixels; (c) mutual orthogonality and unit length of the PC vectors; and (d) zero sum of the PC weights. Formalisms for imposing these constraints in PC space are provided in Sections 3.2.2-3.2.7 and Appendix C. For any uncorrelated field (such as surface reflectance), constraints on the spectral, spatial, or temporal variations of the associated parameters are handled in the same way as is formulated in the original multi-pixel inversion approach [START_REF] Dubovik | Statistically optimized inversion algorithm for enhanced retrieval of aerosol properties from spectral multi-angle polarimetric satellite observations[END_REF] and repeated in Section 3.2.3 and Appendix B. Combining observational and a priori constraints, a system of equations is established at each iteration and then solved to increment the solution (Section 3.2.8) in an inversion process. Error analysis for all retrieved properties is discussed in Section 3.3.

Figure 3.

General structure of the correlated multi-pixel inversion approach. The interpretation of symbols used in the figure can be found in Table A1 of Appendix A.

To preserve generality of the approach, we use the symbol x corr in the notations below to denote a column vector containing an arbitrary set of parameters that have correlation with each other. In the specific implementation of the correlated multi-pixel inversion for AirMSPI aerosol remote sensing (Section 5), we specify the correlation to be among aerosol parameters only, while the surface parameters are uncorrelated with each other or with aerosol properties (though in principle these assumptions could be relaxed [START_REF] Hou | An algorithm for hyperspectral remote sensing of aerosols: 1. Development of theoretical framework[END_REF][START_REF] Hou | An algorithm for hyperspectral remote sensing of aerosols: 2. Information content analysis for aerosol parameters and principal components of surface spectra[END_REF][23]). Utilization of the aerosol correlations enables development of a fast PC-RT model for TOA radiance and polarization calculation over a group of pixels (Section 4).

As the retrieval output, the correlated aerosol fields for pixel p ( x corr, p ) are constructed from pixel-resolved PC weights in vector form ( w p ), spatially effective PC matrix (v), and spatial and temporal mean vector x corr , namely

x corr, p  x corr  v  w p (1)
where x corr is a single column vector and the matrix v is made up of N PC column vectors consisting of correlated aerosol fields, with N PC denoting the number of retrieved principal components. By integrating all PCs into a matrix, we have [START_REF] Twomey | On the numerical solution of Fredholm integral equations of the first kind by the inversion of the linear system produced by quadrature[END_REF] A priori constraints: On the correlated parameters through PCs and on the uncorrelated parameters (Sections 3.2.2-3.2.7 and Appendices B -D) State vector: x q including (x̄c orr , v, w) q for correlated parameters and uncorrelated parameters x uncorr,q (Section 3.1)

Uncorrelated parameters:

x uncorr,q=0

Numerical inversion: Improves fitting f * with iteratively updated solution (Section 3.2.8)

Forward radiative transfer that simulates the observation f(x q ) (Section 4)

Remote sensing definition:

Spectral, angular information, and geometries

x q x retrieved
Remote sensing data: Spectral, angular, and/or polarimetric observations f * and uncertainties f * for a group of pixels

General structure of inversion algorithm

f(x q ), K(x q ) Figure 3. General structure of the correlated multi-pixel inversion approach. The interpretation of symbols used in the figure can be found in Table A1 of Appendix A.

To preserve generality of the approach, we use the symbol x corr in the notations below to denote a column vector containing an arbitrary set of parameters that have correlation with each other. In the specific implementation of the correlated multi-pixel inversion for AirMSPI aerosol remote sensing (Section 5), we specify the correlation to be among aerosol parameters only, while the surface parameters are uncorrelated with each other or with aerosol properties (though in principle these assumptions could be relaxed [START_REF] Hou | An algorithm for hyperspectral remote sensing of aerosols: 1. Development of theoretical framework[END_REF][START_REF] Hou | An algorithm for hyperspectral remote sensing of aerosols: 2. Information content analysis for aerosol parameters and principal components of surface spectra[END_REF][23]). Utilization of the aerosol correlations enables development of a fast PC-RT model for TOA radiance and polarization calculation over a group of pixels (Section 4).

As the retrieval output, the correlated aerosol fields for pixel p (x corr,p ) are constructed from pixel-resolved PC weights in vector form (w p ), spatially effective PC matrix (v), and spatial and temporal mean vector x corr , namely

x corr,p = x corr + v × w p (1)
where x corr is a single column vector and the matrix v is made up of N PC column vectors consisting of correlated aerosol fields, with N PC denoting the number of retrieved principal components. By integrating all PCs into a matrix, we have

v = v 1 v 2 • • • v N PC (2)
where each PC is composed of N corr elements, i.e.,

v k = [v k (1); v k (2); . . . ; v k (N corr )] (3)
where ";" indicates that the elements v k (1), v k (2), . . . and v k (N corr ) are vertically arranged into a column vector. We can further define a PC weight matrix w that consists of N pixel columns, namely,

w = w 1 w 2 • • • w N pixel (4)
where the vector containing PC weights for p-th pixel is

w p = w p (1); w p (2); • • • ; w p (N PC ) (5)
The quantities (x corr , v, w) are contained in the solution vector x from the last iteration of an inversion process.

Inversion

State Vector

Equations ( 1)-( 5) demonstrate the basic operations for constructing correlated fields from PC weight and PC vectors. We define a state vector x state that includes both correlated and uncorrelated fields. The retrieval parameters are arranged in the following order,

x state = x corr ; v k=1 ; v k=2 ; • • • v k=N PC ; w p=1 ; x p=1,uncorr ; w p=2 ; x p=2,uncorr • • • ; w p=N pixel ; x p=N pixel ,uncorr (6)
where the image-effective mean of aerosol fields takes the form of a column vector

x corr = x corr,1 ; x corr,2 ; • • • ; x corr,N TP,corr (7)
and N TP,corr is the total number of types of parameters correlated with others, and each type may have a subset of values. For example, real and imaginary refractive indices are two types of parameters correlated with each other and each of them have a subset of values as a function of wavelength. For the j-th type of correlated parameter with L(j) elements,

x corr,j = [x corr,j (1); x corr,j (2);

• • • x corr,j (L(j))] (8)
Similarly, an arbitrary (k th ) PC vector composed of N TP,corr types of correlated parameters can be arranged into the following column vector,

v k = v k,1 ; v k,2 ; • • • ; v k,N TP,corr (9)
where for the type of correlated parameter with L corr (j) elements,

v k,j = v k,j (1); v k,j (2); . . . ; v k,j (L corr (j)) (10)
The PC weight vector for pixel p is

w p = w p (1); w p (2); • • • w p (N PC ) (11)
For the state vector with N TP,uncorr types of uncorrelated parameters, we have for an arbitrary pixel p,

x p,uncorr = x p,uncorr,1 ; x p,uncorr,2

where the [n TP,uncorr ]-th type of uncorrelated parameter has L uncorr (n TP,uncorr ) elements, namely,

Constraints

Reduction of parameter space is critical for ensuring inversion efficiency and mitigating ill-posedness of large-scale optimization. Moreover, imposition of various types of constraints is important for stabilizing the retrievals and improving accuracy. The correlated multi-pixel inversion algorithm described here uses a total of ten types of constraints:

1.

The first type of constraint, formulated in Section 3.2.1, consists of the observations provided directly by the remote sensing instrument(s).

2.

The second type of constraint consists of a priori values for the retrieval parameters, and is described in Section 3.2.2.

3.

For the uncorrelated fields (such as surface reflection properties), the pixel-resolved properties x uncorr can be subjected to both across-pixel and within-pixel constraints. Imposition of these types of constraints has been incorporated into the original multi-pixel inversion [START_REF] Dubovik | Statistically optimized inversion algorithm for enhanced retrieval of aerosol properties from spectral multi-angle polarimetric satellite observations[END_REF] and is repeated in Section 3.2.3 as well as in Appendix B as the third type of constraint. 4.

When a set of parameters (e.g., aerosol properties) are correlated with each other, their mean can be subjected to smoothness constraints. This is referred to as the fourth type of constraint in Section 3.2.3.

5.

Transformation of smooth variations of aerosol properties from regular aerosol parameters into the PC space forms the fifth type of constraint, discussed in Section 3.2.4.

6.

In PC space, the smoothness constraints can be applied to the PC weights w and vectors v separately. Application to across-pixel weights w, discussed in Section 3.2.3, forms the sixth type of constraint. 7.

Similarly, application of the smoothness constraints to certain type of parameters within a PC vector v k , also discussed in Section 3.2.3, forms the seventh type of constraint. Although it appears that the sixth and seventh types of constraints applied to "w" and "v" separately are redundant with the fifth type of constraint that ensures a smooth variation of overall correlated field constructed from PCs; they are helpful when poor initial guesses of "w" and "v" are provided by a training dataset. 8.

The eighth type of constraint, formulated in Section 3.2.5, imposes a zero sum of the PC weights. 9.

The ninth type of constraint, formulated in Section 3.2.6, imposes mutual orthogonality of the PC vectors. 10. The tenth type of constraint, formulated in Section 3.2.7, imposes unit norm on all PC vectors.

Even with this set of constraints, numerical errors are unavoidable during the iterations. Hence, a post-correction is implemented after each iteration by reapplying PC analysis to the updated PCs. This way the intrinsic properties of PCs are strictly preserved.

With the ten types of constraints (M = 10) introduced above, we now describe the statistical inversion of multi-source (or multi-constraint) data. It involves solving the following system of equations (see [START_REF] Dubovik | Optimization of numerical inversion in photopolarimetric remote sensing[END_REF] for Equations ( 14)-(20) below):

f * i = f i (x) + ∆f * i 1 ≤ i ≤ M (14)
where f * i denotes the i th type of constraint, ∆f * i is the error with this type of constraint, and x = x state , as defined in Equation [START_REF] Dubovik | Statistically optimized inversion algorithm for enhanced retrieval of aerosol properties from spectral multi-angle polarimetric satellite observations[END_REF]. Formally, the statistical independence of different sources of constraints means that the covariance matrix of joint constraint

f * = [f * 1 ; f * 2 ; • • • ; f * M ] has the following structure C f * =       C 1 0 0 0 0 C 2 0 0 • • • • • • . . . • • • 0 0 C M       (15)
where C i indicates the covariance matrix of i-th constraint (f * i ). Following the expressions of f * i and C f * i , the probability distribution function (PDF) of joint data (1 ≤ i ≤ M) can be derived by multiplying PDFs of data from all M sources, namely,

P(f(x)|f * ) = I ∏ i=1 P(f i (x)|f * i ) ∼ exp - 1 2 I ∑ i=1 [f i (x) -f * i ] T (C i ) -1 [f i (x) -f * i ] (16)
Further introducing the weight matrix (W) for multiple (M) types of constraints, the objective cost function to be minimized has a quadratic form, namely,

Ψ total (x) = M ∑ i=1 γ i Ψ i (x) (17)
where

Ψ i (x) = 1 2 [f i (x) -f * i ] T W -1 i [f i (x) -f * i ] (18)
W i = 1 ε 2 i C i ( 19
)
γ i = ε 2 1 ε 2 i ( 20
)
In the above equations, ε 2 i is the first diagonal element of C i (i.e., ε 2 i = {C i } 11 ) and the Lagrange factor γ i weights the contribution of each type of constraint with respect to the first one (γ 1 = 1).

Minimization of Ψ total (x) in Equation ( 17) means its gradient with respect to the solution x approaches zero, such that

∇Ψ total (x) = M ∑ i=1 γ i ∇Ψ i (x) = 0 ( 21
)
which can be ensured by enforcing the gradient of all components approach zero, namely

∇Ψ i (x) = K T i W -1 i (f i (x) -f * i ) = 0 ( 22
)
where K i is the Jacobian matrix containing the derivatives of i-th type of constraint with respect to the retrieval parameters. To solve the above equation iteratively, we replace x by x -∆x in Equation [START_REF] Hou | An algorithm for hyperspectral remote sensing of aerosols: 2. Information content analysis for aerosol parameters and principal components of surface spectra[END_REF] and substitute

f i (x -∆x) = f i (x) -K i ∆x (23)
into it. This results in

∇Ψ i (x -∆x) = K T i W -1 i (f i (x) -K i ∆x -f * i ) = 0 ( 24
)
or equivalently,

(K T i W -1 i K i )∆x = K T i W -1 i (f i (x) -f * i ) = ∇Ψ i (x) (25)
The Jacobian matrix K i in Equations ( 22)-( 25) consists of the derivative of the l-th observational or a priori data with respect to the n-th unknown,

K i,(l,n) = ∂f i,l ∂x n x (26)
More explicit evaluation of f i (x), f * i , and K and W matrices for all ten types of constraints is discussed in the following subsections.

Observational Constraints (i = 1)

Assuming a total number of N pixel pixels, where each pixel has Z observations, then after arranging all pixel data into a column vector we have

f * i=1 (x) = [f * 1,p 1 ; f * 1,p 2 ; • • • ; f * 1,N pixel ] N×1 = [y 1 ; y 2 ; • • • ; y N f ] (27)
where

N f = N pixel × Z.
The parameters of the atmosphere-surface model are adjusted so that the model prediction of radiance and polarization f 1 (x) fit the observational constraints f * 1 (x). The calculation of f 1 (x) is introduced in Section 4. The Jacobian matrix K 1 consists of first order partial derivatives with respect to correlated and uncorrelated parameters in the vicinity of x, namely,

K 1 =       K 1,xcorr K 1,v K 1,w (p 1 ) K 1,uncorr (p 1 ) 0 0 • • • 0 0 K 1,xcorr K 1,v 0 0 K 1,w (p 2 ) K 1,uncorr (p 2 ) • • • 0 0 • • • • • • • • • • • • • • • • • • . . . • • • • • • K 1,xcorr K 1,v 0 0 0 0 • • • K 1,w (N pixel ) K 1,uncorr (N pixel )       (28)
where K 1,x corr is the Jacobian matrix containing derivatives of observations with respect to spatial and temporal mean correlated parameters (total number is N corr ). Variation of x corr impacts observations in all pixels. Then K 1,x corr is evaluated by

K 1,x corr =         ∂y 1 ∂x corr (1)
∂y 1 ∂x corr (2) • • • ∂y 1 ∂x corr (N corr ) ∂y 2 ∂x corr (1)
∂y 2 ∂x corr (2) • • • ∂y 2 ∂x corr (N corr ) • • • • • • . . . • • • ∂y N f ∂x corr (1)
∂y N f ∂x corr (2) • • • ∂y N f ∂x corr (N corr )         (29)
The matrix K 1,v in Equation ( 28) is the Jacobian matrix containing derivatives of observations with respect to all PC elements. It is evaluated as follows:

K 1,v = K 1,v (k = 1) K 1,v (k = 2) • • • K 1,v (k = N PC ) (30)
with

K 1,v (k) =         ∂y 1 ∂v k (1) ∂y 1 ∂v k (2) • • • ∂y 1 ∂v k (N corr ) ∂y 2 ∂v k (1) ∂y 2 ∂v k (2) • • • ∂y 2 ∂v k (N corr ) • • • • • • . . . • • • ∂y N f ∂v k (1)
∂y N f ∂v k (2) • • • ∂y N f ∂v k (N corr )         (31)
The matrix K 1,w in Equation ( 28) is the Jacobian matrix containing derivatives of observations with respect to pixel-resolved PC weights and is evaluated as follows:

K 1,w (p) =         ∂y p,1 ∂w p (1) ∂y p,1 ∂w p (1) • • • ∂y p,1 ∂w p (N PC ) ∂y p,2 ∂w p (2) ∂y p,2 ∂w p (2) • • • ∂y p,2 ∂w p (N PC ) • • • • • • . . . • • • ∂y p,Z ∂w p (1)
∂y p,Z ∂w p (2) • • • ∂y p,Z ∂w p (N PC )         (32)
Moreover, the derivatives with respect to pixel-resolved uncorrelated parameters are evaluated by,

K 1,uncorr (p) =          ∂y p,1 ∂x p,1 ∂y p,1 ∂x p,2
• • • • • . . . • • • ∂y p,Z ∂x 1 ∂y p,Z ∂x p,2 • • • ∂y p,Z ∂x p,Nuncorr          (33)
The covariance matrix for the first constraint (i.e., the observations) assembles the sub-weighting matrix from all pixels, namely,

W 1 =       W 1 (1) 0 • • • 0 0 W 1 (2) • • • 0 • • • • • • . . . • • • 0 0 • • • W 1 (N pixel )       (34) 3.2.2. A Priori Constraints (i = 2)
The a priori constraint is constructed in the same way as in a previous study [START_REF] Dubovik | Optimization of numerical inversion in photopolarimetric remote sensing[END_REF]. Namely, Equation [START_REF] Martonchik | Retrieval of aerosol properties over land using MISR observations[END_REF] becomes

f * i=2 = x a priori = x + ∆x a priori (35)
Then in Equations ( 21)-( 25), K i=2 = I (identity matrix) and

W i=2 = 1 ε 2 a *
C a * . More explicitly, W i=2 can be constructed from estimated range of each parameter relative to the first one,

W i=2 =         1 0 • • • 0 0 (x 2,max -x 2,min ) 2 (x 1,max -x 1,min ) 2 • • • 0 • • • • • • . . . • • • 0 0 • • • (x N,max -x N,min ) 2 (x 1,max -x 1,min ) 2         (36)
3.2.3. Smoothness Constraints in Regular Parameter Space (i = 3, 4, 6, and 7)

The formulation of these types of constraints is similar, so we discuss them together. The third type of constraint reflects the smooth variation of an uncorrelated type of parameter (such as the variation of parameter of the surface bidirectional reflectance distribution function with wavelength). The fourth type of constraint reflects the smooth variation of certain type of parameter in the mean field (e.g., aerosol refractive index as a function of wavelength). The sixth type of constraint reflects the smooth variation of PC weights from pixel to pixel. The seventh equation reflects the smooth variation of certain type of parameters residing in a PC vector (such as the deviation of refractive index from the mean). The strengths of the sixth and seventh types of constraints depends on the rank of PC. Putting all four of these types of constraints together, we have

           f * 3 = 0 * = S 3,m x uncorr + ∆ * g(x uncorr ) f * 4 = 0 * = S 4,m x + ∆ * g(x corr ) f * 6 = 0 * = S 6,m w state + ∆ * g(w) f * 7 = 0 * = S 7,m v state + ∆ * g(v) (37)
where w state and v state are column vectors containing only PC weights and vectors, respectively, and are extracted from the overall state vector expressed in Equation ( 6), and S i,m is the differentiation matrix of m-th order for i-th type of constraint.

When a retrieval parameter varies smoothly as a function of some variable z, it is assumed to be locally approximated by a smooth function g(z), such as a constant, a line, a parabola, etc. With a polynomial form, the m-th derivative approaches zero [START_REF] Dubovik | Statistically optimized inversion algorithm for enhanced retrieval of aerosol properties from spectral multi-angle polarimetric satellite observations[END_REF], such that

S i,m (z)z = d m g(z) dz m = 0 ⇒      g m=1 (z)= const g m=2 (z) = Az + B g m=3 (z) = Az 2 + Bz + C (38)
For a discretized grid of the variable z, the explicit form of S i,m z, is,

S i,m (z j )z j =      dg dz ≈ ∆ 1 g ∆ 1 (z) = g(z j+1 )-g(z j ) ∆ 1 (z j )
, for m= 1

d m g dz m ≈ ∆ m g ∆ m (z) = ∆ m-1 g(z j+1 )/∆ m-1 (z j+1 )-∆ m-1 g(z j )/∆ m-1 (z j ) [∆ m-1 (z j )+∆ m-1 (z j+1 )]/2
, for m ≥ 2 [START_REF] Clark | USGS digital spectral library splib06a[END_REF] Taking the orders of difference m = 1 and 2 as examples, we have

∆ m=1 (z j ) = z j+1 -z j ∆ m=2 (z j ) = [∆ 1 (z j ) + ∆ 1 (z j+1 )]/2 (40)
Application of the above equation to L discretized grids r j (namely, 1 ≤ j ≤ L) leads to

f i (x) = S i,m x (41)
so that by invoking Equation ( 26),

K i = S i,m (42)
where the matrix S i,m is evaluated by,

S i,m=1 =        1 ∆ 1 (1) -1 ∆ 1 (1) 0 • • • 0 0 1 ∆ 1 (2) -1 ∆ 1 (2) • • • 0 0 0 • • • . . . 0 0 0 • • • 1 ∆ 1 (L-1) -1 ∆ 1 (L-1)        (43)
and

S i,m=2 =        2 ∆ 1 (1)[∆ 1 (1)+∆ 1 (2)] -2 ∆ 1 (1)∆ 1 (2) 2 ∆ 1 (2)[∆ 1 (1)+∆ 1 (2)] 0 0 • • • 0 0 2 ∆ 1 (2)[∆ 1 (2)+∆ 1 (3)] -2 ∆ 1 (2)∆ 1 (3) 2 ∆ 1 (3)[∆ 1 (2)+∆ 1 (3)] • • • 0 • • • • • • • • • • • • • • • . . . • • • 0 0 • • • 2 ∆ 1 (L-1)[∆ 1 (L-1)+∆ 1 (L)] -2 ∆ 1 (L-1)∆ 1 (L) • • • 2 ∆ 1 (L)[∆ 1 (L-1)+∆ 1 (L)]        (44)
The same principle applied to higher orders of difference (m > 2) ensures a smooth curve with d m g(z) dz m = 0. Substitution of Equations ( 41), (42), and f * i = 0 * into Equation [START_REF] Xu | Joint retrieval of aerosol and water-leaving radiance from multispectral, multiangular and polarimetric measurements over ocean[END_REF] gives,

Ψ i (x) = 1 2 x T S T i,m W -1 i,m S i,m x (45)
where the weighting matrix W has the following diagonal terms,

{W i,m } jj = 1 ∆ m (z j ) (46)
and ∆ m (z j ) is specified in Equation ( 40) for m = 1 and 2, and can be generalized to an arbitrary higher order.

Substitution of Equations ( 41), (42), and f * i = 0 * into Equation [START_REF] Levenberg | A method for the solution of certain non-linear problems in Least Squares[END_REF] gives

(S T i,m W -1 i,m S m )∆x = S T i,m W -1 i,m (S i,m x) (47)
Further defining

Ω i = S T i,m W -1 i,m S i,m (48)
we have

Ω i ∆x = Ω i x (49)
Note that for equidistant pixels on an imaged grid, the weighting matrix W is a unity matrix (namely W = I). To distinguish the smoothness constraints imposed on different types of parameters, we note here,

Ω i =          Ω uncorr , i = 3 Ω x corr , i = 4 Ω w corr , i = 6 Ω v corr , i = 7 (50) 3.2.4. Smoothness Constraints in PC Space (i = 5)
The fifth type of constraint reflects the smooth variation of correlated parameters in different dimensions. However, unlike the smoothness constraint directly imposed on relevant uncorrelated parameters, it has to be transformed to the PC vectors and weights. For simplicity, we put all PC weights into a single-column vector w state = w k=1 ; w k=2 ; • • • ; w k=N PC and all PC elements for j-th correlated parameter into another single-column vector

v state = v k=1 ; v k=2 ; • • • ; v k=N PC .
The smooth variation of the j-th correlated aerosol parameter is ensured by,

f * i=5 = 0 * = S i,m x wv + ∆ * g([v;w]) (51)
where x wv = v state w state . To evaluate the differentiation matrix S i,m , we recall that a pixel-resolved correlated parameter is constructed from PCs via Equation ( 1). An arbitrary correlated parameter x corr,p (j) associated with p-th pixel is derived as,

x corr,p (j) = x corr,j + N PC ∑ k=1 w p (k)v k (j) = x corr,j + w p (1)v 1 (j) + w p (2)v 2 (j) + . . . + w p (N PC )v N PC (j) (52)
(a) across-pixel smoothness The across-pixel variation of a correlated parameter is contributed by pixel-resolved PC weights, namely,

∆x corr,p = v × (∆w p ) (53)
Taking the first order of difference (m = 1) in the across-pixel variation of the field as an example, Equation (39) becomes,

S m=1 (z j )z j = dg dz ≈ ∆ 1 g ∆ 1 z ≈ N PC ∑ k=1 v k (j) × w j+1 (k) -w j (k) ∆ 1 (z j ) (54)
where z measures the inter-pixel relative distance. As an example, we take the simplest case of a single PC (k = 1), three pixels (p 1 , p 2 , p 3 ) associated with PC weights [w p1 , w p2 , w p3 ], and two correlated parameters associated with PC elements v k=1 (1) and v k=1 [START_REF] Twomey | On the numerical solution of Fredholm integral equations of the first kind by the inversion of the linear system produced by quadrature[END_REF]. The matrix S i,m is constructed in the following way to ensure smooth variation of the correlated parameter across equidistant (e.g., unity spacing) pixels via f i=5 (x wv ) = S i,m x wv :

f i=5 (x wv ) = S i,m x wv =          1 2      w p1 (k = 1) -w p2 (k = 1) 0 v k=1 (1) -v k=1 (1) 0 0 w p1 (k = 1) -w p2 (k = 1) v k=1 (2) -v k=1 (2) 0 w p2 (k = 1) -w p3 (k = 1) 0 0 v k=1 (1) -v k=1 (1) 0 w p2 (k = 1) -w p3 (k = 1) 0 v k=1 (2) -v k=1 (2)               ×        v k=1 (1) v k=1 (2) w p1 (k = 1) w p2 (k = 1) w p3 (k = 1)        (55) (b) within-pixel smoothness
Within-pixel variation of a parameter (such as refractive index) describes its dependence on a non-spatial variable, e.g., wavelength. For an arbitrary pixel p, the within-pixel variation of a correlated parameter is contributed by PC vectors, namely,

∆x corr,p = (∆v) × w p (56)
Taking the first order of difference (m = 1) as an example, the differentiation matrix S i,m is constructed in the following way to ensure smooth variation of the correlated parameter to x(j) as a function of z across which x(j) is smooth (for example, z is wavelength when x(j) is aerosol refractive index),

S i,m=1 (z j )z j = dg dz ≈ ∆ m g ∆ m z ≈ N PC ∑ k=1 w p (k) × [v k (j + 1) -v k (j)] ∆ 1 z j (57)
Taking the simplest case of a single PC (k = 1), three parameters associated with PC elements v k=1 (1), v k=1 (2) and v k=1 (3) as an example,

f i=5 (x wv ) = S i,m=1 x wv =    1 2   w p (k=1) ∆ 1 (z 1 )
-

w p (k=1) ∆ 1 (z 1 ) 0 v k=1 (2)-v k=1 (1) ∆ 1 (z 1 ) 0 w p (k=1) ∆ 1 (z 2 ) - w p (k=1) ∆ 1 (z 2 ) v k=1 (3)-v k=1 (2) ∆ 1 (z 2 )      ×      v k=1 (1) v k=1 (2) v k=1 (3) w p (k = 1)      (58)
Invoking f * i = 0 * , the cost function in the form of Equation ( 45) is derived for both across-pixel and within-pixel smoothness constraints. Further, substitution of f i=5 (x wv ) = S i,m=1 x wv into Equation [START_REF] Levenberg | A method for the solution of certain non-linear problems in Least Squares[END_REF] gives the expression for Jacobian matrix elements,

K i,(l,n) = ∂(S i,m x wv ) l ∂x wv,n x wv (59) in which the matrix S i,m is a function of x (see Equation (57)), therefore K i = S i,m . Intuitively, one can think of f i=5 (x wv ) = S i,m=1
x wv as a model prediction to fit the "observation" f * i , which is a zero vector. Therefore, the above Jacobian matrix is evaluated in the same way as is done for observational constraints by use of finite difference methodology. The finite difference method is used for evaluating Equation (59).

Invoking f i (x wv ) = S i,m x wv and f * i = 0 * , Equation ( 25) becomes

K T i W -1 i K i ∆x wv = K T i W -1 i S i,m x wv (60) Further defining Ω wv corr,1 = K T i W -1 i K i Ω wv corr,2 = K T i W -1 i S i,m (61)
we have Ω corr,1 ∆x wv = Ω corr,2 x wv (62)

Note that if S i,m is independent of x (cases for 3, 4, 6, and 7), then K can be analytically evaluated from Equation (59), namely,

K i = S i,m (63)
so that

Ω wv corr,1 = Ω wv corr,2 = S T i,m W -1 i S i,m (64)
In this manner, the static smoothness matrices for cases for 3, 4, 6, and 7 (Section 3.2.3) are recovered.

Based on above demonstration, the evaluation of S i,m , K i , and Ω i matrices is generalized from the case of regular (non-correlated) retrieval parameters to the case of correlated parameters represented by PCs in Appendix C. The Ω matrix derived above is accounted for later to solve for an increment of the PC terms during optimization.

Zero-sum Constraint on PC Weights (i = 8)

As an intrinsic property of pixel-resolved PC weights, they are required to sum to zero for an arbitrary set of PC vectors. This constraint is ensured by multiplying the PC weights w in the state vector by an O matrix, such that,

f * i=8 = 0 * = f i=8 (w state ) + ∆ * O (65)
where w state is a column vector defined in Section 3.2.4 and contains the weights of all PCs, and

f i=8 (w state ) is expressed as f 8 (w state ) = Ow state (66)
where

O =        1 0 • • • 0 0 1 • • • 0 • • • • • • . . . 0 0 0 • • • 1        N PC ×(N PC ×N pixel ) , 1 = [1, 1, . . . , 1] 1×N pixel ( 67
)
Further substituting f * 8 =0 * and f 8 (w state ) = Ow state into Equation ( 18), the cost function has the following quadratic form,

Ψ i (w state ) = 1 2 (w T state O T Ow state ) (68)
Substitution of Equation (66) into Equation [START_REF] Levenberg | A method for the solution of certain non-linear problems in Least Squares[END_REF] gives the evaluation of Jacobian matrix elements,

K i,(l,n) = ∂(Ow state ) l ∂w state,n w state = O l,n (69)
Further defining a zero-sum matrix to be,

Ω O = Ω i = O T O (70)
and invoking f 8 (w state ) = Ow state and f * 8 =0 * , Equation (25) becomes,

Ω O ∆w state = Ω O w state (71) 3.2.6. Mutual Orthogonality Constraint among PC Vectors (i = 9)
An intrinsic property of PC vectors is that they are mutually orthogonal. This forms the orthogonality constraint imposed on all pairs of PC vectors via f i (i = 9), namely

f * i=9 =0 * = f i=9 (v state ) + ∆ * Γ (72)
where v state is a column vector defined in Section 3.2.4 and contains all PC vectors, and

f i=9 (v state ) is expressed as f 9 (v state ) = Γv state (73)
where,

Γ =                        v T k=2 v T k=1 0 0 • • • 0 0 v T k=3 0 v T k=1 0 • • • 0 0 • • • • • • • • • • • • • • • • • • • • • v T k=N PC 0 0 0 • • • 0 v T k=1 0 v T k=3 v T k=2 0 • • • 0 0 0 v T k=4 0 v T k=2 • • • 0 0 • • • • • • • • • • • • • • • • • • • • • 0 v T k=N PC 0 0 • • • 0 v T k=2 • • • • • • • • • • • • • • • • • • • • • 0 0 0 0 • • • v T k=N PC -1 v T k=N PC                        (74) Further substituting f * 9 =0 * = [ 0 0 • • • 0 ] T { 1 2 [N PC ×(N PC -1)]}×1
and f 9 (v state ) = Γv state into Equation ( 18), the cost function is derived as,

Ψ i (v) = 1 2 v T state Γ T Γv state (75)
Substitution of Equation (73) into Equation ( 26) gives the evaluation of Jacobian matrix elements,

k i,(l,n) = ∂(Γv) l ∂v state,n v state = Γ l,n (76)
Further defining a zero-sum matrix to be,

Ω Γ = Ω i = Γ T Γ (77)
and invoking f 9 (v state ) = Γv state and f * 9 =0 * , Equation ( 25) becomes,

Ω Γ ∆v state = Ω Γ v state (78) 3.2.7. Unity-norm Constraint PC Vectors (i = 10)
Another intrinsic property of each PC vector is that its norm equals unity (i.e.,

v k 2 = v T k v k = 1)
. This forms the unity constraint imposed on all PC vectors via f i (i = 10), such that

f * i=10 =1 * = f i=10 (v state ) + ∆ * U (79) where f 10 (v state ) = Uv state ( 80
)
and U is expressed as,

U =        v T k=1 0 • • • 0 0 v T k=2 • • • 0 • • • • • • . . . 0 0 0 • • • v T k=N PC        N PC ×(N PC ×N corr ) (81)
Then substituting Equation (80) and f * 10

=1 * = [ 1 1 • • • 1 ] T N PC ×1
into Equation ( 18), the cost function is derived as,

Ψ i (v state ) = 1 2 (v T state U T Uv state -N PC ) (82)
Moreover, substitution of Equation (80) into Equation [START_REF] Levenberg | A method for the solution of certain non-linear problems in Least Squares[END_REF] gives the evaluation of Jacobian matrix elements:

K i,(l,n) = ∂(Uv state ) l ∂v state,n v state (83)
Further defining a unity constraining matrix to be,

Ω U,1 = K T i K i Ω U,2 = K T i U (84)
and invoking f 10 (v state ) = Uv state and f * 10 =1 * , Equation ( 25) becomes equivalent to

Ω U,1 ∆v state = Ω U,2 v state -K T i f * i (85)

Construction of overall Equation System

Accounting for the above ten types of above-specified constraints, the solution to minimizing Equation ( 17) is approached iteratively. The solution from iteration q is incremented as follows,

x q+1 = x q -∆x q (86)
where ∆x q is obtained by accounting for both observational and a priori constraints derived in Sections 3.2.1-3.2.7, A q × ∆x q = ∇Ψ total (x q ) (87)

As the explicit form, we have

A q = K T 1,q W -1 1 K 1,q + γΩ total,1 + γ a W -1 a (88) ∇Ψ total (x q ) = K T 1,q W -1 1 f 1 (x q ) -f * 1 + γΩ total,2 x q + γ a W -1 a (x q -x a priori ) +    0 * -γ U K T 10,q f * 10 0 *    (89)
where "γΩ total " incorporates i = 3-10 types of constraints, namely,

γΩ total,1 = γ uncorr Ω Ra uncorr + γ x corr Ω x,Ra corr + γ wv corr Ω wv,Ra corr,1 + γ w corr Ω w,Ra corr + γ v corr Ω v,Ra corr + γ O Ω Ra O + γ Γ Ω Ra Γ + γ U Ω Ra U,1 γΩ total,2 = γ uncorr Ω Ra uncorr + γ x corr Ω x,Ra corr + γ wv corr Ω wv,Ra corr,2 + γ w corr Ω w,Ra corr + γ v corr Ω v,Ra corr + γ O Ω Ra O + γ Γ Ω Ra Γ + γ U Ω Ra U,2 (90)
where the eight terms on the right-hand-side of the above equation are equal to Ω i with i = 3-10, respectively, and the multipliers γ ... control the strength of these constraints. The smoothness matrices

γ uncorr Ω Ra uncorr , γ x corr Ω x,Ra corr , γ O Ω Ra O , γ Γ Ω Ra Γ , γ U Ω Ra U , γ wv
corr Ω wv,Ra corr , and γ w/v corr Ω w/v,Ra corr act on uncorrelated fields, multi-pixel mean of correlated fields, PC weights, PC vectors, PC vectors, combined PC weights and vectors, and separate PC weights and vectors, respectively. They are essentially equal to

γ uncorr Ω uncorr , γ x corr Ω x corr , γ O Ω O ,γ Γ Ω Γ , γ U Ω U , γ wv
corr Ω wv corr , and γ w/v corr Ω w/v corr , respectively, as evaluated in Sections 3.2.3-3.2.7. However, the elements of these matrices are rearranged into the same dimension to add to each other in assembling γΩ total in Equation (90). In performing the rearrangement, specific locations of relevant parameters in the state vector (see Equation ( 6)) have to be accounted for, and zero values have to be filled in to accommodate the parameters not subjected to a specific constraint. For the uncorrelated parameters and fields, the explicit forms of γ uncorr Ω uncorr as well as γ x corr Ω x corr are given in Appendix B. For the correlated parameters and fields, the smoothness constraints imposed on PC vector and weight retrieval incorporates two types: the coupled type (γ wv corr Ω wv corr ), which acts on combined PC weights and vectors, and the decoupled type γ w/v corr Ω w/v corr , which act on PC weights and vectors separately. Based on previous demonstration of principles for deriving Ω i for i = 5, a comprehensive evaluation of γ wv corr Ω wv corr is provided in Appendix C. Comprehensive evaluations of γ w corr Ω w corr and γ v corr Ω v corr are provided in Appendix D. The imposition of γ w corr Ω w corr and γ v corr Ω v corr are effective when pixel-resolved PC weights themselves or PC elements associated with a type of parameter present certain smooth behavior (usually appearing in low order PCs). Our retrieval tests indicate that a combined use of γ wv corr Ω wv corr and γ w/v corr Ω w/v corr stabilize the inversion of PCs. Ideally, a retrieval is deemed successful when the minimization of the cost function is achieved, such that

Ψ total ≤ (N f + N c + N a * -N a )ε 2 f ( 91
)
where N f , N c , N a , and N a* are the total number of observations, total number of constraints imposed on retrieval (including smoothness constraints in different dimensions, zero-sum constraint over PC weights and orthogonality and unity constraints over PC vectors), total number of retrieval parameters, and total number of a priori estimates of parameters, respectively; and ε 2 f is the expected variance due to measurement errors. In practice, forward RT modeling error and other unmodeled effects can impede realization of the required cost function minimization. Therefore, the retrieval is also terminated when the relative difference of fitting residues with solutions from two successive iterations drops below a user-specified threshold value, ε 2 c , namely,

Ψ total (x q+1 ) -Ψ total (x q ) Ψ total (x q ) ≤ ε 2 c . (92) 3.2

.9. Determination of Lagrange Multipliers

The Lagrange multipliers reflect the strength of smoothness constraints for a given parameter to be retrieved. Each multiplier is defined as,

γ i = ε 2 i=1 /ε 2 i ( 93
)
where ε 2 i are the first diagonal elements of the covariance matrices corresponding to i-th type of constraints. To estimate ε 2 i for a given parameter to be retrieved (f = x), which is a function of z, the most unsmooth known solution f us (z) over the target domain is used, as suggested by Dubovik and King [START_REF] Dubovik | A flexible inversion algorithm for retrieval of aerosol optical properties from Sun and sky radiance measurements[END_REF], namely,

ε 2 = z max z min d m [ f us (z)] d m z 2 dz ( 94
)
where z min and z max specify the lower and upper bound of z. For i = 3, 4, 5, γ i is evaluated using least smooth known solution of correlated parameter x over the target domain. For i = 6 and 7, γ i is evaluated using the least smooth known solution of PC weights and PC vectors, respectively, which are derived from a training dataset over the target domain.

In practical implementation of our algorithm, the Lagrange multipliers are modified in the following way:

γ Final i = N f N i ε 2 f ε 2 f γ i (95)
There are two differences between γ Final i and γ i : 1. The multipliers "N f /N i " are introduced to account for possible redundancy of the measured and a priori data. Considering that ε 2 i is a variance of the error in a single measured or estimated a priori value, if we have N values of a similar kind, the total variance increases in proportion to N. Introducing this coefficient ensures that when there are several kinds of data, the data with fewer values are given comparable weight, as the data type for which there is a greater number of available values.

The multiplier ε 2

f /ε 2 f is introduced with ε 2 f estimated as the dynamic fitting residual during iterations:

ε 2 f (x q ) ≈ Ψ total (x q ) (N f + N c + N a * -N a )ε 2 f (96)
With the multiplier ε 2 f /ε 2 f , the fitting residual is used as an estimate of measurement error variance. As a result, during the first few iterations, the contribution of the a priori term is strongest, and its influence decreases as the retrieval progresses. This is done to ensure mostly monotonic convergence, as in the Levenberg-Marquardt procedure [START_REF] Levenberg | A method for the solution of certain non-linear problems in Least Squares[END_REF][START_REF] Marquardt | An algorithm for least-squares estimation of nonlinear parameters[END_REF]. However, the Levenberg-Marquardt approach does not specify a particular scheme for introducing these terms, rather it relies on the implementer's intuition. Our algorithm requires the fitting errors in the initial iterations to be dominated by model linearization errors as opposed to random measurement errors, because at each iterative step the full forward model is replaced by its linear approximation, the "errors of linearization" decrease as convergence toward the final solution progresses, and they practically disappear, so that ε 2 f becomes equal to ε 2 f . As a result of this adjustment of the Lagrange multiplier, the non-linear iteration becomes significantly more monotonic. The Lagrange multiplier that controls the strength of a priori (i = 2), zero sum (i = 8), orthogonality (i = 9), and unity (i = 10) constraints are dynamically updated using Equations ( 95) and (96). When the solution approaches the truth, the constraining effect of smoothness constraints via the Lagrange factor decreases during the optimization process.

Retrieval Error Estiamte

Instrumental errors, forward modeling errors, and physical constraints imposed on the retrievals need to be taken into account when evaluating retrieval errors. Instrumental errors include random and systematic effects, such as absolute bias, band-to-band error, and camera-to-camera error. If all of these error sources are well-characterized, their error propagation to the retrievals can be derived in a comprehensive way [START_REF] Xu | Coupled retrieval of aerosol properties and land surface reflection using the Airborne Multiangle SpectroPolarimetric Imager[END_REF]. If not, then a practical method is to use the square root of the diagonal terms of the following covariance matrix,

∆x j = (C ∆x ) jj (97)
with

C ∆x = (∆x syst ) × (∆x syst ) T + C ∆x,rand (98)
where the covariance matrix C ∆x,rand accounts for the contributions by random errors in the measurement (ε rand ),

C ∆x rand = A(x true ) -1 ε 2 rand ( 99
)
and systematic error ∆x syst is estimated by

∆x syst = [A(x true ) -1 ∇Ψ total (x true ) ( 100
)
where A and ∇Ψ total are calculated by substituting x true into Equations ( 88) and (89), respectively. As the retrieved solution is the closest estimate of the "true" solution x true , we adopt x true = x retrieved for uncertainty estimation. Uncertainties with uncorrelated parameters are available from Equation (97). Uncertainties with pixel-resolved correlated parameters rely on the propagation of retrieval uncertainties of PCs and PC weights, so that,

∆x retrieved p,corr 2 = (∆x retrieved corr ) 2 + N PC ∑ k=1 [w retrieved p (k)] 2 (∆v retrieved k ) 2 + [∆w retrieved p (k)] 2 (v retrieved k ) 2 (101)
To estimate errors for functions of the retrieved parameters (namely f (x retrieved p ), e.g., aerosol single-scattering albedo, (SSA)), the following chain rule is applied,

∆ f (x retrieved p ) = K T 2 (x retrieved p )C ∆x retrieved p K 2 (x retrieved p ) (102)
where K 2 denotes the Jacobian array containing derivatives of f (x retrieved p ) with respect to x retrieved p . When x a priori is not available, we assume x a priori = x q so that the a priori term on the right-hand-side of Equation (89) for ∇Ψ total is neglected. Indeed, instrument errors are already contained in the observation vector. Therefore, by implementing Equations ( 97) and (98), it is possible that errors are double-counted in the case that they bias the solution in the same direction as the modeling errors, resulting in a conservative error estimate.

Most error estimates involving Jacobians assume that the calculation is representative of the whole solution space and that the retrieval error is linear with measurement error. These two assumptions can be problematic in situations where model or observation errors are large. In these cases, closure tests using synthetic data with combined random and systematic errors to obtain improved error estimates are recommended.

Retrieval Options

A priori or first guess of PCs can be obtained from running PCA on a climatology database, field measurements, or model output. In the case of PCs for which confidence in their contents is low, simultaneous retrieval of PC weights and PC vectors has to be performed. For high confidence PCs, a two-step retrieval can be used to speed up the retrieval. That is, we run PC weight retrieval using the predetermined PCs first and then relax PCs during the retrieval to refine the solution. This relaxation is one way to capture unexpected variance that exceeds that represented by the predetermined PCs. In the extreme situation, where PCs are not representative of the actual aerosol properties, the scene model can be reproduced by using N PCs, where N equals to the number of aerosol parameters. In the absence of confidence in an a priori set of PCs, retrieval efficiency is still desired. In this case, one can run the correlated multi-pixel inversion with 2-3 PC vectors to capture the major variation of the correlated fields and then (a) relax more PCs into the retrieval to capture unexpected variance, or (b) return to the original multi-pixel retrieval to refine the solution. These different options have been built into our algorithm. In addition, to stabilize the retrieval of correlated fields varying significantly in magnitude (even in logarithmic space), we built into the algorithm an option of retrieving scaled PC vectors and weights. In this case, the representation of correlated fields by PCs for p-th pixel (namely Equation ( 1)) is modified as,

x corr,p = x corr + c • (v × w p ) (103)
where • in the above equation denotes an operation of element-wise multiplication, and the constant column vector c has N corr elements and is introduced here to balance or weight the contribution of all types of retrieval parameters in the PC analysis. It is calculated from (a) the standard deviation of a correlated field x varying in spatial or temporal scales (σ spatio-temporal, x ), and (b) its uncertainty estimate (σ e, x ), namely,

C x = σ 2 spatio-temporal, x + σ 2 e, x (104)
Such a scaled PC analysis is used to analyze the AERONET retrievals in Figures 1 and2.

To accommodate the constant column vector, the smoothness matrices for the correlated fields and Lagrange multipliers need a slight modification, which is straightforward and is not discussed here.

Radiative Transfer in a Coupled Atmosphere-Surface System

A fast radiative transfer model for a coupled atmosphere-surface (CAS) system plays a critical role in ensuring inversion efficiency. Where correlation among aerosol parameters exists, this section formulates the development of such a model for simulating TOA radiance for a group of pixels simultaneously. Figure 4 shows a CAS system consisting of a land surface and an aerosol/air-molecule mixed layer. The Markov chain model [START_REF] Xu | Markov chain formalism for polarized light transfer in plane-parallel atmospheres, with numerical comparison to the Monte Carlo method[END_REF] is adopted for modeling RT in the atmospheric part of the CAS system. Coupling of atmospheric radiation and surface reflection is performed by use of the adding method [START_REF] Lacis | A parameterization for the absorption of solar radiation in the Earth's atmosphere[END_REF]. As outlined in Figure 5, the principal-component based radiative transfer (PC-RT) model includes two steps: (1) calculating aerosol light scattering properties (aerosol optical depth, single scattering albedo, and phase matrix), and the reflection and transmission matrices of the atmosphere for a group of pixels by utilizing correlations in aerosol properties (Section 4.1.1);

(2) coupling reflection/transmission matrices of the atmosphere with surface reflection pixel-by-pixel to model the TOA observations (Section 4.1.2). For retrieval use, evaluation of the Jacobian matrix that contains the derivatives of TOA observations with respect to all aerosol and surface parameters is described in Section 4.2.

To further ensure modeling efficiency, the correlated aerosol microphysical properties derived from their PCs are input to a light scattering database [START_REF] Dubovik | Application of spheroid models to account for aerosol particle nonsphericity in remote sensing of desert dust[END_REF] for determination of aerosol scattering properties using interpolation. The modeled scattering properties include aerosol optical depth (AOD), single scattering albedo (SSA), and phase matrix. To model surface reflection, the Rahman-Pinty-Verstraete (RPV) model [START_REF] Rahman | Coupled surface-atmosphere reflectance (CSAR) model 2. Semiempirical surface model usable with NOAA advanced very high resolution radiometer data[END_REF] is used to calculate the unpolarized surface bidirectional reflectance distribution function (BRDF). The parameters of the unpolarized part of the surface model include spectral weights and a few angular shape parameters. A modified microfacet model [START_REF] Litvinov | Models for surface reflection of radiance and polarized radiance: Comparison with airborne multi-angle photopolarimetric measurements and implications for modeling top-of-atmosphere measurements[END_REF] is used to calculate the polarized surface BRDF (pBRDF). The pBRDF model parameters include spectral weight, shadowing width, and slope variance of the microfacets. Combination of the two components gives the full surface reflection model used in AirMSPI aerosol retrievals over land [START_REF] Xu | Coupled retrieval of aerosol properties and land surface reflection using the Airborne Multiangle SpectroPolarimetric Imager[END_REF]. With the AOD and single scattering properties, a fast RT model that utilizes the correlation in aerosol fields and orthogonality of PC vectors is used to model the radiometric and polarimetric observations. Remote Sens. 2018, 10, x FOR PEER REVIEW 27 of 49 where in the above equation denotes an operation of element-wise multiplication, and the constant column vector c has N corr elements and is introduced here to balance or weight the contribution of all types of retrieval parameters in the PC analysis. It is calculated from (a) the standard deviation of a correlated field x varying in spatial or temporal scales ( spatio-temporal, x ), and (b) its uncertainty estimate ( e, x ), namely,

c x   spatio-temporal, x 2  e, x 2 (104)
Such a scaled PC analysis is used to analyze the AERONET retrievals in Figures 1 and2. To accommodate the constant column vector, the smoothness matrices for the correlated fields and Lagrange multipliers need a slight modification, which is straightforward and is not discussed here.

Radiative Transfer in a Coupled Atmosphere-Surface System

A fast radiative transfer model for a coupled atmosphere-surface (CAS) system plays a critical role in ensuring inversion efficiency. Where correlation among aerosol parameters exists, this section formulates the development of such a model for simulating TOA radiance for a group of pixels simultaneously. Figure 4 shows a CAS system consisting of a land surface and an aerosol/air-molecule mixed layer. The Markov chain model [START_REF] Xu | Markov chain formalism for polarized light transfer in plane-parallel atmospheres, with numerical comparison to the Monte Carlo method[END_REF] is adopted for modeling RT in the atmospheric part of the CAS system. Coupling of atmospheric radiation and surface reflection is performed by use of the adding method [START_REF] Lacis | A parameterization for the absorption of solar radiation in the Earth's atmosphere[END_REF]. As outlined in Figure 5, the principal-component based radiative transfer (PC-RT) model includes two steps: (1) calculating aerosol light scattering properties (aerosol optical depth, single scattering albedo, and phase matrix), and the reflection and transmission matrices of the atmosphere for a group of pixels by utilizing correlations in aerosol properties (Section 4.1.1); (2) coupling reflection/transmission matrices of the atmosphere with surface reflection pixel-by-pixel to model the TOA observations (Section 4.1.2). For retrieval use, evaluation of the Jacobian matrix that contains the derivatives of TOA observations with respect to all aerosol and surface parameters is described in Section 4.2. To further ensure modeling efficiency, the correlated aerosol microphysical properties derived from their PCs are input to a light scattering database [START_REF] Dubovik | Application of spheroid models to account for aerosol particle nonsphericity in remote sensing of desert dust[END_REF] for determination of aerosol scattering properties using interpolation. The modeled scattering properties include aerosol optical depth (AOD), single scattering albedo (SSA), and phase matrix. To model surface reflection, the Rahman-Pinty-Verstraete (RPV) model [START_REF] Rahman | Coupled surface-atmosphere reflectance (CSAR) model 2. Semiempirical surface model usable with NOAA advanced very high resolution radiometer data[END_REF] is used to calculate the unpolarized surface bidirectional reflectance distribution function (BRDF). The parameters of the unpolarized part of the surface model include spectral weights and a few angular shape parameters. A modified microfacet model [START_REF] Litvinov | Models for surface reflection of radiance and polarized radiance: Comparison with airborne multi-angle photopolarimetric measurements and implications for modeling top-of-atmosphere measurements[END_REF] is used to calculate the polarized surface BRDF (pBRDF). The pBRDF model parameters include spectral weight, shadowing width, and slope variance of the microfacets. Combination of the two components gives the full surface reflection model used in AirMSPI aerosol retrievals over land [START_REF] Xu | Coupled retrieval of aerosol properties and land surface reflection using the Airborne Multiangle SpectroPolarimetric Imager[END_REF]. With the AOD and single scattering properties, a fast RT model that utilizes the correlation in aerosol fields and orthogonality of PC vectors is used to model the radiometric and polarimetric observations.

Fast Multi-Pixel Polarized RT Modeling in the Atmosphere

At a given iteration step, aerosol optical depth (AOD) and light scattering properties, including SSA and phase matrices, are established for all layers. Then, along with the parameters of the surface reflection, RT modeling is performed to evaluate the TOA radiation field. The basic multi-pixel algorithm, which does not take advantage of correlations between aerosol parameters, runs the RT pixel by pixel. Making use of correlation in aerosol fields and orthogonality of PC vectors, the next section describes a fast way to compute radiation in the atmosphere, which is then coupled with surface reflection to derive the TOA radiation field.

Vector radiative transfer model (part 1):

PC-RT model for reflection (R atmos ) and transmission (T atmos ) matrices of atmosphere for a group of pixels (Section 4.1.1)

Retrieval state vector:

x q including (x̄c orr,q , v q , w q ) & x q,surf Forward model* Principal component (PC) terms (x̄c orr , v q , w q ) of aerosol (correlated) properties Simulated observations f(x q ) and Jacobian matrix K q calculation (Section 4.2) Pixel-scale aerosol fields x q,aer calculated from PCs (Section 3.1)

*

Fast Multi-Pixel Polarized RT Modeling in the Atmosphere

At a given iteration step, aerosol optical depth (AOD) and light scattering properties, including SSA and phase matrices, are established for all layers. Then, along with the parameters of the surface reflection, RT modeling is performed to evaluate the TOA radiation field. The basic multi-pixel algorithm, which does not take advantage of correlations between aerosol parameters, runs the RT pixel by pixel. Making use of correlation in aerosol fields and orthogonality of PC vectors, the next section describes a fast way to compute radiation in the atmosphere, which is then coupled with surface reflection to derive the TOA radiation field.

Fast Multiple-Pixel Radiative Transfer Modeling Utilizing Correlation

The spatial and spectral variation of aerosol fields, such as volume concentration as a function of aerosol size, aerosol refractive index, vertical distribution profile, and volume concentration of spherical particles, are captured by a few (N PC ) dominant PC vectors. Taking advantage of orthogonality of the PCs, the quantities (Y), including total AOD (τ aer,tot ), absorption AOD (τ aer,abs ), and reflection and transmission matrices (R and T) for the above-surface atmosphere, can be expanded into Taylor series for multi-variable vector-valued functions v k . Adopting the second order of approximation (justified later in this section) and finite difference method to calculate derivatives, we have,

Y(x p ) ≈ Y(x) + N PC ∑ k=1 Y(x+δ s ×v k )-Y(x-δ s ×v k ) 2δ s w p (k) + Y(x+δ s ×v k )-2Y(x)+Y(x-δ s ×v k ) 2δ 2 s w 2 p (k) (105)
where x is state vector containing multi-pixel mean aerosol properties and δ s is the scale factor that perturbs a PC vector. The scale factor is empirically determined by accounting for: (a) the accuracy of the two Jacobian related terms on the right-hand-side of the above equation; and (b) the sufficiency of expanding the radiation fields to the second order and using finite number of PCs to capture the variation of radiation across the imaging area. For an AirMSPI imaging area 10 km by 10 km, an empirical value of 0.1 is adopted to use the above equation.

The efficiency gain of implementing the above PC based RT modeling is obvious: only (2 × N PC +1) RT runs are necessary for evaluating the radiance fields for all pixels, which is much smaller than required for pixel-by-pixel RT runs. If the variation of aerosol fields across an imaging area is extremely high, we use an alternative method built into the algorithm, which applies PCA to the difference in the RT fields evaluated by RT runs with low and high numbers of streams ("streams" represent the quadrature points for integrating the contribution of light from different directions to an event of scattering in a new direction), namely,

Y(x p ) -Y LS (x p ) ≈ Y HS (x) -Y LS (x) + N PC ∑ k=1 ∆Y(x+δ s ×v k )-∆Y(x-δ s ×v k ) 2δ s w p (k) + ∆Y(x+δ s ×v k )-2∆Y(x)+∆Y(x-δ s ×v k ) 2δ 2 s w 2 p (k) (106)
where "HS" and "LS" denote high and low stream-based RT runs, respectively, and ∆Y = Y HS -Y LS . In this case, RT runs with low order streams are performed for all pixels and then corrected by Equation (106). This methodology was proposed for a hyperspectral RT simulation [START_REF] Natraj | A fast linearized pseudo-spherical two orders of scattering model to account for polarization in vertically inhomogeneous scattering-absorbing media[END_REF] and adopted here for the multi-pixel RT simulation. Compared to directly applying PCA to RT via Equation (105), implementation of Equation ( 106) increases computational cost but captures more variance in the image-scale radiation fields and achieves higher modeling accuracy. Note that the above evaluation assumes the same view and azimuthal angles for all pixels across the whole image, which is not strictly correct, as view angle varies pixel-by-pixel. To account for such an effect, the R and T matrices are calculated for a few view and azimuthal angular grids for each viewed image. Interpolation is then used to obtain the reflection and transmission matrices for an arbitrary pixel before coupling them with the surface reflection (as we assume no correlation between aerosol and surface reflection properties in this paper). Coupling of atmospheric radiation and surface reflection to get the TOA radiance for fitting observation is formulated in the next section.

Using three PCs (N PC = 3), Figure 6 shows a comparison of PC based RT computation of radiance (expressed as bidirectional reflectance factor, BRF) and degree of linear polarization (DOLP) for 30 super-pixels of an AirMSPI image against the RT computation obtained without using PCs. Here, BRF is defined as πI meas d 2 ES /µ 0 E 0 , where I meas is the measured radiance, d ES is the Earth-Sun distance, µ 0 is the cosine of solar zenith angle, and s is the exo-atmospheric solar irradiance. A super-pixel has the resolution of 1 km and is generated by aggregating 100 by 100 original AirMSPI pixels with 10-m resolution. Running RT model over super-pixels can mitigate the errors from the independent pixel approximation [START_REF] Cahalan | The albedo of fractal stratocumulus clouds[END_REF]. For a variation of aerosol loading up to 100%, the error in computed TOA, BRF, and DOLP from PC-based RT computation is found to be less than ~0.5% and ~0.0025, respectively. These values are smaller than typical instrument errors (e.g., ~4% for radiance and 0.005 for DOLP, which are the instrument requirements for AirMSPI [START_REF] Van Harten | Calibration and validation of airborne multiangle spectropolarimetric imager (AirMSPI) polarization measurements[END_REF]).

Note that the PC-RT model formulated here assumes that aerosol properties are correlated. It is possible, however, that some aerosol fields are not highly correlated. In this case, one can use more PCs to capture the variance. In the worst case scenario of no correlation, the RT calculation defaults to the pixel-by-pixel approach.

where x is state vector containing multi-pixel mean aerosol properties and  s is the scale factor that perturbs a PC vector. The scale factor is empirically determined by accounting for: (a) the accuracy of the two Jacobian related terms on the right-hand-side of the above equation; and (b) the sufficiency of expanding the radiation fields to the second order and using finite number of PCs to capture the variation of radiation across the imaging area. For an AirMSPI imaging area 10 km by 10 km, an empirical value of 0.1 is adopted to use the above equation. The viewing geometries are adopted from one of the scenes acquired by AirMSPI during the ACEPOL field campaign. The unpolarized and polarized surface reflectance is calculated from the retrieved BRDF and pBRDF parameters, respectively.

Coupling Atmospheric Radiation with Surface Reflection

Using the deterministic Markov chain RT model for a plane-parallel atmosphere ( [START_REF] Xu | Markov chain formalism for polarized light transfer in plane-parallel atmospheres, with numerical comparison to the Monte Carlo method[END_REF], or other RT models), we calculate two sets of reflection and transmission matrices for the atmosphere, (R atmos , T atmos ) and (R * atmos , T * atmos ), which result from illumination at the top and bottom of the atmosphere, respectively. We denote the surface reflection matrix by R surf , which consists of a BRDF component (R surf, BRDF ) and a pBRDF component (R surf, pBRDF ). In accordance with the adding method [START_REF] Lacis | A parameterization for the absorption of solar radiation in the Earth's atmosphere[END_REF], two matrices Q and S are defined to account for the interaction between the surface and atmosphere, respectively,

S = N max ∑ n=1 Q n (107) Q n = Q 1 Q n-1 (108)
Q 1 = R * atmos R surf ( 109
)
Then the matrices for downwelling and upwelling diffuse light at the atmosphere-land interface are given by

D = T atmos + S exp(- τ atmos µ 0 ) + ST atmos ( 110
)
where τ atmos is the optical thickness of the whole atmosphere, as contributed by all atmospheric constituents, and the diffuse upwelling light from the surface is calculated as

U = R Surf exp(- τ atmos µ 0 ) + R Surf D (111)
The reflection matrix of the full CAS is,

R CAS = R atmos exp(- τ atmos µ )U + T * atmos U (112)
For simplicity in describing the conceptual scheme, the superscript "m" that denotes Fourier series order is not shown in the above expressions. In actuality, the TOA radiation fields for the CAS system are reconstructed from all orders of Fourier terms, namely,

BRF TOA = π ∞ ∑ m=0 (2 -δ 0m )R (m) CAS,11 cos m(φ v -φ 0 ) (113) qBRF TOA = π ∞ ∑ m=0 (2 -δ 0m )R (m) CAS,21 cos m(φ v -φ 0 ) (114) uBRF TOA = π ∞ ∑ m=0 (2 -δ 0m )R (m) CAS,31 cos m(φ v -φ 0 ) (115) vBRF TOA = π ∞ ∑ m=0 (2 -δ 0m )R (m) CAS,41 cos m(φ v -φ 0 ) (116)
where (φ vφ 0 ) is the relative azimuth angle between the view and illumination directions. Equations ( 113)-(116) describe TOA BRFs for total, linearly polarized, and circularly polarized radiation.

When DOLP is used to fit the observations, it is calculated as DOLP ≈ qBRF 2 TOA + uBRF 2 TOA /BRF TOA , where we neglect the minor contribution of circular polarization related term (vBRF TOA ) when it is excluded from polarimetric measurements.

Jacobian Evaluation

In an iterative way, the optimization algorithm adjusts the state vector that parameterizes aerosol and surface properties to approach the solution. At each iterative step, RT calculations are performed to obtain both the modeled radiation fields and their Jacobians, which describe how the radiation fields vary in response to the small perturbation of state vector components and determine the direction and step size to move for the iterative solution to converge. For an observation vector composed of a series of N f measurements (f * 1 = [y 1 ; y 2 ; • • • ; y N f ]) and a state vector of N a components (x = [x 1 ; x 2 ; • • • ; x N a ]), the Jacobian matrix has the following structure,

K =         ∂y 1 ∂x 1 ∂y 1 ∂x 2 • • • ∂y 1 ∂x Na ∂y 2 ∂x 1 ∂y 2 ∂x 2 • • • ∂y 2 ∂x Na • • • • • • . . . • • • ∂y N f ∂x 1 ∂y N f ∂x 2 • • • ∂y N f ∂x Na         N f ×N a (117)
Each matrix element is evaluated by use of the finite difference method. Namely, the derivative of modeled i-th data with respect to n-th state vector component is evaluated by,

∂y i ∂x n = y i (x n + ∆x n ) -y i (x n -∆x n ) 2∆x n (118)
The above equation applies to the calculation of the derivatives of τ aer,tot , τ aer,abs , R, and T with respect to the mean aerosol properties and pixel-resolved PC weights. Then, without applying the finite difference method again, the derivative of quantity Y = {τ aer,tot , τ aer,abs , R, T} at pixel p with respect to an element in k-th PC vector can be derived in a fast manner, namely, from the derivatives with respect to the mean and to the PC weights associated with p-th pixel,

∂Y p ∂v k (i) ≈ w p (k) ∂Y p ∂x(i) (119)
Such a strategy further improves the Jacobian evaluation efficiency when multiple PCs are retrieved.

Inversion of Aerosol and Surface Properties

Section 3 provided a general algorithm formulation to retrieve correlated and uncorrelated parameters. This section describes practical specifics in using AirMSPI data to retrieve aerosol and surface properties. In doing so, we assume correlations in aerosol properties, no correlation in surface properties, and no correlation between aerosol and surface properties.

The first guesses of PCs, PC weights, and mean aerosol fields are derived from a training dataset, which contains all historical AirMSPI retrievals for selected scenes with the criteria of well-calibrated measurements, clear-sky conditions, and collocated ground observations. The aerosol fields include spectrally dependent real and imaginary parts of the aerosol refractive index, volume concentrations of multiple aerosol size components, nonspherical particle fraction, and Gaussian profile-based parameterization of aerosol layer height and standard deviation [START_REF] Xu | Coupled retrieval of aerosol properties and land surface reflection using the Airborne Multiangle SpectroPolarimetric Imager[END_REF]. Constrained by the information content in polarimetric observations, we retrieve aerosol properties by assuming a single layer, and an "effective" set of aerosol properties for the single layer that fits the polarimetric observations is derived.

With these assumptions, the correlated aerosol parameters include volume concentration as a function of five size components (C v 1-N sC with N SC = 5 as adopted in [START_REF] Xu | Coupled retrieval of aerosol properties and land surface reflection using the Airborne Multiangle SpectroPolarimetric Imager[END_REF]), real and imaginary parts of refractive index (n r and n i ), Gaussian distribution-based vertical profile parameterized by central height h a and standard deviation s a , and the volume fraction of spherical particles ( f v, sphere ). Then, the parameter spaces described at the beginning of Section 3 are specified as follows,

x corr = log     x T corr,1 C v1 • • • C v N SC , x T corr,2 n r, λ1 • • • n r, λ7 , x T corr,3 n i, λ1 • • • n i, λ7 , x corr,4 h a , x corr,5 s a , x corr,6 f v, sphere     T (120) v = v 1 v 2 . . . v N PC (121) with v k =    v T k,1 v k,C v1 • • • v k,C v5 , v T k,2 v k,n r, λ1 • • • v k,n r, λ7 , v T k,3 v k,n i, λ1 • • • v k,n i, λ7 , v k,4 v k,h a , v k,5 v k,s a , v k,6 v k, f v,
with N TP,uncorr = 6 types of uncorrelated parameters, including the BRDF spectral weight a λ (j = 1), anisotropy parameter k λ (j = 2), anisotropy parameter g λ (j = 3), pBRDF weight λ (j = 4), shadowing width k γ (j = 5), and slope variance σ s (j = 6), namely,

x p,uncorr, j=1 = a p (λ 1 ); a p (λ 2 ); • • • ; a p (λ 7 ) (125)
x p,uncorr, j=2 = k p (λ 1 );

k p (λ 2 ); • • • ; k p (λ 7 ) (126)
x p,uncorr, j=3 = g p (λ 1 );

g p (λ 2 ); • • • ; g p (λ 7 ) (127)
x p,uncorr, j=4 = ε p (λ 1 );

ε p (λ 2 ); • • • ; ε p (λ 7 ) (128)
x p,uncorr, j=5 = k γ,p (129)

x p,uncorr, j=6 = σ 2 s,p (130)
In Equation ( 120), the natural logarithm is used to ensure non-negativity of the solution after dynamic updates during the iterative optimization process. Though not noted here, the angular shape parameter "g" with the RPV model needs to be offset by a constant before being transformed into logarithmic space. An overview of correlated multi-pixel inversion of aerosol properties and surface reflection algorithm flow is shown in Figure 7. As part of the state vector, the PCs of the correlated aerosol microphysical properties are derived from a training dataset from climatology or other sources. As another part of the state vector, uncorrelated parameters (such as surface reflection properties) are initialized to be static [START_REF] Xu | Coupled retrieval of aerosol properties and land surface reflection using the Airborne Multiangle SpectroPolarimetric Imager[END_REF]. The inversion is stabilized by applying a priori constraints on smooth variations of certain aerosol and surface properties in spatial and spectral dimensions. Iterations repeat until convergence is achieved. For the retrieval tests in this study, it takes five to seven iterations for the solution to converge.

As demonstrated in earlier POLDER and AirMSPI retrievals [START_REF] Dubovik | Statistically optimized inversion algorithm for enhanced retrieval of aerosol properties from spectral multi-angle polarimetric satellite observations[END_REF][START_REF] Xu | Coupled retrieval of aerosol properties and land surface reflection using the Airborne Multiangle SpectroPolarimetric Imager[END_REF], imposition of temporal constraints on the variation of surface reflectance can improve aerosol and surface retrievals. In this case, multiple scenes acquired from revisits of a target have to be used. This functionality is not turned on in the following retrieval tests to simplify the demonstration. An overview of correlated multi-pixel inversion of aerosol properties and surface reflection algorithm flow is shown in Figure 7. As part of the state vector, the PCs of the correlated aerosol microphysical properties are derived from a training dataset from climatology or other sources. As another part of the state vector, uncorrelated parameters (such as surface reflection properties) are initialized to be static [START_REF] Xu | Coupled retrieval of aerosol properties and land surface reflection using the Airborne Multiangle SpectroPolarimetric Imager[END_REF]. The inversion is stabilized by applying a priori constraints on smooth variations of certain aerosol and surface properties in spatial and spectral dimensions. Iterations repeat until convergence is achieved. For the retrieval tests in this study, it takes five to seven iterations for the solution to converge.

As demonstrated in earlier POLDER and AirMSPI retrievals [START_REF] Dubovik | Statistically optimized inversion algorithm for enhanced retrieval of aerosol properties from spectral multi-angle polarimetric satellite observations[END_REF][START_REF] Xu | Coupled retrieval of aerosol properties and land surface reflection using the Airborne Multiangle SpectroPolarimetric Imager[END_REF], imposition of temporal constraints on the variation of surface reflectance can improve aerosol and surface retrievals. In this case, multiple scenes acquired from revisits of a target have to be used. This functionality is not turned on in the following retrieval tests to simplify the demonstration.

AirMSPI Datasets

The retrieval algorithm is designed to retrieve column aerosol and surface reflection properties from observations by AirMSPI, which flies on NASA's ER-2 aircraft at an altitude of 20 km and operates in eight spectral channels: 355, 380, 445, 470*, 555, 660*, 865*, and 935 nm, with the asterisk denoting polarimetric bands, in which the Stokes parameters Q and U are measured in addition to radiance I. Images of each targeted area were obtained at 9 viewing angles: 0° (nadir), ±29°, ±48°, ±59°, and ±66° in AirMSPI's step-and-stare mode. At nadir, the imaged area covers a 10 km x 11 km region and the data are mapped to a 10-m spatial grid. Without using the water-vapor influenced band (935 nm), a total of 117 signals per pixel are used, which include 63 radiances (transformed to logarithmic space in retrieval) at nine angles and seven spectral bands, and 27 signals of q=Q/I and Forward RT model Jacobian matrix (K q ) Construct constraints on a): (x̄q ,corr , w q , v q ) and x q,surf regarding smooth aerosol and surface property variations b): sum-to-zero of PC weights and orthogonality&unity of PCs x q = [x q,aer ; x q,surf ], where x q,aer is derived from (x̄q ,corr , w q , v q ) via Eq.( 1)

x q+1 = x q -t q ×x q q=0: Initialized PC vectors v q =[v q,1 ; v q,2 ; …; v q,N_{PC} ], PC weights w q =[w q,1 , w q,2 , .. w q,N_{PC} ], mean fields (x̄q ,corr ) Initial estimate of uncorrelated surface parameters (x q,surf )

PCA of correlated fields in a training dataset from climatology or other sources

Update x̄q +1,corr ,v q+1 and w q+1 in x q+1 by reapplying PCA to the constructed correlated fields Solve Eq.( 87) for x q Spectral, angular & polarimetric observations from a groups of pixels

x q = x q+1 t q =1 Evaluate cost function  q+1 2 Evaluate modeled data and cost function  q 2 No Yes Output x retrieved = x q+1  q+1 2 <  q 2 - q+1 2 / q 2 < Forward RT model t q = t q /2
Yes Evaluate uncertainty for all retrieval output

AirMSPI Datasets

The retrieval algorithm is designed to retrieve column aerosol and surface reflection properties from observations by AirMSPI, which flies on NASA's ER-2 aircraft at an altitude of 20 km and operates in eight spectral channels: 355, 380, 445, 470*, 555, 660*, 865*, and 935 nm, with the asterisk denoting polarimetric bands, in which the Stokes parameters Q and U are measured in addition to radiance I. Images of each targeted area were obtained at 9 viewing angles: 0 • (nadir), ±29 • , ±48 • , ±59 • , and ±66 • in AirMSPI's step-and-stare mode. At nadir, the imaged area covers a 10 km × 11 km region and the data are mapped to a 10-m spatial grid. Without using the water-vapor influenced band (935 nm), a total of 117 signals per pixel are used, which include 63 radiances (transformed to logarithmic space in retrieval) at nine angles and seven spectral bands, and 27 signals of q = Q/I and another 27 signals of u = U/I in the three polarimetric bands. Retrievals for all pixels of a surface area viewed from all 9 angles are performed simultaneously. The measurement errors are adopted as 4% for radiance (to account for angle-to-angle and pixel-to-pixel uncertainties) and 0.005 for DOLP.

A wide range of atmospheric conditions and terrestrial environments have been covered by AirMSPI during over a hundred flights in several airborne campaigns. In this paper, we use AirMSPI data from the Polarimeter Definition Experiment (PODEX) (January to February 2013), Studies of Emissions, Atmospheric Composition, Clouds, and Climate Coupling by Regional Surveys (SEAC 4 RS) (August to September 2013), CalWater-2 (January to March 2015), and Imaging Polarimetric Assessment and Characterization of Tropospheric Particulate Matter (ImPACT-PM) (July 2016) campaigns. From these, 27 AirMSPI step-and-stare data collection sequences were identified to be cloud-free and collocated with AERONET sun photometers for retrieval validation. Locations of these AERONET sites and AirMSPI/AERONET measurement times can be found in a previous study [START_REF] Xu | Coupled retrieval of aerosol properties and land surface reflection using the Airborne Multiangle SpectroPolarimetric Imager[END_REF]. To control

Retrieval Validation against AERONET Products

The retrieved AirMSPI AOD, SSA, size distribution, and refractive index are validated against AERONET Level 1.5 aerosol products. We choose Level 1.5 AERONET product (cloud-screened and quality controlled), as it reports more fields to validate our retrievals. Though the Level 2.0 AERONET product (quality-assured) has higher confidence, the SSA, refractive index, and aerosol size distribution fields were not all generated in 27 test cases. However, for the fields reported by both versions, negligible differences were observed. As a check, a retrieval is performed using AirMSPI observations acquired on September 9, 2013, over the AERONET Baskin, Louisiana, site which is located at longitude = -91.738 • and latitude = 32.282 • . The left image in Figure 8a shows nadir radiance using the spectral band combination of 445, 555, and 660 nm, while the right image displays DOLP at 470, 660, and 865 nm. The retrievals are performed over the area viewed in common at all 9 AirMSPI view angles, outlined by the yellow box. Figure 8b shows TOA BRF at 445, 555, and 660 nm in the left, middle, and right panels, respectively, for the retrieval area with spatial resolution ~1.0 km. Figure 8c shows maps of retrieved AOD, SSA, and surface albedo (A surf ) at 555 nm in the left, middle, and right panels, respectively. The retrieved AOD, SSA, and volume weighted aerosol size distribution for the atmosphere above the super-pixel closest to the Baskin site is compared to the AERONET reference data in the left, middle, and right panels of Figure 8d. Good agreement (quantified below) is obtained for all of these quantities, except for the coarse particle size distribution, likely due to the lack of bands longward of 1000 nm in AirMSPI. Generally, the difference between AirMSPI retrievals of AOD, SSA, and size distribution and AERONET reference data is within their retrieval uncertainties. The AirMSPI uncertainties plotted in Figure 8d are estimated as the root mean square of the retrieval uncertainties of these aerosol quantities and the standard deviation of their variations over the whole image. The AERONET uncertainties consist of two parts: temporal variation within the ±~1-h window centered on the AirMSPI nadir overpass time, and aerosol measurement and retrieval error [START_REF] Dubovik | Accuracy assessment of aerosol optical properties retrieval from AERONET sun and sky radiance measurements[END_REF]. A temporally closest AERONET reference data was identified compare to AirMSPI retrieval at the spatially closest pixel. To account for airmass change during the measurements, the ±~1-h window centered on AirMSPI nadir overpass time is used to calculate the AERONET uncertainty from temporal variation. Figure 9 shows a comparison of retrieved pixel-resolved AODs, SSAs, and aerosol size distributions (dV/dln(r), in μm 3 /μm 2 ) from the correlated multi-pixel inversion with those retrieved using original multi-pixel algorithm adapted for AirMSPI [START_REF] Xu | Coupled retrieval of aerosol properties and land surface reflection using the Airborne Multiangle SpectroPolarimetric Imager[END_REF]. The AOD results for seven AirMSPI spectral bands are plotted in different colors: pink (355 nm), purple (380 nm), dark blue (445 nm), light blue (470 nm), green (555 nm), red (660 nm), and brown (865 nm). Linear regression is performed to obtain slope a, intercept b, as well as the coefficient of determination R 2 . The mean Figure 9 shows a comparison of retrieved pixel-resolved AODs, SSAs, and aerosol size distributions (dV/dln(r), in µm 3 /µm 2 ) from the correlated multi-pixel inversion with those retrieved using original multi-pixel algorithm adapted for AirMSPI [START_REF] Xu | Coupled retrieval of aerosol properties and land surface reflection using the Airborne Multiangle SpectroPolarimetric Imager[END_REF]. The AOD results for seven AirMSPI spectral bands are plotted in different colors: pink (355 nm), purple (380 nm), dark blue (445 nm), light blue (470 nm), green (555 nm), red (660 nm), and brown (865 nm). Linear regression is performed to obtain slope a, intercept b, as well as the coefficient of determination R 2 . The mean absolute difference (MAD) in AERONET and AirMSPI results is also calculated to measure the overall deviation. Strong correlation and low bias (R 2 ≥ 0.88, a ~0.90, b ≤ 0.05, and MAD < 0.02) are observed. It can also be observed that a variation of aerosol loading by ~30% around the mean (namely in the range 0.28 ≤ AOD 445nm ≤ 0.40 with mean value 0.32) across the retrieval area is captured by the correlated multi-pixel inversion. Implementation of our using several datasets with even higher (~90%) variation of aerosol loading over several smoke scenes acquired by AirMSPI during the recent Aerosol Characterization from Polarimeter and Lidar (ACEPOL) campaign showed the algorithm to be capable of capturing this variation. The regression in Figure 9b shows correlations and low bias of SSAs (R 2 > 0.40, a > 0.60, b < 0.030, and MAD ≤ 0.004) from the two inversions as well. Figure 9c shows basic consistency in the retrieved aerosol size distributions: both algorithms find the peaks of fine and coarse mode aerosol size to be around ~0.15 µm and ~2 µm, respectively. Due to the lower sensitivity of AirMSPI's longest wavelength 865 nm to coarse mode aerosols, some differences in coarse mode aerosol size can be observed. This indicates the impact of insufficient observational information about certain aerosol properties. Comparisons of pixel-scale AOD, SSA, and size distribution at other retrieval cases show a similar quality of agreement. (namely in the range 0.28 ≤ AOD 445nm ≤ 0.40 with mean value 0.32) across the retrieval area is captured by the correlated multi-pixel inversion. Implementation of our approach using several datasets with even higher (~90%) variation of aerosol loading over several smoke scenes acquired by AirMSPI during the recent Aerosol Characterization from Polarimeter and Lidar (ACEPOL) campaign showed the algorithm to be capable of capturing this variation. The regression in Figure 9b shows correlations and low bias of SSAs ( R 2 > 0.40, a > 0.60, b < 0.030, and MAD ≤ 0.004) from the two inversions as well. Figure 9c shows basic consistency in the retrieved aerosol size distributions: both algorithms find the peaks of fine and coarse mode aerosol size to be around ~0.15μm and ~2μm, respectively. Due to the lower sensitivity of AirMSPI's longest wavelength 865 nm to coarse mode aerosols, some differences in coarse mode aerosol size can be observed. This indicates the impact of insufficient observational information about certain aerosol properties. Comparisons of pixel-scale AOD, SSA, and size distribution at other retrieval cases show a similar quality of agreement.

(a) against the previous retrievals using multi-pixel inversion (MPI, [START_REF] Xu | Coupled retrieval of aerosol properties and land surface reflection using the Airborne Multiangle SpectroPolarimetric Imager[END_REF]). The AirMSPI dataset used for retrieval is the same as in Figure 5 against the previous retrievals using multi-pixel inversion (MPI, [START_REF] Xu | Coupled retrieval of aerosol properties and land surface reflection using the Airborne Multiangle SpectroPolarimetric Imager[END_REF]). The AirMSPI dataset used for retrieval is the same as in Figure 5 over the AERONET Baskin site. The results for seven AirMSPI spectral bands are plotted in different colors: pink (355 nm), purple (380 nm), dark blue (445 nm), light blue (470 nm), green (555 nm), red (660 nm), and brown (865 nm). (b) The same as Figure 9a but for SSA. (c) Comparison of image-mean volume-weighted aerosol size distribution retrieved from correlated multi-pixel inversion and multi-pixel inversion.

Figures 10 and 11 compare AirMSPI and AERONET retrievals of AOD and SSA respectively. Since the temporal variation of aerosol loading and properties is not constant, non-symmetric AERONET error bars can be observed in Figures 10 and11. Under the circumstance of no AERONET reference data before or after AirMSPI measurement, only AERONET measurement/retrieval error is reported. Figures 12 and13 compare real and imaginary parts of aerosol refractive index, respectively. The AERONET spectral aerosol products were linearly interpolated in wavelength to match the AirMSPI band centers. Comparisons of fine and coarse mode effective radii are shown in Figure 14. To facilitate the comparison of aerosol size, an effective radius was calculated for fine and coarse mode aerosols from AirMSPI retrievals using Equation [START_REF] Dubovik | GRASP: A versatile algorithm for characterizing the atmosphere[END_REF] in a previous study [START_REF] Xu | Coupled retrieval of aerosol properties and land surface reflection using the Airborne Multiangle SpectroPolarimetric Imager[END_REF]. For AOD, linear regression is performed. Values of regression related parameters (a, b, R 2 ) and MAD are indicated in all panels. The AOD regression shows a spectral means of coefficient of determination 0.91, slope 0.93, and intercept 0.03, reflecting high retrieval quality. While SSA and refractive index in Figures 111213show relatively larger differences between the AirMSPI and AERONET retrievals, the differences are generally within their respective uncertainties, which in turn depend on AirMSPI and AERONET observation errors and the sensitivities of the respective retrieval algorithms. Figure 14 shows a maximum difference of 30% between AirMSPI and AERONET retrieved fine mode aerosol size, whereas larger differences (up to 80%) are observed in coarse mode aerosol size. As noted above, shortwave infrared spectral bands, which AirMSPI lacks, are necessary to constrain the coarse mode aerosol size.

Remote Sens. 2018, 10, x FOR PEER REVIEW 42 of 49 Figure 14. To facilitate the comparison of aerosol size, an effective radius was calculated for fine and coarse mode aerosols from AirMSPI retrievals using Equation ( 35) in a previous study [START_REF] Xu | Coupled retrieval of aerosol properties and land surface reflection using the Airborne Multiangle SpectroPolarimetric Imager[END_REF]. For AOD, linear regression is performed. Values of regression related parameters (a, b, R 2 ) and MAD are indicated in all panels. The AOD regression shows a spectral means of coefficient of determination 0.91, slope 0.93, and intercept 0.03, reflecting high retrieval quality. While SSA and refractive index in show relatively larger differences between the AirMSPI and AERONET retrievals, the differences are generally within their respective uncertainties, which in turn depend on AirMSPI and AERONET observation errors and the sensitivities of the respective retrieval algorithms. Figure 14 shows a maximum difference of 30% between AirMSPI and AERONET retrieved fine mode aerosol size, whereas larger differences (up to 80%) are observed in coarse mode aerosol size. As noted above, shortwave infrared spectral bands, which AirMSPI lacks, are necessary to constrain the coarse mode aerosol size. Table 6 summarizes MAD in several key aerosol properties from correlated multi-pixel inversion approach (this paper) and from the original multi-pixel inversion approach adapted to AirMSPI [START_REF] Xu | Coupled retrieval of aerosol properties and land surface reflection using the Airborne Multiangle SpectroPolarimetric Imager[END_REF]. These properties include AOD, SSA, real and imaginary parts of refractive index ( n r , n i ), and effective radii of fine and coarse particles ( r eff, fine , r eff, coarse ). The bias is evaluated by taking the mean of the absolute difference between AirMSPI and AERONET retrievals at collocated pixels. The correlated multi-pixel inversion has slightly higher MAD than that of original multi-pixel retrieval, namely by ~0.01, 0.015, 0.001, 0.002, and 0.05 for AOD and SSA in visible, n r , n i , r eff, fine , and r eff, coarse , respectively. The deviation of both correlated and original multi-pixel retrievals from AERONET reference data are mostly within the retrieval uncertainties from our algorithms (see Section 3.3 and Section 4.2. of the previous study [START_REF] Xu | Coupled retrieval of aerosol properties and land surface reflection using the Airborne Multiangle SpectroPolarimetric Imager[END_REF]) and those estimated for AERONET, as observed in Figures 10-14 of this work, as well as Figures 45678in the previous study [START_REF] Xu | Coupled retrieval of aerosol properties and land surface reflection using the Airborne Multiangle SpectroPolarimetric Imager[END_REF].

Note that to enhance the retrieval efficiency, two PCs are adopted to perform the retrieval and achieve the retrieval quality in Table 6. Including more PCs in our retrieval would allow capturing more spatial variations of aerosol properties across the imaged area and improve accuracy, but at a cost of decreased retrieval efficiency. In correlated multi-pixel inversion, the retrieval efficiency is Table 6 summarizes MAD in several key aerosol properties from correlated multi-pixel inversion approach (this paper) and from the original multi-pixel inversion approach adapted to AirMSPI [START_REF] Xu | Coupled retrieval of aerosol properties and land surface reflection using the Airborne Multiangle SpectroPolarimetric Imager[END_REF]. These properties include AOD, SSA, real and imaginary parts of refractive index ( n r , n i ), and effective radii of fine and coarse particles ( r eff, fine , r eff, coarse ). The bias is evaluated by taking the mean of the absolute difference between AirMSPI and AERONET retrievals at collocated pixels. The correlated multi-pixel inversion has slightly higher MAD than that of original multi-pixel retrieval, namely by ~0.01, 0.015, 0.001, 0.002, and 0.05 for AOD and SSA in visible, n r , n i , r eff, fine , and r eff, coarse , respectively. The deviation of both correlated and original multi-pixel retrievals from AERONET reference data are mostly within the retrieval uncertainties from our algorithms (see Section 3.3 and Section 4.2. of the previous study [START_REF] Xu | Coupled retrieval of aerosol properties and land surface reflection using the Airborne Multiangle SpectroPolarimetric Imager[END_REF]) and those estimated for AERONET, as observed in Figures 10-14 of this work, as well as Figures 45678in the previous study [START_REF] Xu | Coupled retrieval of aerosol properties and land surface reflection using the Airborne Multiangle SpectroPolarimetric Imager[END_REF].

Note that to enhance the retrieval efficiency, two PCs are adopted to perform the retrieval and achieve the retrieval quality in Table 6. Including more PCs in our retrieval would allow capturing more spatial variations of aerosol properties across the imaged area and improve accuracy, but at a cost of decreased retrieval efficiency. In correlated multi-pixel inversion, the retrieval efficiency is Table 6 summarizes MAD in several key aerosol properties from correlated multi-pixel inversion approach (this paper) and from the original multi-pixel inversion approach adapted to AirMSPI [START_REF] Xu | Coupled retrieval of aerosol properties and land surface reflection using the Airborne Multiangle SpectroPolarimetric Imager[END_REF]. These properties include AOD, SSA, real and imaginary parts of refractive index (n r , n i ), and effective radii of fine and coarse particles (r eff, fine , r eff, coarse ). The bias is evaluated by taking the mean of the absolute difference between AirMSPI and AERONET retrievals at collocated pixels. The correlated multi-pixel inversion has slightly higher MAD than that of original multi-pixel retrieval, namely by ~0.01, 0.015, 0.001, 0.002, and 0.05 for AOD and SSA in visible, n r , n i , r eff, fine , and r eff, coarse , respectively. The deviation of both correlated and original multi-pixel retrievals from AERONET reference data are mostly within the retrieval uncertainties from our algorithms (see Sections 3.3 and 4.2 of the previous study [START_REF] Xu | Coupled retrieval of aerosol properties and land surface reflection using the Airborne Multiangle SpectroPolarimetric Imager[END_REF]) and those estimated for AERONET, as observed in Figures 10-14 of this work, as well as Figures 45678in the previous study [START_REF] Xu | Coupled retrieval of aerosol properties and land surface reflection using the Airborne Multiangle SpectroPolarimetric Imager[END_REF].

Note that to enhance the retrieval efficiency, two PCs are adopted to perform the retrieval and achieve the retrieval quality in Table 6. Including more PCs in our retrieval would allow capturing more spatial variations of aerosol properties across the imaged area and improve accuracy, but at a cost of decreased retrieval efficiency. In correlated multi-pixel inversion, the retrieval efficiency is gained from two aspects: (a) reduction of parameter space by retrieving PCs of correlated fields; (2) use of PC-based RT model. For a case study with 22 correlated aerosol parameters and 30 uncorrelated (assumed) surface parameters per super-pixel, the correlated multi-pixel saves 50% CPU time by using two PCs. The speed gain can be greater if we further capitalize the correlation in surface parameters, as well as the correlation between aerosol and surface parameters. It is anticipated that correlated multi-pixel inversion will save 90% CPU time if an image of 100 × 100 pixels is retrieved simultaneously and all 52 parameters are correlated with each other. Table 6. Mean absolute bias (MAD) in key aerosol properties from correlated multi-pixel inversion approach using two PCs (this paper) and from original multi-pixel inversion approach adapted to AirMSPI [START_REF] Xu | Coupled retrieval of aerosol properties and land surface reflection using the Airborne Multiangle SpectroPolarimetric Imager[END_REF]. The bias is calculated by taking the mean of the absolute difference between AirMSPI and AERONET retrievals at collocated pixels. A total of 27 AirMSPI datasets are used. Interpolation is used to calculate some AERONET parameters at central wavelengths of the AirMSPI spectral bands. By performing regression against AERONET reference AOD, the coefficients of determination (R 2 ) are given in columns 3 and 5. The ranges of non-AOD parameters are too small and the sample size is limited. Therefore a reliable regression analysis is not established and R 2

Summary and Outlook

Without utilizing correlations among aerosol parameters, optimization-based retrievals are confronted with a high-dimensional parameter space. However, certain types of aerosols or certain combinations of aerosol fields generally prevail within a targeted area, and consequently some aerosol parameters are correlated with each other (in other words, have high linear dependency, as captured by PC analysis). Due to the lack of accurate physical models, however, it is hard to quantify the physical processes and accurately quantify the correlations between all parameters. To mitigate the influence of model assumptions, a priori information about aerosol correlation informed by ground or other types of measurements is helpful. This motivates our development of PC-based aerosol inversion approach to improve the inversion stability and efficiency by reducing the number of retrieval parameters. The algorithm makes use of multiple types of constraints, including across-pixel smoothness constraints transformed to be imposed on the PCs to retrieve zero sum of the PC weights, and orthogonality and unit norm conditions on PC weights and PC vectors, respectively. While applying smoothness constraints to the PCs instead of the individual aerosol parameters, the regularization (smoothness) remains faithful to both aerosol correlations and smooth spatial and spectral variation within the scene. To accelerate the multi-pixel retrieval, a PC-based RT model is developed, which capitalizes on the aerosol correlations and mutual orthogonality of the PC vectors. The retrieval methodology was tested by comparing aerosol retrievals from 27 AirMSPI datasets acquired between 2012 and 2016 with collocated aerosol reference data reported by AERONET. Mean absolute differences between AirMSPI and AERONET retrievals are found to be ~0.029 and 0.038 for AOD and SSA, respectively.

The correlated multi-pixel inversion established in this paper is informed by prior estimates of aerosol properties within the retrieval in order to generate initial guesses for the PC vectors and to calculate the co-variance of a priori. Potential sources of such information can include AERONET climatologies, chemical transport model results, or satellite-based aerosol inversion output obtained, for example, from the MISR operational aerosol retrieval [START_REF] Martonchik | Techniques for the retrieval of aerosol properties over land and ocean using multiangle data[END_REF][START_REF] Martonchik | Retrieval of aerosol properties over land using MISR observations[END_REF], and the near-real time POLDER aerosol retrieval [START_REF] Dubovik | GRASP: A versatile algorithm for characterizing the atmosphere[END_REF]. Though not implemented in this study, it would be interesting to further impose on retrieval the correlations in surface properties and the correlations between aerosol and surface properties. To derive a priori of the correlations in surface properties, PC analysis can be performed over a surface reflectance dataset, such as the one based on the Multi-Angle Implementation of Atmospheric Correction (MAIAC) algorithm [START_REF] Lyapustin | Multiangle implementation of atmospheric correction (MAIAC): 1. Radiative transfer basis and look-up tables[END_REF][START_REF] Lyapustin | Multiangle implementation of atmospheric correction (MAIAC): 2. Aerosol algorithm[END_REF], the ASTER spectral library [START_REF] Baldridge | The ASTER spectral library version 2.0. Remote Sens[END_REF], and the U.S. Geological Survey (USGS) digital spectral library [START_REF] Clark | USGS digital spectral library splib06a[END_REF]. To impose the correlation between aerosol and surface, however, one might encounter a complex surface condition over a targeted area (e.g., with co-existence of multiple types of surfaces such as soil, water, snow, forest, etc.). Under such circumstances, surface pixel classification can be performed to identify surface types first, either from existing surface climatology or from various indices for vegetation, snow, soil, water, etc. based on their spectral difference. Then one can assign surface-type dependent sets of PCs to pixels, each set containing a small number of PCs. Finally, these surface-type dependent sets of PCs are retrieved simultaneously with the imposition of the constraints regarding smooth variations of relevant aerosol/surface properties in spatial, spectral and/or temporal directions. To account for dependent sets of PCs in a state vector, some modifications are necessary in formulating the smoothness constraints. Using such a strategy is expected to gain more retrieval efficiency than if one directly relaxes more PCs in retrievals to capture strong spatial variations of surface properties.

While the correlated multi-pixel approach developed here allows a retrieval of both PC vectors and PC weights, another way to capture regionally limited variability in aerosol type is to use a traditional lookup table (LUT) based aerosol retrieval. This approach is implemented in some operational aerosol retrievals employed by, e.g., Multi-angle Imaging SpectroRadiometer (MISR) [START_REF] Martonchik | Techniques for the retrieval of aerosol properties over land and ocean using multiangle data[END_REF][START_REF] Martonchik | Retrieval of aerosol properties over land using MISR observations[END_REF] and Moderate resolution Imaging Spectroradiometer (MODIS) [START_REF] Levy | Second-generation operational algorithm: Retrieval of aerosol properties over land from inversion of moderate resolution imaging spectroradiometer spectral reflectance[END_REF], and has extremely high retrieval efficiency. With reliable estimates of aerosol properties from other sources as noted above, a "smart" LUT can be generated, which then serves as the basis for a reliable set of PC vectors. If there is high confidence in the representativeness of these PCs, the retrieval could be confined to the PC weights only, which will be faster than the combined inversion of PC vectors and PC weights. Such an approach would compensate for the weakness in traditional LUT approach, namely that aerosol mixtures are confined to a discretized aerosol parameter space. Further testing of these ideas are planned using multi-angle satellite observations from MISR [START_REF] Diner | Multiangle imaging spectro-radiometer (MISR) description and experiment overview[END_REF]. The vector consisting of correlated aerosol properties -calculated from the solution at q-th iteration

x q,surf The vector consisting of uncorrelated surface reflection properties -containing in the solution at q-th iteration x retrieved Retrieved column state vector ∆x syst Systematic error in retrieval x uncorr,p Column vector containing uncorrelated parameters (fields) for p-th pixel x true Column state vector associated with true solution x wv Column vector including PC weights and vectors (∆x) j

The retrieval error in j-th parameter

Appendix B Smoothness Matrix to Constrain Uncorrelated Parameter Retrieval

To explain construction of the smoothness matrix for PC vectors and weights, we start by describing the smoothness matrix used in the original multi-pixel algorithm formulated by Dubovik et al. [START_REF] Dubovik | Statistically optimized inversion algorithm for enhanced retrieval of aerosol properties from spectral multi-angle polarimetric satellite observations[END_REF]. This appendix forms the basis for extension to PC-based smoothness matrix, described in Appendices C and D. Two major classes of constraints are imposed on the PC retrieval: across-pixel (spatial) constraints and within-pixel (e.g., spectral) constraints. The following matrix incorporates both across-pixel and within-pixel constraints for a set of uncorrelated aerosol parameters:

γ uncorr Ω uncorr = γ uncorr,∆ Ω uncorr,∆ +       γ uncorr, Ω uncorr,1 0 • • • 0 0 γ uncorr, Ω uncorr,2 • • • 0 • • • • • • . . . • • • 0 0 • • • γ uncorr, Ω uncorr,N pixel       (A1)
where γ uncorr,∆ Ω uncorr,∆ is a block matrix that includes across-pixel smoothness constraints, γ uncorr, Ω uncorr, ,p is a block matrix that includes within-pixel smoothness constraints over the parameters associated with p-th pixel, and 0 is the zero block matrix. Pixel resolved matrices γ uncorr, Ω uncorr, ,p do not interfere with each other, so they are arranged along the diagonal axis of the large matrix on the right-hand-side of the above equation. The across-pixel constraints that ensure smooth variation of a parameter in space and time is expressed by,

γ uncorr,∆ Ω uncorr,∆ = γ uncorr,u Ω uncorr,u + γ uncorr,v Ω uncorr,v + γ uncorr,t Ω uncorr,t (A2)
To simplify the explanation, consider an image of N u pixels along the horizontal spatial dimension u, two pixels along the horizontal spatial dimension v (N v = 2) and two successive measurements in time (N t = 2) of some parameter x. The construction of a smoothness matrix along the vertical direction is neglected for simplicity. In this case, the single-column state vector is arranged as,

x uncorr =                         x(u 1 ; v 1 ; t 1 )
x(u 2 ; v 1 ; t 1 ) x(u 3 ; v 1 ; t 1 ) . . .

x(u N u ; v 1 ; t 1 ) x(u 1 ; v 2 ; t 1 ) x(u 2 ; v 2 ; t 1 ) x(u 3 ; v 2 ; t 1 ) . . .

x(u N u ; v 2 ; t 1 ) to continue

                                                continued x(u 1 ; v 1 ; t 2 ) x(u 2 ; v 1 ; t 2 ) x(u 3 ; v 1 ; t 2 ) . . . x(u N u ; v 1 ; t 2 ) x(u 1 ; v 2 ; t 2 )
x(u 2 ; v 2 ; t 2 ) x(u 3 ; v 2 ; t 2 ) . . .

x(u N u ; v 2 ; t 2 )

                        (A3)
The corresponding smoothness matrix constraining the horizontal and temporal variation of x is given by,

where

s (muncorr,u) uncorr,u =       d (muncorr,u) 1 (1) d (muncorr,u) 1 (2) • • • d (muncorr,u) 1 (muncorr,u + 1) 0 • • • 0 0 d (muncorr,u) 2 (1) d (muncorr,u) 2 (2) • • • d (muncorr,u) 2 (muncorr,u + 1) • • • 0 • • • • • • • • • • • • • • • • • • • • • 0 0 • • • d (muncorr,u) Nu-muncorr,u (1) d
(muncorr,u) Nu-muncorr,u (2) • • • d (muncorr,u) Nu-muncorr,u (muncorr,u + 1)       (A6)
As an example, the differentiation array for the first order of difference in above equation is,

d (m uncorr,u =1) j = 1 ∆ 1 (j) , -1 ∆ 1 (j) (A7)
where ∆ 1 (j) accounts for the distance between neighboring pixels in the direction u and is evaluated by Equation ( 40), and for the second order of difference,

d (m uncorr,u =2) j = 2 ∆ 1 (j)[∆ 1 (j)+∆ 1 (j+1)] , - 2 ∆ 1 (j)∆ 1 (j+1) , 2 ∆ 1 (j+1)[∆ 1 (j)+∆ 1 (j+1)] (A8)
Moreover, the smoothness matrix in the direction v is,

√ γ uncorr,v S (m uncorr,v ) uncorr,v = √ γ uncorr,v s (m uncorr,v ) uncorr,v 0 0 s (m uncorr,v ) uncorr,v (A9)
where

s (muncorr,v=1) uncorr,v =        d (muncorr,v) 1 (1) 0 • • • 0 d (muncorr,v) 1 (2) 0 • • • 0 0 d (muncorr,v) 1 (1) • • • 0 0 d (muncorr,v) 1 (2) • • • 0 • • • • • • . . . • • • • • • • • • . . . • • • 0 0 • • • d (muncorr,v) 1 (1) 0 0 • • • d (muncorr,v) 1 (2)        (A10)
where we assume m uncorr,v = 1 and the calculation of the differentiation array d accounts for the distance between neighboring pixels in the direction v and is evaluated by Equation ( 40). The smoothness matrices in the direction t is,

s (muncorr,t=1) uncorr,t =       d (muncorr,t) 1 (1) 0 • • • 0 0 0 • • • 0 d (muncorr,t) 1 (2) 0 0 • • • 0 0 • • • 0 0 d (muncorr,t) 1 (1) • • • 0 0 0 • • • 0 0 d (muncorr,t) 1 (2) 0 • • • 0 0 • • • 0 d (muncorr,t) 1 (1) 0 • • • 0 0 0 • • • 0 d (muncorr,t) 1 (2) 0 0 • • • 0 0 • • • 0 0 d (muncorr,t) 1 (1) • • • 0 0 0 • • • 0 0 d (muncorr,t) 1 (2) 0 • • • 0 0 • • • 0       (A11)
where we assume m uncorr,t = 1 and the calculation of the differentiation array d accounts for the temporal gap between successive measurements and is evaluated by Equation [START_REF] Levy | Second-generation operational algorithm: Retrieval of aerosol properties over land from inversion of moderate resolution imaging spectroradiometer spectral reflectance[END_REF]. Note that the above formalism applies to all uncorrelated but smoothly varying parameters. The incorporation of these smoothness matrices for different parameters into an overall matrix γ uncorr,∆ Ω uncorr,∆ is designed to account for the locations of these parameters in a retrieval state vector. (muncorr, (j) + 1)

• • • 0 • • • • • • • • • • • • • • • . . . • • • 0 0 • • • d (muncorr, (j))
TP,uncorr(j)-muncorr, (j) (1) d

(muncorr, (j))

TP,uncorr(j)-muncorr, (j) (2)

• • • d (muncorr, (j))
TP,uncorr(j)-muncorr, (j) (muncorr, (j) + 1)

       (A15)
where we assume j-th parameter belongs to [TP corr (j)]-th type of correlated parameter, the calculation of the differentiation array d accounts for the distance between neighboring wavelength and proceeds in the way as expressed in Equations (A7) and (A8). The weight matrix in Equation (A13) is expressed as, where

W (muncorr, ) =         W (muncorr, ( 1
W (m uncorr, (i)) (i) =         1 ∆ muncorr, (i) (1) 0 • • • 0 0 1 ∆ muncorr, (i) (2) • • • 0 • • • • • • . . . • • • 0 0 • • • 1 ∆ muncorr, (i) (L(i)-m uncorr, (i))         (A17)
where ∆ m (j) is evaluated by Equation [START_REF] Levy | Second-generation operational algorithm: Retrieval of aerosol properties over land from inversion of moderate resolution imaging spectroradiometer spectral reflectance[END_REF].

To generate γ x corr Ω x corr (see Equation (50) for Ω x corr ), the evaluation of smoothness matrix of the correlated fields in the column vector x containing multi-pixel mean fields proceeds in a similar way as formulated for the within-pixel constraint.

Appendix C Smoothness Matrix for Correlated Parameters

As shown in Appendix B, the original multi-pixel inversion algorithm imposes smoothness constraints directly on the aerosol fields. In the correlated multi-pixel inversion, the correlated parameters are not directly retrieved; the retrieval involves the PC vectors, PC weights, and multi-pixel mean fields. Therefore, the smoothness constraints must be imposed on PCs and PC weights. Like the implementation for uncorrelated parameter retrievals, the overall smoothness matrix includes both across-pixel (γ wv corr,∆ Ω wv corr,∆ ) and within-pixel constraints (γ wv corr, Ω wv corr, ), namely,

γ
The explicit forms of γ wv corr,∆ Ω wv corr,∆ and γ wv corr, Ω wv corr, will be given in the Appendices C.1 and C.2.

Appendix C.1 Across-Pixel Smoothness Constraints

To simplify the discussion, the formalism given in this and the next section assumes that the smoothness constraint applied to a parameter x occurs in the dimension u. It can be easily extended to enable construction of smoothness matrices in other spatial dimensions.

To constrain across-pixel variation of a certain parameter, the γ wv corr,∆ Ω wv corr,∆ matrix in Equation (A18) is given by

            
0 0 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 0 • • • • • • 0 C (m
        0 0 • • • 0 • • • • • • • • • • • • v 1 (i) v 2 (i) • • • v N PC (i) • • • • • • . . . • • • 0 0 • • • 0         (A23)
Moreover, Equation ( 59) is implemented to evaluate d(S (m wv corr,∆ ) wv,i x i )/dx i in Equation (A19).

Appendix C.2 Within-Pixel Smoothness Matrix

The overall within-pixel smoothness matrix γ wv corr, Ω wv corr, in Equation (A18) is expressed as, where the column vector x i has the same length as that of the state vector and contains all-pixel PC weights and the i-th type of correlated fields included in all PC vectors (all rest elements of x i are zero), and the weight matrix W that has a similar structure as Equations (A16) and (A17) (but for the i-th type of correlated parameters and is expanded to account for the number of pixels). To account for the smooth variation of a type of correlated parameters, the contributions of PC vectors and weights are coupled in the following way, (m wv corr, + 1)

γ
• • • 0 0 0 0 0 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 0 0 0 • • • • • • • • • • • • d (m wv corr, ) L(j)-m wv corr, (1) d

Appendix D Decoupled Smoothness Constraints

Smooth variations are often directly observed in PC weights and elements in a PC vector associated with a certain type of parameter. Direct imposition of smoothness helps stabilize the first few iterations when a priori information about the PC vectors and weights is insufficient. In an integrated form, the across-pixel and within-PC smoothness matrix is expressed as

γ w/v

Figure 1 .

 1 Figure 1. (a) Top left panel: AERONET inversion of aerosol fields and properties (in natural logarithm space) including total volume concentration ( C v,tot ), fraction of fine mode aerosols ( f fine ),

Figure 1 .

 1 Figure 1. (a) Top left panel: AERONET inversion of aerosol fields and properties (in natural logarithm space) including total volume concentration (C v,tot ), fraction of fine mode aerosols ( f fine ), effective radii of fine (r eff, fine ) and coarse (r eff, coarse ) mode aerosols, column effective real part of refractive indices (n r, 1-4 ) at 0.439, 0.675, 0.870, and 1.018 µm, respectively, imaginary part of refractive indices (n i, 1-4 ), and spherical particle volume concentration (C v,sphere ). Though not specified by use of legend, each color is associated with an independent set of retrieval parameters. Top right panel: spatial and temporal mean of all retrievals and the first four principal components (PCs). Bottom left panel: percentage variance of the aerosol fields captured by PCs, indicating that 41%, 62%, 77%, and 85% of variance is captured by the first one, two, three, and four PCs, respectively. Bottom right panel: Aerosol optical depth (AOD) and single scattering albedo (SSA) at 440 nm that are reported in the AERONET retrievals analyzed here. (b) Regression of the derived aerosol properties from four PC against input AERONET values. The scatter plot is colored by the density of the points. Each panel corresponds to a specific aerosol property indicated in the title.

Figure 2 .

 2 Figure 2. Same as Figure 1 but analysis was performed for a circular domain of 2000 km around AERONET site in Namibe, Angola. The first four PCs capture 56%, 74%, 83%, and 90% variance of the aerosol fields.

Figure 2 .

 2 Figure 2. Same as Figure 1 but analysis was performed for a circular domain of 2000 km around AERONET site in Namibe, Angola. The first four PCs capture 56%, 74%, 83%, and 90% variance of the aerosol fields.

Figure 4 .

 4 Figure 4. Depiction of the coupled atmosphere-surface (CAS) system model. The sun illuminates the top-of-atmosphere with solar zenith angle  0 and azimuthal plane  0 . The sensor views the atmosphere at view zenith angle  v and azimuthal angle  v . A Gaussian vertical distribution profile for aerosols. The Markov chain model is used for computing polarized RT in the atmosphere. It is then coupled with surface reflection using an adding method.

Figure 4 .

 4 Figure 4. Depiction of the coupled atmosphere-surface (CAS) system model. The sun illuminates the top-of-atmosphere with solar zenith angle θ 0 and azimuthal plane φ 0 . The sensor views the atmosphere at view zenith angle θ v and azimuthal angle φ v . A Gaussian vertical distribution profile for aerosols. The Markov chain model is used for computing polarized RT in the atmosphere. It is then coupled with surface reflection using an adding method.

Figure 5 .

 5 Figure 5. The scheme of PC-based forward radiative transfer modeling of remote sensing observations from an airborne or spaceborne sensor. The interpretation of symbols used in the figure can be found in Table A1 of Appendix A.

Figure 5 .

 5 Figure 5. The scheme of PC-based forward radiative transfer modeling of remote sensing observations from an airborne or spaceborne sensor. The interpretation of symbols used in the figure can be found in Table A1 of Appendix A.

Figure 6 .

 6 Figure 6. (a) Pixel-resolved spectral aerosol optical depth (AOD) used in simulating the radiance and degree of linear polarization (DOLP) in (b). The AODs for seven Airborne Multiangle SpectroPolarimetric Imager (AirMSPI) spectral bands are plotted in different colors: pink (355 nm), purple (380 nm), dark blue (445 nm), light blue (470 nm), green (555 nm), red (660 nm), and brown (865 nm). b) The principal component based RT (PC-RT) simulation of radiance using three PCs (upper left panel) and DOLP (lower left panel) at AirMSPI's nine viewing angles within (-66 • , +66 • ) range around nadir, and the relative error of PC based RT simulation of radiance (top right) and absolute error of DOLP (bottom right) as compared to direct RT computation pixel-by-pixel. Errors are estimated by 100 × (Y PC-RT -Y Direct-RT )/Y Direct-RT for Y = radiance and by (Y PC-RT -Y Direct-RT ) for Y = DOLP.The viewing geometries are adopted from one of the scenes acquired by AirMSPI during the ACEPOL field campaign. The unpolarized and polarized surface reflectance is calculated from the retrieved BRDF and pBRDF parameters, respectively.

Figure 7 .

 7 Figure 7. Algorithm flowchart for correlated multi-pixel retrieval of aerosol and land surface reflection properties. The interpretation of symbols used in the figure can be found in Table A1 of Appendix A.

Figure 7 .

 7 Figure 7. Algorithm flowchart for correlated multi-pixel retrieval of aerosol and land surface reflection properties. The interpretation of symbols used in the figure can be found in Table A1 of Appendix A.

Figure 8 .

 8 Figure 8. (a) High resolution AirMSPI nadir imagery of Baskin, Louisiana, acquired on September 9, 2013. The left image is of radiance at 445, 555, and 660 nm. The location of the Baskin AERONET site is marked. The right image displays DOLP in the three polarimetric bands (470, 660, and 865 nm). The yellow box indicates the viewed at all 9 AirMSPI view angles, and where data were used for retrieval algorithm testing. (b) Lower-resolution imagery (~1.0 km) of the retrieval area after pixel aggregation. The left, middle, and right panels are images of BRF at 445, 555, and 660 nm, respectively. (c) Retrieved AOD, SSA, and surface albedo (A) maps at 555 nm in the left, middle, and right panels, respectively. (d) The AirMSPI retrieved AOD, SSA, and volume weighted aerosol size distribution at the pixel closest to the Baskin AERONET site, compared to the AERONET-derived values.

Figure 8 .

 8 Figure 8. (a) High resolution AirMSPI nadir imagery of Baskin, Louisiana, acquired on September 9, 2013. The left image is of radiance at 445, 555, and 660 nm. The location of the Baskin AERONET site is marked. The right image displays DOLP in the three polarimetric bands (470, 660, and 865 nm). The yellow box indicates the area viewed at all 9 AirMSPI view angles, and where data were used for retrieval algorithm testing. (b) Lower-resolution imagery (~1.0 km) of the retrieval area after pixel aggregation. The left, middle, and right panels are images of BRF at 445, 555, and 660 nm, respectively. (c) Retrieved AOD, SSA, and surface albedo (A) maps at 555 nm in the left, middle, and right panels, respectively. (d) The AirMSPI retrieved AOD, SSA, and volume weighted aerosol size distribution at the pixel closest to the Baskin AERONET site, compared to the AERONET-derived values.

Figure 9 .Figure 9 .

 99 Figure 9. Cont.

  Figure 9. (a) Regression of pixel-scale AOD retrieved from correlated multi-pixel inversion (CMPI)against the previous retrievals using multi-pixel inversion (MPI,[START_REF] Xu | Coupled retrieval of aerosol properties and land surface reflection using the Airborne Multiangle SpectroPolarimetric Imager[END_REF]). The AirMSPI dataset used for retrieval is the same as in Figure5over the AERONET Baskin site. The results for seven AirMSPI spectral bands are plotted in different colors: pink (355 nm), purple (380 nm), dark blue (445 nm), light blue (470 nm), green (555 nm), red (660 nm), and brown (865 nm). (b) The same as Figure9abut for SSA. (c) Comparison of image-mean volume-weighted aerosol size distribution retrieved from correlated multi-pixel inversion and multi-pixel inversion.

Figures 10 and 11

 11 Figures 10 and 11 compare AirMSPI and AERONET retrievals of AOD and SSA respectively, whereas Figures 12 and 13 compare real and imaginary parts of aerosol refractive index, respectively. The AERONET spectral aerosol products were linearly interpolated in wavelength to match the AirMSPI band centers. Comparisons of fine and coarse mode effective radii are shown in

Figure 9 .

 9 Figure9. (a) Regression of pixel-scale AOD retrieved from correlated multi-pixel inversion (CMPI) against the previous retrievals using multi-pixel inversion (MPI,[START_REF] Xu | Coupled retrieval of aerosol properties and land surface reflection using the Airborne Multiangle SpectroPolarimetric Imager[END_REF]). The AirMSPI dataset used for retrieval is the same as in Figure5over the AERONET Baskin site. The results for seven AirMSPI spectral bands are plotted in different colors: pink (355 nm), purple (380 nm), dark blue (445 nm), light blue (470 nm), green (555 nm), red (660 nm), and brown (865 nm). (b) The same as Figure9abut for SSA. (c) Comparison of image-mean volume-weighted aerosol size distribution retrieved from correlated multi-pixel inversion and multi-pixel inversion.

Figure 10 .Figure 10 .

 1010 Figure 10. Regression of AirMSPI retrieved aerosol optical depth (AOD) against AERONET measured values. The upper three panels are for 355, 445, and 470 nm, and the lower three panels are for 555, 660, and 865 nm. Linear interpolation is used to obtain AERONET AOD values at the AirMSPI wavelengths. The AERONET uncertainties are from the ~1-hour window around the time of AirMSPI overflight plus measurement uncertainties (0.01), while the AirMSPI uncertainties are the root mean square of the retrieval uncertainties and standard deviation of pixel-resolved AODs over the whole image. Linear regression analysis yields values of slope a, intercept b, coefficient of determination R 2 , and mean absolute difference (MAD). Values of each are indicated in all panels.

Figure 11 .

 11 Figure 11. Comparison of AirMSPI retrieved single scattering albedo (SSA) against AERONET retrievals. The upper left and right panels are for 445 and 555 nm and the lower left and right panels are for 660 and 865 nm. Linear interpolation is used to obtain AERONET SSA values at the AirMSPI wavelengths. The AirMSPI errors are computed from statistics obtained over the whole image plus the errors evaluated using the method in Section 3.3. Values of MAD are shown.

Figure 12 .

 12 Figure 12. Same as Figure 11 but for the real part of aerosol refractive index.

Figure 11 .Figure 11 .

 1111 Figure 11. Comparison of AirMSPI retrieved single scattering albedo (SSA) against AERONET retrievals. The upper left and right panels are for 445 and 555 nm and the lower left and right panels are for 660 and 865 nm. Linear interpolation is used to obtain AERONET SSA values at the AirMSPI wavelengths. The AirMSPI errors are computed from statistics obtained over the whole image plus the errors evaluated using the method in Section 3.3. Values of MAD are shown.

Figure 12 .

 12 Figure 12. Same as Figure 11 but for the real part of aerosol refractive index.

Figure 12 .

 12 Figure 12. Same as Figure 11 but for the real part of aerosol refractive index.

Figure 13 .

 13 Figure 13. Same as Figure 11 but for the imaginary part of aerosol refractive index.

Figure 14 .

 14 Figure 14. Same as Figure 11 but for the effective radii of fine and coarse mode aerosols.

Figure 13 .Figure 13 .

 1313 Figure 13. Same as Figure 11 but for the imaginary part of aerosol refractive index.

Figure 14 .

 14 Figure 14. Same as Figure 11 but for the effective radii of fine and coarse mode aerosols.

Figure 14 .

 14 Figure 14. Same as Figure 11 but for the effective radii of fine and coarse mode aerosols.

Appendix B. 2

 2 Within-Pixel Smoothness MatrixCertain types of parameters, such as the spectral weight of the microfacet model-based pBRDF function and spectral shape parameters in the RPV model-based BRDF function as a function of wavelength, are subjected to inherent smoothness. The smooth variation of such a type of within-pixel parameters is ensured byγ uncorr, Ω uncorr, = , Ω uncorr, ) p2 • • • 0 (γ uncorr, Ω uncorr, ) N pixel C(m uncorr, (N TP,uncorr )) p (N TP,uncorr ) submatrix for j-th type of correlated parameters,

  muncorr, (N TP,uncorr )) (N TP,uncorr )

  wv corr,∆ ) Npixel-m wv corr,∆ (1, i) 0uncorr,N pixel-m wv corr,∆C (m wv corr,∆ ) Npixel-m wv corr,∆ (2, i) 0uncorr,N pixel-m wv corr,∆ +1 • • • C (m wv corr,∆ ) Npixel-m wv corr,∆ (m wv corr,∆ + 1, i) 0uncorr,N pixelwhere the first column of zero matrices accommodate the multi-pixel mean of correlated parameters, and the explicit form of B matrix is expressed as, corr,∆ (i))p (i, w ... (1)) B (m wv corr,∆ (i)) p (i, w ... (2)) • • • B (m wv corr,∆ (i)) p (i,w ... (N PC )) w...(k)) = γ wv corr,∆ (i) the i-th correlated parameter. The matrix C is a null matrix except its i-th row has fill-in values, namely C (m wv corr,∆ ) p (m , i) = γ wv corr,∆ (i)d (m wv corr,∆ (i)) p (m )

  γ wv corr, (N TP,corr )S (m wv corr, (N TP,corr ),k) v (N TP,corr )

  corr, (NTP,corr),k=1) wp (NTP,corr γ wv corr, (NTP,corr)S (m wv corr, (NTP,corr),k=2) wp (NTP corr ) • • • γ wv corr, (NTP,corr)S (m wv corr, (NTP,corr),k=NPC) wp(NTP,corr) vector for j-th type of correlated parameter described by,

  (59) is implemented to evaluate d(S (m wv corr,∆ ) wv x i )/dx i in Equation (A24).

Table 1 .

 1 Correlation coefficient of AERONET retrieved aerosol properties analyzed in Figure1(for 2000 km circular around Fresno, California). The correlation is calculated for all parameters in logarithmic space.

		C v,tot	f fine	r eff,fine r eff,coarse n r,1	n r,2	n r,3	n r,4	n i,1	n i,2	n i,3	n i,4	C v,sphere
	C v,tot	1.00	-0.34 -0.41	0.08	-0.12 -0.01 0.01	0.02	-0.01 -0.16 -0.23 -0.25	0.06
	f fine	-	1.00	0.18	0.30	-0.34 -0.36 -0.34 -0.35 0.22	0.32	0.31	0.30	0.67
	r eff,fine	-	-	1.00	-0.01 -0.17 -0.35 -0.41 -0.43 -0.16 -0.08 -0.01 -0.01	-0.06
	r eff,coarse	-	-	-	1.00	-0.01 -0.01 -0.01 -0.01 0.25	0.23	0.21	0.20	0.25
	n r,1	-	-	-	-	1.00	0.94	0.88	0.84	0.33	0.31	0.28	0.28	-0.29
	n r,2	-	-	-	-	-	1.00	0.98	0.96	0.28	0.22	0.18	0.17	-0.28
	n r,3	-	-	-	-	-	-	1.00	0.99	0.24	0.16	0.12	0.11	-0.26
	n r,4	-	-	-	-	-	-	-	1.00	0.21	0.12	0.08	0.08	-0.28
	n i,1	-	-	-	-	-	-	-	-	1.00	0.95	0.90	0.87	0.23
	n i,2	-	-	-	-	-	-	-	-	-	1.00	0.98	0.96	0.29
	n i,3	-	-	-	-	-	-	-	-	-	-	1.00	0.99	0.25
	n i,4	-	-	-	-	-	-	-	-	-	-	-	1.00	0.23
	C v,sphere	-	-	-	-	-	-	-	-	-	-	-	-	1.00

Table 2 .

 2 Correlation coefficient of AERONET retrieved aerosol properties analyzed in Figure2(for 2000 km circular domain around Namibe, Angola).

		C v,tot	f fine	r eff,fine r eff,coarse n r,1	n r,2	n r,3	n r,4	n i,1	n i,2	n i,3	n i,4	C v,sphere
	C v,tot	1.00	-0.32 -0.37	-0.05 -0.43 -0.41 -0.42 -0.40 -0.52 -0.54 -0.54 -0.54	0.39
	f fine	-	1.00	0.23	0.42	-0.12 -0.06 0.02	0.09	0.62	0.59	0.57	0.56	0.45
	r eff,fine	-	-	1.00	-0.07	0.20	0.18	0.17	0.14	0.26	0.28	0.27	0.26	-0.01
	r eff,coarse	-	-	-	1.00	0.14	0.19	0.25	0.33	0.36	0.41	0.41	0.42	0.26
	n r,1	-	-	-	-	1.00	0.97	0.92	0.86	0.51	0.58	0.58	0.58	-0.31
	n r,2	-	-	-	-	-	1.00	0.98	0.94	0.57	0.61	0.61	0.61	-0.28
	n r,3	-	-	-	-	-	-	1.00	0.98	0.61	0.63	0.63	0.62	-0.22
	n r,4	-	-	-	-	-	-	-	1.00	0.62	0.63	0.63	0.62	-0.17
	n i,1	-	-	-	-	-	-	-	-	1.00	0.96	0.94	0.92	0.19
	n i,2	-	-	-	-	-	-	-	-	-	1.00	0.99	0.99	0.12
	n i,3	-	-	-	-	-	-	-	-	-	-	1.00	1.00	0.09
	n i,4	-	-	-	-	-	-	-	-	-	-	-	1.00	0.08
	C v,sphere	-	-	-	-	-	-	-	-	-	-	-	-	1.00

  Remote Sens. 2018, 10, x FOR PEER REVIEW 40 of 49 absolute difference (MAD) in AERONET and AirMSPI results is also calculated to measure the overall deviation. Strong correlation and low bias (with a ~ 0.90, b ≤ 0.05, R 2 ≥ 0.88 and MAD < 0.02) are observed. It can also be observed that a variation of aerosol loading by ~30% around the mean

  is reported only for AOD in the table.

	Parameters	Correlated Multi-pixel Inversion	Original Multi-pixel Inversion
		MAD	R 2	MAD	R 2
	AOD 355 nm	0.058	0.917	0.040	0.925
	AOD 445 nm	0.035	0.927	0.024	0.955
	AOD 470 nm	0.030	0.933	0.021	0.960
	AOD 555 nm	0.020	0.933	0.016	0.960
	AOD 660 nm	0.016	0.922	0.013	0.959
	AOD 865 nm	0.015	0.815	0.014	0.851
	SSA 445 nm	0.035	-	0.030	-
	SSA 555 nm	0.036	-	0.030	-
	SSA 660 nm	0.040	-	0.032	-
	SSA 865 nm	0.041	-	0.035	-
	n r,445 nm	0.052	-	0.039	-
	n r,555 nm	0.051	-	0.039	-
	n r,660 nm	0.046	-	0.036	-
	n r,865 nm	0.038	-	0.037	-
	n i,445 nm	0.004	-	0.004	-
	n i,555 nm	0.004	-	0.004	-
	n i,660 nm	0.005	-	0.004	-
	n i,865 nm	0.005	-	0.005	-
	r eff,fine (µm)	0.024	-	0.022	-
	r eff,coarse (µm)	1.050	-	0.993	-

Table A1 .

 A1 Cont.

	Symbol/Abbreviation	Description
	N corr	Number of correlated parameters (fields)
	N f	Total number of observations (all pixels are accounted)
	N i	Total number of i-th type of constraint
	N PC	Total number of principle components
	N pixel	Total number of pixels
	N SC	Number of aerosol size components
	N TP,uncorr	Total number of types of uncorrelated parameters
	O	Constraining matrix that reflects zero-sum of PC weights
	pBRDF	Polarized BRDF
	P	Probability distribution function (PDF)
	PC	Principal component
	PCA	Principal component analysis
	POLDER	Polarization and Directionality of Earth's Reflectance
	q	Ratio of Stokes components Q and I
	q	Iterative step during optimization
	qBRF TOA	Top-of-atmosphere BRF associated with Stokes vector component Q
	Q	Second Stokes vector component
	r	Radius of aerosol
	r eff,coarse	Effective radius of coarse mode aerosols
	r eff,fine	Effective radius of fine mode aerosols
	R	Reflection matrix
	R 2	Coefficient of determination
	R atmos R * atmos	Reflection matrix for atmosphere associated with light illumination from top of the atmosphere Reflection matrix for atmosphere associated with light illumination from bottom of the atmosphere
	R CAS	Reflection matrix for the coupled atmosphere-surface system (CAS)
	RPV	Rahman-Pinty-Verstraete (surface BRDF model)
	R surf	Surface reflection matrix
	R surf, BRDF	Depolarizing part of surface reflection matrix
	R surf, pBRDF	Polarizing part of surface reflection matrix
	RT	Radiative transfer
	s a	Standard deviation of aerosol vertical profile (constrained by Gaussian profile)
	S i,m	Differentiation matrix of m-th order for i-th type of constraint
	SSA	Single scattering albedo
	t	Temporal coordinate
	T	Transmission matrix
	T atmos	Transmission matrix for atmosphere associated with light illumination from top of the atmosphere
	T * atmos	Transmission matrix for atmosphere associated with light illumination from bottom of the atmosphere
	TOA	Top-of-atmosphere
	u	Cartesian coordinate in the direction u
	u	Ratio of Stokes components U and I
	uBRF TOA	Top-of-atmosphere BRF associated with Stokes vector component U
	U	Third Stokes vector component
	U	Constraining matrix that reflects the unit length of a PC
	USGS	U.S. Geological Survey
	v	Cartesian coordinate in the direction v
	v	PC matrix containing N PC columns PC vectors
	vBRF TOA	Top-of-atmosphere BRF associated with Stokes vector component V
	v k	The k-th PC vector
	v state	Column vector containing all PC vectors
	w	PC weight matrix containing N pixel column vectors containing PC weights
	w i	Weight matrix for i-th type of constraint
	w p	Column vector containing PC weights for p-th pixel
	w state	Column vector containing all PC weights
	x	Column state vector including all retrieval parameters
	x a priori	a priori of state vector
	x corr	Column vector containing spatial/temporal mean of correlated parameters (fields)
	x corr,p	Column vector containing correlated parameters (fields) for p-th pixel
	x q,aer

Table A1 .

 A1 Cont.

	Symbol/Abbreviation	Description
	y i Y HS	i-th observational signal Output of RT calculation with high stream approximation
	Y LS	Output of RT calculation with low stream approximation
	z	Variable of a smooth function
	z min	Lower bound of z
	z max	Upper bound of z
	Z	Number of observations per pixel
	δ	Kronecker delta
	δ s	Scale factor that perturbs a PC vector
	ε rand	Random error in measurements
	ε λ ε 2 i ε 2 c ε 2 f θ 0	Spectral weight of pBRDF First diagonal element of C i User-specified threshold value to diagnose the convergence of optimization Expected variance due to measurement errors Solar zenith angle
	θ v	View zenith angle
	λ	Wavelength
	µ 0	Cosine of solar zenith angle
	γ i	Lagrange factor for i-th type of constraint
	φ 0	Solar azimuthal angle
	Ψ i	Objective cost function for i-th type of constraint
	Ψ total	Overall objective cost function
	∇Ψ i	Gradient of the objective cost function for i-th type of constraint
	∇Ψ total	Gradient of the overall objective cost function
	Ω i Ω Ra ...	The smoothness matrix associated with i-th type of constraints The rearranged smoothness matrix from Ω ...
	Ω uncorr	The smoothness matrix imposed on uncorrelated parameters (fields)
	Ω x corr Ω w corr Ω v corr Ω wv corr	The smoothness matrix imposed on spatial/temporal mean of correlated parameters (fields) Smoothness matrix imposed on pixel resolved PC weights Smoothness matrix imposed on a PC vector Smoothness matrix imposed on correlated parameters (fields)
	σ spatio-temporal, x	Standard deviation of a correlated field x
	σ e, x	Uncertainty estimate of a correlated field x
	σ s	Slope variance of polarized BRDF
	τ aer,tot	Total aerosol optical depth
	τ aer,abs	Total absorbing aerosol optical depth
	τ atmos	Atmospheric optical depth
	Γ	Constraining matrix that reflects the mutual orthogonality in PCs

  (u 1 ; v i ; t j ); x(u 2 ; v i ; t j ); • • • ; x(u N u ; v i ; t j ) , and x(u k ; v i ; t j ) is a vector that contains uncorrelated parameters for the pixel (u k , v i ) observed at temporal point t j . Evaluations of γ uncorr,∆ Ω uncorr,∆ and γ uncorr, Ω uncorr, ,p are discussed in Appendices B.1 and B.2, respectively.

	Appendix B.1 Across-Pixel Smoothness Matrix

To facilitate the use of γ uncorr Ω uncorr calculated via Equation (A1), the uncorrelated parameters are grouped together in the order of x uncorr = [x(t 1 ); x(t 2 ); • • • ; x(t N t )], where x(t j ) = x(v 1 ; t j ); x(v 2 ; t j ); • • • ; x(v N v ; t j ) , x(v i ; t j ) =

x

  where the column vector x i has the same length as that of the state vector and contains all-pixel PC weights and the i-th elements of all PC vectors (all rest elements of x i are zero), and γ wv corr,∆ S two matrix components B and C that account for the contributions by PCs and PC weights, respectively, in the following form,

					γ wv corr,∆ Ω wv corr,∆,1 = γ wv corr,∆ Ω wv corr,∆,2 =	N corr ∑ i=1 N corr ∑ i=1	γ wv corr,∆ γ wv corr,∆	d(S d(S	(m wv corr,∆ ) wv,i dx i (m wv corr,∆ ) wv,i dx i	x i ) x i )	T T	d(S wv,i (m wv corr,∆ ) wv,i dx i corr,∆ S γ wv corr,∆ γ wv corr,∆ ) (m wv	x i )	(A19)
											(m wv corr,∆ )
											wv,i
								γ wv corr,∆ S (m wv corr,∆ ) wv,i	=
	       	0 B (m wv w, Npixel-m wv corr,∆ corr,∆ ) • • • • • • 0 B (m wv corr,∆ ) w, p1 (i) 0 corr,∆ ) w, p2 (i) B (m wv	(i)	C (m wv corr,∆ ) p1	(1, i) 0uncorr,p1 C (m wv corr,∆ ) p1	(2, i) 0uncorr,p2

includes

  v corr represent the matrix form of smoothness constraints imposed on PC weights and PC vectors, respectively. The explicit forms of γ w corr Ω w corr and γ v corr Ω v corr are given in the Appendices D.1 and D.2, respectively.

	Remote Sens. 2019, 11, 746	48 of 51
	where γ w corr Ω w corr and γ v corr Ω
	corr Ω w/v corr = γ w corr Ω w corr + γ v corr Ω v corr	(A32)

Remote Sens. 2018, 10, x FOR PEER REVIEW

Acknowledgments:

We are grateful to Carol J. Bruegge, Michael J. Garay, Gerard van Harten, Veljko M. Jovanovic, Olga V. Kalashnikova, Brian E. Rheingans, Felix C. Seidel, Irina N. Tkatcheva, and Mika Tosca of JPL for their support of AirMSPI data acquisition and calibration in multiple NASA field campaigns. We also thank Xu Liu at NASA Langley Research Center and Vijay Natraj at JPL for helpful discussions on principal component analysis to accelerate radiative transfer modeling and atmospheric remote sensing.

the strength of multiple types of constraints, the initial values of Lagrange multipliers for the two PCs used in our retrieval are provided in Tables 345. Note that the difference between Tables 3 and4 and Table 5 is that Tables 3 and4 are for constructing the smoothness constraints over the correlated fields constructed from the combined set of PCs and PC weights. Table 5 is for constructing the smoothness constraints over PC weights and PC vectors separately. Tables 3-5 also list the first guesses of the relevant state vector components and the order of difference for imposing the smoothness constraints on these components. Table 3. Initial guess of image-effective (multi-pixel mean) aerosol parameters and uncorrelated pixel-resolved surface parameters, and the order of difference and Lagrange multipliers for imposing within-pixel smoothness constraints.

Range

First Guess

Order of Finite

Conflicts of Interest:

The authors declare no conflict of interest.

Appendix A Symbols and Abbreviations

Table A1. List of symbols and abbreviations.

Symbol/Abbreviation Description

Appendix D.1 Across-Pixel Smoothness Constraints on PC Weights

The smoothness constraints imposed directly on PC weights is expressed as

where

where the first two columns of zero matrices are to accommodate the multi-pixel mean of correlated parameters and PC vectors, and C(i, k) is calculated by

where the delta function δ(k, k ) = 1 for k = k and 0 for k = k , respectively, and the differentiation arrays are expressed similarly as Equations (A7) and (A8), except that an extra subscript "k" is introduced to the differentiation arrays d to allow different strength and orders of difference of smoothness constraints to impose on a specific (k-th) principal component.

Appendix D.2 Within-PC Smoothness Constraints

The matrix form of within-PC smothness constraints imposed directly upon a PC is evaluated by,

where

where the first zero matrix accommodates the multi-pixel mean of correlated parameters, the last two zero matrices accommodate the PC weights and uncorrelated parameters, and

where where d depends on the type of correlated parameters and is expressed in Equations (A7) and (A8). It can also vary as a function of a specific principal component. The weights matrix W (m v corr,k ) in Equation (A36) is evaluated in the same way as for within-pixel constraints for uncorrelated parameters (see Equation (A16)).