李巍　范文義　毛學剛　王蘭霞

摘要對反距離加權(quán)法、樣條函數(shù)法、普通克里格法及協(xié)同克里格法等幾種常用的降雨空間插值方法的優(yōu)缺點進行了分析和比較，考慮到高程對降雨量影響較大，在協(xié)同克里格法的基礎(chǔ)上將高程作為第2類影響因素引入降雨量的空間插值方法中，并提出了引入高程信息的協(xié)同克里格法。將4種方法（反距離加權(quán)法、樣條函數(shù)法、普通克里格法及協(xié)同克里格法）用于大興安嶺地區(qū)降雨量的插值計算，結(jié)果表明，考慮高程信息的協(xié)同克里格法的插值效果明顯優(yōu)于其他3種方法。

關(guān)鍵詞降雨量；空間插值方法；協(xié)同克里格法；高程信息

中圖分類號S161.6文獻標識碼A文章編號0517-6611（2014）12-03667-03

基金項目“十二五”國家科技支撐項目（2011BAD08B01）；黑龍江省青年基金（QC2012C102）；黑龍江省普通高校重點實驗室項目（KJKF1203）；黑龍江省教育廳科學技術(shù)研究項目（12511490，12531599）。

作者簡介李?。?982- ），女，黑龍江哈爾濱人，講師，在讀博士，從事3S技術(shù)研究及3S在土壤侵蝕、土地規(guī)劃方面研究。*通訊作者，教授，博士后，博士生導(dǎo)師，從事3S技術(shù)研究及應(yīng)用。

氣象資料是進行地學研究的一個基本而重要的參數(shù)，氣象資料可以反映出區(qū)域的森林植被分布、土壤侵蝕分布、水文生態(tài)的規(guī)律[1]。由于經(jīng)濟條件和技術(shù)手段的制約，大部分區(qū)域內(nèi)設(shè)置的氣象站點觀測密度均不高，這就使得研究者們無法在研究區(qū)內(nèi)得到連續(xù)有序的氣象空間數(shù)據(jù)分布。因此利用現(xiàn)有的氣象觀測站的數(shù)據(jù)通過空間插值對觀測數(shù)據(jù)進行補充就尤為重要。筆者利用大興安嶺地區(qū)2000～2010年降雨量資料，采用4種常用的降雨空間插值方法對降雨量插值進行比較分析。

1資料與方法

選取大興安嶺地區(qū)2000～2010年的多年平均降雨量資料，采用反距離加權(quán)法、樣條函數(shù)法、普通克里格法及協(xié)同克里格法4種方法對大興安嶺地區(qū)降雨量進行插值分析，并對插值結(jié)果進行檢驗，找出最優(yōu)插值方法。

1.1插值方法

1.1.1反距離加權(quán)法?，F(xiàn)階段，對地學研究最常用的插值方法就是反距離加權(quán)插值方法（IDW）[2]，它是從地理學第一定律的角度推出，也叫距離倒數(shù)乘方法。1972年，美國國家氣象局提出了反距離加權(quán)法。反距離加權(quán)法就是根據(jù)數(shù)據(jù)點間的空間距離遠近加權(quán)插字段進行的插值方法。距離中心越近的點，其估算值越受影響，中心點的影響隨著離它的距離越遠而減小。反距離加權(quán)插值的計算公式是：z（x0）=ni=1z（xi）（di0）/ni=11（di0）p，式中，z（x0）是插值點的預(yù)估值，z（xi）（i=1，2，3，…，n）是實測樣本值，n表示參與插值的實測樣本數(shù)；di0表示插值點和第i0個站點之間的距離；p是距離的冪，插值的結(jié)果被p影響，p越大，內(nèi)插后的效果越平滑，反之，p越小，內(nèi)插后的效果越尖銳，一般根據(jù)最小平均絕對誤差的大小來確定p。在日常使用中，經(jīng)常出現(xiàn)反距離平方加權(quán)法，就是令p=2的反距離加權(quán)法。如果p=3，就進行了反向距離立方插值。在ArcGIS中，p的默認值為2。此次降雨量插值時p的值設(shè)為2，站點搜索的范圍設(shè)為周邊6個，即n=6。

