降雨量空間插值方法比較研究

李巍　范文義　毛學剛　王蘭霞

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

關鍵詞降雨量；空間插值方法；協同克里格法；高程信息

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

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

作者簡介李巍（1982- ），女，黑龍江哈爾濱人，講師，在讀博士，從事3S技術研究及3S在土壤侵蝕、土地規劃方面研究。*通訊作者，教授，博士后，博士生導師，從事3S技術研究及應用。

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

1資料與方法

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

1.1插值方法

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

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

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

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

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

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

2結果與分析

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

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

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

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

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

3結論

（1）從4種方法對多年年均降雨數據插值結果來看，在所有的插值方法中，優劣排序為OCK

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

（3）基于普通協克里格方法對大興安嶺地區4個及其周邊14個站點進行多年平均降雨量插值，得到研究區多年平均降雨空間分布特征，大興安嶺地區的年降雨量在358～462 mm，強降水中心在呼瑪縣南部，其最大值為462 mm，降水的低值區漠河縣西部，年降水量為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] 侯景儒，黃競先.地質統計學的理論與方法[M].北京：地質出版社，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，FLINT 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] 王家華，高海余，周葉.克里金地質繪圖技術[M].北京：石油工業出版社，1999：157-179.

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