基于地理加權(quán)回歸張力樣條函數(shù)的廣西PM2.5濃度插值

謝劭峰，魏朋志，黃良珂*，張 偉，黎峻宇

(1.桂林理工大學(xué)測(cè)繪地理信息學(xué)院，桂林 541006；2.廣西空間信息與測(cè)繪重點(diǎn)實(shí)驗(yàn)室，桂林 541006)

大氣霾污染因其對(duì)人體健康、生態(tài)環(huán)境和氣候變化的影響而成為全球關(guān)注的嚴(yán)重環(huán)境問(wèn)題，PM2.5是霾污染的主要原因[1]。中外學(xué)者對(duì)PM2.5預(yù)測(cè)研究熱度也在不斷提高。周體鵬[2]基于克里金插值法對(duì)昆明市PM2.5濃度變化進(jìn)行了估算；盧月明等[3]基于局部加權(quán)線性回歸模型提出了一種引入正則化項(xiàng)的空間插值方法，該方法提高了PM2.5插值效果；謝劭峰等[4]對(duì)南寧市PM2.5濃度與氣象因素的關(guān)系進(jìn)行了探討并結(jié)合多種氣象因素運(yùn)用多元線性回歸模型對(duì)其濃度變化進(jìn)行了預(yù)測(cè)；王德冬等[5]利用時(shí)空回歸克里金法對(duì)區(qū)域PM2.5進(jìn)行了時(shí)空建模及插值；王娟[6]利用灰色關(guān)聯(lián)度、多元回歸分析等方法定性定量分析了中國(guó)30個(gè)代表城市的污染程度及污染規(guī)律；胡穩(wěn)等[7]利用普通克里金(ordinary Kriging，OK)法進(jìn)行空間插值獲取PM2.5、PM10分布特征并比較了6種半變異函數(shù)模型的適用性；Masood等[8]建立了基于機(jī)器學(xué)習(xí)方法的PM2.5預(yù)測(cè)模型；焦利民等[9]基于土地利用回歸模型進(jìn)行了武漢市PM2.5濃度高分辨率空間分布模擬研究；李爽等[10]將主成分分析與逐步多元線性回歸相結(jié)合,提出了一種改進(jìn)的土地利用回歸模型模擬大區(qū)域PM2.5濃度空間分布的方法;劉妍月等[11]運(yùn)用多種插值方法對(duì)長(zhǎng)沙市大氣中PM2.5濃度分布進(jìn)行比較研究，發(fā)現(xiàn)基于反距離加權(quán)的克里金插值方法效果較好；車?yán)诘萚12]運(yùn)用一種基于多尺度最小二乘支持向量機(jī)優(yōu)化的克里金插值方法對(duì)青島市PM2.5濃度變化進(jìn)行插值計(jì)算，效果優(yōu)于傳統(tǒng)克里金模型插值效果；趙陽(yáng)陽(yáng)等[13]運(yùn)用協(xié)同時(shí)空地理加權(quán)回歸PM2.5濃度估算法對(duì)京津冀地區(qū)進(jìn)行實(shí)驗(yàn)，實(shí)驗(yàn)效果較傳統(tǒng)單一核函數(shù)時(shí)空地理加權(quán)回歸模型有所提升；陳輝等[14]利用地理加權(quán)回歸模型進(jìn)行全國(guó)區(qū)域PM2.5遙感估算，其效果優(yōu)于多元線性回歸模型；鄧悅等[15]以北京市為實(shí)驗(yàn)區(qū)域，在地理加權(quán)回歸模型基礎(chǔ)上加入了貝葉斯先驗(yàn)信息以降低弱數(shù)據(jù)對(duì)回歸模型的影響。

上述方法在一定條件下都取得了較好的預(yù)測(cè)效果，但這些方法并沒有對(duì)模型的回歸殘差進(jìn)行很好的處理。另外PM2.5存在著很強(qiáng)的空間異質(zhì)性和空間非平穩(wěn)性，上述模型難以處理或同時(shí)處理這兩個(gè)PM2.5分布特征。現(xiàn)以廣西地區(qū)49個(gè)空氣質(zhì)量監(jiān)測(cè)站點(diǎn)和18個(gè)氣象監(jiān)測(cè)站點(diǎn)2018年監(jiān)測(cè)數(shù)據(jù)年均值為數(shù)據(jù)基礎(chǔ)，建立地理加權(quán)回歸張力樣條函數(shù)(geographically weighted regression-tension spline function, GWR-TSF)組合模型進(jìn)行PM2.5濃度插值，并與克里金和地理加權(quán)回歸模型進(jìn)行對(duì)比。

1 數(shù)學(xué)模型

1.1 克里金插值

