黃土高原降雨量空間插值精度比較——KRIGING與 TPS法

萬龍，馬芹，張建軍，付艷玲，張曉萍

(西北農林科技大學資源環境學院;中國科學院水利部水土保持研究所;黃土高原土壤侵蝕與旱地農業國家重點實驗室:712100，陜西楊凌)

降雨是地表水資源的主要補給來源，認識降雨時空演變規律對科學利用和調控有限地表水資源尤為重要。近年來，隨著GIS在多種領域的廣泛應用，很多學者對于降雨量空間插值方法的優劣性做了深入研究，認為各種插值方法都有自己的優點，但也有明顯的缺陷［1-3］。如泰森多邊形方法(Thiessen)適用于站點密集地區，要求地形大致相同，且不考慮高程影響［4］，反距離加權法(IDW)通過權重調整空間插值等值線結構，但該方法也不考慮地形因素對降雨的影響［5］?；诮y計插值技術的克里金插值方法(KRIGING)利用半方差函數來反映插值對象的空間相關性，被認為是空間數據的最佳內插法［6］。薄板光滑樣條函數方法(Thin Plate smoothing Spine-TPS)是樣條函數方法的一種，兼顧了插值曲面的平滑性與精確度［7］。除了插值方法外，站點數量及分布均勻程度、時間尺度變化和柵格單元大小等因素，均會引起空間信息的不確定性［8-9］。

黃土高原地區年降雨量在300～650 mm之間，東南季風、地理位置及地勢構架因素的影響，空間上呈現出降雨量由東南向西北遞減的趨勢，時間上則表現為年際變化大，月際分配極不均勻的特點。受地面強烈切割和起伏的影響，稀疏覆被地表加強了夏季熱力對流作用，降雨形式多暴雨和冰雹［10］。在氣象站點觀測值基礎上，如何準確估算站點外區域的氣象信息，是科技和管理人員實踐中經常面對的問題。KRIGING插值方法是國內外文獻中通過站點數據獲取空間信息常用的方法［11］。為氣象數據曲面擬合而編寫的專用軟件ANUSPLIN是基于樣條插值理論，允許引進多元協變量的線性子模型，近年來在國際上得到了廣泛應用［12］。針對黃土高原降雨時空變異性大的特征，哪種方法更能準確反應氣候要素的空間信息，少有文獻探討。筆者擬利用黃土高原的降雨數據對這2種方法在多年平均、年和月尺度上進行插值精度比較，以期為土壤侵蝕模型的建立提供參考。

1 研究區概況

內蒙古河口鎮—陜西龍門區間(河龍區間)是黃河中游處于山西—陜西交界的區域，面積約11.3萬km2，以黃土丘陵溝壑地貌為主，海拔 450～2 079 m，切割深度在100～500 m之間，≥25°陡坡面積比例10% ～40%，地表割裂度30% ～70%。東部以呂梁山為界，西南鄰北洛河，西北為毛烏素沙地。屬于干旱半干旱氣候，年內降雨集中于5—9月，暴雨強度大，因引起劇烈侵蝕，每年輸入黃河的泥沙量超過三門峽以上總輸沙量的70%［13］。河龍區間植被以草原為主，北部為典型草原和荒漠草原，南部為灌叢草原、草甸草原和森林草原。區間內黃綿土廣泛分布，由南向北逐漸粗化。截至1996年底，河龍區間總人口約800萬，其中農業人口約697.2萬，土地利用不合理，水土流失嚴重。

2 數據資料及研究方法

2.1 數據資料

黃河中游河龍區間及毗鄰區共包括50個氣象站點，其中27個站點分布于河龍區間(圖1)，站點分布經緯度范圍為 E106.5°～113.5°，N35°～40.5°。50個氣象站點1981—2000年的月降雨數據從國家氣象站獲得。各站點月降雨量數據累加，計算多年平均降雨量值。在ArcGIS支持下使用KRIGING方法和采用ANUSPLIN專用軟件的TPS方法，對多年平均、年及月降雨量進行插值，利用河龍區間內27個站點數據進行交叉驗證和對比，比較2種方法插值精度的優劣。插值設置投影參數為Albers等積投影，第1條緯線為N 25°，第2條緯線為N 47°，中央經線為 E 110°，插值柵格為500 m ×500 m。研究區氣象站點平均分布密度為4200 km2/站。

