基于GTWRK的PM2.5濃度空間插值方法

謝劭峰，張 偉，潘夢清，黃良珂*，劉立龍

(1.桂林理工大學(xué) 測繪地理信息學(xué)院，廣西 桂林 541006；2.河南省地震局 鶴壁地震監(jiān)測中心站，河南 鶴壁 456200；3.桂林市測繪研究院，廣西 桂林 541002)

0 引言

隨著我國經(jīng)濟的快速發(fā)展和城市化進程的加快，大氣細顆粒物PM2.5已經(jīng)成為影響我國大氣環(huán)境污染的主要因素之一[1]。目前PM2.5的監(jiān)測以定點的地面監(jiān)測為主，由于監(jiān)測站點分散，僅能觀測有限空間范圍內(nèi)的空氣污染情況，無法獲取區(qū)域內(nèi)連續(xù)的PM2.5空間分布格局[2]。

為了獲取連續(xù)的PM2.5空間分布，文獻[3-5]采用克里金方法進行PM2.5濃度空間插值。然而由于PM2.5具有很強的空間異質(zhì)性，因此僅僅依靠普通克里金插值獲得PM2.5空間分布，效果不太理想[2]。地理加權(quán)回歸(Geographically Weighted Regression，GWR)模型[6-8]可以在一定程度上解決PM2.5濃度估算的空間異質(zhì)性問題，但常規(guī)GWR模型忽略了時間非平穩(wěn)性的影響。針對這個問題，文獻[9-11]將時間因素加入到GWR模型中，采用時空地理加權(quán)回歸(Geographically and Temporally Weighted Regression,GTWR)方法研究PM2.5濃度的時空分布。盡管GTWR模型具有較高的擬合優(yōu)度，但GTWR模型難以解決PM2.5分布的復(fù)雜非線性問題，對回歸殘差也未進行進一步的處理。

土地利用回歸(Land Use Regression,LUR)模型是PM2.5時空分布模擬的有效方法。文獻[12-14]基于PM2.5濃度與土地利用、地形、氣候、交通和人口密度等因素之間的相關(guān)關(guān)系，使用LUR模型對PM2.5濃度的空間分布進行反演，較好地反映了區(qū)域PM2.5濃度分布情況，但LUR模型難以解決PM2.5分布的時空非平穩(wěn)性問題。還有研究人員采用機器學(xué)習方法進行PM2.5濃度預(yù)測或反演：文獻[15]采用隨機森林方法估計PM2.5濃度；文獻[16]基于氣溶膠光學(xué)厚度(Aerosol Optical Depth，AOD)產(chǎn)品和氣象再分析資料，采用隨機森林等方法反演北京市近地面PM2.5濃度；文獻[17]利用卷積神經(jīng)網(wǎng)絡(luò)預(yù)報模型進行京津冀區(qū)域PM2.5濃度預(yù)報；文獻[18]構(gòu)建了基于卷積神經(jīng)網(wǎng)絡(luò)和長短時記憶的深度學(xué)習預(yù)報模型。這些機器學(xué)習方法也取得了較好的預(yù)測結(jié)果，但同樣無法有效處理PM2.5時空分布的非平穩(wěn)性問題。

上述研究基于不同方法、不同數(shù)據(jù)進行區(qū)域PM2.5濃度預(yù)測或模擬，取得了較好的效果。但這些方法難以同時有效處理PM2.5分布的時空非平穩(wěn)性、空間自相關(guān)性問題。因此，本文基于GTWR方法處理時空非平穩(wěn)性問題和克里金方法處理空間自相關(guān)性問題的良好能力，建立了顧及多元因子影響的時空地理加權(quán)回歸克里金(Geographically and Temporally Weighted Regression Kriging,GTWRK)模型，進行廣西PM2.5濃度預(yù)測，并對比分析了模型的預(yù)測精度，繪制了2015—2019年廣西PM2.5濃度空間分布圖。

1 數(shù)據(jù)與方法

1.1 研究區(qū)概況

