于志英，張福浩，仇阿根，趙陽陽

(中國測(cè)繪科學(xué)研究院，北京 100830)

地理加權(quán)回歸(geographically weighted regression,GWR)是地學(xué)領(lǐng)域常用的空間分析方法，被廣泛用于空間非平穩(wěn)性探測(cè)。其基本思想是將采樣點(diǎn)的空間位置嵌入回歸參數(shù)中，利用局部加權(quán)最小二乘方法逐點(diǎn)進(jìn)行參數(shù)估計(jì)。通過構(gòu)建回歸模型，探究空間模式背后的影響因素或預(yù)測(cè)空間現(xiàn)象的發(fā)展情況。

諸多學(xué)者對(duì)地理加權(quán)回歸模型進(jìn)行了深入研究。文獻(xiàn)[1]提出了混合地理加權(quán)回歸模型，將回歸模型中隨空間位置變化的參數(shù)稱為局域參數(shù)，不受空間位置影響的參數(shù)稱為全域參數(shù)，并采用迭代的方法對(duì)參數(shù)進(jìn)行近似估計(jì)。文獻(xiàn)[2]在此基礎(chǔ)上，對(duì)混合地理加權(quán)回歸模型進(jìn)行推導(dǎo)，得到了全域參數(shù)和局域參數(shù)的精確表達(dá)。文獻(xiàn)[3]將時(shí)間要素融入地理加權(quán)回歸模型，提出了時(shí)空地理加權(quán)回歸。文獻(xiàn)[4]對(duì)時(shí)空地理加權(quán)回歸模型估計(jì)方法、核函數(shù)選擇、因子選擇、多重共線性檢驗(yàn)和參數(shù)估計(jì)過程進(jìn)行了詳細(xì)介紹。文獻(xiàn)[5]提出了基于半監(jiān)督的地理加權(quán)回歸方法，利用有標(biāo)記樣本訓(xùn)練無標(biāo)記樣本，選擇置信度高的結(jié)果擴(kuò)充有標(biāo)記樣本，解決了地理加權(quán)回歸樣本量較少情況下模型精度不高的問題。文獻(xiàn)[6]針對(duì)現(xiàn)有模型無法充分?jǐn)M合復(fù)雜非線性關(guān)系的問題，提出了地理時(shí)空神經(jīng)網(wǎng)絡(luò)加權(quán)回歸模型，有效擬合時(shí)間鄰近和空間鄰近的非線性融合作用。另有學(xué)者從中間過程出發(fā)，對(duì)地理加權(quán)回歸模型進(jìn)行改進(jìn)。文獻(xiàn)[7]將非線性主成分用于回歸模型變量選擇，既避免了變量間的共線性，又保留了原始影響因素的主要信息。文獻(xiàn)[8]用絕對(duì)值交叉驗(yàn)證(absolute value cross-validation，ACV)替代交叉驗(yàn)證(cross-validation,CV)[9]進(jìn)行了最優(yōu)帶寬選取，避免了CV法中二次誤差標(biāo)準(zhǔn)放大離群值的影響。

通過以上研究現(xiàn)狀發(fā)現(xiàn)，各方法從不同角度對(duì)地理加權(quán)回歸模型進(jìn)行改進(jìn)，但缺乏異常值檢測(cè)與處理過程。因此，本文提出基于IGGⅢ的地理加權(quán)回歸模型，將IGGⅢ權(quán)函數(shù)用于地理加權(quán)回歸參數(shù)估計(jì)過程，對(duì)鄰域觀測(cè)值進(jìn)行降權(quán)或剔除。首先對(duì)方法原理進(jìn)行詳細(xì)介紹，并介紹方法流程，最后將模擬數(shù)據(jù)和真實(shí)數(shù)據(jù)與GWR、ACV-GWR進(jìn)行對(duì)比試驗(yàn)，利用均方誤差(mean square error，MSE)、平均絕對(duì)誤差(mean absolute error，MAE)和R2作為指標(biāo)進(jìn)行評(píng)價(jià)。

