粵港澳大灣區秋冬季節PM2.5濃度時空模擬分析

李天翔, 蔣 超, 龔建周*

(1.廣州蘢騰園林景觀設計有限公司, 廣東 廣州 510520; 2.廣州大學 地理科學與遙感學院, 廣東 廣州 510006)

伴隨國家層面賦予了大灣區新的發展內涵，大灣區正成為我國現代區域發展新的增長極，但也需關注增長給大灣區生態環境安全帶來的壓力[1-2]。廣東省政府聯合港、澳特區政府于2014年建立了我國首個區域空氣監測網絡——粵港澳珠江三角洲區域空氣監測網絡，網絡監測數據顯示：①2015年PM2.5年均濃度值超過了國家《環境空氣質量標準》(GB3095-2012)規定的二級質量濃度限值35 μg·m-3;②2016-2018年連續3年的年均濃度均未超過二級質量濃度限值[3-5];③2020年各子站相關達標率在98.6%及以上，大灣區空氣質量正逐年改善(資料來源于廣東省生態環境廳,網址為http://gdee.gd.gov.cn/kqjc/content/post_2725979.html)。

PM2.5能直接進入呼吸道深部和肺泡，危害人體健康。同時，由于PM2.5對可見光具有消光作用，使能見度降低，導致城市環境健康水平下降。隨著城市高質量發展和空氣污染治理進入精細化階段，通過協調大灣區的規劃與管理降低超標期PM2.5濃度將是下一階段空氣質量全面達標的重點內容，因此，需要加強針對空氣污染超標期、典型區的相關研究。在濃度監測方面，因氣溶膠與PM2.5有相近的物理機制，衛星遙感亦可大空間、多時相監測氣溶膠[6-8]，因此，衛星遙感監測氣溶膠間接估測PM2.5濃度被廣泛應用在大氣環境研究中[9]，相較于地表覆被、人類活動等要素[10-11]，溫度[12]、濕度[13-14]和大氣壓[15]等氣象條件對單一城市尺度或城市群的PM2.5濃度影響更加顯著[16-17]。

2015年PM2.5年均濃度值超過了《環境空氣質量標準》(GB3095-2012)規定的二級質量濃度限值，若對全年的逐日濃度進行數據分析，可獲取2015年逐日濃度超標發生的日期，即可作為本次研究的超標期。在此基礎上，本文將利用2015年超標期的氣溶膠光學厚度與氣象要素數據構建多元回歸模型，采用逐步回歸模擬2015年日均濃度超標時的PM2.5濃度空間特征。此外，濃度序列數據可識別濃度演變趨勢，對精準防控具有重要意義，序列數據已揭示研究區秋冬月份為PM2.5濃度管控的重點時段[4]。在各類序列數據模擬方法中，ARIMA模型(AutoRegressive Integrated Moving Average model)能根據時序數據變量的自身特征構建模型以實現短期模擬，具有精度高、易操作的優勢[18]。彭斯俊等[19]、牟敬鋒等[20]運用ARIMA模型分別對PM2.5濃度和空氣質量指數(AQI)進行分時段的短期預測，發現其預測模擬的結果均優于灰色GM(1,1)模型和全年時間序列模型。對大灣區PM2.5濃度季節變化規律的研究表明，冬季和秋季是濃度最高的重點時段[21-22]，本文將利用2019年以來秋冬季節的濃度序列數據，通過構建ARIMA模型模擬2021年11～12月的濃度狀況，進而分析2019-2021年重點時段的濃度值演變、同期數值增減情況。從綜上兩方面識別超標期下大灣區PM2.5濃度的時空演變，為粵港澳聯防、聯控區域大氣污染研究提供參考資料。

1 研究區域

粵港澳大灣區(111°21′～114°53′E，21°27′～24°24′N)大部分位于北回歸線以南，屬熱帶、亞熱帶氣候，年均溫在21 ℃～23 ℃，年平均降水量在1 300～2 500 mm，全年溫暖濕潤[23]。灣區城市群由廣東省9市及2個特別行政區組成，總面積5.6萬km2(圖1)。截止2020年底，大灣區共有人口約7 000萬，區域經濟總量達1.4萬億美元，是我國最具經濟活力、創新氛圍最濃、人口素質較高的地區之一。

圖1 研究區區位示意圖

2 數據及方法

2.1 數據及預處理

本文使用的數據包括PM2.5逐時濃度數據集、MODIS AOD和氣象數據。

(1)PM2.5逐時濃度數據集來源于中國環境監測總站空氣質量監測發布平臺，該數據涵蓋了粵港澳區域空氣監測網絡56個子站的常規空氣質量監測指標。本研究所用為2015年以來56個子站的PM2.5逐時濃度數據，并于建模前采用濕度訂正的方法進行PM2.5濕度訂正以提高模擬精度。濕度訂正的公式如下所示[24]：

(1)

(2)MODIS AOD是空氣污染研究中最為常用的AOD產品數據之一，在全球、區域大氣污染的宏觀分布方面具有很大潛力[14]。本研究所使用的MODIS數據L1B產品可從NASA官方遙感影像數據下載(https://ladsweb.modaps.eosdis.nasa.gov/search/)，并利用AERONET(AErosol RObotic NETwork)提供的監測數據作為真實值，檢驗AOD反演數據的精度[24]。

(3)PM2.5濃度相關的氣象要素[7]，包括各站點氣壓(P)、氣溫(T)、相對濕度(U)、風向風速(W)，以上要素通過克里金插值以生產與反演AOD數據空間分辨率一致的柵格數據集。氣象數據由中國科學院資源環境科學數據中心提供(http://www.resdc.cn/Default.aspx)。

2.2 研究方法

2.2.1 AOD數據反演

我國華南區域四季植被覆蓋度較高，適用“暗像元算法”反演氣溶膠光學厚度數據[17]，本文利用MODIS L1B數據反演大灣區氣溶膠光學厚度，包括輻射校正、幾何糾正、重采樣、角度數據的幾何校正、去云處理和氣溶膠光學厚度反演等流程，采用絕對百分比對反演AOD數據進行誤差檢驗[25]。

2.2.2 多元逐步回歸

逐步回歸的思路：將自變量逐個引入方程，每添加一個變量后，對已選入的原變量進行顯著性檢驗，若已引入變量因新引入變量而對被解釋變量的影響不再顯著時，將剔除已引入變量。在方程中引入新變量或剔除不顯著變量被稱為逐步回歸的一個步驟。對逐步回歸的每一個步驟反復進行F檢驗，確保顯著變量被引入，不顯著變量被刪除。

相較于其他多元線性回歸分析，逐步回歸分析優勢表現為在多元線性回歸的基礎上，模型逐步引入、剔除變量，從而建立最優的回歸模型，避免模型中變量的多重共線性，變量量綱不同也會影響分析結果，本研究使用歸一法對數據進行標準化處理，消除數據量綱[26]。

2.2.3 基于ARIMA(AutoRegressive Integrated Moving Average)時序模擬模型

ARIMA模型是基于Box-Jenkins模型提出的一種針對時間序列的模擬模型，由自回歸模型(AR模型)和移動平均模型(MA模型)構成，即ARIMA(p,d,q)模型[27]。ARIMA模型的建模步驟如下：①序列平穩性檢驗；②根據自相關和偏相關圖對模型進行定階，確定p,q值；③在保證SIC與AC最小的原則下，選擇最佳的p,q值。在此基礎上對所建模型進行白噪音檢驗，利用所建模型檢驗其對原序列的擬合效果，最后利用模型ARIMA(p,d,q)進行未來短期時間尺度模擬，并用雷達圖分析2019-2021年濃度高值月份下PM2.5濃度時序特征。

3 結果與分析

3.1 超標期確定與AOD反演

表1顯示了2015年內大灣區PM2.5日濃度超標詳情，濃度最大值出現在2月11日，最大濃度為107.48 μg·m-3，超標出現月份為1月、2月與12月。由于每月出現超標的時間相對集中，將超標時間進行合并，然后下載超標時間內的遙感影像(T1、T2、T3)，利用絕對百分比誤差檢驗本研究反演AOD數據，誤差通過檢驗，可開展下一步研究。

表1 研究區PM2.5污染超標詳情及AOD反演誤差檢驗

3.2 2015年PM2.5濃度超標時的空間模擬

利用ArcGIS的數據提取工具在各變量圖層分別提取建模所需變量數據，構建多元逐步回歸模型，模型匯總見表2。逐步回歸分析逐次引入變量(大氣壓)X2，(AOD)X1，(風速)X5，(相對濕度)X3，隨著模型中變量數量的增加，調整后R2逐漸增加，模型中變量X2，X1，X5，X3能共同解釋因變量變異的68.2%。

表2 模型匯總統計表

依據《環境空氣質量標準》(GB3095-2012)，將濃度低于35 μg·m-3(優)定義為第一級；鑒于污染等級為“良”的范圍較大，為進一步細化濃度值結構特征，將濃度介于(35～75]細分為4級，步長為10，各濃度區間依次為(35～45]、(45～55]、(55～65]和(65～75]，濃度等級依次記為：良(1)、良(2)、良(3)和良(4)；將濃度高于75～115 μg·m-3(輕度污染)定義為第6級，以上共包括6個濃度等級。2015年超標期內PM2.5的濃度狀況模擬結果見圖2。

圖2 超標期內PM2.5濃度值空間模擬結果/(μg·m-3)

由圖2可知，濃度出現超標時，模擬結果顯示區域的極小值為40.23 μg·m-3，極大值為107.75 μg·m-3，3次的極差分別為52.15 μg·m-3、62.74 μg·m-3和52.14 μg·m-3，超標情況下，全域濃度值差異較大，大部分介于50～60 μg·m-3之間；3次濃度超標中，濃度極大值與極差最大均發生在T2。濃度超標時，空間上整體呈現“中心-外圍”濃度值遞減的特征，廣佛交界，中山以北為高濃度典型中心，此外，東莞等靠近該中心的區域會形成次級高值區；結合各濃度等級范圍占全域面積比例(表3)，可知當中心濃度值顯著超標發生時(超標期T2)，濃度等級在良(4)及以上約占了全域3/4的范圍，此類情況下，次級高值區將會向外擴散至更廣范圍(圖2(b))。

表3 各濃度等級占全域面積比例

3.3 時序濃度模擬結果

以廣州市建模過程為例，具體說明建模過程，濃度時序數據一階差分后序列的自相關圖(ACF)和偏自相關圖(PACF)見圖3。據圖可知序列基本位于置信區間，說明序列是平穩的。由圖判斷自相關函數在K=8截尾，偏相關函數在K=3截尾，則p、q值分別為3和8。

圖3 一階差分后自相關、偏自相關函數

根據p、q值初步建立的模型，模型ARIMA(3,1,8)的模型擬合與模擬圖見圖4，ARIMA(3,1,8)R2為0.41，利用2021年10月觀測值進行模擬精度檢驗，10月觀測值為24.22 μg·m-3，模擬值為21.47 μg·m-3，相對誤差為0.11，精度檢驗滿足模擬要求，說明ARIMA的短期模擬精度較高。模擬結果顯示，2021年11月和12月濃度值分別為25.88 μg·m-3和26.36 μg·m-3。

圖4 2019-2021年濃度高值月份濃度值擬合與模擬結果

繪制雷達圖來分析2019-2021年濃度高值易發月份的濃度值時序的演變與同期數值增減規律見圖5。根據圖5的結果，各子區域顯示1～2月的濃度值在全年較高，且1月濃度值明顯高于2月份濃度值，但進入10月以后，相較于上一年同期，各市濃度值均呈現遞減的態勢。

圖5 粵港澳各市及特區2019-2021年濃度高值月份月均濃度狀況/(μg·m-3)

以廣州市為例具體分析，2019年1月和2021年1月濃度值分別為48.30 μg·m-3和46.58 μg·m-3，2020年1月濃度值為 34.20 μg·m-3，增減趨勢不穩定，2月份濃度值變化幅度較小。 進入10月份，2019年濃度值為40.51 μg·m-3，而2020年與2021年濃度值已降至21.86 μg·m-3和24.22 μg·m-3，變化幅度較大，據觀測數據，2020年11～12月濃度，相較2019年同期，呈現遞減的變化，變化幅度在5～10 μg·m-3之間，而根據模擬數據，2021年11～12月的濃度將繼續延續這一規律，濃度值繼續減小。

東莞、深圳和惠州濃度特征在時序上也呈現年初濃度高，增減幅度小的特征，但進入10月后濃度值顯著減小，且降幅也強于廣州市。與珠三角各市相比，香港與澳門年初的濃度值也較高，在30 μg·m-3左右，10～12月數值穩定保持在20 μg·m-3左右，低于同一時段的珠三角各市域。

4 結論與討論

本文以粵港澳大灣區超標期的PM2.5濃度特征為研究對象，一方面，利用氣溶膠與氣象要素數據模擬超標期時大灣區的PM2.5濃度空間特征；另一方面，結合ARIMA模型分析了研究區2019-2021年秋冬季節的濃度值演變情況，得出了以下結果：

在2015年超標期下，PM2.5濃度空間上呈現“中心-外圍”濃度值遞減的特征，全域濃度值介于50～60 μg·m-3之間。當中心濃度值顯著超標發生時，高值范圍會向外擴散。

經過檢驗，ARIMA模型模擬短期城市群PM2.5濃度的精度較高，2019-2021年的觀測與模擬結果顯示，近年來，1～2月的濃度值增減不明顯，數值保持穩定；1月的濃度值為全年最高，需在該時段重點管控PM2.5濃度，各市10～12月濃度值均呈現遞減的趨勢。地域上相鄰的單元其濃度值增減規律具有共性，如廣、肇、佛等市，又如莞、深、惠等市。

由此得到以下結論：相較于研究區全年濃度值在空間上呈“東南-西北”的三級梯度分異[27]，濃度發生超標時，研究區濃度值呈現“中心-外圍”遞減，2019-2021年秋冬季節，研究區僅1月PM2.5濃度值相對較高，需要開展整治專項行動，如綜合運用熱點網格、移動監測車等科技執法手段，在濃度重點防控區域逐步建立“監測-通報-溯源-整改”的長效機制。

盡管如此，由于研究認知水平的局限，本文的研究還有待進一步完善：①本研究選取了與PM2.5濃度具有關聯性的因子構建多元回歸模型，但地形條件、時間序列變化等因素也會影響濃度值[28]。同時，由于粵港澳大灣區空間監測站較少且分布不均衡，將此數據與AOD數據結合開展PM2.5濃度的反演可能會存在不確定性；②空間統計模型適用于中長期的濃度空間靜態評估，而短期的濃度可能會受到各種因素影響，仍用空間統計模型來模擬濃度，可能會與真實情況出現偏差，研究結果對制定濃度管控措施的參考價值有限，因此，對PM2.5濃度空間模擬尚需長期、深入研究。