Aerosol Retrieval Study from a Particulate Observing Scanning Polarimeter Onboard Gao-Fen 5B without Prior Surface Knowledge, Based on the Optimal Estimation Method

Abstract

To meet the demand for the aerosol detection of single-angle and multi-band polarization instrument containing short-wave infrared bands, an inversion algorithm that makes full use of multi-band intensity and polarization information is proposed based on optimal estimation theory. This method uses the polarization information in the short-wave infrared band to perform surface and atmosphere decoupling without a prior information on the surface. This obtains the initial value of the aerosol, and then it uses the scalar information to obtain the final result. Moreover, the multi-band information of the instrument is used for decoupling the surface and atmospheric information, which avoids the inversion error caused by the untimely update of the surface reflectance database and the error of spatio-temporal matching. The measured data of the Particulate Observing Scanning Polarimeter (POSP) are used to test the proposed algorithm. Firstly, to verify the effectiveness of the algorithm under different surface conditions, four regions with large geographical differences (Beijing, Hefei, Baotou, and Taiwan) are selected for aerosol optical depth (AOD) inversion, and they are compared with the aerosol robotic network (AERONET) products of the nearby stations. The validation against the AERONET products produces high correlation coefficients of 0.982, 0.986, 0.718, and 0.989, respectively, which verifies the effectiveness of the algorithm in different regions. Further, we analyzed the effectiveness of the proposed algorithm under different pollution conditions. Regions with AOD >0.7 and AOD < 0.7 are screened by using the AOD products of the Moderate-Resolution Imaging Spectroradiomete (MODIS), and the AOD of the corresponding region is inverted using POSP data. It was found to be spatially consistent with the MODIS products. The correlation coefficient and root mean square error (RMSE) in the AOD high region were 0.802 and 0.217, respectively, and 0.944 and 0.022 in the AOD low region, respectively, which verified the effectiveness of the proposed algorithm under different pollution conditions.

1. Introduction

Atmospheric aerosol, solid or liquid particles suspended in the atmosphere, including those from natural and anthropogenic sources, are an important component of the Earth’s atmosphere and they play an important role in the climate system and atmospheric environment system. In the climate system, aerosol not only affects the radiative energy balance by interacting with solar and terrestrial radiation [1,2,3], but they also change the characteristics and lifetimes of clouds [4,5,6,7,8]. In terms of the environment, aerosols are considered to be the main cause of air pollution, and they have a significant negative impact on human health [9,10,11,12,13,14,15]. Therefore, it is of great significance to carry out the remote sensing retrieval of aerosol and particulate matter concentration.

According to different remote sensing platforms, aerosol remote sensing can be divided into ground-based, airborne, and satellite remote sensing. Satellite remote sensing has the advantages of timeliness, periodicity, and large scale, and it is an effective method for detecting the optical properties and spatial distributions of aerosols [16,17,18,19,20]. The inversion of aerosol properties using satellite data has been carried out for more than 40 years. According to the difference of remote sensing data and the method of surface-atmosphere decoupling, the following methods are mainly used: dark target method using the statistical relationship between the apparent reflectance of short-wave infrared channels and the surface reflectance of red and blue channels [21,22,23]; deep blue algorithm for bright surface area [24,25]; in the Multi-angle Imaging SpectroRadiometer (MISR) operational aerosol retrieval algorithm, a mathematical method, empirical orthogonal functions, is used to obtain the surface contribution, which gives a general solution for the surface reflectance estimation in multi-angle aerosol remote sensing [26,27,28]; Inversion of aerosol optical thickness using polarization information [29,30,31,32,33,34].

Although for the aerosol inversion, lookup table (LUT)-based methods are widely applied, the method turned out to be incapable of fully exploiting the information implicit in the satellite measurements [35]. So, the optimal estimation theory is introduced into atmospheric research. The optimal estimation method was applied early to ground-based observations [36]. After that Dubovik et al. developed the generalized retrieval of an aerosol and surface properties algorithm which can be applied to a multi-angle observation satellite such as Polarization and Directionality of the Earth’s Reflectance (POLDER) [37]. There are some inversion framework based on optimal estimation theory that have been established [38,39,40,41]. Nevertheless, further studies based on optimal estimation theory need to be developed to apply to different satellite measurements such as single-angle multi-band polarization measurements, including short-wave infrared band.

