劉 寧，鄒 濱，張鴻輝

1.中南大學(xué)地球科學(xué)與信息物理學(xué)院，湖南 長(zhǎng)沙 410083; 2.廣東國(guó)地規(guī)劃科技股份有限公司，廣東 廣州 510075; 3.湖南師范大學(xué)資源與環(huán)境科學(xué)學(xué)院，湖南 長(zhǎng)沙 410012

地理加權(quán)回歸(geographically weighted regression,GWR)是通過(guò)假設(shè)回歸點(diǎn)僅與空間一定范圍觀測(cè)點(diǎn)有關(guān)而建立起的一種局部加權(quán)最小二乘模型[1]。建模過(guò)程中，首先基于修正赤池信息量準(zhǔn)則或留一交叉驗(yàn)證準(zhǔn)則搜索空間帶寬[2-3]，確定影響回歸點(diǎn)的臨近觀測(cè)點(diǎn)；然后構(gòu)建距離衰減權(quán)重函數(shù)對(duì)臨近觀測(cè)點(diǎn)定權(quán)后，通過(guò)局部加權(quán)最小二乘估算回歸點(diǎn)系數(shù)值與預(yù)測(cè)值[1]。該模型在探索自然地理及社會(huì)人文現(xiàn)象的空間分布、空間非平穩(wěn)性及空間尺度轉(zhuǎn)換等多個(gè)領(lǐng)域有著廣泛的應(yīng)用。文獻(xiàn)[4]采用GWR模型估算了全國(guó)PM2.5濃度空間分布。文獻(xiàn)[5]基于GWR模型揭示了上海房?jī)r(jià)的空間分異特征及其顯著影響因素。文獻(xiàn)[6]通過(guò)GWR模型將遙感MODIS降水量產(chǎn)品從990 m分辨率提升至90 m精細(xì)尺度。

然而受空間樣本稀疏與局部共線(xiàn)性等因素影響，GWR建模結(jié)果不確定性往往呈現(xiàn)空間異質(zhì)性[7]。文獻(xiàn)[8]在住宅地價(jià)估算的修正研究中發(fā)現(xiàn)，回歸樣本點(diǎn)密度降低時(shí)，GWR模型估算精度不可避免隨之下降。文獻(xiàn)[9—11]發(fā)現(xiàn)，當(dāng)局部區(qū)域兩個(gè)或多個(gè)自變量存在共線(xiàn)性問(wèn)題時(shí)，GWR模型回歸預(yù)測(cè)結(jié)果常出現(xiàn)奇異值，回歸系數(shù)符號(hào)甚至偏離預(yù)期。

針對(duì)上述問(wèn)題，文獻(xiàn)[12]嘗試將主成分分析(principle component analysis，PCA)方法引入GWR建模過(guò)程，即預(yù)先用PCA方法對(duì)GWR自變量進(jìn)行降維轉(zhuǎn)換，通過(guò)僅保留自變量主要信息的方式弱化其相關(guān)性達(dá)到解決局部共線(xiàn)性問(wèn)題的目的，但轉(zhuǎn)換后的自變量常失去原有物理含義。同時(shí)也有學(xué)者通過(guò)對(duì)自變量單獨(dú)構(gòu)建空間帶寬[13-14]、引入嶺回歸思路[15]減弱局部共線(xiàn)性影響，或引入條件數(shù)[16]、方差膨脹系數(shù)[17]等多個(gè)指標(biāo)甄別共線(xiàn)性問(wèn)題，但這些方法與指標(biāo)仍都無(wú)法表征由樣本稀疏造成的GWR模型估算精度衰減。

對(duì)此，本文擬在依據(jù)協(xié)方差傳播定律構(gòu)建后驗(yàn)標(biāo)準(zhǔn)差評(píng)估指標(biāo)基礎(chǔ)上[18-19]，提出一種地理加權(quán)回歸建模結(jié)果不確定性度量與約束方法，并基于衛(wèi)星遙感估算地表PM2.5濃度實(shí)例與傳統(tǒng)條件數(shù)約束方法開(kāi)展對(duì)比，驗(yàn)證該方法的可靠性。

1 后驗(yàn)標(biāo)準(zhǔn)差約束下的GWR模型

