黃土高原中游降雨量的空間插值方法研究

曹曉敏，劉志紅，張曉萍

（1.青海省氣象局，西寧810001；2.成都信息工程學院 資源環境學院，成都610225；3.中國科學院 水利部 保持研究所，陜西 楊陵712100）

黃土高原位于中國中部偏北，黃河中游地區。總面積40萬km2。這一地區自然資源尤其是煤炭和天然氣資源豐富。在氣候變化和人類活動的共同作用下，該區已成為世界上水土流失較嚴重的地區之一。降水量作為干旱半干旱地區植被生長的主要限制性因子，直接決定著多沙粗沙區植被恢復重建成果，進而影響到對水土流失的防治作用。獲取精確的降雨量空間分布特征的方法之一是建立高密度的氣象站點。然而，由于經濟、技術等原因，氣象站點的數量是有限的，而定點觀測到的數據大多不能直接用于其他地點，更不能代替某一較大面積上的平均值。在實際工作中，研究人員將地統計學方法與地理信息系統相結合，通過對已知站點氣象數據進行空間插值，實現由點數據到面數據的轉化，生成研究區氣象要素的空間分布圖，則是一種有效的解決方法。因此，本研究利用黃土高原中游已知的氣象站點的降水資料進行區域降水量的空間插值，探討氣象站點密度與布局對插值結果的影響，為黃土高原環境變化研究及區域環境治理提供參考和指導。

1 資料來源與研究方法

數據來自黃土高原中游地區58個氣象臺（站），數據年限為28a（部分觀測站數據為20a），1981－2008年各氣象臺（站）28a的月降雨量值。用Arc－Map 9.2地理信息軟件將空間數據轉換成Albers等積投影，投影參數為：第1條緯線，25°N；第2條緯線，47°N；中央經線，110°E。圖1為黃土高原中游多沙粗沙區及周邊觀測臺（站）的分布狀況。反映黃土高原地形與氣象資料同步的數據字段包括：經緯度、高程。

目前，基于降雨量資料空間插值的方法很多，本研究主要采用反距離加權插值法、徑向基函數插值法、普通克里格法和協同克里格插值等方法。

圖1 黃土高原中游觀測站分布狀況

2 基于降雨量的各種插值方法研究

傳統氣候要素空間插值方法包括反距離加權法（Inverse Distance Weighting）、徑向基函數插值（Radial Basis Function Interpolation）、普通克里格插值（Kriging）等。

2.1 反距離加權插值

反距離加權插值法（IDW）是基于相近相似的原理：即兩個物體離得近，它們的性質就越相似。反之，離得越遠則相似性越小。它以插值點與樣本點間的距離為權重進行加權平均，離插值點越近的樣本點賦予的權重越大。其公式如式（1）。

式中：Z（S0）——S0處的預測值；N——預測計算過程中要使用的預測點周圍樣點的數量；λi——預測計算過程中使用的各樣點的權重，該值隨著樣點與預測點之間距離的增加而減少；Z（Si）——在Si處獲得的測量值。確定權重的計算公式如式（2）。

式中：p——指數值；di0——點S0與各已知點Si之間的距離。

樣點在預測點值的計算過程中所占權重的大小受參數p的影響，也就是說，隨采樣點與預測點之間距離的增加，標準樣點對預測點影響的權重按指數規律減少。

2.2 徑向基函數插值

徑向基函數插值法（RBF）如同將一個軟膜插入并經過各個已知樣點，同時又使表面的總曲率最小，是精確插值方法。所謂精確插值方法就是指表面必須經過每一個已知樣點。徑向基函數插值法適用于對大量點數據進行插值計算，同時要求獲得平滑表面的情況。將徑向基函數應用于表面變化平緩的表面，如表面上平緩的點高程插值，能得到令人滿意的結果。而在一段較短的水平距離內，表面值發生較大的變化，或無法確定采樣點數據的準確性，或采樣點數據具有很大的不確定性時，徑向基函數插值的方法并不適用。

2.3 普通克里格插值

地統計學就是針對區域化變量而發展的空間統計理論，主要研究那些分布于空間中并顯示出一定結構性和隨機性的自然現象。它在考慮樣本點位置方向和彼此之間距離的基礎上，直接測定空間結構的相關性和依賴性，研究具有一定隨機性和一定結構性的各種變量的空間分布及變異規律。

