鄧　悅，劉　洋，劉紀(jì)平，徐勝華，羅　安

(1. 中國測繪科學(xué)研究院，北京 100830； 2. 中南大學(xué)地理信息系，湖南 長沙 410083)

據(jù)統(tǒng)計，2015年中國南方地區(qū)共出現(xiàn)5次區(qū)域性暴雨過程，造成全國15個省(自治區(qū)、直轄市)123.6萬人次受災(zāi)，近3000間房屋倒塌，1.2萬間不同程度損壞，農(nóng)作物受災(zāi)面積5.86×104hm2。因此，洪澇和地質(zhì)災(zāi)害已成為中國公眾及政府部門重點關(guān)注的問題。但目前我國山區(qū)的氣象站和雨量站多設(shè)在河谷低處，在海拔較高的地方雨量測點很少。利用如今分布的站點觀測資料計算山區(qū)的面雨量，是很不可靠的[1]。

目前，已有大量對降雨量空間分布模擬的研究，主要有以下兩個方向：①數(shù)學(xué)建模方法：傅抱璞[1]提出了一個表示降水與地形、海拔高度關(guān)系的數(shù)學(xué)模型；MarquíNez等[2]利用回歸模型分析了西班牙北部區(qū)域降水和一系列地形變量的關(guān)系；Naoum等[3]建立了希臘Crete島的年降雨多元線性回歸模型；Boni等[4]采用多元線性回歸方法分析了降雨與氣候、地形因子的關(guān)系；彭曉芬等[5]通過研究對比發(fā)現(xiàn)回歸模型模擬方法適用于云南地區(qū)的空間模擬，但線性回歸模型并不適應(yīng)于降雨量這種隨時間呈某種變化趨勢的非平穩(wěn)隨機過程[6]。②空間插值方法：白江濤等[7]發(fā)現(xiàn)協(xié)同克里金插值方法更適合于降雨量空間分布模擬；何紅艷等[8]采用5種插值方法對青藏高原及其周邊地區(qū)降雨量進(jìn)行了模擬；Delbari等[9]利用多種插值方法對伊朗地區(qū)進(jìn)行了降雨量空間分布研究；SHEN Zhenyao等[10]采用了兩種克里金法對降雨量進(jìn)行空間插值方法比較。上述方法共有的缺點在于誤差較大，不同地區(qū)所適用的插值方法不同，甚至同一地區(qū)也無法用同一種插值模型進(jìn)行模擬。

Brunsdon和Fotheringham等提出的地理加權(quán)回歸(geographic weighted regression,GWR)模型是一種局部加權(quán)回歸模型[11-12]，可以模擬地理現(xiàn)象局部特征，但往往會因為異常值或“弱數(shù)據(jù)”的出現(xiàn)而對局部參數(shù)估計產(chǎn)生較大的影響；LeSage提出了貝葉斯地理加權(quán)回歸(Bayesian geographic weighted regression,BGWR)模型[14-15]，可以彌補GWR模型無法處理異常觀測值和“弱數(shù)據(jù)”的不足；Nirmal Subedi等證實了BGWR模型可以有效改進(jìn)傳統(tǒng)的地理加權(quán)回歸，有助于提高樹木或植物群落的空間結(jié)構(gòu)預(yù)測精度[16]；崔長彬等利用BGWR對河北省136個縣市經(jīng)濟影響因素進(jìn)行了回歸分析[17]；丁剛等利用BGWR對政府效能對創(chuàng)新能力建設(shè)的影響進(jìn)行了客觀評價[18]。但鮮有學(xué)者將貝葉斯地理加權(quán)回歸模型應(yīng)用于降雨量預(yù)測。

因此，本文提出了一種結(jié)合吉布斯采樣的貝葉斯地理加權(quán)回歸方法(GBGWR)，進(jìn)行降雨量空間分布模擬。該方法利用GWR模型的空間非平穩(wěn)性探測能力，解決了如OLS等線性模型對局部地理現(xiàn)象模擬不精確的問題，又將最大后驗概率作為先驗信息引入GWR模型，彌補GWR模型無法處理異常觀測值和“弱數(shù)據(jù)”的不足。相較于傳統(tǒng)回歸分析，本文提出的GBGWR方法能夠更好地進(jìn)行降雨量空間分布模擬。

