Industrial Heat Source-Related PM 2.5 Concentration Estimates and Analysis Using New Three-Stage Model in the Beijing–Tianjin–Hebei Region

: The prevalent high-energy, high-pollution and high-emission economic model has led to significant air pollution challenges in recent years. The industrial sector in the Beijing–Tianjin–Hebei (BTH) region is a notable source of atmospheric pollutants, with industrial heat sources (IHSs) being primary contributors to this pollution. Effectively managing emissions from these sources is pivotal for achieving air pollution control goals in the region. A new three-stage model using multi-source long-term data was proposed to estimate atmospheric, delicate particulate matter (PM 2.5 ) concentrations caused by IHS. In the first stage, a region-growing algorithm was used to identify the IHS radiation areas. In the second and third stages, based on a seasonal trend decomposition procedure based on Loess (STL), multiple linear regression, and U-convLSTM models, IHS-related PM 2.5 concentrations caused by meteorological and anthropogenic conditions were removed using long-term data from 2012 to 2021. Finally, this study analyzed the spatial and temporal variations in IHS-related PM 2.5 concentrations in the BTH region. The findings reveal that PM 2.5 concentrations in IHS radiation areas were higher than in background areas, with approximately 33.16% attributable to IHS activities. A decreasing trend in IHS-related PM 2.5 concentrations was observed. Seasonal and spatial analyses indicated higher concentrations in the industrially dense southern region, particularly during autumn and winter. Moreover, a case study in Handan’s She County demonstrated dynamic fluctuations in IHS-related PM 2.5 concentrations, with notable reductions during periods of industrial inactivity. Our results aligned closely with previous studies and actual IHS operations, showing strong positive correlations with related industrial indices. This study’s outcomes are theoretically and practically significant for understanding and addressing the regional air quality caused by IHSs, contributing positively to regional environmental quality improvement and sustainable industrial development.


Introduction
Particulate matter with a diameter of 2.5 µm or less, commonly called PM 2.5 , has significant adverse effects on the atmospheric environment, climate change, and human health [1][2][3].Industrial emissions from industrial heat sources (IHS) have been identified as one of the primary contributors to regional air pollution [4].Moreover, these emissions' distribution and transformation processes profoundly impact the spatial and temporal variations in PM 2.5 concentrations [5,6].Data show that the intensity of industrial pollution emissions in North China with dense IHSs is four times the national average, with atmospheric pollutant concentrations still exceeding the national average [7,8].Therefore, accurately estimating and analyzing the spatial and temporal variations in PM 2.5 concentrations caused by IHSs is crucial for optimizing emission control policies and combating air pollution.
Fuel consumption during industrial production causes many PM 2.5 emissions [9].Estimating the contribution of industrial emissions to PM 2.5 requires complex and timeconsuming emission source-apportionment studies [10].Within the last few decades, significant efforts have been made by researchers to quantify the impact of industry on PM 2.5 concentrations, including emission inventories [11,12], source apportionment [13,14], and numerical models [15,16], as well as statistical models [17,18].Many studies have predominantly aggregated industrial emission sources at a regional level, thereby assessing their collective impact on air quality [19,20].This approach is grounded in analyzing emission inventories and the statistical data delineated in annual reports.Duan et al. quantified the contribution rates of the iron and steel industry to PM 2.5 emissions in the BTH region based on the emission inventory specific to this sector, primarily utilizing data derived from field measurements and emission factors [21].Analogously, An et al. estimated the contribution of industrial sectors to PM 2.5 emissions in the Yangtze River Delta, utilizing an updated pollutant emission inventory and local measurement data [22].However, these studies have relied on large-scale data for the spatial distribution of emissions and need help to locate emission sources accurately [23].Source apportionment and numerical simulation methods elucidate the physical mechanisms behind pollutant formation and transport.These methodologies facilitate a nuanced assessment of PM 2.5 concentration contributions by various industries [24,25].Dong et al. utilized the positive matrix factorization (PMF) model to intricately dissect the origins of PM 2.5 in Changchun, primarily employing data from air quality monitoring stations [26].Their analysis revealed that secondary sources are the principal contributors to the city's PM 2.5 concentrations.Liu et al. employed the CAMx-PSAT model to assess the seasonal variations and sources of PM 2.5 in Beijing, identifying the industrial sector as one of the primary sources of PM 2.5 concentrations [27].However, such models have required extensive a priori knowledge and intricate parameterization.Furthermore, these models need more generalizability, making them difficult to extend [28].
Statistical models bypass the intricate atmospheric physical and chemical processes, focusing instead on analyzing pollutant concentrations based on the data's inherent features [29].Yuen conducted an analysis of the relationships between air quality and weather conditions using multiple regression techniques [30].Ding et al. applied the seasonal-trend decomposition procedure based on Loess (STL) method to isolate meteorological fluctuations, subsequently positing that changes in ozone levels in Tianjin are largely driven by variations in precursor emissions [31].However, traditional regression methods must be revised to address high-dimensional and nonlinear data, consequently constraining their performance.Therefore, these methods are frequently integrated with other models to enhance analytical capabilities [32].The advancement of deep learning has provided a novel perspective for monitoring and decomposing PM 2.5 concentrations.Deep learning methods, which integrate non-linear patterns and network models, outperform simple statistical models regarding fitting accuracy [33].Among various neural networks, long short-term memory (LSTM) networks, with their unique gate structures, have shown superiority in predicting time-series data [34].Concurrently, as popular deep learning models, convolutional neural networks (CNN) effectively extract spatial features from data [35].Therefore, some researchers have attempted to combine CNN with LSTM to study PM 2.5 concentrations from both temporal and spatial dimensions.Huang et al. proposed an APNet algorithm based on CNN-LSTM for predicting PM 2.5 concentrations in the upcoming hour, primarily utilizing data from air quality monitoring stations [36].Xu et al. developed a weight-sharing deep learning framework and multi-objective optimization network based on the CNN-LSTM model to predict PM 2.5 concentrations in Xi'an [33].The down-sampling structure in the U-Net model strengthens the extraction of spatial information, and its integration with CNN-LSTM has been shown to improve PM 2.5 forecasting

Study Area
The BTH region (36.0 • N-42.6 • N, 113.5 • E-119.8 • E) is situated on the North China Plain (Figure 1), encompassing Beijing, Tianjin, and Hebei province [42].The former two are municipalities, and Hebei Province comprises 11 prefecture-level cities.The Yan and Taihang mountain ranges in the north form a semi-enclosed structure around the region, resulting in low wind speeds and increased atmospheric stability.The BTH region hosts over 8% of China's population and is pivotal in driving economic growth in northern China [43].Additionally, this region stands as the most extensive energy industrial base and a significant steel hub in China.The region's complex topographical and adverse meteorological conditions accentuate severe atmospheric pollution caused by substantial energy consumption emissions.In 2022, the annual average PM 2.5 concentration in the BTH and surrounding areas was 44 µg/m 3 , exceeding the national average by 16 µg/m 3 .This persistent environmental issue has sustained attention from government entities, the media, and various societal sectors [44].

Data Source
The data required for this study encompass PM2.5 concentration data, natural geographic data, and socio-economic data.More information on the data can be seen in Table 1.

Data Source
The data required for this study encompass PM 2.5 concentration data, natural geographic data, and socio-economic data.More information on the data can be seen in Table 1.

