Estimation of Ultrahigh Resolution PM2.5 Mass Concentrations Based on Mie Scattering Theory by Using Landsat8 OLI Images over Pearl River Delta

The aerosol optical depth (AOD), retrieved by satellites, has been widely used to estimate ground-level PM2.5 mass concentrations, due to its advantage of large-scale spatial continuity. However, it is difficult to obtain urban-scale pollution patterns from the coarse resolution retrieval results (e.g., 1 km, 3 km, or 10 km) at present, and little research has been conducted on PM2.5 mass concentration retrieval from high resolution remote sensing data. In this study, a physical model is proposed based on Mie scattering theory to evaluate the PM2.5 mass concentrations by using Landsat8 Operational Land Imager (OLI) images. First, the Second Simulation of the Satellite Signal in the Solar Spectrum (6S) model (which can simulate the transmission process of solar radiation in the Earth-atmosphere system and calculate the radiance at the top of the atmosphere) is used to build a lookup table to retrieve the AOD of the coast and blue bands based on the improved deep blue (DB) method. Then, the Angstrom formula is used to obtain the AOD of the green and red bands. Second, the dry near-surface AOD of four bands (coast, blue, green, red) is obtained through vertical correction and humidity correction. Third, aerosol particles are divided into four types based on the standard radiation atmosphere (SRA) model, and the optical properties of different aerosol types are analyzed to derive the volume distribution of aerosol particles. Finally, the relationship between the dry near-surface AOD of each band and the volume distribution of four aerosol particles is correlated, based on Mie scattering theory, and a physical model is established between the AOD and PM2.5 mass concentrations. Then, the distribution of PM2.5 mass concentrations is obtained. The retrieval results show that the distribution of AOD and PM2.5 at the urban scale in detail. The AOD results show that a reasonable relationship with a correlation coefficient (R2) of 0.66 and root mean square error (RMSE) of 0.1037 between Landsat8 OLI AOD and MODO4 DB AOD at 550 nm. The PM2.5 retrieval results are compared with the PM2.5 values measured by ground monitoring stations. The RMSEs for a certain day in different years, including 2017, 2018, 2019, and 2020, are 11.9470 μg/m³, 11.9787 μg/m³, 7.4217 μg/m³, and 5.4723 μg/m³, respectively. The total RMSE is 10.0224 μg/m³. The ultrahigh resolution PM2.5 results can provide pollution details at the urban scale and support better decisions on urban atmospheric environmental governance.

Delta (PRD) region was retrieved from Landsat8 Operational Land Imager (OLI) coast and blue spectral band data using the improved deep blue (DB) method [11], and then the Angstrom formula was used to obtain the AOD of the green and red bands. In the process, both the original Landsat8 OLI data and the surface reflectance data were from the United States Geological Survey (USGS) website. An aerosol model was determined by statistical analysis of long-term ground monitoring CE318 data. Second, the dry near-surface AOD of each band was obtained from RH and PBLH data by humidity correction and vertical correction. These data were from the China Environmental Monitoring Station (CEMS) and European Centre for Medium-Range Weather Forecasts (ECMWF). Third, to avoid directly obtaining FMF, the aerosol particles were divided into four types based on the standard radiation atmosphere (SRA) model, and the optical properties of different aerosol types were analyzed to derive the volume distribution of aerosol particles. Then, the relationship between the dry near-surface AOD of each band and the volume distribution of four aerosol particles was correlated based on Mie scattering theory, and a physical model was established between AOD and PM 2.5 mass concentrations. To verify our proposed method, the distribution of AOD and PM 2.5 is presented at the urban scale. The AOD results were compared with MODO4 DB AOD products. The PM 2.5 results were verified by ground monitoring sites. Finally, the performance and limitations of the model are also discussed.
This paper is organized as follows: Section 2 introduces materials and related methods; Section 3 mainly shows the distribution results of AOD and PM 2.5 mass concentrations; Section 4 discusses some advantages and limitations of this method by comparing it with other studies; and Section 5 is the conclusion of this paper.

