Combining Data from Multiple Sources to Evaluate Spatial Variations in the Economic Costs of PM2.5-Related Health Conditions in the Beijing–Tianjin–Hebei Region

Fine particulate matter, known as PM2.5, is closely related to a range of adverse health outcomes and ultimately imposes a high economic cost on the society. While we know that the costs associated with PM2.5-related health outcomes are not uniform geographically, a few researchers have considered the geographical variations in these costs because of a lack of high-resolution data for PM2.5 and population density. Satellite remote sensing provides highly precise, high-resolution data about how PM2.5 and population density vary spatially, which can be used to support detailed health-related assessments. In this study, we used high-resolution PM2.5 concentration and population density based on remote sensing data to assess the effects of PM2.5 on human health and the related economic costs in the Beijing–Tianjin–Hebei (BTH) region in 2016 using exposure-response functions and the relationship between health and economic costs. The results showed that the PM2.5-related economic costs were unevenly distributed and as with the population density, the costs were mainly concentrated in urban areas. In 2016, the economic costs of PM2.5-related health endpoints amounted to 4.47% of the total gross domestic product in the BTH region. Of the health endpoints, the cost incurred by premature deaths accounted for more than 80% of the total economic costs associated with PM2.5. The results of this study provide new and detailed information that could be used to support the implementation of national and regional policies to reduce air pollution.


Introduction
Air pollution has significant impacts on health. The World Health Organization (WHO) estimated that, worldwide, three million people died prematurely because of air pollution in 2012 [1]. As a major air pollutant, fine particulate matter with an aerodynamic diameter of ≤2.5 µm, known as PM 2.5 [2,3], remains suspended in the air for a long time. With a small particle size and large surface area, it is very active and carries various poisonous substances into the body, mainly through the respiratory tract, and penetrates deep into the tracheobronchial and alveolar regions [4,5]. Exposure to PM 2.5 is associated with various health conditions such as respiratory diseases [6], cardiovascular diseases [7,8], lung cancer [9][10][11], nervous system diseases [12], and congenital heart defects [13]. These negative effects on health inevitably increase the economic costs of health, which then increase the burden on economic development [14]. Therefore, it is important to assess how air pollution impacts on society and economic development and ascertain the health-related economic costs of PM 2.5 .
To date, researchers have estimated the economic cost of health impacts caused by PM 2.5 in several cities and countries worldwide [15][16][17]. The economic loss caused by PM 2.5 was estimated to account for 4.7% of the GDP in the Beijing-Tianjin-Hebei (BTH) area of China in 2009 [18]. The health-related modeled from the DMSP/OLS NTL data. Then, the spatial variations in the health-related economic costs of PM 2.5 in the BTH region were assessed quantitatively by using the spatially highly resolved remote sensing data, the exposure-response function, and the health economic-loss relationship adjusted to the BTH region. The results provide new information that can be used to formulate effective measures to reduce the economic burden associated with air pollution and, at the same time, offer effective technological support for large-scale social activities such as the Winter Olympic Games.