PM 2.5 Concentration Data
The PM 2.5 concentration dataset was derived from the CHAP dataset (available at https://weijing-rs.github.io/product.html,accessed on 30 June 2023).The CHAP dataset is known for its encapsulation of long-term, full-coverage, high-resolution, and high-quality data on ground-level air pollutants within China [39][40][41].Many significant data sources, including ground-based measurements, satellite remote sensing products, atmospheric reanalysis, and model simulations, were incorporated to generate the dataset.It boasted a temporal resolution of one day and a spatial resolution of one kilometer.Additionally, the daily PM 2.5 estimates exhibit high accuracy, evidenced by an average cross-validation coefficient of determination of 0.92.Daily PM 2.5 data spanning from 2012 to 2021 were employed as the primary data source in this study.The focus was exploring the spatial and temporal trends of PM 2.5 concentrations in the BTH regions.

Natural Geographic Data
Meteorological data were obtained from the ECMWF Reanalysis v5 (ERA5) dataset, the fifth-generation atmospheric reanalysis product of global climate produced by the European Centre for Medium-Range Weather Forecasts (ECMWF) (available at https://earthengine. google.com/,accessed on 30 June 2023).Compared to its predecessor, the ERA-Interim, ERA5 provided atmospheric analysis data with superior temporal and spatial resolutions [45].The ERA5 DAILY data, downloaded from the Google Earth Engine (GEE) platform, comprised a temporal resolution of one day and a spatial resolution of 0.1 • × 0.1 • .This study extended the dataset period from 1 January 2012 to 31 December 2021.
The following parameters were selected to correlate meteorological variables with PM 2.5 concentrations: 2 m air temperature (T2M), total precipitation (TP), surface pressure (SP), 10 m u wind component (10U), and 10 m v wind component (10V).Relative humidity (RH) was then calculated using the t2m and 2 m dew point temperature.These six variables were crucial in isolating meteorological factors from PM 2.5 concentrations [46].
The digital elevation model (DEM) data were sourced from the Copernicus DEM, a product derived from global radar satellite observations (available at https://panda.copernicus.eu/panda,accessed on 30 June 2023) [47].These observations were processed through interferometric techniques to yield the TanDEM-X DEM product.Subsequently, a semi-automated editing and sampling process was employed to refine this data into a 30 m grid, resulting in the Copernicus 30 m resolution DEM dataset.

Socio-Economic Data
The dataset for IHS was derived from research conducted by Ma et al. (available at https://doi.org/10.57760/sciencedb.j00001.00430,accessed on 30 June 2023) [48].The k-means method was employed to identify heat source objects generated through the fire points' spatial density segmentation.Subsequently, a random forest algorithm was employed to identify and analyze IHSs.The results were validated and revised using Google Earth, other high-resolution satellite imagery, points of interest (POI) information, and field surveys.This process culminated in identifying 474 IHSs in the BTH region from 2012 to 2021.Furthermore, the dataset included attribute information on the operational status of the IHS, providing a precise reference for constructing the radiation area of the IHS.
The carbon emissions dataset was sourced from the China City Greenhouse Gas Working Group (CCG) (available at http://www.cityghg.com/toCauses?id=4, accessed on 30 June 2023).Drawing from internationally established and widely applied urban CO 2 emissions accounting methods and considering China's urban realities, the CCG calculated all direct and indirect emissions from cities purchasing electricity outside their administrative boundaries.Industrial indices, environmental emissions, and raw material production were sourced from the statistical yearbooks of Beijing, Tianjin, and Hebei provinces (available at https://www.gotohui.com,accessed on 30 June 2023).For this study, data from 2012 to 2021 were selected as key factors for measuring the correlation between IHS-related PM 2.5 concentrations and these indices.

Method
The technical route of this study is illustrated in Figure 2.This study proposed a new three-stage model using multi-source long-term data to estimate IHS-related PM 2.5 concentrations.In the first stage, the IHS radiation areas were established based on a region-growing algorithm using IHS and DEM data to identify the IHS radiation areas.In the second stage, pixel-by-pixel STL decomposition and multiple linear regression methods were employed to eliminate the influence of meteorological conditions.In the third stage, the U-ConvLSTM model was adopted to extract the background concentrations within the IHS radiation areas, thereby removing interference from other anthropogenic factors.Finally, the concentration and contribution of IHS-related PM 2.5 in the BTH Region were analyzed in three parts: the entire study area, all IHS radiation areas, and the typical individual IHS radiation area.The detailed steps of this procedure will be elaborated on in the following sections.

Data Preparation
PM2.5 concentration data and meteorological data were transformed to the WGS_1984 projection coordinate system in this study.Subsequently, the target area data were clipped based on the administrative boundaries of the BTH region.Due to variations in data sources, the temporal resolution of all datasets was harmonized to a daily scale, while the spatial resolution was consistently set to 1 km × 1 km.To ensure the integrity and reliability of the PM2.5 data, outliers and missing values were systematically excluded.Anomalies within the data were identified and removed using a quartile-based method.Lastly, the Kriging interpolation technique was employed to fill in the gaps in PM2.5 concentration and meteorological datasets.

Stage 1: Determining IHS Radiation Areas Based on Region-Growing Algorithm Using IHS and DEM Data
The IHS dataset provides spatial and temporal distributions of IHS [49].Emissions and pollution from these sources impact not just the site but also the surrounding areas [50].The region-growing algorithm using IHS and DEM data was implemented to define the IHS radiation areas [51].The detailed steps of the region-growing method are as follows: (1) The IHS data were transformed from vector to raster format.The DEM and IHS data were resampled to align with the PM2.5 concentration grid as a reference layer.(2) The region-growing algorithm was initiated from the central grid representing the IHSs, identified in Figure 3b.This grid's elevation was compared with the surrounding eight pixels.If the adjacent pixels displayed lower elevations, they were amalgamated into the IHS grid.

Data Preparation
PM 2.5 concentration data and meteorological data were transformed to the WGS_1984 projection coordinate system in this study.Subsequently, the target area data were clipped based on the administrative boundaries of the BTH region.Due to variations in data sources, the temporal resolution of all datasets was harmonized to a daily scale, while the spatial resolution was consistently set to 1 km × 1 km.To ensure the integrity and reliability of the PM 2.5 data, outliers and missing values were systematically excluded.Anomalies within the data were identified and removed using a quartile-based method.Lastly, the Kriging interpolation technique was employed to fill in the gaps in PM 2.5 concentration and meteorological datasets.

Stage 1: Determining IHS Radiation Areas Based on Region-Growing Algorithm Using IHS and DEM Data
The IHS dataset provides spatial and temporal distributions of IHS [49].Emissions and pollution from these sources impact not just the site but also the surrounding areas [50].The region-growing algorithm using IHS and DEM data was implemented to define the IHS radiation areas [51].The detailed steps of the region-growing method are as follows: (1) The IHS data were transformed from vector to raster format.The DEM and IHS data were resampled to align with the PM 2.5 concentration grid as a reference layer.(2) The region-growing algorithm was initiated from the central grid representing the IHSs, identified in Figure 3b.This grid's elevation was compared with the surrounding eight pixels.If the adjacent pixels displayed lower elevations, they were amalgamated into the IHS grid.(3) Repeat step (2) with all newly added grids serving as new growth points.The process persisted until either no lower elevation grids were identified or the predefined maximum impact distance was reached; at this point, the expansion ceased.Consequently, all incorporated grids were classified as the IHS radiation areas.The maximum distance for each IHS pollution impact was set to 10 km, according to the recommended method for the risk assessment of environmental emergencies in administrative regions (Figure 3c) [52].(4) Replicating steps 2 and 3 was repeated to delineate all IHS radiation areas of each IHSs, as shown in Figure 3a.
Atmosphere 2024, 15, x FOR PEER REVIEW 8 of 27 (3) Repeat step ( 2) with all newly added grids serving as new growth points.The process persisted until either no lower elevation grids were identified or the predefined maximum impact distance was reached; at this point, the expansion ceased.Consequently, all incorporated grids were classified as the IHS radiation areas.The maximum distance for each IHS pollution impact was set to 10 km, according to the recommended method for the risk assessment of environmental emergencies in administrative regions (Figure 3c) [52].(4) Replicating steps 2 and 3 was repeated to delineate all IHS radiation areas of each IHSs, as shown in Figure 3a.

Stage 2: Exclusion of Meteorological Factors Based on Pixel-by-Pixel STL Decomposition and Multiple Linear Regression
The relationship between meteorological elements and air quality exhibits regional variability attributable to the distinct topographical features and climatic conditions.To objectively estimate the impact of IHSs on PM2.5 concentrations, it is necessary to exclude the influence of meteorological factors on the short-term and long-term variation in pollutant concentrations.This study used STL decomposition to remove the components influenced by short-term meteorological changes.Concurrently, a multiple linear regression approach was employed to negate the influence of long-term meteorological conditions on the PM2.5 concentration series.The specific process of elimination is depicted in Figure 4 (using the value of a single pixel as an example).The specific steps are as follows:

Stage 2: Exclusion of Meteorological Factors Based on Pixel-by-Pixel STL Decomposition and Multiple Linear Regression
The relationship between meteorological elements and air quality exhibits regional variability attributable to the distinct topographical features and climatic conditions.To objectively estimate the impact of IHSs on PM 2.5 concentrations, it is necessary to exclude the influence of meteorological factors on the short-term and long-term variation in pollutant concentrations.This study used STL decomposition to remove the components influenced by short-term meteorological changes.Concurrently, a multiple linear regression approach was employed to negate the influence of long-term meteorological conditions on the PM 2.5 concentration series.The specific process of elimination is depicted in Figure 4 (using the value of a single pixel as an example).The specific steps are as follows:

Exclusion of Short-term Meteorological Factors Based on STL Decomposition
STL decomposition is a non-parametric statistical method that consists of an inner loop for fitting trend and seasonal components and an outer loop for calculating robustness weights to reduce the impact of outliers [53].Compared to other time-series decomposition methods, it offers greater precision and configurability.STL decomposition was utilized to separate the PM2.5 concentration of each grid pixel into short-term, seasonal, and long-term components while excluding short-term fluctuations present in the original data.The short-term component is primarily influenced by short-term meteorological changes and, to a lesser extent, by short-term variations in pollutant emissions [54].Seasonal and long-term components are mainly attributed to changes in human activity-related emissions of air pollutants and the seasonality and inter-annual variations in meteorological conditions.For the STL decomposition applied to our daily PM2.5 concentration dataset spanning from 2012 to 2021, the seasonal window was set to 365 days, aiming to accurately capture the annual fluctuation patterns of PM2.5 concentrations.After decomposition, the expression for the composition of atmospheric acceptable particulate matter concentrations for each pixel can be calculated as follows: where PM .(t), T(t), S(t), and R(t) represent the observed values of atmospheric particulate matter concentration series, long-term fluctuation, seasonal fluctuation, and shortterm fluctuation at time t, respectively, and N denotes the time series length.
The baseline concentration, obtained by excluding the short-term component, can be expressed as:

Exclusion of Short-term Meteorological Factors Based on STL Decomposition
STL decomposition is a non-parametric statistical method that consists of an inner loop for fitting trend and seasonal components and an outer loop for calculating robustness weights to reduce the impact of outliers [53].Compared to other time-series decomposition methods, it offers greater precision and configurability.STL decomposition was utilized to separate the PM 2.5 concentration of each grid pixel into short-term, seasonal, and long-term components while excluding short-term fluctuations present in the original data.The shortterm component is primarily influenced by short-term meteorological changes and, to a lesser extent, by short-term variations in pollutant emissions [54].Seasonal and long-term components are mainly attributed to changes in human activity-related emissions of air pollutants and the seasonality and inter-annual variations in meteorological conditions.For the STL decomposition applied to our daily PM 2.5 concentration dataset spanning from 2012 to 2021, the seasonal window was set to 365 days, aiming to accurately capture the annual fluctuation patterns of PM 2.5 concentrations.After decomposition, the expression for the composition of atmospheric acceptable particulate matter concentrations for each pixel can be calculated as follows: where PM 2.5 (t), T(t), S(t), and R(t) represent the observed values of atmospheric particu- late matter concentration series, long-term fluctuation, seasonal fluctuation, and short-term fluctuation at time t, respectively, and N denotes the time series length.
The baseline concentration, obtained by excluding the short-term component, can be expressed as: where PM 2.5bl (t),T(t)and S(t) represent the baseline concentration of atmospheric particu- late matter, long-term fluctuation, and seasonal fluctuation at time t, respectively, and N is the time series length.

Exclusion of Long-term Meteorological Factors Based on Multiple linear regression
PM 2.5bl (t) is typically associated with long-term variations in meteorological condi- tions and emission sources.The research has indicated that meteorological factors can account for 31.96% of PM 2.5 concentration variations in the BTH region [55].In particular, temperature and relative humidity have been found to exert the most significant influence on air quality in the area.At the same time, precipitation, wind direction, wind speed, and surface pressure also play a role to some extent.As a commonly used statistical method, multiple linear regression facilitates the quantitative elucidation of the relationship between meteorological variables and PM 2.5 concentrations.This methodology effectively discerns the extent of meteorological factors' impact on air quality [56].The advantage of multiple linear regression lies in accurately describing the relationships between variables by adjusting the coefficients of each independent variable.In this study, a multiple linear regression model was constructed to eliminate the influence of long-term meteorological conditions on the PM 2.5 concentration series.The specific implementation steps of model construction and analysis are as follows: The STL decomposition method was employed to break down the pre-processed meteorological parameters into long-term, seasonal, and short-term fluctuations.Shortterm fluctuations caused by meteorological changes and human activities were excluded, thereby retaining the baseline components of each meteorological parameter.
Multiple linear regression was performed for each meteorological variable using the daily averages or cumulative values of all meteorological elements, expressed as: where met bl i (t)(iϵ [1, 6]) represents the baseline component of meteorological variable i at t, with i = 1, 2, 3, 4, 5, 6 corresponding to t2m, tp, sp, 10u, 10v and rh, respectively.a i is the regression coefficient for the corresponding meteorological variable, a 0 is the intercept, and N represents the length of the time series.
The impacts of seasonal and inter-annual variations in meteorological conditions on atmospheric fine particulate matter concentrations were eliminated.The expression is as follows: PM 2.5bl tem (t) = PM2.5 bl (t) − PM2.5 bl mlr (t), t = 1, 2, . . ., N where PM 2.5bl tem (t) represents the concentration of PM 2.5 at t after excluding meteorologi- cal factors.

Stage 3: Exclusion of Other Anthropogenic Factors Based on U-ConvLSTM Model
In this research, the U-ConvLSTM model was employed to extract the spatial and temporal features of PM 2.5 concentrations influenced by other anthropogenic factors in background areas, aiming to predict the PM 2.5 concentrations in IHS radiation areas affected by these same factors.The exclusion process is shown in Figure 5, and the specific architecture of the model is depicted in Figure 5. Combining CNN with LSTM, the ConvL-STM model enhances traditional LSTM by incorporating convolutional operations within its structure [57].This innovation lets the model grasp the data's temporal associations and spatial details.Developed initially for precipitation forecasting, it has extensive application in meteorological predictions.In our study, the model effectively discerned and encoded the local spatiotemporal characteristics of particulate matter concentrations influenced by various non-meteorological factors.The essential equations related to ConvLSTM can be calculated as follows: where i t denotes the output of the input gate, f t denotes the output of the forget gate, C t represents the state of the cell at the current moment, o t corresponds to the gate's output, and H t denotes the predicted output of updating the current sequence index.
data [38].This feature is beneficial for identifying localized PM2.5 concentration characteristics within background regions and precisely predicting PM2.5 concentrations in IHS radiation areas influenced by other anthropogenic factors.The encoder comprised four layers, each incorporating a 3 × 3 convolution, 2 × 2 max pooling, and a ReLU activation function.Following downsampling, the encoder generated feature maps enriched with semantic information.The decoder then fused these acquired spatiotemporal features across various scales, employing upsampling and skip connections, ultimately producing the forecast through a 1 × 1 convolution.In this model, PM .concentration images were inputted, with the PM2.5 pixel concentrations in background areas considered to be influenced by anthropogenic factors.These data were the training set for predicting PM2.5 concentrations in IHS radiation areas impacted by similar factors.Ultimately, the expression for IHS-related PM2.5 concentrations was formulated as follows: where PM .represents the IHS-related PM2.5 concentrations, PM .represents the PM2.5 concentrations influenced by other anthropogenic factors.Furthermore, the encoder-decoder structure of the U-net maintained the spatial structure of images, facilitating the distinction and reconstruction of fine details in the data [38].This feature is beneficial for identifying localized PM 2.5 concentration characteristics within background regions and precisely predicting PM 2.5 concentrations in IHS radiation areas influenced by other anthropogenic factors.The encoder comprised four layers, each incorporating a 3 × 3 convolution, 2 × 2 max pooling, and a ReLU activation function.Following downsampling, the encoder generated feature maps enriched with semantic information.The decoder then fused these acquired spatiotemporal features across various scales, employing upsampling and skip connections, ultimately producing the forecast through a 1 × 1 convolution.In this model, PM 2.5bl tem concentration images were inputted, with the PM 2.5 pixel concentrations in background areas considered to be influenced by anthropogenic factors.These data were the training set for predicting PM 2.5 concentrations in IHS radiation areas impacted by similar factors.Ultimately, the expression for IHS-related PM 2.5 concentrations was formulated as follows: where PM 2.5ihs represents the IHS-related PM 2.5 concentrations, PM 2.5bl oem represents the PM 2.5 concentrations influenced by other anthropogenic factors.

Statistical Analysis of IHS-Related PM 2.5 Concentrations
This research utilized various methods, such as the Mann-Whitney U test, correlation analysis, the Mann-Kendall test, and spatial autocorrelation analysis, to examine the spatiotemporal features and trends of IHS-related PM 2.5 concentrations in the study region.The methodologies are detailed below: Mann-Whitney U test: A Mann-Whitney U test was conducted between the PM 2.5 concentrations in background areas and those in the IHS radiation areas to determine if IHS activities had a quantifiable impact on PM 2.5 levels.The Mann-Whitney U test, a non-parametric statistical method, compares the differences between two independent groups [58].This test is beneficial when the data deviate from a normal distribution as a reliable substitute for the independent samples t-test.In this study, the Mann-Whitney U statistic can be computed as follows: where n 1 and n 2 , respectively, denote the sample sizes of the PM 2.5 concentration data from the background area and the IHS radiation area.R 1 is the sum of the ranks for the background area.The lower the U statistic, the more substantial the evidence of a significant difference between the groups.Correlation analysis: Correlation analysis is a statistical method used to assess the strength and direction of the linear relationship between two or more variables [59].This study primarily employed correlation analysis to measure the association between IHSrelated PM 2.5 concentrations and indicators such as the number of operational IHSs and carbon emissions.The correlation is quantified using Pearson's correlation coefficient, which ranges from −1 to +1.A coefficient value closer to 1 or −1 indicates a stronger linear relationship between the variables.The Pearson correlation coefficient is calculated using the formula: where x i , y i represent the observed values of IHS-related PM 2.5 concentrations and other related variables, respectively, while x and y denote their mean values.
Mann-Kendall test: The Mann-Kendall test, a non-parametric statistical tool, is utilized in this study for trend analysis in data series [60].This method is apt for identifying monotonic trends, either increasing or decreasing, in time series data without necessitating a specific data distribution.This study employed the p-values from the Mann-Kendall test to ascertain the trend in the contribution of IHS-related PM 2.5 concentrations.It provides both the slope of the contribution change and its 95% confidence interval.A trend is deemed significant if the p-value is below 0.05.Additionally, a z-value above 0 indicates an upward trend, whereas a z-value below 0 indicates a downward trend.
Spatial autocorrelation analysis: Spatial autocorrelation analysis is a statistical approach used to evaluate whether the spatial distribution of a variable shows systematic clustering or scattering [61].This study applied spatial autocorrelation analysis to determine the spatial distribution patterns of IHS-related PM 2.5 concentrations.Typically, Moran's I index is employed for this analysis, serving as an indicator of spatial correlation.Moran's I index ranges from −1 to +1, with values near +1 indicating positive spatial correlation (clustering), those closer to −1 suggesting negative correlation (dispersion), and values around 0 implying a random distribution.The Moran's I index calculation is as follows: where n represents the number of pixels, x i and x j denote the IHS-related PM 2.5 concentration values at locations i and j, x is the average of all concentration values, and ω ij is the spatial weight between locations i and j.

Overall Spatial and Temporal
Trends of PM 2.5 in the BTH Region 3.1.1.Distribution and Variation of PM 2.5 Concentrations in the BTH Region Figure 6 illustrates the spatial and temporal distribution of PM 2.5 concentrations throughout the BTH region.From 2012 to 2021, the area displayed a spatial pattern of higher concentrations in the southern parts and lower in the northern regions.The mean concentration per pixel varied between 26.89 µg/m 3 and 93.28 µg/m 3 , with an overall average of 54.96 µg/m 3 , as depicted in Figure 6a.Spatial analysis revealed that the average PM 2.5 concentration in IHS radiation areas was 84.90 µg/m 3 , reaching a peak of 281.74 µg/m 3 and a standard deviation of 39.84 µg/m 3 .Conversely, the background areas showed an average PM 2.5 concentration of 21.97 µg/m 3 , characterized by lower variability (standard deviation of 18.41 µg/m 3 ) and a peak concentration of 145.61 µg/m 3 , with less frequent occurrences of high concentration levels compared to the IHS radiation areas, as shown in Figure 6b.Statistical analysis using the Mann-Whitney U test indicated a significant difference (p < 0.01) in the average PM 2.5 concentrations between IHS radiation and background areas, indicating a discernible impact of IHS operations on PM 2.5 concentrations.
Atmosphere 2024, 15, x FOR PEER REVIEW 13 of 27   7b).Correlational analysis, considering the annual operational status of IHS, revealed a correlation coefficient of 0.94, highlighting a significant link between IHS operations and PM2.5 levels (Figure 7c).When combining background and IHS contributions, approximately 33.16% of the total PM2.5 concentration in these areas was attributed to IHSs (Figure 7d).Further, a Mann-Kendall trend analysis  7b).Correlational analysis, considering the annual operational status of IHS, revealed a correlation coefficient of 0.94, highlighting a significant link between IHS operations and PM 2.5 levels (Figure 7c).When combining background and IHS contributions, approximately 33.16% of the total PM 2.5 concentration in these areas was attributed to IHSs (Figure 7d).Further, a Mann-Kendall trend analysis of the contributions of IHS-related PM 2.5 indicated a decreasing trend over time (p < 0.05), signifying a year-over-year reduction in the influence of IHS production activities on PM 2.5 concentrations and indicating some success in industrial reform efforts.

Overall Spatial and Temporal Trends of PM2.5 in the BTH Region
Atmosphere 2024, 15, x FOR PEER REVIEW 14 of 27 of the contributions of IHS-related PM2.5 indicated a decreasing trend over time (p < 0.05), signifying a year-over-year reduction in the influence of IHS production activities on PM2.5 concentrations and indicating some success in industrial reform efforts.

Temporal Trends and Seasonal Fluctuations on IHS-Related PM2.5 Concentrations
The BTH region experienced considerable variability and a degree of periodicity in overall IHS-related PM2.5 concentrations from 2012 to 2021 (Figure 8a).The Mann-Kendall trend test indicated a declining trend in the annual mean IHS-related PM2.5 concentrations (S = −0.007,p < 0.05), signifying statistical significance.
Regarding yearly variations, the peak average concentration occurred in 2013 at 42.42 μg/m 3 , while the lowest was recorded in 2021 at 18.14 μg/m 3 .The average concentration represented a 52.05% reduction from 2012 to 2021.Notably, there was a 7.2% increase in  The BTH region experienced considerable variability and a degree of periodicity in overall IHS-related PM 2.5 concentrations from 2012 to 2021 (Figure 8a).The Mann-Kendall trend test indicated a declining trend in the annual mean IHS-related PM 2.5 concentrations (S = −0.007,p < 0.05), signifying statistical significance.
2012, each season in 2021 showed improvements, with varying degrees of reduction (Figure 8d).Summer saw the most significant decrease at 64.20%, an annual average reduction of 1.89 μg/m 3 , followed by winter at 54.33% (3.40 μg/m 3 annually), autumn at 52.23% (2.07 μg/m 3 annually), and spring with the most minor decrease at 39.01% (1.45 μg/m 3 annually).In particular after 2014, reductions were noted in each season, potentially reflecting the impact of regional environmental policies.From 2012 to 2021, seasonal variations in IHS-related PM 2.5 concentrations in the BTH region showed a distinct pattern (Figure 8c).The trend displayed a 'low in springsummer and high in autumn-winter' pattern, with winter concentrations being the highest, followed by autumn, spring, and summer.Spring concentrations averaged 25.50 µg/m 3 , summer at 18.46 µg/m 3 , autumn at 28.77 µg/m 3 , and winter at 44.64 µg/m 3 .Compared to 2012, each season in 2021 showed improvements, with varying degrees of reduction (Figure 8d).Summer saw the most significant decrease at 64.20%, an annual average reduction of 1.89 µg/m 3 , followed by winter at 54.33% (3.40 µg/m 3 annually), autumn at 52.23% (2.07 µg/m 3 annually), and spring with the most minor decrease at 39.01% (1.45 µg/m 3 annually).In particular after 2014, reductions were noted in each season, potentially reflecting the impact of regional environmental policies.

Spatial Distribution and Analysis of IHS-Related PM 2.5 Concentrations
The spatial distribution of IHS-related PM 2.5 concentrations partially mirrors the industrial structure and the efficacy of environmental policies in distinct areas.Across the BTH region from 2012 to 2021, there was a consistent spatial trend of 'lower in the north, higher in the south' in IHS-related PM 2.5 concentrations, with a stable distribution over the study period.The Moran's I index for these concentrations stood at 0.98, reflecting a solid correlation with the density of industrial heat sources.
Notably, the pixel-level concentration map (Figure 9a) revealed that areas with the highest average PM 2.5 values, mainly in Tangshan and Handan, showed peak concentrations up to 78.40 µg/m 3 .These areas were known for their dense industrial heat sources.In contrast, regions with fewer IHSs showed lower PM 2.5 concentrations, indicative of reduced industrial emission densities.At the city level, the average IHS-related PM 2.5 concentrations varied, not aligning with IHS density, which can be attributed to different industrial pollution control standards in each city (Figure 9b).Handan had the highest average concentration (37.60 µg/m 3 ), followed by Cangzhou (37.45 µg/m 3 ), with Zhangjiakou presenting the lowest (18.16 µg/m 3 ).This pattern suggests that the southern parts of the BTH region are inclined towards industries with higher pollution and emissions, highlighting the need for more intensive industrial restructuring and enforcement of environmental regulations.The spatial distribution of IHS-related PM2.5 concentrations partially mirrors the industrial structure and the efficacy of environmental policies in distinct areas.Across the BTH region from 2012 to 2021, there was a consistent spatial trend of 'lower in the north, higher in the south' in IHS-related PM2.5 concentrations, with a stable distribution over the study period.The Moran's I index for these concentrations stood at 0.98, reflecting a solid correlation with the density of industrial heat sources.
Notably, the pixel-level concentration map (Figure 9a) revealed that areas with the highest average PM2.5 values, mainly in Tangshan and Handan, showed peak concentrations up to 78.40 μg/m 3 .These areas were known for their dense industrial heat sources.In contrast, regions with fewer IHSs showed lower PM2.5 concentrations, indicative of reduced industrial emission densities.At the city level, the average IHS-related PM2.5 concentrations varied, not aligning with IHS density, which can be attributed to different industrial pollution control standards in each city (Figure 9b).Handan had the highest average concentration (37.60 μg/m 3 ), followed by Cangzhou (37.45 μg/m 3 ), with Zhangjiakou presenting the lowest (18.16 μg/m 3 ).This pattern suggests that the southern parts of the BTH region are inclined towards industries with higher pollution and emissions, highlighting the need for more intensive industrial restructuring and enforcement of environmental regulations.

Heterogeneity and Evolution of IHS-Related PM 2.5 Concentrations in the BTH
As Figure 10 and Table 2 illustrate, the spatial and temporal evolution of IHS-related PM 2.5 concentrations in the 13 cities of the BTH region exhibit notable heterogeneity.In 2012, concentrations varied from 14.33 µg/m 3 to 55.14 µg/m 3 , with contributions ranging between 26% and 55% (Figure 10a, Table 2).Regions with active IHS clusters, such as Cangzhou, Tangshan, Tianjin, and Handan, had noticeably higher concentrations and contributions than other areas.By 2021, all 13 cities significantly improved concentration levels and contribution rates of IHS-related PM 2.5 (Figure 10b).Concentrations ranged from 7.39 µg/m 3 to 23.83 µg/m 3 , with contributions between 22% and 46%.Between 2012 and 2021, the decrease in IHS-related PM 2.5 concentrations ranged from 36.09% to 57.23%, with Xingtai experiencing the most significant drop and Qinhuangdao the least (Table 2).Correspondingly, the decrease in contribution rates varied from 9.41% to 35.62%, with Tianjin showing the most substantial decline and Qinhuangdao the least.The observations from 2012, 2017, and 2021, as well as the average over the decade, were consistent, indicating a gradual improvement, although not uniform, in the impact of IHSs on air pollution.There were significant disparities both within and between cities, necessitating continued attention to fine-tuning and implementing localized industrial emission reduction strategies.A variance analysis further determined the differences in IHS-related PM 2.5 concentrations among cities from 2012 to 2021 and the extent of fluctuations in each city.Figure 11 displays the distribution of IHS-related PM 2.5 concentrations across three critical years (2012, 2017, and 2021) and over the decade.Generally, Cangzhou had the highest median concentration in all measured years, with a wide interquartile range (IQR), indicating substantial variability in IHS-related PM 2.5 concentrations.In contrast, Zhangjiakou exhibited the lowest median concentration and a narrow IQR, suggesting consistent pollution levels with fewer outliers (Figure 11d).The observations from 2012, 2017, and 2021, as well as the average over the decade, were consistent, indicating a gradual improvement, although not uniform, in the impact of IHSs on air pollution.There were significant disparities both within and between cities, necessitating continued attention to fine-tuning and implementing localized industrial emission reduction strategies.The annual IHS-related PM 2.5 concentration in this IHS radiation area varied.From 2012 to 2013, there was an uptrend with an average increase of 7.48%.Post the launch of the 'Air Pollution Prevention and Control Action Plan' in 2013, a consistent decline in concentrations was observed from 2013 to 2016, with the rate of decrease growing each year: 7.38%, 15.93%, and 34.12%, respectively.With Xilin Technology commencing in 2017, the area's industrial heat source count rose to two.Between 2016 and 2019, the annual mean IHS-related PM 2.5 concentration showed limited fluctuations.Impacted by the COVID-19 pandemic and the dual-carbon policy, a reduction of 13.79% in IHS-related PM 2.5 concentration occurred in 2020 compared to 2019, with a continued decrease of 1.86% in 2021.
As shown in Figure 12c, the daily contributions of IHS-related PM 2.5 concentrations ranged from 7.61% to 55.06% between 2012 and 2021.The trend line indicated an overall decrease in contribution rates, with the annual average contribution in 2021 being 29.21%, a decrease of 8.39% from 37.60% in 2012.Notably, during a partial shutdown of Jinyu Cement Factory from November 2016 to March 2017, the contribution rate of IHS-related PM 2.5 fell to merely 20.95%, a substantial decline of 14.94% compared to the same period in the previous year.In sum, this observed temporal decrease and consistent industrial contributions highlight the importance of continuous monitoring of industrial emissions.

Interrelations of Factors Affecting IHS-Related PM 2.5 Concentrations
Table 3 elucidates the intricate dynamics between IHS-related PM 2.5 concentrations and various indicators signaling industrial activity, energy usage, and environmental emissions.Predominantly, metrics associated with operational IHSs, industrial indices, and environmental emissions showcased a strong positive correlation with IHS-related PM 2.5 concentrations; however, specific raw material production correlations were less pronounced.Specifically, (1) The quantity of operational IHSs across the three provinces exhibited a potent positive correlation with IHS-related PM 2.5 concentrations.Considering the direct emission of particulates from these IHSs during production processes, the number and active status of these IHSs have been identified as significant factors influencing ambient PM 2.5 concentrations.(2) Metrics like energy consumption level, industry scale, and secondary sector gross domestic product (GDP) were highly correlated with IHS-related PM 2.5 concentrations, especially in Hebei (0.93, 0.88, and 0.95, respectively).This implied that the existence of IHSs and their operational magnitudes substantially impact air quality.Significantly, regional differences were evident in the correlations between IHS-related PM 2.5 concentrations and industrial indices.Industry scale in Tianjin and secondary sector GDP in Beijing with IHS-related PM 2.5 concentrations were relatively modest.This pattern indicated that higher economic output was not necessarily synonymous with increased pollution levels, potentially signifying a shift in the economic composition of these regions towards cleaner technologies and services.(3) The correlations between raw material production and IHS-related PM 2.5 concentrations exhibited notable variations.These differences suggested that distinct production processes, varying degrees of technological adoption, and the effectiveness of pollution control measures significantly impact PM 2.5 emissions.Steel and cement production in Beijing and Tianjin moderately correlated with IHS-related PM 2.5 concentrations, indicating contributions to particulate levels, yet subject to policy-induced production moderation and response adjustments.The exceedingly high correlation with raw coal production (0.98) in Hebei signaled the region's historical dependence on coal and high-polluting industrial processes.Also, it highlighted the province's pace of coal industry reform, directly affecting PM 2.5 levels.Conversely, the negative correlation with steel production (−0.86) reflected a transition to modernized steel production technologies, ensuring environmental cleanliness while boosting production efficiency.
(4) SO 2 and NO x emissions across all three provinces were strongly correlated with IHS-related PM 2.5 concentrations, notably at 0.94 for SO 2 and 0.93 for NO x in Hebei.These correlations highlight the need for improved combustion efficiency and better desulfurization and denitrification processes.Additionally, transitioning to cleaner energy sources is essential for substantially reducing PM 2.5 levels.The high correlation between industrial wastewater discharges and carbon emissions reflects broader environmental management practices within industries.Adequate wastewater treatment and low-carbon energy sources reduce IHS-related PM 2.5 concentrations.In a word, the correlation between economic activities, industrial practices, and environmental emissions profoundly influenced IHS-related PM 2.5 concentration levels.The above correlations between IHS-related PM 2.5 concentrations and industrial activities highlighted the necessity for stricter emission controls and cleaner production processes.Furthermore, the sustained strong correlations also illustrated the environmental latency in response to industrial practices and policy changes, underscoring the necessity for continuous policy assessment and refinement.

Comparison with Other Previous Studies
The previous literature has predominantly focused on industry-specific source apportionment of PM 2.5 .However, these studies have varied in their methods of industrial categorization and analytical models.IHSs typically include cement and steel plants, oil refineries and exploration fields, chemical processing plants, and power generation facilities [48].Wang, Li, and Wang utilized different numerical simulation models-RegAEM-APSA, CAMx, and CMAQ-for the source apportionment of PM 2.5 in the BTH region during January and February 2013 [62][63][64].Their results in industrial source decomposition were comparable to the IHS-related PM 2.5 concentrations in this study (Table 4).The contribution rates of IHS-related PM 2.5 concentrations in our study were observed aligning with the industrial proportions reported in these previous studies.However, our study reported higher proportions in certain regions, possibly due to variations in classification criteria and regional delineations.IHS-related concentrations in our study, including contributions from some power plants and combustion sources, did not entirely match the industrial categories in atmospheric simulation methods.Furthermore, this study concentrated on IHS radiation areas, inherently exhibiting higher average PM 2.5 levels than the broader regional analysis in source apportionment methods.Liu used an emissions inventory and extensions-particulate matter source apportionment technology (PSAT) for decomposing PM 2.5 levels in Beijing in July 2012 [27].Their categorization of steel mills, cement plants, and power plants can be compared with the IHSs in this study.The findings showed that the contribution sources in this study were slightly lower than those in Liu (Table 4).This could be due to the classification of oil and gas development and coal chemical-related IHS contributions under other industrial sectors in Liu's analysis.Compared with earlier studies, the contribution rates of IHS-related PM 2.5 concentrations in our research on an urban scale align well with source apportionment methodologies for the same period.However, the literature has yet to offer a similar analysis for specific IHS areas, thus limiting direct comparisons.This study's advantages over previous research include: (1) Long Time Series with High Temporal Resolution: Unlike past studies that have mainly analyzed PM 2.5 compositions over shorter monthly periods, our study provided a long-term daily time series of IHS-related PM 2.5 concentrations.This supported the identification of daily anomalies and long-term trends in industrial production.(2) Focus on individual IHS targets: In advancing past the methodologies of the previous research that have predominantly utilized regional emission inventories, this study uniquely pinpointed the exact locations of each IHS.This precision facilitated a more granular analysis of the specific impact exerted by individual IHS units on PM 2.5 concentrations.Constructing a 1 km x 1 km resolution radiation area around each IHS accurately extracted IHS-related PM 2.5 concentrations for the impacted areas, providing a scientific basis for assessing the specific impact of factories on air quality.(3) Model Portability and Computational Efficiency: The model employed in this study not only demonstrated computational efficiency but also offered adaptability across different industrial and geographical contexts.This versatility enhanced the model's applicability in varied environmental research scenarios.

Significance and Uncertainties of the Study
This study extracted IHS-related PM 2.5 concentrations in the BTH region from 2012 to 2021 and analyzed their spatial and temporal trends and influencing factors.These findings are crucial for comprehending the dynamics of regional air pollution and the impact of industrial emissions on the environment.This study effectively separated PM 2.5 concentrations influenced by IHSs through the exclusion of meteorological and other non-industrial factors, offering a relatively precise delineation of industrial air pollution.Notably, this research pioneered a model tailored for extracting IHS-related PM 2.5 concentrations, enabling monitoring of the impact of individual IHSs on air quality.These innovative approaches have enriched the methods and theories for assessing the effectiveness of industrial air pollution control policies, providing valuable data and analytical tools for atmospheric policy evaluation.Furthermore, the outcomes offer vital scientific insights for developing and optimizing air quality management strategies, enhancing regional environmental quality, and fostering sustainable industrial development.
Nevertheless, this research is characterized by inherent uncertainties.Limitations in data acquisition and processing could introduce errors in the models and analytical methods employed, such as in the quality of source data and the definition of IHS radiation areas.Potential biases in the CHAP dataset, such as varying calibration standards and satellite data discrepancies, could impact PM 2.5 concentration representation.Similarly, the ERA5 dataset may have uncertainties influenced by urban heat islands and specific meteorological conditions in the BTH region, affecting temperature and precipitation accuracy.Moreover, while this study endeavored to isolate non-IHS-related PM 2.5 concentrations thoroughly, residues of influences from biomass combustion and vehicular emissions were only partially extricated.Additionally, the methodology cannot disaggregate specific PM 2.5 concentrations for each heat source in areas with multiple IHSs.
Future research should focus on refining methodologies to assess the impact of IHSs on air quality more accurately.A key area for further study involves effectively separating non-IHS-related PM concentrations, particularly those contributions from biomass burning and traffic emissions.This would likely entail integrating more diverse data types, such as traffic flow statistics and spatial distribution data of biomass burning activities, and considering using more complex statistical models and machine learning algorithms to enhance analytical precision.

Conclusions
This research developed a new three-stage model using multi-source long-term sequential data to estimate IHS-related PM 2.5 concentrations, focusing on the BTH region from 2012 to 2021.The new three-stage model based on region-growing, STL decomposition, multiple regression methods, and U-ConvLSTM algorithms can identify the IHS radiation areas and eliminate the influence of meteorological and anthropogenic factors on IHS-related PM 2.5 concentrations.Then, the daily IHS-related PM 2.5 concentrations in the BTH region were obtained.The key findings are as follows: (1) Within the study area, the average PM 2.5 concentrations in IHS radiation areas were significantly higher than in background areas, with approximately 33.16% of PM 2.5 concentration attributable to IHS activities.Furthermore, a year-over-year decline in the contribution of IHS-related PM 2.5 was observed, indicating the effectiveness of industrial reform measures.(2) The annual mean IHS-related PM 2.5 concentration in the BTH region exhibited a general downward trend with a 5.78% average annual reduction.Seasonal analysis revealed a pronounced "low in spring-summer, high in autumn-winter" pattern, with the highest concentrations in winter.Spatial distribution analysis showed that IHSrelated PM 2.5 concentrations in the southern, industrially dense areas were significantly higher than in the north, and the 13 cities within the region displayed varied temporal and spatial trends in IHS-related PM 2.5 concentrations.These findings underscore the importance of industrial activities and regional environmental policies in air pollution control.(3) In the specific industrial area of She County, Handan, two IHSs contributed an average of 19.20 µg/m 3 to the IHS-related PM 2.5 concentration.From 2012 to 2021, these concentrations fluctuated dynamically, peaking in 2013 and notably decreasing during partial shutdowns of IHS operations.This highlights the significant impact of the operational status of IHSs on local air quality.
The results not only offer a novel perspective on regional air pollution dynamics and the environmental impact of industrial emissions on a regional scale but also demonstrate the value of the developed IHS-related PM 2.5 concentration extraction model.This model enables monitoring of air quality impacts from individual IHS perspectives, providing valuable scientific support for formulating and optimizing air quality management strategies.Future studies should focus on refining the separation of IHS-related PM 2.5 from other sources like biomass burning and traffic emissions, potentially employing advanced statistical models and diverse data integration for enhanced accuracy.

Figure 2 .
Figure 2. Technology route of IHS-related PM2.5 concentration estimates and analyses in the BTH region using a three-stage model.

Figure 2 .
Figure 2. Technology route of IHS-related PM 2.5 concentration estimates and analyses in the BTH region using a three-stage model.

Figure 3 .
Figure 3. IHS radiation areas in the BTH region.(a) IHS radiation areas.(b) Location of IHSs.(c) Location of IHS radiation areas.

Figure 3 .
Figure 3. IHS radiation areas in the BTH region.(a) IHS radiation areas.(b) Location of IHSs.(c) Location of IHS radiation areas.

Figure 4 .
Figure 4. Process of meteorological and anthropogenic influence removal for estimating IHS impact on PM2.5 concentrations (illustrated using a single pix value).

Figure 4 .
Figure 4. Process of meteorological and anthropogenic influence removal for estimating IHS impact on PM 2.5 concentrations (illustrated using a single pix value).

Figure 5 .
Figure 5. Structure of the U-convlstm model.2.3.5.Statistical Analysis of IHS-Related PM2.5 Concentrations This research utilized various methods, such as the Mann-Whitney U test, correlation analysis, the Mann-Kendall test, and spatial autocorrelation analysis, to examine the

3. 1 . 1 .
Figure 6 illustrates the spatial and temporal distribution of PM2.5 concentrations throughout the BTH region.From 2012 to 2021, the area displayed a spatial pattern of higher concentrations in the southern parts and lower in the northern regions.The mean concentration per pixel varied between 26.89 μg/m 3 and 93.28 μg/m 3 , with an overall average of 54.96 μg/m 3 , as depicted in Figure 6a.Spatial analysis revealed that the average PM2.5 concentration in IHS radiation areas was 84.90 μg/m 3 , reaching a peak of 281.74 μg/m 3 and a standard deviation of 39.84 μg/m 3 .Conversely, the background areas showed an average PM2.5 concentration of 21.97 μg/m 3 , characterized by lower variability (standard deviation of 18.41 μg/m 3 ) and a peak concentration of 145.6 1μg/m 3 , with less frequent occurrences of high concentration levels compared to the IHS radiation areas, as shown in Figure 6b.Statistical analysis using the Mann-Whitney U test indicated a significant difference (p < 0.01) in the average PM2.5 concentrations between IHS radiation and background areas, indicating a discernible impact of IHS operations on PM2.5 concentrations.

Figure 6 .Figure 7
Figure 6.The spatial and temporal distribution of PM2.5 concentrations in the BTH region from 2012 to 2021.(a) Spatial distribution of average PM2.5 concentrations between 2012 and 2021.(b) Temporal distribution during 2012 and 2021 on background areas and IHS radiation areas.3.1.2.Analysis of IHS-Related PM2.5 Concentrations and Influencing Factors in the BTH Region Figure 7 displays the spatial and temporal distribution of IHS-related PM2.5 concentrations and background concentrations, along with the factors influencing and the extent of the contributions of IHS-related PM2.5.The background concentration, influenced by meteorological and other anthropogenic factors, varied between −3.92 μg/m 3 and 85.36 μg/m 3 with an average of 55.63 μg/m 3 in the BTH area from 2012 to 2021, as shown in Figure 7a.Excluding the background concentration, the IHS-related PM2.5 concentration in radiation areas averaged 29.27 μg/m 3 (Figure7b).Correlational analysis, considering the annual operational status of IHS, revealed a correlation coefficient of 0.94, highlighting a significant link between IHS operations and PM2.5 levels (Figure7c).When combining background and IHS contributions, approximately 33.16% of the total PM2.5 concentration in these areas was attributed to IHSs (Figure7d).Further, a Mann-Kendall trend analysis

Figure 6 .Figure 7
Figure 6.The spatial and temporal distribution of PM 2.5 concentrations in the BTH region from 2012 to 2021.(a) Spatial distribution of average PM 2.5 concentrations between 2012 and 2021.(b) Temporal distribution during 2012 and 2021 on background areas and IHS radiation areas.3.1.2.Analysis of IHS-Related PM 2.5 Concentrations and Influencing Factors in the BTH Region Figure 7 displays the spatial and temporal distribution of IHS-related PM 2.5 concentrations and background concentrations, along with the factors influencing and the extent of the contributions of IHS-related PM 2.5 .The background concentration, influenced by meteorological and other anthropogenic factors, varied between −3.92 µg/m 3 and 85.36 µg/m 3 with an average of 55.63 µg/m 3 in the BTH area from 2012 to 2021, as shown in Figure 7a.Excluding the background concentration, the IHS-related PM 2.5 concentration in radiation areas averaged 29.27 µg/m 3 (Figure7b).Correlational analysis, considering the annual operational status of IHS, revealed a correlation coefficient of 0.94, highlighting a significant link between IHS operations and PM 2.5 levels (Figure7c).When combining background and IHS contributions, approximately 33.16% of the total PM 2.5 concentration in these areas was attributed to IHSs (Figure7d).Further, a Mann-Kendall trend analysis

Figure 7 .
Figure 7.The spatial and temporal distribution of IHS-related PM2.5 concentrations and background concentrations in the BTH region from 2012 to 2021.(a) The average PM2.5 background concentration in the BTH region.(b) The temporal distribution of IHS-related PM2.5 concentrations and background concentrations.(c) Correlation of the number of IHSs and IHS-related PM2.5 concentrations on different years.(d) PM2.5 contribution related to IHS.

Figure 7 .
Figure 7.The spatial and temporal distribution of IHS-related PM 2.5 concentrations and background concentrations in the BTH region from 2012 to 2021.(a) The average PM 2.5 background concentration in the BTH region.(b) The temporal distribution of IHS-related PM 2.5 concentrations and background concentrations.(c) Correlation of the number of IHSs and IHS-related PM 2.5 concentrations on different years.(d) PM 2.5 contribution related to IHS.

3. 2 .
Analysis of the Spatial and Temporal Variations on IHS-Related PM 2.5 Concentrations in the BTH Region 3.2.1.Temporal Trends and Seasonal Fluctuations on IHS-Related PM 2.5 Concentrations

Figure 8 .
Figure 8. Temporal variations in IHS-related PM 2.5 concentrations in the BTH region during 2012 and 2021.(a) Temporal trends of IHS-related PM 2.5 concentrations from 2012 to 2021.(b) Rates of change in IHS-related PM 2.5 concentrations over two adjacent years.(c) Seasonal interannual variability of IHS-related PM 2.5 concentrations.(d) Seasonal trends in IHS-related PM 2.5 concentrations.Regarding yearly variations, the peak average concentration occurred in 2013 at 42.42 µg/m 3 , while the lowest was recorded in 2021 at 18.14 µg/m 3 .The average concentration represented a 52.05% reduction from 2012 to 2021.Notably, there was a 7.2% increase in IHS-related PM 2.5 concentrations from 2012 to 2013.However, from 2014 to 2021, a continual decrease was observed, particularly in 2014-2015 and 2017-2018, with reductions of 16.8% and 19.2%, respectively, as shown in Figure 8b.These declines correlated with enforcing the Air Pollution Prevention and Control Action Plan and subsequent supply-side reforms.Post-2018, a continued decrease in IHS-related PM 2.5 concentrations was observed, with an increasing annual reduction rate, culminating in an 11.6% decrease in 2021 compared to 2020, following the dual-carbon policy implementation.From 2012 to 2021, seasonal variations in IHS-related PM 2.5 concentrations in the BTH region showed a distinct pattern (Figure8c).The trend displayed a 'low in springsummer and high in autumn-winter' pattern, with winter concentrations being the highest, followed by autumn, spring, and summer.Spring concentrations averaged 25.50 µg/m 3 , summer at 18.46 µg/m 3 , autumn at 28.77 µg/m 3 , and winter at 44.64 µg/m 3 .Compared to 2012, each season in 2021 showed improvements, with varying degrees of reduction (Figure8d).Summer saw the most significant decrease at 64.20%, an annual average

Figure 9 .
Figure 9. Spatial distribution of average IHS-related PM 2.5 concentrations between 2012 and 2021.(a) The pixel-level average IHS-related PM 2.5 concentrations.(b) The city-level average IHS-related PM 2.5 concentrations.

Atmosphere 2024 ,
15, x FOR PEER REVIEW 18 of 27 respectively.Cities like Baoding, Beijing, Hengshui, Langfang, Tianjin, Xingtai, and others experienced significant reductions in IQR fluctuations, with Hengshui showing the most notable decrease, dropping from 25.09 μg/m 3 in 2017 to 15.60 μg/m 3 in 2021.Chengde, Qinhuangdao, and Tangshan saw slight increases in IQR fluctuations, with Tangshan having the largest increase of 31.32%.Other areas maintained relatively stable IQR fluctuations.

Figure 10 .
Figure 10.The contribution of IHS-related PM2.5 concentrations in the BTH region.(a) The annual average contributions of IHS-related PM2.5 in cities of the BTH region from 2012 to 2021.(b) Annual mean IHS-related PM2.5 concentration and contribution rate changes in cities of the BTH region from 2012 to 2021.

Figure 10 .
Figure 10.The contribution of IHS-related PM 2.5 concentrations in the BTH region.(a) The annual average contributions of IHS-related PM 2.5 in cities of the BTH region from 2012 to 2021.(b) Annual mean IHS-related PM 2.5 concentration and contribution rate changes in cities of the BTH region from 2012 to 2021.

3. 3 .
Analysis of IHS-Related PM 2.5 Concentrations Focusing on Typical IHS In the southwest of She County in Handan (113 • 35 ′ 00 ′′ E-113 • 42 ′ 00 ′′ E, 36 • 28 ′ 00 ′′ N-36 • 35 ′ 00 ′′ N), the IHS radiation area included two adjacent IHSs: Jinyu Cement Factory and Xilin Technology Company (Figure 12a).While Jinyu Cement Factory has been continuously operational, Xilin Technology commenced operations in 2017.From 2012 to 2021, the IHS-related PM 2.5 concentration in this industrial area fluctuated with an average concentration of 19.20 µg/m 3 .A clear peak in IHS-related PM 2.5 concentration occurred in 2013 (Figure 12b).Despite these short-term fluctuations, a declining trend was evident (S = −0.005,p < 0.05), signifying the effectiveness of long-term emission reduction initiatives.Atmosphere 2024, 15, x FOR PEER REVIEW 20 of 27 a decrease of 8.39% from 37.60% in 2012.Notably, during a partial shutdown of Jinyu Cement Factory from November 2016 to March 2017, the contribution rate of IHS-related PM2.5 fell to merely 20.95%, a substantial decline of 14.94% compared to the same period in the previous year.In sum, this observed temporal decrease and consistent industrial contributions highlight the importance of continuous monitoring of industrial emissions.

Figure 12 .
Figure 12. Analysis of IHS-related PM 2.5 concentrations focusing on southwest of She county in Handan.(a) The location of IHS radiation area.(b) Temporal trends of IHS-related PM 2.5 concentrations from 2012 to 2021.(c) PM 2.5 contribution related to IHSs.

Table 1 .
Data description used in this study.

Table 1 .
Data description used in this study.

Table 2 .
IHS-related PM 2.5 concentrations (µg/m 3 ) in cities of the BTH region from 2012 to 2021.

Table 3 .
Correlation coefficients between IHS-related PM 2.5 concentrations and associated indicators.

Table 4 .
Comparison of the contribution of this study with previous studies of IHS-related PM 2.5 concentration.