The Particulate Observing Scanning Polarimeter (POSP) on board the Gaofen-5B satellite obtains atmospheric information in three dimensions: spectrum, observation angle and polarization through time synchronization and spatial coverage, and it obtains atmospheric parameters over the imaging area. The unique advantage of POSP is that it has nine detection channels from 380 to 2250 nm, all of which have polarization detection capabilities. Although its observation mode is single-angle observation, its advantages of multiple detection bands, wide band range and polarization detection capability will help with the inversion of aerosol optical depth (AOD).

According to the characteristics of POSP instruments, the intensity polarization joint inversion algorithm is proposed on the basis of optimal inversion theory, which makes full use of multi-band information. Firstly, the polarization reflectance in the 2250 nm band is used as the initial value of the surface polarization reflectance, and the optimal aerosol column concentration is obtained via the optimization iterative method. Then, according to the preliminary obtained AOD, the scalar signal is subjected to atmospheric correction to obtain the initial value of surface reflectance. Scalar information is used for the optimization iteration to obtain the final AOD.

This paper is divided into four parts: firstly, the characteristics of POSP data and the theory of optimization inversion framework are introduced. Secondly, the implementation details of the intensity polarization joint inversion algorithm are introduced in detail. Then, the effectiveness of the algorithm is tested in different regions and under different pollution conditions. Finally, the relevant conclusions are given.

2. Data and Optimization Estimate Framework

2.1. POSP Data Introduction

POSP is a high-precision polarization scanner developed by Anhui Institute of Optics and Fine Mechanics, Hefei Institute of Physical Sciences, Chinese Academy of Sciences, which is an onboard satellite for earth observation and can obtain multispectral and polarized radiation information on the surface and in the atmosphere. POSP adopts an on-board calibration to ensure the measurement accuracy and accuracy of the POSP operation in orbit [42]. POSP provides the first three components (I, Q, and U) of the Stokes vector corresponding to the band, and the basic parameters are shown in

Table 1. The basic parameters of POSP.

Parameter	Value
Central wavelength/nm	380, 410, 443, 490, 670, 865, 1380, 1610, 2250
Bandwidth/nm	20, 20, 20, 20, 20, 40, 40, 60, 80
Stokes parameters	I, Q, U
Quantized digit	14 bit
Radiation calibration error	≤5%
Polarization calibration error	≤0.5%

.

Since the 380 nm band is greatly influenced by atmospheric molecules, the 1380 nm band is mainly used for cirrus detection. So, seven bands except for 380 nm and 1380 nm are used for AOD inversion in this paper.

2.2. Optimization Estimate Framework

The optimal inversion framework in this paper is based on the previous information analysis of aerosol [43] plus the solution of the cost function. Let x be a state vector composed of n parameters (such as aerosol volume, refractive index, particle size distribution, etc.). The state vector will be retrieved from the observation vector y containing m measured elements, such as spectral scalar and polarization reflectivity. The relationship between y and x is shown in Formula (1):

$$y=F\left(x\right)+\in$$

(1)

where the forward model F describes the forward model of the physical relationship between y and x. ∈ represents the combined error of observation uncertainty and forward model uncertainty.

For satellite remote sensing, the forward model F is a nonlinear function, and the optimized solution of Equation (1) needs to go through multiple Gaussian hypothesis processes. According to the maximum likelihood method, the state vector of the t’th iteration corresponds to the solution that minimizes the quadratic cost function as follows:

$$J\left(x^t\right)=\frac{1}{2}{[y-F\left(x^t\right)]}^T{S}_{ϵ}^{-1}\left[y-F\left(x^t\right)\right]+\frac{1}{2}{r}_a{(x^t-x_a)}^T{S}_a^{-1}\left(x^t-x_a\right)$$

(2)

