曹曉敏，劉志紅，張曉萍

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

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

1 資料來(lái)源與研究方法

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

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

圖1 黃土高原中游觀測(cè)站分布狀況

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

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

2.1 反距離加權(quán)插值

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

式中：Z（S0）——S0處的預(yù)測(cè)值；N——預(yù)測(cè)計(jì)算過(guò)程中要使用的預(yù)測(cè)點(diǎn)周圍樣點(diǎn)的數(shù)量；λi——預(yù)測(cè)計(jì)算過(guò)程中使用的各樣點(diǎn)的權(quán)重，該值隨著樣點(diǎn)與預(yù)測(cè)點(diǎn)之間距離的增加而減少；Z（Si）——在Si處獲得的測(cè)量值。確定權(quán)重的計(jì)算公式如式（2）。

式中：p——指數(shù)值；di0——點(diǎn)S0與各已知點(diǎn)Si之間的距離。

樣點(diǎn)在預(yù)測(cè)點(diǎn)值的計(jì)算過(guò)程中所占權(quán)重的大小受參數(shù)p的影響，也就是說(shuō)，隨采樣點(diǎn)與預(yù)測(cè)點(diǎn)之間距離的增加，標(biāo)準(zhǔn)樣點(diǎn)對(duì)預(yù)測(cè)點(diǎn)影響的權(quán)重按指數(shù)規(guī)律減少。

2.2 徑向基函數(shù)插值

徑向基函數(shù)插值法（RBF）如同將一個(gè)軟膜插入并經(jīng)過(guò)各個(gè)已知樣點(diǎn)，同時(shí)又使表面的總曲率最小，是精確插值方法。所謂精確插值方法就是指表面必須經(jīng)過(guò)每一個(gè)已知樣點(diǎn)。徑向基函數(shù)插值法適用于對(duì)大量點(diǎn)數(shù)據(jù)進(jìn)行插值計(jì)算，同時(shí)要求獲得平滑表面的情況。將徑向基函數(shù)應(yīng)用于表面變化平緩的表面，如表面上平緩的點(diǎn)高程插值，能得到令人滿意的結(jié)果。而在一段較短的水平距離內(nèi)，表面值發(fā)生較大的變化，或無(wú)法確定采樣點(diǎn)數(shù)據(jù)的準(zhǔn)確性，或采樣點(diǎn)數(shù)據(jù)具有很大的不確定性時(shí)，徑向基函數(shù)插值的方法并不適用。

2.3 普通克里格插值

地統(tǒng)計(jì)學(xué)就是針對(duì)區(qū)域化變量而發(fā)展的空間統(tǒng)計(jì)理論，主要研究那些分布于空間中并顯示出一定結(jié)構(gòu)性和隨機(jī)性的自然現(xiàn)象。它在考慮樣本點(diǎn)位置方向和彼此之間距離的基礎(chǔ)上，直接測(cè)定空間結(jié)構(gòu)的相關(guān)性和依賴性，研究具有一定隨機(jī)性和一定結(jié)構(gòu)性的各種變量的空間分布及變異規(guī)律。

普通克里格插值（Kriging）又稱空間局部插值法，是以變異函數(shù)理論和結(jié)構(gòu)分析為基礎(chǔ)，在有限區(qū)域內(nèi)對(duì)區(qū)域化變量進(jìn)行無(wú)偏最優(yōu)估計(jì)的一種方法，克里格方法考慮距離，而且通過(guò)變異函數(shù)和結(jié)構(gòu)分析，考慮了已知樣本點(diǎn)的空間分布及與未知樣點(diǎn)的空間方位關(guān)系，是地統(tǒng)計(jì)學(xué)的主要內(nèi)容之一。

普通克里格是區(qū)域化變量的線性估計(jì)，它假設(shè)數(shù)據(jù)變化成正態(tài)分布，認(rèn)為區(qū)域化變量Z的期望值是未知的。插值過(guò)程類似于加權(quán)滑動(dòng)平均，權(quán)重值的確定來(lái)自于空間數(shù)據(jù)分析?？杀硎救缡剑?）。