Fine Particulate Matter (PM 2.5 ) Data
To date, a variety of empirical models [36,37] and advanced models [28,30,38] have been developed to estimate PM 2.5 concentrations from satellite-derived AOD products for different parts of the world. In this study, we estimated the PM 2.5 concentrations of the BTH region in 2016 at a 3-km resolution using multiple remote sensing images and the linear mixed-effects (LME) model [31,39]. The LME model can be used to account for daily variations in the AOD-PM 2.5 relationship. Daily variations in the relationship between additional meteorological variables and PM 2.5 were not considered in our model because they led to overfitting of the model [31]. The model structure was expressed as: where PM s,t (µg/m 3 ) is the averaged observed PM 2.5 concentration at a grid cell s on day t; b 0 and b 0,t are the fixed and day-specific random intercepts, respectively; b i (i = 1, 2, . . . . . . ) represents the fixed slopes for the predictors; b 1,t is the random slope for AOD for each day; AOD s,t refers to the MODIS AOD value (unitless) at grid cell s on day t; and PBLH s,t , RH s,t , WS s,t , and PS s,t are the meteorological parameters at grid cell s on day t. The meteorological parameters included the planetary boundary layer height (PBLH, m), relative humility (RH, %), wind speed (WS, m/s), and surface pressure (PS, hPa). The data used in the model included hourly PM 2.5 concentrations from ground monitoring sites ( Figure 1) and the MODIS 3-km AOD and assimilated meteorological data. We obtained PM 2.5 data for the BTH region for 2016 from China's Air Quality Online Monitoring and Analysis Platform [40] and the Beijing Municipal Environmental Monitoring Center [41]. We also obtained MODIS Collection 6 Level 2 aerosol products with a spatial resolution of 3 km for 2016 [42]. The data from the 3-km AOD products agreed well with the ground-based data from six Aerosol Robotic Network (AERONET) sites in China [43]. In this study, Aqua (overpass at~13:30 LocalTime) and Terra (overpass at 10:30 LocalTime) MODIS AOD data at 550 nm were first combined and averaged to improve the spatial coverage [31]. The 2016 meteorological data were obtained from the Goddard Earth Observing System Data Assimilation Forward Processing (GEOS-5 FP) system. The NASA Global Modeling and ASSimilation Office (GMAO) GEOS-5 FP products, available since April 2012, is the latest version of the GEOS-5 meteorological data, and has a spatial resolution of 0.25 • (latitude) × 0.3125 • (longitude) in a nested grid that covers China [44]. The hourly PM 2.5 concentrations and the GEOS-5 FP data from 10:00 to 14:00 local time were averaged as the daytime PM 2.5 concentrations and meteorological parameters to match the overpass times of the MODIS Terra and Aqua satellites. All data were reprojected using the Albers Equal Area Conic Projection and were resampled to a spatial resolution of 3 km. The data were clipped to match the administrative boundary of the study area. Finally, all the variables for 2016 were matched by grid cell and date for model fitting.
We used the coefficient of determination (R 2 ) to indicate the degree of agreement between the predicted and observed PM 2.5 concentrations. We also applied a 10-fold cross-validation (CV) method to check for potential model overfitting [45] and used R 2 as an indicator of the agreement. 0.3125° (longitude) in a nested grid that covers China [44]. The hourly PM2.5 concentrations and the GEOS-5 FP data from 10:00 to 14:00 local time were averaged as the daytime PM2.5 concentrations and meteorological parameters to match the overpass times of the MODIS Terra and Aqua satellites. All data were reprojected using the Albers Equal Area Conic Projection and were resampled to a spatial resolution of 3 km. The data were clipped to match the administrative boundary of the study area. Finally, all the variables for 2016 were matched by grid cell and date for model fitting. We used the coefficient of determination (R 2 ) to indicate the degree of agreement between the predicted and observed PM2.5 concentrations. We also applied a 10-fold cross-validation (CV) method to check for potential model overfitting [45] and used R 2 as an indicator of the agreement.

The Spatial Distribution of the Population Density
In this study, the spatial distribution of population density was determined with the spatial lag regression model (SLM) based on DMSP/OLS NTL and land cover data. Traditional regression models assume that the observation is mutually independent, which is not valid because of the spatial structure of the data [46]. The SLM method estimates the spatial dependency in regression analysis and provides information about the spatial relationships among the parameters involved in the model [47]. The population density model is expressed as: where POPi denotes the initial population of the ith pixel; DNi denotes the digital number (DN) of the DMSP/OLS NTL image; Lik denotes the kth land class area of the ith pixel; ρ is the regression coefficient of the spatial lag variable ωPOPi; u and ak are the regression coefficients