克里金(Kriging)是一種依據(jù)協(xié)方差函數(shù)對(duì)隨機(jī)過(guò)程或者隨機(jī)場(chǎng)進(jìn)行空間建模和插值的回歸算法，在特定的隨機(jī)過(guò)程，例如固有平穩(wěn)過(guò)程中，克里金法能夠給出最優(yōu)線性無(wú)偏估計(jì)，因此在地統(tǒng)計(jì)學(xué)中也被稱為空間最優(yōu)無(wú)偏估計(jì)器，該方法不僅考慮被估點(diǎn)位置與已知數(shù)據(jù)位置的相互關(guān)系，而且還考慮已知點(diǎn)位置之間的相互聯(lián)系，因此更能反映客觀地質(zhì)規(guī)律，估值精度相對(duì)較高，該方法的適用條件為區(qū)域化變量存在空間相關(guān)性，其原理[16]可表示為

(1)

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

地理加權(quán)回歸(GWR)是一種空間分析技術(shù)。空間數(shù)據(jù)一般具有空間非平穩(wěn)性的特征，用一般線性回歸模型來(lái)擬合空間數(shù)據(jù)，其分析結(jié)果不能全面反映空間數(shù)據(jù)的真實(shí)特征。GWR是一種相對(duì)簡(jiǎn)單而又有效的探測(cè)空間非平穩(wěn)性的方法，屬于局域空間分析模型。它允許不同的地理空間存在不同的空間關(guān)系，其結(jié)果是局域而不是全域的參數(shù)估計(jì)，因此能夠探測(cè)到空間數(shù)據(jù)的空間非平穩(wěn)性。GWR通過(guò)建立空間范圍內(nèi)每個(gè)點(diǎn)處的局部回歸方程，來(lái)探索研究對(duì)象在某一尺度下的空間變化及相關(guān)驅(qū)動(dòng)因素。由于它考慮到了空間對(duì)象的局部效應(yīng)，因此其優(yōu)勢(shì)是具有較好的準(zhǔn)確性，其基本原理[16]可表示為

(2)

式(2)中：(ui,vi)為第i個(gè)采樣點(diǎn)坐標(biāo)；βk(ui,vi)為第i個(gè)采樣點(diǎn)上的第k個(gè)回歸參數(shù)；xik為第i個(gè)觀測(cè)點(diǎn)的第k個(gè)影響變量；p為影響變量個(gè)數(shù)；εi為回歸殘差。

1.3 地理加權(quán)回歸張力樣條函數(shù)插值

地理加權(quán)回歸張力樣條函數(shù)(GWR-TSF)插值是一種將地理加權(quán)回歸(GWR)與張力樣條函數(shù)(TSF)插值結(jié)合起來(lái)的綜合分析方法。GWR模型是對(duì)普通線性回歸模型的擴(kuò)展；張力樣條函數(shù)是徑向基函數(shù)插值法的一種，該方法插值速度快以及估測(cè)大小的范圍不局限，張力樣條函數(shù)的基本原理[17]可表示為

(3)

式(3)中：S(x,y)為插值結(jié)果；a為趨勢(shì)函數(shù)；N為插值區(qū)域點(diǎn)的個(gè)數(shù)；λj為通過(guò)求解線性方程組獲得的系數(shù)；rj為點(diǎn)(x,y)到第j個(gè)點(diǎn)的距離；φ為權(quán)重參數(shù)；k0()為修正貝塞爾函數(shù)；c為常數(shù)，c≈0.577 215。

GWR-TSF模型利用張力樣條函數(shù)對(duì)GWR模型得到的回歸殘差ε進(jìn)行空間插值，然后將得到的殘差插值結(jié)果與GWR回歸估計(jì)值進(jìn)行疊加，從而獲得GWR-TSF模型估算值，其原理可表示為

GWR-TSFPM2.5=GWRPM2.5+TSFGWRr

(4)

式(4)中：GWR-TSFPM2.5為GWR-TSF模型對(duì)PM2.5的濃度估算值；GWRPM2.5為GWR模型對(duì)PM2.5的濃度估算值；TSFGWRr為TSF插值法對(duì)GWR模型估算PM2.5的濃度值后產(chǎn)生的回歸殘差進(jìn)行區(qū)域插值得到的結(jié)果。

2 實(shí)驗(yàn)數(shù)據(jù)

