基于時空指示克里格的PM2.5不確定性分布

梅 楊,張文婷,楊 勇*,趙 玉,李露露 (.寧波市鎮海規劃勘測設計研究院,浙江 寧波 35200；2.華中農業大學資源與環境學院,湖北 武漢 430070；3.農業部長江中下游耕地保育重點實驗室,湖北 武漢 430070)

近年來,隨著國民經濟的持續快速發展,我國PM2.5污染有常態化趨勢[1],研究表明[2-5],若長時間暴露于高濃度 PM2.5中,不僅會引發各種呼吸系統病變,同時還會增加癌癥的發病率,嚴重威脅人類身體健康.當前,對于 PM2.5的研究,主要有以下4個方面：(1)從時間角度出發,運用相關方法研究單個監測站點的PM2.5數據在時間維度上的變化特性,如李梓銘等[6]通過研究北京城區 PM2.5在不同時間尺度上的周期性得出,北京城區 PM2.5濃度存在多個明顯的周期性變化;(2)從空間角度出發,運用相關方法研究不同監測站點的 PM2.5數據在空間上的變化特性,如 Gehrig等[7]研究了瑞士1998～2001年期間 PM2.5與 PM10變化特征,發現PM2.5和PM10存在高相關性,且冬季PM2.5濃度高于夏季;莊欣等[8]研究了珠三角大氣污染的空間分布特征,結果表明在珠三角地區大氣污染存在明顯的區域性特征;(3)從 PM2.5自身角度出發,研究PM2.5的化學組成成分及來源,如周甜等[9]研究了華北平原夏季 PM2.5的化學組成成分及來源,結果指出PM2.5中二次無機離子的含量達到60%,其來源也受到工業源和塵源的影響;(4)從時空的角度出發,運用時空地統計學理論,研究PM2.5數據在時空維度的變化特征,如梅楊等[10]研究了山東省PM2.5的時空演變特征,結果發現PM2.5具有明顯的時間污染特征和空間污染特征;Christakos等[11]在時空克里格的基礎上,提出了一種新的時空降維預測方法,并以山東省 PM2.5為研究對象進行預測并對結果進行分析,證明了該方法優于時空克里格.此外,其他學者也從上述4個方面對PM2.5的時空濃度分布特征進行研究[12-15],但大多數均著力于對 PM2.5的濃度值進行分析,屬于確定性分析范疇,而對于 PM2.5達到某種空氣質量等級的概率性分析,國內外乏善可陳;同時,PM2.5作為大氣污染物的一種,其變化受溫度、濕度、風力等氣象因素影響較大,若僅對PM2.5的濃度值做確定性分析,不僅會使得分析結果不夠嚴謹,同時也不利于環保部門對空氣污染來源的探索以及制定完善的空氣污染防治措施.基于此,本研究擬以山東省 2014年PM2.5日均質量監測濃度為數據源,根據2012年我國環境保護部門頒布的《環境空氣質量標準》[16]和環境保護環境監測司公布的空氣質量指數對人體健康的影響狀況[17]為閾值,擴展空間領域的指示克里格法,提出時空指示克里格法,并對山東省2014年PM2.5日均質量濃度進行不確定性分析.

1 材料與方法

1.1 數據來源及預處理