Study Area and Datasets
The study area in this paper is part of the PRD, located in Guangdong Province of China, bordering the South China Sea (Figure 1 Left). The longitude of the study area is between 112 • 23 and 114 • 54 , and the latitude is between 21 • 49 and 24 • 09 . With the development of industrialization in the PRD over recent decades, regional air pollution has gradually increased and spread [35,36]. Therefore, remote sensing retrieval of PM 2.5 mass concentrations over the PRD region is significant.
As shown in Table 1, this paper used the following data: (a) Landsat8 OLI image data and reflectance data from the USGS (https://glovis.usgs.gov/, accessed on 9 October 2020), which have 11 bands (the spatial resolution of band 1-7, 9-11 is 30 m, and band 8 is a panchromatic band with 15m resolution) and a round trip period of about 16 days. The surface reflectance data was from atmospherically-corrected datasets from Landsat8 OLI sensor. Since only the visible bands of Landsat8 OLI were used to retrieve AOD in this paper, Table 2 mainly displays the detailed characteristics of these bands. All Landsat8 OLI image data have been preprocessed before beginning the AOD retrieval process. The preprocesses include geometric correction and radiometric calibration. The accuracy of the geometric correction is maintained at 1~2 pixels. We selected a certain day of the image with better data quality as experimental data in different years (including 2017-2020); (b) Hourly PM 2.5 site data and RH data from the weather stations of CEMS (http://data.cma.cn/, accessed on 25 October 2020). The distribution of weather stations in the study area is shown in Figure 1. The stations can provide average measured PM 2. 5 and RH values within one hour; thus, we select the hourly average PM 2.5 values close to imaging time as validation data. RH data is mainly used to study the moisture absorption characteristics of aerosol particles, which helps to obtain the extinction characteristics of dry aerosol particles, and then establish a relationship with near-surface PM 2.5 more accurately; (c) PBLH data from the ECMWF (http://www.ecmwf.int/ accessed on 18 November 2020) at a nominal spatial resolution (0.125 • × 0.125 • ), which can approximately represent the atmospheric height of the aerosol particle distribution over PRD region. We collected monthly PBLH data over the PRD region from 2017 to 2020 for vertical correction. All PBLH data have been re-projected and cropped, consistent with Landsat8 experimental data (Figure 1 Right);  17 40 E). The instrument provides long-term observations of aerosols over Guangzhou city over 2011. Multiyear observation results (January 2017 to December 2019) from Guangzhou city were used to obtain AOD and analyze the aerosol model, which can generate the lookup table suitable for the PRD region; (e) MODO4 DB AOD products derived from MODIS sensors onboard Terra and Aqua, are provided by DB algorithm at a coarse spatial resolution of 3 km × 3 km (Level-1 and Atmosphere Archive & Distribution System Distribution Active Archive Center: https://ladsweb.modaps.eosdis.nasa.gov/, accessed on 2 January 2021). Compared with the Dark Target (DT) algorithm, it can provide more accurate AOD results in urban areas of the PRD. In addition, the spatial coverage of MODIS DB retrieval can reach up to 100% over PRD region in cloudless conditions, which is suitable for the verification of Landsat8 OLI retrieval results; and (f) industrial zone and road data, which were obtained using crawler technology, from the Baidu map (https://map.baidu.com/, accessed on 2 January 2021), these geographic data are mainly used to analyze pollution characteristics in urban areas. All data (except CE318 data) was selected in accordance with the time to that Landsat8 satellite transits over the PRD region.

Methods
The workflow of our proposed method is shown in Figure 2. First, Landsat8 OLI data preprocessing was performed, and the input parameters of the Second Simulation of the Satellite Signal in the Solar Spectrum (6S) model were set to generate a lookup table for the coast and blue bands. Second, according to the generated lookup table, the surface reflectance and apparent reflectance data were brought into the radiation transmission model to obtain the AOD of the coast and blue bands, and then the Angstrom formula was used to derive the AOD of the green and red bands. Third, vertical correction and humidity correction were performed using PBLH and RH data based on the entire atmospheric AOD to obtain dry near-surface AOD. Fourth, a physical model between the dry near-surface AOD and PM2.5 mass concentrations was established based on Mie scattering theory. Fifth, the PM2.5 mass concentration was solved by using multiband dry near-surface AOD. Finally, result analysis and error estimation were performed on the AOD and PM2.5 retrieval results.

Methods
The workflow of our proposed method is shown in Figure 2. First, Landsat8 OLI data preprocessing was performed, and the input parameters of the Second Simulation of the Satellite Signal in the Solar Spectrum (6S) model were set to generate a lookup table for the coast and blue bands. Second, according to the generated lookup table, the surface reflectance and apparent reflectance data were brought into the radiation transmission model to obtain the AOD of the coast and blue bands, and then the Angstrom formula was used to derive the AOD of the green and red bands. Third, vertical correction and humidity correction were performed using PBLH and RH data based on the entire atmospheric AOD to obtain dry near-surface AOD. Fourth, a physical model between the dry nearsurface AOD and PM 2.5 mass concentrations was established based on Mie scattering theory. Fifth, the PM 2.5 mass concentration was solved by using multiband dry nearsurface AOD. Finally, result analysis and error estimation were performed on the AOD and PM 2.5 retrieval results.

AOD Retrieval Algorithm
This paper used radiative transfer theory to retrieve AOD from satellite remote sensing imagery. Spectral information received by the satellite was the combined effect of the extinction of all substances in the Earth's atmosphere and the reflection on the ground. The reflectance at the top of the atmosphere (TOA) received by the satellite sensor can be expressed by the following equation [37]: In the formula, θ is the solar zenith angle, is the view zenith angle, is the relative azimuth angle, λ is the corresponding sensor band, represents the TOA, representes the atmospheric path radiation reflectance, is the hemisphere albedo in the atmosphere, (θ ) and ( ) represent the atmospheric upward transmittance and downward transmittance, respectively, and represents the surface reflectance. To solve Equation (1), this study adopted the 6S model to simulate the transmission process of light in the atmosphere. By setting the imaging time, view geometry, specific band,

AOD Retrieval Algorithm
This paper used radiative transfer theory to retrieve AOD from satellite remote sensing imagery. Spectral information received by the satellite was the combined effect of the extinction of all substances in the Earth's atmosphere and the reflection on the ground. The reflectance at the top of the atmosphere (TOA) received by the satellite sensor can be expressed by the following equation [37]: In the formula, θ s is the solar zenith angle, θ v is the view zenith angle, φ is the relative azimuth angle, λ is the corresponding sensor band, ρ TOA λ represents the TOA, ρ α λ representes the atmospheric path radiation reflectance, S λ is the hemisphere albedo in the atmosphere, T(θ s ) and T(θ v ) represent the atmospheric upward transmittance and downward transmittance, respectively, and ρ S λ represents the surface reflectance. To solve Equation (1), this study adopted the 6S model to simulate the transmission process of light in the atmosphere. By setting the imaging time, view geometry, specific band, atmospheric model, aerosol type, and surface reflectivity, the model simulates the radiative transmission process in the atmosphere at the time input, and then ρ TOA λ under different AODs is simulated by three key parameters: ρ α λ , S λ and T λ (T(θ s )T(θ v )). Finally, the AOD can be obtained by comparing the real and simulated TOA.

Aerosol Model Determination
When using the 6S model to generate a lookup table, there were two important input parameters: The surface reflectance (the surface reflectance problem was solved by downloading it from the USGS) and the aerosol model. Related research shows that the aerosol model has a significant impact on the retrieval of AOD [38]. In this study, three years (2017-2019) of Sun photometer data from the CE318 instrument were used to determine the aerosol model. The site was located on the South Campus of Sun Yat-sen University at a longitude of 113 • 06 12 E and a latitude of 23 • 06 12 N. Through simple data processing and classification, AODs at 550 nm in different seasons (spring: March-May; summer: June-August; autumn: September-November; winter: December-February) were obtained. Figure 3a shows that the proportion of AOD between 0 and 0.3 in each season was relatively low throughout the year, while the proportion of AOD between 0.3 and 0.5 was highest in summer and lowest in spring. The proportion of AOD between 0.5 and 0.8 was relatively high and was approximately 30% in each season. The proportions of AOD in the intervals from 0.8-1.0 and 1-1.7 were highest in spring and lowest in summer, and the proportions of AOD in these intervals in summer were lower than the annual averages. This also reflected that the AOD in Guangzhou was higher in spring, reached a lower level in summer, and rebounded in autumn and winter. Then, aerosol particles were divided into four types (water-soluble, dust, marine, soot) based on the SRA model in this paper. To obtain the volume ratio of different aerosol particles, a Monte Carlo forward modeling method was adopted to determine the aerosol model by using processed AOD data in various bands [39,40], which mainly consisted of two steps. Firstly, we randomly generated the volume ratio of three particles were ω 1 , ω 2 , ω 3 (all values were between 0 and 1), and the volume ratio of the other particle was ω 4 = 1 − ω 1 − ω 2 − ω 3 (the value of ω 4 must be between 0 and 1), and then substituted them into the Equation (19) to calculate the simulated AOD of different bands. Secondly, a comparison was made between the simulated AOD and real AOD retrieved by the satellite. When their differences meet a certain accuracy range, ω 1 , ω 2 , ω 3 , and ω 4 can be considered as the correct volume ratio of aerosol particles, then the aerosol model will be finally determined. The distribution of average particle volumes in different seasons is shown in Figure 3b. The structural ratios of aerosol components in different seasons are shown in Table 3. Generally, the volumetric proportions of various components of aerosol particles were stable and did not change with the season, indicating a relatively fixed pollution source, and the particle volume distribution showed a bimodal structure that was mainly composed of water-soluble and soot particles with sizes less than 0.1 µm. According to the definitions of the SRA model, the aerosol model in the PRD is a mixed structure of urban industrial and ocean type. Therefore, when inputting the aerosol model parameters based on the 6S model, according to the transit season of the satellite image, an accurate AOD can be obtained by inputting the proportions of various components of aerosol particles corresponding to Table 3.

Multiband AOD Retrieval
After the aerosol model and surface reflectance were determined, the intervals and ranges of some input parameters (such as view geometry, AOD, atmospheric model, band) were set to generate a lookup table by using the 6S model. The details of these input parameters are listed in Table 4. The generated lookup table of the coast and blue band contains seven columns of values ( , , , θ , , , AOD). For the AOD retrieval of these two bands, the lookup table of the corresponding band was selected, and the corresponding parameters ( , , ) were derived by obtaining the view geometric information (θ , , ) of each pixel. Then, it was substituted into equation (1) to obtain the simulated TOA. The real AOD can be interpolated by comparing the real and simulated TOA. In this way, the coast and blue bands of AOD can be obtained accurately in this study, and then the green and red bands of AOD can be obtained by the Angstrom equation [41]: where τ(λ) is the AOD at different wavelengths, is the wavelength index, and represents the Angstrom turbidity coefficient. If the AOD of two wavelengths was obtained, and can be expressed as: = τ( ) = τ( ) Substituting and into Equation (2), the AOD of other wavelengths can be derived.

Multiband AOD Retrieval
After the aerosol model and surface reflectance were determined, the intervals and ranges of some input parameters (such as view geometry, AOD, atmospheric model, band) were set to generate a lookup table by using the 6S model. The details of these input parameters are listed in Table 4. The generated lookup table of the coast and blue band contains seven columns of values (ρ α λ , S λ , T λ , θ s , θ v , φ, AOD). For the AOD retrieval of these two bands, the lookup table of the corresponding band was selected, and the corresponding parameters (ρ α λ , S λ , T λ ) were derived by obtaining the view geometric information (θ s , θ v , φ) of each pixel. Then, it was substituted into equation (1) to obtain the simulated TOA. The real AOD can be interpolated by comparing the real and simulated TOA. In this way, the coast and blue bands of AOD can be obtained accurately in this study, and then the green and red bands of AOD can be obtained by the Angstrom equation [41]: where τ(λ) is the AOD at different wavelengths, α is the wavelength index, and β represents the Angstrom turbidity coefficient. If the AOD of two wavelengths was obtained, α and β can be expressed as: λ 2 −α Substituting α and β into Equation (2), the AOD of other wavelengths can be derived. In addition, to demonstrate the feasibility of AOD retrieval for various bands of Landsat8 OLI, the changes in the sensitivity of the TOA to the surface reflectance at different AOD values were analyzed by using the 6S model. In the simulation, the view geometry was set as θ s = 30 • , θ v = 30 • , and φ = 120 • , and the aerosol model and atmospheric model were set as custom and mid-latitude winter model respectively. Figure 4 shows that the TOAs of the four bands were all sensitive to the surface reflectance, and there was a threshold point at a surface reflectance of 0.2. When the surface reflectance was less than 0.2, the TOA increased with increasing AOD. When the surface reflectance was greater than 0.2, the TOA decreased with increasing AOD. When the surface reflectance was fixed, it can be seen from Figure 4 that the TOA was gradually insensitive to changes in AOD as the wavelength increased, especially in the red band, the TOA changed very little with the growth of AOD. Therefore, it was reasonable to choose the coast and blue bands to retrieve AOD in this paper.  In addition, to demonstrate the feasibility of AOD retrieval for various bands of Landsat8 OLI, the changes in the sensitivity of the TOA to the surface reflectance at different AOD values were analyzed by using the 6S model. In the simulation, the view geometry was set as θ = 30°, = 30°, and = 120°, and the aerosol model and atmospheric model were set as custom and mid-latitude winter model respectively. Figure 4 shows that the TOAs of the four bands were all sensitive to the surface reflectance, and there was a threshold point at a surface reflectance of 0.2. When the surface reflectance was less than 0.2, the TOA increased with increasing AOD. When the surface reflectance was greater than 0.2, the TOA decreased with increasing AOD. When the surface reflectance was fixed, it can be seen from Figure 4 that the TOA was gradually insensitive to changes in AOD as the wavelength increased, especially in the red band, the TOA changed very little with the growth of AOD. Therefore, it was reasonable to choose the coast and blue bands to retrieve AOD in this paper.  The AOD of four bands over the entire atmospheric column (AODcolumn), representing the extinction of all aerosol particles in the vertical direction, is obtained through the above steps. However, the result measured by the ground monitoring site is the PM2.5 value near

PM 2.5 -AOD Model Building
The AOD of four bands over the entire atmospheric column (AOD column ), representing the extinction of all aerosol particles in the vertical direction, is obtained through the above steps. However, the result measured by the ground monitoring site is the PM 2.5 value near the surface under dry conditions. In order to establish the relationship with the PM 2.5 mass concentration monitored on the ground, the vertical correction and humidity correction should be performed on AOD column to obtain the near-surface dry AOD. The vertical correction, is achieved by assuming that the aerosol vertical distribution follows a negative exponential form, then the near-surface AOD can be derived using the following equation [42,43]: where σ λ represents the near-surface AOD at wavelength λ (in units of µm −1 ), z is the vertical height (in units of km), and H is the aerosol scale height (in units of km), which can be approximated by the PBLH [44]. Regarding the humidity correction, it is necessary to obtain the extinction properties of dry particles and the aerosol particle hydroscopic growth function f r (RH) (with no unit), which have been studied for a long time [45][46][47]. The aerosol particle hydroscopic growth function can be expressed by the following equation: where RH is the relative humidity of the atmosphere at ground level (with no unit) and r and r dry represent the dry and wet radii, respectively. The values of a and b can be found in the related literature [48]. Then, the optical hydroscopic growth function f o (RH) can be expressed as: Then the near-surface dry AOD can be expressed by the following equation [49,50]: where σ dry,λ (in units of µm −1 ) is the near-surface dry AOD at wavelength λ. and then the σ dry,λ of the four bands can be obtained. Assuming that all particles in the atmosphere are spherical particles, according to Mie scattering theory, the near-surface dry AOD (σ dry,λ ) can be expressed as [51][52][53]: where Q(λ, r, m) is the extinction coefficient of a single particle (with no unit), r is the particle radius (in units of µm), r M is the upper limit radius of extinction particles, λ is the wavelength (in units of µm), m represents the particle complex refractive index, which is composed of real and imaginary parts, and n(r) is the aerosol particle size distribution (in units of µm −4 ). Defining K n (λ, r, m) = πr 2 Q(λ, r, m), Equation (8) can be rewritten as: K n (λ, r, m) was defined as the aerosol particle size distribution weight function (in units of µm 2 ). For the volume distribution of aerosol particles, based on v(r) = 4 3 πr 3 ·n(r), K v (λ, r, m) = πr 2 Q(λ,r,m) ume distribution weight function (in units of µm −1 ). Then, Equation (9) can be changed to: where v(r) is the aerosol particle volume distribution (in units of µm −1 ). The complex normal logarithmic volume distribution, which can describe various aerosol modes, was chosen to simulate the particle volume distribution in this study [54]. Aerosol particles are divided into four types based on the SRA model: Water-soluble, dust, marine, and soot. Each aerosol type is composed of aerosol components in different proportions. The specific expression is as follows: where V i (i = 1, 2, 3, 4) represents the water-soluble, dust, marine, and soot particle volumes, σ i is the standard deviation of different particle radii, and r vi is the mean volume median radius (in units of µm). The values of parameters are shown in Table 5. To solve Equation (10), the weight function K v (λ, r, m) must be obtained. The key is to know the extinction coefficient Q(λ, r, m), which can be calculated by the following equation based on Mie theory [56]: where α = 2πr λ , representing the scale parameter, and a n and b n are Mie scattering parameters that are related to the complex refractive index of the particle. The extinction coefficient of different aerosol particles is derived by obtaining the complex refractive index of these four particles [55]. Figure 5 shows that the extinction coefficient of different aerosol components varies with particle radius at the center wavelength of different bands. The trends of different bands were the same, but the peaks gradually moved to the shortwave direction for different bands. With increasing particle radius, the extinction coefficient oscillated and converged to 2. Then, the volume weight function K v (λ, r, m) can be determined based on the previous definition.
If the σ dry,λ at different bands are known, the volume distribution v(r) can be solved by specific methods. Then, the mass concentration of PM 2.5 can be expressed as: where ρ represents the dry mass density of fine particles near the ground (in units of g/cm 3 ). The value of ρ in the PRD can be found in the relevant literature [57]. aerosol components varies with particle radius at the center wavelength of different bands. The trends of different bands were the same, but the peaks gradually moved to the shortwave direction for different bands. With increasing particle radius, the extinction coefficient oscillated and converged to 2. Then, the volume weight function ( , , ) can be determined based on the previous definition. If the , at different bands are known, the volume distribution ( ) can be solved by specific methods. Then, the mass concentration of PM2.5 can be expressed as: where represents the dry mass density of fine particles near the ground (in units of g/cm 3 ). The value of in the PRD can be found in the relevant literature [57].

Volume Distribution Solution
Combining Equations (10) and (11) yields the following equation: where K vi (λ, r, m) is the volume weight function of different aerosol particles (i = 1, 2, 3, 4 represent water-soluble, dust, marine, and soot). Then, the definite integral is discretized, and the complex trapezoidal formula is taken: where h = y−x n ,r j = x + (j − 1)h, j = 1, 2 · · · n. The values of a and b depend on the extinction radius range of different aerosol particles, and Figure 5 shows the range is from 0.01 µm to 10 µm because when the radius is less than 0.01, the extinction coefficient tends to 0, and when the radius is greater than 10, the extinction coefficient oscillates and converges to 2. In addition, n is the degree of discretization with a value of 1000 in this study, so x = 0.01, y = 10, n = 1000. Assuming that the total volume of aerosol particles per column is V and that the volumetric proportion of four aerosol particles is ω i and 2, 3, 4 represent water-soluble, dust, marine, and soot), the volume distribution of different particles can be expressed as: Equation (17) is defined as: Substituting Equation (17) into Equation (15): Equation (19) is defined as: δ i (λ) can be calculated by these known parameters, denoted as |δ i (λ)|, which is a 4*1 matrix. Then, Equation (19) can be converted as: There are four variables in Equation (18) because ω i satisfied ∑ 4 i=1 ω i = 1. Equation (20) will be solved by the σ dry,λ of the four bands obtained earlier and then substituted into Equation (16) to obtain v(r). Finally, the PM 2.5 mass concentrations can be obtained through Equation (13).