1.1 GWR模型原理

空間位置為(μi,vi)的GWR局部回歸模型[1]為

(1)

由式(1)可知，系數(shù)βj(μi,vi)與β0(μi,vi)為空間位置的函數(shù)。位置(μi,vi)的回歸系數(shù)列向量β(μi,vi)=[β0(μi,vi);β1(μi,vi);…;βm(μi,vi)]由局部最小加權(quán)二乘估計(jì)

(2)

式中，y=[y1y2…yn]為n個(gè)點(diǎn)位下的觀測(cè)列矢量；X為自變量矩陣為

(3)

式中，W(μi,vi)為點(diǎn)位(μi,vi)的對(duì)角權(quán)陣，其對(duì)角陣元素為[w1,iw2,i…wn,i]，wk,i代表了k點(diǎn)位對(duì)i點(diǎn)位影響的權(quán)重值，權(quán)重由隨距離衰減的核函數(shù)構(gòu)造，常用核函數(shù)有Guassian與Bisquare等[1,20]

(4)

回歸系數(shù)β(μi,vi)由式(2)估計(jì)得出后，點(diǎn)位(μi,vi)的GWR估計(jì)值為

(5)

式中，xi為自變量矩陣第i行向量。

1.2 后驗(yàn)標(biāo)準(zhǔn)差約束下GWR優(yōu)化

設(shè)Ri=(XTW(μi,vi)X)-1XTW(μi,vi)、Si=xiRi，則式(2)與式(5)可重寫(xiě)為

(6)

(7)

(8)

(9)

圖1 GWR估算值在均值為標(biāo)準(zhǔn)差分別為σ1、σ2、σ3時(shí)的正態(tài)分布曲線(xiàn)

本文擬通過(guò)剔除后驗(yàn)標(biāo)準(zhǔn)差較大的GWR估算值提升GWR模型的總體精度。具體而言，①基于式(5)與式(7)計(jì)算每個(gè)樣本點(diǎn)的十折交叉驗(yàn)證值[22]及其后驗(yàn)SD值；②設(shè)置SD閾值，剔除后驗(yàn)SD值大于閾值的樣本點(diǎn)；③比較剩余樣本的十折交叉驗(yàn)證值與觀測(cè)值，評(píng)估十折交叉驗(yàn)證精度；④精度不滿(mǎn)足要求則降低SD閾值重復(fù)②—④步，滿(mǎn)足要求則以當(dāng)前SD閾值約束GWR估算值。算法流程圖如圖2所示。

圖2 后驗(yàn)標(biāo)準(zhǔn)差約束下的GWR模型建模流程

2 試驗(yàn)設(shè)計(jì)與結(jié)果分析

2.1 數(shù)據(jù)源與預(yù)處理

本文以2019年基于地理加權(quán)回歸模型的中國(guó)區(qū)域近地表PM2.5濃度遙感制圖試驗(yàn)進(jìn)行驗(yàn)證。模型因變量為來(lái)自中國(guó)環(huán)境監(jiān)測(cè)總站發(fā)布的1356個(gè)地面站點(diǎn)PM2.5小時(shí)濃度數(shù)據(jù)(http:∥www.cnemc.cn/)。模型自變量包括遙感氣溶膠光學(xué)厚度(aerosol optical depth，AOD)、氣象數(shù)據(jù)及增強(qiáng)型植被指數(shù)(enhanced vegetation index，EVI)。遙感AOD與EVI數(shù)據(jù)均為NASA Terra/Aqua MODIS傳感器(https:∥ladsweb.modaps.eosdis.nasa.gov/)所反演的一天兩次的1 km產(chǎn)品，其兩次的過(guò)境時(shí)刻分別為當(dāng)?shù)貢r(shí)間10:30 a.m.與13:30 p.m.左右。ERA5氣象數(shù)據(jù)來(lái)自第五代ECMWF大氣再分析全球氣候產(chǎn)品(https:∥cds.climate.copernicus.eu/)，選取的氣象變量分別為邊界層高度、相對(duì)濕度、地表氣壓、2 m氣溫、經(jīng)向風(fēng)速、緯向風(fēng)速。ERA5氣象數(shù)據(jù)的空間分辨率為0.1°～0.25°，時(shí)間分辨率為1 h。