In Formula (2), ${\mathrm{x}}^{\mathrm{t}}$ represents the state vector at the t iteration, the two terms on the right side represent the constraints on the inversion results from observation and a prior estimation, respectively. The observation error covariance ${\mathrm{S}}_{ϵ}$ represents the uncertainty of the model and measurements. The prior error covariance ${\mathrm{S}}_{\mathrm{a}}$ represents the uncertainty of the prior estimate ${\mathrm{x}}_{\mathrm{a}}$.${\mathrm{r}}_{\mathrm{a}}=\mathrm{n}/2$, where n is the number of parameters to be inverted [44].

The process of finding the minimum value of the cost function needs to go through multiple iterations, using the LM (Levenberg-Marquardt) iterative method [45]. The state vector at the (t+1)-th iteration can be expressed as:

$$x^{t+1}=x^t-{[\left(r_a+\gamma \right)S_a^{-1}+K_t^T{S}_{ϵ}^{-1}{K}_t]}^{-1}\left\{K_t^T{S}_{ϵ}^{-1}\left[F\left(x^t\right)-y\right]+S_a^{-1}\left[x^t-x_a\right]\right\}$$

(3)

In Formula (3), $\mathsf{\gamma }$ is the LM parameter, which is initially set to 10, and is adjusted according to Rodgers’ update strategy after each iteration of x. K is called the weighting function matrix or Jacobian matrix, and it consists of the partial derivative of each forward model element with respect to each state vector element, i.e., ${\mathrm{K}}_{\mathrm{i},\mathrm{j}}=\partial {\mathrm{F}}_{\mathrm{i}}/\partial {\mathrm{x}}_{\mathrm{j}},\left(\mathrm{i}=1,\dots \mathrm{m};\mathrm{j}=1,\dots \mathrm{n}\right)$(K_i,j denotes the element in the i-th row and j-th column of matrix K, F_i is the i-th simulated measurement, and x_j is the j-th element in the state vector x).

3. Optimal Estimation Inversion Algorithm

3.1. A Priori Information on Aerosols and Surface

In the framework of optimal inversion, a priori estimation is an important factor affecting the inversion results. Reasonable priori estimation can reduce the number of inversion iterations, improve the time efficiency of the algorithm and avoid unreasonable convergence results. A priori estimation mainly comes from the a priori knowledge of state vectors and non-state vectors, such as historical meteorological data, ground measured data, etc. The influence of prior estimation on inversion is manifested in two aspects: first, the a priori estimation of the state parameter is input into the inversion program as the initial value for iteration, and the prior term that constitutes the cost function constrains the final inversion result. Secondly, the reasonable values of non-state parameters can reduce the model error and improve the accuracy of the inversion results.

3.1.1. Aerosol Model

In many studies, the bimodal lognormal distribution function is commonly used to characterize the aerosol particle size distribution, as follows:

$$n\left(r\right)={\displaystyle \sum }_{i=1}^2\frac{C_i}{r{\sigma }_i\sqrt{2\pi }}exp\left(-\frac{{\left(lnr-ln{r}_{m,i}\right)}^2}{2{\sigma }_i^2}\right)$$

(4)

In Formula (4), C denotes the volume concentration, σ is the standard deviation, and $r_m$ is the mean radius. The bimodal log normal distribution requires six parameters to represent.

Under the assumption of a bimodal log-normal distribution, five pairs of parameters can be used to describe the aerosol model: the particle column concentration $v_0^f$ and $v_0^c$, the effective radius $r_{eff}^f$ and $r_{eff}^c$, the effective variance $v_{eff}^f$ and $v_{eff}^c$, the real part of the refractive index $m_r^f$ and $m_r^c$, and the imaginary part of the refractive index $m_i^f$ and $m_i^c$. The superscripts f and c represent fine mode particles and coarse mode particles, respectively. A priori information about the aerosol models can use the results of previous studies [46], and the aerosol particle size distribution is shown in Figure 1.