1.1.2樣條函數(shù)法。樣條函數(shù)（SPLINE）通過一個使表面整體曲率最小的數(shù)學函數(shù)來估計單元值，所得表面較為光滑。樣條函數(shù)法主要是針對一些特征點，對估計方差進行控制，適用多項式擬合的方法產(chǎn)生平滑的插值曲線。在降雨量空間插值中，樣條函數(shù)就是用一個最小曲率面來充分逼近各降雨量觀測點，進而推算出整個研究區(qū)的降雨量分布。樣條函數(shù)的計算公式為：Z=ni=1λiR（γi）+T（x，y），式中，Z為降雨量的預(yù)測值；n為參與此次插值的觀測站點數(shù)；λi是線性方程組求解確定的系數(shù)；γi是預(yù)估測點到第i點的距離，R（γi）和T（x，y）的表達式分別為：R（γi）=γ24[ln[γ2π]+c-1]+τ2[k0[γτ]+c]+ln[r2π]2π、T（x，y）=a1+a2x+a3y，式中，τ2為權(quán)重系數(shù)，γ為一直點與采樣點之間的距離，k0是修改后的貝塞爾函數(shù)，c為常數(shù)，a為線性方程的系數(shù)。此次插值所用的樣條函數(shù)為規(guī)則樣條函數(shù)，指數(shù)選擇0.1，點點搜索范圍為臨近的6個觀測站點。

1.1.3普通克里格插值法。普通克里金插值法（OK）是區(qū)域化變量的線性估計，是由法國的地理學家Matheron和南非礦山工程師Krige提出[3]，普通克里格法最初用于礦山勘探。它假設(shè)數(shù)據(jù)變化呈正態(tài)分布，通過對數(shù)據(jù)的空間分析獲取權(quán)重值，插值的整個過程相當于在未知區(qū)域化變量Z的期望值時對樣點進行加權(quán)滑動求取平均值的過程。普通克里格方法的公式為：Z（x0）=ni=1λiZ（xi）、ni=1λi=1，式中，Z（x0）為待插值點的估計值，Z（xi）為第i個樣本點的實測值，n為參與計算的觀測站樣本個數(shù)，λi為第i個樣本點的權(quán)重系數(shù)。權(quán)重λi的選擇必須保證Z（x0）能進行無偏估計，且估計的方差小于其他線性組合生成的方差?？死锔癫逯捣椒ǖ年P(guān)鍵問題是選擇合適的變異函數(shù)，目前經(jīng)常使用的變異函數(shù)包括高斯、考克斯（冪）、指數(shù)、球面及線性模型等。此次研究選擇考克斯變換模型作為基本變異函數(shù)模型，采用離預(yù)估點最臨近的6個觀測站點的數(shù)據(jù)進行插值。

1.1.4協(xié)同克里格插值法。Dirks等發(fā)現(xiàn)當站網(wǎng)密度較高時，普通克里金方法的插值效果與其他常用方法相比并無多大優(yōu)勢[4]。Borgam等也曾得到過同樣的結(jié)論[5]。但隨著高程的增加，降雨量有增加的趨勢。Hevesi等研究了年平均降雨量與高程的相關(guān)性，并得到了其相關(guān)系數(shù)達0.175的結(jié)果[6]。因此，在此采用協(xié)同克里格方法（OCK），并將高程作為第2影響因素引入到對降雨量的空間插值中來[7-8]。協(xié)同克里格法的公式為：Z（x0）=ni=1λiz（xi）+λ[y（x0）-my+mz]，式中，Z（x0）為x0點的預(yù)測值；z（xi）是第i站點的觀測值；y（x0）是x0點的高程；n為降雨量觀測站點的個數(shù)；my和mz為海拔高程和降雨量的全局平局值；λi和λ為協(xié)同克里格插值的權(quán)重系數(shù)。嘗試使用集中變異函數(shù)，結(jié)果顯示指數(shù)模型的插值效果最好，站點搜索范圍仍然為臨近6個站點。