采用最鄰近插值以及樣條函數(shù)插值方法將EVI和氣象等要素值降尺度至遙感AOD網(wǎng)格，再將EVI、氣象、AOD等網(wǎng)格化自變量值與地面站點(diǎn)PM2.5濃度值進(jìn)行時(shí)空匹配。剔除因云雨天氣及積雪覆蓋等導(dǎo)致的無(wú)效值后，共計(jì)生成236 206個(gè)樣本數(shù)據(jù)。此外，為使GWR模型估算的回歸系數(shù)不受不同變量之間量級(jí)差異的影響，在GWR建模前對(duì)所有自變量與因變量均基于z-score方法開(kāi)展了標(biāo)準(zhǔn)化處理。表1為因變量地面站點(diǎn)PM2.5濃度數(shù)據(jù)與自變量AOD/EVI等數(shù)據(jù)的統(tǒng)計(jì)信息。在此基礎(chǔ)上，本文在365 d Terra/Aqua過(guò)境時(shí)刻共建立了730個(gè)模型。通過(guò)建立的GWR模型，估算近地表PM2.5遙感制圖結(jié)果及其不確定性，最后再將近地表PM2.5遙感制圖結(jié)果反標(biāo)準(zhǔn)化為濃度值。

表1 近地表日均PM2.5濃度遙感制圖樣本數(shù)據(jù)統(tǒng)計(jì)信息

2.2 試驗(yàn)設(shè)計(jì)

為充分檢驗(yàn)后驗(yàn)標(biāo)準(zhǔn)差約束下GWR模型的有效性，選擇不同參數(shù)構(gòu)建4類(lèi)GWR模型，全面探究不確定性約束前后GWR模型精度表現(xiàn)。模型A，自適應(yīng)類(lèi)型帶寬、Bisquare核函數(shù)及CV帶寬擇優(yōu)指標(biāo)；模型B，自適應(yīng)類(lèi)型帶寬、Bisquare核函數(shù)及AICc帶寬擇優(yōu)指標(biāo)；模型C，固定類(lèi)型帶寬、Guassian核函數(shù)及CV帶寬擇優(yōu)指標(biāo)；模型D，固定類(lèi)型帶寬、Guassian核函數(shù)及AICc帶寬擇優(yōu)指標(biāo)。因條件數(shù)(condition number，CN)通?？捎脕?lái)診斷GWR估算結(jié)果的局部共線(xiàn)性問(wèn)題[23]，本文還進(jìn)一步深入系統(tǒng)評(píng)估了GWR模型精度、奇異值甄別效果和回歸系數(shù)正負(fù)比例隨不同后驗(yàn)SD與CN閾值約束的敏感性。

GWR模型精度檢驗(yàn)過(guò)程中，本文采用基于樣本、站點(diǎn)及區(qū)域的3種十折交叉驗(yàn)證方法[24]。基于樣本的交叉驗(yàn)證是指將本文中所有236 206個(gè)樣本點(diǎn)的時(shí)空位置隨機(jī)打亂分成十組，換言之，同一個(gè)PM2.5站點(diǎn)上的樣本數(shù)據(jù)某些時(shí)刻可用于精度檢驗(yàn)，某些時(shí)刻可用于建模。基于站點(diǎn)的交叉驗(yàn)證則是將1356個(gè)地面站點(diǎn)隨機(jī)打亂分成十組，驗(yàn)證站點(diǎn)在所有時(shí)刻的樣本全用于模型檢驗(yàn)?；趨^(qū)域的交叉驗(yàn)證是將省級(jí)行政區(qū)隨機(jī)打亂分成十組，落在用于驗(yàn)證的省級(jí)行政區(qū)內(nèi)所有站點(diǎn)在所有時(shí)刻的樣本點(diǎn)均用于模型檢驗(yàn)。模型精度檢驗(yàn)指標(biāo)采用擬合優(yōu)度(coefficient of determination，R2)與均方根誤差(root mean square error，RMSE)。

