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

李巍　范文義　毛學(xué)剛　王蘭霞

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

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

中圖分類號(hào)S161.6文獻(xiàn)標(biāo)識(shí)碼A文章編號(hào)0517-6611（2014）12-03667-03

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

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

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

1資料與方法

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

1.1插值方法

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

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

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

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

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

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

2結(jié)果與分析

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

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

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

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

2.3降雨量空間插值精度驗(yàn)證分析利用4種方法對(duì)多年

3結(jié)論

（1）從4種方法對(duì)多年年均降雨數(shù)據(jù)插值結(jié)果來(lái)看，在所有的插值方法中，優(yōu)劣排序?yàn)镺CK

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

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

參考文獻(xiàn)

[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] 侯景儒，黃競(jìng)先.地質(zhì)統(tǒng)計(jì)學(xué)的理論與方法[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] 朱會(huì)義，賈紹鳳.降雨信息空間插值的不確定性分析[J].地理科學(xué)進(jìn)展，2004，23（2）：34-42.