1 基于IGGⅢ的地理加權(quán)回歸原理

1.1 方法原理

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

(1)

式中，y為因變量；x為自變量；(u,v)為采樣點(diǎn)位置；βk(ui,vi)為第i個(gè)采樣點(diǎn)的第k個(gè)參數(shù)，β的取值與采樣點(diǎn)位置有關(guān)；ε為隨機(jī)誤差，符合正態(tài)分布。

未知參數(shù)β的估計(jì)采用最小二乘方法實(shí)現(xiàn)。當(dāng)數(shù)據(jù)中存在離群值時(shí)，離群值參與回歸點(diǎn)的參數(shù)估計(jì)過程，其二次誤差會(huì)主導(dǎo)殘差平方和的值，影響參數(shù)估計(jì)結(jié)果。為削弱鄰域點(diǎn)中的離群值對(duì)回歸點(diǎn)參數(shù)估計(jì)的影響，本文在地理加權(quán)回歸模型中引入基于權(quán)函數(shù)的粗差處理方法。常用的權(quán)函數(shù)有Huber法[10]、Hample法、Turkey法、Danish法、IGG方案[11]、IGGⅢ方案等。考慮Turkey權(quán)函數(shù)為有界連續(xù)函數(shù)；Danish法實(shí)質(zhì)上為淘汰法，沒有抗差上的論證[11]；IGG法為有淘汰區(qū)的M估計(jì)，權(quán)因子變化平緩[12]，性能優(yōu)于Huber法、Hample法[12-13]；IGG為跳躍函數(shù)，IGGⅢ權(quán)函數(shù)為連續(xù)函數(shù)。因此，本文選擇IGGⅢ方案[14]中的權(quán)函數(shù)用于加權(quán)最小二乘參數(shù)估計(jì)過程。該函數(shù)采用三段法進(jìn)行權(quán)重定義，對(duì)正常段的觀測(cè)采用最小二乘估計(jì)，對(duì)可用觀測(cè)采用權(quán)因子降權(quán)，權(quán)因子取0～1之間的變值，對(duì)達(dá)到淘汰界的離群值進(jìn)行剔除。權(quán)函數(shù)如下所示

(2)

將其繪制成圖直觀展示，如圖1所示。

將空間距離權(quán)重和觀測(cè)點(diǎn)可靠性權(quán)重同時(shí)納入地理加權(quán)回歸參數(shù)估計(jì)模型，模型表達(dá)為

(3)

式中，wij為空間距離權(quán)重，在地理加權(quán)回歸中常采用Gauss和Bi-square兩種核函數(shù)進(jìn)行計(jì)算。由于不同核函數(shù)對(duì)模型參數(shù)估計(jì)的影響相差不大，最優(yōu)帶寬的選取對(duì)參數(shù)估計(jì)結(jié)果影響較大[15]，因此，選擇最優(yōu)帶寬值以確定合適的空間權(quán)重矩陣對(duì)模型解算至關(guān)重要。常用的帶寬選取方法有Akaike信息量準(zhǔn)則(Akaike information criterion，AIC)、貝葉斯信息準(zhǔn)則(Bayesian information criterion,BIC)及交叉驗(yàn)證等方法，AIC準(zhǔn)則和BIC準(zhǔn)則通過極大似然估計(jì)計(jì)算，CV法采用二次誤差標(biāo)準(zhǔn)進(jìn)行計(jì)算。

1.2 方法流程

方法流程如圖2所示，流程說明：

(1) 構(gòu)建地理加權(quán)回歸模型。利用相關(guān)性分析和共線性分析選擇變量，并構(gòu)建地理加權(quán)回歸模型。

(2) 計(jì)算最優(yōu)帶寬。利用AIC準(zhǔn)則、BIC準(zhǔn)則、CV或ACV進(jìn)行最優(yōu)帶寬選取。