根據(jù)不確定性約束的敏感性分析結(jié)果，擇優(yōu)選擇GWR模型與后驗(yàn)SD閾值實(shí)現(xiàn)地表PM2.5濃度制圖，并對(duì)比不確定性約束前后制圖結(jié)果，實(shí)證GWR建模結(jié)果不確定性度量與約束方法的有效性。

2.3 不確定性約束前GWR模型精度分析

圖3中4類(lèi)GWR模型帶寬時(shí)間序列值對(duì)比可知，自適應(yīng)型帶寬優(yōu)化結(jié)果對(duì)CV/AICc優(yōu)化指標(biāo)的選取并不敏感，模型A與模型B的帶寬均值相近；但固定型帶寬優(yōu)化結(jié)果在不同優(yōu)化指標(biāo)的選取下相差較大，模型D的帶寬均值明顯低于模型C。與此同時(shí)，表2中不確定性約束前模型A與模型B的擬合、樣本驗(yàn)證與站點(diǎn)驗(yàn)證精度明顯優(yōu)于模型C與模型D，說(shuō)明采用自適應(yīng)型帶寬確保局部回歸建模樣本數(shù)量相同可有效提升GWR建模精度。但對(duì)于固定型帶寬的模型C與模型D，因較低空間帶寬易導(dǎo)致某些區(qū)域局部回歸缺少充足建模樣本，因此模型D預(yù)測(cè)結(jié)果常產(chǎn)生較多奇異值，使得總體擬合與驗(yàn)證精度嚴(yán)重偏低。不僅如此，表2中基于區(qū)域的交叉驗(yàn)證結(jié)果顯示，4類(lèi)GWR模型的區(qū)域驗(yàn)證R2與RMSE結(jié)果均不理想，表明未經(jīng)不確定性約束的GWR模型預(yù)測(cè)結(jié)果可信度低。

圖3 4類(lèi)GWR模型帶寬時(shí)間序列

2.4 不確定性約束下GWR模型的敏感性分析

2.4.1 模型精度穩(wěn)定性分析

圖4繪制了不同后驗(yàn)SD閾值約束下，GWR模型擬合及十折交叉驗(yàn)證的R2與RMSE結(jié)果，以及相應(yīng)剔除率的變化情況。結(jié)果表明，4類(lèi)GWR模型擬合、十折交叉驗(yàn)證RMSE隨著后驗(yàn)SD閾值的降低而降低，后驗(yàn)SD約束可較好提升GWR模型精度。同時(shí)，除模型D外，其他3類(lèi)GWR模型的擬合與十折交叉驗(yàn)證R2分別呈現(xiàn)出隨后驗(yàn)SD閾值降低而下降和升高，且兩者之間差異逐步變小的結(jié)果，一定程度上展現(xiàn)了后驗(yàn)SD約束在緩解GWR模型過(guò)擬合現(xiàn)象中的作用。相比較表2中約束前的擬合與驗(yàn)證結(jié)果，在后驗(yàn)SD閾值僅設(shè)為30 μg/m3時(shí)，4類(lèi)GWR模型精度仍均有提高，且精度提升比例在模型D上尤為突出。模型D的擬合、樣本驗(yàn)證、站點(diǎn)驗(yàn)證及區(qū)域驗(yàn)證結(jié)果分別僅在剔除0.56%、2.89%、2.96%及21.15%的樣本后，RMSE提升比例分別可達(dá)92.99%、96.06%、97.10%及99.83%，說(shuō)明未經(jīng)后驗(yàn)SD約束的GWR預(yù)測(cè)結(jié)果中存在的少量奇異值嚴(yán)重影響了GWR模型總體精度。

表2 后驗(yàn)SD與CN約束前各類(lèi)GWR模型的擬合、樣本驗(yàn)證、站點(diǎn)驗(yàn)證及區(qū)域驗(yàn)證的R2/RMSE

圖4 不同后驗(yàn)標(biāo)準(zhǔn)差閾值約束下GWR模型擬合及十折交叉驗(yàn)證的R2、RMSE與剔除率統(tǒng)計(jì)