2.2 KRIGING插值法

KRIGING插值法是由法國地理學家Georges Martheron和南非采礦工程師D.G.Krige創立的地質統計學的最佳內插方法［14］，近年來被廣泛應用于地學研究中。區域化變量和變異函數是克里金插值的基石。

圖1 河龍區間氣象站點分布圖Fig.1 Location of the meteorological stations in the He-Long Section

當空間點x在一維x軸上變化時，把區域化變量在x與x+h處的Z(x)與Z(x+h)差的方差之半定義為區域化變量Z(x)在x軸方向上的變異函數，即半方差函數，記為γ(h)，在滿足二階平穩假設條件下，有

式中:N(h)是研究區內距離為h的點對數，Z(x)和Z(x+h)分別為x及x+h點的變量值。

圖2描述了半方差函數與站點點對距離h之間的關系，隨著距離h的增加，觀測點的半方差增大，相關性也逐漸降低。當h達到一定程度，即h=a(變程)時，γ(h)逐漸趨近于定值。理論上γ(h)在h=0時γ(h)=0，但是，有時會在原點附近出現不連續的現象。這種現象稱為塊金效應(Nugget Effect)，這時的γ(h)就叫做塊金效應值(C0)，而(C+C0)為基臺值。當C0/(C+C0)比較小時，各采樣點的加權值變化比較大，數據的空間自相關性較好。隨著C0/(C+C0)的增大，空間數據相關性變差，各個采樣點的加權值變化減小，數據的空間自相關程度降低。塊金效應值等于基臺值時，即純塊金效應，說明采樣數據沒有包含空間相關性的任何量化信息，這種情況的估計結果就完全是采樣點的算術平均值，各個采樣點的加權值都一樣大［15］。

圖2 半方差函數變異曲線Fig.2 Semi-variogram variation curve

由于計算得到的實驗半方差函數是離散的且變化也往往不規則，在克里金估計中難以直接應用于指導克里金插值，所以要用理論模型去擬合。半方差理論模型常采用的模型有球面模型，模型參數通過實驗數據用最小二乘法擬合得到。文中變程設置為各點對的平均距離50 967 m，用離估計點最近的12個站點的數據來進行插值。

擬合半方差重要的用途是確定局部內插需要的權重因子，其值按采樣點數據的半方差圖的統計分析原理計算，即:

式中:Zi(xi)為已知第i個站點的觀測值;λi為第i個站點的權重;^Z(x0)為待估點x0的預測值。λi權重的選擇應使 ^Z(x0)是無偏估計。

2.3 TPS方法

薄板樣條函數(TPS)插值法，具有連續、光滑的數學特性，是自然樣條函數在多維空間的推廣。它除普通的樣條自變量外允許引入線性協變量子模型，如溫度，海拔等。由于黃土高原降雨與海拔無明顯的相關關系;故僅采用以經度和緯度為雙變量薄板光滑樣條函數進行插值，沒有引入海拔值作為協變量因子［16］，且文中樣條次數采用通常的二次樣條函數。

TPS 法插值原理［17］如下:

式中:zi為位于空間i點的因變量;xi為d維樣條獨立變量矢量;f(xi)為要估算的關于xi的未知光滑函數;yi為p維獨立協變量矢量;b為yi的p維系數;ei為具有期望值為0和方差為wiσ2的自變量隨機誤差，wi為作為權重的已知相對誤差方差;σ2為所有數據點上為一常數的誤差方差 ，但通常未知;N為觀測值個數。

函數f和系數b通過下述最小二乘法估計確定:

式中:Jm(f)為函數f(xi)的粗糙度測度函數，定義為函數f的m階偏導(在ANUSPLIN中稱為樣條次數，也叫糙度次數);ρ為正的光滑參數，在數據保真度與曲面的粗糙度之間起平衡作用，本文中由廣義交叉驗證GCV(generalized cross validation)的最小化來確定［18］。ANUSPLIN軟件在日志文件(Log file)中提供了用于判別誤差來源和插值質量的統計參數，可用于確定擬合過程中的最優光滑參數。

2.4 誤差分析方法

采用交叉驗證的方法對數據進行比較，依次移去一個點，用所有其他站點插值，得到被移去站點的預測值，該預測值與實測值的誤差為該站點的交叉驗證誤差。常用的評價插值精度的方法有標準均方根誤差(Root Mean Squared Interpolation Error，E)，可以反映估計值與實測值誤差的絕對量［19］。為了消除對異常值的敏感程度，同時利用一致性指標(the index of agreement，A)來評價KRIGING方法和TPS方法插值精度的優劣［20］。公式如下:

式中:Oi為第i個站點的實際觀測值;Pi為第i個站點的交叉驗證預測值;為n個站點的實際觀測值的平均值;n為用于參與驗證的站點數。A值在0～1之間，1表示預測值與實測值非常吻合，0表示不吻合。

2.5 全域型Moran’s I指數

計算空間自相關的方法有許多種，最為知名也最為常用的是Moran’s I方法。用全域型Moran’s I指數來表示研究區要素空間自相關程度。全域型Moran’s I指數是基于統計學相關系數的變異數和共變量關系推算得來。假設數組 xi全域型的Moran’s I的公式如下:

式中:I為moran’s I指數值，Wij為研究范圍內每一個空間單元 i與空間單元 j(i、j=1，2，…，n)的空間相鄰權重矩陣，以1表示i與j相鄰，以0表示i與j不相鄰;n為數組xi的元素個數。

I值標準化Z(I)的公式為

式中:E(I)為I的期望值;Var(I)為I的方差。對I值進行顯著性確定時，在5%顯著水平下，Z(I)大于1.96時，表示研究范圍內某現象的分布有顯著的空間自相關性。

3 結果與分析

3.1 多年平均降雨量插值結果比較

受東南季風的影響，河龍區間多年(1981—2000年)平均降雨量由東南向西北遞減，其空間變異情況為:平均值424.0 mm，極值比1.66，標準差52.9 mm，離差系數0.125，Moran’s I為0.57，統計值 Z=8.0，空間上具有極顯著的自相關性。KRIGING與TPS方法的插值結果均能很好地反映多年平均降雨量的空間變化趨勢，具有大致相同的空間變化梯度，TPS方法插值法降雨量的等值線更加平滑和柔和(圖3)。

對河龍區間內27個站點進行KRIGING與TPS方法交叉驗證，標準均方根誤差E分別為21.2和21.7 mm，一致性指標A分別為0.775和0.765。這2種方法對河龍區間多年平均降雨量的插值精度差別很小，E差值為多年平均降雨量的0.12%，A相差0.01。KRIGING方法(圖3(a))要略微優于TPS方法(圖3(b))。從數據點對的半方差云圖(圖3(c))可知，數據的空間相關性較好，無塊金效應，插值精度較高。

圖4顯示了KRIGING與TPS插值方法對27個站點多年平均降雨量的交叉驗證誤差，可見，這2種方法交叉驗證誤差具有極顯著的相關性，近85%的誤差集中在-30～30 mm之間。交叉驗證最大誤差站點均出現在興縣，且均大于40 mm。插值誤差較大的區域集中在山西河曲、五寨、興縣、離石、隰縣等站點分布較稀疏的地區。

3.2 年尺度插值結果比較

圖3 河龍區間1981—2000年多年平均降雨量插值結果Fig.3 Interpolation maps and semi-variogram cloud of the average annual precipitation of 1981—2000