1.2檢驗方法以上4種插值方法利用交叉檢驗的方法對插值的結(jié)果進行對比分析。交叉檢驗方法是首先假設(shè)一個站點的觀測值是未知的，用這個站點周圍的其他站點的觀測值來進行插值獲取估計值，再假設(shè)另一個站點的觀測值未知，用周圍其他站點的觀測值進行插值獲取估計值，以此輪換將17個站點分別作為未知站點，獲取17個站點的估計值，然后計算這17個站點的估計值和觀測值之間的誤差值，根據(jù)誤差值的大小來判斷4種插值方法效果的好壞。

利用平均相對誤差指標可以判斷出估計值和實測值之間的誤差大小，利用平均絕對誤差指標可以估算出插值獲取的估計值的誤差范圍，因此，在此利用MRE（平均相對誤差）、MAE（平均絕對誤差）這2個指標作為誤差的檢驗指標。MRE、MAE的表達式分別是：MRE=1nni=1Zai-ZeiZai、MAE=1nni=1|Zai-Zei|，式中，Zai是第i個站點的觀測值；Zei是第i個站點的插值估計值；n為作為交叉檢驗觀測站的數(shù)量。

2結(jié)果與分析

2.1多年平均降雨量插值為了更加直觀地看到各種插值方法的效果，圖1給出了大興安嶺地區(qū)不同方法得到的2000～2010年的多年平均降雨量插值結(jié)果。從圖1可知，各種方法得到的插值結(jié)果明顯不同，總體上大興安嶺地區(qū)的年降雨量在337～475 mm，強降水中心在呼瑪縣南部。

2.2不同插值模型的插值精度分析

2.2.1年尺度插值結(jié)果比較。從4種方法對多年年均降雨數(shù)據(jù)插值結(jié)果（圖2）來看，在所有的插值方法中，優(yōu)劣排序為OCK

2.2.2月尺度插值結(jié)果比較。利用2007年全年12個月的降雨量數(shù)據(jù)進行月尺度插值結(jié)果對比，從4種方法對2007年的月降雨數(shù)據(jù)插值結(jié)果（圖3）來看，6～9月的插值結(jié)果的不確定性遠遠小于11月～次年3月。

2.3降雨量空間插值精度驗證分析利用4種方法對多年

3結(jié)論

（1）從4種方法對多年年均降雨數(shù)據(jù)插值結(jié)果來看，在所有的插值方法中，優(yōu)劣排序為OCK

（2）從4種方法對2007年的月降雨數(shù)據(jù)插值結(jié)果來看，6～9月的插值結(jié)果的不確定性遠遠小于11月～次年3月。

（3）基于普通協(xié)克里格方法對大興安嶺地區(qū)4個及其周邊14個站點進行多年平均降雨量插值，得到研究區(qū)多年平均降雨空間分布特征，大興安嶺地區(qū)的年降雨量在358～462 mm，強降水中心在呼瑪縣南部，其最大值為462 mm，降水的低值區(qū)漠河縣西部，年降水量為358 mm。

參考文獻

[1] MOHAMED A S.Reliabilty estimation of rainfallrunoff models[D].New York：State University of New York，1999.

[2] LAMN.Spatial interpolation methods：a review[J].The Amercian Cartographer，1983，10（2）：129-149.

[3] 侯景儒，黃競先.地質(zhì)統(tǒng)計學的理論與方法[M].北京：地質(zhì)出版社，1990：69-78.

[4] DIRKS K N，HAY J E，STOW C D，et al.Highresolution studies of rainfall on Norfolk Island.Part：interpolation of rainfall data[J].J Hydrol，1998，208（3/4）：187-193.

[5] BORGAM，VIZZACCARO A.On the interpolation of hydrologic variables：formal equivalence of multiquadratic surface fitting and Kriging[J].J Hydrol，1997，195：160-171.

[6] HEVESI J A，F(xiàn)LINT A L，ISTO J D.Precipitation estimation in mountainous terrain using multivariate geostatistics.part I：structuralanalysis[J].J Appl Meteor，1992，31：661-676.

[7] 王家華，高海余，周葉.克里金地質(zhì)繪圖技術(shù)[M].北京：石油工業(yè)出版社，1999：157-179.

[8] 朱會義，賈紹鳳.降雨信息空間插值的不確定性分析[J].地理科學進展，2004，23（2）：34-42.