石家庄市PM2.5时空特征及其对土地利用变化的响应

赵文斐，于占江，王让会,3*，高祺，彭擎，杨鹏

1. 南京信息工程大学应用气象学院，江苏 南京 210044；2. 石家庄市气象局，河北 石家庄 050081；3. 江苏省农业气象重点实验室，江苏 南京 210044

大气污染是中国城市化快速发展过程中亟需解决的问题，大气颗粒物已成为影响中国空气质量的首要污染物（Huang et al.，2012），尤其大气细颗粒物（PM2.5）污染问题更加严重。PM2.5不仅会危及人体健康，引起呼吸系统甚至心血管系统疾病（Mirabelli et al.，2016），还会影响城市宜居性（刘圆等，2017）以及区域生态系统平衡（王让会，2018）。2012年新颁布的国家空气质量标准将PM2.5设为常规监测指标，PM2.5已经成为中国大气污染研究的热点问题。目前对PM2.5的研究大多集中在来源解析（Huang et al.，2014）、时空特征（王飞等，2017）、空间分布模拟（许刚等，2016）和健康风险评价（李晓宝等，2020）等方面。

现有研究表明影响PM2.5浓度的因素有生物化学因素（如气溶胶光学厚度 AOD）（焦利民等，2015）、气象因子（如温度、风速）（刘寿东等，2018）、人类活动（Johansson et al.，2009）等。此外，作为大气下垫面的土地的类型、规模和景观格局也会直接或间接影响 PM2.5污染程度，一方面不同地类作为“源-汇”景观直接影响PM2.5浓度（Liu et al.，2018），另一方面则会改变区域水热条件从而产生间接影响（王方等，2016）。目前关于土地利用变化和 PM2.5关系的研究大概分为三方面：一是基于相关性分析法来定性研究土地利用类型和 PM2.5浓度的时空关联性，有研究单一地类与PM2.5的关系（杨昆等，2017），也有以研究区土地利用/覆被为整体与PM2.5浓度做相关性研究（唐新明等，2015）；二是以土地利用作为影响因素，构建土地利用回归（LUR）模型来探究对PM2.5浓度的影响（杨婉莹等，2019），以及实现区域PM2.5空间分布的模拟（满旺等，2019），在此基础上还有许多研究利用地理加权回归模型、主成分分析法对该模型进行优化（刘炳杰等，2018；Olvera et al.，2012）；此外，还有研究基于景观生态学原理，采用景观格局分析法研究土地利用变化对PM2.5浓度的影响，并表明景观类型的边缘密度（ED）、斑块密度（PD）等景观格局指数均与其密切相关（谢舞丹等，2017）。然而，大多研究的 PM2.5数据来源于站点监测，研究多针对短时间内的城镇区域，缺乏其他区域和长时间序列的研究，其时空连续性均存在一定的局限性。

《2018生态环境状况公报》表明，石家庄市是环境空气质量相对较差的城市之一（中华人民共和国生态环境部，2019），虽然近年来其空气质量有所好转，但大气污染问题依然严重。同时石家庄市也属于“2+26”京津冀大气污染传输通道，如何防治大气污染仍是石家庄市的当务之急。因此，该研究以石家庄市为研究区域，基于高精度遥感反演数据分析其2000—2017年PM2.5时空变化特征，并进一步探讨 PM2.5浓度对土地利用变化的响应，以期在合理控制人类活动和严格控制污染物排放的基础上，合理配置土地资源，通过优化土地利用方式促进大气污染防治，更好地打好蓝天保卫战，实现“绿水青山就是金山银山”的理念。

1 材料与方法

1.1 研究区概况

石家庄市（37°27′— 38°47′N ， 113°30′—115°20′E）位于河北省中南部（图1），地跨太行山和华北平原两大地貌单元，地势西高东低，西部为太行山中段，东部为滹沱河冲积平原。石家庄市地处环渤海湾经济区，是京津冀地区重点发展城市之一，城镇化发展十分迅速，同时还拥有丰富的能源和矿产资源，是河北省的工业大市，但伴随而来的便是该地区大气污染问题突出，常年位居全国城市空气质量排名的末位。近年来，虽然石家庄市地区“煤烟型”空气污染问题有所改善，但其空气质量状况仍然不容乐观。