廣西地處中國華南地區(qū)西部，位于20°54′～26°24′N、104°26′～112°04′E，陸地總面積約23.67萬平方千米，南瀕熱帶海洋，北為南嶺山地，西延云貴高原，境內(nèi)河流縱橫，地理環(huán)境比較復(fù)雜。廣西氣候類型多樣，夏長冬短，年均溫在16～23 ℃。從全國來看，盡管廣西不是霧霾高發(fā)地區(qū)，但部分地市出現(xiàn)霧霾的天數(shù)比較高，尤其是梧州、桂林和柳州3市霧霾天數(shù)出現(xiàn)的頻率超過20%[19]。

1.2 數(shù)據(jù)來源及預(yù)處理

PM2.5濃度數(shù)據(jù)來源于中國環(huán)境監(jiān)測總站(http:∥www.cnemc.cn/)，時間跨度為2015年1月1日—2019年12月31日。得到的數(shù)據(jù)為逐時監(jiān)測的站點數(shù)據(jù)，在計算得到日均值的基礎(chǔ)上進一步得到建模使用的月均值和年均值。

氣象數(shù)據(jù)來源于中國氣象數(shù)據(jù)網(wǎng)(http:∥data.cma.cn/)，獲取的數(shù)據(jù)主要包括降水量(Precipitation，PRE)、風速(Wind Speed，WS)、氣壓(Pressure，PRS)、溫度(Temperature，TEM)、相對濕度(Relative humidity，RHU)和水汽壓(Vapor Pressure，VP)。

由于氣象站點和空氣質(zhì)量監(jiān)測站點的空間位置不同且數(shù)量也不同，本文對氣象數(shù)據(jù)進行反距離加權(quán)插值(Inverse Distance Weighting，IDW)到空氣質(zhì)量監(jiān)測站點位置處。氣象站點共有17個，其中百色市有3個，河池市、柳州市和來賓市各有2個，桂林市、賀州市、梧州市、玉林市、欽州市、南寧市、防城港市和崇左市各有1個。空氣質(zhì)量監(jiān)測站點共有50個，其中南寧市有8個，柳州市有6個，桂林市、梧州市、貴港市和北海市各有4個，河池市、欽州市、防城港市和玉林市各有3個，百色市、崇左市、賀州市和來賓市各有2個。

AOD數(shù)據(jù)來源于美國國家航天局(https:∥ladsweb.modaps.eosdis.nasa.gov/)，數(shù)據(jù)為MCD19A2，空間分辨率為1 km，時間分辨率為1 d。

全球?qū)Ш叫l(wèi)星系統(tǒng)(Global Navigation Satellite System，GNSS)天頂對流層延遲(Zenith Tropospheric Delay,ZTD)數(shù)據(jù)來源于中國地震局GNSS數(shù)據(jù)產(chǎn)品服務(wù)平臺(ftp:∥ftp.cgps.ac.cn/products)。由于陸態(tài)網(wǎng)在廣西僅有6個基準站，為了便于研究，本文將與廣西相鄰(包括湖南、云南、貴州、廣東和海南)的共計44個GNSS基準站的ZTD數(shù)據(jù)分月份對廣西進行反距離加權(quán)插值，最后得到廣西空氣質(zhì)量監(jiān)測站點月均的ZTD值。

高程數(shù)據(jù)(Digital Elevation Model,DEM)來源于中國科學(xué)院地理科學(xué)與資源研究所(http:∥www.resdc.cn/Default.aspx)；植被指數(shù)(Normalized Difference Vegetation Index,NDVI)來源于美國國家航天局(https:∥ladsweb.modaps.eosdis.nasa.gov/)；人口數(shù)據(jù)來源于Worldpop網(wǎng)站(https:∥www.worldpop.org/)發(fā)布的1 km的柵格數(shù)據(jù)集，時間范圍為2015—2019年。

由于AOD數(shù)據(jù)的空間分辨率為1 km×1 km，所以將所有數(shù)據(jù)全部轉(zhuǎn)化成相同的空間分辨率。為了保持所有數(shù)據(jù)的空間范圍一致，本文將所有柵格數(shù)據(jù)集的坐標系都轉(zhuǎn)化為CGCS2000-Alberts投影坐標系，用以準確獲得空氣質(zhì)量監(jiān)測站點上的各個特征變量。