1　研究區(qū)和數(shù)據(jù)

1.1　研究區(qū)概況

湖南省地處我國腹地，全省地勢東高西低，處于云貴高原與江漢平原的過渡地帶，地貌形態(tài)復(fù)雜多樣，屬于典型的亞熱帶季風(fēng)氣候區(qū)，處在東南季風(fēng)和西南季風(fēng)相交的地帶，降水分布不均勻，年際變化較大，旱澇災(zāi)害發(fā)生頻率高。因此，本文選擇湖南省作為試驗區(qū)域，對其進(jìn)行降水空間分布進(jìn)行模擬。

1.2　數(shù)　據(jù)

本文采用的試驗數(shù)據(jù)分別為：①中國氣象科學(xué)數(shù)據(jù)共享網(wǎng)(data.cma.cn)提供的湖南省35個監(jiān)測站點1985—2015年的月值氣象觀測資料，對于部分缺失數(shù)據(jù)，首先采用MarquíNez等的數(shù)據(jù)不完整站點和完整站點之間的回歸方程法，對降雨數(shù)據(jù)進(jìn)行插補延長，最終整理得到湖南省降雨量及其相關(guān)因子30年的持續(xù)觀測數(shù)據(jù)，時間分辨率為每月。表1 給出了降雨量以及相關(guān)變量說明。

圖1　研究區(qū)站點分布

變量名稱類型單位降雨量因變量mm平均最高溫自變量0.1℃平均氣溫自變量0.1℃平均風(fēng)速自變量0.1m/s平均水汽壓自變量0.1hPa高程自變量m

2　模型及算法過程

2.1　基于吉布斯采樣的貝葉斯地理加權(quán)回歸模型

GBGWR模型相較于GWR模型的優(yōu)點在于不僅可以將地理位置納入模型之中，而且能夠通過將參數(shù)的先驗信息融入模型來消除或減少異常值對模型的影響，表述為

Wiy=WiXβi+Wiεi

(1)

(2)

2.2　吉布斯采樣

參數(shù)估計是貝葉斯地理加權(quán)回歸方法中的一個重要步驟，運用吉布斯采樣計算參數(shù)的后驗概率是最常用的方法之一。它從間接的概率分布生成隨機樣本，而不必計算密度本身，GBGWR的吉布斯采樣程序可以歸納如下：

(1) 從參數(shù)β(Si)0、σ(Si)0、δ0、V(Si)0、γ0的任意值開始。

(2) 對每一個觀測值i=1，2，…，n：

a. 從p(β(Si)|δ0，σ(Si)0，V(Si)0，γ0)中抽樣出一個β(Si)1；

b. 從p(σ(Si)|δ0，V(Si)0，β(Si)1，γ0)中抽樣出一個σ(Si)1；

c. 從p(V(Si)|δ0，β(Si)1，σ(Si)1，γ0)中抽樣出一個V(Si)1。

(3) 使用上述n次循環(huán)中的每一次采樣值β(Si)1，i=1，2，…，n來更新γ0、γ1的值。

(5) 返回步驟(1)，利用更新后的值β(Si)1、σ(Si)1、δ1、V(Si)1、γ1繼續(xù)進(jìn)行抽樣。

由此可見，該方法主要是通過對每一個觀測值的單獨迭代來模擬模型參數(shù)，通過捕獲較大方差得到異常值，并且在估計系數(shù)的同時使用加權(quán)回歸方法和方差的導(dǎo)數(shù)作為權(quán)重。

2.3　最優(yōu)帶寬的選擇