1.2 数据来源及处理

1.2.1 PM2.5浓度数据

选用PM2.5浓度的年均值作为分析指标，PM2.5空间分布数据来源于大气成分分析组织（Atmospheric Composition Analysis Group，ACAG）的中国区域估算数据集（V4.CH.02），该数据集基于中国 PM2.5监测值进行了修正及地理加权回归校准，空间分辨率为0.01°，目前已被国内许多研究所引用（刘海猛等，2018；Luo et al.，2017）。

本文利用ArcGIS 10.2软件对2000—2017年全国PM2.5空间分布数据进行重投影和裁剪，提取出石家庄市各年的PM2.5数据，然后利用分区统计工具得到不同县市和不同地类区域的PM2.5年均浓度。

1.2.2 土地利用数据

利用遥感影像进行监督分类解译得到石家庄市土地利用分类数据。遥感数据来源于地理空间数据云，网址为http://www.gscloud.cn/，影像数据集为Landsat TM和Landsat OLI，选用2000年和2017年云量低于10%的遥感影像进行处理使用，具体信息见表1。

1.3 研究方法

1.3.1 热点分析

图1 研究区位置Fig. 1 Location of study area

表1 遥感数据参数Table 1 Parameters of remote sensing data

利用基于Getis-Ord Gi*的热点分析法来探究研究区内PM2.5浓度局部聚集特征。热点分析是被广泛应用于多领域的的局部自相关统计，该方法是通过局部区域的所有值来识别具有统计显著性的低值和高值的空间聚类，即冷点和热点。其计算模式如下（Ord et al.，1995）：

式中，n为空间数据单元总数；xj是数据单元j的属性值；wi,j是数据单元i和j之间的空间权重，表示i和j的空间邻接关系。对PM2.5浓度进行Getis-OrdGi*统计会得到一个z得分和p值，z得分（即Gi*）绝对值大于1.65时才具有统计意义，其值越高代表聚集程度越高。

1.3.2 广义相加模型（GAM）

广义相加模型（GAM）是由数据驱动的非参数回归模型（南洋等，2020），GAM模型不需要提前设立参数模型，可以根据各解释变量的光滑函数建立因变量与解释变量之间的非线性关系。研究建立了石家庄市2 km×2 km的渔网数据，共产生了3600个网格单元，通过R语言4.0.0版本中mgcv包的GAM函数构建非线性模型：

式中，g(μi,j)为连接函数；μi,j代表网格的PM2.5质量浓度（μg·m-3）；解释变量有耕地（cropland）、林地（forest land）、草地（grassland）和建设用地（construction land）；s为针对解释变量的光滑函数，该函数是通过p样条、三次样条等方法对散点图进行光滑处理得到的；β0为截距；ε为随机变量。

1.3.3 地理加权回归分析

地理加权回归（GWR）是量化空间异质性的有效工具，通过建立空间范围内不同位置处的局部回归方程，来分析存在空间非平稳性的问题，使得结果更符合实际情况（Propastin，2012）。GWR是对OLS模型的扩展，将要素的空间位置嵌入到回归方程中，其形式为：

式中，i代表空间单元要素；εi为误差项，(ui,vi)为要素i的地理坐标；βk(ui,vi)是要素i上的第k个回归参数，回归参数通过空间权函数的方法得到，常见的有距离反比法、距离阈值法和Gauss函数法等。将研究区划分为2 km×2 km和3 km×3 km的渔网，分别产生了3600个和1636个网格单元，以每个网格区域内2000—2017年PM2.5浓度变化量为因变量，以耕地、林地、草地和建设用地面积变化量为解释变量，构建GWR模型来分析PM2.5浓度对地类变化的响应情况。

2 结果分析与讨论

2.1 PM2.5时空变化特征