1.3 研究方法

1.3.1 地理加權(quán)回歸

GWR是一種改進的空間線性回歸模型，其系數(shù)不是固定的，而是隨著空間位置變化，并且將局部特征作為權(quán)重。GWR模型可以表示為[20]：

(1)

式中，yi為在位置i處的因變量值；xik(k=1,2,…,p)為位置i處的自變量值；(ui,νi)為回歸分析點i的坐標；β0(ui,νi)為截距項；βk(ui,νi) (k=1,2,…,p)為回歸系數(shù)；εi為回歸殘差。

GWR模型的回歸系數(shù)采用加權(quán)最小二乘法求解[20]：

(2)

式中，W(ui,νi)為m×m的空間權(quán)重對角矩陣；X為m×(n+1)自變量的設(shè)計矩陣；Y為m×1因變量向量。

空間權(quán)重采用Bi-square函數(shù)計算[20]：

(3)

式中，wij為空間已知點j去估計未知點i時的權(quán)重；dij為點i與j之間的歐氏距離；h為帶寬。

帶寬采用最小AICc準則(corrected Akaike Information Criterion)進行判斷[20]：

AICc=-nln(σ′)+nln(2π)+

n{[n+tr(S)]/[n-2-tr(S)]}，

(4)

式中，n為數(shù)據(jù)點的數(shù)量；σ′為誤差項(估計標準偏差)；tr(S)為帽子矩陣S的跡，是帶寬h的函數(shù)。

1.3.2 時空地理加權(quán)回歸

將時間效應(yīng)整合到GWR模型中，就成為了GTWR模型，可以用來同時處理時間和空間非平穩(wěn)問題。GTWR可以表示為[11]：

(5)

式中，(ui,vi,ti)為第i個樣本點的時空坐標。

與GWR模型相似，回歸系數(shù)采用加權(quán)最小二乘法求解：

(6)

引入時空權(quán)重矩陣W(ui,νi,ti)來衡量樣本i在估計樣本j方面相對于時間和空間的權(quán)重。時空距離dST為空間距離dS與時間距離dT的線性組合[11]：

(7)

式中，ti和tj為觀測時間i和j；λ和μ為用于協(xié)調(diào)空間距離和時間之間不同單位的影響權(quán)重。

1.3.3 克里金法

克里金(Kriging)插值是建立在變異函數(shù)理論及結(jié)構(gòu)分析基礎(chǔ)上，利用區(qū)域化變量的原始數(shù)據(jù)和變異函數(shù)的結(jié)構(gòu)性，對未知點區(qū)域化變量的取值進行線性無偏最優(yōu)估計的一種方法?？死锝鸱ㄊ且环N空間局部插值方法，適用于具有空間相關(guān)性的區(qū)域化變量研究。克里金插值可表示為[20]：

(8)

式中，z*(x0)為插值點x0處的估計值；z(xi)為n個樣本點中xi處的觀測值；λi為對應(yīng)的權(quán)重。

1.3.4 時空地理加權(quán)回歸克里金

GTWRK是在GTWR的基礎(chǔ)上，對其回歸殘差進行普通克里金(Ordinary Kriging，OK)插值，最后將OK殘差插值結(jié)果和GTWR預(yù)測值相加而得。GTWRK可表示為：

yGTWRK(μi,νi,ti)=yGTWR(μi,νi,ti)+εOK(μi,νi,ti)，

(9)

式中，yGTWRK(μi,νi,ti)為GTWRK的預(yù)測值；yGTWR(μi,νi,ti)為GTWR的預(yù)測值；εOK(μi,νi,ti)為GTWR回歸殘差進行普通克里金插值的結(jié)果。

2 模型構(gòu)建與精度評估

2.1 相關(guān)性分析與特征變量選擇