式中：Z（x0）——未知樣點(diǎn)的值；Z（xi）——未知樣點(diǎn)周圍的已知樣本點(diǎn)的值；λi——第i個(gè)已知樣本點(diǎn)對(duì)未知樣點(diǎn)的權(quán)重；n——已知樣本點(diǎn)的個(gè)數(shù)。

2.4 協(xié)同克里格插值

協(xié)同克里格法（CoKriging）是普通克里格的單個(gè)區(qū)域化變量向多個(gè)區(qū)域化變量的一種拓展，理論上并無(wú)本質(zhì)的區(qū)別，因此可以用推導(dǎo)普通克里格法的過(guò)程，推導(dǎo)協(xié)同克里格法。假設(shè)研究區(qū)域上有K個(gè)變量構(gòu)成協(xié)同區(qū)域化變Zk（x）＝（k＝1，2，…，k）滿足二階平穩(wěn)假設(shè)和本征假設(shè)，那么，可以確定交叉協(xié)方差函數(shù)和交叉變異函數(shù)存在，并定義如式（4）。

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

式中：aj和bj分別是兩個(gè)區(qū)域化變量的權(quán)重值。則整理后的協(xié)同克里格線性方程組表達(dá)式如式（6）。

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

3.1 站點(diǎn)數(shù)量與插值方法的插值精度比較

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

采用交叉驗(yàn)證法（cross－validation）來(lái)選擇氣象因子的最優(yōu)插值方法。用相對(duì)平均誤差（RME）作為驗(yàn)證幾種空間插值方法的標(biāo)準(zhǔn)。相對(duì)平均誤差是絕對(duì)平均誤差（RME）與站點(diǎn)實(shí)測(cè)值的平均值的百分比，見表1。

表1 不同站點(diǎn)數(shù)的相對(duì)平均誤差 %

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

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

3.2 反距離加權(quán)插值、徑向基函數(shù)插值、普通克里格插值比較

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

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

3.3 普通克里格插值、協(xié)同克里格插值比較

為了檢驗(yàn)插值精度以及比較兩種方法的優(yōu)劣，從58個(gè)樣本點(diǎn)中隨機(jī)抽取5個(gè)觀測(cè)站作為檢驗(yàn)點(diǎn)，進(jìn)行交叉驗(yàn)證，其中對(duì)多年平均降雨量的驗(yàn)證結(jié)果如表2所示，結(jié)果顯示插值精度較好。從總體上比較，普通克里格插值結(jié)果要好于協(xié)同克里格插值［3］。

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

圖3 反距離加權(quán)插值柵格圖

圖4 徑向基函數(shù)插值柵格圖

圖5 多年平均降雨量協(xié)同克里格柵格圖

圖5為協(xié)同克里格插值的結(jié)果，在總體空間格局上，圖2和圖5是相近的。協(xié)同克里格的空間插值中加入了高程因子，表現(xiàn)了因高程的影響而導(dǎo)致的局域變化。但是由于協(xié)同克里格考慮了高程的影響，各個(gè)氣象站點(diǎn)的海拔高度差異導(dǎo)致局域的變化明顯，空間變化相對(duì)突破單一的帶狀過(guò)渡，局域性、復(fù)雜性增強(qiáng)。

表2 多年平均降雨量?jī)煞N插值結(jié)果的交叉驗(yàn)證

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

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

4 結(jié)語(yǔ)

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

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

［1］ 朱會(huì)義，賈紹鳳.降雨信息空間插值的不確定性［J］.地球科學(xué)進(jìn)展，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］ 岳文澤，徐建華，徐麗華.基于地統(tǒng)計(jì)方法的氣候要素空間插值研究［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變量薄盤光滑樣條函數(shù)的區(qū)域蒸發(fā)量空間插值［J］.中國(guó)水土保持科學(xué)，2006，4（6）：23－30.

［8］ 劉志紅，劉文兆，李銳.基于3S技術(shù)的區(qū)域蒸散研究進(jìn)展［J］.中國(guó)水土保持科學(xué)，2006，4（3）：117－122.

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

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

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