图2为石家庄市2000—2017年的PM2.5年均浓度变化图，中国《环境空气质量标准》（GB 3095—2012）（环境保护部，2012）规定PM2.5年平均浓度一级浓度限值为 15 μg·m-3，二级浓度限值为 35 μg·m-3，可以看出，研究区自2000年以来PM2.5浓度年均值始终高于二级浓度限值，属于空气质量不达标，可见研究区的大气污染形势严峻。2000—2017年石家庄市 PM2.5年均浓度总体呈升高的态势，PM2.5年均浓度由2000年的46.70 μg·m-3增加到了2017年的67.49 μg·m-3，增幅达44.5%。其中，2000—2006年呈波动增加，并于2006年达到最高值83.81 μg·m-3，这是由于该期间石家庄市工业迅速发展，能源消费量增加造成的。2007—2017年呈波动下降趋势，其中2008年PM2.5年均浓度呈大幅度下降，这是因为2007年国务院批准了《第29届奥运会北京空气质量保障措施》，各省市实施了严格的污染治理和减排措施；2013年国务院发布了《大气污染防治计划》，石家庄市制定了相关的攻坚行动计划来防治大气污染，PM2.5年均浓度自2013年开始持续下降，可见石家庄市大气污染防治工作取得了一定的成效。

图2 石家庄市PM2.5年均浓度变化图Fig. 2 Variation trend of PM2.5 in Shijiazhuang

从空间分布来看，石家庄市PM2.5空间分布呈西北部低、中部和东部高的格局，结合图3可发现研究区2000—2017年来均保持这一空间分布格局。石家庄市西部地区 PM2.5多年平均浓度最低值为45.04 μg·m-3，最高值为 65.38 μg·m-3；东部地区PM2.5多年平均浓度最低值为 77.85 μg·m-3，平均浓度最高值为 80.75 μg·m-3；石家庄市中部市区的PM2.5浓度为最高值，PM2.5多年平均浓度最高达88.66 μg·m-3。此外，其余区县均呈现城镇区域高、周边区域低的分布情况。究其原因，一方面是由于城镇区域人口密集，城市建设和交通排放等人类活动带来的扬尘和汽车尾气使得大气颗粒物浓度升高，导致市区和城镇的PM2.5浓度远高于农村和郊区；另一方面，城镇区域内各种的工业活动和能源使用也会造成大气污染物密集排放，致使其大气颗粒物浓度高于周边郊区（Zhao et al.，2009）。

利用热点分析法探究石家庄市PM2.5浓度空间聚集特征（图4），可以看出2000年和2017年石家庄市的PM2.5浓度冷热点分布情况大致保持不变，但不同等级的冷点和热点区域的面积和数量发生了一定的变化。PM2.5浓度的热点主要集中在研究区中部和东部地区，包括石家庄市中心城区辛集市的大部分区域，以及晋州市、赵县等的城镇区域；2017年热点区在原有基础上面积有所扩大，如栾城区、藁城区、高邑县等，而晋州市、辛集市南部和赵县东部的热点区相较于 2000年有一定的缩小，但热点区整体分布相较于 2000年呈南移趋势。热点区所在地城镇化水平高，地势较平坦，交通及工厂污染排放较多；加之各区县之间联系紧密，各城区PM2.5的流动和积累互相影响，进而形成大气污染热点区域。此外，研究区西部井陉矿区由于以工矿用地为主，空气污染也较重，因而也表现为热点区域。冷点区则在灵寿县、平山县、井陉县和赞皇县的西部地区，分布区域较稳定，2017年核心冷点区的空间范围则有所扩大。冷点区所在区地处太行山地，城镇分布较少，污染排放较少，同时林地、草地地类占比大，因而空气质量相对较好。

图3 石家庄市不同年度PM2.5年均浓度分布图Fig. 3 Distribution of PM2.5 in different years in Shijiazhuang

图4 石家庄市PM2.5年均浓度热点分布图Fig. 4 Hotspot distribution of PM2.5 in Shijiazhuang

2.2 土地利用动态变化特征

由图5和表2可以看出，石家庄市的土地利用类型主要为耕地、林地、草地和建设用地，四种用地面积占石家庄市总面积的97%以上，耕地类型面积占全市面积50%左右。

表2 2000年和2017年石家庄市土地利用类型面积及所占比例Table 2 Area and proportion of land use types in 2000 and 2017 in Shijiazhuang