Model Evaluation
We used two kinds of data to verify the AOD and PM 2.5 evaluation results. One was MODO4 DB AOD (550 nm) products, which were used to verify the Landsat8 OLI AOD (550 nm) retrieval results synchronously. Only if the retrieval results of AOD (550 nm) are reliable can the evaluation of PM 2.5 be accurate. The second was ground monitoring site PM 2.5 data; there were 52 ground monitoring stations in the study area, and the hourly average values were used to verify the estimated PM 2.5 results. In addition, the root mean square error (RMSE) was used as an evaluation index for the AOD and PM 2.5 results. The overall uncertainty of the model was evaluated using error propagation theory [58,59], and the errors of PM 2.5 can be written as:

Retrieved Results from the Proposed Model
Using the physical model proposed above, the four-year Landsat8 OLI images of the PRD area in spring and winter were selected as the retrieval target because more rainfall leads to poor data in summer and autumn. The cloud coverage of all images was less than 10%, which was convenient for better retrieval results. In addition, water and land separation and cloud mask extraction were carried out before inversion, and no inversion was performed in the areas of water and clouds. Figure 6 shows the temporal and spatial distributions of Landsat8 OLI PM 2.5 in the PRD for a certain day in different years. The In addition, the pollution sources of several major industrial cities in the PRD can be identified through the ultrahigh resolution retrieval results of PM2.5 concentrations. Therefore, the areas with high PM2.5 values in three years (except for 14 November 2019) were used to analyze the pollution status of Guangzhou, Foshan, and Dongguan cities. Figure 7 shows that the PM2.5 mass concentrations in downtown Guangzhou and Panyu District were relatively high, and the PM2.5 mass concentrations for different days in In addition, the pollution sources of several major industrial cities in the PRD can be identified through the ultrahigh resolution retrieval results of PM 2.5 concentrations. Therefore, the areas with high PM 2.5 values in three years (except for 14 November 2019) were used to analyze the pollution status of Guangzhou, Foshan, and Dongguan cities. Figure 7 shows that the PM 2.5 mass concentrations in downtown Guangzhou and Panyu District were relatively high, and the PM 2.5 mass concentrations for different days in the three years reached 80 µg/m 3 or more, while the Zengcheng and Conghua Districts in northern Guangzhou were relatively lightly polluted. The areas with high PM 2.5 mass concentrations in Foshan city were mainly distributed in Nanhai District and Sanshui District. The PM 2.5 mass concentrations in Nanhai District exceeded 100 µg/m 3 on 23 October 2017. The overall pollution area in Dongguan was relatively large and unevenly distributed. The main pollution areas were concentrated in western, central, and southeastern Dongguan. Relative to other parts of the study area, areas adjacent to Guangzhou and Shenzhen were seriously polluted. Through the above ultrahigh resolution PM 2.5 retrieval results in different cities, the range of urban scale pollution can be clearly known.  Figure 8 presents the spatial and temporal distributions of Landsat8 OLI AOD (left) and MODO4 DB AOD (right) in the PRD for a certain day in different years. We compared MODO4 DB AOD and Landsat8 OLI AOD and found that Landsat8 OLI AOD tended to  Figure 8 presents the spatial and temporal distributions of Landsat8 OLI AOD (left) and MODO4 DB AOD (right) in the PRD for a certain day in different years. We compared MODO4 DB AOD and Landsat8 OLI AOD and found that Landsat8 OLI AOD tended to be slightly higher than MODO4 DB AOD, which may be caused by different resolutions or a slight difference in transit time. In addition, MODO4 DB AOD products sometimes lack too much area and cannot provide abundant details, due to their coarse resolution (3 km); nevertheless, Landsat8 OLI AOD at a spatial resolution of 30 m can satisfy the needs of urban scale monitoring. To further verify the correlation between the MODO4 DB AOD and Landsat8 OLI AOD, we resampled Landsat8 OLI AOD to the same resolution as MODO4 DB AOD and then conditionally picked verification points in the area that had MODO4 DB AOD data for comparison. In the process of selecting verification points, to avoid areas with high MODO4 DB AOD values caused by clouds and MODO4 DB AOD areas above waters, we eliminated verification points with high AOD values (greater than 1.5), so that they better match with Landsat8 AOD for verification. Figure 9 shows a reasonable relationship with a correlation coefficient (R 2 ) of 0.66 and RMSE of 0.1037 between Landsat8 OLI AOD and MODO4 DB AOD at 550 nm. The slope (0.7893) was less than 1, indicating that Landsat8 OLI AOD was generally higher than MODO4 AOD. The MODO4 DB AOD was slightly larger than the Landsat8 OLI AOD only when the Landsat8 OLI AOD was small. areas above waters, we eliminated verification points with high AOD values (greater than 1.5), so that they better match with Landsat8 AOD for verification. Figure 9 shows a reasonable relationship with a correlation coefficient (R 2 ) of 0.66 and RMSE of 0.1037 between Landsat8 OLI AOD and MODO4 DB AOD at 550 nm. The slope (0.7893) was less than 1, indicating that Landsat8 OLI AOD was generally higher than MODO4 AOD. The MODO4 DB AOD was slightly larger than the Landsat8 OLI AOD only when the Landsat8 OLI AOD was small.