圖4 KRIGING與TPS方法對河龍區間多年平均降雨量插值交叉驗證誤差對比Fig.4 Comparison of the cross-validation error for average annual precipitation interpolation between KRIGING and TPS methods at stations within He-Long Section

分別采用KRIGING與TPS方法對1981—2000年降雨量進行插值，對河龍區間內27個站點進行交叉驗證，計算平均E及A值，結果見圖5。結果表明，KRIGING和TPS方法對河龍區間1981—2000年降雨量插值結果的標準均方根誤差E(圖5(a))近95%的年份均介于40～70 mm之間，其均值分別為56.2和57.0 mm，一致性指標A值(圖5(b))均主要介于0.4～0.8之間，平均值分別為0.582和0.576。這2種方法的插值精度非常接近，平均E差值為多年均值的0.19%，平均A值相差0.006，KRIGING方法仍略微優于TPS方法。各年份2種方法A值相差基本在±0.05范圍內，與多年平均降雨量插值比較，2種方法的年插值精度平均A值均降低約0.19(19%)。

從圖5中選取E值較小、A值較大的1987年和E值較大、A值較小的1996年進行具體比較分析。1987年河龍區間暴雨次數較少，降雨強度較低，各站點降雨量差別較小，均在510 mm以下;而1996年一些站點降雨次數及暴雨強度較大，降雨量在各站點差別大，個別站點如興縣、安塞、清澗等地區降雨量甚至超過600 mm。研究區這2年的降雨量空間變異情況統計結果見表1?？芍?，盡管1987年各站點年降雨極值比、離差系數較1996年大，但1987年各站點年降雨量空間自相關性要較1996年好得多，也就是說，1996年各站點年降雨量空間變異性較1987年大得多。

圖5 河龍區間1981—2000年年均降雨量KRIGING和TPS插值方法的交叉驗證對比Fig.5 Cross-validation for annual precipitation using KRIGING and TPS methods

表1 河龍區間1987和1996年降雨量空間變異統計Tab.1 Spatial variation of annual precipitation in 1987 and 1996

對于不同空間變異的具體年份，KRIGING和TPS方法插值結果差異較大(圖6)。1987年(圖6(a))降雨量空間自相關性較好，KRIGING和TPS方法插值精度均較高，交叉驗證平均 A值分別為0.711和0.671，而1996年(圖6(b))降雨量的空間自相關性較差，2種方法插值的精度均較低，交叉驗證平均A值分別為0.407和0.471。與KRIGING方法相比，TPS方法的平滑作用可以改變局部精度而影響到全局效果。如1987年東南部插值的平滑作用，使全局交叉驗證A值提高了0.04，E值減少了6.9 mm;但1996年局部平滑使全局A值降低了0.06，E值增加了1.6 mm。

圖6 河龍區間1987、1996年降雨量KTIGING和TPS方法插值效果圖和半方差云圖比較Fig.6 Interpolation maps of annual precipitation in 1987 and 1996 in the He-Long Section using KRIGING and TPS methods

3.3 月尺度插值結果比較

河龍區間降雨量月際分配極不均勻，空間變異非常大。研究時段內，春季4月份約占區間多年平均降雨量的3% ～8%，平均降雨量10～35 mm，從南向北逐漸減少且梯度平緩。夏季7月份約占多年平均降雨量的15% ～30%，平均降雨量60～120 mm，從南向北逐漸減少，有多個強降雨中心。采用KRIGING和TPS方法對降雨較平穩的4月份，以及空間變異大、暴雨強度高的7月份的降雨數據進行27個氣象站點插值精度的交叉驗證，結果見圖7。