从土地利用类型的动态变化来看，2000—2017年石家庄市耕地面积呈减少趋势，其占比由2000年的57.81%减少到2017年的48.69%；林地面积则有所增加，主要表现在石家庄市东南部地区的林地类型用地增加；建设用地面积呈增加趋势，其占比由9.46%增加到16.45%；草地、水域和未利用土地的面积则变化不大。

图5 石家庄市土地利用分类Fig. 5 Land use types in Shijiazhuang

表3 2000—2017年石家庄市土地利用类型转移矩阵Table 3 Area conversion matrix for land use change in 2000-2017 in Shijiazhuang

石家庄市 2000—2017年土地利用类型转移矩阵如表3所示，可以看出石家庄市的土地利用类型变化状况较为显著。其中，耕地的转出较为显著，转出占比为12.41%。耕地主要转出的用地类型为建设用地，占耕地总转出的63.87%，发生转移的地区多为石家庄市各市县城镇的周边区域，表明了其城镇化发展趋势；其次，耕地还转移为林地和草地类型，分别占总转出的18.27%和1.51%，尤其是石家庄市东南部地区表现明显，这与退耕还林政策有关。林地类型主要转移为草地类型，转出了林地面积的3%，1.34%和1.11%的林地分别转移为建设用地和耕地。建设用地则有10.81%转变为耕地，转为其他类型的用地则变化不大。草地、水域和未利用土地部分地区也有所增减，但由于转入转出占比相近，其面积基本没有太大的变化。

2.3 PM2.5对不同地类面积的非线性响应

表4为石家庄市各土地利用类型的PM2.5浓度最值和均值，结果表明建设用地的PM2.5浓度均值最高，2000年和2017年分别为59.31 μg·m-3和77.18 μg·m-3，远高于其他土地利用类型；2000年林地的PM2.5浓度最低，为 27.58 μg·m-3，2017 年则是未利用土地 PM2.5值最低，其次为林地和草地。综合来看，各土地利用类型的PM2.5浓度值大小特征为：林地、未利用土地＜草地＜水域＜耕地＜建设用地。

表4 各土地利用类型区域的PM2.5特征情况Table 4 Characteristics of PM2.5 in various land use types μg·m-3

基于石家庄市2 km×2 km渔网数据，构建GAM模型分析 PM2.5浓度对不同地类面积变化的非线性响应，由于石家庄市水域和未利用土地面积占比很小，故本文不予考虑。以通过单因素分析显著性检验的4个地类为解释变量，PM2.5浓度为因变量，构建GAM模型。PM2.5浓度和不同地类面积的GAM拟合结果如表5所示，2000年和2017年各解释变量均通过显著性检验，置信水平为0.001，模型方差解释率分别为 78.6%和 77.6%，调整后R2分别为0.785和0.774，表明模型拟合结果良好且四种地类面积均与 PM2.5浓度有显著的非线性关系。F值反映模型各解释变量对因变量的相对重要程度（胡成媛等，2019），可以发现2000年建设用地、草地的F值分别为105.69、95.67，是影响PM2.5浓度的主要因素，耕地的影响程度次之，F值为40.08；2017年的拟合结果与2000年类似，建设用地（F=105.85）和草地（F=85.58）为影响PM2.5浓度的主导因子。

表5 PM2.5浓度和不同地类面积的GAM拟合结果Table 5 GAM fitting results of PM2.5 and the area of different land use types

图6为2000年和2017年GAM模型中不同地类面积对PM2.5浓度的影响效应，图中实线表示解释变量对PM2.5浓度的平滑拟合；虚线间阴影部分为置信区间。横坐标为不同土地利用类型的面积，纵坐标为解释变量对 PM2.5浓度的平滑拟合值，纵坐标括号中为自由度，自由度值大于1时为非线性关系，其值越大，非线性关系越显著。