Comparison Landsat8 OLI PM2.5 with Ground Measurements
PM2.5 data from ground monitoring sites were used to verify the accuracy of the retrieval results. The time difference between the ground monitoring sites and satellite transit time was less than 30 minutes, and the average PM2.5 value of 3*3 pixels around the ground monitoring site was compared with the PM2.5 value of the ground monitoring site. Figure 10 shows the PM2.5 absolute error diagram between the estimated results and measured results. The total RMSE was 10.0224 μg/m³. The RMSEs of a certain day in different years, including 2017, 2018, 2019, and 2020, were 11.9470 μg/m³, 11.9787 μg/m³, 7.4217 μg/m³, and 5.4723 μg/m³, respectively. When the PM2.5 value was approximately 50 μg/m³, the relative error was the smallest. When the PM2.5 value was greater than 60 μg/m³, the error increased. Among the evaluation results of the four years, the results from 14 November 2019 and 18 February 2020 were more accurate; in contrast, on 23 October 2017 and 12 February 2018, severe air pollution led to poor accuracy. The overall error of AOD and PM2.5 was small, which reflected the feasibility of the algorithm. 3.3. Comparison Landsat8 OLI PM 2.5 with Ground Measurements PM 2.5 data from ground monitoring sites were used to verify the accuracy of the retrieval results. The time difference between the ground monitoring sites and satellite transit time was less than 30 minutes, and the average PM 2.5 value of 3 × 3 pixels around the ground monitoring site was compared with the PM 2.5 value of the ground monitoring site. Figure 10 shows the PM 2.5 absolute error diagram between the estimated results and measured results. The total RMSE was 10.0224 µg/m 3 . The RMSEs of a certain day in different years, including 2017, 2018, 2019, and 2020, were 11.9470 µg/m 3 , 11.9787 µg/m 3 , 7.4217 µg/m 3 , and 5.4723 µg/m 3 , respectively. When the PM 2.5 value was approximately 50 µg/m 3 , the relative error was the smallest. When the PM 2.5 value was greater than 60 µg/m 3 , the error increased. Among the evaluation results of the four years, the results from 14 November 2019 and 18 February 2020 were more accurate; in contrast, on 23 October 2017 and 12 February 2018, severe air pollution led to poor accuracy. The overall error of AOD and PM 2.5 was small, which reflected the feasibility of the algorithm.