可知，對于4月份的插值結果(圖7(a)、(b))，KRIGING和TPS方法交叉驗證E值范圍均在2～18 mm之間，平均為6.2和6.6 mm。A值介于0.40～0.90之間，平均為0.629和0.631。對于7月份的插值結果(圖7(c)、(d))，KRIGING和TPS方法 E值范圍均在20～60 mm之間，平均為34.3和35.0 mm，A值介于0.10～0.70之間，平均為0.419和0.421。結果表明:無論降雨過程平穩的4月還是空間變異性較大的7月，KRIGING和TPS方法插值精度平均結果差別較小;這2種方法對月雨量插值交叉驗證平均一致性指標A值均相差0.002，但在具體年份有較大差別，2種方法A值絕對值差異基本在±0.05范圍內波動;在降雨變異較大的7月，2種方法插值的平均精度均低于年值插值結果約0.16(16%)，更低于多年平均插值結果約0.35(35%)，但對于降雨過程較平穩的4月，2種方法的插值精度均高于年值插值精度約0.05(5%)，低于多年平均插值結果約0.14(14%)。

由于7月降雨量插值平均精度較差，各年份間差異較大，故對插值精度較高的1992年7月和精度較低的1989、1996年7月插值表面、半方差云圖及站點交叉驗證誤差進行深入對比，其空間變異統計結果見表2，插值結果比較見圖8和圖9。

圖7 用KRIGING與TPS方法對河龍區間1981—2000年4和7月降雨量插值交叉驗證誤差對比Fig.7 Comparison of cross-validation in April and July precipitation using KRIGING and TPS methods

表2 河龍區間1989、1992、1996年7月降雨量空間變異統計Tab.2 Spatial variation of July precipitation in 1989，1992 and 1996

1992年(圖8(b))7月，要素在空間上呈正自相關，并且自相關性很顯著。2種方法插值表面光滑，TPS方法更為平滑，半方差隨著站點距離的增大而增大，TPS方法平滑參數大小適中，誤差判斷GCV較小(2.17)。由于有顯著的空間自相關性，所以，2種插值方法的精度較高，2種方法交叉驗證預測值與實測值一致性較好(圖9)，交叉驗證有效系數 ＞0.7。

而1989(圖8(a))和1996年(圖8(c))7月，要素空間自相關性很差，2種方法插值降雨量空間分布不一致，KRIGING方法插值表面不光滑。1989年，半方差甚至隨站點距離增大基本不變，導致KRIGING方法插值產生純塊金效應，交叉驗證預測值僅為相鄰一些站點的平均計算結果，與實測值沒有相關關系(圖9)，站點交叉驗證A＜0.2。TPS方法對1989和1996年插值平滑參數過大和過小，數據存在短相關，找不到合適的平滑表面，GCV較大，分別為3.01和5.50，交叉驗證A＜0.2。

總之，在1989和1996年的7月，由于空間上較低的相關性或呈現負的相關性，2種方法插值均不能正確反映降雨量的實際空間分布狀況，尤其在出現暴雨的站點，交叉驗證預測值與實際值可相差至100 mm以上，交叉驗證的預測值與實際值誤差較大，僅為各站點平均結果，基本沒有相關關系，2種方法插值結果均不能客觀地表現該時期的空間信息。對于這樣的空間分布狀況，引入協變量可能會改善其精度。

4 結論

圖8 河龍區間1989、1992和1996年7月降雨量KTIGING和TPS方法插值圖和半方差云圖比較Fig.8 Interpolation maps of July precipitation of 1989，1992 and 1996 in the He-Long Section using KRIGING and TPS methods

圖9 河龍區間1989、1992和1996年7月降雨量實測值與交叉驗證預測值相關關系Fig.9 Correlation between the predicted and actual values of July precipitation in 1989，1992 and 1996 at the stations within the He-Long Section

1)無論多年平均、還是年和月尺度，KRIGING和TPS方法插值結果都能正確反映河龍區間降雨量的空間分布趨勢，不同時間尺度上，2種方法的交叉驗證區域時段交叉驗證平均一致性指標A值相差均在±0.01范圍內。不同時間尺度2種方法的面平均插值精度均沒有顯著性和實質性差異。隨著尺度變小，降雨量空間分布變異性增大，2種方法插值精度均顯著降低，不同年份間一致性指標A值絕對值差別基本在±0.05范圍內波動。僅從數值上可以看出，多年平均和年尺度，KRIGING方法要略優于TPS法。