由图 6可知，PM2.5浓度随耕地面积的增加缓慢上升，且 2017年耕地面积增加到一定程度后PM2.5浓度上升趋势开始加剧，说明耕地面积和对PM2.5浓度的非线性影响表现为正效应。众多研究表明，耕地对大气颗粒物浓度有“源”和“汇”两方面的作用，在农作物生长旺盛时，可以吸附大气颗粒物；而在播种和收割季节，化肥施放、土地裸露、秸秆燃烧等均会导致耕地区域 PM2.5浓度升高（Li et al.，2014）。在年均尺度下，耕地的双向季节效应有所削弱，整体表现为耕地面积对PM2.5浓度有正向促进作用，这是由于石家庄市地区耕地主要分布在中东部平原地区，这一区域工业发达，污染物排放量大，加之农田连片分布导致的化肥大量使用，使得耕地的“源”作用更明显。2000年和2017年均表现为 PM2.5浓度随林地面积的增加而波动下降，且林地面积增加到一定程度时（约3.5 km2），PM2.5浓度下降加剧，这说明林地规模达到一定规模时，对PM2.5的滞留效应会更加明显。许多研究也表明，林地密度和规模越大，其滞尘能力越强，大气颗粒物浓度越小（许珊等，2015）。草地对PM2.5浓度的影响同样表现为反向抑制，随着草地面积的增加，PM2.5浓度波动下降，但相比于林地，草地面积对PM2.5浓度的影响效应较小。这是由于草地覆盖能够有效缩小地面尘源，从而降低PM2.5浓度，但由于其植被类型较矮小，无法吸附空中的大气颗粒物，因而其滞尘能力略低于林地（彭威等，2014）。在不同年份里，建设用地的面积对PM2.5浓度的影响均表现为正向促进，且当建设用地面积大于一定值（2 km2）时，PM2.5浓度增长趋势也由缓变陡。一方面，建设用地面积会增加人类活动以及工业、交通污染排放，这些均是PM2.5的主要来源；另一方面，相较于其他地类，建设用地受人类活动扰动更加频繁，汽车尾气等颗粒物更容易扩散到空气中，从而导致PM2.5污染加重。有研究表明，建设用地最大斑块面积和面积比等指标与颗粒物浓度呈正相关，且最大斑块面积越大，对大气颗粒物浓度影响越显著，这与GAM模型得出的结果相一致（崔岩岩，2013）。

图6 PM2.5浓度对不同地类面积的非线性响应Fig. 6 Nonlinear response of PM2.5 to area of different land use types

2.4 PM2.5对土地利用变化的响应

基于研究区 2 km和 3 km的渔网数据，以2000—2017年PM2.5浓度变化作为因变量，以这一期间石家庄市主要土地利用类型的变化情况作为解释变量，在通过OLS建模预测后，利用地理加权回归（GWR）模型进行构建拟合模型，从而分析PM2.5空间变化对土地利用变化的响应。

图7为地理加权回归模型的拟合结果，其中图7a、b是以2 km渔网化的石家庄市数据进行地理加权回归拟合的结果；图7c、d则是以3 km渔网进行地理加权回归拟合的结果。2 km渔网的模型结果的R2为0.87，调整后R2为0.86，模型拟合结果较好；图7a为模型拟合的标准化残差（StdResid），其值绝对值小于 2的结果有 94.42%，表明模型拟合较稳定；图7b为地理加权回归的LocalR2分布情况，结合石家庄市土地利用类型变化情况来看，LocalR2低值区为土地利用类型未发生显著变化的地区，而LocalR2高值区则是发生了明显土地利用类型变化的地区，如石家庄市中部和东部城镇扩展区和西南部果林种植区。3 km渔网的模型结果的R2为0.82，调整后R2为0.80，模型拟合结果同样不错；图7c的标准化残差（StdResid）结果有94.5%处于[-2, 2]范围内，模型稳定性较好；图7d中LocalR2低值区和高值区分布与图 7b基本一致，这表明石家庄市PM2.5变化对土地利用类型的改变有着明显的响应。杨伟等（2020）基于GWR模型探究华北地区大气细颗粒物对土地利用/覆被变化的响应，其地理加权回归计算结果同样表明土地利用/覆被变化的区域，其LocalR2为高值区。

由 2000—2017年石家庄市 PM2.5整体趋势可知，石家庄市 PM2.5浓度的升高主要是由于人为活动导致的排放量增加而造成的，为去除该影响，先分析未发生土地利用变化区域的 PM2.5的变化情况，并以此为基础探究土地利用变化情况下对PM2.5的影响。如表6所示，不同土地利用类型在未发生地类转化时，其PM2.5增幅均在21%左右。