帶寬的選擇對于GBGWR至關(guān)重要，傳統(tǒng)的GBGWR采用CV交叉驗證法來確定帶寬。但CV交叉驗證及AIC信息準(zhǔn)則下的最優(yōu)帶寬并不能適應(yīng)數(shù)據(jù)點疏密不均的情況。因此，本文在傳統(tǒng)GWR模型的基礎(chǔ)上對帶寬選擇方法進(jìn)行了改進(jìn)，采用自適應(yīng)權(quán)函數(shù)來確立帶寬，從而使得模型不受數(shù)據(jù)點空間分布形式的影響。其中Rij代表數(shù)據(jù)點j的序號，第k個臨近數(shù)據(jù)點(Rij=k)的距離dij為帶寬。

(3)

2.4　算　法

GBGWR模型有兩個關(guān)鍵步驟：①采用吉布斯采樣方法計算變量的后驗概率，并通過貝葉斯公式計算回歸變量的先驗信息，然后對局部空間內(nèi)的各個變量進(jìn)行加權(quán)，如若檢測出異常觀測值，則對該觀測值進(jìn)行降權(quán)以降低其對參數(shù)估計的影響；②對于“弱數(shù)據(jù)”問題，通過在貝葉斯模型中引入顯式參數(shù)平滑關(guān)系(如式(2))，對參數(shù)變化的空間性質(zhì)施加限制，最后結(jié)合主觀先驗信息求出系數(shù)權(quán)重矩陣。其算法流程如圖2所示。

圖2　算法流程

3　試驗分析

本文分別采用GWR、BGWR及GBGWR 3種方法對湖南省35個監(jiān)測站點1985—2015年的降雨量觀測數(shù)據(jù)集進(jìn)行回歸，分別從模型系數(shù)、模型性能及空間分布模擬3個方面展開對比分析。建模前利用VIF進(jìn)行多重共線性診斷，結(jié)果表明解釋變量之間不存在多重共線性，最終建立了如下模型：其中i為數(shù)據(jù)點編號，β為自變量系數(shù)，εi為殘差項，R代表降雨量，AMT代表平均最高溫，AT代表平均氣溫，AWS代表平均風(fēng)速，AVP平均水汽壓，E代表高程。

Ri=β1i+β2i(AMTi)+β3i(ATi)+β4i(AWSi)+

β5i(AVPi)+β6i(Ei)+εi

對該模型進(jìn)行最小二乘估計，主要用于對模型的初步探測，保證后續(xù)回歸的可行性，同時與GWR、BGWR及GBGWR 3種方法得到的結(jié)果進(jìn)行對比分析，最小二乘估計結(jié)果見表2。

表2　最小二乘估計結(jié)果

3.1　系數(shù)對比分析

本文所用的GWR、BGWR及GBGWR是根據(jù)式(3)中提出的指數(shù)距離衰減函數(shù)產(chǎn)生的，使用式(2)中提出的一階鄰接先驗調(diào)整關(guān)系進(jìn)行參數(shù)調(diào)整，以此減少估計偏差。兩種方法得到的回歸系數(shù)(如圖3 所示)，通過在單一圖表中表示的3種方法對35個監(jiān)測站點的估計值可以看出：①貝葉斯先驗信息對空間平滑的影響主要是通過限制每個監(jiān)測站點的系數(shù)與鄰近監(jiān)測站點的系數(shù)相似。②如果35個監(jiān)測站點的GBGWR估計值及BGWR估計值與GWR估計值有很大的偏差，那么表明貝葉斯先前的信息在估計中引入了大量的偏差。但從圖中可以看出，3種方法得到的常數(shù)項值估計值均在-890左右變化，平均最高溫在1.1左右變化，平均氣溫在-4.9左右變化，平均風(fēng)速在3.8左右變化，平均水汽壓在16.3左右變化，高程在0.45左右變化，即均接近于最小二乘估計值，GBGWR估計和BGWR估計只是GWR估計的平滑版本，則表明先前的信息有助于強化對異常值的估計和非恒定方差。③基于吉布斯采樣的貝葉斯地理加權(quán)回歸結(jié)果與貝葉斯地理加權(quán)回歸結(jié)果整體上一致，表明在小樣本的情況下,不考慮樣本規(guī)模時,吉布斯抽樣能得到較為正確的后驗分布。