The legends in Figure 1 represent different types of aerosol models. LH-LP: Low Humidity-Low Pollution Aerosol, CBC: Carbonaceous-Black Carbon Aerosol, DC: Dust Urban Aerosol, DD: Sand Dust-Desert Aerosol, LH-HP: Low Humidity-High Pollution type aerosol, HH-LP: High Humidity-Low Pollution type aerosol, HH-Hp: High Humidity-High Pollution type aerosol, CBrC: Carbonaceous-Brown Carbon type aerosol.

3.1.2. Processing of A Priori Information about Surface

At present, the methods for obtaining a priori surface information mainly include the surface model method [29], the prior surface reflectivity database [24], etc. In this paper, according to the polarization reflectancenot being sensitive to band changes [47,48], the polarization reflectance of the shortwave infrared band is used as the initial value of the surface polarization reflections of other bands, and then the aerosol parameters to be retrieved are obtained according to an iterative strategy. According to the aerosol parameters obtained via polarization inversion, the atmospheric transmittance, hemispheric albedo and atmospheric path radiation are obtained using Unified Linearized Vector Radiative Transfer Model (UNL-VRTM) [47], and then the empirical value of the surface reflectance is obtained using Formula (5):

$${\rho }_s=\frac{\rho -{\rho }_0}{\left(\rho -{\rho }_0\right)\ast S+T}$$

(5)

In Formula (5), $\rho$, ${\rho }_0$, and ${\rho }_s$ are the normalized reflectances of the atmospheric top, atmospheric path radiation, and surface, respectively. S and T are the atmospheric hemisphere albedo and the total atmospheric transmittance, respectively.

3.2. Intensity Polarization Joint Optimization Inversion Algorithm

3.2.1. Satellite Observation Model

The radiation signals obtained by satellites can be characterized by the apparent reflectance. There are multiple interaction processes between the surface and the atmosphere, which are complicated by different surface reflectances and different atmospheric conditions. Therefore, the apparent reflectance model is as follows [49]:

$$\begin{gathered} {\mathrm{R}}^{\mathrm{TOA}}\left({\mathsf{\mu }}_{\mathrm{s}},{\mathsf{\mu }}_{\mathrm{v}},\mathsf{\phi }\right)={\mathrm{R}}^{\mathrm{Atm}}\left({\mathsf{\mu }}_{\mathrm{s}},{\mathsf{\mu }}_{\mathrm{v}},\mathsf{\phi }\right)+\frac{{\mathrm{T}}_{\mathrm{s}}{\mathrm{T}}_{\mathrm{v}}\mathsf{\rho }\left({\mathsf{\mu }}_{\mathrm{s}},{\mathsf{\mu }}_{\mathrm{v}},\mathsf{\phi }\right)}{1-\mathrm{s}\mathsf{\rho }\left({\mathsf{\mu }}_{\mathrm{s}},{\mathsf{\mu }}_{\mathrm{v}},\mathsf{\phi }\right)} \\ {\mathrm{R}}^{\mathrm{TOA}}\left({\mathsf{\mu }}_{\mathrm{s}},{\mathsf{\mu }}_{\mathrm{v}},\mathsf{\phi }\right)=\frac{\mathsf{\pi }\mathrm{I}\left({\mathsf{\mu }}_{\mathrm{s}},{\mathsf{\mu }}_{\mathrm{v}},\mathsf{\phi }\right)}{{\mathsf{\mu }}_{\mathrm{s}}{\mathrm{E}}_0} \end{gathered}$$

(6)

In Formula (6), ${\mathrm{R}}^{\mathrm{TOA}}$ and ${\mathrm{R}}^{\mathrm{Atm}}$ are the top-of-atmosphere (TOA) reflectance and the atmospheric path radiation, ${\mathrm{T}}_{\mathrm{s}}$ and ${\mathrm{T}}_{\mathrm{v}}$ are atmospheric downward and upward transmittances, respectively, $\mathsf{\rho }$ is the surface albedo and $\mathrm{s}$ is the atmospheric hemispheric albedo. ${\mathsf{\mu }}_{\mathrm{s}},{\mathsf{\mu }}_{\mathrm{v}}$, and $\mathsf{\phi }$ are the cosine of the solar zenith angle, the cosine of the observed zenith angle, and the relative azimuth angle. ${\mathrm{E}}_0$ is the solar constant (the solar spectral irradiance outside the atmosphere). The surface albedo can be obtained if the apparent reflectance, atmospheric path radiation, atmospheric transmittance, and atmospheric hemispheric albedo are known.

