PM2.5 Exposure and Health Risk Assessment Using Remote Sensing Data and GIS

Assessing personal exposure risk from PM2.5 air pollution poses challenges due to the limited availability of high spatial resolution data for PM2.5 and population density. This study introduced a seasonal spatial-temporal method of modeling PM2.5 distribution characteristics at a 1-km grid level based on remote sensing data and Geographic Information Systems (GIS). The high-accuracy population density data and the relative exposure risk model were used to assess the relationship between exposure to PM2.5 air pollution and public health. The results indicated that the spatial-temporal PM2.5 concentration could be simulated by MODIS images and GIS method and could provide high spatial resolution data sources for exposure risk assessment. PM2.5 air pollution risks were most serious in spring and winter, and high risks of environmental health hazards were mostly concentrated in densely populated areas in Shanghai-Hangzhou Bay, China. Policies to control the total population and pollution discharge need follow the principle of adaptation to local conditions in high-risk areas. Air quality maintenance and ecological maintenance should be carried out in low-risk areas to reduce exposure risk and improve environmental health.


Introduction
Air pollution has been considered a global health priority regarding several Sustainable Development Goals (SDGs), such as Goal 3 (Ensure healthy lives and promote well-being for all at all ages) [1,2]. Fine particulate matter (a diameter of less than 2.5 µm, PM 2.5 ) pollution is a common type of air pollution in recent years. According to the China Ecological Environmental Bulletin [3][4][5][6][7], 40~75% cities in China exceeded the standard for ambient air quality (PM 2.5 < 75 µg/m 3 ), and on 45~80% of days, PM 2.5 as the primary pollutant exceeded the standard. One study by Son et al. found that short-term exposure to PM 2.5 was positively associated with increased risk of mortality [8]. Other studies indicated that the long-term chronic effects of PM 2.5 may cause cardiovascular diseases such as lung cancer, myocardial infarction, and myocardial ischemia [9][10][11], and it was an important cause of acute triggering of common respiratory system diseases such as asthma, bronchitis, rhinitis, and upper and lower respiratory tract infections [12][13][14]. Long-term exposure to PM 2.5 pollution may lead to slow growth, slow neurological development, and brain dysfunction in children [15,16]. In addition, it may lead to depression and pessimism, and even suicidal behavior [17][18][19]. Due to its critical impact on health, exposure to PM 2.5 and health risk assessment has been a critical concern for ensuring healthy lives and promoting well-being.
There are various data and methods for quantifying ambient air pollution in recent years. Remote sensing (RS) and Geographic Information Systems (GIS) methods have been increasingly used in environment and health research, with corresponding improvements in the use and accessibility of multi-temporal satellite-derived environmental data [20,21]. Some studies estimated PM 2.5 concentration and spatial-temporal distribution from air quality monitoring stations using spatial interpolation techniques from GIS [22][23][24], but this approach ignored the uneven distribution of monitoring stations and the rough spatial resolution of the data. The rapid development of RS technology and advanced satellite sensors (with particulate matter detection instruments) addressed this issue. Many satellites (e.g., GOES, Terra, Aqua, METOP, PARASOL) equipped with multifunctional sensors (e.g., MODIS, AVHRR, SeaWiFS, POLDER) greatly promoted the development of Aerosol Optical Depth (AOD, one of the most important parameters of aerosols) inversion by remote sensing imagery interpretation and processing [25][26][27][28][29]. Kaufman introduced the dark target algorithm of AOD inversion based on band relationships [30], and Levy et al. improved the accuracy of AOD inversion of the algorithm [31]. Tanré and Holben developed the structural function method by using atmospheric transmittance to obtain aerosol information [32,33]. Hsu et al. directly adopted surface reflectance data on visibles band to retrieve AOD using the deep blue algorithm [34,35]. Lyapustin et al. proposed a new correction algorithm for AOD inversion, the multi-angle atmospheric correction algorithm, using time series analysis and image processing technology for atmospheric correction and AOD inversion [36,37].
Many researchers found a specific relationship between AOD and PM 2.5 that could be estimated by some methods from AOD to PM 2.5 [38][39][40]. The model correction method used different models to correct various influencing factors and simulate the proportion relationship between AOD and PM 2.5 , such as the atmospheric chemical transport model (GEOS-Chem) [41][42][43] and the atmospheric boundary layer model (RAMS) [44]. However, this method ignored the physical mechanism between AOD and PM 2.5 . Both aerosol type and vertical distribution lead to differences in scattering extinction. The mechanism correction method solved this problem and obtained the stable aerosol extinction coefficient by vertical correction and scattering extinction correction for PM 2.5 estimation [45][46][47][48]. However, this method was highly dependent on setting physical mechanism parameters, and these parameters are different in different areas. In order to monitor and estimate PM 2.5 in real time, it is necessary to update these physical mechanism parameters in time. In addition, the statistical model method established linear or nonlinear statistical models between the AOD and PM 2.5 based on various meteorological or environmental elements (wind speed, direction, position, humidity, height, etc.), such as multiple linear regression model (GLM) [49,50], generalized summation model (GAM) [51,52], and geographically weighted regression model (GWR) [53][54][55][56][57]. Some researchers combined several statistical models to construct multilevel statistical models to estimate PM 2.5 concentration. In recent years, the machine learning (ML) method has been widely used to associate AOD with PM 2.5 , incorporating big geographic spatialtemporal data as examples and with self-supervision and training functions [58,59]. Although this method had high accuracy in estimating results and could deal with the complex relationship between AOD and PM 2.5 , it required multistep processing of training samples in advance and relevant operations on physical and chemical mechanisms, which increased the difficulty of use to a certain extent [60]. Many methods have associated AOD with PM 2.5 , each with its characteristics.
To explore the relationship between PM 2.5 air pollution and health exposure risk, one study by Zou et al. found a spatial pattern of population exposure to air pollution by constructing a relative exposure risk assessment model [61]. Tong et al. showed that population density data and the relative exposure risk assessment model could more reasonably represent the relationship between PM 2.5 pollution and environmental health [62]. Lu et al. developed a personal mobility model to quantify long-term air pollution exposure for individuals [63]. Park developed an air dispersion model based on a large sample of travel-activity diary data to assess personal exposure to PM 2.5 [64]. These studies confirmed the applicability to and efficiency of remote sensing data and GIS for exposure to PM 2.5 air pollution and environmental health. However, risk assessment results from these studies were not sufficiently accurate owing to the limited availability of high spatial resolution data for PM 2.5 and population density.
Toward this end, this study presents a framework that incorporates high spatial resolution remote-sensing images and geographic spatial-temporal data to retrieve AOD and estimate PM 2.5 based on the GIS platform. The Enhanced Dark Target Algorithm (EDTA) is introduced to retrieve daily AOD from MODIS images, improving the spatial resolution at 1 km grid level. Spatial-temporal seasonal models are used to estimate PM 2.5 concentration using geospatial-temporal data and GIS spatial analysis methods. Then, we assess the exposure risk to PM 2.5 pollution using high-accuracy population density data and the relative exposure risk models. Finally, we discuss the usage and deficiency of assessment results and offer our expectations for future application. We look forward to provide scientific reference for improving urban atmospheric pollution and living environment.