表6 未发生土地利用变化区域的PM2.5浓度变化情况Table 6 Change of PM2.5 concentration in the area without land use change μg·m-3

图7 地理加权回归模型拟合结果Fig. 7 Results of GWR model fitting

表7 发生土地利用变化区域的PM2.5浓度变化情况Table 7 Change of PM2.5 concentration in the land use change area μg·m-3

以未发生土地利用类型变化的PM2.5浓度变化幅度（表 6）作为参考标准，结合发生地类变化的PM2.5浓度统计值，分析发生土地利用类型转化区域的PM2.5变化趋势（表7）。结果表明，当用地类型由林地转化为耕地、草地和建设用地时，PM2.5浓度呈升高趋势，其中，PM2.5浓度对林地转建设用地的响应最大，升高了5.12 μg·m-3，林地转化为耕地时浓度升高了2.84 μg·m-3，林地转化为草地时浓度变化量相对较小，升高了1.21 μg·m-3。草地与林地相似，转变为耕地和建设用地时同样会造成PM2.5浓度不同程度的升高，浓度增量分别为2.86 μg·m-3和1.20 μg·m-3；草地转变为林地时，PM2.5浓度下降了0.71 μg·m-3，其变化幅度并不明显。当地类由建设用地转变为林地、草地时，PM2.5浓度则有明显下降的趋势，转变为林地时下降最多，为4.34 μg·m-3，转变为草地下降了2.73 μg·m-3。耕地转变为建设用地时 PM2.5浓度呈增加趋势，增加量为 3.54 μg·m-3；而当耕地转变为林地和草地时，PM2.5浓度呈下降趋势，且对转变为林地的响应最为明显，下降了6.54 μg·m-3，而耕地转变为草地时浓度下降了 1.32 μg·m-3。综上，耕地、林地和草地均属于不同类型的生态用地（龙花楼等，2015），当生态用地向建设用地转化时，会导致 PM2.5浓度显著上升，而建设用地转化为绿地时，PM2.5污染显著下降，可见在城市建设时增加绿地规划可在一定程度上缓解大气污染；而耕地转变为绿地引起的 PM2.5浓度下降也表明石家庄市落实的退耕还林政策对于其大气污染防治有着不错的成效。

3 结论

本文利用石家庄市2000—2017年的PM2.5年均浓度数据和遥感影像数据，基于GIS和RS技术分析了石家庄市 PM2.5年均浓度的时空变化特征和土地利用动态变化特征，在此基础上结合GAM模型、GWR模型等方法探究了 PM2.5浓度对土地利用变化的响应，得出了如下结论：

（1）石家庄市PM2.5浓度空间分布差异显著，呈西北低、中东部高的空间格局，空间聚集特征表现为热点区增多并且有向南扩张的趋势，冷点区范围基本保持稳定。时间上来看，2000—2017年PM2.5年均浓度总体呈升高态势，2006年达到污染峰值，PM2.5浓度值为 83.81 μg·m-3，之后 PM2.5浓度波动下降，大气污染防治工作取得一定成效。

（2）石家庄市土地利用类型以耕地、林地、草地和建设用地为主，2000—2017年土地利用变化情况表现为建设用地和林地的增加以及耕地的大量减少，草地转入转出面积相近，因而总面积变化不大，水域和未利用土地面积基本没变。

（3）GAM模型多因素拟合结果表明，PM2.5浓度对各地类（耕地、林地、草地、建设用地）的面积存在显著的非线性响应，具体表现为耕地、建设用地面积增加对PM2.5浓度升高呈正效应，而林地、草地面积增加对PM2.5浓度升高呈负效应。

（4）GWR模型拟合结果表明，LocalR2较高区域为土地利用类型变化明显的区域，说明PM2.5浓度变化对地类变化有着显著响应。具体表现为林地、草地和耕地等生态用地转变为建设用地时PM2.5浓度升高，建设用地和耕地向林地、草地转换时，PM2.5浓度呈下降趋势。