普通克里格插值（Kriging）又稱空間局部插值法，是以變異函數理論和結構分析為基礎，在有限區域內對區域化變量進行無偏最優估計的一種方法，克里格方法考慮距離，而且通過變異函數和結構分析，考慮了已知樣本點的空間分布及與未知樣點的空間方位關系，是地統計學的主要內容之一。

普通克里格是區域化變量的線性估計，它假設數據變化成正態分布，認為區域化變量Z的期望值是未知的。插值過程類似于加權滑動平均，權重值的確定來自于空間數據分析。可表示如式（3）。

式中：Z（x0）——未知樣點的值；Z（xi）——未知樣點周圍的已知樣本點的值；λi——第i個已知樣本點對未知樣點的權重；n——已知樣本點的個數。

2.4 協同克里格插值

協同克里格法（CoKriging）是普通克里格的單個區域化變量向多個區域化變量的一種拓展，理論上并無本質的區別，因此可以用推導普通克里格法的過程，推導協同克里格法。假設研究區域上有K個變量構成協同區域化變Zk（x）＝（k＝1，2，…，k）滿足二階平穩假設和本征假設，那么，可以確定交叉協方差函數和交叉變異函數存在，并定義如式（4）。

假設空間估計值Z＊x由兩個區域化變量Z（xi）和Y（xj）共同決定，下面給出兩個變量的協同克里格空間插值計算公式如式（5）。

式中：aj和bj分別是兩個區域化變量的權重值。則整理后的協同克里格線性方程組表達式如式（6）。

3 黃土高原中游多年平均降雨量的空間插值

3.1 站點數量與插值方法的插值精度比較

為了檢驗插值精度以及比較各種方法的優劣，選取相鄰的不同氣象站點數目進行反距離加權插值法、徑向基函數插值法和普通克里格插值法，選用3，5，10，15，20，25，30個站點數據，分別對年降雨量進行空間插值比較分析。從中選取基于降雨量的各種空間插值最合適的站點數目。

采用交叉驗證法（cross－validation）來選擇氣象因子的最優插值方法。用相對平均誤差（RME）作為驗證幾種空間插值方法的標準。相對平均誤差是絕對平均誤差（RME）與站點實測值的平均值的百分比，見表1。

表1 不同站點數的相對平均誤差 %

由表1看出：反距離加權插值法、徑向基函數插值法和普通克里格插值法對年降雨量適宜空間插值的氣象站點數分別為10個、15個和10個，其相對平均誤差分別為12.32%、11.32%和11.2%，對年降雨量插值的較優程度是OK＞RBF＞IDW。而普通克里格法對不同數目氣象站點的插值相對誤差較穩定。

從以上還可以得到：1）不同的插值方法對同一氣象要素適合空間插值的站點數目不一定相同；2）插值的站點數目不同，其插值結果的精度不同；該結論與朱會義等［1］的研究相一致。3）空間插值時，選擇合適的相鄰站點數目，插值結果的精度最高，而不是站點越多，插值精度越高［2］。該結論與朱會義等［1］的結論不同。

3.2 反距離加權插值、徑向基函數插值、普通克里格插值比較

根據半變異函數理論模型，利用ArcMap 9.2中的地統計學模塊對黃土高原中游1981－2000年20a年平均降雨量進行插值，得到的插值結果如圖2、圖3、圖4所示。

1km×1km的年降雨量的柵格空間數據庫圖2－4分別給出了不同方法插值出的年平均降雨量。從3幅圖中可看出，降雨量受地形的影響較大。特別是對位于長城西北側和東南側邊坡地區的陜西省的影響尤為明顯。而我們的樣本點，即各個氣象站點的海拔高度差異顯著。3種插值均能反映出我國黃土高原中游年降雨量的空間分布特征是從東南沿海到西北地區年降雨量逐漸減少，呈東高西低的梯度變化。

3.3 普通克里格插值、協同克里格插值比較

為了檢驗插值精度以及比較兩種方法的優劣，從58個樣本點中隨機抽取5個觀測站作為檢驗點，進行交叉驗證，其中對多年平均降雨量的驗證結果如表2所示，結果顯示插值精度較好。從總體上比較，普通克里格插值結果要好于協同克里格插值［3］。

圖2 普通克里格插值柵格圖