由圖5中CN閾值約束結(jié)果可知，雖然GWR模型基于樣本驗(yàn)證及站點(diǎn)驗(yàn)證的R2與RMSE結(jié)果隨著CN閾值降低均分別呈現(xiàn)出上升及下降趨勢(shì)，但擬合及基于區(qū)域驗(yàn)證的R2與RMSE結(jié)果卻呈現(xiàn)相反的變化趨勢(shì)。這一結(jié)果表明，傳統(tǒng)CN閾值約束或許不是有效提升GWR模型精度的手段。此外因?yàn)檫^(guò)多樣本剔除，即樣本點(diǎn)偏少易導(dǎo)致R2與RMSE等指標(biāo)計(jì)算不穩(wěn)定。圖5表明，當(dāng)CN閾值由50降至5時(shí)，樣本剔除率可由30%升至90%。

圖5 不同條件數(shù)閾值約束下GWR模型擬合及十折交叉驗(yàn)證的RMSE與剔除率統(tǒng)計(jì)

2.4.2 奇異值甄別效果分析

圖6與圖7分別展示了采用自適應(yīng)型帶寬與固定型帶寬的4類(lèi)GWR模型擬合及交叉驗(yàn)證散點(diǎn)圖，結(jié)果表明4類(lèi)GWR模型的預(yù)測(cè)結(jié)果均存在奇異值，但采用自適應(yīng)型帶寬的GWR模型奇異值在量級(jí)與數(shù)量上均要小于采用固定型帶寬的GWR模型。對(duì)于采用自適應(yīng)型帶寬的模型A與模型B，在地面站點(diǎn)PM2.5濃度較高時(shí)存在少量低估的GWR預(yù)測(cè)值，隨著后驗(yàn)SD閾值從30 μg/m3降至10 μg/m3時(shí)，這些低估的預(yù)測(cè)值可被逐步剔除。對(duì)于采用固定型帶寬的模型C與模型D，其GWR預(yù)測(cè)的奇異值卻可達(dá)1500 μg/m3及-40 000 μg/m3，說(shuō)明局部建模樣本不充足是產(chǎn)生奇異值的主要原因。后驗(yàn)SD對(duì)甄別奇異值仍具有較好效果，當(dāng)后驗(yàn)SD閾值降為10 μg/m3時(shí)奇異值均可被剔除。

注：紅實(shí)線(xiàn)為回歸線(xiàn)，黑虛線(xiàn)為1∶1參考線(xiàn)。

注：紅實(shí)線(xiàn)為回歸線(xiàn)，黑虛線(xiàn)為1∶1參考線(xiàn)。

2.4.3 模型局部共線(xiàn)性?xún)?yōu)化

在地表PM2.5濃度建模制圖過(guò)程中，由于建模變量在空間臨近范圍內(nèi)常表現(xiàn)出較強(qiáng)的相似性，因此易發(fā)生局部共線(xiàn)性，進(jìn)而導(dǎo)致GWR回歸系數(shù)符號(hào)與實(shí)際情況存在偏離。圖8顯示了4類(lèi)模型擬合結(jié)果在不同后驗(yàn)SD閾值/CN閾值約束下，符號(hào)為正的AOD回歸系數(shù)與符號(hào)為負(fù)的EVI回歸系數(shù)占總回歸系數(shù)比例的變化趨勢(shì)。由圖8可知，隨著后驗(yàn)SD閾值從30 μg/m3降至2 μg/m3時(shí)，4類(lèi)GWR模型正符號(hào)的AOD回歸系數(shù)比例提升范圍分別為75.67%～82.11%、75.24%～78.60%、90.09%～94.12%與78.22%～85.34%，而負(fù)符號(hào)的EVI回歸系數(shù)提升范圍分別為62.51%～66.38%、62.44%～67.43%、71.99%～78.51%與62.73%～67.40%。與之類(lèi)似，CN閾值約束結(jié)果同樣呈現(xiàn)了相似規(guī)律，當(dāng)CN閾值逐漸從100降低至10時(shí)，正符號(hào)的AOD回歸系數(shù)比例與負(fù)符號(hào)的EVI回歸系數(shù)比例也有所提升。這些改進(jìn)結(jié)果更加符合PM2.5濃度地理要素作用機(jī)制，即在地表PM2.5濃度關(guān)聯(lián)分析中，遙感AOD常呈現(xiàn)較強(qiáng)的正相關(guān)性[25]，遙感EVI則因森林葉面可吸附并捕獲PM2.5微粒常呈現(xiàn)負(fù)相關(guān)性[26]。綜上可知，后驗(yàn)SD與CN約束均能緩解局部共線(xiàn)性帶來(lái)的GWR建模過(guò)程中的系數(shù)符號(hào)偏離問(wèn)題，但考慮圖5中CN約束易導(dǎo)致過(guò)多樣本被剔除，后驗(yàn)SD約束方案更加合理。