The Spatial Distribution of the Population Density
In this study, the spatial distribution of population density was determined with the spatial lag regression model (SLM) based on DMSP/OLS NTL and land cover data. Traditional regression models assume that the observation is mutually independent, which is not valid because of the spatial structure of the data [46]. The SLM method estimates the spatial dependency in regression analysis and provides information about the spatial relationships among the parameters involved in the model [47]. The population density model is expressed as: where POP i denotes the initial population of the ith pixel; DN i denotes the digital number (DN) of the DMSP/OLS NTL image; L ik denotes the kth land class area of the ith pixel; ρ is the regression coefficient of the spatial lag variable ωPOP i ; u and a k are the regression coefficients of DN i and L ik , respectively; ω denotes the spatial weight matrix; and µ denotes the vector of error terms that are independent, but not necessarily identically distributed. Because there were errors associated with the polynomial regression function, there were discrepancies between the initial population estimate and the census data. A proportionality coefficient was built to reduce the errors [48][49][50], as follows: where k is the proportionality coefficient of the BTH region, Census is all the census data for the BTH region, and POP is the modeled initial total population of the BTH region. The modeled initial population of all the pixels was multiplied by the proportionality coefficient, and the spatialized population of each pixel was obtained, as follows: where EstPOP i is the spatialized population of the ith pixel. As in a previous study [51], the mean relative error (MRE) was adopted to evaluate the overall accuracy of the population density estimates as follows: where γ j is relative error of the jth county, POP j is the estimated population of the jth county, POP j is the census data of the jth county, and n is the number of counties. The data sources included three aspects, namely, DMSP/OLS stable NTL data, MODIS land use data, and county-level data from the 6th Census ( Figure 2). Because the data from the 6th Census were collected in 2010, we also used NTL data and land use data from 2010. The DMSP/OLS stable NTL data (1-km resolution) values represent the average of all light intensity values observed on cloud-free days during a calendar year from temporally stable sources such as cities, towns, villages, industrial sites, and persistent gas flares [52]; ephemeral events (e.g., fires) were discarded. The DN values of the stable NTL data ranged from 1 to 63 and the background and noise were identified and replaced with values of zero [53]. For land cover, we used the MODIS MCD12Q1 products with a 500-m spatial resolution and the International Geosphere Biosphere Program (IGBP) global vegetation classification scheme. We adopted the IGBP scheme as it is the primary land cover classification scheme used in data from the Land Processes Distribution Active Archive Center of the United States Geological Survey (USGS) [54]. The land cover types of the MCD12Q1 products were amalgamated into five land cover classes, namely, urban land (built-up), forest, grass, water, and bare land (soil) [55,56]. The DMSP/OLS stable NTL data and MCD12Q1 from 2010 were reprojected to the Albers Conical Equal Area projection and were then co-registered and resampled to a pixel size of 1 km.