以廣西地區(qū)2018年的氣象和空氣質(zhì)量參數(shù)年均值為實(shí)驗(yàn)數(shù)據(jù)，廣西地區(qū)氣象監(jiān)測(cè)站點(diǎn)為18個(gè)，采集其站點(diǎn)2018年年均氣溫、風(fēng)速、氣壓和降水量等氣象數(shù)據(jù)，數(shù)據(jù)來(lái)源為中國(guó)氣象科學(xué)數(shù)據(jù)共享服務(wù)網(wǎng)，空氣質(zhì)量監(jiān)測(cè)站點(diǎn)為49個(gè)，選擇其中7個(gè)站點(diǎn)作為模型驗(yàn)證集，42個(gè)站點(diǎn)為模型訓(xùn)練集，采集其站點(diǎn)2018年年均PM2.5、CO、SO2、NO2和O3等大氣因子濃度數(shù)據(jù)，數(shù)據(jù)來(lái)源為環(huán)境專業(yè)知識(shí)服務(wù)系統(tǒng)http://envi.ckcest.cn/environment/，站點(diǎn)分布如圖1所示。

圖1 廣西氣象站與空氣質(zhì)量監(jiān)測(cè)站分布

3 實(shí)驗(yàn)結(jié)果分析與討論

3.1 克里金模型

以廣西地區(qū)空氣質(zhì)量監(jiān)測(cè)訓(xùn)練集42個(gè)站點(diǎn)2018年P(guān)M2.5濃度年均值數(shù)據(jù)為基礎(chǔ)，依據(jù)式(1)建立廣西地區(qū)克里金模型，得到的克里金插值結(jié)果如圖2所示。

圖2 克里金法PM2.5插值結(jié)果

從圖2可以看出，廣西地區(qū)的PM2.5濃度分布情況主要為東高西低，最嚴(yán)重的區(qū)域主要集中在柳州市和來(lái)賓市一帶。

3.2 地理加權(quán)回歸模型

因?yàn)镚WR模型為回歸模型，具有參考多個(gè)解釋變量建模的優(yōu)勢(shì)，由于空氣質(zhì)量監(jiān)測(cè)站點(diǎn)只能得到各類大氣污染物含量數(shù)據(jù)，并不能有效獲取其站點(diǎn)位置的氣象數(shù)據(jù)，而廣西地區(qū)的氣象監(jiān)測(cè)站數(shù)量較少，因此為了有效得到各個(gè)空氣質(zhì)量監(jiān)測(cè)站點(diǎn)的氣象數(shù)據(jù)，采用反距離加權(quán)插值法對(duì)氣象站點(diǎn)氣溫、風(fēng)速、氣壓以及降水量進(jìn)行空間插值，并將其對(duì)應(yīng)的插值結(jié)果提取至各個(gè)空氣質(zhì)量監(jiān)測(cè)站點(diǎn)，各氣象參數(shù)插值結(jié)果如圖3所示。

圖3 氣象參數(shù)反距離加權(quán)插值結(jié)果

在得到各氣象參數(shù)插值結(jié)果后，將其值提取到各個(gè)空氣質(zhì)量監(jiān)測(cè)站點(diǎn)，得到各站點(diǎn)的氣象參數(shù)數(shù)據(jù)。由于GWR模型不能對(duì)具有多重共線性的變量進(jìn)行建模，所以找到合適的變量組合是完成模型的前提條件和關(guān)鍵因素，因此在建模之前還應(yīng)當(dāng)對(duì)空氣質(zhì)量監(jiān)測(cè)站點(diǎn)的各類數(shù)據(jù)進(jìn)行共線性診斷，得到的結(jié)果如表1所示。

從表1中方差比例來(lái)看，除第7維度中有兩個(gè)變量(CO和NO2)方差比例同時(shí)高于50%以外，其余維度中均最多只有一個(gè)變量方差比例高于50%，因此不具備多變量存在多重共線性的判斷條件；而從條件指數(shù)來(lái)看，條件指數(shù)為最大的主成分與當(dāng)前主成分比值的平方根，從第6維度到第9維度的條件指數(shù)均大于30，即存在多個(gè)維度條件指數(shù)大于30，證明此處用于建模所選的8類變量之間存在著多重共線性，而多重共線性是指回歸模型中的解釋變量之間由于存在精確相關(guān)關(guān)系或高度相關(guān)關(guān)系而使模型估計(jì)失真或難以估計(jì)準(zhǔn)確，因此需要剔除掉一些多余的變量才能進(jìn)行建模實(shí)驗(yàn)，經(jīng)過(guò)反復(fù)比較實(shí)驗(yàn)，最后得到的變量組合為CO、SO2、NO2和風(fēng)速，變量組合共線性診斷結(jié)果如表2所示。

表1 變量共線性診斷結(jié)果

從表2中數(shù)據(jù)可以看出，所有維度條件指數(shù)均小于30，方差比例也沒有出現(xiàn)某一維度具有多個(gè)高于50%的變量，說(shuō)明各變量間不存在強(qiáng)多重共線性，可以用于模型構(gòu)建。確定變量組合后，以廣西地區(qū)42個(gè)空氣質(zhì)量監(jiān)測(cè)站點(diǎn)作為訓(xùn)練集，7個(gè)空氣質(zhì)量監(jiān)測(cè)站點(diǎn)作為驗(yàn)證集建立GWR模型，模型解釋變量為CO、SO2、NO2和風(fēng)速，得到的GWR驗(yàn)證集結(jié)果如圖4所示，其模型殘差結(jié)果如圖5所示。