圖8 不同后驗(yàn)SD閾值與CN閾值下，4類(lèi)GWR模型擬合結(jié)果AOD系數(shù)正數(shù)比例趨勢(shì)圖與EVI系數(shù)負(fù)數(shù)比例趨勢(shì)圖

2.5 地表PM2.5濃度制圖結(jié)果對(duì)比分析

在綜合考慮4類(lèi)GWR模型后驗(yàn)SD約束的建模效果基礎(chǔ)上，本文設(shè)定后驗(yàn)SD閾值為10 μg/m3，基于模型B模擬生成了GWR模型約束前/后研究區(qū)地表PM2.5濃度。圖9(a)、(b)分別為2019年2月6日與2019年9月25日Aqua衛(wèi)星PM2.5濃度模擬制圖結(jié)果。由圖9可知，未經(jīng)后驗(yàn)SD約束的PM2.5濃度GWR模擬估算結(jié)果中，監(jiān)測(cè)站點(diǎn)稀疏的西藏地區(qū)部分區(qū)域PM2.5濃度均高于100 μg/m3，與西藏地區(qū)地勢(shì)高且污染排放量極少的事實(shí)相違背；相比，后驗(yàn)SD約束后的PM2.5濃度則未出現(xiàn)這些奇異值。與此同時(shí)，圖10也展示了后驗(yàn)SD約束前/后PM2.5濃度GWR模擬估算結(jié)果中負(fù)數(shù)比例的時(shí)間序列及其頻數(shù)分布。圖10表明，約束后不合理PM2.5濃度負(fù)值比例大幅度降低，負(fù)數(shù)比例頻數(shù)分布由約束前多數(shù)天的5%～10%下降至了約束后的0%附近，春冬兩季更為明顯。此外，2019年站點(diǎn)PM2.5濃度年均值分布、后驗(yàn)SD約束后Terra/Aqua衛(wèi)星PM2.5濃度GWR模擬估算的年均制圖結(jié)果、及后驗(yàn)SD約束前后站點(diǎn)PM2.5濃度與GWR模擬估算制圖結(jié)果偏差頻數(shù)分布表明，后驗(yàn)SD約束后GWR模擬估算PM2.5濃度與地面站點(diǎn)PM2.5濃度偏差總體減小，污染分布空間特征識(shí)別基本一致，污染嚴(yán)重區(qū)域主要集中在華北平原、四川盆地及新疆塔克拉瑪干沙漠區(qū)域。

圖9 2019年2月6日至2019年9月25日不確定約束前后地表PM2.5濃度制圖結(jié)果

圖10 不確定性約束前后地表PM2.5濃度制圖結(jié)果中負(fù)數(shù)比例時(shí)間序列及其頻數(shù)分布

3 總 結(jié)

針對(duì)樣本空間稀疏程度及局部共線(xiàn)性等因素導(dǎo)致GWR模型精度呈現(xiàn)空間異質(zhì)性的問(wèn)題，本文提出了一種GWR建模結(jié)果不確定性度量與約束方法，并基于地表PM2.5濃度遙感制圖案例對(duì)其進(jìn)行了驗(yàn)證。結(jié)果表明，后驗(yàn)標(biāo)準(zhǔn)差約束下的GWR建模方法能有效甄別GWR模擬結(jié)果中的奇異值，提升模型的精度與模型模擬的可靠性。未來(lái)研究中，可進(jìn)一步專(zhuān)注GWR建模中間過(guò)程，利用后驗(yàn)標(biāo)準(zhǔn)差約束GWR模型參數(shù)的訓(xùn)練(如空間帶寬)，直接獲得更加穩(wěn)健的GWR模擬估算結(jié)果。