Error Estimation
This study was based on multiband AOD to establish a physical model of PM2.5 mass concentration retrieval. Therefore, the main error comes from the multiband AOD and the parameters in the physical model. The model parameterization error mainly comes from ( ). The average uncertainty in the hygroscopic growth item was estimated to be 24%, according to Table 6. Water and land separation was performed based on the Landsat8 OLI datasets, so the uncertainty of urban or mixed aerosols was reduced to 12%, as shown in Table 6. The uncertainty of AOD was determined by the comparison between Landsat8 OLI AOD and MODO4 DB AOD at 550 nm. Section 3.2 shows that the relative error between the two was 0.1251. In addition, the fine particle density can also vary in different regions. This study adopted 1.5 g/cm³ as the fine particle density in the PRD region, which produced an error of approximately 28% compared with the national average fine particle mass density of 2.09 g/cm³ [59]. In summation, the total uncertainty of PM2.5 was approximately 33% following equation (21). The uncertainty of parameters, such as PBLH and RH, depends greatly on the sensors and retrieval algorithm and needs improvement in future work.

Error Estimation
This study was based on multiband AOD to establish a physical model of PM 2.5 mass concentration retrieval. Therefore, the main error comes from the multiband AOD and the parameters in the physical model. The model parameterization error mainly comes from f o (RH). The average uncertainty in the hygroscopic growth item was estimated to be 24%, according to Table 6. Water and land separation was performed based on the Landsat8 OLI datasets, so the uncertainty of urban or mixed aerosols was reduced to 12%, as shown in Table 6. The uncertainty of AOD was determined by the comparison between Landsat8 OLI AOD and MODO4 DB AOD at 550 nm. Section 3.2 shows that the relative error between the two was 0.1251. In addition, the fine particle density ρ can also vary in different regions. This study adopted 1.5 g/cm 3 as the fine particle density in the PRD region, which produced an error of approximately 28% compared with the national average fine particle mass density of 2.09 g/cm 3 [59]. In summation, the total uncertainty of PM 2.5 was approximately 33% following equation (21). The uncertainty of parameters, such as PBLH and RH, depends greatly on the sensors and retrieval algorithm and needs improvement in future work. Table 6. Hygroscopic growth functions on particle radius (r) and particle extinction (o). The error column shows the difference between two hygroscopic growth functions f r 3 (RH) and f o (RH) [29].