(3) 計(jì)算空間權(quán)重矩陣。根據(jù)最優(yōu)帶寬，利用Gauss核函數(shù)、Bi-square核函數(shù)或自適應(yīng)核函數(shù)計(jì)算空間權(quán)重矩陣。

(4) 估計(jì)回歸模型參數(shù)。構(gòu)建地理加權(quán)回歸參數(shù)估計(jì)模型，根據(jù)最小化損失函數(shù)原則進(jìn)行參數(shù)求解。

(5) 計(jì)算殘差。依據(jù)步驟(4)計(jì)算所得的參數(shù)值結(jié)果，計(jì)算因變量y的估計(jì)值。因變量觀測(cè)值與估計(jì)值的差值即為殘差。

(6) 計(jì)算評(píng)價(jià)指標(biāo)。計(jì)算標(biāo)準(zhǔn)化殘差及MSE、MAE、R2等評(píng)價(jià)指標(biāo)。

(7) 判斷是否存在離群值。通過標(biāo)準(zhǔn)化殘差判斷模型中是否存在離群值。若存在，利用IGGⅢ計(jì)算權(quán)因子，構(gòu)建基于IGGⅢ的地理加權(quán)回歸模型，迭代進(jìn)行模型求解，直到模型中不存在離群值時(shí)，結(jié)束迭代；若不存在，輸出結(jié)果。

2 試 驗(yàn)

2.1 模擬數(shù)據(jù)試驗(yàn)

2.1.1 數(shù)據(jù)生成

本文根據(jù)地理加權(quán)回歸模型特性設(shè)計(jì)模擬數(shù)據(jù)。其中，自變量和因變量滿足線性回歸關(guān)系，系數(shù)與采樣點(diǎn)空間位置有關(guān)，具體公式見表1。

表1 模擬數(shù)據(jù)生成公式

表1中，x1、x2為自變量，服從(0,1)均勻分布；u、v為位置變量，服從[0,20]均勻分布；ε為隨機(jī)誤差。此外，向模擬數(shù)據(jù)中添加高斯白噪聲。

2.1.2 對(duì)比試驗(yàn)設(shè)置

本文將GWR、ACV-GWR和IGGⅢ-GWR進(jìn)行對(duì)比試驗(yàn)，采用Gauss核函數(shù)計(jì)算空間權(quán)重矩陣，GWR和IGGⅢ-GWR采用CV確定最優(yōu)帶寬，ACV-GWR采用ACV確定最優(yōu)帶寬。ACV計(jì)算方法為

(4)

2.1.3 試驗(yàn)結(jié)果分析

利用模擬數(shù)據(jù)對(duì)以上3種方法分別試驗(yàn)40次，表2列舉部分試驗(yàn)結(jié)果。表3展示IGGⅢ-GWR較ACV-GWR、GWR各指標(biāo)性能平均提升情況。

表2 部分試驗(yàn)結(jié)果

表3 各指標(biāo)平均提升情況

從MSE、MAE、R2性能提升百分比來看，IGGⅢ-GWR比GWR性能分別提升51.14%、23.77%、28.4%，比ACV-GWR分別提升49.96%、22.57%、27.1%。

2.2 真實(shí)數(shù)據(jù)試驗(yàn)

2.2.1 試驗(yàn)數(shù)據(jù)

本文選用2016年1月至2018年3月北京地區(qū)空氣質(zhì)量及其影響因素作為試驗(yàn)數(shù)據(jù)進(jìn)行分析驗(yàn)證。計(jì)算CO、NO2、O3、PM10、PM2.5、SO2與空氣質(zhì)量指數(shù)的Pearson相關(guān)系數(shù)(見表4)，CO、NO2、PM2.5、SO2與空氣質(zhì)量指數(shù)顯著正相關(guān)，O3與空氣質(zhì)量指數(shù)呈較強(qiáng)的負(fù)相關(guān)，PM10與空氣質(zhì)量指數(shù)呈較強(qiáng)的正相關(guān)。

表4 各影響因素與AQI間的Pearson相關(guān)系數(shù)