PM2.5 Health Risk Assessment
We compared various studies published in China and elsewhere in recent years that reported the health impacts caused by exposure to PM2.5 [18,24,[59][60][61][62][63][64] and established the exposure-response correlation between PM2.5 pollution and health effects in the BTH region. We estimated the PM2.5 health impacts in this study from mortality data, hospital admission rates for respiratory and cardiovascular diseases, outpatient visits to internal medicine and pediatric practitioners, and co-occurrences of chronic bronchitis, acute bronchitis, and asthma.
Disease and death in the population are probabilistic events that correspond to the Poisson distribution. The linear exposure-response functions (Equations (7) and (8)) were used [65,66]: Highly accurate population density data of 2010 with a spatial resolution of 1 km were then produced. The population density distribution in 2016 was produced using the annual population growth rate of the cities from 2010 [57,58].

PM 2.5 Health Risk Assessment
We compared various studies published in China and elsewhere in recent years that reported the health impacts caused by exposure to PM 2.5 [18,24,[59][60][61][62][63][64] and established the exposure-response correlation between PM 2.5 pollution and health effects in the BTH region. We estimated the PM 2.5 health impacts in this study from mortality data, hospital admission rates for respiratory and cardiovascular diseases, outpatient visits to internal medicine and pediatric practitioners, and co-occurrences of chronic bronchitis, acute bronchitis, and asthma.
Disease and death in the population are probabilistic events that correspond to the Poisson distribution. The linear exposure-response functions (Equations (7) and (8)) were used [65,66]: where β i is the exposure-response coefficient of health endpoint i, C is the ambient concentration of PM 2.5 (µg/m 3 ), C 0 is the threshold level at which no health effects are yet assumed (µg/m 3 ), E i is the actual morbidity or mortality of the health endpoint i under pollution level C, E 0i is the morbidity or mortality of a certain health endpoint i at the threshold concentration level, HI ij refers to the health impact i in city j under pollution level C, and Pop j is the exposed population (person) in city j. In this study, considering the potential of results in China, we used a baseline concentration (C 0 ) of 35 µg/m 3 , the threshold suggested by China's Environmental Protection Ministry [18,24,67,68]. The concentration-response coefficients for short-term mortality and morbidity are based on the previously published values listed in Table 1. To maximize the regional accuracy, we prioritized studies carried out in China when developing the exposure-response parameters. The mortality, respiratory hospital admission, and cardiovascular hospital admission parameters mainly refer to the meta-analysis results from Xie et al. (2009) [63] and Aunan and Pan (2004) [61]. The outpatient visits to pediatrics and internal medicine data are from the research about association of air pollution with hospital outpatient visits in Beijing [59]. The data of acute bronchitis, chronic bronchitis, and asthma attack are from the relevant domestic research. Table 2 shows the mortality and morbidity for the health endpoints in the BTH region, making reference mainly to the statistical and research results in the BTH region or across the nation. This information was used to calculate the patient and excess death numbers at health effect endpoints caused by PM 2.5 .

Quantitative Estimation of the Economic Costs
In this study, the value of a statistical life (VOSL) and the cost of illness (COI) methods were used to estimate the external costs, based on the health impacts. The VOSL is the value of a small change in the risk associated with the death of an unnamed member of a large group [73]. The VOSL represents society's collective willingness to pay (WTP) [74] for a small reduction in the annual mortality risk of death [75]. In this study, the VOSL method was applied to assess the economic losses in the context of excess deaths and chronic diseases caused by PM 2.5 [65]. We used data from a contingent valuation (CV) method conducted in Chongqing [65,70,75,76], where the VOSL of a local resident was about USD 34,458 and could increase by USD 14,434 for an annual income increase of USD 144.58. The VOSL of residents in the BTH region in 2016 was calculated as follows: where VOSL j and VOSL Cq are the VOSL of the city j and Chongqing, respectively; I j and I Cq are the personal income of the city j and Chongqing in 2016, respectively; and e is the elastic coefficient of WTP, which was assumed to be 1.0 in this study. The cost of chronic bronchitis was estimated as 0.055 of the VOSL, as suggested by Hammitt and Zhou (2006) [77]. The COI method calculates the cost of a disease from the costs of medical treatment and hospitalization, and productivity loss [78]. Based on the health effects calculated above, the associated economic costs were obtained using the following equation: where EC ij is the economic loss from health impact i in city j and Cost ij is its economic cost per case. The values used for the economic costs per case of different health impacts and the literature on which the estimates were based are provided in Table 3. The exchange rate in 2016 was 6.6423 CNY to 1 USD [79].

Spatial Variations in PM 2.5 and Population Density
For model fitting, the overall value of R 2 between the predicted and observed PM 2.5 concentrations was 0.7495, and the Root Mean Squard Error (RMSE) and Maximum Permissible Error (MPE) were 20.91 µg/m 3 and 14.55, respectively. The CV of the R 2 value of the model was 0.7077, which was 0.04 less than that for the model fitting results. The CV of the RMSE value for the model was about 1 µg/m 3 higher than the model fitting results, which suggests that the model was not substantially overfitted.
The spatial distribution of PM 2.5 and the mean distributions of PM 2.5 for 2016 for two seasons across the BTH region are shown in Figure 3. Because there were limited satellite data for winter (January, February, and December), the data were grouped into two seasons, namely, warm for the period from 15 April to 14 October and cold for the period from 15 October to 14 April. The blank areas in Tianjin in Figure 3 reflect that there is no predicted PM 2.5 due to the lack of satellite AOD data. Spatially, there was a southeast-to-northwest gradient in the predicted concentrations of PM 2.5 . The concentrations were high in Beijing and in the south of Hebei province, and were low in the north of Hebei province. There were distinct seasonal variations in PM 2.5 in the BTH region, and, apart from Beijing, the concentrations were the highest in the cold season and the lowest in the warm season. The abnormal seasonal variation in Beijing areas is due to the lack of satellite AOD. Both Terra and Aqua AOD data are not available in some days in cold season, especially in the urban center areas of Beijing. The predicted PM 2.5 concentrations in the warm season and the cold season ranged from 8 to 190 µg/m 3 , and from 2 to 300 µg/m 3 , respectively. This seasonal variation reflects both the meteorological conditions and local emission sources.  The MRE calculated for the validation population result was about 34.62%, which proved that the estimates of the population density were relatively high. The simulated population density is presented in Figure 4. The BTH region is one of the most densely populated regions in China, with population densities exceeding 500 people/km 2 . The population densities were very obviously higher in the urban areas (generally exceeded 1000 people/km 2 ) than in the rural areas (generally less than 250 people/km 2 ). The MRE calculated for the validation population result was about 34.62%, which proved that the estimates of the population density were relatively high. The simulated population density is presented in Figure 4. The BTH region is one of the most densely populated regions in China, with population densities exceeding 500 people/km 2 . The population densities were very obviously higher in the urban areas (generally exceeded 1000 people/km 2 ) than in the rural areas (generally less than 250 people/km 2 ).
The MRE calculated for the validation population result was about 34.62%, which proved that the estimates of the population density were relatively high. The simulated population density is presented in Figure 4. The BTH region is one of the most densely populated regions in China, with population densities exceeding 500 people/km 2 . The population densities were very obviously higher in the urban areas (generally exceeded 1000 people/km 2 ) than in the rural areas (generally less than 250 people/km 2 ).

Health Effects Linked to PM 2.5
The health impacts calculated from a baseline concentration of 35 µg/m 3 are shown in Table 4. Our results suggest that a PM 2.5 pollution event of this magnitude would result in 89,005 acute deaths from all causes in the BTH region, with 15,934 deaths in Beijing, 11,004 deaths in Tianjin, and 62,067 deaths in Hebei. Of the health endpoints, the number of cases of outpatient visits to internal medicine was the highest in Beijing and Tianjin, followed by acute bronchitis and outpatient visits to pediatrics. For Hebei, the number of cases of acute bronchitis was the highest of all health endpoints, followed by outpatient visits to internal medicine and pediatrics. These effects probably reflect both the PM 2.5 pollution and population distribution. While the mean PM 2.5 concentrations of the three regions were similar (about 70 µg/m 3 ), the number of cases for each health impact was much greater in Hebei than in Beijing and Tianjin, because the total population of Hebei (about 7.47 million) was much larger than the population of either Beijing (about 2.173 million) or Tianjin (about 1.562 million). Generally, as the PM 2.5 pollution and population increased, the health effects in the BTH region also increased.
The spatial distribution of the calculated health impacts of PM 2.5 in the BTH region are shown in Figure 5. The blank areas of the figure represent the areas where the PM 2.5 concentration is lower than 35 µg/m 3 and there are no health impacts. From the figure, we can see there was considerable spatial heterogeneity, and that more residents were affected by the PM 2.5 pollution in central Beijing, Tianjin, and the southern districts of Hebei. The health impacts were also directly proportional to the population size. In the cold season, PM 2.5 pollution was high to the northwest of Zhangjiakou and Chengde in the Hebei province (Figure 3), but the health impacts in these districts were smaller than those in central Beijing because of the differences in population distribution across these districts (Figure 4). The spatial patterns of health impacts in different seasons, based on the same population density (shown in Figure 5a-c), were similar to the patterns in the PM 2.5 concentrations (Figure 3).

Economic Costs of PM2.5
The economic costs calculated for all health endpoints and all-cause mortality for a baseline PM2.5 concentration of 35 µg/m 3 are shown in Table 5. Of the three areas, the external costs of PM2.5 were the highest in Hebei and lowest in Tianjin, where it was around 23% of the cost in Hebei. Economic losses due to all-cause mortality accounted for about 81% of the total external costs of PM2.5 pollution. Chronic bronchitis accounted for most of the cost of morbidity health impacts.
The ratio of the economic cost to GDP (Ecal/GDP) describes the economic burden of the health impacts linked to PM2.5 pollution and thus, it is an important parameter for formulating effective policies to maintain sustainable economic development [76]. The current economic development means that the values of Ecal/GDP are relatively high for this study. The economic cost of the health effects related to PM2.5 pollution in Beijing was 15.16 billion USD, around 3.9% of the GDP, while the economic costs of the health effects of PM2.5 pollution were 2.5% and 6% of the GDP in Tianjin and Hebei, respectively. These results suggest that the economic cost of health impacts related to PM2.5 has increased the economic burden in the BTH region.

Economic Costs of PM 2.5
The economic costs calculated for all health endpoints and all-cause mortality for a baseline PM 2.5 concentration of 35 µg/m 3 are shown in Table 5. Of the three areas, the external costs of PM 2.5 were the highest in Hebei and lowest in Tianjin, where it was around 23% of the cost in Hebei. Economic losses due to all-cause mortality accounted for about 81% of the total external costs of PM 2.5 pollution. Chronic bronchitis accounted for most of the cost of morbidity health impacts.
The ratio of the economic cost to GDP (Ecal/GDP) describes the economic burden of the health impacts linked to PM 2.5 pollution and thus, it is an important parameter for formulating effective policies to maintain sustainable economic development [76]. The current economic development means that the values of Ecal/GDP are relatively high for this study. The economic cost of the health effects related to PM 2.5 pollution in Beijing was 15.16 billion USD, around 3.9% of the GDP, while the economic costs of the health effects of PM 2.5 pollution were 2.5% and 6% of the GDP in Tianjin and Hebei, respectively. These results suggest that the economic cost of health impacts related to PM 2.5 has increased the economic burden in the BTH region. The spatial distribution of the total economic cost associated with PM 2.5 in 2016 and changes in the spatial distribution over the seasons are presented in Figure 6. The blank areas of the figure represent the areas where there was no economic cost linked to PM 2.5 . The blank areas in Tianjin in Figure 6b reflect the lack of PM 2.5 data for cold season, which we have not discussed. The economic losses reflected the population density ( Figure 4) and were concentrated in the urban areas and limited in the suburbs. The urban areas of Beijing, Tianjin, Shijiazhuang, Baoding, Xingtai, and Handan had high health-related economic losses because of PM 2.5 . The ranking was similar in different seasons and highly impacted areas were concentrated in the south of the BTH region, while the lower-ranked areas were mainly concentrated in the northern part of the BTH region. In the north of Hebei province, there was no economic cost in the warm season and the economic costs were low in the cold season (blank areas in Figure 6). In Chengde and Qinhuangdao in the southeast, the economic costs because of PM 2.5 were lower in the cold season than in the warm season. The seasonal variation in the spatial distribution of the economic costs associated with PM 2.5 mainly reflects the seasonal variations in PM 2.5 . In general, PM 2.5 pollution and population distribution together account for the variable economic costs, and as the PM 2.5 and population density increased, the economic costs also increased. areas in Tianjin in Figure 6(b) reflect the lack of PM2.5 data for cold season, which we have not discussed. The economic losses reflected the population density ( Figure 4) and were concentrated in the urban areas and limited in the suburbs. The urban areas of Beijing, Tianjin, Shijiazhuang, Baoding, Xingtai, and Handan had high health-related economic losses because of PM2.5. The ranking was similar in different seasons and highly impacted areas were concentrated in the south of the BTH region, while the lower-ranked areas were mainly concentrated in the northern part of the BTH region. In the north of Hebei province, there was no economic cost in the warm season and the economic costs were low in the cold season (blank areas in Figure 6). In Chengde and Qinhuangdao in the southeast, the economic costs because of PM2.5 were lower in the cold season than in the warm season. The seasonal variation in the spatial distribution of the economic costs associated with PM2.5 mainly reflects the seasonal variations in PM2.5. In general, PM2.5 pollution and population distribution together account for the variable economic costs, and as the PM2.5 and population density increased, the economic costs also increased.

Discussion
We assessed the spatial variation in the economic costs of PM2.5-related health impacts in 2016 in the BTH region. Overall, the annual economic cost of these illnesses in the BTH region is high, especially in the central and southern areas of Beijing, and in the urban areas of

Discussion
We assessed the spatial variation in the economic costs of PM 2.5 -related health impacts in 2016 in the BTH region. Overall, the annual economic cost of these illnesses in the BTH region is high, especially in the central and southern areas of Beijing, and in the urban areas of Tianjin, Shijiazhuang, and Baoding. The high costs in these areas reflect pollution emissions, meteorological conditions, and population density. The areas with the most difference between the cold season and warm season were concentrated in Zhangjiakou, Chengde, and the southern parts of Hebei. The variability in the spatial distribution between the seasons reflects the seasonal variations in the PM 2.5 concentrations.
The results of this study show that more than 2.3 million people suffered from health impacts related to PM 2.5 . These health impacts resulted in more than 15,000 premature deaths, representing about 0.07% of the 21.73 million residents in Beijing in 2016. In Tianjin, 1.6 million people suffered health impacts and 11,000 or 0.07% of the 15.62 million residents died prematurely. In Hebei, 7.8 million people suffered health impacts and 62,000 or 0.08% of the 74.70 million residents died prematurely. The most common health endpoints were outpatient visits to internal medicine and acute bronchitis. These results are generally consistent with those of previous studies; for example, Xie et al. (2011) reported that about 0.08% of the total population in the Pearl River Delta region in China died prematurely from PM 2.5 -related illnesses [80], while Apte et al. (2015) from their assessment of the health impacts of PM 2.5 on people of all ages in China found that premature deaths accounted for 0.09% of the total population in 2013 [81].
We calculated that the total external costs of PM 2.5 -related health impacts in Beijing, Tianjin, and Hebei were 15.16, 6.8, and 29 billion USD, which was equivalent to around 3.9%, 2.5%, and 6% of the GDP of these regions in 2016, respectively. These findings are similar to those from other studies. For example, Lv and Li (2016) reported that the economic losses from health impacts as a result of PM 2.5 pollution accounted for 2.16% (VSL) of the regional GDP in the BTH region in 2013 [67], while other researchers reported that the economic costs of PM 2.5 pollution accounted for 4.98% and 3.79% of Beijing's GDP in 2008 and 2012, respectively [66]. These estimates of the economic losses may vary depending on the PM 2.5 concentration levels, baseline levels, health impact assessment models, exposure-response coefficients, and the choice of health endpoints.
We used high-resolution population data and satellite-derived PM 2.5 data to explore the spatial distribution of economic costs from PM 2.5 -related health impacts. Because of the limited availability of high-resolution PM 2.5 and population data, previous studies only estimated the economic costs at the district scale. However, the population density and PM 2.5 concentrations vary considerably within cities and districts. Our study, therefore, is significant because we estimated the district-level and within-district economic costs related to air pollution in the BTH region. Our results build on those obtained from previous studies and provide more detailed information for policymakers about areas that should be given priority when developing strategies to control and prevent air pollution in this era of increased environmental protection in China.
However, there are some limitations in this analysis. First, it is difficult to consider all the possible health endpoints related to PM 2.5 and obtain the associated data and information [24]. Furthermore, the biological evidence for some health impacts is still not certain [82]. Hence, we did not include some important potential PM 2.5 health impacts such as diabetes and mental and behavioral disorders, which could impact on a person's ability to work or engage in routine daily activities, and could represent significant additional health and economic burdens on social services. Second, in our study, Aqua and Terra MODIS AOD data were combined to improve spatial coverage, but neither Terra nor Aqua AOD data are available for some days because of cloud contamination, especially in the cold season. The predicted PM 2.5 concentrations are affected by the missing values of the satellite AOD data, which require more accurate PM 2.5 spatial distributions in the study area. Third, the population density data also has its limitations. Although the estimates of the population density in 2010 were relatively high, the usage of annual population growth rate of the cities from 2010 to 2016 to produce the population density distribution in 2016 may introduce some errors, which may affect the results.
Fourth, we used a baseline PM 2.5 concentration of 35 µg/m 3 , as suggested by the Chinese Environmental Protection Ministry, and considered a sensitive parameter for both health impacts and estimating external costs [24]. In future studies, the WHO-recommended baseline of 10 µg/m 3 could also be used to estimate the external costs and the results calculated from both the baseline concentrations could be compared.
Despite these limitations, the results from this study can still be used to show the levels of air pollution and the district-level and within-district spatial variations of economic losses from PM 2.5 -related health impacts in the BTH region. However, to curb the health-related economic costs related to PM 2.5 , the PM 2.5 concentrations need to be reduced further.

Conclusions
We combined simulated high-resolution population density data and satellite-derived PM 2.5 concentration data and developed maps of the spatial distributions of economic losses related to multiple health impacts caused by PM 2.5 pollution in 2016 in the BTH region. The results showed that the PM 2.5 -related economic losses of health impacts in the BTH region were enormous, and were close to 4.47% of the BTH's regional GDP in 2016. We also found that the economic costs were unevenly distributed and as observed for population density, the costs were concentrated in urban areas and low in the suburbs. The losses were the greatest in Beijing, Tianjin, and cities in the southern Hebei province such as Shijiazhuang, Baoding, Xingtai, and Handan, and reflect the spatial distribution of the PM 2.5 pollution. Priority should be given to controlling the PM 2.5 pollution in these central and southern areas. There are some limitations associated with the economic losses estimated in this study because of data and research limitations. For instance, we did not include all PM 2.5 -related health impacts in the assessment and we used a static population exposure model in the study. These limitations should be addressed in future studies. Regardless of these limitations, the results from this study can still be used to enhance our fundamental understanding of the effects of PM 2.5 pollution on our health and economy. This study highlights important information that can be used to reduce the health risks and economic burdens associated with PM 2.5 . It should also be useful for policymakers and should support the design of new policies to prevent and control air pollution.