A Correlated Multi-Pixel Inversion Approach for Aerosol Remote Sensing

Abstract

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.

1. 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 [1,2,3,4,5]), 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 [6]. 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.

Figure 1 and Figure 2 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. Correlation coefficient of AERONET retrieved aerosol properties analyzed in Figure 1 (for 2000 km circular around Fresno, California). The correlation is calculated for all parameters in logarithmic space.

	Cv,tot	ffine	reff,fine	reff,coarse	nr,1	nr,2	nr,3	nr,4	ni,1	ni,2	ni,3	ni,4	Cv,sphere
Cv,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
ffine	-	1.00	0.18	0.30	−0.34	−0.36	−0.34	−0.35	0.22	0.32	0.31	0.30	0.67
reff,fine	-	-	1.00	−0.01	−0.17	−0.35	−0.41	−0.43	−0.16	-0.08	−0.01	−0.01	−0.06
reff,coarse	-	-	-	1.00	−0.01	−0.01	−0.01	−0.01	0.25	0.23	0.21	0.20	0.25
nr,1	-	-	-	-	1.00	0.94	0.88	0.84	0.33	0.31	0.28	0.28	−0.29
nr,2	-	-	-	-	-	1.00	0.98	0.96	0.28	0.22	0.18	0.17	−0.28
nr,3	-	-	-	-	-	-	1.00	0.99	0.24	0.16	0.12	0.11	−0.26
nr,4	-	-	-	-	-	-	-	1.00	0.21	0.12	0.08	0.08	−0.28
ni,1	-	-	-	-	-	-	-	-	1.00	0.95	0.90	0.87	0.23
ni,2	-	-	-	-	-	-	-	-	-	1.00	0.98	0.96	0.29
ni,3	-	-	-	-	-	-	-	-	-	-	1.00	0.99	0.25
ni,4	-	-	-	-	-	-	-	-	-	-	-	1.00	0.23
Cv,sphere	-	-	-	-	-	-	-	-	-	-	-	-	1.00

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. Correlation coefficient of AERONET retrieved aerosol properties analyzed in Figure 2 (for 2000 km circular domain around Namibe, Angola).

	Cv,tot	ffine	reff,fine	reff,coarse	nr,1	nr,2	nr,3	nr,4	ni,1	ni,2	ni,3	ni,4	Cv,sphere
Cv,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
ffine	-	1.00	0.23	0.42	−0.12	−0.06	0.02	0.09	0.62	0.59	0.57	0.56	0.45
reff,fine	-	-	1.00	−0.07	0.20	0.18	0.17	0.14	0.26	0.28	0.27	0.26	−0.01
reff,coarse	-	-	-	1.00	0.14	0.19	0.25	0.33	0.36	0.41	0.41	0.42	0.26
nr,1	-	-	-	-	1.00	0.97	0.92	0.86	0.51	0.58	0.58	0.58	−0.31
nr,2	-	-	-	-	-	1.00	0.98	0.94	0.57	0.61	0.61	0.61	−0.28
nr,3	-	-	-	-	-	-	1.00	0.98	0.61	0.63	0.63	0.62	−0.22
nr,4	-	-	-	-	-	-	-	1.00	0.62	0.63	0.63	0.62	−0.17
ni,1	-	-	-	-	-	-	-	-	1.00	0.96	0.94	0.92	0.19
ni,2	-	-	-	-	-	-	-	-	-	1.00	0.99	0.99	0.12
ni,3	-	-	-	-	-	-	-	-	-	-	1.00	1.00	0.09
ni,4	-	-	-	-	-	-	-	-	-	-	-	1.00	0.08
Cv,sphere	-	-	-	-	-	-	-	-	-	-	-	-	1.00

. Moreover, comparison of Table 1 and Table 2 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.

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 [7,8,9]. 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 [10,11,12]. 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 [13,14]. The similarity in the angular shapes of surface bidirectional reflectance factors among MISR’s four spectral bands in the visible and near-infrared [15] 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 [6]), while “correlated” captures the extra correlation constraint. In contrast to most optimization approaches that retrieve the individual aerosol fields (see reviews in [16,17]) 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 [18,19], 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 [20] 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.