對(duì)包含CO、NO2、O3、PM10、PM2.5、SO2在內(nèi)的6種污染物進(jìn)行多重共線性分析，以方差膨脹因子小于2且條件索引小于10為限定條件進(jìn)行因子選取，最終選擇O3、PM2.5、SO2作為影響指標(biāo)構(gòu)建回歸模型。北京地區(qū)空氣質(zhì)量監(jiān)測(cè)站點(diǎn)分布如圖3所示。

2.2.2 試驗(yàn)結(jié)果分析

分別采用GWR、ACV-GWR、IGGⅢ-GWR進(jìn)行試驗(yàn)，以MSE、MAE和R2作為指標(biāo)對(duì)試驗(yàn)結(jié)果進(jìn)行評(píng)價(jià),見表5。

表5 真實(shí)數(shù)據(jù)試驗(yàn)結(jié)果

從MSE、MAE、R2指標(biāo)性能上看，IGGⅢ-GWR較GWR分別提升12.65%、7.44%、0.37%，較ACV-GWR分別提升11.85%、6.96%、0.34%。

隨機(jī)選取任意月份繪制空氣質(zhì)量分布和標(biāo)準(zhǔn)化殘差分布，比較各模型結(jié)果，本文以2017年1月為例進(jìn)行結(jié)果展示，如圖4—圖7所示。

圖4—圖7展示了空氣質(zhì)量觀測(cè)值和不同模型估計(jì)結(jié)果。從圖4—圖7可以看出，2017年1月北京市空氣質(zhì)量指數(shù)介于88～202之間，GWR估算結(jié)果介于29～150之間，ACV-GWR估算結(jié)果介于39～182之間，IGGⅢ-GWR估計(jì)結(jié)果介于87～208之間。從估算結(jié)果上看，IGGⅢ-GWR模型估算得到的空氣質(zhì)量指數(shù)更符合真實(shí)情況。

從空氣狀況空間分布來看，2017年1月北京市南部地區(qū)空氣質(zhì)量相對(duì)較差，GWR估計(jì)結(jié)果顯示中部地區(qū)空氣質(zhì)量相對(duì)較差，ACV-GWR估計(jì)結(jié)果規(guī)律性不明顯，IGGⅢ-GWR估計(jì)得到的空氣質(zhì)量空間分布情況與觀測(cè)值更吻合。

繪制不同回歸模型計(jì)算所得的標(biāo)準(zhǔn)化殘差分布，如圖8—圖10所示。

從圖8—圖10可以看出，GWR標(biāo)準(zhǔn)化殘差計(jì)算結(jié)果介于0～2.05之間，ACV-GWR標(biāo)準(zhǔn)化從殘差計(jì)算結(jié)果介于0～2.04之間，IGGⅢ-GWR標(biāo)準(zhǔn)化殘差計(jì)算結(jié)果介于0～1.98之間。從整體上講，IGGⅢ-GWR估計(jì)效果更好一些。從局部來看，IGGⅢ-GWR對(duì)北京中南部地區(qū)空氣質(zhì)量估計(jì)效果優(yōu)于其他兩種方法。

3 結(jié) 語

本文提出了基于IGGⅢ的地理加權(quán)回歸模型，將IGGⅢ方案應(yīng)用于地理加權(quán)回歸，降低了離群值對(duì)參數(shù)估計(jì)的影響，提高了地理加權(quán)回歸模型對(duì)離群值的抵抗能力。通過模擬數(shù)據(jù)和真實(shí)數(shù)據(jù)與GWR、ACV-GWR進(jìn)行對(duì)比試驗(yàn)，以MSE、MAE、R2作為評(píng)價(jià)指標(biāo)進(jìn)行驗(yàn)證。試驗(yàn)結(jié)果表明，IGGⅢ-GWR可用于空間非平穩(wěn)性表達(dá)與未知量預(yù)測(cè)，當(dāng)數(shù)據(jù)中存在離群值時(shí)，基于IGGⅢ的地理加權(quán)回歸模型擬合效果更好。