The polarization signal obtained by the satellite is characterized by the TOA polarized reflectance. The polarization contribution caused by the multiple interactions between the surface and the atmosphere is very small. After the attenuation of the atmosphere, the contribution to the observation at the top of the atmosphere can be ignored. Therefore, the TOA polarized reflectance modeling only needs to consider the upward polarization scattering of the atmosphere and the single reflection of the surface [50], and the TOA polarized reflectance can be written as follows:

$$\begin{gathered} {\mathrm{R}}_{\mathrm{P}}^{\mathrm{TOA}}\left({\mathsf{\mu }}_{\mathrm{s}},{\mathsf{\mu }}_{\mathrm{v}},\mathsf{\phi }\right)={\mathrm{R}}_{\mathrm{P}}^{\mathrm{Atm}}\left({\mathsf{\mu }}_{\mathrm{s}},{\mathsf{\mu }}_{\mathrm{v}},\mathsf{\phi }\right)+{\mathrm{T}}_{\mathrm{s}}{\mathsf{\rho }}_{\mathrm{P}}\left({\mathsf{\mu }}_{\mathrm{s}},{\mathsf{\mu }}_{\mathrm{v}},\mathsf{\phi }\right){\mathrm{T}}_{\mathrm{v}} \\ {\mathrm{R}}_{\mathrm{P}}^{\mathrm{TOA}}\left({\mathsf{\mu }}_{\mathrm{s}},{\mathsf{\mu }}_{\mathrm{v}},\mathsf{\phi }\right)=\frac{{\mathsf{\pi }\mathrm{I}}_{\mathrm{p}}\left({\mathsf{\mu }}_{\mathrm{s}},{\mathsf{\mu }}_{\mathrm{v}},\mathsf{\phi }\right)}{{\mathsf{\mu }}_{\mathrm{s}}{\mathrm{E}}_0} \end{gathered}$$

(7)

In Formula (7), ${\mathrm{R}}_{\mathrm{P}}^{\mathrm{TOA}}$ is the TOA polarized reflectance, ${\mathrm{R}}_{\mathrm{P}}^{\mathrm{Atm}}$ is the polarized reflectance of the atmosphere (which can be calculated by setting the surface polarized reflectance to zero), and ${\mathsf{\rho }}_{\mathrm{P}}$ is the surface polarized reflectance. ${\mathrm{I}}_{\mathrm{p}}=\sqrt{{\mathrm{Q}}^2+{\mathrm{U}}^2}$ is the polarized radiance at the top of the atmosphere.

In summary, when using the intensity information and polarization information for the optimization solution, due to different observation models and different assumptions of the surface, the state vector x and observation vector y are shown in

Table 2. State vector x and observation vector y settings.

Category	Polarization Setting	Intensity Setting
Observation vector	$y={\left[I_{p_{{\lambda }_1}},\dots ,I_{p_{{\lambda }_6}}\right]}^T$	$y={\left[I_{{\lambda }_1},\dots ,I_{{\lambda }_7}\right]}^T$
State vector	$x={\left[V_f,V_c,{\rho }_p\right]}^T$	$x={\left[V_f,V_c,{\rho }_{{\lambda }_1},\dots ,{\rho }_{{\lambda }_7}\right]}^T$

:

In Table 2, the superscript T represents the transpose of the vector, and ${\lambda }_1-{\lambda }_7$ are the corresponding center wavelengths of POSP: 410 nm, 443 nm, 490 nm, 670 nm, 865 nm, 1610 nm and 2250 nm, respectively. $V_f$ and $V_c$ are the aerosol fine-mode and coarse-mode column concentrations, respectively, and $\mathsf{\rho }$ represents the surface reflectance. Since it is assumed that the surface polarized reflectance is independent of the wavelength, and the TOA polarized reflectance in the 2250 nm band is used as the initial value of the surface polarized reflectance, in the polarization inversion, the state vector x has only three parameters, and the observation vector has 6 parameters. In the intensity inversion, the state vector has nine parameters and the observation vector has seven parameters.