2. 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_{\mathrm{corr}}\times N_{\mathrm{pixel}}$ to $(N_{\mathrm{PC}}+1)\times N_{\mathrm{corr}}+N_{\mathrm{PC}}\times N_{\mathrm{pixel}}$, where $N_{\mathrm{pixel}}$, $N_{\mathrm{corr}}$, and $N_{\mathrm{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. List of symbols and abbreviations.

Symbol/Abbreviation	Description
a	Slope of a regression line
$\overline{a}$	Mean slope of a set of regression lines
$a_{\lambda }$	Spectral weight of surface BRDF
A	Coefficient for a basic smooth function shape, e.g., “g(z) = Az 2 + Bz + C” for a parabola
A	Fischer matrix
AERONET	Aerosol Robotic Network
AirMSPI	Airborne Multiangle SpectroPolarimetric Imager
AOD	Aerosol optical depth
ASTER	Advanced Spaceborne Thermal Emission and Reflection Radiometer
Asurf	Surface albedo
b	Intercept of a regression line
$\overline{b}$	Mean intercept of a set of regression lines
B	Coefficient for a basic smooth function shape, e.g., “g(z) = Az 2 + Bz + C” for a parabola
BRDF	Bidirectional reflectance distribution function
BRF	Bidirectional reflectance factor (BRF) associated with Stokes vector component I
c	Constant column vector that weight correlated fields in PC analysis
C	Coefficient for a basic smooth function shape, e.g., “g(z) = Az 2 + Bz + C” for a parabola
${\mathbf{C}}_i$	Covariance matrix of i-th type of constraint
${\mathbf{C}}_{\Delta \mathbf{x},\mathrm{rand}}$	Covariance matrix of random errors in the measurements
Cv,sphere	Column volume concentration of spherical aerosols
Cv,tot	Total column volume concentration of aerosols
d	Differentiation array
$d_{\mathrm{ES}}$	Earth-Sun distance
DOLP	Degree of linear polarization
dV/dln(r)	Volume weighted aerosol size distribution
$E_0$	Exo-atmospheric solar irradiance
${\mathbf{f}}_i^{\ast }$	Column vector of i-th type of constraint
${\mathbf{f}}_i$	Column vector that contains model prediction to fit i-th type of constraint
$\Delta {\mathbf{f}}_i^{\ast }$	Column vector that contains the errors of i-th type of constraint
ffine	Fine mode aerosol fraction
fv,sphere	Volume fraction of spherical aerosols
g(z)	Smooth function with variable z
$g_{\lambda }$	Spectral anisotropy parameter of surface BRDF
GRASP	Generalized Retrieval of Aerosol and Surface ProperEes
h	Cartesian coordinate in the direction h
ha	Central height of aerosol vertical profile (constrained by Gaussian profile)
I	First Stokes vector component
I	Identity matrix
Imeas	Measured radiance
$k_{\lambda }$	Spectral anisotropy parameter of surface BRDF
$k_{\gamma }$	Shadowing width of polarized BRDF
$\mathbf{K}$	Jacobian matrix
${\mathbf{K}}_i$	The Jacobian matrix associated with i-th type of constraint
L(j)	Length of j-th type of correlated parameters (fields)
m	Order of difference in constructing smoothness matrix
M	Number of types of constraints imposed on retrieval
MAD	Mean absolute difference
MAIAC	MultiAngle Implementation of Atmospheric Correction
MISR	Multi-angle Imaging SpectroRadiometer
MODIS	Moderate Resolution Imaging Spectroradiometer
nr, j	Real part of refractive index at j-th wavelength
ni, j	Imaginary part refractive index at j-th wavelength
$N_{\mathrm{a}}$	Total number of retrieval parameters
$N_{{\mathrm{a}}^{\ast }}$	Total number of a priori estimate of parameters
$N_{\mathrm{c}}$	Total number of constraints imposed on retrieval
$N_{\mathrm{corr}}$	Number of correlated parameters (fields)
$N_{\mathrm{f}}$	Total number of observations (all pixels are accounted)
$N_i$	Total number of i-th type of constraint
$N_{\mathrm{PC}}$	Total number of principle components
$N_{\mathrm{pixel}}$	Total number of pixels
$N_{\mathrm{SC}}$	Number of aerosol size components
$N_{\mathrm{TP},\mathrm{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
qBRFTOA	Top-of-atmosphere BRF associated with Stokes vector component Q
Q	Second Stokes vector component
r	Radius of aerosol
reff,coarse	Effective radius of coarse mode aerosols
reff,fine	Effective radius of fine mode aerosols
R	Reflection matrix
R2	Coefficient of determination
${\mathbf{R}}_{\mathrm{atmos}}$	Reflection matrix for atmosphere associated with light illumination from top of the atmosphere
${\mathbf{R}}_{\mathrm{atmos}}^{\ast }$	Reflection matrix for atmosphere associated with light illumination from bottom of the atmosphere
${\mathbf{R}}_{\mathrm{CAS}}$	Reflection matrix for the coupled atmosphere-surface system (CAS)
RPV	Rahman-Pinty-Verstraete (surface BRDF model)
${\mathbf{R}}_{\mathrm{surf}}$	Surface reflection matrix
${\mathbf{R}}_{\mathrm{surf}, \mathrm{BRDF}}$	Depolarizing part of surface reflection matrix
${\mathbf{R}}_{\mathrm{surf}, \mathrm{pBRDF}}$	Polarizing part of surface reflection matrix
RT	Radiative transfer
sa	Standard deviation of aerosol vertical profile (constrained by Gaussian profile)
${\mathbf{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
${\mathbf{T}}_{\mathrm{atmos}}$	Transmission matrix for atmosphere associated with light illumination from top of the atmosphere
${\mathbf{T}}_{\mathrm{atmos}}^{\ast }$	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
uBRFTOA	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
$\mathbf{v}$	PC matrix containing NPC columns PC vectors
vBRFTOA	Top-of-atmosphere BRF associated with Stokes vector component V
${\mathbf{v}}_k$	The k-th PC vector
${\mathbf{v}}_{\mathrm{state}}$	Column vector containing all PC vectors
w	PC weight matrix containing Npixel column vectors containing PC weights
${\mathbf{w}}_i$	Weight matrix for i-th type of constraint
${\mathbf{w}}_p$	Column vector containing PC weights for p-th pixel
${\mathbf{w}}_{\mathrm{state}}$	Column vector containing all PC weights
$\mathbf{x}$	Column state vector including all retrieval parameters
${\mathbf{x}}^{a priori}$	a priori of state vector
${\overline{\mathbf{x}}}_{\mathrm{corr}}$	Column vector containing spatial and temporal mean of correlated parameters (fields)
${\mathbf{x}}_{\mathrm{corr},p}$	Column vector containing correlated parameters (fields) for p-th pixel
${\mathbf{x}}_{q,\mathrm{aer}}$	The vector consisting of correlated aerosol properties – calculated from the solution at q-th iteration
${\mathbf{x}}_{q,\mathrm{surf}}$	The vector consisting of uncorrelated surface reflection properties – containing in the solution at q-th iteration
${\mathbf{x}}^{\mathrm{retrieved}}$	Retrieved column state vector
$\Delta {\mathbf{x}}_{\mathrm{syst}}$	Systematic error in retrieval
${\mathbf{x}}_{\mathrm{uncorr},p}$	Column vector containing uncorrelated parameters (fields) for p-th pixel
${\mathbf{x}}^{\mathrm{true}}$	Column state vector associated with true solution
${\mathbf{x}}_{\mathrm{wv}}$	Column vector including PC weights and vectors
${(\Delta x)}_j$	The retrieval error in j-th parameter
$y_i$	i-th observational signal
${\mathbf{Y}}^{\mathrm{HS}}$	Output of RT calculation with high stream approximation
${\mathbf{Y}}^{\mathrm{LS}}$	Output of RT calculation with low stream approximation
z	Variable of a smooth function
zmin	Lower bound of z
zmax	Upper bound of z
Z	Number of observations per pixel
$\delta$	Kronecker delta
${\delta }_{\mathrm{s}}$	Scale factor that perturbs a PC vector
${\epsilon }_{\mathrm{rand}}$	Random error in measurements
${\epsilon }_{\lambda }$	Spectral weight of pBRDF
${\epsilon }_i^2$	First diagonal element of Ci
${\epsilon }_{\mathrm{c}}^2$	User-specified threshold value to diagnose the convergence of optimization
${\epsilon }_{\mathrm{f}}^2$	Expected variance due to measurement errors
${\theta }_0$	Solar zenith angle
${\theta }_{\mathrm{v}}$	View zenith angle
$\lambda$	Wavelength
${\mu }_0$	Cosine of solar zenith angle
${\gamma }_i$	Lagrange factor for i-th type of constraint
${\varphi }_0$	Solar azimuthal angle
${\Psi }_i$	Objective cost function for i-th type of constraint
${\Psi }_{\mathrm{total}}$	Overall objective cost function
$\nabla {\Psi }_i$	Gradient of the objective cost function for i-th type of constraint
$\nabla {\Psi }_{\mathrm{total}}$	Gradient of the overall objective cost function
${\mathbf{\Omega }}_i$	The smoothness matrix associated with i-th type of constraints
${\mathbf{\Omega }}_{\dots }^{\mathrm{Ra}}$	The rearranged smoothness matrix from ${\mathbf{\Omega }}_{\dots }$
${\mathbf{\Omega }}_{\mathrm{uncorr}}$	The smoothness matrix imposed on uncorrelated parameters (fields)
${\mathbf{\Omega }}_{\mathrm{corr}}^{{\displaystyle \overline{\mathrm{x}}}}$	The smoothness matrix imposed on spatial and temporal mean mean of correlated parameters (fields)
${\mathbf{\Omega }}_{\mathrm{corr}}^{\mathrm{w}}$	Smoothness matrix imposed on pixel resolved PC weights
${\mathbf{\Omega }}_{\mathrm{corr}}^{\mathrm{v}}$	Smoothness matrix imposed on a PC vector
${\mathbf{\Omega }}_{\mathrm{corr}}^{\mathrm{wv}}$	Smoothness matrix imposed on correlated parameters (fields)
${\sigma }_{\mathrm{spatio}-\mathrm{temporal}, x}$	Standard deviation of a correlated field x
${\sigma }_{\mathrm{e}, x}$	Uncertainty estimate of a correlated field x
${\sigma }_{\mathrm{s}}$	Slope variance of polarized BRDF
${\tau }_{\mathrm{aer},\mathrm{tot}}$	Total aerosol optical depth
${\tau }_{\mathrm{aer},\mathrm{abs}}$	Total absorption aerosol optical depth
${\tau }_{\mathrm{atmos}}$	Atmospheric optical depth
$\mathbf{\Gamma }$	Constraining matrix that reflects the mutual orthogonality in PCs

). 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 Section 3.2.2, Section 3.2.3, Section 3.2.4, Section 3.2.5, Section 3.2.6 and Section 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 [6] 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.

To preserve generality of the approach, we use the symbol ${\mathbf{x}}_{\mathrm{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 [21,22,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 (${\mathbf{x}}_{\mathrm{corr},p}$) are constructed from pixel-resolved PC weights in vector form (${\mathbf{w}}_p$), spatially and temporally effective PC matrix (v), and spatial and temporal mean vector ${\overline{\mathbf{x}}}_{\mathrm{corr}}$, namely

$${\mathbf{x}}_{\mathrm{corr},p}={\overline{\mathbf{x}}}_{\mathrm{corr}}+\mathbf{v}\times {\mathbf{w}}_p$$

(1)

where ${\overline{\mathbf{x}}}_{\mathrm{corr}}$ is a single column vector and the matrix v is made up of $N_{\mathrm{PC}}$ column vectors consisting of correlated aerosol fields, with $N_{\mathrm{PC}}$ denoting the number of retrieved principal components. By integrating all PCs into a matrix, we have

$$\mathbf{v}=\left[\begin{array}{cccc}{\mathbf{v}}_1& {\mathbf{v}}_2& \cdots & {\mathbf{v}}_{N_{\mathrm{PC}}}\end{array}\right]$$

(2)

where each PC is composed of $N_{\mathrm{corr}}$ elements, i.e.,

$${\mathbf{v}}_k=\left[v_k(1);v_k(2);\dots ;v_k(N_{\mathrm{corr}})\right]$$

(3)

where “;” indicates that the elements $v_k(1)$, $v_k(2)$, … and $v_k(N_{\mathrm{corr}})$ are vertically arranged into a column vector. We can further define a PC weight matrix w that consists of N_pixel columns, namely,

$$\mathbf{w}=\left[\begin{array}{cccc}{\mathbf{w}}_1& {\mathbf{w}}_2& \cdots & {\mathbf{w}}_{N_{\mathrm{pixel}}}\end{array}\right]$$

(4)

where the vector containing PC weights for p-th pixel is

$${\mathbf{w}}_p=\left[w_p(1);w_p(2);\cdots ;w_p(N_{\mathrm{PC}})\right]$$

(5)

The quantities (${\overline{\mathbf{x}}}_{\mathrm{corr}}$, v, w) are contained in the solution vector x from the last iteration of an inversion process.

3. Inversion

3.1. State Vector

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

$${\mathbf{x}}_{\mathrm{state}}=\left[{\overline{\mathbf{x}}}_{\mathrm{corr}};{\mathbf{v}}_{k=1};{\mathbf{v}}_{k=2}; \cdots {\mathbf{v}}_{k=N_{\mathrm{PC}}}; {\mathbf{w}}_{p=1};{\mathbf{x}}_{p=1,\mathrm{uncorr}}; {\mathbf{w}}_{p=2};{\mathbf{x}}_{p=2,\mathrm{uncorr}}\cdots ; {\mathbf{w}}_{p=N_{\mathrm{pixel}}};{\mathbf{x}}_{p=N_{\mathrm{pixel}},\mathrm{uncorr}}\right]$$

(6)

where the mean of aerosol fields takes the form of a column vector

$${\overline{\mathbf{x}}}_{\mathrm{corr}}=\left[{\overline{\mathbf{x}}}_{\mathrm{corr},1};{\overline{\mathbf{x}}}_{\mathrm{corr},2};\cdots ;{\overline{\mathbf{x}}}_{\mathrm{corr},N_{\mathrm{TP},\mathrm{corr}}}\right]$$

(7)

and $N_{\mathrm{TP},\mathrm{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,

$${\overline{\mathbf{x}}}_{\mathrm{corr},j}=[{\overline{\mathbf{x}}}_{\mathrm{corr},j}(1);{\overline{\mathbf{x}}}_{\mathrm{corr},j}(2);\cdots {\overline{\mathbf{x}}}_{\mathrm{corr},j}(L(j))]$$

(8)

Similarly, an arbitrary ($k^{\mathrm{th}}$) PC vector composed of $N_{\mathrm{TP},\mathrm{corr}}$ types of correlated parameters can be arranged into the following column vector,

$${\mathbf{v}}_k=\left[{\mathbf{v}}_{k,1};{\mathbf{v}}_{k,2};\cdots ; {\mathbf{v}}_{k,N_{\mathrm{TP},\mathrm{corr}}}\right]$$

(9)

where for the type of correlated parameter with $L_{\mathrm{corr}}(j)$ elements,

$${\mathbf{v}}_{k,j}=\left[v_{k,j}(1);v_{k,j}(2);\dots ; v_{k,j}(L_{\mathrm{corr}}(j))\right]$$

(10)

The PC weight vector for pixel p is

$${\mathbf{w}}_p=\left[w_p(1);w_p(2);\cdots w_p(N_{\mathrm{PC}})\right]$$

(11)

For the state vector with $N_{\mathrm{TP},\mathrm{uncorr}}$ types of uncorrelated parameters, we have for an arbitrary pixel p,

$${\mathbf{x}}_{p,\mathrm{uncorr}}=\left[{\mathbf{x}}_{p,\mathrm{uncorr},1};{\mathbf{x}}_{p,\mathrm{uncorr},2};\cdots ;{\mathbf{x}}_{p,\mathrm{uncorr},N_{\mathrm{TP},\mathrm{uncorr}}}\right]$$

(12)

where the [$n_{\mathrm{TP},\mathrm{uncorr}}$]-th type of uncorrelated parameter has $L_{\mathrm{uncorr}}(n_{\mathrm{TP},\mathrm{uncorr}})$ elements, namely,

$${\mathbf{x}}_{\mathrm{uncorr},p,n_{\mathrm{TP},\mathrm{uncorr}}}=\left[x_{\mathrm{uncorr},p,n_{\mathrm{TP},\mathrm{uncorr}}}(1);x_{\mathrm{uncorr},p,n_{\mathrm{TP},\mathrm{uncorr}}}(2);\cdots ;x_{\mathrm{uncorr},p n_{\mathrm{TP},\mathrm{uncorr}}}(L_{\mathrm{uncorr}}(n_{\mathrm{TP},\mathrm{uncorr}}))\right]$$

(13)

3.2. 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:

• The first type of constraint, formulated in Section 3.2.1, consists of the observations provided directly by the remote sensing instrument(s).
• The second type of constraint consists of a priori values for the retrieval parameters, and is described in Section 3.2.2.
• For the uncorrelated fields (such as surface reflection properties), the pixel-resolved properties ${\mathbf{x}}_{\mathrm{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 [6] and is repeated in Section 3.2.3 as well as in Appendix B as the third type of constraint.
• 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.
• 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.
• 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.
• Similarly, application of the smoothness constraints to certain type of parameters within a PC vector ${\mathbf{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.
• The eighth type of constraint, formulated in Section 3.2.5, imposes a zero sum of the PC weights.
• The ninth type of constraint, formulated in Section 3.2.6, imposes mutual orthogonality of the PC vectors.
• 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 [24] for Equations (14)–(20) below):

$${\mathbf{f}}_i^{\ast }={\mathbf{f}}_i(\mathbf{x})+\Delta {\mathbf{f}}_i^{\ast },\hspace{1em}1\le i\le M$$

(14)

where ${\mathbf{f}}_i^{\ast }$ denotes the $i^{\mathrm{th}}$ type of constraint, $\Delta {\mathbf{f}}_i^{\ast }$ is the error with this type of constraint, and x = ${\mathbf{x}}_{\mathrm{state}}$, as defined in Equation (6). Formally, the statistical independence of different sources of constraints means that the covariance matrix of joint constraint ${\mathbf{f}}^{\ast }=[{\mathbf{f}}_1^{\ast };{\mathbf{f}}_2^{\ast };\cdots ;{\mathbf{f}}_M^{\ast }]$ has the following structure

$${\mathbf{C}}_{{\mathbf{f}}^{\ast }}=\left[\begin{array}{cccc}{\mathbf{C}}_1& \mathbf{0}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& {\mathbf{C}}_2& \mathbf{0}& \mathbf{0}\\ \cdots & \cdots & \ddots & \cdots \\ \mathbf{0}& \mathbf{0}& & {\mathbf{C}}_{\mathrm{M}}\end{array}\right]$$

(15)

where ${\mathbf{C}}_i$ indicates the covariance matrix of i-th constraint (${\mathbf{f}}_i^{\ast }$). Following the expressions of ${\mathbf{f}}_i^{\ast }$ and ${\mathbf{C}}_{{\mathbf{f}}_i^{\ast }}$, the probability distribution function (PDF) of joint data (1 ≤ i ≤ M) can be derived by multiplying PDFs of data from all M sources, namely,

$$\mathrm{P}(\mathbf{f}(\mathbf{x})|{\mathbf{f}}^{\ast })={\displaystyle \prod _{i=1}^M\mathrm{P}({\mathbf{f}}_i(\mathbf{x})|{\mathbf{f}}_i^{\ast })}~\exp \left\{-\frac{1}{2}{\displaystyle \sum _{i=1}^M{\left[{\mathbf{f}}_i(\mathbf{x})-{\mathbf{f}}_i^{\ast }\right]}^{\mathrm{T}}{({\mathbf{C}}_i)}^{-1}\left[{\mathbf{f}}_i(\mathbf{x})-{\mathbf{f}}_i^{\ast }\right]}\right\}$$

(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,

$${\Psi }_{\mathrm{total}}(\mathbf{x})={\displaystyle \sum _{i=1}^M{\gamma }_i{\Psi }_i(\mathbf{x})}$$

(17)

where

$${\Psi }_i(\mathbf{x})=\frac{1}{2}{\left[{\mathbf{f}}_i(\mathbf{x})-{\mathbf{f}}_i^{\ast }\right]}^{\mathrm{T}}{\mathbf{W}}_i^{-1}\left[{\mathbf{f}}_i(\mathbf{x})-{\mathbf{f}}_i^{\ast }\right]$$

(18)

$${\mathbf{W}}_i=\frac{1}{{\epsilon }_i^2}{\mathbf{C}}_i$$

(19)

$${\gamma }_i=\frac{{\epsilon }_1^2}{{\epsilon }_i^2}$$

(20)

In the above equations, ${\epsilon }_i^2$ is the first diagonal element of ${\mathbf{C}}_i$ (i.e., ${\epsilon }_i^2={\left\{{\mathbf{C}}_i\right\}}_{11}$) and the Lagrange factor ${\gamma }_i$ weights the contribution of each type of constraint with respect to the first one (${\gamma }_1$ = 1).

Minimization of ${\Psi }_{\mathrm{total}}(\mathbf{x})$ in Equation (17) means its gradient with respect to the solution x approaches zero, such that

$$\nabla {\Psi }_{\mathrm{total}}(\mathbf{x})={\displaystyle \sum _{i=1}^M{\gamma }_i\nabla {\Psi }_i(\mathbf{x})}=0$$

(21)

which can be ensured by enforcing the gradient of all components approach zero, namely

$$\nabla {\Psi }_i(\mathbf{x})={\mathbf{K}}_i^{\mathrm{T}}{\mathbf{W}}_i^{-1}({\mathbf{f}}_i(\mathbf{x})-{\mathbf{f}}_i^{\ast })=0$$

(22)

where ${\mathbf{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 (22) and substitute

$${\mathbf{f}}_i(\mathbf{x}-\Delta \mathbf{x})={\mathbf{f}}_i(\mathbf{x})-{\mathbf{K}}_i\Delta \mathbf{x}$$

(23)

into it. This results in

$$\nabla {\Psi }_i(\mathbf{x}-\Delta \mathbf{x})={\mathbf{K}}_i^{\mathrm{T}}{\mathbf{W}}_i^{-1}({\mathbf{f}}_i(\mathbf{x})-{\mathbf{K}}_i\Delta \mathbf{x}-{\mathbf{f}}_i^{\ast })=0$$

(24)

or equivalently,

$$({\mathbf{K}}_i^{\mathrm{T}}{\mathbf{W}}_i^{-1}{\mathbf{K}}_i)\Delta \mathbf{x}={\mathbf{K}}_i^{\mathrm{T}}{\mathbf{W}}_i^{-1}({\mathbf{f}}_i(\mathbf{x})-{\mathbf{f}}_i^{\ast })=\nabla {\Psi }_i(\mathbf{x})$$

(25)

The Jacobian matrix ${\mathbf{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)}={\frac{\partial {\mathbf{f}}_{i,l}}{\partial {\mathbf{x}}_n}|}_{\mathbf{x}}$$

(26)

More explicit evaluation of ${\mathbf{f}}_i(\mathbf{x})$, ${\mathbf{f}}_i^{\ast }$, and K and W matrices for all ten types of constraints is discussed in the following subsections.

3.2.1. Observational Constraints (i = 1)

Assuming a total number of $N_{\mathrm{pixel}}$ pixels, where each pixel has Z observations, then after arranging all pixel data into a column vector we have

$${\mathbf{f}}_{i=1}^{\ast }(\mathbf{x})={[{\mathbf{f}}_{1,p_1}^{\ast };{\mathbf{f}}_{1,p_2}^{\ast };\cdots ;{\mathbf{f}}_{1,N_{\mathrm{pixel}}}^{\ast }]}_{N\times 1}=[y_1;y_2;\cdots ;y_{N_{\mathrm{f}}}]$$

(27)

where $N_{\mathrm{f}}$ = $N_{\mathrm{pixel}}$× Z.

The parameters of the atmosphere-surface model are adjusted so that the model prediction of radiance and polarization ${\mathbf{f}}_1(\mathbf{x})$ fit the observational constraints ${\mathbf{f}}_1^{\ast }(\mathbf{x})$. The calculation of ${\mathbf{f}}_1(\mathbf{x})$ is introduced in Section 4. The Jacobian matrix ${\mathbf{K}}_1$ consists of first order partial derivatives with respect to correlated and uncorrelated parameters in the vicinity of x, namely,

$${\mathbf{K}}_1=\left[\begin{array}{ccccccccc}{\mathbf{K}}_{1,{\overline{\mathbf{x}}}_{\mathrm{corr}}}& {\mathbf{K}}_{1,\mathbf{v}}& {\mathbf{K}}_{1,\mathbf{w}}(p_1)& {\mathbf{K}}_{1,\mathrm{uncorr}}(p_1)& \mathbf{0}& \mathbf{0}& \cdots & \mathbf{0}& \mathbf{0}\\ {\mathbf{K}}_{1,{\overline{\mathbf{x}}}_{\mathrm{corr}}}& {\mathbf{K}}_{1,\mathbf{v}}& \mathbf{0}& \mathbf{0}& {\mathbf{K}}_{1,\mathbf{w}}(p_2)& {\mathbf{K}}_{1,\mathrm{uncorr}}(p_2)& \cdots & \mathbf{0}& \mathbf{0}\\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \ddots & \cdots & \cdots \\ {\mathbf{K}}_{1,{\overline{\mathbf{x}}}_{\mathrm{corr}}}& {\mathbf{K}}_{1,\mathbf{v}}& \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{0}& \cdots & {\mathbf{K}}_{1,\mathbf{w}}(N_{\mathrm{pixel}})& {\mathbf{K}}_{1,\mathrm{uncorr}}(N_{\mathrm{pixel}})\end{array}\right]$$

(28)

where ${\mathbf{K}}_{1,{\overline{\mathbf{x}}}_{\mathrm{corr}}}$ is the Jacobian matrix containing derivatives of observations with respect to spatial and temporal mean correlated parameters (total number is $N_{\mathrm{corr}}$). Variation of ${\overline{\mathbf{x}}}_{\mathrm{corr}}$ impacts observations in all pixels. Then ${\mathbf{K}}_{1,{\overline{\mathbf{x}}}_{\mathrm{corr}}}$ is evaluated by

$${\mathbf{K}}_{1,{\overline{\mathbf{x}}}_{\mathrm{corr}}}=\left[\begin{array}{cccc}\frac{\partial y_1}{\partial {\overline{x}}_{\mathrm{corr}}(1)}& \frac{\partial y_1}{\partial {\overline{x}}_{\mathrm{corr}}(2)}& \cdots & \frac{\partial y_1}{\partial {\overline{x}}_{\mathrm{corr}}(N_{\mathrm{corr}})}\\ \frac{\partial y_2}{\partial {\overline{x}}_{\mathrm{corr}}(1)}& \frac{\partial y_2}{\partial {\overline{x}}_{\mathrm{corr}}(2)}& \cdots & \frac{\partial y_2}{\partial {\overline{x}}_{\mathrm{corr}}(N_{\mathrm{corr}})}\\ \cdots & \cdots & \ddots & \cdots \\ \frac{\partial y_{N_f}}{\partial {\overline{x}}_{\mathrm{corr}}(1)}& \frac{\partial y_{N_f}}{\partial {\overline{x}}_{\mathrm{corr}}(2)}& \cdots & \frac{\partial y_{N_f}}{\partial {\overline{x}}_{\mathrm{corr}}(N_{\mathrm{corr}})}\end{array}\right]$$

(29)

The matrix ${\mathbf{K}}_{1,\mathbf{v}}$ in Equation (28) is the Jacobian matrix containing derivatives of observations with respect to all PC elements. It is evaluated as follows:

$${\mathbf{K}}_{1,\mathbf{v}}=\left[\begin{array}{cccc}{\mathbf{K}}_{1,\mathbf{v}}(k=1)& {\mathbf{K}}_{1,\mathbf{v}}(k=2)& \cdots & {\mathbf{K}}_{1,\mathbf{v}}(k=N_{\mathrm{PC}})\end{array}\right]$$

(30)

with

$${\mathbf{K}}_{1,\mathbf{v}}(k)=\left[\begin{array}{cccc}\frac{\partial y_1}{\partial v_k(1)}& \frac{\partial y_1}{\partial v_k(2)}& \cdots & \frac{\partial y_1}{\partial v_k(N_{\mathrm{corr}})}\\ \frac{\partial y_2}{\partial v_k(1)}& \frac{\partial y_2}{\partial v_k(2)}& \cdots & \frac{\partial y_2}{\partial v_k(N_{\mathrm{corr}})}\\ \cdots & \cdots & \ddots & \cdots \\ \frac{\partial y_{N_f}}{\partial v_k(1)}& \frac{\partial y_{N_f}}{\partial v_k(2)}& \cdots & \frac{\partial y_{N_f}}{\partial v_k(N_{\mathrm{corr}})}\end{array}\right]$$

(31)

The matrix ${\mathbf{K}}_{1,\mathbf{w}}$ in Equation (28) is the Jacobian matrix containing derivatives of observations with respect to pixel-resolved PC weights and is evaluated as follows:

$${\mathbf{K}}_{1,\mathbf{w}}(p)=\left[\begin{array}{cccc}\frac{\partial y_{p,1}}{\partial w_p(1)}& \frac{\partial y_{p,1}}{\partial w_p(2)}& \cdots & \frac{\partial y_{p,1}}{\partial w_p(N_{\mathrm{PC}})}\\ \frac{\partial y_{p,2}}{\partial w_p(1)}& \frac{\partial y_{p,2}}{\partial w_p(2)}& \cdots & \frac{\partial y_{p,2}}{\partial w_p(N_{\mathrm{PC}})}\\ \cdots & \cdots & \ddots & \cdots \\ \frac{\partial y_{p,Z}}{\partial w_p(1)}& \frac{\partial y_{p,Z}}{\partial w_p(2)}& \cdots & \frac{\partial y_{p,Z}}{\partial w_p(N_{\mathrm{PC}})}\end{array}\right]$$

(32)

Moreover, the derivatives with respect to pixel-resolved uncorrelated parameters are evaluated by,

$${\mathbf{K}}_{1,\mathrm{uncorr}}(p)=\left[\begin{array}{cccc}\frac{\partial y_{p,1}}{\partial x_p(1)}& \frac{\partial y_{p,1}}{\partial x_p(2)}& \cdots & \frac{\partial y_{p,1}}{\partial x_p(L(1)+L(2)+\dots +L(N_{\mathrm{uncorr}}))}\\ \frac{\partial y_{p,2}}{\partial x_p(1)}& \frac{\partial y_{p,2}}{\partial x_p(2)}& \cdots & \frac{\partial y_{p,2}}{\partial x_p(L(1)+L(2)+\dots +L(N_{\mathrm{uncorr}}))}\\ \cdots & \cdots & \ddots & \cdots \\ \frac{\partial y_{p,Z}}{\partial x_p(1)}& \frac{\partial y_{p,Z}}{\partial x_p(2)}& \cdots & \frac{\partial y_{p,Z}}{\partial x_p(L(1)+L(2)+\dots +L(N_{\mathrm{uncorr}}))}\end{array}\right]$$

(33)

The covariance matrix for the first constraint (i.e., the observations) assembles the sub-weighting matrix from all pixels, namely,

$${\mathbf{W}}_1=\left[\begin{array}{cccc}{\mathbf{W}}_1(1)& \mathbf{0}& \cdots & \mathbf{0}\\ \mathbf{0}& {\mathbf{W}}_1(2)& \cdots & \mathbf{0}\\ \cdots & \cdots & \ddots & \cdots \\ \mathbf{0}& \mathbf{0}& \cdots & {\mathbf{W}}_1(N_{\mathrm{pixel}})\end{array}\right]$$

(34)

3.2.2. A Priori Constraints (i = 2)

The a priori constraint is constructed in the same way as in a previous study [24]. Namely, Equation (14) becomes

$${\mathbf{f}}_{i=2}^{\ast }={\mathbf{x}}^{a priori}=\mathbf{x}+\Delta {\mathbf{x}}^{a priori}$$

(35)

Then in Equations (21)–(25), ${\mathbf{K}}_{i=2}$ = I (identity matrix) and ${\mathbf{W}}_{i=2}=\frac{1}{{\epsilon }_{a\ast }^2}{\mathbf{C}}_{a\ast }$. More explicitly, ${\mathbf{W}}_{i=2}$ can be constructed from estimated range of each parameter relative to the first one,

$${\mathbf{W}}_{i=2}=\left[\begin{array}{cccc}1& 0& \cdots & 0\\ 0& \frac{{(x_{2,\max }-x_{2,\min })}^2}{{(x_{1,\max }-x_{1,\min })}^2}& \cdots & 0\\ \cdots & \cdots & \ddots & \cdots \\ 0& 0& \cdots & \frac{{(x_{N,\max }-x_{N,\min })}^2}{{(x_{1,\max }-x_{1,\min })}^2}\end{array}\right]$$

(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

$$\{\begin{array}{l}{\mathbf{f}}_3^{\ast }={\mathbf{0}}^{\ast }={\mathbf{S}}_{3,\mathrm{m}}{\mathbf{x}}_{\mathrm{uncorr}}+{\Delta }_{{\mathrm{g}(\mathrm{x}}_{\mathrm{uncorr}})}^{\ast }\\ {\mathbf{f}}_4^{\ast }={\mathbf{0}}^{\ast }={\mathbf{S}}_{4,\mathrm{m}}\overline{\mathbf{x}}+{\Delta }_{{\mathrm{g}(\overline{\mathrm{x}}}_{\mathrm{corr}})}^{\ast }\\ {\mathbf{f}}_6^{\ast }={\mathbf{0}}^{\ast }={\mathbf{S}}_{6,\mathrm{m}}{\mathbf{w}}_{\mathrm{state}}+{\Delta }_{\mathrm{g}\left(\mathrm{w}\right)}^{\ast }\\ {\mathbf{f}}_7^{\ast }={\mathbf{0}}^{\ast }={\mathbf{S}}_{7,\mathrm{m}}{\mathbf{v}}_{\mathrm{state}}+{\Delta }_{\mathrm{g}\left(\mathrm{v}\right)}^{\ast }\end{array}$$

(37)

where ${\mathbf{w}}_{\mathrm{state}}$ and ${\mathbf{v}}_{\mathrm{state}}$ are column vectors containing only PC weights and vectors, respectively, and are extracted from the overall state vector expressed in Equation (6), and ${\mathbf{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 [6], such that

$$S_{i,m}(z)z=\frac{{\mathrm{d}}^{m}g(z)}{\mathrm{d}{z}^m}=0\Rightarrow \{\begin{array}{l}{g}_{m=1}(z)=\mathrm{const}\\ g_{m=2}(z)=Az+B\\ g_{m=3}(z)=A{z}^2+Bz+C\end{array}$$

(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=\{\begin{array}{l}\frac{\mathrm{d}g}{\mathrm{d}z}\approx \frac{{\Delta }^{1}g}{{\Delta }_1(z)}=\frac{g(z_{j+1})-g(z_j)}{{\Delta }_1(z_j)},\hspace{1em}\hspace{1em}\hspace{1em}\mathrm{for} m=1\\ \frac{{\mathrm{d}}^{m}g}{\mathrm{d}{z}^m}\approx \frac{{\Delta }^{m}g}{{\Delta }_m(z)}=\frac{{\Delta }^{m-1}g(z_{j+1})/{\Delta }_{m-1}(z_{j+1})-{\Delta }^{m-1}g(z_j)/{\Delta }_{m-1}(z_j)}{[{\Delta }_{m-1}(z_j)+{\Delta }_{m-1}(z_{j+1})]/2},\hspace{1em}\mathrm{for} m\ge 2\end{array}$$

(39)

Taking the orders of difference m = 1 and 2 as examples, we have

$$\{\begin{array}{l}{\Delta }_{m=1}(z_j)=z_{j+1}-z_j\\ {\Delta }_{m=2}(z_j)=[{\Delta }_1(z_j)+{\Delta }_1(z_{j+1})]/2\end{array}$$

(40)

Application of the above equation to L discretized grids r_j (namely, 1 ≤ j ≤ L) leads to

$${\mathbf{f}}_i(\mathbf{x})={\mathbf{S}}_{i,m}\mathbf{x}$$

(41)

so that by invoking Equation (26),

$${\mathbf{K}}_i={\mathbf{S}}_{i,m}$$

(42)

where the matrix ${\mathbf{S}}_{i,m}$ is evaluated by,

$${\mathbf{S}}_{i,m=1}=\left[\begin{array}{ccccc}\frac{1}{{\Delta }_1(1)}& -\frac{1}{{\Delta }_1(1)}& 0& \cdots & 0\\ 0& \frac{1}{{\Delta }_1(2)}& -\frac{1}{{\Delta }_1(2)}& \cdots & 0\\ 0& 0& \cdots & \ddots & 0\\ 0& 0& \cdots & \frac{1}{{\Delta }_1(L-1)}& -\frac{1}{{\Delta }_1(L-1)}\end{array}\right]$$

(43)

and

$${\mathbf{S}}_{i,m=2}=\left[\begin{array}{ccccccc}\frac{2}{{\Delta }_1(1)[{\Delta }_1(1)+{\Delta }_1(2)]}& \frac{-2}{{\Delta }_1(1){\Delta }_1(2)}& \frac{2}{{\Delta }_1(2)[{\Delta }_1(1)+{\Delta }_1(2)]}& 0& 0& \cdots & 0\\ 0& \frac{2}{{\Delta }_1(2)[{\Delta }_1(2)+{\Delta }_1(3)]}& \frac{-2}{{\Delta }_1(2){\Delta }_1(3)}& \frac{2}{{\Delta }_1(3)[{\Delta }_1(2)+{\Delta }_1(3)]}& & \cdots & 0\\ \cdots & \cdots & \cdots & \cdots & \cdots & \ddots & \cdots \\ 0& 0& \cdots & \frac{2}{{\Delta }_1(L-1)[{\Delta }_1(L-1)+{\Delta }_1(L)]}& \frac{-2}{{\Delta }_1(L-1){\Delta }_1(L)}& \cdots & \frac{2}{{\Delta }_1(L)[{\Delta }_1(L-1)+{\Delta }_1(L)]}\end{array}\right]$$

(44)

The same principle applied to higher orders of difference (m > 2) ensures a smooth curve with $\frac{{\mathrm{d}}^{m}g(z)}{\mathrm{d}{z}^m}=0$. Substitution of Equations (41), (42), and ${\mathbf{f}}_i^{\ast }={\mathbf{0}}^{\ast }$ into Equation (18) gives,

$${\Psi }_i(\mathbf{x})=\frac{1}{2}{\mathbf{x}}^{\mathrm{T}}{\mathbf{S}}_{i,m}^{\mathrm{T}}{\mathbf{W}}_{i,m}^{-1}{\mathbf{S}}_{i,m}\mathbf{x}$$

(45)

where the weighting matrix W has the following diagonal terms,

$${\left\{{\mathbf{W}}_{i,m}\right\}}_{jj}=\frac{1}{{\Delta }_m(z_j)}$$

(46)

and ${\Delta }_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 ${\mathbf{f}}_i^{\ast }={\mathbf{0}}^{\ast }$ into Equation (26) gives

$$({\mathbf{S}}_{i,m}^{\mathrm{T}}{\mathbf{W}}_{i,m}^{-1}{\mathbf{S}}_m)\Delta \mathbf{x}={\mathbf{S}}_{i,m}^{\mathrm{T}}{\mathbf{W}}_{i,m}^{-1}({\mathbf{S}}_{i,m}\mathbf{x})$$

(47)

Further defining

$${\mathbf{\Omega }}_i={\mathbf{S}}_{i,m}^{\mathrm{T}}{\mathbf{W}}_{i,m}^{-1}{\mathbf{S}}_{i,m}$$

(48)

we have

$${\mathbf{\Omega }}_i\Delta \mathbf{x}={\mathbf{\Omega }}_i\mathbf{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,

$${\mathbf{\Omega }}_i=\{\begin{array}{ll}{\mathbf{\Omega }}_{\mathrm{uncorr}},& i=3\\ {\mathbf{\Omega }}_{\mathrm{corr}}^{{\displaystyle \overline{\mathrm{x}}}},& i=4\\ {\mathbf{\Omega }}_{\mathrm{corr}}^{\mathrm{w}},& i=6\\ {\mathbf{\Omega }}_{\mathrm{corr}}^{\mathrm{v}},& i=7\end{array}$$

(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 ${\mathbf{w}}_{\mathrm{state}}=\left[\begin{array}{cccc}{\mathbf{w}}_{k=1};& {\mathbf{w}}_{k=2};& \cdots ;& {\mathbf{w}}_{k=N_{\mathrm{PC}}}\end{array}\right]$ and all PC elements for j-th correlated parameter into another single-column vector ${\mathbf{v}}_{\mathrm{state}}=\left[\begin{array}{cccc}{\mathbf{v}}_{k=1};& {\mathbf{v}}_{k=2};& \cdots ;& {\mathbf{v}}_{k=N_{\mathrm{PC}}}\end{array}\right]$. The smooth variation of the j-th correlated aerosol parameter is ensured by,

$${\mathbf{f}}_{i=5}^{\ast }={\mathbf{0}}^{\ast }={\mathbf{S}}_{i,m}{\mathbf{x}}_{\mathrm{wv}}+{\Delta }_{\mathrm{g}\left(\right[\mathrm{v};\mathrm{w}\left]\right)}^{\ast }$$

(51)

where ${\mathbf{x}}_{\mathrm{wv}}=\left[\begin{array}{l}{\mathbf{v}}_{\mathrm{state}}\\ {\mathbf{w}}_{\mathrm{state}}\end{array}\right]$. To evaluate the differentiation matrix ${\mathbf{S}}_{i,m}$, we recall that a pixel-resolved correlated parameter is constructed from PCs via Equation (1). An arbitrary correlated parameter $x_{\mathrm{corr},p}(j)$ associated with p-th pixel is derived as,

$$x_{\mathrm{corr},p}(j)={{\displaystyle \overline{x}}}_{\mathrm{corr},j}+{\displaystyle \sum _{k=1}^{N_{\mathrm{PC}}}{w}_p(k)v_k(j)}={{\displaystyle \overline{x}}}_{\mathrm{corr},j}+\left[w_p(1)v_1(j)+w_p(2)v_2(j)+\dots +w_p(N_{\mathrm{PC}})v_{N_{\mathrm{PC}}}(j)\right]$$

(52)

(a) across-pixel smoothness

The across-pixel variation of a correlated parameter is contributed by pixel-resolved PC weights, namely,

$$\Delta {\mathbf{x}}_{\mathrm{corr},p}=\mathbf{v}\times (\Delta {\mathbf{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=\frac{\mathrm{d}g}{\mathrm{d}z}\approx \frac{{\Delta }^{1}g}{{\Delta }_{1}z}\approx \frac{{\displaystyle \sum _{k=1}^{N_{\mathrm{PC}}}\left\{v_k(j)\times \left[w_{j+1}(k)-w_j(k)\right]\right\}}}{{\Delta }_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}$(2). The matrix ${\mathbf{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 ${\mathbf{f}}_{i=5}({\mathbf{x}}_{\mathrm{wv}})={\mathbf{S}}_{i,m}{\mathbf{x}}_{\mathrm{wv}}$:

$${\mathbf{f}}_{i=5}({\mathbf{x}}_{\mathrm{wv}})={\mathbf{S}}_{i,m}{\mathbf{x}}_{\mathrm{wv}}=\left\{\frac{1}{2}\left[\begin{array}{ccccc}{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)\end{array}\right]\right\}\times \left[\begin{array}{c}{v}_{k=1}(1)\\ v_{k=1}(2)\\ w_{p1}(k=1)\\ w_{p2}(k=1)\\ w_{p3}(k=1)\end{array}\right]$$

(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,

$$\Delta {\mathbf{x}}_{\mathrm{corr},p}=(\Delta \mathbf{v})\times {\mathbf{w}}_p$$

(56)

Taking the first order of difference (m = 1) as an example, the differentiation matrix ${\mathbf{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=\frac{\mathrm{d}g}{\mathrm{d}z}\approx \frac{{\Delta }^{m}g}{{\Delta }_{m}z}\approx \frac{{\displaystyle \sum _{k=1}^{N_{\mathrm{PC}}}\left\{w_p(k)\times \left[v_k(j+1)-v_k(j)\right]\right\}}}{{\Delta }_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,

$${\mathbf{f}}_{i=5}({\mathbf{x}}_{\mathrm{wv}})={\mathbf{S}}_{i,m=1}{\mathbf{x}}_{\mathrm{wv}}=\left\{\frac{1}{2}\left[\begin{array}{cccc}\frac{w_p(k=1)}{{\Delta }_1(z_1)}& -\frac{w_p(k=1)}{{\Delta }_1(z_1)}& 0& \frac{v_{k=1}(2)-v_{k=1}(1)}{{\Delta }_1(z_1)}\\ 0& \frac{w_p(k=1)}{{\Delta }_1(z_2)}& -\frac{w_p(k=1)}{{\Delta }_1(z_2)}& \frac{v_{k=1}(3)-v_{k=1}(2)}{{\Delta }_1(z_2)}\end{array}\right]\right\}\times \left[\begin{array}{c}{v}_{k=1}(1)\\ v_{k=1}(2)\\ v_{k=1}(3)\\ w_p(k=1)\end{array}\right]$$

(58)

Invoking ${\mathbf{f}}_i^{\ast }={\mathbf{0}}^{\ast }$, the cost function in the form of Equation (45) is derived for both across-pixel and within-pixel smoothness constraints. Further, substitution of ${\mathbf{f}}_{i=5}({\mathbf{x}}_{\mathrm{wv}})={\mathbf{S}}_{i,m=1}{\mathbf{x}}_{\mathrm{wv}}$ into Equation (26) gives the expression for Jacobian matrix elements,

$$K_{i,(l,n)}={\frac{\partial {({\mathbf{S}}_{i,m}{\mathbf{x}}_{\mathrm{wv}})}_l}{\partial x_{\mathrm{wv},n}}|}_{{\mathbf{x}}_{\mathrm{wv}}}$$

(59)

in which the matrix ${\mathbf{S}}_{i,m}$ is a function of x (see Equation (57)), therefore ${\mathbf{K}}_i\ne {\mathbf{S}}_{i,m}$. Intuitively, one can think of ${\mathbf{f}}_{i=5}({\mathbf{x}}_{\mathrm{wv}})={\mathbf{S}}_{i,m=1}{\mathbf{x}}_{\mathrm{wv}}$ as a model prediction to fit the “observation” ${\mathbf{f}}_i^{\ast }$, 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 ${\mathbf{f}}_i({\mathbf{x}}_{\mathrm{wv}})={\mathbf{S}}_{i,m}{\mathbf{x}}_{\mathrm{wv}}$ and ${\mathbf{f}}_i^{\ast }={\mathbf{0}}^{\ast }$, Equation (25) becomes

$$\left({\mathbf{K}}_i^{\mathrm{T}}{\mathbf{W}}_i^{-1}{\mathbf{K}}_i\right)\Delta {\mathbf{x}}_{\mathrm{wv}}={\mathbf{K}}_i^{\mathrm{T}}{\mathbf{W}}_i^{-1}{\mathbf{S}}_{i,m}{\mathbf{x}}_{\mathrm{wv}}$$

(60)

Further defining

$$\{\begin{array}{l}{\mathbf{\Omega }}_{\mathrm{corr},1}^{\mathrm{wv}}={\mathbf{K}}_i^{\mathrm{T}}{\mathbf{W}}_i^{-1}{\mathbf{K}}_i\\ {\mathbf{\Omega }}_{\mathrm{corr},2}^{\mathrm{wv}}={\mathbf{K}}_i^{\mathrm{T}}{\mathbf{W}}_i^{-1}{\mathbf{S}}_{i,m}\end{array}$$

(61)

we have

$${\mathbf{\Omega }}_{\mathrm{corr},1}\Delta {\mathbf{x}}_{\mathrm{wv}}={\mathbf{\Omega }}_{\mathrm{corr},2}{\mathbf{x}}_{\mathrm{wv}}$$

(62)

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

$${\mathbf{K}}_i={\mathbf{S}}_{i,m}$$

(63)

so that

$${\mathbf{\Omega }}_{\mathrm{corr},1}^{\mathrm{wv}}={\mathbf{\Omega }}_{\mathrm{corr},2}^{\mathrm{wv}}={\mathbf{S}}_{i,m}^{\mathrm{T}}{\mathbf{W}}_i^{-1}{\mathbf{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 ${\mathbf{S}}_{i,m}$, ${\mathbf{K}}_i$, and ${\mathbf{\Omega }}_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 $\mathbf{\Omega }$ matrix derived above is accounted for later to solve for an increment of the PC terms during optimization.

3.2.5. 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,

$${\mathbf{f}}_{i=8}^{\ast }={\mathbf{0}}^{\ast }={\mathbf{f}}_{i=8}({\mathbf{w}}_{\mathrm{state}})+{\Delta }_{\mathrm{O}}^{\ast }$$

(65)

where ${\mathbf{w}}_{\mathrm{state}}$ is a column vector defined in Section 3.2.4 and contains the weights of all PCs, and ${\mathbf{f}}_{i=8}({\mathbf{w}}_{\mathrm{state}})$ is expressed as

$${\mathbf{f}}_8({\mathbf{w}}_{\mathrm{state}})=\mathbf{O}{\mathbf{w}}_{\mathrm{state}}$$

(66)

where

$$\mathbf{O}={\left[\begin{array}{cccc}{\displaystyle \stackrel{\rightharpoonup }{1}}& \stackrel{\rightharpoonup }{0}& \cdots & \stackrel{\rightharpoonup }{0}\\ \stackrel{\rightharpoonup }{0}& {\displaystyle \stackrel{\rightharpoonup }{1}}& \cdots & \stackrel{\rightharpoonup }{0}\\ \cdots & \cdots & \ddots & \stackrel{\rightharpoonup }{0}\\ \stackrel{\rightharpoonup }{0}& \stackrel{\rightharpoonup }{0}& \cdots & {\displaystyle \stackrel{\rightharpoonup }{1}}\end{array}\right]}_{N_{\mathrm{PC}}\times (N_{\mathrm{PC}}\times N_{\mathrm{pixel}})},{\displaystyle \stackrel{\rightharpoonup }{1}}={[1, 1, \dots , 1]}_{1\times N_{\mathrm{pixel}}}$$

(67)

Further substituting ${\mathbf{f}}_8^{\ast } ={\mathbf{0}}^{\ast }$ and ${\mathbf{f}}_8({\mathbf{w}}_{\mathrm{state}})=\mathbf{O}{\mathbf{w}}_{\mathrm{state}}$ into Equation (18), the cost function has the following quadratic form,

$${\Psi }_i({\mathbf{w}}_{\mathrm{state}})=\frac{1}{2}({\mathbf{w}}_{\mathrm{state}}^{\mathrm{T}}{\mathbf{O}}^{\mathrm{T}}\mathbf{O}{\mathbf{w}}_{\mathrm{state}})$$

(68)

Substitution of Equation (66) into Equation (26) gives the evaluation of Jacobian matrix elements,

$$K_{i,(l,n)}={\frac{\partial {(\mathbf{O}{\mathbf{w}}_{\mathrm{state}})}_l}{\partial w_{\mathrm{state},n}}|}_{{\mathbf{w}}_{\mathrm{state}}}=O_{l,n}$$

(69)

Further defining a zero-sum matrix to be,

$${\mathbf{\Omega }}_{\mathrm{O}}={\mathbf{\Omega }}_i={\mathbf{O}}^{\mathrm{T}}\mathbf{O}$$

(70)

and invoking ${\mathbf{f}}_8({\mathbf{w}}_{\mathrm{state}})=\mathbf{O}{\mathbf{w}}_{\mathrm{state}}$ and ${\mathbf{f}}_8^{\ast } ={\mathbf{0}}^{\ast }$, Equation (25) becomes,

$${\mathbf{\Omega }}_{\mathrm{O}}\Delta {\mathbf{w}}_{\mathrm{state}}={\mathbf{\Omega }}_{\mathrm{O}}{\mathbf{w}}_{\mathrm{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 ${\mathbf{f}}_i$ (i = 9), namely

$${\mathbf{f}}_{i=9}^{\ast } ={\mathbf{0}}^{\ast }={\mathbf{f}}_{i=9}({\mathbf{v}}_{\mathrm{state}})+{\Delta }_{\Gamma }^{\ast }$$

(72)

where ${\mathbf{v}}_{\mathrm{state}}$ is a column vector defined in Section 3.2.4 and contains all PC vectors, and ${\mathbf{f}}_{i=9}({\mathbf{v}}_{\mathrm{state}})$ is expressed as

$${\mathbf{f}}_9({\mathbf{v}}_{\mathrm{state}})=\mathbf{\Gamma }{\mathbf{v}}_{\mathrm{state}}$$

(73)

where,

$$\mathbf{\Gamma }=\left[\begin{array}{ccccccc}{\mathbf{v}}_{k=2}^{\mathrm{T}}& {\mathbf{v}}_{k=1}^{\mathrm{T}}& \stackrel{\rightharpoonup }{0}& \stackrel{\rightharpoonup }{0}& \cdots & \stackrel{\rightharpoonup }{0}& \stackrel{\rightharpoonup }{0}\\ {\mathbf{v}}_{k=3}^{\mathrm{T}}& \stackrel{\rightharpoonup }{0}& {\mathbf{v}}_{k=1}^{\mathrm{T}}& \stackrel{\rightharpoonup }{0}& \cdots & \stackrel{\rightharpoonup }{0}& \stackrel{\rightharpoonup }{0}\\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ {\mathbf{v}}_{k=N_{\mathrm{PC}}}^{\mathrm{T}}& \stackrel{\rightharpoonup }{0}& \stackrel{\rightharpoonup }{0}& \stackrel{\rightharpoonup }{0}& \cdots & \stackrel{\rightharpoonup }{0}& {\mathbf{v}}_{k=1}^{\mathrm{T}}\\ \stackrel{\rightharpoonup }{0}& {\mathbf{v}}_{k=3}^{\mathrm{T}}& {\mathbf{v}}_{k=2}^{\mathrm{T}}& \stackrel{\rightharpoonup }{0}& \cdots & \stackrel{\rightharpoonup }{0}& \stackrel{\rightharpoonup }{0}\\ \stackrel{\rightharpoonup }{0}& {\mathbf{v}}_{k=4}^{\mathrm{T}}& \stackrel{\rightharpoonup }{0}& {\mathbf{v}}_{k=2}^{\mathrm{T}}& \cdots & \stackrel{\rightharpoonup }{0}& \stackrel{\rightharpoonup }{0}\\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ \stackrel{\rightharpoonup }{0}& {\mathbf{v}}_{k=N_{\mathrm{PC}}}^{\mathrm{T}}& \stackrel{\rightharpoonup }{0}& \stackrel{\rightharpoonup }{0}& \cdots & \stackrel{\rightharpoonup }{0}& {\mathbf{v}}_{k=2}^{\mathrm{T}}\\ \cdots & \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\ \stackrel{\rightharpoonup }{0}& \stackrel{\rightharpoonup }{0}& \stackrel{\rightharpoonup }{0}& \stackrel{\rightharpoonup }{0}& \cdots & {\mathbf{v}}_{k=N_{\mathrm{PC}}-1}^{\mathrm{T}}& {\mathbf{v}}_{k=N_{\mathrm{PC}}}^{\mathrm{T}}\end{array}\right]$$

(74)

Further substituting ${\mathbf{f}}_9^{\ast } ={\mathbf{0}}^{\ast }={\left\{{[\begin{array}{cccc}0& 0& \cdots & 0\end{array}]}^{\mathrm{T}}\right\}}_{\left\{\frac{1}{2}[N_{\mathrm{PC}}\times (N_{\mathrm{PC}}-1)]\right\}\times 1}$ and ${\mathbf{f}}_9({\mathbf{v}}_{\mathrm{state}})=\mathbf{\Gamma }{\mathbf{v}}_{\mathrm{state}}$ into Equation (18), the cost function is derived as,

$${\Psi }_i(\mathbf{v})=\frac{1}{2}{\mathbf{v}}_{\mathrm{state}}^{\mathrm{T}}{\mathbf{\Gamma }}^{\mathrm{T}}\mathbf{\Gamma }{\mathbf{v}}_{\mathrm{state}}$$

(75)

Substitution of Equation (73) into Equation (26) gives the evaluation of Jacobian matrix elements,

$$k_{i,(l,n)}={\frac{\partial {(\mathbf{\Gamma }\mathbf{v})}_l}{\partial v_{\mathrm{state},n}}|}_{{\mathbf{v}}_{\mathrm{state}}}={\mathbf{\Gamma }}_{l,n}$$

(76)

Further defining a zero-sum matrix to be,

$${\mathbf{\Omega }}_{\Gamma }={\mathbf{\Omega }}_i={\mathbf{\Gamma }}^{\mathrm{T}}\mathbf{\Gamma }$$

(77)

and invoking ${\mathbf{f}}_9({\mathbf{v}}_{\mathrm{state}})=\mathbf{\Gamma }{\mathbf{v}}_{\mathrm{state}}$ and ${\mathbf{f}}_9^{\ast } ={\mathbf{0}}^{\ast }$, Equation (25) becomes,

$${\mathbf{\Omega }}_{\Gamma }\Delta {\mathbf{v}}_{\mathrm{state}}={\mathbf{\Omega }}_{\Gamma }{\mathbf{v}}_{\mathrm{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.,${\Vert {\mathbf{v}}_k\Vert }^2={\mathbf{v}}_k^{\mathrm{T}}{\mathbf{v}}_k=1$). This forms the unity constraint imposed on all PC vectors via ${\mathbf{f}}_i$ (i = 10), such that

$${\mathbf{f}}_{i=10}^{\ast } ={\mathbf{1}}^{\ast }={\mathbf{f}}_{i=10}({\mathbf{v}}_{\mathrm{state}})+{\Delta }_{\mathrm{U}}^{\ast }$$

(79)

where

$${\mathbf{f}}_{10}({\mathbf{v}}_{\mathrm{state}})=\mathbf{U}{\mathbf{v}}_{\mathrm{state}}$$

(80)

and U is expressed as,

$$\mathbf{U}={\left[\begin{array}{cccc}{\mathbf{v}}_{k=1}^{\mathrm{T}}& \stackrel{\rightharpoonup }{0}& \cdots & \stackrel{\rightharpoonup }{0}\\ \stackrel{\rightharpoonup }{0}& {\mathbf{v}}_{k=2}^{\mathrm{T}}& \cdots & \stackrel{\rightharpoonup }{0}\\ \cdots & \cdots & \ddots & \stackrel{\rightharpoonup }{0}\\ \stackrel{\rightharpoonup }{0}& \stackrel{\rightharpoonup }{0}& \cdots & {\mathbf{v}}_{k=N_{\mathrm{PC}}}^{\mathrm{T}}\end{array}\right]}_{N_{\mathrm{PC}}\times (N_{\mathrm{PC}}\times N_{\mathrm{corr}})}$$

(81)

Then substituting Equation (80) and ${\mathbf{f}}_{10}^{\ast } ={\mathbf{1}}^{\ast }={\left\{{[\begin{array}{cccc}1& 1& \cdots & 1\end{array}]}^{\mathrm{T}}\right\}}_{N_{\mathrm{PC}}\times 1}$ into Equation (18), the cost function is derived as,

$${\Psi }_i({\mathbf{v}}_{\mathrm{state}})=\frac{1}{2}({\mathbf{v}}_{\mathrm{state}}^{\mathrm{T}}{\mathbf{U}}^{\mathrm{T}}\mathbf{U}{\mathbf{v}}_{\mathrm{state}}-N_{\mathrm{PC}})$$

(82)

Moreover, substitution of Equation (80) into Equation (26) gives the evaluation of Jacobian matrix elements:

$$K_{i,(l,n)}={\frac{\partial {(\mathbf{U}{\mathbf{v}}_{\mathrm{state}})}_l}{\partial v_{\mathrm{state},n}}|}_{{\mathbf{v}}_{\mathrm{state}}}$$

(83)

Further defining a unity constraining matrix to be,

$$\{\begin{array}{l}{\mathbf{\Omega }}_{\mathrm{U},1}={\mathbf{K}}_i^{\mathrm{T}}{\mathbf{K}}_i\\ {\mathbf{\Omega }}_{\mathrm{U},2}={\mathbf{K}}_i^{\mathrm{T}}\mathbf{U}\end{array}$$

(84)

and invoking ${\mathbf{f}}_{10}({\mathbf{v}}_{\mathrm{state}})=\mathbf{U}{\mathbf{v}}_{\mathrm{state}}$ and ${\mathbf{f}}_{10}^{\ast } ={\mathbf{1}}^{\ast }$, Equation (25) becomes equivalent to

$${\mathbf{\Omega }}_{\mathrm{U},1}\Delta {\mathbf{v}}_{\mathrm{state}}={\mathbf{\Omega }}_{\mathrm{U},2}{\mathbf{v}}_{\mathrm{state}}-{\mathbf{K}}_i^{\mathrm{T}}{\mathbf{f}}_i^{\ast }$$

(85)

3.2.8. 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,

$${\mathbf{x}}_{q+1}={\mathbf{x}}_q-\Delta {\mathbf{x}}_q$$

(86)

where $\Delta {\mathbf{x}}_q$ is obtained by accounting for both observational and a priori constraints derived in Section 3.2.1, Section 3.2.2, Section 3.2.3, Section 3.2.4, Section 3.2.5, Section 3.2.6 and Section 3.2.7,

$${\mathbf{A}}_q\times \Delta {\mathbf{x}}_q=\nabla {\Psi }_{\mathrm{total}}({\mathbf{x}}_q)$$

(87)

As the explicit form, we have

$${\mathbf{A}}_q={\mathbf{K}}_{1,q}^{\mathrm{T}}{\mathbf{W}}_1^{-1}{\mathbf{K}}_{1,q}+\gamma {\mathbf{\Omega }}_{\mathrm{total},1}+{\gamma }_{\mathrm{a}}{\mathbf{W}}_{\mathrm{a}}^{-1}$$

(88)

$$\nabla {\Psi }_{\mathrm{total}}({\mathbf{x}}_q)={\mathbf{K}}_{1,q}^{\mathrm{T}}{\mathbf{W}}_1^{-1}\left[{\mathbf{f}}_1({\mathbf{x}}_q)-{\mathbf{f}}_1^{\ast }\right]+\gamma {\mathbf{\Omega }}_{\mathrm{total},2}{\mathbf{x}}_q+{\gamma }_{\mathrm{a}}{\mathbf{W}}_{\mathrm{a}}^{-1}({\mathbf{x}}_q-{\mathbf{x}}^{a priori})+\left[\begin{array}{c}{\mathbf{0}}^{\ast }\\ -{\gamma }_{\mathrm{U}}{\mathbf{K}}_{10,q}^{\mathrm{T}}{\mathbf{f}}_{10}^{\ast }\\ {\mathbf{0}}^{\ast }\end{array}\right]$$

(89)

where “$\gamma {\mathbf{\Omega }}_{\mathrm{total}}$” incorporates i = 3–10 types of constraints, namely,

$$\{\begin{array}{l}\gamma {\mathbf{\Omega }}_{\mathrm{total},1}={\gamma }_{\mathrm{uncorr}}{\mathbf{\Omega }}_{\mathrm{uncorr}}^{\mathrm{Ra}}+{\gamma }_{\mathrm{corr}}^{{\displaystyle \overline{\mathrm{x}}}}{\mathbf{\Omega }}_{\mathrm{corr}}^{{\displaystyle \overline{\mathrm{x}}},\mathrm{Ra}}+{\gamma }_{\mathrm{corr}}^{\mathrm{wv}}{\mathbf{\Omega }}_{\mathrm{corr},1}^{\mathrm{wv},\mathrm{Ra}}+{\gamma }_{\mathrm{corr}}^{\mathrm{w}}{\mathbf{\Omega }}_{\mathrm{corr}}^{\mathrm{w},\mathrm{Ra}}+{\gamma }_{\mathrm{corr}}^{\mathrm{v}}{\mathbf{\Omega }}_{\mathrm{corr}}^{\mathrm{v},\mathrm{Ra}}+{\gamma }_{\mathrm{O}}{\mathbf{\Omega }}_{\mathrm{O}}^{\mathrm{Ra}}+{\gamma }_{\Gamma }{\mathbf{\Omega }}_{\Gamma }^{\mathrm{Ra}}+{\gamma }_{\mathrm{U}}{\mathbf{\Omega }}_{\mathrm{U},1}^{\mathrm{Ra}}\\ \gamma {\mathbf{\Omega }}_{\mathrm{total},2}={\gamma }_{\mathrm{uncorr}}{\mathbf{\Omega }}_{\mathrm{uncorr}}^{\mathrm{Ra}}+{\gamma }_{\mathrm{corr}}^{{\displaystyle \overline{\mathrm{x}}}}{\mathbf{\Omega }}_{\mathrm{corr}}^{{\displaystyle \overline{\mathrm{x}}},\mathrm{Ra}}+{\gamma }_{\mathrm{corr}}^{\mathrm{wv}}{\mathbf{\Omega }}_{\mathrm{corr},2}^{\mathrm{wv},\mathrm{Ra}}+{\gamma }_{\mathrm{corr}}^{\mathrm{w}}{\mathbf{\Omega }}_{\mathrm{corr}}^{\mathrm{w},\mathrm{Ra}}+{\gamma }_{\mathrm{corr}}^{\mathrm{v}}{\mathbf{\Omega }}_{\mathrm{corr}}^{\mathrm{v},\mathrm{Ra}}+{\gamma }_{\mathrm{O}}{\mathbf{\Omega }}_{\mathrm{O}}^{\mathrm{Ra}}+{\gamma }_{\Gamma }{\mathbf{\Omega }}_{\Gamma }^{\mathrm{Ra}}+{\gamma }_{\mathrm{U}}{\mathbf{\Omega }}_{\mathrm{U},2}^{\mathrm{Ra}}\end{array}$$

(90)

where the eight terms on the right-hand-side of the above equation are equal to ${\mathbf{\Omega }}_i$ with i = 3–10, respectively, and the multipliers ${\gamma }_{\dots }$ control the strength of these constraints. The smoothness matrices ${\gamma }_{\mathrm{uncorr}}{\mathbf{\Omega }}_{\mathrm{uncorr}}^{\mathrm{Ra}}$, ${\gamma }_{\mathrm{corr}}^{{\displaystyle \overline{\mathrm{x}}}}{\mathbf{\Omega }}_{\mathrm{corr}}^{{\displaystyle \overline{\mathrm{x}}},\mathrm{Ra}}$, ${\gamma }_{\mathrm{O}}{\mathbf{\Omega }}_{\mathrm{O}}^{\mathrm{Ra}}$, ${\gamma }_{\Gamma }{\mathbf{\Omega }}_{\Gamma }^{\mathrm{Ra}}$, ${\gamma }_{\mathrm{U}}{\mathbf{\Omega }}_{\mathrm{U}}^{\mathrm{Ra}}$, ${\gamma }_{\mathrm{corr}}^{\mathrm{wv}}{\mathbf{\Omega }}_{\mathrm{corr}}^{\mathrm{wv},\mathrm{Ra}}$, and ${\gamma }_{\mathrm{corr}}^{\mathrm{w}/\mathrm{v}}{\mathbf{\Omega }}_{\mathrm{corr}}^{\mathrm{w}/\mathrm{v},\mathrm{Ra}}$ 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 ${\gamma }_{\mathrm{uncorr}}{\mathbf{\Omega }}_{\mathrm{uncorr}}$, ${\gamma }_{\mathrm{corr}}^{{\displaystyle \overline{\mathrm{x}}}}{\mathbf{\Omega }}_{\mathrm{corr}}^{{\displaystyle \overline{\mathrm{x}}}}$, ${\gamma }_{\mathrm{O}}{\mathbf{\Omega }}_{\mathrm{O}}$,${\gamma }_{\Gamma }{\mathbf{\Omega }}_{\Gamma }$, ${\gamma }_{\mathrm{U}}{\mathbf{\Omega }}_{\mathrm{U}}$, ${\gamma }_{\mathrm{corr}}^{\mathrm{wv}}{\mathbf{\Omega }}_{\mathrm{corr}}^{\mathrm{wv}}$, and ${\gamma }_{\mathrm{corr}}^{\mathrm{w}/\mathrm{v}}{\mathbf{\Omega }}_{\mathrm{corr}}^{\mathrm{w}/\mathrm{v}}$, respectively, as evaluated in Section 3.2.3, Section 3.2.4, Section 3.2.5, Section 3.2.6 and Section 3.2.7. However, the elements of these matrices are rearranged into the same dimension to add to each other in assembling $\gamma {\mathbf{\Omega }}_{\mathrm{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 ${\gamma }_{\mathrm{uncorr}}{\mathbf{\Omega }}_{\mathrm{uncorr}}$ as well as ${\gamma }_{\mathrm{corr}}^{{\displaystyle \overline{\mathrm{x}}}}{\mathbf{\Omega }}_{\mathrm{corr}}^{{\displaystyle \overline{\mathrm{x}}}}$ 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 (${\gamma }_{\mathrm{corr}}^{\mathrm{wv}}{\mathbf{\Omega }}_{\mathrm{corr}}^{\mathrm{wv}}$), which acts on combined PC weights and vectors, and the decoupled type ${\gamma }_{\mathrm{corr}}^{\mathrm{w}/\mathrm{v}}{\mathbf{\Omega }}_{\mathrm{corr}}^{\mathrm{w}/\mathrm{v}}$, which act on PC weights and vectors separately. Based on previous demonstration of principles for deriving ${\mathbf{\Omega }}_i$ for i = 5, a comprehensive evaluation of ${\gamma }_{\mathrm{corr}}^{\mathrm{wv}}{\mathbf{\Omega }}_{\mathrm{corr}}^{\mathrm{wv}}$ is provided in Appendix C. Comprehensive evaluations of ${\gamma }_{\mathrm{corr}}^{\mathrm{w}}{\mathbf{\Omega }}_{\mathrm{corr}}^{\mathrm{w}}$ and ${\gamma }_{\mathrm{corr}}^{\mathrm{v}}{\mathbf{\Omega }}_{\mathrm{corr}}^{\mathrm{v}}$ are provided in Appendix D. The imposition of ${\gamma }_{\mathrm{corr}}^{\mathrm{w}}{\mathbf{\Omega }}_{\mathrm{corr}}^{\mathrm{w}}$ and ${\gamma }_{\mathrm{corr}}^{\mathrm{v}}{\mathbf{\Omega }}_{\mathrm{corr}}^{\mathrm{v}}$ 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 ${\gamma }_{\mathrm{corr}}^{\mathrm{wv}}{\mathbf{\Omega }}_{\mathrm{corr}}^{\mathrm{wv}}$ and ${\gamma }_{\mathrm{corr}}^{\mathrm{w}/\mathrm{v}}{\mathbf{\Omega }}_{\mathrm{corr}}^{\mathrm{w}/\mathrm{v}}$ stabilize the inversion of PCs.

Ideally, a retrieval is deemed successful when the minimization of the cost function is achieved, such that

$${\Psi }_{\mathrm{total}}\le (N_{\mathrm{f}}+N_{\mathrm{c}}+N_{{\mathrm{a}}^{\ast }}-N_{\mathrm{a}}){\epsilon }_{\mathrm{f}}^2$$

(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 ${\epsilon }_{\mathrm{f}}^2$ 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, ${\epsilon }_{\mathrm{c}}^2$, namely,

$$\frac{|{\Psi }_{\mathrm{total}}({\mathbf{x}}_{q+1})-{\Psi }_{\mathrm{total}}({\mathbf{x}}_q)|}{{\Psi }_{\mathrm{total}}({\mathbf{x}}_q)}\le {\epsilon }_{\mathrm{c}}^2.$$

(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,

$${\gamma }_i={\epsilon }_{i=1}^2/{\epsilon }_i^2$$

(93)

where ${\epsilon }_i^2$ are the first diagonal elements of the covariance matrices corresponding to i-th type of constraints. To estimate ${\epsilon }_i^2$ for a given parameter to be retrieved (f = x), which is a function of z, the most unsmooth known solution $f^{\mathrm{us}}(z)$ over the target domain is used, as suggested by Dubovik and King [25], namely,

$${\epsilon }^2={{\displaystyle {\int }_{z_{\min }}^{z_{\max }}\left(\frac{{\mathrm{d}}^m[f^{\mathrm{us}}(z)]}{{\mathrm{d}}^{m}z}\right)}}^2\mathrm{d}z$$

(94)

where $z_{\min }$ and $z_{\max }$ specify the lower and upper bound of z. For i = 3, 4, 5, ${\gamma }_i$ is evaluated using least smooth known solution of correlated parameter x over the target domain. For i = 6 and 7, ${\gamma }_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:

$${\gamma }_i^{\mathrm{Final}}=\frac{N_{\mathrm{f}}}{N_i}\frac{{\tilde{\epsilon }}_{\mathrm{f}}^2}{{\epsilon }_{\mathrm{f}}^2}{\gamma }_i$$

(95)

There are two differences between ${\gamma }_i^{\mathrm{Final}}$ and ${\gamma }_i$:

1. The multipliers “$N_f$/$N_i$” are introduced to account for possible redundancy of the measured and a priori data. Considering that ${\epsilon }_i^2$ 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.

2. The multiplier ${\tilde{\epsilon }}_{\mathrm{f}}^2$/${\epsilon }_{\mathrm{f}}^2$ is introduced with ${\tilde{\epsilon }}_{\mathrm{f}}^2$ estimated as the dynamic fitting residual during iterations:

$${\tilde{\epsilon }}_{\mathrm{f}}^2({\mathbf{x}}_q)\approx \frac{{\Psi }_{\mathrm{total}}({\mathbf{x}}_q)}{(N_{\mathrm{f}}+N_{\mathrm{c}}+N_{{\mathrm{a}}^{\ast }}-N_{\mathrm{a}}){\epsilon }_{\mathrm{f}}^2}$$

(96)

With the multiplier ${\tilde{\epsilon }}_{\mathrm{f}}^2$/${\epsilon }_{\mathrm{f}}^2$, 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 [26,27]. 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 ${\tilde{\epsilon }}_{\mathrm{f}}^2$ becomes equal to ${\epsilon }_{\mathrm{f}}^2$. 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.

3.3. 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 [19]. If not, then a practical method is to use the square root of the diagonal terms of the following covariance matrix,

$$\Delta x_j=\sqrt{{({\mathbf{C}}_{\Delta \mathbf{x}})}_{jj}}$$

(97)

with

$${\mathbf{C}}_{\Delta \mathbf{x}}=(\Delta {\mathbf{x}}_{\mathrm{syst}})\times {(\Delta {\mathbf{x}}_{\mathrm{syst}})}^{\mathrm{T}}+{\mathbf{C}}_{\Delta \mathbf{x},\mathrm{rand}}$$

(98)

where the covariance matrix ${\mathbf{C}}_{\Delta \mathbf{x},\mathrm{rand}}$ accounts for the contributions by random errors in the measurement (${\epsilon }_{\mathrm{rand}}$),

$${\mathbf{C}}_{\Delta {\mathrm{x}}_{\mathrm{rand}}}={\left[\mathbf{A}({\mathbf{x}}^{\mathrm{true}})\right]}^{-1}{\epsilon }_{\mathrm{rand}}^2$$

(99)

and systematic error $\Delta {\mathbf{x}}_{\mathrm{syst}}$ is estimated by

$$\Delta {\mathbf{x}}_{\mathrm{syst}}={\mathbf{\left[}\mathbf{A}({\mathbf{x}}^{\mathrm{true}})\right]}^{-1}\nabla {\Psi }_{\mathrm{total}}({\mathbf{x}}^{\mathrm{true}})$$

(100)

where A and $\nabla {\Psi }_{\mathrm{total}}$ are calculated by substituting ${\mathbf{x}}^{\mathrm{true}}$ into Equations (88) and (89), respectively.

As the retrieved solution is the closest estimate of the “true” solution ${\mathbf{x}}^{\mathrm{true}}$, we adopt ${\mathbf{x}}^{\mathrm{true}}$ = ${\mathbf{x}}^{\mathrm{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,

$${\left(\Delta {\mathbf{x}}_{p,\mathrm{corr}}^{\mathrm{retrieved}}\right)}^2={(\Delta {\overline{\mathbf{x}}}_{\mathrm{corr}}^{\mathrm{retrieved}})}^2+{\displaystyle \sum _{k=1}^{N_{\mathrm{PC}}}\left\{{[w_p^{\mathrm{retrieved}}(k)]}^2{(\Delta {\mathbf{v}}_k^{\mathrm{retrieved}})}^2+{[\Delta w_p^{\mathrm{retrieved}}(k)]}^2{({\mathbf{v}}_k^{\mathrm{retrieved}})}^2\right\}}$$

(101)

To estimate errors for functions of the retrieved parameters (namely $f({\mathbf{x}}_p^{\mathrm{retrieved}})$, e.g., aerosol single-scattering albedo, (SSA)), the following chain rule is applied,

$$\Delta f({\mathbf{x}}_p^{\mathrm{retrieved}})=\sqrt{{\mathbf{K}}_2^{\mathrm{T}}({\mathbf{x}}_p^{\mathrm{retrieved}}){\mathbf{C}}_{\Delta {\mathbf{x}}_p^{\mathrm{retrieved}}}{\mathbf{K}}_2({\mathbf{x}}_p^{\mathrm{retrieved}})}$$

(102)

where ${\mathbf{K}}_2$ denotes the Jacobian array containing derivatives of $f({\mathbf{x}}_p^{\mathrm{retrieved}})$ with respect to ${\mathbf{x}}_p^{\mathrm{retrieved}}$.

When ${\mathbf{x}}^{a priori}$ is not available, we assume ${\mathbf{x}}^{a priori}$ = ${\mathbf{x}}_q$ so that the a priori term on the right-hand-side of Equation (89) for $\nabla {\Psi }_{\mathrm{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.

3.4. 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,

$${\mathbf{x}}_{\mathrm{corr},p}={\overline{\mathbf{x}}}_{\mathrm{corr}}+\mathbf{c}\circ (\mathbf{v}\times {\mathbf{w}}_p)$$

(103)

where $\circ$ in the above equation denotes an operation of element-wise multiplication, and the constant column vector c has $N_{\mathrm{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 (${\sigma }_{\mathrm{spatio}-\mathrm{temporal}, x}$), and (b) its uncertainty estimate (${\sigma }_{\mathrm{e}, x}$), namely,

$$C_x=\sqrt{{\sigma }_{\mathrm{spatio}-\mathrm{temporal}, x}^2+{\sigma }_{\mathrm{e}, x}^2}$$

(104)

Such a scaled PC analysis is used to analyze the AERONET retrievals in Figure 1 and Figure 2. 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.

4. 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 [28] 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 [29]. 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 [30] 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 [31] 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 [32] 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 [19]. 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.

4.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.

4.1.1. 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_{\mathrm{PC}}$) dominant PC vectors. Taking advantage of orthogonality of the PCs, the quantities (Y), including total AOD (${\tau }_{\mathrm{aer},\mathrm{tot}}$), absorption AOD (${\tau }_{\mathrm{aer},\mathrm{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 ${\mathbf{v}}_k$. Adopting the second order of approximation (justified later in this section) and finite difference method to calculate derivatives, we have,

$$\mathbf{Y}({\mathbf{x}}_p)\approx \mathbf{Y}(\overline{\mathbf{x}})+{\displaystyle \sum _{k=1}^{N_{\mathrm{PC}}}\left[\frac{\mathbf{Y}(\overline{\mathbf{x}}+{\delta }_{\mathrm{s}}\times {\mathbf{v}}_k)-\mathbf{Y}(\overline{\mathbf{x}}-{\delta }_{\mathrm{s}}\times {\mathbf{v}}_k)}{2{\delta }_{\mathrm{s}}}{w}_p(k)+\frac{\mathbf{Y}(\overline{\mathbf{x}}+{\delta }_{\mathrm{s}}\times {\mathbf{v}}_k)-2\mathbf{Y}(\overline{\mathbf{x}})+\mathbf{Y}(\overline{\mathbf{x}}-{\delta }_{\mathrm{s}}\times {\mathbf{v}}_k)}{2{\delta }_{\mathrm{s}}^2}{w}_p^2(k)\right]}$$

(105)

where $\overline{\mathbf{x}}$ is state vector containing multi-pixel mean aerosol properties and ${\delta }_{\mathrm{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_{\mathrm{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,

$$\begin{array}{ll}\mathbf{Y}({\mathbf{x}}_p)-{\mathbf{Y}}^{\mathrm{LS}}({\mathbf{x}}_p)\approx & \left[{\mathbf{Y}}^{\mathrm{HS}}(\overline{\mathbf{x}})-{\mathbf{Y}}^{\mathrm{LS}}(\overline{\mathbf{x}})\right]+\\ & {\displaystyle \sum _{k=1}^{N_{\mathrm{PC}}}\left[\frac{\Delta \mathbf{Y}(\overline{\mathbf{x}}+{\delta }_{\mathrm{s}}\times {\mathbf{v}}_k)-\Delta \mathbf{Y}(\overline{\mathbf{x}}-{\delta }_{\mathrm{s}}\times {\mathbf{v}}_k)}{2{\delta }_{\mathrm{s}}}{w}_p(k)+\frac{\Delta \mathbf{Y}(\overline{\mathbf{x}}+{\delta }_{\mathrm{s}}\times {\mathbf{v}}_k)-2\Delta \mathbf{Y}(\overline{\mathbf{x}})+\Delta \mathbf{Y}(\overline{\mathbf{x}}-{\delta }_{\mathrm{s}}\times {\mathbf{v}}_k)}{2{\delta }_{\mathrm{s}}^2}{w}_p^2(k)\right]}\end{array}$$

(106)

where “HS” and “LS” denote high and low stream-based RT runs, respectively, and $\Delta \mathbf{Y}={\mathbf{Y}}^{\mathrm{HS}}-{\mathbf{Y}}^{\mathrm{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 [8] 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_{\mathrm{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 $\mathsf{\pi }{I}_{\mathrm{meas}}{d}_{\mathrm{ES}}^2/{\mu }_0{E}_0$, where $I_{\mathrm{meas}}$ is the measured radiance, $d_{\mathrm{ES}}$ is the Earth-Sun distance, ${\mu }_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 [20]. 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 [33]).

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.

4.1.2. Coupling Atmospheric Radiation with Surface Reflection

Using the deterministic Markov chain RT model for a plane-parallel atmosphere ([28], or other RT models), we calculate two sets of reflection and transmission matrices for the atmosphere, (${\mathbf{R}}_{\mathrm{atmos}}$, ${\mathbf{T}}_{\mathrm{atmos}}$) and (${\mathbf{R}}_{\mathrm{atmos}}^{\ast }$, ${\mathbf{T}}_{\mathrm{atmos}}^{\ast }$), which result from illumination at the top and bottom of the atmosphere, respectively. We denote the surface reflection matrix by ${\mathbf{R}}_{\mathrm{surf}}$, which consists of a BRDF component (${\mathbf{R}}_{\mathrm{surf}, \mathrm{BRDF}}$) and a pBRDF component (${\mathbf{R}}_{\mathrm{surf}, \mathrm{pBRDF}}$). In accordance with the adding method [29], two matrices Q and S are defined to account for the interaction between the surface and atmosphere, respectively,

$$\mathbf{S}={\displaystyle \sum _{n=1}^{N_{\max }}{\mathbf{Q}}_n}$$

(107)

$${\mathbf{Q}}_{\mathrm{n}}={\mathbf{Q}}_1{\mathbf{Q}}_{n-1}$$

(108)

$${\mathbf{Q}}_1={\mathbf{R}}_{\mathrm{atmos}}^{\ast }{\mathbf{R}}_{\mathrm{surf}}$$

(109)

Then the matrices for downwelling and upwelling diffuse light at the atmosphere-land interface are given by

$$\mathbf{D}={\mathbf{T}}_{\mathrm{atmos}}+\mathbf{S}\exp (-\frac{{\tau }_{\mathrm{atmos}}}{{\mu }_0})+\mathbf{S}{\mathbf{T}}_{\mathrm{atmos}}$$

(110)

where ${\tau }_{\mathrm{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

$$\mathbf{U}={\mathbf{R}}_{\mathrm{Surf}}\exp (-\frac{{\tau }_{\mathrm{atmos}}}{{\mu }_0})+{\mathbf{R}}_{\mathrm{Surf}}\mathbf{D}$$

(111)

The reflection matrix of the full CAS is,

$${\mathbf{R}}_{\mathrm{CAS}}={\mathbf{R}}_{\mathrm{atmos}}\exp (-\frac{{\tau }_{\mathrm{atmos}}}{\mu })\mathbf{U}+{\mathbf{T}}_{\mathrm{atmos}}^{\ast }\mathbf{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,

$${\mathrm{BRF}}_{\mathrm{TOA}}=\mathsf{\pi }{\displaystyle \sum _{m=0}^{\infty }(2-{\delta }_{0m}){\mathbf{R}}_{\mathrm{CAS},11}^{\left(\mathrm{m}\right)}\cos m({\varphi }_{\mathrm{v}}-{\varphi }_0)}$$

(113)

$${\mathrm{qBRF}}_{\mathrm{TOA}}=\mathsf{\pi }{\displaystyle \sum _{m=0}^{\infty }(2-{\delta }_{0m}){\mathbf{R}}_{\mathrm{CAS},21}^{\left(\mathrm{m}\right)}\cos m({\varphi }_{\mathrm{v}}-{\varphi }_0)}$$

(114)

$${\mathrm{uBRF}}_{\mathrm{TOA}}=\mathsf{\pi }{\displaystyle \sum _{m=0}^{\infty }(2-{\delta }_{0m}){\mathbf{R}}_{\mathrm{CAS},31}^{\left(\mathrm{m}\right)}\cos m({\varphi }_{\mathrm{v}}-{\varphi }_0)}$$

(115)

$${\mathrm{vBRF}}_{\mathrm{TOA}}=\mathsf{\pi }{\displaystyle \sum _{m=0}^{\infty }(2-{\delta }_{0m}){\mathbf{R}}_{\mathrm{CAS},41}^{\left(\mathrm{m}\right)}\cos m({\varphi }_{\mathrm{v}}-{\varphi }_0)}$$

(116)

where $({\varphi }_{\mathrm{v}}-{\varphi }_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 $\mathrm{DOLP}\approx \sqrt{{\mathrm{qBRF}}_{\mathrm{TOA}}^2+{\mathrm{uBRF}}_{\mathrm{TOA}}^2}/{\mathrm{BRF}}_{\mathrm{TOA}}$, where we neglect the minor contribution of circular polarization related term (${\mathrm{vBRF}}_{\mathrm{TOA}}$) when it is excluded from polarimetric measurements.

4.2. 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 (${\mathbf{f}}_1^{\ast }=[y_1;y_2;\cdots ;y_{N_f}]$) and a state vector of $N_a$ components ($\mathbf{x}=[x_1;x_2;\cdots ;x_{N_a}]$), the Jacobian matrix has the following structure,

$$\mathbf{K}={\left[\begin{array}{cccc}\frac{\partial y_1}{\partial x_1}& \frac{\partial y_1}{\partial x_2}& \cdots & \frac{\partial y_1}{\partial x_{N_a}}\\ \frac{\partial y_2}{\partial x_1}& \frac{\partial y_2}{\partial x_2}& \cdots & \frac{\partial y_2}{\partial x_{N_a}}\\ \cdots & \cdots & \ddots & \cdots \\ \frac{\partial y_{N_f}}{\partial x_1}& \frac{\partial y_{N_f}}{\partial x_2}& \cdots & \frac{\partial y_{N_f}}{\partial x_{N_a}}\end{array}\right]}_{N_f\times 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,

$$\frac{\partial y_i}{\partial x_n}=\frac{y_i(x_n+\Delta x_n)-y_i(x_n-\Delta x_n)}{2\Delta x_n}$$

(118)

The above equation applies to the calculation of the derivatives of ${\tau }_{\mathrm{aer},\mathrm{tot}}$, ${\tau }_{\mathrm{aer},\mathrm{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 = {${\tau }_{\mathrm{aer},\mathrm{tot}}$, ${\tau }_{\mathrm{aer},\mathrm{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,

$$\frac{\partial {\mathbf{Y}}_p}{\partial v_k(i)}\approx w_p(k)\frac{\partial {\mathbf{Y}}_p}{\partial \overline{x}(i)}$$

(119)

Such a strategy further improves the Jacobian evaluation efficiency when multiple PCs are retrieved.

5. 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 [19]. 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_{{\mathrm{v}}_{1-N_{\mathrm{S}\mathrm{C}}}}$ with $N_{\mathrm{SC}}$ = 5 as adopted in [19]), real and imaginary parts of refractive index ($n_{\mathrm{r}}$ and $n_{\mathrm{i}}$), Gaussian distribution-based vertical profile parameterized by central height $h_{\mathrm{a}}$ and standard deviation $s_{\mathrm{a}}$, and the volume fraction of spherical particles ($f_{\mathrm{v}, \mathrm{sphere}}$). Then, the parameter spaces described at the beginning of Section 3 are specified as follows,

$${\overline{\mathbf{x}}}_{\mathrm{corr}}=\log {\left[\stackrel{{\overline{\mathbf{x}}}_{\mathrm{corr},1}^{\mathrm{T}}}{\overbrace{{\overline{C}}_{\mathrm{v}1} \cdots {\overline{C}}_{{\mathrm{v}}_{N_{\mathrm{SC}}}}}},\stackrel{{\overline{\mathbf{x}}}_{\mathrm{corr},2}^{\mathrm{T}}}{\overbrace{{\overline{n}}_{\mathrm{r}, \lambda 1}\cdots {\overline{n}}_{\mathrm{r}, \lambda 7}}},\stackrel{{\overline{\mathbf{x}}}_{\mathrm{corr},3}^{\mathrm{T}}}{\overbrace{{\overline{n}}_{\mathrm{i}, \lambda 1}\cdots {\overline{n}}_{\mathrm{i}, \lambda 7}}},\stackrel{{\overline{\mathbf{x}}}_{\mathrm{corr},4}}{\overbrace{{\overline{h}}_{\mathrm{a}}}},\stackrel{{\overline{\mathbf{x}}}_{\mathrm{corr},5}}{\overbrace{{\overline{s}}_{\mathrm{a}}}}, \stackrel{{\overline{\mathbf{x}}}_{\mathrm{corr},6}}{\overbrace{{\overline{f}}_{\mathrm{v}, \mathrm{sphere}}}}\right]}^{\mathrm{T}}$$

(120)

$$\mathbf{v}=\left[{\mathbf{v}}_1 {\mathbf{v}}_2 \dots {\mathbf{v}}_{N_{\mathrm{PC}}}\right]$$

(121)

with

$${\mathbf{v}}_k={\left[\stackrel{{\mathbf{v}}_{k,1}^{\mathrm{T}}}{\overbrace{v_{k,C_{\mathrm{v}1}} \cdots v_{k,C_{\mathrm{v}5}}}},\stackrel{{\mathbf{v}}_{k,2}^{\mathrm{T}}}{\overbrace{v_{k,n_{\mathrm{r}, \lambda 1}}\cdots v_{k,n_{\mathrm{r}, \lambda 7}}}},\stackrel{{\mathbf{v}}_{k,3}^{\mathrm{T}}}{\overbrace{v_{k,n_{i, \lambda 1}}\cdots v_{k,n_{i, \lambda 7}}}},\stackrel{{\mathbf{v}}_{k,4}}{\overbrace{v_{k,h_a}}}, \stackrel{{\mathbf{v}}_{k,5}}{\overbrace{v_{k,s_a}}}, \stackrel{{\mathbf{v}}_{k,6}}{\overbrace{v_{k,f_{\mathrm{v}, \mathrm{sphere}}}}}\right]}^{\mathrm{T}}$$

(122)

$$\mathbf{w}=\left[{\mathbf{w}}_1 {\mathbf{w}}_2 \dots {\mathbf{w}}_{N_{\mathrm{pixel}}}\right], \mathrm{with} {\mathbf{w}}_p=\left[\begin{array}{cccc}{w}_p(1);& w_p(2);& \cdots ;& w_p(N_{\mathrm{PC}})\end{array}\right]$$

(123)

$${\mathbf{x}}_{p,\mathrm{uncorr}}=\left[{\mathbf{x}}_{p,\mathrm{uncorr},1};{\mathbf{x}}_{p,\mathrm{uncorr},2};\cdots ;{\mathbf{x}}_{p,\mathrm{uncorr},N_{\mathrm{TP},\mathrm{uncorr}}}\right]$$

(124)

with $N_{\mathrm{TP},\mathrm{uncorr}}$= 6 types of uncorrelated parameters, including the BRDF spectral weight $a_{\lambda }$ (j = 1), anisotropy parameter $k_{\lambda }$ (j = 2), anisotropy parameter $g_{\lambda }$ (j = 3), pBRDF weight ${ϵ}_{\lambda }$ (j = 4), shadowing width $k_{\gamma }$ (j = 5), and slope variance ${\sigma }_s$ (j = 6), namely,

$${\mathbf{x}}_{p,\mathrm{uncorr}, j=1}=\left[a_p({\lambda }_1);a_p({\lambda }_2);\cdots ; a_p({\lambda }_7)\right]$$

(125)

$${\mathbf{x}}_{p,\mathrm{uncorr}, j=2}=\left[k_p({\lambda }_1);k_p({\lambda }_2);\cdots ; k_p({\lambda }_7)\right]$$

(126)

$${\mathbf{x}}_{p,\mathrm{uncorr}, j=3}=\left[g_p({\lambda }_1);g_p({\lambda }_2);\cdots ; g_p({\lambda }_7)\right]$$

(127)

$${\mathbf{x}}_{p,\mathrm{uncorr}, j=4}=\left[{\epsilon }_p({\lambda }_1);{\epsilon }_p({\lambda }_2);\cdots ; {\epsilon }_p({\lambda }_7)\right]$$

(128)

$${\mathbf{x}}_{p,\mathrm{uncorr}, j=5}=\left[k_{\gamma ,p}\right]$$

(129)

$${\mathbf{x}}_{p,\mathrm{uncorr}, j=6}=\left[{\sigma }_{\mathrm{s},p}^2\right]$$

(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 [19]. 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 [6,19], 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.

5.1. 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⁴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 [19]. To control the strength of multiple types of constraints, the initial values of Lagrange multipliers for the two PCs used in our retrieval are provided in

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 Difference for Spectral Smoothness Constraints $\left({\mathit{m}}_{\mathbf{uncorr},\diamond }\right)$	Lagrange Regularization Factor $\left({\mathit{\gamma }}_{\mathbf{uncorr},\diamond }\right)$
Aerosol parameters (scene-averaged)
Volume concentration of size components (Cv, 1-5, μm3/μm2)	[1.0 × 10−6,5]	0.002	-	-
Central height of aerosol distribution profile (ha, km)	[0.05,10]	1	-	-
Standard deviation of aerosol distribution profile (sa)	[0.5,2.5]	0.75	-	-
Real part of refractive index (nr(λ))	[1.33,1.60]	1.50	1	0.1
Imaginary part of refractive index (ni(λ))	[5.0 × 10−7,5.0 × 10−1]	0.005	2	0.01
Spherical particle volume fraction (fv,sphere)	[0.5,1.0]	0.95	-	-
Surface parameters (pixel-resolved)
BRDF spectral weight (aλ)	[0,0.7]	0.015–0.1	3	0.1
Anisotropy parameter (kλ)	[0,1]	0.6	1	0.5
Anisotropy parameter (gλ)	[−1,1]	0.1	1	0.5
pBRDF weight (ελ)	[0,10]	0.01	-	-
Shadowing width (kγ)	[0,1]	0.75	1	0.1
Slope variance (σs )	[0.05,0.5]	0.075	-	-

,

Table 4. Initial guess and the order of difference and Lagrange multipliers for imposing within-pixel and across-pixel smoothness constraints on the correlated aerosol fields through the first two PCs.

	Initial Guess (PC 1)	Range (PC 1)	Initial Guess (PC 2)	Range (PC 2)	Order of Finite Difference for Spectral Smoothness Constraints $\left({\mathit{m}}_{\mathbf{corr},\diamond }^{\mathbf{wv}}\right)$	Lagrange Multiplier $\left({\mathit{\gamma }}_{\mathbf{corr},\diamond }^{\mathbf{wv}}\right)$	Order of Finite Difference for Spatial Smoothness Constraints $\left({\mathit{m}}_{\mathbf{corr},\mathbf{\Delta }}^{\mathbf{wv}}\right)$	Lagrange Multiplier $\left({\mathit{\gamma }}_{\mathbf{corr},\mathbf{\Delta }}^{\mathbf{wv}}\right)$
log(Cv, 1-5)	0.1–0.6	[−0.75,+0.75]	−3 × 10−1–7 × 10−1	[−1,+1]	-	-	1	1
log(ha)	≈−1 × 10−1	[−0.9,+0.9]	≈5 × 10−1	[−1,+1]	-	-	1	0.01
log(sa)	≈−2 × 10−2	[−0.4,+0.4]	≈7 × 10−2	[−0.4,+0.4]	-	-	1	0.01
log(nr(λ))	≈3 × 10−3	[−0.1,+0.1]	≈5 × 10−3	[−0.1,+0.1]	1	0.1	1	10
log(ni(λ))	≈−2 × 10−2	[−0.1,+0.1]	≈2 × 10−2	[−0.1,+0.1]	2	0.01	1	1
log(fv,sphere)	≈6 × 10−3	[−0.05,+0.05]	≈4 × 10−3	[−0.05,+0.05]	-	-	1	0.1

and

Table 5. The order of difference and Lagrange multipliers for imposing within-PC constraints on the first two PC vectors and for imposing across-pixel constraint on the PC weights.

	Initial Guess	Range (in log-space)	Order of Finite Difference for Spectral Smoothness Constraints on First (Second) PC Vectors $\left({\mathit{m}}_{\mathbf{c}\mathbf{o}\mathbf{r}\mathbf{r},\diamond }^{\mathbf{v}}\right)$	Lagrange Multiplier on First (Second) PC Vectors $\left({\mathit{\gamma }}_{\mathbf{c}\mathbf{o}\mathbf{r}{r},\diamond }^{\mathbf{v}}\right)$	Order of Finite Difference for Spatial Smoothness Constraints on First (Second) PC weight $\left({\mathit{m}}_{\mathbf{c}\mathbf{o}\mathbf{r}\mathbf{r},\mathbf{\Delta }}^{\mathbf{w}}\right)$	Lagrange Multiplier on PC Weights $\left({\mathit{\gamma }}_{\mathbf{c}\mathbf{o}\mathbf{r}\mathbf{r},\mathbf{\Delta }}^{\mathbf{w}}\right)$
log(Cv, 1-5)	−0.3–0.7	[−1,+1]	-	-	-	-
log(ha)	≈5 × 10−1	[−1,+1]	-	-	-	-
log(sa)	≈7 × 10−2	[−0.4,+0.4]	-	-	-	-
log(nr(λ))	≈5 × 10−3	[−0.1,+0.1]	1(1)	0.01(0.001)	-	-
log(ni(λ))	≈2.5 × 10−2	[−0.01,+0.01]	2(2)	0.001(0.00001)	-	-
log(fv,sphere)	≈4 × 10−3	[−0.05,+0.05]	-	-	1(1)	0.05(0.005)
wp	0	[−10,10]	-	-	1(1)	0.05(0.005)

. Note that the difference between Table 3 and Table 4 and Table 5 is that Table 3 and Table 4 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. Table 3, Table 4 and Table 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.

5.2. 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 first 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 [34]. 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³/μm²) from the correlated multi-pixel inversion with those retrieved using original multi-pixel algorithm adapted for AirMSPI [19]. 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². The mean absolute difference (MAD) in AERONET and AirMSPI results is also calculated to measure the overall deviation. Strong correlation and low bias (R² ≥ 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 ≤ ${\mathrm{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 (${\mathrm{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.

Figure 10 and Figure 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 Figure 10 and Figure 11. Under the circumstance of no AERONET reference data before or after AirMSPI measurement, only AERONET measurement/retrieval error is reported. Figure 12 and Figure 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 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 [19]. For AOD, linear regression is performed. Values of regression related parameters (a, b, R²) 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 Figure 11, Figure 12 and Figure 13 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. 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 [19]. 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²) 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² is reported only for AOD in the table.

Parameters	Correlated Multi-pixel Inversion		Original Multi-pixel Inversion
	MAD	R2	MAD	R2
AOD355 nm	0.058	0.917	0.040	0.925
AOD445 nm	0.035	0.927	0.024	0.955
AOD470 nm	0.030	0.933	0.021	0.960
AOD555 nm	0.020	0.933	0.016	0.960
AOD660 nm	0.016	0.922	0.013	0.959
AOD865 nm	0.015	0.815	0.014	0.851
SSA445 nm	0.035	-	0.030	-
SSA555 nm	0.036	-	0.030	-
SSA660 nm	0.040	-	0.032	-
SSA865 nm	0.041	-	0.035	-
nr,445 nm	0.052	-	0.039	-
nr,555 nm	0.051	-	0.039	-
nr,660 nm	0.046	-	0.036	-
nr,865 nm	0.038	-	0.037	-
ni,445 nm	0.004	-	0.004	-
ni,555 nm	0.004	-	0.004	-
ni,660 nm	0.005	-	0.004	-
ni,865 nm	0.005	-	0.005	-
reff,fine (μm)	0.024	-	0.022	-
reff,coarse (μm)	1.050	-	0.993	-

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 [19]. These properties include AOD, SSA, real and imaginary parts of refractive index ($n_{\mathrm{r}}$, $n_{\mathrm{i}}$), and effective radii of fine and coarse particles ($r_{\mathrm{eff}, \mathrm{fine}}$, $r_{\mathrm{eff}, \mathrm{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_{\mathrm{r}}$, $n_{\mathrm{i}}$, $r_{\mathrm{eff}, \mathrm{fine}}$, and $r_{\mathrm{eff}, \mathrm{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 [19]) and those estimated for AERONET, as observed in Figure 10, Figure 11, Figure 12, Figure 13 and Figure 14 of this work, as well as Figure 4, Figure 5, Figure 6, Figure 7 and Figure 8 in the previous study [19].

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.

6. 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 [13,14], and the near-real time POLDER aerosol retrieval [35]. 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 [36,37], the ASTER spectral library [38], and the U.S. Geological Survey (USGS) digital spectral library [39]. 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) [13,14] and Moderate resolution Imaging Spectroradiometer (MODIS) [40], 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 [41].