Discussion
This paper presented a new physical model based on multiband AOD for PM 2.5 estimation. The results showed ultrahigh spatial resolution and abundant spatial details, which can be applied to urban scale atmospheric monitoring and can be combined with other geographic data to analyze pollution characteristics in different cities. Figure 11a-i shows the pollution characteristics of the three major industrial cities (Guangzhou, Foshan, and Dongguan) in the PRD of a certain day in different years. The pollution in Guangzhou (Figure 11a-c) and Foshan (Figure 11d-f) was mainly industrial pollution, and the denser industrial zones had higher PM 2.5 values. The pollution characteristics of Dongguan were linearly distributed along the roads, indicating that traffic pollution was the main cause, which was very similar to the pollution situation in northern Poland [60]. With the previous coarse resolution PM 2.5 retrieval results, which pay more attention to PM 2.5 estimation at a large scale, it is difficult to obtain pollution characteristics at the urban scale. The ultrahigh resolution retrieval results from this approach can be provided to the urban environmental department to support better regulatory measures.

Conclusions
In this study, we used Landsat8 OLI images to retrieve the AOD of the coast and blue bands based on the 6S model, and then the Angstrom formula was used to obtain the AOD of the green and red bands. Then, the dry near-surface AOD of visible bands was obtained through vertical correction and humidity correction, and the aerosol particles in Figure 11. Analysis of pollution characteristics of the three major industrial cities (Guangzhou: (a-c); FoShan: (d-f); DongGuan: (g-i)) in PRD for a certain day in different years.
However, there were still some limitations and deficiencies in this study. The first was that long-term ground monitoring data were needed to determine the aerosol mode in the study area; otherwise, it was difficult to obtain accurate AOD, which had a great impact on the subsequent PM 2.5 retrieval. The second was that the simple aerosol vertical distribution function had difficulty representing the vertical distribution of aerosols in the entire PRD area, and the uniform use of the negative exponential distribution function inevitably caused errors in the retrieval results. Third, this study used Mie scattering theory to simulate the extinction characteristics of different aerosol particles, which varied from the extinction characteristics of real aerosol particles in the atmosphere. Measurement experiments of extinction characteristics should be performed on different types of aerosols in the study area to obtain their true extinction characteristics, which will make the physical model of PM 2.5 -AOD more suitable for the real conditions in the study area.