本研究 PM2.5數據來源于中國空氣質量在線監測分析平臺(http：//www.aqistu dy.cn/),該平臺在山東省分布有99個監測點位(圖1),每小時監測一次,并同步進行數據信息更新.數據監測時間為2014年1月1日0：00時至2014年12月31日23：00.對于獲取的監測數據,由于數據更新周期短,數據量極大,不利于對數據進行時空建模分析,因此對各站點的監測數據按天取均值處理.圖 2為各監測站點PM2.5日均質量濃度基本統計信息.

圖1 山東省空氣質量監測點位分布Fig.1 Spatial distribution of air quality monitoring stations in Shandong Province

圖2 山東省PM2.5質量濃度基本統計特征Fig.2 Statistical characteristics of the PM2.5 concentrations for Shandong Province

1.2 多閾值時空指示克里格方法

作為一種大氣污染物,PM2.5在不同尺度下其變量之間的自相關程度相差較大,且隨著監測站點時間與空間距離的增大,其時空變異函數的隨機成分也逐漸增大,若僅對PM2.5運用單一閾值的時空指示克里格(STIK)進行分析,即只分析 PM2.5在某一閾值下的時空結構特征,不僅掩蓋了其在大尺度下的結構特征,同時也不利于分析其在小尺度下的宏觀結構.因此,本研究采用多閾值時空指示克里格法進行分析.

1.2.1 多閾值時空指示變量與時空指示經驗變異函數 多閾值時空指示克里格(MTSTIK),是指在多個時空閾值下,分別運用時空指示克里格法對時空變量進行時空指示預測分析[18].在MTSTIK中,定義D為空間域,T為時間域,若在研究區域存在指示值(閾值)Zc(c=1,2,…,n),則對于時空域D×T上的每一點Z(s,t)(s=(x,y)∈D,t∈T),其時空指示變量值I c (s , t; Z c )為：

對于各閾值下的時空指示變量,其時空指示經驗變異函數與時空普通克里格經驗變異函數

[10]類似,即在本征假設的條件下,定義時空指示變異函數為：

式中：hS,hT分別為空間和時間間隔變量,N(hS,hT)為樣點中符合所定義空間和時間間隔的點對數.

1.2.2 時空理論變異函數模型 時空理論變異函數模型通常分為時空分離模型和時空非分離模型[19],兩者的差別在于時空分離模型中存在獨立的時間和空間部分,而時空非分離模型中則將時間與空間統一考慮.在本次研究中,PM2.5質量濃度變化受時間和地域的共同影響,若將 PM2.5質量濃度變化分別在時間維度和空間維度進行分解,不僅與實際情況不符,其分析結果也不具備科學性和嚴謹性.因此,本研究采用 Cressie等[20]提出的時空協方差函數模型(式 3)對各閾值下的時空經驗指示變異函數進行擬合.

1.2.3 時空指示克里格預測 對于各閾值下的時空指示變量,其時空指示克里格插值理論基礎同普通指示克里格插值算法類似,即構造克里格權重計算矩陣A?λ= b,然后通過時空變異函數計算A、b矩陣,求取時空權重,進而利用求取待估點小于閾值Zc的概率[19].

1.3 PM2.5時空分布的不確定性分析

根據各閾值對應的時空預測結果,繪制相應的時空立方體.并統計各空間點位對應的時空預測結果,分析各時空立方體結構特征以及計算各空間點位對應的年均概率值、月均概率值和部分城市城區范圍內日均概率值,實現對山東省PM2.5時空分布的不確定性分析.

2 結果與分析

2.1 時空指示克里格閾值

表1 PM2.5濃度值對應的空氣質量級別及對健康的影響狀況Table 1 Air quality level and health effects to the concentration of PM2.5

MTSTIK關鍵在于如何確定多個既能反映研究區整體的分布特征,又能揭示局部時空結構特征的時空閾值.在本文中,根據我國環境保護部2012年頒布的《環境空氣質量標準》中規定的日均 PM2.5濃度標準值和我國環境保護部環境監測司2012年公布的空氣質量指數對人體健康的影響情況,依次選擇 35,75,150和 250μg/m3作為本次研究的時空指示值(表1).

2.2 時空指示經驗變異函數與模型擬合

按照公式(1)和對應的時空指示值,將山東省PM2.5質量濃度指示化.同時,設置空間步長為5km,最大空間計算距離為100km,時間步長為1d,最大時間計算距離為 14d,運用公式(2)計算各閾值下的時空指示經驗變異函數,并根據公式(3)對各時空模型進行擬合(表2,圖3).

表2 4種時空指示理論變異函數擬合參數Table 2 Parameters in 4 different spatio-temporal semivariogram models

由表2和圖3可知,在空間距離大于100km時,各時空指示克里格變異函數值仍有增大趨勢,說明山東省PM2.5濃度變化的空間自相關范圍大于 100km;在時間距離在 3d時,各變異函數值已趨于穩定,說明山東省 PM2.5濃度變化的時間自相關范圍為3d左右.

圖3 4種時空指示理論模型擬合Fig.3 Spatio-temporal indicator semivariograms under 4thresholds

2.3 多閾值時空指示克里格預測立方體

設置MTSTIK預測格網大小為2km×2km×1d,運用時空指示克里格方法,對各閾值下的PM2.5指示值進行時空插值.并基于插值結果,在三維XYZ立方體中,以XOY平面為空間坐標,Z軸為時間坐標,繪制出各閾值對應的時空指示克里格預測立方體(圖4).

圖4 4種閾值下時空指示預測立方體Fig.4 The estimated cubes obtained by MTSTIK under 4thresholds

表3 四種時空立方體基本統計特征Table 3 Statistical characteristics of spatiotemporal indicators of PM2.5 estimated cubes under 4thresholds

對于各時空立方體,進行基本統計特征分析.由表3和圖5可知,對于閾值Zc=35μg/m3,時空立方體數據值的平均值為 0.1309,四分位距為0.0593,數據值小于 0.2 的概率(ρ)為 81%,數據值大于0.8的概率為7%,說明其時空立方體數據值小,數據離散程度低,境內全年空氣質量以超過0.8的概率達到空氣質量優級別的時空占比為7%;對于閾值Zc=75μg/m3,時空立方體數據值的平均值為0.4945,四分位距為0.8782,數據值小于0.2的概率和大于0.8的概率均為34%,其余3個統計層次的概率值均為 10%左右,說明時空立方體數據值分布較為對稱,離散程度高,境內全年空氣質量以超過0.8的概率達到和超過輕度污染級別的時空占比均為 34%;對于閾值Zc=150μg/m3和Zc=250μg/m3,其時空立方體數據值的平均值分別為0.8111和0.9545,四分位距依次為0.2728和0,數據值大于0.8的概率均依次為74%和93%,說明其時空立方體數據值大,數據離散程度低,境內全年空氣質量以超過0.8的概率超過嚴重污染級別的時空占比為1%.

圖5 4種時空立方體日均概率百分比Fig.5 Percentages of daily average probability in estimated cubes under 4thresholds

2.4 PM2.5時空分布不確定性分析

2.4.1 山東省年均概率不確定性分析 根據各閾值下的時空預測立方體數據,獲取各空間點位的日概率值,取其年平均值,繪制相應的山東省PM2.5年均概率分布(圖6所示).

由圖6可知,對于閾值Zc=35μg/m3,山東省境內東部沿海地域年均概率為 0.3～0.4,局部地域大于0.5,其余地域均小于 0.2,說明山東省大部分地域空氣質量達到優級別的概率小于 0.2;對于閾值Zc=75μg/m3,東部沿海地域年均概率為0.6～0.9,其余地域為0.4～0.6,說明東部沿海空氣質量達到輕度污染及以上級別的概率均大于0.6,其余地域大于0.4;對于閾值Zc=150μg/m3和Zc=250μg/m3,全境年均概率均大于0.8,部分地域接近于1,說明山東省境內空氣質量超過嚴重污染級別的概率小于0.1.

此外,各尺度東部沿海空氣質量達到各空氣質量級別的概率均大于中部和西部,結合莊欣等[8]在對珠三角 PM2.5污染情況分析以及山東省地理位置和地理環境可知,受大氣氣流運輸的影響,西部地區與河北、河南接壤,而陸雅靜等[21]研究發現,河北省空氣質量污染嚴重,從而使得山東省西部PM2.5污染較為嚴重;而東部則臨海,距離河南河北等省較遠,外省氣流抵達東部區域時,空氣中的顆粒物經過自然沉降后,氣流攜帶顆粒物濃度降低,造成東部受到外省污染情況的影響遠小于西部.

2.4.2 山東省月均概率不確定性分析 由上述分析可知,雖在不同閾值下山東省年均概率分布不確定性結構特征不盡相同,但在后續對各閾值下PM2.5月均濃度不確定性結構特征進行分析后發現,各閾值對應概率的時空分布模式大體相似,同時根據表 1空氣質量分指數對健康的影響情況,在避免重復分析各閾值下時空不確定性特征的條件下,本文選取閾值Zc=75μg/m3的時空預測結果來分析山東省PM2.5月均概率和部分城市日均概率的時空分布特征(圖7).

由圖 7可知,從空間上,各月份達到 PM2.5≤75μg/m3的最小概率均位于西部區域(如菏澤、聊城、濟寧等),尤其是在1月份,部分城市概率小于0.1;各月份概率最大位于東部沿海區域(青島、煙臺、威海等),絕大部分月份概率超過 0.6;中部為過渡區域,不同月份概率變化較大整體上呈現出從西至東,概率值逐漸增大,說明山東省空氣質量從西至東,污染程度逐漸減輕,具有明顯的空間分布特征,這與梅楊[10]在研究山東省 PM2.5質量濃度時空分布時得到的結論相同.

從時間上,各區域達到 PM2.5≤75μg/m3的概率值最大為7、8月份,大部分地域概率值大于0.6,部分地域概率超過 0.9;概率值最小為 1月份,大部分地域概率小于0.2,局部地域介于0.5～0.6之間.整體上,從 1～12 月,概率值先增大后減小,說明山東省空氣污染程度從 1～12月,先減輕后加劇,具有明顯的時間分布特征,這與康桂紅等[22]對山東省2014～2015年間,不同月份PM2.5污染程度分析后得到的結果,以及成亞利對上海市 2014年PM2.5污染情況分析得到的結果相一致[23].

2.4.3 山東省部分城市單一閾值日均概率不確定性分析 根據山東省行政區劃和城市經濟排名,依次在山東省北部、西部、中部、南部、中東部以及東部沿海各選擇一個城市作為PM2.5日均質量濃度不確定性分析觀測點(北部選擇東營市,西部選擇濟寧市,中部選擇濟南市,南部選擇臨沂市,中東部選擇濰坊市,東部沿海選擇煙臺市),獲取各城市城區范圍在閾值 Zc=75μg/m3的日均概率值,并繪制相應的百分比(圖8).

圖6 多閾值時空指示克里格年均概率分布Fig.6 Spatial distribution of annual average probability by MTSTIK under various thresholds

圖7 PM2.5≤75μg/m3下時空指示克里格月均概率分布特征Fig.7 Spatial distribution of monthly average probabilities for PM2.5≤75μg/m3

圖8 山東省部分城市日均概率百分比Fig.8 Percentages of daily average probabilities for hot-spot cities in Shandong province

由圖 8 可知,對于 PM2.5≤75μg/m3,6 個城市全年日均概率大于 0.75的占比中,最高為煙臺市,超過70%,最低為濟南市,僅為36%,說明煙臺、東營、濟寧、臨沂、濰坊和濟南6個城市全年的空氣污染程度以0.75的概率達到輕度污染的天數最高為煙臺市,260d左右,最低為濟南市,僅有130d;全年日均概率小于0.25的占比中,除煙臺市(13%)以外,其余五個城市比例均高于 30%(東營市最高,為 41%),說明全年空氣污染程度以大于0.75的概率超過輕度污染的天數依次為48,150,139,130,105d.各城市空氣污染程度依次為：煙臺＜東營＜濰坊＜濟南＜濟寧≈臨沂,大致符合從西至東,空氣污染逐漸減輕的特點.此外,根據山東省《2014山東省環境狀況公報》[24]顯示,6個城市年均 PM2.5濃度大小順序依次為：煙臺＜濰坊＜東營＜濟寧＜臨沂＜濟南[24],說明濟南市PM2.5日均濃度相較于其他5個城市,其波動最為明顯.

3 結論

3.1 山東省境內 PM2.5的空間自相關范圍超過100km,時間自相關范圍為3d左右.

3.2 對于閾值 Zc=35μg/m3,其時空預測立方體數據值偏小、離散程度低,各空間點位以超過0.8的概率達到空氣質量優級別的時空占比為 7%;對于閾值 Zc=75μg/m3,其時空立方體數據值呈對稱分布,離散程度高,其以超過 0.8的概率達到空氣質量輕度污染級別的時空占比為 34%;對于閾值 Zc=150μg/m3和 Zc=250μg/m3,其時空立方體數據值較大、離散程度低,其以超過0.8的概率超過嚴重污染級別的時空占比為1%.

3.3 在年均概率分布上,山東省大部分地域空氣質量達到優級別的概率介于 0～0.1;達到輕度污染級別的概率為東部大于 0.6,其余地區大于0.4;超過嚴重污染級別的概率均小于0.1.

3.4 從空間上,各月份達到 PM2.5≤75μg/m3的最小概率均位于西部區域,概率最大位于東部沿海區域,中部為過渡區域,.整體上呈現出從西至東,概率值逐漸增大;從時間上,月均概率值最大為7、8月份,概率值最小為1月份.整體上,從1～12月,概率值先增大后減小.

3.5 在濟寧、濟南、臨沂、東營、濰坊和煙臺6個城市中,全年空氣污染程度以大于0.75的概率值達到和優于輕度污染的天數中,煙臺市最高(超過260d),濟南市最低(130d左右);以小于0.75的概率值超過輕度污染的天數中,最高為東營(150d左右),最低為煙臺(48d).

[1]王庚辰,王普才.中國 PM2.5污染現狀及其對人體健康的危害[J]. 科技導報, 2014,32(26)：72-78.

[2]戴海夏,宋偉民.大氣 PM2.5的健康影響 [J]. 環境衛生學雜志,2001,28(5)：299-303.

[3]張文麗,崔九思,徐東群,等.大氣中細顆粒物的污染特征及其生物效應 [J]. 中國環境衛生, 2003,(1)：90-102.

[4]覃知還.關于PM2.5的十個問答 [J]. 證券導刊, 2011,(46)：93-96.

[5]Weichenthal S, Villeneuve P J, Burnett R T, et al. Long-term exposure to fine particulate matter： association with nonaccidental and cardiovascular mortality in the agricultural health study cohort[J]. Environmental health perspectives, 2014,122(6)：609-15.

[6]李梓銘,孫兆彬,邵 勰,等.北京城區PM2.5不同時間尺度周期性研究 [J]. 中國環境科學, 2017,37(2)：407-415.

[7]Gehrig R, Buchmann B. Characterising seasonal variations and spatial distribution of ambient PM10and PM2.5concentrations based on long-term Swiss monitoring data [J]. Atmospheric Environment, 2003,37(19)：2571-2580.

[8]莊 欣,黃曉鋒,陳多宏,等.基于日變化特征的珠江三角洲大氣污染空間分布研究 [J]. 中國環境科學, 2017,37：2001-2006.

[9]周 甜,閆才青,李小瀅,等.華北平原城鄉夏季PM2.5組成特征及來源研究 [J]. 中國環境科學, 2017,37(9)：3227-3236.

[10]梅 楊,黨麗娜,楊 勇,等.基于時空克里格的PM2.5時空預測及分析 [J]. 環境科學與技術, 2016,39(7)：157-163.

[11]Christakos G, Yang Y, Wu J P, et al. Improved space-time mapping of PM2.5distribution using a domain [J]. Ecological Indicators, 2018,85：1273-1279.

[12]Begum B A, Biswas S K, Hopke P K. Temporal variations and spatial distribution of ambient PM2.5and PM10concentrations in Dhaka, Bangladesh [J]. Scie of the Total Envi [J], 2006,358(1-3)：36—45.

[13]Yu H L, Chiting C, Lin S D, et al. Spatiotemporal analysis and mapping of oral cancer risk in Changhua County (Taiwan)： an application of generalized Bayesian maximum entropy method.[J]. Annals of Epidemiology, 2010,20(2)：99.

[14]Hu J, Wang Y, Ying Q, et al. Spatial and temporal variability of PM2.5, and PM10, over the North China Plain and the Yangtze River Delta, China [J]. Atmospheric Environment, 2014,95(1)：598-609.

[15]Kloog I, Chudnovsky A A, Just A C, et al. A New Hybrid Spatio-Temporal Model For Estimating Daily Multi-Year PM2.5Concentrations Across Northeastern USA Using High Resolution Aerosol Optical Depth Data [J]. Atmospheric Environment, 2014,95(1)：581-590.

[16]中國環境科學研究院.環境空氣質量標準 [M]. 北京：中國環境科學出版社, 2012.

[17]環境保護部環境監測司,中國.“十二五”環境監測工作手冊 [M].中國環境科學出版社, 2012.

[18]Yang Y, Christakos G. Uncertainty assessment of heavy metal soil contamination mapping using spatiotemporal sequential indicator simulation with multi-temporal sampling points [J].Environmental Monitoring & Assessment, 2015,187(9)：571.

[19]楊 勇,梅 楊,張楚天,等.基于時空克里格的土壤重金屬時空建模與預測 [J]. 農業工程學報, 2014,30(21)：249-255.

[20]Cressie N, Huang H C. Classes of Nonseparable, Spatio-Temporal Stationary Covariance Functions [J]. Publications of the American Statistical Association, 1999,94(448)：1330-1339.

[21]陸雅靜,王 輝,倪爽英,等.河北省 2013～2014年環境空氣質量現狀對比分析 [J]. 安徽農業科學, 2015,(14)：271-274.

[22]康桂紅,孫興池,韓永清,等.山東省大氣污染時空分布特征分析[J]. 山東氣象, 2016,36(1)：13-17.

[23]成亞利,王 波,上海市 PM2.5的時空分布特征及污染評估 [J].計算機與應用化學, 2014,31(10)：1189-1192.

[24]山東省環保保護廳.2014年山東省環境狀況公報 [EB/OL].http：//www.zhb.gov.cn/gkml/hbb/qt/201506/t20150604_302942.html.