石家莊市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濃度呈下降趨勢。