基于地理加權回歸張力樣條函數的廣西PM2.5濃度插值

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

(1.桂林理工大學測繪地理信息學院，桂林 541006；2.廣西空間信息與測繪重點實驗室，桂林 541006)

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

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

1 數學模型

1.1 克里金插值

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

(1)

1.2 地理加權回歸

地理加權回歸(GWR)是一種空間分析技術。空間數據一般具有空間非平穩性的特征，用一般線性回歸模型來擬合空間數據，其分析結果不能全面反映空間數據的真實特征。GWR是一種相對簡單而又有效的探測空間非平穩性的方法，屬于局域空間分析模型。它允許不同的地理空間存在不同的空間關系，其結果是局域而不是全域的參數估計，因此能夠探測到空間數據的空間非平穩性。GWR通過建立空間范圍內每個點處的局部回歸方程，來探索研究對象在某一尺度下的空間變化及相關驅動因素。由于它考慮到了空間對象的局部效應，因此其優勢是具有較好的準確性，其基本原理[16]可表示為

(2)

式(2)中：(ui,vi)為第i個采樣點坐標；βk(ui,vi)為第i個采樣點上的第k個回歸參數；xik為第i個觀測點的第k個影響變量；p為影響變量個數；εi為回歸殘差。

1.3 地理加權回歸張力樣條函數插值

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

(3)

式(3)中：S(x,y)為插值結果；a為趨勢函數；N為插值區域點的個數；λj為通過求解線性方程組獲得的系數；rj為點(x,y)到第j個點的距離；φ為權重參數；k0()為修正貝塞爾函數；c為常數，c≈0.577 215。

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

GWR-TSFPM2.5=GWRPM2.5+TSFGWRr

(4)

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

2 實驗數據

以廣西地區2018年的氣象和空氣質量參數年均值為實驗數據，廣西地區氣象監測站點為18個，采集其站點2018年年均氣溫、風速、氣壓和降水量等氣象數據，數據來源為中國氣象科學數據共享服務網，空氣質量監測站點為49個，選擇其中7個站點作為模型驗證集，42個站點為模型訓練集，采集其站點2018年年均PM2.5、CO、SO2、NO2和O3等大氣因子濃度數據，數據來源為環境專業知識服務系統http://envi.ckcest.cn/environment/，站點分布如圖1所示。

圖1 廣西氣象站與空氣質量監測站分布

3 實驗結果分析與討論

3.1 克里金模型

以廣西地區空氣質量監測訓練集42個站點2018年PM2.5濃度年均值數據為基礎，依據式(1)建立廣西地區克里金模型，得到的克里金插值結果如圖2所示。

圖2 克里金法PM2.5插值結果

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

3.2 地理加權回歸模型

因為GWR模型為回歸模型，具有參考多個解釋變量建模的優勢，由于空氣質量監測站點只能得到各類大氣污染物含量數據，并不能有效獲取其站點位置的氣象數據，而廣西地區的氣象監測站數量較少，因此為了有效得到各個空氣質量監測站點的氣象數據，采用反距離加權插值法對氣象站點氣溫、風速、氣壓以及降水量進行空間插值，并將其對應的插值結果提取至各個空氣質量監測站點，各氣象參數插值結果如圖3所示。

圖3 氣象參數反距離加權插值結果

在得到各氣象參數插值結果后，將其值提取到各個空氣質量監測站點，得到各站點的氣象參數數據。由于GWR模型不能對具有多重共線性的變量進行建模，所以找到合適的變量組合是完成模型的前提條件和關鍵因素，因此在建模之前還應當對空氣質量監測站點的各類數據進行共線性診斷，得到的結果如表1所示。

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

表1 變量共線性診斷結果

從表2中數據可以看出，所有維度條件指數均小于30，方差比例也沒有出現某一維度具有多個高于50%的變量，說明各變量間不存在強多重共線性，可以用于模型構建。確定變量組合后，以廣西地區42個空氣質量監測站點作為訓練集，7個空氣質量監測站點作為驗證集建立GWR模型，模型解釋變量為CO、SO2、NO2和風速，得到的GWR驗證集結果如圖4所示，其模型殘差結果如圖5所示。

表2 CO、SO2、NO2、風速共線性診斷結果

圖4 GWR驗證集結果

圖5 GWR殘差

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

3.3 地理加權回歸張力樣條模型

在完成GWR模型對該地區PM2.5濃度估算以后，對GWR殘差進行空間自相關分析，其結果如表3所示。

表3中，MoranI指數為正，說明殘差呈空間正相關，其值越大則空間相關性越明顯，Z得分和P值分別表示標準差的倍數和空間分析中產生隨機事件的概率，GWR殘差Z得分為3.429，P值為0.001，則表示隨機產生此聚類模式的可能性小于1%，結果可信，可根據GWR殘差的空間自相關性運用張力樣條函數對其進行空間插值運算，結果如圖6所示。

表3 GWR殘差空間自相關分析

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

完成殘差插值計算后得到的張力樣條函數插值結果提取值至驗證集站點得到新的殘差結果，接著將GWR預測值與插值處理后的殘差值進行疊加運算便得到了GWR-TSF模型的最終插值結果，對7個驗證集站點由3種模型得到的插值結果與真實值對比，結果如表4所示。

表4 驗證集數據對比

3.4 插值結果評估

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

(5)

(6)

3種模型插值精度統計結果如表5所示。

表5 插值精度對比

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

4 結論

主要以廣西地區空氣質量監測站點和氣象站點監測數據為基礎，運用了3種模型對廣西地區進行了PM2.5濃度插值分析，證明了GWR-TSF組合模型效果相較于傳統的克里金模型和地理加權回歸模型精度更好，更適用于廣西地區PM2.5濃度插值研究，不過由于獲取的氣象監測站點較少，氣象數據插值后的結果多重共線性現象嚴重，因此對模型的精度產生了一定的影響，該組合模型仍然有很大的改進空間。