3.2.2. Cost Function Solution Method

The optimal estimation theory provides the cost function J(x) for finding the optimal inversion result, and gives the specific form. J(x) is a nonlinear function, and Newton’s iterative method is a common method for solving nonlinear functions. In this paper, the LM algorithm is used to solve the minimum value of J(x), and in the LM program module, the forward radiated transmission simulation module and the calculation module of the cost function are combined to establish the program framework for optimal iterative inversion. Under the framework of the program, the optimal estimation of the state vector is found after several iterations. The convergence condition is as follows:

$$\Delta =\frac{{\left|\left|{\mathrm{x}}^{\mathrm{t}+1}-{\mathrm{x}}^{\mathrm{t}}\right|\right|}_2}{{\left|\left|{\mathrm{x}}^{\mathrm{t}}\right|\right|}_2}\le \mathsf{\epsilon }$$

(8)

In Formula (8), where || ||2 is the L2 norm, $\mathsf{\epsilon }$ is the convergence threshold for judging the end of the iteration, which is an empirical value, and is taken as 0.1% in this paper.

In the actual operation process, when the number of iterations exceeds the set maximum number of iterations of six, or when the element value in x exceeds the boundary range, the iteration stops and an invalid value is output.

3.2.3. Algorithm Implementation

The LM program module, UNL-VRTM, and the cost function calculation module are combined to make full use of the multi-band polarization information of POSP to carry out the optimal inversion of AOD. The detailed flow chart is shown in Figure 2. After reading the POSP data, the proposed algorithm can be divided into the following two stages.

In the first stage, the optimal inversion is performed using the polarization information to obtain the aerosol column concentration. There are three main steps at this stage: (1) Obtain the initial value of the surface polarization reflectivity of other bands through the apparent polarization reflectivity of 2250 nm; (2) Use the radiation transfer software UNL-VRTM to calculate the Jacobian matrix of the state vector; (3) Obtain a new state vector according to the update strategy of Formula (3), and then after several iterations, when the iterative convergence threshold is met, the aerosol column concentration parameter is output.

In the second stage, the final aerosol column concentration parameters are obtained by optimizing the inversion using the intensity information. First, input the aerosol column concentration parameters obtained in the previous stage into UNL-VRTM and run it, and then perform simple atmospheric correction according to Equation (5) to obtain the initial values of surface reflectance in different bands. Then, according to (2)–(3) of the first stage, the column concentration parameters of the aerosol are finally obtained.

The aerosol volume column concentration parameters are obtained via inversion according to the above process, and the spectral AOD can be further calculated by combining with the aerosol model parameters, but this calculation process does not need to be carried out separately. In the iterative process, when using the UNL-VRTM software for forward simulation, AOD can be obtained at the same time, and the calculation result can be output when the last step satisfies the iterative convergence condition.

4. Result and Discussion

4.1. Algorithm Iteration Process

According to the prior information acquisition and retrieval strategy described above, the POSP data from November 2021 to July 2022 are selected for the retrieval of AOD. An inversion case of 8 June 2022, over the Beijing site is selected to illustrate the iteration and convergence processes of the algorithm, as shown in Figure 3.

Figure 3a,c shows the fitting between the simulated values of the model and the measured values of POSP with the increase in the number of iterations in the process of inversion using polarization information and intensity information, respectively. After the last iteration, the simulated values using polarization information agree well with the observed values in the first five bands, and there is only a slight gap in the 1610 nm band. The simulated values using intensity information are in good agreement with the observed values in all bands, indicating that in single-angle observation, the intensity contains more information about the atmosphere than the polarization measurement. In Figure 3b,d of the cost function varies with the number of iterations in the process of polarization inversion and intensity inversion, respectively. It is also reflected that the intensity contains more information than polarization from (b) and (d). When convergence, the cost function value of intensity information inversion is less than that of polarization information inversion. After five iterations, the cost function almost does not change with the increase in the number of iterations, so that in order to save the time of the algorithm, the maximum number of iterations can be set to 6.

