鄧　悅，劉　洋，劉紀平，徐勝華，羅　安

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

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

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

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

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

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

1.1　研究區(qū)概況

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

1.2　數(shù)　據(jù)

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

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

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

2　模型及算法過程

2.1　基于吉布斯采樣的貝葉斯地理加權回歸模型

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

Wiy=WiXβi+Wiεi

(1)

(2)

2.2　吉布斯采樣

參數(shù)估計是貝葉斯地理加權回歸方法中的一個重要步驟，運用吉布斯采樣計算參數(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ù)進行抽樣。

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

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

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

(3)

2.4　算　法

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

圖2　算法流程

3　試驗分析

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

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

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

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

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

3.1　系數(shù)對比分析

本文所用的GWR、BGWR及GBGWR是根據(jù)式(3)中提出的指數(shù)距離衰減函數(shù)產(chǎn)生的，使用式(2)中提出的一階鄰接先驗調(diào)整關系進行參數(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左右變化，平均風速在3.8左右變化，平均水汽壓在16.3左右變化，高程在0.45左右變化，即均接近于最小二乘估計值，GBGWR估計和BGWR估計只是GWR估計的平滑版本，則表明先前的信息有助于強化對異常值的估計和非恒定方差。③基于吉布斯采樣的貝葉斯地理加權回歸結(jié)果與貝葉斯地理加權回歸結(jié)果整體上一致，表明在小樣本的情況下,不考慮樣本規(guī)模時,吉布斯抽樣能得到較為正確的后驗分布。

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

3.2　模型性能對比分析

表3　GWR和GBGWR性能對比分析

3.3　空間分布特征分析

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

4　結(jié)　論

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

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

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

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

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

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

參考文獻：

[1] 傅抱璞.地形和海拔高度對降水的影響[J].地理學報，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è)大學學報，2010，30(5):25-28.

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

[8] 何紅艷，郭志華，肖文發(fā),等.利用GIS和多變量分析估算青藏高原月降水[J].生態(tài)學報，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]覃文忠,王建梅,劉妙龍.地理加權回歸分析空間數(shù)據(jù)的空間非平穩(wěn)性[J].遼寧師范大學學報(自然科學版),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]武文娟，徐京華，時進，等.基于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)濟影響因素的空間差異分析——基于貝葉斯地理加權回歸方法[J].經(jīng)濟地理，2012，32(2)：39-45.

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

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