Study Area
SHB is located in eastern China and the north Pacific coastal area (28 • 51 ~31 • 53 N, 118 • 21 ~123 • 25 E) and is an important part of the integrated development of the Yangtze River Delta, China. From north to south, it contains seven cities, including Shanghai, Jiaxing, Huzhou, Hangzhou, Shaoxing, Ningbo and Zhoushan. Its total land area is about 52,300 km 2 , and its population is about 54.41 million (2018 Statistical Yearbook of China). Belonging to the northern subtropical region, SHB is characterized by a mild and humid climate with abundant rainfall. The terrain is higher in the southwest and lower in the northeast; plains, hills and mountains cover the most part (

Vector and Elevation Data
The data on the region boundaries and elevation were from the Resource and Environmental Sciences and Data Center, China (ESDC, http://www.resdc.cn/Default.aspx, accessed on 18 May 2019).

Remote Sensing Data
Terra MODIS L1B data (Level 1, with Global Certification Level) were applicable and reliable for China [65,66] at a spatial resolution of 1 km ×1 km. Excluding a period of no data from 19 to 27 February, there were 1 to 3 image data every day, for a total of 827 images of original data from SHB from 1 January to 31 December 2016. Data were downloaded from NASA LAADS Web to retrieve AOD (LAADS, https://ladsweb.nascom.nasa.gov/, accessed on 4 December 2019).

Ground-Based AOD Observation Data
The observation error for ground-based AOD from AERONET is 0.01~0.02, which can be used directly for the correction of estimated AOD [67]. The Version 3 and Level 1.5 data (data after cloud processing) in 2016 were downloaded from the Aerosol Robotic Network (AERONET, http://aeronet.gsfc.nasa.gov/, accessed on 18 May 2019) for verifying the accuracy of retrieved AOD. We finally chose Level-1.5 AOD data from two stations (SONET_Shanghai and SONET_Zhoushan) according to the actual situation in SHB (Table 1).

Ground-Level PM 2.5 Observation Data
Based on China national standards and specifications for ambient air quality assessment, daily observation data on PM 2.5 concentrations (unit: µg/m 3 ) at air quality monitoring stations were derived from the official website of the China Environmental Monitoring Center for exploring the optimal relationship between retrieved AOD and PM 2.5 observation data (CEMC, http://106.37.208.233:20035/, accessed on 23 December 2019). There were 41 monitoring stations in SHB, and specific information is in Table 2.

Population Density Data
The population density data came from the open population data set platform Worldpop (https://www.worldpop.org/, accessed on 23 December 2019), with a spatial resolution of 1 km × 1 km. Based on remote sensing images and geospatial-temporal big data, this data set used the random forest algorithm to simulate the spatial-temporal distribution of the population density, comprehensively considering land cover and land use, residential areas, roads, buildings, public facilities, night lights, vegetation, geographic and geomorphic conditions, and so on [68]. The dark target algorithm (DTA) is an operational algorithm with high inversion accuracy and maturity. Kaufman found a good linear relationship between surface reflectivity in the near-infrared band (2.1 µm) and the visible red band (0.66 µm) and visible blue band (0.47 µm) [69]. The near-infrared band is affected very little by aerosols, and its apparent reflectivity can be approximated as surface reflectivity [26]. The surface reflectivity of the visible red and blue bands can be estimated with a linear calculation equation: where ρ r is the surface reflectivity of the visible red band (0.66 µm); ρ b is the surface reflectivity of the visible blue band (0.47 µm); and ρ n is the surface reflectivity of nearinfrared band (2.1 µm). Apparent reflectivity is the result of surface reflectivity and atmospheric reflectivity, as shown in Equation (3). The apparent reflectivity information of the visible red and blue bands can be obtained from satellite remote sensing data, and the surface reflectivity of the visible red and blue bands can be obtained after removing the linear estimation, which gives the atmospheric reflectivity [70]: where ρ * is apparent reflectivity; µ is the cosine of the satellite zenith angle, µ 0 is s the cosine of the solar zenith angle, ϕ is the satellite azimuth angle, ϕ 0 is the solar azimuth angle; ρ a is atmospheric reflectivity; T(µ 0 ) is the total transmittivity from the sun to the earth's atmosphere; T(µ) is the total transmittivity from the earth's surface to the satellite's atmosphere; ρ is surface reflectivity; and s is the spherical surface albedo of the atmosphere. Based on the basic dark target algorithm (DTA), we introduced the enhanced dark target algorithm (EDTA) to retrieve the AOD according to the actual situation in SHB. The basic inversion process included radiometric correction, geometric correction, resampling, composition and clipping, cloud detection and elimination, Lookup Table (LUT) setting, and accuracy verification of the inversion results ( Figure 2). We improved the inversion algorithm, especially in building LUT.
Geometry correction. MODIS L1B 1 km data (MOD02_1KM) contains emissivity and reflectivity files and angle data (sensor zenith and sensor azimuth of the satellite, solar zenith and solar azimuth of the sun), which are obviously different types. We employed ENVI software to resample the row and column numbers of angle data set from 271 × 406 to 1354 × 2030 (e.g., emissivity and reflectivity files) before correction. The HDF file was used to generate ground control points (GCPs). The correction model was Triangulation, and the resampling method was Bilinear.
Band operation and cloud detection. The Layer Stacking tool was adopted for synthesizing geometry files and angle data after correction. The stacked order of emissivity and reflectivity files (as well as angle data) affected the results. The reflectivity file must be placed in up and the emissivity file in down. In fact, angle data had been expanded 100 times in HDF files, so it should have been multiplied by 0.01 during band operation. In order to remove the influence of cloud reflection, absorption and scattering noise, the Cloud Detection tool was installed.
LUT Setting and AOD retrieval. The accuracy of LUT determines the accuracy of AOD retrieval to a certain extent. A 6S Radiation Transmission model was used to distinguish various surface types and observation bands. Different parameters of atmospheric aerosol and observation parameters were preset for radiation transmission calculation to obtain inversion results. After geometric correction, band operation, band clipping, cloud detection, and invalid value elimination, the results for emissivity, reflectance, and angle data set were combined with the LUT to perform aerosol inversion calculation, and we obtained AOD inversion values. According to geographical features and the actual situation in SHB, we improved the related model and enhanced the LUT setting in the geometric parameters (Table 3) to produce the enhanced dark target algorithm (EDTA).  Validation. The ground-based AOD from AERONET data was used to verify the accuracy of the retrieved AOD. We set the value of AOD at 550 nm wavelength (0, 0.25, 0.50, 1.00, 1.50, 1.95) in LUT, but the AERONET data were available in 340 nm, 380 nm, 440 nm, 500 nm, 675 nm, 870 nm, 936 nm, 1020 nm, and 1640 nm wavelength channels (no AOD observation value in 550 nm). Only AOD values of the same observation band could be compared, so the AOD values of AERONET between different observation band channels needed to be converted. We used the Angstrom formula to convert [71]: where λ is the wavelength; τ λ is AOD of the wavelength λ; β is the atmospheric turbidity index; and α is the wavelength index.

AOD-PM 2.5 Spatial-Temporal Regression Models
When observing the relationships between AOD samples and PM 2.5 samples, we found that PM 2.5 tended to increase gradually with increasing AOD. Therefore, AOD samples were taken as the independent variable and PM 2.5 samples as the dependent variable. In order to infer the relationship between AOD (independent variable) and PM 2.5 (dependent variable), we mainly used 6 linear and nonlinear regression models, including linear, logarithmic, exponential, power, quadratic polynomial, and cubic polynomial function regression models (Table 4). By combining terrain, landscape, and other geospatial information, we applied and tested each model pixel by pixel based on seasonal characteristics and the actual situation in SHB. Table 4. Regression models for AOD-PM 2.5 relationship prediction.

Regression Model Equation
Linear y = a 0 + a 1 x Logarithmic y = a 0 + a 1 ln(x) Exponential y = a 0 × e a1x Power y = a 0 (x a1 ) Quadratic Polynomial y = a 0 + a 1 x + a 2 x 2 Cubic Polynomial y = a 0 + a 1 x + a 2 x 2 + a 3 x 3 x for independent variable. y for dependent variable. a 0 , a 1 , a 2 , a 3 for relevant parameters.
The first part of the AOD samples (estimated from the MODIS images) and the PM 2.5 samples (observed from the monitoring stations) was used for model modeling, and the second part was used for model testing. The optimal model was determined based on model fit (determinant coefficient, R 2 ) and error results (root mean square error, RMSE) (Equations (5) and (6)) to build a spatial-temporal estimation model of AOD-PM 2.5 suitable for SHB and to simulate the spatial-temporal pattern of PM 2.5 concentration based on the grid: where R 2 is the determine coefficient; RMSE is the root mean square error; n is the sample number;ŷ i is the value of independent variable i; y i is the value of dependent variable i; and y is the mean value of the dependent variable.

Pearson's and Spearman's Rank Correlation Coefficients
Pearson's correlation coefficient (R) was applied to the correlation analysis between observation AOD (ground-based from AERONET) and inversion AOD (estimated from MODIS images), seeing by Equation (7): where R is the Pearson correlation coefficient; n is the sample number; z i is the observation AOD value of i; z is the average value of observation AOD; u i is the inversion AOD value of i; and u is the average value of inversion AOD. Spearman's rank correlation coefficient (ρ) was used to indicate the correlation between estimated AOD and observed PM 2.5 as in Equation (8): where ρ is the Spearman's rank correlation coefficient, G i is the rank of estimated AOD; H i is the rank of observed PM 2.5 value; is the rank difference of estimated AOD and observed PM 2.5 value.

Relative Exposure Risk Model
Based on the estimated PM 2.5 concentration and population density data at 1-km grid level, this study used the relative exposure risk model to assess the resident exposure level to PM 2.5 air pollution [61,72], as in Equation (9): where i is the grid number; Q i is the relative population exposure risk of i; P i is the population density of i (unit: person/km 2 ); M i is the PM 2.5 concentration of i (unit: µg/m 3 ); and n is the total number of grids.

Spatial Autocorrelation Analysis
Moran's I is used to represent the spatial autocorrelation, including the global Moran's I and the local Moran's I. The global Moran's I is for the spatial autocorrelation of variables in the study area as a whole. When the value of I approaches 1, the correlation of variables in the spatial distribution is more significant; when the value of I approaches 0, the correlation is weaker. The local Moran's I refers to the correlation degree between the local area and the surrounding area, and its results can be shown in the LISA agglomeration figure. Its I i value is calculated by Equations (10) and (11): where n is the number of grids, i = j, x i , and x j are values of variables of i and j, x is the average value, and w i,j is the weight matrix for the proximity between i and j. On the significance test, the significance level is α = 0.05. The Results of I i are divided into five types for the spatial agglomeration characteristics: HH (high-high) is the high-value agglomeration phenomenon, LL (lowlow) is the low-value agglomeration phenomenon, LH (low-high) and HL (high-low) are high values alternate with low values, and NS (not significant) is no obvious agglomeration features.

AOD Inversion Results
The value range of AOD inversion results was from 0 to 1.95, which was a dimensionless value. The larger the AOD was, the greater the aerosol optical thickness was, indicating the lower atmospheric transmittance. The 0 value indicated the existence of no aerosol particles, which was the best atmospheric condition and indicated that solar radiation was not reduced when it passed through the atmosphere.
Because the remote sensing images were large and susceptible to cloud influence, the AOD inversion process and results would have been affected by large areas of cloud in those time periods. We screened MODIS images from all of the year of 2016 one by one, and we excluded images that were seriously blurred by clouds. Finally, we obtained monthly AOD results for SHB (from average daily AOD) for February, March, April, May, June, July, August, September, November, and December (10 months in total) and seasonal AOD results (from average monthly AOD) for spring (March, April, May), summer (June, July, August), autumn (September, November), and winter (February, December) in SHB.

Monthly AOD Results
In the monthly AOD results in Figure 3, higher AOD values were concentrated in Shanghai, Jiaxing, northwest of Hangzhou, Shaoxing, and northern Ningbo, and lower AODs were concentrated in southwest and south of SHB, namely, southwest of Hangzhou, south of Shaoxing and Ningbo, and Zhoushan. Because some MODIS images in a certain time were affected by cloud interference, a small portions of them were eliminated and then showed null values. However, this did not affect the inversion results for other regions, so they continued to participate in the calculation. In 2016, the AOD began increasing steadily beginning in February through May and gradually decreased in June and July, reached the lowest in August; then they gradually increased from September to November and reached the peak in December. Therefore, AOD showed a time variation tendency of "high-low-high" in SHB.

Seasonal AOD Results
In the seasonal AOD results in Figure 4, inversion AOD still presented the distribution characteristics of higher at the coast and lower in the south in spatial scale. In terms of time scale, the overall AOD was lowest in summer and highest in winter, with little difference between spring and autumn, and seasonal AOD still presented a high-low-high tendency in SHB. For higher-value coastal regions, the distribution of low and middle values was interphase in summer, the distribution of middle and high values was interphase in winter, and the distribution of middle values was uniform in spring and autumn, indicating that AODs in the coastal regions were greatly different in summer and winter. Moreover, Zhoushan had the lowest AOD throughout all four seasons, indicating the least PM 2.5 pollution and the best air quality in 2016.

Verification Result
Using Level 1.5 AERONET AOD data (from ground-based AOD observation stations) of "SONET_Shanghai" and "SONET_Zhoushan", we verified the accuracy of the inversion AODs. In 2016, there were 29 days of Level 1.5 data available for SHB. AOD inversion values on those days were extracted through software, and then invalid AOD values (disturbed by clouds) were removed to obtain valid values. Finally, values of inversion AOD and observation AOD could be used for verification, a total of 21 groups of data (Table 5).   The verification results in Table 5 show that inversion AOD and observation AOD in terms of time trend were approximately the same. From the Pearson's correlation analysis, M and SD for the inversion AOD were higher than those for the observation AODs, and R was 0.781 at the 0.01 significance level, showing that the observation AODs were strongly related to the inversion AODs. The AOD result was at a spatial resolution of 1 km × 1 km. Therefore, we verified that the inversion AODs obtained by MODIS data and EDTA were relatively accurate, and the inversion results could be used for the estimation of PM 2.5 concentrations at higher accuracy and resolution.   Values for inversion AOD and observation PM 2.5 were processed by the min-max normalization method using a total of 410 groups ( Table 6). The monthly correlation analysis results showed inversion AOD had the best correlation with observation PM 2.5 in May, with ρ of 0.631 (at 0.01 confidence level), followed by in August, with ρ of 0.607. Since there was a big difference between inversion AOD and observation PM 2.5 values every month in 2016, the correlations between them were calculated by season. In the seasonal correlation analysis in Table 6, we found that the best correlation was in summer, with ρ of 0.684 (at 0.01 confidence level), followed by spring, with ρ of 0.538. Overall, PM 2.5 concentration increased gradually with the increase in AOD. The correlation in each season was better than that in each month. Therefore, seasonal modeling was more effective for PM 2.5 estimation in SHB.

Seasonal Model Building and Verification
Based on the correlation analysis of inversion AOD and observation PM 2.5 , we obtained 123 groups of data in spring, 123 groups in summer, 82 groups in autumn, and 82 groups in winter. A part of each season's data was used for model building, and another part was used for model verification. AOD was the independent variable and PM 2.5 the dependent variable. R 2 represented the model fit degree, and F is for the significance of the model. The larger these values, the more significant and suitable the model is.
In spring, the first 82 groups were used for model building, and the last 35 groups were used for model verification (6 groups of abnormal values eliminated). The seasonal model building results ( Figure 7A) show that the power model was the best in spring, with R 2 = 0.511 and F = 77.209. In summer, the first 84 groups were used for model building, and the last 35 groups were used for model verification (4 groups of abnormal values eliminated), and the exponential model was the best in summer, with R 2 = 0.551 and F = 127.519 ( Figure 7B). In autumn, the first 41 groups were used for model building, and the last 41 groups were used for model verification, and the power model was the best, with R 2 = 0.524 and F = 34.180 ( Figure 7C). In winter, the first 41 groups were used for model building, and the last 41 groups were used for model verification, and the power model was the best, with R 2 = 0.504 and F = 39.556 ( Figure 7D).
The seasonal model verification results in Table 7 show that the power model had better fit and the fewest errors in spring, with R 2 = 0.513 and RMSE = 6.204; the exponential model showed the best fit and the smallest error in summer, with R 2 = 0.640 and RMSE = 3.979; the fit of the power model was the best and the error was smaller in autumn, with R 2 = 0.520 and RMSE = 7.893; and the power model had better fit and minimum error in winter, with R 2 = 0.540 and RMSE = 7.392. Overall, the power model was the best in spring, autumn, and winter and the exponential model was the most suitable for summer, showing that both played significant roles in estimating PM 2.5 concentration, especially the power model. Based on these optimal seasonal models, we combined observation PM 2.5 concentrations and the geographic big data according to the actual situation and features in SHB to produce seasonal spatial-temporal models for PM 2.5 estimation.

PM 2.5 Estimation Results
Based on these seasonal spatial-temporal models, we obtained seasonal PM 2.5 estimation results for SHB ( Figure 8). Spatially, the PM 2.5 concentrations in each season showed the same distribution characteristic of higher in the coastal areas and lower in the mountainous areas. Zhoushan was the city with the lowest PM 2.5 concentration in all seasons, which was always lower than 40 µg/m 3 . Shanghai and Jiaxing were both above 40 µg/m 3 in all seasons, making them the two most polluted cities. From the perspective of time, the average PM 2.5 concentrations over the four seasons was winter > spring > autumn > summer, and the PM 2.5 concentration showed a high-low-high tendency in a year.
From the annual average PM 2.5 results in Figure 8E, we see that the PM 2.5 concentration range was 0~68.16 µg/m 3 , and its spatial resolution was 1 km × 1 km. Concentrations were high (50~70 µg/m 3 ) in Shanghai, Jiaxing, northwest of Hangzhou, Shaoxing, and the northern part of Ningbo; medium (30~50 µg/m 3 ) in the vicinity of high concentration; and low (0~30 µg/m 3 ) mainly in the southwest of Hangzhou and Zhoushan. Overall, the spatial distribution of PM 2.5 concentrations was higher in the northeast and lower in the southwest, indicating that the area along the Hangzhou Bay was seriously affected by fine particles, and the air quality in the mountainous area with higher altitude was better. This correlated with the elevation characteristics of SHB (Figure 1). The coastal areas are of low altitude, with dense populations, convenient transportation, and developed industry and commerce, which is not conducive to haze diffusion. However, the mountainous areas are of high altitude, with high vegetation coverage and abundant rain, so there is less PM 2.5 pollution and good air quality.  In this work, the spatial distribution result of observation PM 2.5 concentration was obtained by the Kriging interpolation of GIS, and its optimal spatial scale could only reach 12 km × 12 km (because of the small number of monitoring stations) ( Figure 8F). Compared with estimates of PM 2.5 concentrations from remote sensing data (1 km × 1 km), although PM 2.5 concentrations are high in similar areas (Shanghai, the junctional zone of Jiaxing, Huzhou, and Hangzhou), the estimated PM 2.5 concentrations are more continuous and even more accurate in some areas, such as east and west of Huzhou and south and north of Shaoxing. Eastern Huzhou has a larger population and a more developed economy, so the impacts of PM 2.5 population on eastern areas should be greater than that on western areas. The same is true in the northern and southern areas of Shaoxing. Therefore, PM 2.5 concentration estimates from remote sensing images are closer to the real situation in SHB, reflecting the necessity of remote sensing to estimate PM 2.5 for assessing fine pollution levels at a regional scale.

Exposure Risk Assessment
Based on population density data from Worldpop (1-km grid level) and the relative exposure risk model, the exposure risk to PM 2.5 pollution was evaluated using GIS. Figure 9A shows that each level of average annual exposure risk was unevenly distributed around the cities, indicating that PM 2.5 air pollution degree differs greatly around SHB. The range of exposure risk values was from 0 to 133, with a large difference between the lowest value and the highest value. In addition, about 2/3 of the land area was a safe zone, and the rest was in a danger zone (all low-elevation urban areas). By city, the average annual exposure risk was Shanghai (3.57) > Jiaxing (1.24) > Ningbo (0.77) > Shaoxing (0.51) > Hangzhou (0.50) > Zhoushan (0.46) > Huzhou (0.44). Except for Shanghai, other cities were in the safe range, with residents of Huzhou city the least exposed to PM 2.5 pollution. Moreover, the danger level of exposure risk in each city was concentrated in the coastal areas around SHB, and the higher elevations (in the south) were all safe, with the lowest PM 2.5 pollution degree. The highest danger level was mainly concentrated in the main urban areas, with Shanghai and Hangzhou as the most serious. On the global Moran's I for SHB, the values were 0.6436 in spring, 0.6174 in summer, 0.6365 in autumn, and 0.6351 in winter. Values of p for all four seasons were all above 0.01, indicating strong positive spatial correlations and high significance of PM 2.5 exposure risk in SHB. The results of the LISA spatial aggregation ( Figure 9B) indicated that the annual average PM 2.5 exposure risk in SHB had a strong spatial autocorrelation and an obvious spatial aggregation. The high-value agglomeration (HH) was distributed in Shanghai, northeastern Hangzhou, central Jiaxing, and central Ningbo. The lowvalue agglomeration (LL) was mainly distributed in the areas with higher elevations, Shanghai and Zhoushan City in the Chongming District. Overall, the coastal areas were dominated by high-value clustering, and mountainous areas were dominated by low-value clustering.

Discussion
Due to the high cost of construction and maintenance of air quality monitoring stations, the number of stations is small in China. Monitoring stations are unevenly distributed, mostly concentrated in relatively developed areas. The observation PM 2.5 concentrations could only reflect the local PM 2.5 situations but could not truly reflect the actual characteristics of spatial distribution. In addition, there are few available, suitable, and high spatial resolution PM 2.5 products, so we use daily MODIS images and introduce the EDTA method to retrieve high-resolution AODs in space and time to create high-precision and highresolution PM 2.5 concentrations at a 1-km grid level based on a seasonal spatial-temporal estimation model. In this work, we addressed problems of the spatial resolution and the spatial-temporal continuity of PM 2.5 concentration data.
Previous studies that adopt spatial distribution analysis of PM 2.5 using monitoring station data and GIS interpolation methods only account for local area situations or analyze variations in PM 2.5 characteristics and time trends using annual average data. Studies using monitoring stations to obtain PM 2.5 did not take suitable spatial resolution into account; studies using interpolation methods to derive PM 2.5 did not focus on seasonal differences in PM 2.5 pollution. Both are not precise enough in time and space.
By introducing the EDTA method and seasonal spatial-temporal models, we estimate the seasonal PM 2.5 concentration in a high spatial resolution of 1 km × 1 km. The verification accuracy (R 2 ) of estimation PM 2.5 concentration reached 0.513 in spring, 0.640 in summer, 0.520 in autumn, and 0.540 in winter, and the estimation error (RMSE) was in the range of 3.979~7.893 µg/m 3 . We showed the feasibility and reliability of retrieving AOD and estimating PM 2.5 from MODIS remote sensing images. According to the characteristics of different seasons, we also constructed a corresponding estimation model that has seasonal applicability. Today, 1 km-grid level of PM 2.5 concentration is a higher-resolution data source, which could allow for assessing fine PM 2.5 pollution at small and medium scale. These data sets and research results are useful for policymakers in air pollution control administration to plan a more sustainable living environment. SHB is seriously affected by the southeast monsoon and plum rain season in summer. Despite the high time resolution of the MODIS images (daily), these images are often influenced by massive cloud cover in AOD inversion and PM 2.5 estimation. leading to poor estimated results on some days. Therefore, the replacement or supplement of these data (seriously disturbed by clouds) will be helpful for improving the accuracy of AOD inversion and PM 2.5 estimation. While the inversion AOD accuracy is high by the method of EDTA, the basic DTA method has certain requirements for the vegetation coverage in the region. In high surface reflectance areas, this method is only moderately effective.
Combining more suitable inversion methods according to different surface types could increase the accuracy.
Future investigations would be helpful in the following aspects: changes in particulate matter at a micro scale would improve the accuracy of AOD-PM 2.5 estimation models; problems of data loss and signal-noise ratio (caused by cloud and rainfall interference) could be anticipated in the development of multi-source spatial-temporal data fusion; for instance, we are trying to solve the single-source data in the product resolution and spatial-temporal coverage and improve the accuracy and reliability of estimated PM 2.5 data; improved estimation PM 2.5 data could be used to evaluate the population exposure risk to PM 2.5 , health economic losses, and early warning and prevention, in order to provide scientific reference for policymakers for improving urban atmospheric pollution and living environments.

Conclusions
This study developed a framework to improve the spatial resolution of AOD and PM 2.5 dataset and present the health risk assessment from PM 2.5 pollution: first, we retrieve the daily AOD using MODIS remote sensing images, AERONET AOD data, and the Enhanced Dark Target Algorithm (EDTA). Then we mapped the monthly and seasonal AOD results at a spatial resolution of 1 km × 1 km. Second, we inferred the optimal relationship between retrieved AOD and observed PM 2.5 in four seasons, and spatial-temporal seasonal models were developed to estimate PM 2.5 concentration. According to geographical features and seasonal characteristics in the study area, we obtain seasonal PM 2.5 at 1-km grid level by GIS platform. Third, we assessed health risk from PM 2.5 pollution using highaccuracy population density data and the relative exposure risk model. Last, the usage and deficiency of PM 2.5 dataset and risk assessment results were discussed. Therefore, reasonable assessments on health risk from PM 2.5 pollution are important for improving public health and living environment.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
Data is available on request to the corresponding author.

Acknowledgments:
We thank all the participants involved in the project for their contribution to our research data.

Conflicts of Interest:
The authors declare no conflict of interest.