4.2. Evaluation Index of the Inversion Result

In the evaluation of the inversion results, quantitative indexes such as correlation coefficient (R), root mean square error (RMSE), expected error (EE), and the percentage of except error (Scope) are used. Taking AOD as an example, the specific definitions of each index are as follows:

$$\mathrm{R}=\frac{\mathrm{cov}\left({\mathsf{\tau }}_{\mathrm{retrieval}},{\mathsf{\tau }}_{\mathrm{true}}\right)}{\mathrm{std}\left({\mathsf{\tau }}_{\mathrm{retrieval}}\right)\times \mathrm{std}\left({\mathsf{\tau }}_{\mathrm{true}}\right)}$$

(9)

$$\mathrm{RMSE}=\sqrt{\frac{1}{\mathrm{N}}{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{N}}{\left({\mathsf{\tau }}_{\mathrm{retrieval}}-{\mathsf{\tau }}_{\mathrm{true}}\right)}^2}$$

(10)

$$\Delta =\mathrm{EE}=\pm 0.15{\mathsf{\tau }}_{\mathrm{true}}\pm 0.05$$

(11)

$$\mathrm{Scope}={\mathrm{N}}_{\mathrm{valid}}/\mathrm{N}\times 100\%$$

(12)

In Formulas (9)–(12), ${\mathsf{\tau }}_{\mathrm{retrieval}}$ is the AOD result of inversion using POSP data, and ${\mathsf{\tau }}_{\mathrm{true}}$ is the data for verification. Cov() and std() represent the covariance and standard deviation calculations, respectively. Scope is the proportion of the amount of data within the range of ${\mathsf{\tau }}_{\mathrm{true}}+\Delta$ for the inversion AOD to the total data. N and ${\mathrm{N}}_{\mathrm{valid}}$ are the amount of data and the total amount of data that fall into the error range, respectively.

4.3. Validation against Ground-Based Data

The AERONET level 2.0 category is used for ground-based validation in this paper. AERONET is a ground-based network of sun photometers that provides high-quality aerosol measurements including AOD products [36,51].

The spectral AOD are retrieved using the proposed algorithm. With regard to validation, the AOD at 670 nm is selected. Correspondingly, the AOD at 675 nm in the AERONET products is used.

In order to verify the effectiveness of the inversion algorithms in different regions, four regions (Beijing, Hefei, Baotou, and Taiwan) with considerable geographical differences were selected. Additionally these four areas have AERONET ground-based long-term observation stations, namely: Beijing area (Beijing station), Baotou area (AOE_Baotou station), Hefei area (Hefei station) and Taiwan area (Chen-Kung_Univ station), the corresponding latitude and longitude are shown in

Table 3. AERONET data for validation.

AERONET Sites	Longitude	Latitude	Date Range
Beijing	116.3814	39.9769	2021.11–2022.7
AOE_Baotou	109.6288	40.8517	2021.11–2022.7
Hefei	117.1622	31.9047	2021.11–2022.7
Chen-Kung_Univ	120.2047	22.9934	2021.11–2022.7

.

The location and observation time of the AERONET sites are matched to extract the observation data corresponding to POSP from November 2021 to July 2022, as shown in Figure 4.

The validation against the corresponding AEROENT product is shown in Figure 5a–d are the verification results in Beijing, Baotou, Hefei and Taiwan, respectively.

It can be found from Figure 5 that the inversion results based on the algorithm proposed in this paper are basically consistent with AOD products of the AERONET sites, and the results are good. However, due to the influences of different surface types and aerosol types, the Scope is different. The order from largest to smallest is Hefei, Taiwan, Beijing, and Baotou. The reason for this may be that the surface vegetation in the Baotou area is sparse and is often considered to be a bright surface, while the Hefei site is mostly covered by evergreen vegetation, and the surface reflectance value is relatively small, so that the impact of the surface on the apparent reflectance is relatively smaller than the bright surface.