圖3　地理加權(quán)回歸模型、貝葉斯地理加權(quán)回顧模型和基于吉布斯采樣的貝葉斯地理加權(quán)回歸模型估計

3.2　模型性能對比分析

表3　GWR和GBGWR性能對比分析

3.3　空間分布特征分析

本文分別選取2015年夏季即6—8月湖南省35個監(jiān)測站點的月均值降雨量數(shù)據(jù)，采用反距離加權(quán)進(jìn)行插值得到降雨量空間分布圖。圖4分別為這3個月湖南省月度降雨量實測值與模擬值空間分布對比圖。實測值與模擬值總體呈一致狀態(tài)，6、7、8三個月份的降雨量空間分布也總體程一致狀態(tài)。從圖中可以看出，降雨量在空間上主要呈現(xiàn)北高南低，從北向南依次降低的特征，中部地區(qū)有一條明顯的綠色條帶將湖南省分為南中北3個降水量區(qū)域；其次，湖南省的降雨分布主要與大氣環(huán)流和地形有關(guān)，邵陽、衡陽等地降雨量明顯低于其他地區(qū)，主要原因在于衡邵盆地6—8月受西太平洋副熱帶高壓的控制，表現(xiàn)出長時間的晴熱干燥天氣，同時由于其處于南嶺北側(cè)，為東南季風(fēng)的背風(fēng)區(qū)，易產(chǎn)生焚風(fēng)效應(yīng)；從時間上來看，每年夏季6月為一年中平均降水量最大的月份，之后依次遞減，但總體趨勢呈現(xiàn)一致的狀態(tài)。

4　結(jié)　論

本文使用GWR模型、BGWR模型及GBGWR對湖南省35個監(jiān)測站點30年的降雨量數(shù)據(jù)進(jìn)行空間分布模擬，相較于GWR和BGWR，GBGWR很好地結(jié)合了貝葉斯理論與地理加權(quán)回歸二者的優(yōu)點，尤其是對于降雨量監(jiān)測站點稀疏、空間分布不均、降雨量的不穩(wěn)定性有著良好的適用性。試驗表明，該方法對于降雨量的空間分布模擬具有更好的穩(wěn)健性和更佳的模擬效果。綜上所述，得出如下結(jié)論：

圖4　2015年湖南省降雨量空間分布實際與模擬對比

(1) 湖南省夏季的降雨量整體上表現(xiàn)為北高南低、由南往北依次減少的分布特征，且降雨量越大時空間分布上的差異越大。

按WoS分類方法，全部論文分布于160種學(xué)科中，其中20個學(xué)科發(fā)文超過了100篇，是創(chuàng)業(yè)研究的主要學(xué)科載體。其中，商學(xué)以3193篇的發(fā)文優(yōu)勢力壓群雄，占據(jù)霸主地位，無疑是創(chuàng)業(yè)研究的最重要學(xué)科載體；管理學(xué)以微弱的劣勢屈居第二，發(fā)文3150篇，同樣顯示了其在創(chuàng)業(yè)研究中的重要地位。位于第三的是顯著落后于前兩者，但以較大優(yōu)勢領(lǐng)先于居第四位的規(guī)劃管理學(xué)的經(jīng)濟學(xué)，發(fā)文1831篇。這三大學(xué)科的絕對發(fā)文優(yōu)勢彰顯了其在創(chuàng)業(yè)研究中最為重要的學(xué)科載體地位。在學(xué)科分布中，中國與國際基本同步，管理學(xué)、商學(xué)和經(jīng)濟學(xué)依次為中國發(fā)文較多的前三種學(xué)科。

(2) GBGWR模型的平均擬合優(yōu)度達(dá)到了0.87，遠(yuǎn)高于GWR模型的0.73和BGWR模型的0.81，主要原因在于在小樣本的情況下,不考慮樣本規(guī)模時,吉布斯抽樣能得到較為正確的后驗分布。同時GBGWR模型估計系數(shù)對異常值敏感。當(dāng)鄰域中存在異常值時，GWR方法需要通過剔除異常值重新建立模型。相比之下，GBGWR可以自動識別異常值，并且通過降低異常值權(quán)重來減少異常值的影響，并且得到的估計系數(shù)是穩(wěn)健的。