2)KRIGING和TPS方法在不同時間尺度插值精度均與全域要素自相關程度顯著相關。2種方法多年平均降雨量的插值精度均優于4月，其次是年值、7月。與多年平均降雨量插值平均一致性指標A值比較，4月插值精度均降低約14%，年值均降低約19%，而7月均降低約35%。

3)TPS方法插值降雨量的等值線更加平滑和柔和，但河龍區間降雨空間變異性較大，與KRIGING方法相比，TPS方法的局部平滑作用可以改善或削弱全局插值精度。

［1］蔡福，于慧波，矯玲玲，等.降水要素空間插值精度的比較［J］.資源科學，2006，28(6):73-79

［2］李新，程國棟，盧玲.空間內插方法比較［J］.地球科學進展，2000，15(3):260-265

［3］林忠輝，莫興國，李紅軒，等.中國陸地區域氣象要素的空間插值［J］.地理學報，2002，57(1):47-56

［4］王淑英，陳守煜.加權平均的權重優選算法及其應用［J］.水利學報，2003(12):1-7

［5］何紅艷，郭志華，肖文發.降水空間插值技術的研究進展［J］.生態學雜志，2005，24(10):1187-1191

［6］李軍龍，張劍，張叢，等.氣象要素空間插值方法的比較分析［J］.草業科學，2006，23(8):6-11

［7］劉志紅，Tim R McVicar，Li Lingtao，等.基于 5 變量局部薄盤光滑樣條函數的蒸發空間插值［J］.中國水土保持科學，2006，4(6):23-30

［8］許家琦，舒紅.降水數據空間插值的時間尺度效應［J］.測繪信息與工程，2009，34(3):29-30

［9］朱會義，賈紹鳳.降雨信息空間插值的不確定性分析［J］.地理科學進展，2004，23(2):34-42

［10］張漢雄.黃土高原的暴雨特性及其分布規律［J］.地理學報，1983，38(4):416-425

［11］Bargaoui Z K，Chebbi A.Comparison of two kriging interpolation methods applied to spatiotemporal rainfall［J］.Journal of Hydrology，2009，365(1-2):56-73

［12］Tiat A，Henderson R，Turner R，et al.Thin plate smoothing spline interpolation of daily rainfall for New Zealand using a climatological rainfall surface［J］.International Journal of Climatology，2006，26(14):2097-2115

［13］冉大川，柳林旺，趙力儀，等.黃河中游河口鎮至龍門區間水土保持與水沙變化［M］.鄭州:黃河水利出版社，2000

［14］王廣德，過常齡.“Krige”空間內插技術在地理學中的應用［J］.地理學報，1987，42(4):366-375

［15］劉永社，印興耀，賀維勝.空間相關分析因素對儲層建模中克里金估計結果的影響［J］.石油大學學報，2004，28(2):24-27

［16］劉志紅，Tim R McVicar，Li Lingtao，等.基于 ANUSPLIN的時間序列氣象要素空間插值［J］.西北農林科技大學學報，2008，36(10):227-234

［17］Hutchinson M F.ANUSPLIN Version 4.3 User Guide［M］.Canberra:The Australia National University，Center for Resource and Environment Studies，2004

［18］Hancock P A，Hutchinson M F.Spatial interpolation of large climate data sets using bivariate thin plate smoothing splines［J］.Environmental Modelling ＆ Software，2006，21(12):1684-1694

［19］Newlands N K，Davidson A，Howard A，et al.Validation and inter-comparison of three methodologies for interpolating daily precipitation and temperature across Canada［J］.Environmetrics，2010，22(2):205-223

［20］張曉萍，張櫓，穆興民，等.黃河中游河口一龍門區間多年平均流域水平衡特征［J］.地理學報，2007，62(7):753-763