為了研究廣西地區(qū)影響PM2.5濃度變化的因素，根據(jù)相關(guān)研究選擇氣象數(shù)據(jù)(包括降水量(PRE)、風速(WS)、氣壓(PRS)、溫度(TEM)、相對濕度(RHU)和水汽壓(VP))、高程數(shù)據(jù)(DEM)、氣溶膠數(shù)據(jù)(AOD)、人口密度(POP)、植被指數(shù)(NDVI)和GNSS ZTD數(shù)據(jù)為特征變量，使用SPSS軟件的雙變量分析功能分析各個特征變量與PM2.5濃度之間的相關(guān)性，結(jié)果如表1所示。

表1 各變量與PM2.5濃度的相關(guān)性

從表1可知，NDVI，PRS，DEM與PM2.5濃度的相關(guān)性不顯著；AOD，POP與PM2.5濃度呈顯著正相關(guān)，AOD的相關(guān)系數(shù)達0.37；PRE，WS，TEM，RHU，VP，ZTD與PM2.5濃度呈顯著負相關(guān)，其中VP的相關(guān)系數(shù)達-0.74。

然而，并不是所有相關(guān)性顯著的變量都適用于GTWRK模型的建模預(yù)測，因此這些變量還需要經(jīng)過共線性診斷。方差膨脹因子(Variance Inflation Factor,VIF)診斷法常用來檢驗各個變量之間是否存在多重共線性。VIF值越大，變量之間的多重共線性越大。本文以7.5為界限，VIF值高于7.5的變量可以判斷為存在共線性的變量。

為了檢驗變量之間的共線性，以相關(guān)性較高的變量為自變量、PM2.5濃度為因變量，在ArcGIS中進行普通最小二乘(Ordinary Least Square，OLS)建模，OLS建模結(jié)果能夠檢驗出變量之間是否存在共線性。剔除共線性較高的變量前、后結(jié)果分別如表2和表3所示。

表2 剔除變量前OLS的建模結(jié)果

表3 剔除變量后OLS的建模結(jié)果

由表2和表3可知,剔除共線性較高的變量后各變量的VIF值均小于7.5，說明變量之間已經(jīng)不存在多重共線性，表3中所有變量可以用來建模。變量系數(shù)的正負也反映了變量對PM2.5濃度的影響方向，AOD和POP與PM2.5濃度呈正相關(guān)，PRE，WS，TEM和ZTD與PM2.5濃度呈負相關(guān)，與上述的皮爾遜相關(guān)系數(shù)分析的影響方向相同，且變量系數(shù)的大小也說明了對PM2.5濃度影響的強弱關(guān)系。從表3還可以看出，對PM2.5濃度正向影響較大的是AOD，負向影響較大的是WS。

2.2 模型構(gòu)建

GWR模型在ArcGIS軟件工具箱里空間關(guān)系建模中的地理加權(quán)回歸模塊實現(xiàn)，GTWR模型在ArcGIS軟件的GTWR插件實現(xiàn)。將POP，AOD，PRE，WS，TEM和ZTD選為解釋變量，PM2.5選為因變量。2個模型的建模結(jié)果如表4所示。

表4 GWR和GTWR模型建模結(jié)果

表4中的R2表示擬合優(yōu)度，值范圍0～1，值越接近于0模型擬合越差，反之表明模型的擬合效果越好。一般情況下，模型參考的是調(diào)整后的R2，GTWR模型調(diào)整后的R2更接近1,說明GTWR模型擬合效果更好；GTWR模型的AICc的值比GWR模型小，說明GTWR模型比GWR模型更優(yōu)；GTWR模型的殘差平方和明顯小于GWR模型的殘差平方和；Sigma表示剩余平方和除以殘差的有效自由度的平方根，該統(tǒng)計值越小越好，表中的GTWR模型的Sigma值也小于GWR模型。綜上，GTWR模型結(jié)果的各項參數(shù)均優(yōu)于GWR模型，說明GTWR模型更適用于PM2.5濃度預(yù)測。

2.3 精度評估

