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

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

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

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

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

1 研究區(qū)概況

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

2 數(shù)據(jù)資料及研究方法

2.1 數(shù)據(jù)資料

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

2.2 KRIGING插值法

KRIGING插值法是由法國地理學(xué)家Georges Martheron和南非采礦工程師D.G.Krige創(chuàng)立的地質(zhì)統(tǒng)計學(xué)的最佳內(nèi)插方法［14］，近年來被廣泛應(yīng)用于地學(xué)研究中。區(qū)域化變量和變異函數(shù)是克里金插值的基石。

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

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

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

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

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

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

擬合半方差重要的用途是確定局部內(nèi)插需要的權(quán)重因子，其值按采樣點數(shù)據(jù)的半方差圖的統(tǒng)計分析原理計算，即:

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

2.3 TPS方法

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

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

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

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

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

2.4 誤差分析方法

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

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

2.5 全域型Moran’s I指數(shù)

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

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

I值標(biāo)準(zhǔn)化Z(I)的公式為

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

3 結(jié)果與分析

3.1 多年平均降雨量插值結(jié)果比較

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

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

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

3.2 年尺度插值結(jié)果比較

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

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

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

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

表1 河龍區(qū)間1987和1996年降雨量空間變異統(tǒng)計Tab.1 Spatial variation of annual precipitation in 1987 and 1996

對于不同空間變異的具體年份，KRIGING和TPS方法插值結(jié)果差異較大(圖6)。1987年(圖6(a))降雨量空間自相關(guān)性較好，KRIGING和TPS方法插值精度均較高，交叉驗證平均 A值分別為0.711和0.671，而1996年(圖6(b))降雨量的空間自相關(guān)性較差，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 河龍區(qū)間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 月尺度插值結(jié)果比較

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

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

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

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

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

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

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

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

4 結(jié)論

圖8 河龍區(qū)間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 河龍區(qū)間1989、1992和1996年7月降雨量實測值與交叉驗證預(yù)測值相關(guān)關(guān)系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方法插值結(jié)果都能正確反映河龍區(qū)間降雨量的空間分布趨勢，不同時間尺度上，2種方法的交叉驗證區(qū)域時段交叉驗證平均一致性指標(biāo)A值相差均在±0.01范圍內(nèi)。不同時間尺度2種方法的面平均插值精度均沒有顯著性和實質(zhì)性差異。隨著尺度變小，降雨量空間分布變異性增大，2種方法插值精度均顯著降低，不同年份間一致性指標(biāo)A值絕對值差別基本在±0.05范圍內(nèi)波動。僅從數(shù)值上可以看出，多年平均和年尺度，KRIGING方法要略優(yōu)于TPS法。

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

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

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

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

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

［4］王淑英，陳守煜.加權(quán)平均的權(quán)重優(yōu)選算法及其應(yīng)用［J］.水利學(xué)報，2003(12):1-7

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

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

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

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

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

［10］張漢雄.黃土高原的暴雨特性及其分布規(guī)律［J］.地理學(xué)報，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］冉大川，柳林旺，趙力儀，等.黃河中游河口鎮(zhèn)至龍門區(qū)間水土保持與水沙變化［M］.鄭州:黃河水利出版社，2000

［14］王廣德，過常齡.“Krige”空間內(nèi)插技術(shù)在地理學(xué)中的應(yīng)用［J］.地理學(xué)報，1987，42(4):366-375

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

［16］劉志紅，Tim R McVicar，Li Lingtao，等.基于 ANUSPLIN的時間序列氣象要素空間插值［J］.西北農(nóng)林科技大學(xué)學(xué)報，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］張曉萍，張櫓，穆興民，等.黃河中游河口一龍門區(qū)間多年平均流域水平衡特征［J］.地理學(xué)報，2007，62(7):753-763