由于降雨量隨著地形和季節(jié)變化的波動性較大[19]，因此下一步的研究工作主要集中在GBGWR模型的基礎(chǔ)上納入地形因素和時間因素，從貝葉斯時空地理加權(quán)的角度進(jìn)一步優(yōu)化模型。

參考文獻(xiàn)：

[1] 傅抱璞.地形和海拔高度對降水的影響[J].地理學(xué)報，1992(4):302-314.

[3] NAOUM S，TSANIS I K.A Multiple Linear Regression GIS Module Using Spatial Variables to Model Orographic Rainfall[J].Journal of Hydroinformatics，2004，6(1):39-56.

[4] BONI G，PARODI A，SICCARDI F.A New Parsimonious Methodology of Mapping the Spatial Variability of Annual Maximum Rainfall in Mountainous Environments[J].Journal of Hydrometeorology，2007，9(3):492-506.

[5] 彭曉芬，黃甫則，周汝良.云南省年均降雨量空間插值模擬方法比較[J].西南林業(yè)大學(xué)學(xué)報，2010，30(5):25-28.

[7] 白江濤，白建軍，王磊,等.基于GIS的關(guān)中——陜南地區(qū)降雨量空間插值分析[J].安徽農(nóng)業(yè)科學(xué)，2011，39(33):20872-20876.

[8] 何紅艷，郭志華，肖文發(fā),等.利用GIS和多變量分析估算青藏高原月降水[J].生態(tài)學(xué)報，2005，25(11):2933-2938.

[9] DELBARI M，AFRASIAB P,JAHANI S，et al.Spatial Interpolation of Monthly and Annual Rainfall in Northeast of Iran[J].Meteorology & Atmospheric Physics，2013，122(1-2):103-113.

[10]SHEN Z，CHEN L，LIAO Q，et al.Impact of Spatial Rainfall Variability on Hydrology and Nonpoint Source Pollution Modeling[J].Journal of Hydrology，2012，s 472-473(24):205-215.

[11]覃文忠,王建梅,劉妙龍.地理加權(quán)回歸分析空間數(shù)據(jù)的空間非平穩(wěn)性[J].遼寧師范大學(xué)學(xué)報(自然科學(xué)版),2005,28(4):476-479.

[12]BRUNSDON C.Geographically Weighted Reg Ression:A Method for Exploring Spatial Nonstation Arity[J].Geographical Analysis，1996，28(4)：281-298.

[13]武文娟，徐京華，時進(jìn)，等.基于GWR的四川省醫(yī)院床位數(shù)時空分布及其影響因素研究[J].測繪通報，2016(4)：49-53.

[14]LESAGE J P.A Spatial Econometric Examination of China’s Economic Growth[J].Geographic Information Sciences，1999，5(2)：143-153.

[15]LESAGE J P.A Family of Geographically Weighted Regression Models[C]∥Advances in Spatial Econometrics.Berlin Heidelberg：Springer，2004：241-264.

[16]SUBEDI N.Local Modeling of Tree Crown Area Using Bayesian Geographically Weighted Regression[D].[S.l.]：State University of New York College of Environmental Science and Forestry，2012.

[17]崔長彬，姜石良，張正河.河北縣域經(jīng)濟影響因素的空間差異分析——基于貝葉斯地理加權(quán)回歸方法[J].經(jīng)濟地理，2012，32(2)：39-45.

[18]丁剛.基于BGWR模型的區(qū)域創(chuàng)新能力建設(shè)中政府效能測評[J].技術(shù)經(jīng)濟，2013，32(4)：20-26.

[19]葉金桃，王俊超，馮光柳，等.降雨型滑坡聯(lián)合監(jiān)測及預(yù)警預(yù)報系統(tǒng)研究[J].測繪通報，2012(2)：31-34.