In addition, the aerosols in Baotou are mainly composed of coarse particles, so when polarization inversion is performed, the AOD obtained will have a large error, and then a large error will be introduced when the initial value of surface reflectance is provided based on the results of polarization inversion, resulting in relatively poor final inversion results.

4.4. The Validation against MODIS Products

In order to verify the effectiveness of the intensity polarization joint inversion algorithm under different pollution conditions, MODIS’s AOD products at 550 nm are first used to screen the time and longitude ranges of AOD high value (AOD > 0.7) and low value (AOD < 0.7), and then filter out the corresponding POSP observation data. The specific latitude and longitude ranges and times are shown in

Table 4. Data information used for validation.

AOD Range	Longitude Range	Latitude Range	Date
AOD>0.7	114.6–115.9	37.5–38.5	2022.6.9
AOD<0.7	116.6–118.2	31.57–33.46	2022.5.4

.

The comparison between the inversion results in the high-value area and the MODIS product is shown in Figure 6. The blank area in the figure indicates no inversion results, because they are cloud pixels or no data.

Comparing the (a) and (b) figures in Figure 6, it can be found that the inversion results of POSP are consistent with the MODIS product as a whole, but that the inversion value of POSP in some areas is smaller than the AOD value in MODIS (prismatic area in the figure), and that the AOD value of POSP is greater than that of MODIS in some areas (circular area in the figure). In order to further analyze the difference between the AOD inversion results of POSP and the MODIS products, the inverted AOD and MODIS products were statistically analyzed and linear regression analysis was performed, and the results are shown in Figure 7.

According to the statistical results of Figure 7, a total of 125 high-value AOD data are obtained, of which 78.4% of the data fall within the error range, and the correlation coefficient R is 0.802.

Figure 8 shows the comparison of the AOD inversion results of POSP in the low AOD value area with the MODIS product. Comparing (a) and (b), it can be found that the POSP inversion results are consistent with the MODIS product as a whole. In some areas, the inversion result of POSP is smaller than MODIS (the prismatic area in the figure). In order to further analyze the effectiveness of the inversion algorithm, the inversion results and MODIS products were statistically analyzed, and the results are shown in Figure 9.

As can be seen from Figure 9, a total of 670 data points are obtained, 99.6% of the data are within the error range, and the correlation coefficient R is 0.944, indicating that the inversion results are in good agreement with the MODIS.

In summary, comparing the inversion results of high-value AOD and low-value AOD, it can be found that the inversion results of low-value AOD are significantly better than those in high-value AOD regions. This may be due to the assumption that the 2250 nm band is not sensitive to the atmosphere when performing polarization inversion, so that the TOA polarization reflectance of 2250 nm is used as the surface polarization reflectance of other bands. When the AOD is low, this assumption may have little effect on the inversion results, but in the high AOD value region, due to the increased probability of collision between light and particles, the 2250 nm band will carry more aerosol information, which will have a certain impact on the inversion results.

5. Conclusions

In this paper, based on the observation characteristics of POSP single-angle multi-band polarization measurements, an intensity polarization joint optimization inversion algorithm based on optimization theory is proposed, which can effectively perform ground-air decoupling without prior knowledge of the surface. Additionally, the effectiveness of the algorithm is verified in different regions and under different pollution conditions. It is found that in the vegetation coverage area (Hefei) and the urban underlying surface (Beijing and Taiwan), the correlation coefficient between the inversion AOD at 670 nm from POSP and the ground-based AOD at 675 nm from AERONET can reach more than 0.98. Good results can also be obtained on the bright surface (Baotou), the correlation coefficient can reach 0.71, and the proportion of data within the error range can reach 85.5%. In addition, the inversion results of the high-value area and the low-value area of AOD are also analyzed. It is found that the AOD at a 550 nm correlation coefficient between POSP and MODIS in the high-value area is 0.802, and in the low-value area, it is 0.994, indicating that the performance of the algorithm in the low-value area is better than that of the high-value area.

This algorithm uses the insensitivity of the short-wave infrared band to the atmosphere to obtain the initial value of the surface polarized reflectance. Therefore, in the high-value area of aerosol or under the condition coarse mode dominated aerosol, it will have a certain impact on the inversion results.