表2 CO、SO2、NO2、風(fēng)速共線性診斷結(jié)果

圖4 GWR驗(yàn)證集結(jié)果

圖5 GWR殘差

從圖4可以看出，基于GWR模型得到的驗(yàn)證集站點(diǎn)PM2.5濃度的大小與分布規(guī)律同克里金插值法所得到的插值結(jié)果大致相同，以柳州市和貴港市一帶區(qū)域數(shù)值最高。由圖5的殘差結(jié)果分析可知，殘差絕對(duì)值較大的區(qū)域也主要分布于廣西中東部地區(qū)。

3.3 地理加權(quán)回歸張力樣條模型

在完成GWR模型對(duì)該地區(qū)PM2.5濃度估算以后，對(duì)GWR殘差進(jìn)行空間自相關(guān)分析，其結(jié)果如表3所示。

表3中，MoranI指數(shù)為正，說(shuō)明殘差呈空間正相關(guān)，其值越大則空間相關(guān)性越明顯，Z得分和P值分別表示標(biāo)準(zhǔn)差的倍數(shù)和空間分析中產(chǎn)生隨機(jī)事件的概率，GWR殘差Z得分為3.429，P值為0.001，則表示隨機(jī)產(chǎn)生此聚類模式的可能性小于1%，結(jié)果可信，可根據(jù)GWR殘差的空間自相關(guān)性運(yùn)用張力樣條函數(shù)對(duì)其進(jìn)行空間插值運(yùn)算，結(jié)果如圖6所示。

表3 GWR殘差空間自相關(guān)分析

從圖6中對(duì)于GWR殘差的張力樣條函數(shù)插值結(jié)果可知，廣西地區(qū)GWR殘差值的分布規(guī)律為中東部區(qū)域數(shù)值為正，其絕對(duì)值較大的區(qū)域主要位于柳州、來(lái)賓、貴港、梧州、賀州和桂林六市；而廣西南部區(qū)域殘差多呈負(fù)值，其絕對(duì)值較大區(qū)域以南寧、崇左、防城港、欽州、北海和玉林六市為主。

完成殘差插值計(jì)算后得到的張力樣條函數(shù)插值結(jié)果提取值至驗(yàn)證集站點(diǎn)得到新的殘差結(jié)果，接著將GWR預(yù)測(cè)值與插值處理后的殘差值進(jìn)行疊加運(yùn)算便得到了GWR-TSF模型的最終插值結(jié)果，對(duì)7個(gè)驗(yàn)證集站點(diǎn)由3種模型得到的插值結(jié)果與真實(shí)值對(duì)比，結(jié)果如表4所示。

表4 驗(yàn)證集數(shù)據(jù)對(duì)比

3.4 插值結(jié)果評(píng)估

為了更直觀地評(píng)價(jià)各模型的性能，采用均方根誤差(root mean square error, RMSE)以及平均絕對(duì)誤差(mean absolute error, MAE)對(duì)插值結(jié)果進(jìn)行精度評(píng)定，相應(yīng)的計(jì)算公式為

(5)

(6)

3種模型插值精度統(tǒng)計(jì)結(jié)果如表5所示。

表5 插值精度對(duì)比

從表5可以看出，GWR-TSF組合模型插值精度最高，GWR-TSF組合模型均方根誤差的值較前兩種模型分別提高了20.68%和25.71%，而平均絕對(duì)誤差的值也分別提高了20.22%和11.62%，提升幅度基本都在20%左右，說(shuō)明該組合模型在區(qū)域性插值PM2.5這一類空間異質(zhì)性較強(qiáng)的變化因素時(shí)，效果要優(yōu)于傳統(tǒng)的單一模型。

4 結(jié)論

主要以廣西地區(qū)空氣質(zhì)量監(jiān)測(cè)站點(diǎn)和氣象站點(diǎn)監(jiān)測(cè)數(shù)據(jù)為基礎(chǔ)，運(yùn)用了3種模型對(duì)廣西地區(qū)進(jìn)行了PM2.5濃度插值分析，證明了GWR-TSF組合模型效果相較于傳統(tǒng)的克里金模型和地理加權(quán)回歸模型精度更好，更適用于廣西地區(qū)PM2.5濃度插值研究，不過(guò)由于獲取的氣象監(jiān)測(cè)站點(diǎn)較少，氣象數(shù)據(jù)插值后的結(jié)果多重共線性現(xiàn)象嚴(yán)重，因此對(duì)模型的精度產(chǎn)生了一定的影響，該組合模型仍然有很大的改進(jìn)空間。