謝劭峰，張 偉，潘夢(mèng)清，黃良珂*，劉立龍

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

0 引言

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

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

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

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

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

1.1 研究區(qū)概況

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

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

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

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

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

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

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

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

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

1.3 研究方法

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

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

(1)

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

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

(2)

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

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

(3)

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

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

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

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

(4)

式中，n為數(shù)據(jù)點(diǎn)的數(shù)量；σ′為誤差項(xiàng)(估計(jì)標(biāo)準(zhǔn)偏差)；tr(S)為帽子矩陣S的跡，是帶寬h的函數(shù)。

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

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

(5)

式中，(ui,vi,ti)為第i個(gè)樣本點(diǎn)的時(shí)空坐標(biāo)。

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

(6)

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

(7)

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

1.3.3 克里金法

克里金(Kriging)插值是建立在變異函數(shù)理論及結(jié)構(gòu)分析基礎(chǔ)上，利用區(qū)域化變量的原始數(shù)據(jù)和變異函數(shù)的結(jié)構(gòu)性，對(duì)未知點(diǎn)區(qū)域化變量的取值進(jìn)行線(xiàn)性無(wú)偏最優(yōu)估計(jì)的一種方法。克里金法是一種空間局部插值方法，適用于具有空間相關(guān)性的區(qū)域化變量研究?？死锝鸩逯悼杀硎緸閇20]：

(8)

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

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

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

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

(9)

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

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

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

為了研究廣西地區(qū)影響PM2.5濃度變化的因素，根據(jù)相關(guān)研究選擇氣象數(shù)據(jù)(包括降水量(PRE)、風(fēng)速(WS)、氣壓(PRS)、溫度(TEM)、相對(duì)濕度(RHU)和水汽壓(VP))、高程數(shù)據(jù)(DEM)、氣溶膠數(shù)據(jù)(AOD)、人口密度(POP)、植被指數(shù)(NDVI)和GNSS ZTD數(shù)據(jù)為特征變量，使用SPSS軟件的雙變量分析功能分析各個(gè)特征變量與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ù)達(dá)0.37；PRE，WS，TEM，RHU，VP，ZTD與PM2.5濃度呈顯著負(fù)相關(guān)，其中VP的相關(guān)系數(shù)達(dá)-0.74。

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

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

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

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

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

2.2 模型構(gòu)建

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

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

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

2.3 精度評(píng)估

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

為了驗(yàn)證模型的預(yù)測(cè)精度，采用平均絕對(duì)誤差(Mean Absolute Error，MAE)、均方根誤差(Root Mean Square Error，RMSE)和決定系數(shù)(R2)3個(gè)指標(biāo)來(lái)評(píng)價(jià)各個(gè)模型的預(yù)測(cè)精度，結(jié)果如表5所示。

表5 模型精度對(duì)比

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

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

(a) GWR模型

(b) GTWR模型

(c) GWRK模型

(d) GTWRK模型

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

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

(a) 2015年

(b) 2016年

(c) 2017年

(d) 2018年

(e) 2019年

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

4 結(jié)束語(yǔ)

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

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

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

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

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