Conclusions
In this study, we used Landsat8 OLI images to retrieve the AOD of the coast and blue bands based on the 6S model, and then the Angstrom formula was used to obtain the AOD of the green and red bands. Then, the dry near-surface AOD of visible bands was obtained through vertical correction and humidity correction, and the aerosol particles in the PRD were divided into four different types according to the SRA model. The extinction characteristics of different types of aerosol particles were simulated by the Mie scattering theory. Finally, a physical model was established between PM 2.5 and dry near-surface AOD to obtain the spatial and temporal distribution characteristics of PM 2.5 in the PRD. The resolution of the retrieved results was the same as that of the original images, and the distribution of PM 2.5 at the urban scale was analyzed by using ultrahigh resolution remote sensing results. The retrieval results showed the distribution of PM 2.5 mass concentrations in the PRD for a certain day in different years, and the cities with higher degrees of industrialization (such as Guangzhou, Foshan, and Dongguan) were analyzed in detail. The urban areas of Guangzhou and Panyu District had high pollution except on 14 November 2019, and the PM 2.5 mass concentrations exceeded 80 µg/m 3 . The pollution in Foshan was mainly in Nanhai District, while the PM 2.5 mass concentrations in Dongguan were mainly distributed in the western, central, and southeastern regions. In addition, ultrahigh resolution retrieval results were also used to analyze the pollution characteristics of different cities on the day; the industrial pollution was dominant in Guangzhou and Foshan, while traffic pollution was the main factor in Dongguan which was consistent with other studies on pollution characteristics of industrial cities [61,62]. Finally, the retrieval results were verified by using MODO4 DB AOD products and ground monitoring station data. The AOD at 550 nm retrieved in this study had high consistency with MODO4 DB AOD products. The RMSE of the PM 2.5 mass concentration retrieval results and ground monitoring station results was only 10.0224 µg/m 3 , which indicated that the retrieval results were reliable.
It should be mentioned that the resolution of parameters, such as PBLH and RH, used in this study was 0.125 • × 0.125 • , which was different from the resolution of Landsat8 OLI data. Corresponding errors are introduced when making the humidity correction and vertical correction. In addition, the 30 m AOD retrieval result was resampled to 3 km for verification with the 3 km MODO4 DB AOD products, which introduces certain verification errors. Despite some limitations, our proposed physical model can still provide some information support to urban environmental departments. In the future, we will carry out a quantitative comparison based on the vertical profile of the dissipation coefficient by Landsat8 and CALIPSO [63] and the generation of PBLH images covering the original images [64,65] or use Lidar to perform multipoint measurements of the study area to obtain an accurate vertical distribution of aerosols, which will greatly improve the vertical correction [66]. In addition, we will also conduct different types of aerosol extinction coefficient measurements in the study area to improve the model.