為了檢驗?zāi)Ｐ偷念A(yù)測精度，隨機選擇了80%的站點數(shù)據(jù)進行訓(xùn)練，20%的站點數(shù)據(jù)進行驗證。為了保證所選擇的訓(xùn)練樣本和驗證樣本的科學(xué)性，根據(jù)均勻分布原則選取建模樣本和驗證樣本，即每個月選擇了不同城市的不同站點作為訓(xùn)練集和驗證集。

為了驗證模型的預(yù)測精度，采用平均絕對誤差(Mean Absolute Error，MAE)、均方根誤差(Root Mean Square Error，RMSE)和決定系數(shù)(R2)3個指標來評價各個模型的預(yù)測精度，結(jié)果如表5所示。

表5 模型精度對比

從表5可知，加入時間效應(yīng)的GTWR模型和GTWRK模型優(yōu)于對應(yīng)的未加入時間效應(yīng)的GWR模型和GWRK模型；殘差經(jīng)過克里金插值修正后預(yù)測精度優(yōu)于未進行殘差修正的模型；GTWRK模型的各項評價指標均優(yōu)于其他模型，表明GTWRK比其他模型具有更好的預(yù)測精度。

用Matlab軟件繪制4個模型的實測PM2.5濃度與預(yù)測PM2.5濃度的散點圖，如圖1所示。由圖1可知，相比于其他模型，GTWRK模型的預(yù)測與實測散點圖更接近于一條直線，該模型的R2=0.814，說明實測濃度與預(yù)測濃度擬合程度高，預(yù)測值更接近實測值，進一步證明了GTWRK模型的預(yù)測能力。

(a) GWR模型

(b) GTWR模型

(c) GWRK模型

(d) GTWRK模型

3 廣西PM2.5濃度空間分布特征

為了直觀地顯示廣西PM2.5濃度的時空分布特征，利用ArcGIS中的柵格計算器功能，以年為時間尺度進行PM2.5濃度年均值計算，將GTWRK模型的PM2.5濃度預(yù)測結(jié)果進行空間化表達，得到2015—2019年廣西PM2.5濃度的空間分布，結(jié)果如圖2所示。

(a) 2015年

(b) 2016年

(c) 2017年

(d) 2018年

(e) 2019年

由圖2可知，2015年高濃度PM2.5主要分布在百色、柳州和桂林，其他地區(qū)PM2.5濃度相對較低。2016年高濃度PM2.5主要分布在柳州和桂林，百色的高濃度PM2.5分布面積減少。2017年高濃度PM2.5主要分布在柳州和桂林，但是分布面積較上年有所減少。2018年高濃度PM2.5分布面積明顯減少，主要還是集中在柳州和桂林；來賓較上年有所增加。2019年高濃度PM2.5主要位于柳州、桂林、來賓，但高濃度PM2.5分布面積顯著減少；防城港、北海、欽州、崇左、南寧、河池等市PM2.5濃度較低?？傮w來看，2015—2019年廣西高濃度PM2.5主要位于柳州和桂林，但分布面積逐年減少，說明廣西環(huán)境治理初見成效。

4 結(jié)束語

通過構(gòu)建顧及多因子影響的GTWRK模型進行PM2.5濃度預(yù)測和空間分布制圖，得到以下結(jié)論：

① 廣西地區(qū)影響PM2.5濃度變化顯著性較高的變量有AOD,POP,PRE,WS,TEM,RHU,VP和ZTD，正向影響最大的變量是AOD，負向影響最大的變量是VP。

② GTWRK模型更適用于區(qū)域PM2.5濃度預(yù)測，其預(yù)測結(jié)果的MAE,RMSE,R2分別為4.69 μg/m3,6.44 μg/m3,0.81，全部優(yōu)于GWR,GWRK,GTWR,GTWRK模型。

③ 2015—2019年廣西高濃度PM2.5分布面積逐年遞減，2015年分布面積最大，2019年分布面積最??；高濃度PM2.5主要分布在桂林和柳州，低濃度主要分布于防城港和北海。

由于本文氣象參數(shù)和ZTD等數(shù)據(jù)獲取站點數(shù)量不一，只能通過插值法進行區(qū)域數(shù)據(jù)獲取，因此精度還可以進一步優(yōu)化。