圖3 反距離加權插值柵格圖

圖4 徑向基函數插值柵格圖

圖5 多年平均降雨量協同克里格柵格圖

圖5為協同克里格插值的結果，在總體空間格局上，圖2和圖5是相近的。協同克里格的空間插值中加入了高程因子，表現了因高程的影響而導致的局域變化。但是由于協同克里格考慮了高程的影響，各個氣象站點的海拔高度差異導致局域的變化明顯，空間變化相對突破單一的帶狀過渡，局域性、復雜性增強。

表2 多年平均降雨量兩種插值結果的交叉驗證

由表2看出：通過克里格空間插值得到了黃土高原中游各個網格的多年平均降雨量數值。結果顯示：多年平均降雨量在空間上都呈現明顯的梯度變化，且空間變化幅度比較大。降雨量從東南部向西北部逐漸減少，從東南向西北逐漸增加。

上述兩種插值方法雖然可以反映氣候要素的總體分布格局，但插值精度并不高，因為是在黃土高原中游范圍內，只有58個樣本點，而且這些樣本點在空間分布上也是不均勻的。此外，相鄰的兩個樣本點，海拔相差達上百米。觀測設備以及觀測人員水平的參差不齊，觀測精度也不一致。所以在小范圍內的插值精度仍然較差，特別是樣本點較少或缺少樣本點的地區，例如西北端的地區。所以除了精確的插值方法外，合理的空間抽樣和一致的觀測精度對精確插值也是必不可少的。

4 結語

（1）基于地統計的插值方法，根據實驗方差最小的原理，選擇合適的半變異函數理論模型，進行空間插值，能夠很好地模擬區域化變量的空間連續分布格局。克里格插值由于充分考慮了區域化變量的特性，經過檢驗發現較傳統方法，其插值精度大大提高。另外，對比普通克里格法和協同克里格法，后者增加了高程對降水量的影響因素，理論上說在空間上更為合理，但此次插值的精度表現為普通克里格法要略高于協同克里格法。

（2）采用地統計方法雖然在總體上能夠較好地反映氣候要素的空間分布格局，但是由于理論模型的擬合優度、研究區域和數據（研究區域范圍過大，地形變化過于復雜、形狀不規則、樣本數量有限、空間分布不均勻、樣本測量精度）等問題，導致上述兩種方法空間插值的精度還不高。所以進行地統計插值結果檢驗時，除了模型本身外，要充分考慮研究區域的特征，還可以通過增加樣本數量，提高觀測精度以及改進半變異函數理論模型等方式進一步提高插值精度。

［1］ 朱會義，賈紹鳳.降雨信息空間插值的不確定性［J］.地球科學進展，2004，23（3）：34－21.

［2］ Booth T H.Mapping regions climatically suitablefor particular tree species at the global scale［J］.Forest E－cology and Management，1990，36：172601.

［3］ 岳文澤，徐建華，徐麗華.基于地統計方法的氣候要素空間插值研究［J］.高原氣象，2005，24（6）：974－979.

［4］ Lamn.Spatial interpolation methods：a review［J］.The Amercian Cartographer，1983，10（2）：129－149.

［5］ Mat heron G.Principles of Geostatistics［J］.Economic Geoogy，1963，58：1246－1266.

［6］ Issaks E H R M Srivastava.An introduction to applied geostatistics［M］.New York：Oxford Univ.Press，1989.

［7］ 劉志紅，Tim R，McVicar，等.基于5變量薄盤光滑樣條函數的區域蒸發量空間插值［J］.中國水土保持科學，2006，4（6）：23－30.

［8］ 劉志紅，劉文兆，李銳.基于3S技術的區域蒸散研究進展［J］.中國水土保持科學，2006，4（3）：117－122.

［9］ 劉志紅，Tim等.專用氣候數據空間插值軟件ANUSPLIN及其應用［J］.氣象，2008，34（2）：18－23.

［10］ 劉志紅，Tim K，Mcvicar，等.基于ANUSPLIN的時間序列氣象要素空間插值［J］.西北農林科技大學學報：自然科學版，2008，36（6）：227－234..

［11］ 汪幫穩，楊勤科，劉志紅.基于DEM和GIS的修正通用土壤流失方程地形因子值的提取［J］.中國水土保持科學，2007，5（2）：18－23.