基于GIS的黃土高原氣象要素空間插值方法*

孟慶香，劉國彬，楊勤科

(1.河南農(nóng)業(yè)大學資源與環(huán)境學院，鄭州450002;2.中國科學院水利部水土保持研究所陜西楊陵712100)

黃土高原位于中國西北部，包括山西、甘肅、寧夏全部和陜西、河南、內(nèi)蒙、青海部分地區(qū)，共計7省區(qū)287個縣，總面積 62.4萬 km2。溫度、降水等氣候因素是黃土高原農(nóng)牧業(yè)發(fā)展和生態(tài)環(huán)境惡化的重要制約因子，在氣候變化和人類活動的共同作用下，該區(qū)成為世界上水土流失較嚴重的地區(qū)之一。因此，研究氣候空間分布特征成為促進農(nóng)牧業(yè)發(fā)展和保護生態(tài)環(huán)境的必須，而獲取精確氣候數(shù)據(jù)的方法之一就是建立高密度的氣象站點。然而，由于經(jīng)濟、技術等原因，氣象站點是有限的，而定點觀測到的數(shù)據(jù)大多不能直接用于其他地點，更不能代替某一較大面積上的平均值。在實際工作中，研究人員將地統(tǒng)計學方法與地理信息系統(tǒng)相結合，通過對已知站點氣象數(shù)據(jù)進行空間插值，實現(xiàn)由點數(shù)據(jù)到面數(shù)據(jù)的轉化，生成研究區(qū)氣象要素的空間分布圖，則是一種有效的解決方法。目前，國外對此已進行了大量研究[1-5]，而國內(nèi)研究還較少。因此，本研究利用黃土高原已知的102個氣象站點氣候資料，進行區(qū)域年均溫的空間插值，探討不同插值方法對研究黃土高原氣候空間分布特征的影響，對黃土高原環(huán)境變化研究以及環(huán)境治理也有十分重要的指導作用。

1 數(shù)據(jù)與方法

1.1 數(shù)據(jù)基礎

氣象數(shù)據(jù)來源于國家氣象局102個氣象站點1950-2000年觀測值[6]。102個氣象站點的空間分布見圖1。以氣象站點的地理坐標和海拔高程作為自變量，由日降水量計算的多年年均降水量和由月均溫計算的多年年均氣溫用于插值。同時，以1∶50萬黃土高原行政區(qū)圖和1∶25萬黃土高原數(shù)字高程圖(DEM)為參照。

圖1 參與空間插值的102個氣象站點分布圖

1.2 軟件工具

軟件工具采用美國環(huán)境系統(tǒng)研究所(ERSI)的GIS桌面平臺系統(tǒng)ArcGIS 8.3以及目前常用的微軟電子表格工具Excel 2003。在ArcGIS 8.3的三大模塊中，ArcMap主要用于空間數(shù)據(jù)顯示、編輯、查詢、檢索、統(tǒng)計報表生成、空間分析和高級制圖等;ArcCatalog主要用于空間數(shù)據(jù)的瀏覽;ArcToolbox主要用于空間數(shù)據(jù)格式轉換、疊加處理、緩沖區(qū)生成和坐標轉換等[7]。電子表格工具Excel 2003主要用于原始數(shù)據(jù)的輸入和數(shù)據(jù)的簡單計算等。

1.3 數(shù)據(jù)處理轉換

上述研究成果多以微軟電子表格xls文件形式保存，由于目前ARCGIS尚不能很好地支持xls文件格式(必須加載 EXCEL_TOOLS工具才能實現(xiàn))，因此首先把含有各氣象站點編號、經(jīng)緯度、海拔高程和年均降水量的.xls文件保存為.dbf格式，然后通過ArcMap中的轉換工具將此文件生成黃土高原氣象站點圖。

為了與其他圖件(如行政區(qū)圖和DEM)進行疊加，上述數(shù)據(jù)均被統(tǒng)一到統(tǒng)一的坐標系和投影下。所采用的投影為等面積割圓錐投影(ALBERS)，中央經(jīng)線為東經(jīng)109.30°，雙標準緯線分別為北緯36.30°和北緯 37.10°，所采用的橢球體為 KRASOVSKY橢球體。

1.4 研究方法

目前，用于氣象要素空間插值的方法有很多，本研究主要采用反距離權重法、多項式插值法、徑向基函數(shù)插值法和克里金等方法。這里主要介紹反距離加權插值法、克里金法和協(xié)克里金法。

1.4.1 反距離加權插值法(IDW) 該法是20世紀60年代末提出的方法。IDW實際是一種加權移動平均方法，將計算區(qū)域劃分成若干矩形網(wǎng)格，每個網(wǎng)格的寬度和長度分別為Δx和Δ y，網(wǎng)格格點處的降水量xj可用其周圍鄰近的氣象站實測資料按距離平方的倒數(shù)插值求得，即:

式中:xj——第 j個格點處插得的降水量;pi——第j個格點鄰近的第i個氣象站點的實測降水量;di——第j個格點到其周圍鄰近的第i個氣象站點距離;m——第j個格點鄰近的氣象站點個數(shù)。反距離加權插值法的缺點是其計算值易受數(shù)據(jù)點集的影響，從而使計算結果常出現(xiàn)一種孤立點數(shù)據(jù)明顯高于周圍數(shù)據(jù)點的情況。

1.4.2 克里金方法(Kriging) 克里金方法是建立在地統(tǒng)計學基礎上的一種插值方法，最早由法國地理學家Matheron和南非礦山工程師Krige提出[8]，并用于礦山勘探?？死锝鸱椒ǔ浞治樟说乩斫y(tǒng)計的思想，其原理是假設某種屬性的空間變化既不是完全隨機也不是完全確定。反之，空間變化可能包括三個影響因素:空間相關因素，代表區(qū)域變量的變化;偏移或結構，代表趨勢;還有隨機誤差。偏移出現(xiàn)與否和對區(qū)域變量的解釋導致了用于空間插值的不同克里金方法的出現(xiàn)，主要有簡單克里金方法(Simple Kriging)和普通克里金方法(Ordinary Kriging)。簡單克里金方法假設不存在偏移，關注與空間相關因素。普通克里金方法假設除了已知點之間的空間關系外，空間變量在z值上還有偏移或有結構因素，表現(xiàn)為一個趨勢?？傮w而言，克里金方法認為，當空間變量的結構性成分確定后，剩余的差異變化屬于同質變化，不同位置之間的差異僅是距離的函數(shù)，可以表示為

式中:z(x0)——x0處的估計值;z(xi)——xi處的觀測值;λi——克里金權重系數(shù);n——觀測點個數(shù)。

1.4.3 引入高程信息的協(xié)克里金方法(Co-Kriging) Dirks等[3]發(fā)現(xiàn)，當站網(wǎng)密度較高時，普通克里金方法的插值效果與其他常用方法相比并無多大優(yōu)勢。Borga等[4]也曾得到同樣的結論。但是，隨著高程的增加，降水量有增加的趨勢。Hevesi等[8]研究了年平均降水量與高程的相關性，得到了其相關系數(shù)達到0.75的結果。因此，本文采用協(xié)克里金方法，并將高程作為第二影響因素引入到對降水量的空間插值中來。對多個具有空間相關性的空間變量進行估計的克里金方法可以歸類為協(xié)克里金方法。借助這類方法，可以利用幾個空間變量之間的相關性，對其中的一個變量或多個變量進行空間估計，以提高估計的精度和合理性。雖然協(xié)克里金方法的基本原理早已為人們所熟知，但至今仍不經(jīng)常使用，究其原因，主要在于其符號和算法比較復雜，交叉協(xié)方差函數(shù)和交叉變異函數(shù)的求取比較困難。

當研究區(qū)域內(nèi)的第2類信息(如高程)隨處可知且變化平穩(wěn)時，可將這類信息作為第2類影響因素引入?yún)f(xié)克里金方法?？紤]高程影響的協(xié)克里金估值方法可表示為

式中 :z(x0)，z(xi)，λi含義同式(2);y(x)— —x 點處的高程;my，mz——分別為高程和降水量的全局平均值;n——觀測點個數(shù)。

與克里金方法一樣，引入高程信息的協(xié)克里金方法也有簡單克里金方法(Simple Kriging)和普通克里金方法(Ordinary Kriging)之分。

2 插值結果與分析

2.1 不同插值方法比較

本研究采用克里金(Kiging)和協(xié)克里金(Co-Kriging)方法進行插值(選用Spherical模型)，每種方法分別采用簡單插值和普通插值。為了與常規(guī)插值方法進行比較，同時運用了反距離加權插值法(IDW)、全局多項式插值法(GPI)、局部多項式插值法(LPI)和徑向基函數(shù)插值法(RBF)進行插值。運用Cross-Validation進行交叉驗證，并對各參數(shù)進行修正，以期得到最好的插值結果[9-13]。

2.1.1 常規(guī)方法與地統(tǒng)計學方法的比較 將以上插值方法所得擬合值與實測值進行比較，計算誤差均值(MEAN)和誤差均方根(RMS)，結果見表1。一般來說，各種插值方法的MEAN和RMS總體最小者，具有較好的插值效果，尤其是RMS越小越好。從表1可知，對于降水量空間插值結果，誤差均方根排序為普通克里金插值法＜局部多項式插值法＜徑向基函數(shù)插值法＜普通協(xié)克里金插值法＜簡單克里金插值法＜反距離加權插值法＜簡單協(xié)克里金插值法＜全局多項式插值法;誤差均值排序為普通協(xié)克里金插值法＜全局多項式插值法＜簡單協(xié)克里金插值法＜普通克里金插值法＜反距離加權插值法＜徑向基函數(shù)插值法＜簡單克里金插值法＜局部多項式插值法。綜合來看，普通克里金和普通協(xié)克里金效果較好，進一步選擇插值方法，可以通過不同克里金插值的廣義交叉驗證得到。

同理，對于年均溫空間插值，誤差均方根以簡單協(xié)克里金插值法最小，普通協(xié)克里金插值法次之;且簡單協(xié)克里金插值法誤差均值最小，普通協(xié)克里金插值法次之。局部多項式插值法的誤差均值最大，徑向基函數(shù)插值法的誤差均值也大于普通協(xié)克里金和簡單協(xié)克里金。綜合來看，普通協(xié)克里金和簡單協(xié)克里金效果較好，尤其是簡單協(xié)克里金效果更好。從插值的MEAN和RMS總體最小看氣候插值結果，無論是降水量還是年均溫，地統(tǒng)計學方法均優(yōu)于常規(guī)方法。

表1 基于102個氣象站點不同插值方法檢驗比較的結果

2.1.2 不同克里金插值方法的比較 地統(tǒng)計學中通常采用以下5個指標評價其預測精度，即:預測誤差的均值(MEAN)、預測誤差的均方根(RMS)、平均預測標準差(ASE)、平均標準差(MS)和標準均方根預測誤差(RMSS)。其中，前4個指標越小越好，主要是RMS越小越好，第5個指標越大越好。

表2 基于102個氣象站點不同克里金插值方法檢驗比較的結果

通過廣義交叉驗證，從表2可以看出，對于年均溫空間插值，從絕對值看，預測誤差的均值排序為簡單克里金插值法＞普通克里金插值法＞普通協(xié)克里金插值法＞簡單協(xié)克里金插值法;預測誤差的均方根排序為簡單克里金＞普通克里金插值法＞普通協(xié)克里金插值法＞簡單協(xié)克里金插值法;平均預測標準差排序為簡單協(xié)克里金插值法＞簡單克里金＞普通克里金插值法＞普通協(xié)克里金插值法;平均標準差排序為簡單克里金＞普通克里金插值法＞普通協(xié)克里金插值法＞簡單協(xié)克里金插值法;標準均方根預測誤差排序為普通克里金插值法＞普通協(xié)克里金插值法＞簡單克里金＞簡單協(xié)克里金插值法。雖然，從標準均方根預測誤差看，普通克里金插值法較好，從平均預測標準差看，普通協(xié)克里金插值法較好，但是，從預測誤差的均值、預測誤差的均方根和平均標準差看，簡單協(xié)克里金插值法較好，尤其是簡單協(xié)克里金插值法預測誤差的均方根最小，綜合5項指標，選用102個氣象站點，對黃土高原年均溫進行插值，簡單協(xié)克里金插值法效果最好。其原因主要是簡單協(xié)克里金法將海拔高程作為第二影響因素引入到對年均溫的空間插值中來，利用地理位置、海拔高程和氣溫等空間變量的相關性，對溫度進行空間估計，大大提高了插值精度和合理性。

從表2可以看出，對于降水量空間插值結果，普通克里金與普通協(xié)克里金相比，RMS和ASE小，雖然MEAN和MS略大，但是RMSS也略大，因此，普通克里金法是降水量空間插值的最佳方法。同理，簡單協(xié)克里金法插值效果最好。

2.2 黃土高原氣象要素的空間分布

2.2.1 黃土高原降水量空間插值結果分析 通過以上分析，降水量采用普通克里金方法進行空間插值效果最好。黃土高原降水量空間分析就是采用基于102個氣象站點普通克里金方法插值實現(xiàn)點數(shù)據(jù)到面數(shù)據(jù)的轉換。用黃土高原行政區(qū)劃圖對102個站點插值得到的柵格圖進行截邊，得到黃土高原降水量等值線圖。黃土高原年均降水量的空間分布見附圖1。

由附圖1可知，黃土高原降水量偏少，多年平均降水量變化在117～721 mm。降水量分布存在明顯的地區(qū)差異和年季變化。黃土高原降水量總體分布呈現(xiàn)西北低、東南高的態(tài)勢，界限很明顯。內(nèi)蒙古的河套地區(qū)到寧夏銀川一線的西北部分降水量不足200 mm，屬于干旱地區(qū)，植被類型是草原化荒漠，植被稀疏。神木-靖邊-會寧一線以北的大部分地區(qū)降水為200～400 mm，屬于半干旱地區(qū)，其中，包頭-鹽池-同心一線以北為重半干旱區(qū)，以荒漠草原植被為主，以南屬于輕半干旱區(qū)，植被類型為典型草原。但青海部分地區(qū)比較復雜，年均降水量為400～500 mm。區(qū)域內(nèi)山西、河南以及西安-天水以北的陜甘地區(qū)降水量為400～600 mm，此線以南降水為600～800 mm，屬于半濕潤地區(qū)。400 mm不僅是半干旱-半濕潤氣候的分界線，亦是畜牧業(yè)區(qū)和農(nóng)耕區(qū)的分界線，該線以北適于牧業(yè)生產(chǎn)，以南適于農(nóng)業(yè)生產(chǎn)。

從降水量分布面積上看，在黃土高原，年均降水量小于200 mm的干旱地區(qū)占總面積的7.02%，年均降水量為200～400 mm的半干旱地區(qū)占總面積的35.46%，年均降水量為400～800 mm的半濕潤地區(qū)占總面積的57.52%。其中，在半濕潤地區(qū)，年均降水量400～600 mm占半濕潤區(qū)面積的96.87%，年均降水量600～800 mm占3.13%。黃土高原半濕潤和半干旱面積占整個區(qū)域面積的4/5以上，因此，黃土高原以半濕潤-半干旱氣候為主。

黃土高原天然降水量少，而年蒸散量多為700～1 000 mm，水分虧缺嚴重，不僅使農(nóng)業(yè)生產(chǎn)受到限制，而且常常又與氣候干旱互為因果關系，招致大風所導致的風蝕危害，使區(qū)域環(huán)境進一步惡化。

2.2.2 黃土高原年均溫空間插值結果分析 通過分析，年均溫空間插值采用簡單協(xié)克里金方法效果最好。黃土高原年均溫空間插值分析就是以高程為協(xié)變量的簡單協(xié)克里金方法。用黃土高原行政區(qū)劃圖對102個站點插值得到的柵格圖進行截邊，得到黃土高原年均溫等值線圖(見附圖2)。

黃土高原地區(qū)年均溫變化在1.0～14.0℃，由附圖2可知，黃土高原年均溫總體分布呈現(xiàn)西北低、東南高的態(tài)勢，界限很明顯。青海省西寧-同仁一線以西年均溫低于4.0℃，氣候寒冷，植被稀疏。永登-蘭州-隴西一線以西地區(qū)以及土默特左旗-原平-五臺以北年均溫為4.0～7.0℃，氣候也是相對寒冷，不利于植物生長。長武-志丹-陽泉一線以北的大部分地區(qū)溫度在7.0～10.0℃，很多植被可以生長，農(nóng)業(yè)牧業(yè)都有所發(fā)展。此線以南的陜西關中、山西中南部以及河南西北部年均溫在10.0～14.0℃，大部分植被生長茂盛，農(nóng)業(yè)發(fā)達[13-14]。

從年均溫分布面積上看，在黃土高原，年均溫小于4.0℃的寒冷干旱地區(qū)占總面積的3.20%，年均溫為 4.0～7.0℃的較寒冷地區(qū)占總面積的11.47%，年均溫為7.0～10.0℃的較溫暖地區(qū)占總面積的54.63%，年均溫10.0～14.0℃占總面積的30.71%。黃土高原7.0～14.0℃面積占整個區(qū)域面積的85%，因此，黃土高原氣候適宜農(nóng)牧業(yè)發(fā)展。

3 結論與討論

黃土高原降水量和年均溫空間插值表明，傳統(tǒng)的反距離權重插值、多項式插值和徑向基函數(shù)插值方法，插值誤差比地統(tǒng)計學中的克里金方法大。對比各種插值方法，對于降水量宜采用普通克里金法進行空間插值，年均溫宜采用簡單協(xié)克里金法。這是因為降水量受高程影響不大，而年均溫受高程影響較大，因此，將高程作為第2影響因素引入到空間插值中的協(xié)克里金方法，利用地理位置、高程和年均溫的空間變量相關性，進行空間估計，大大提高了插值精度和合理性。因此，不能籠統(tǒng)地說哪種插值方法最好，必須在特定條件下，針對特定區(qū)域的特定項目選擇最佳插值方法。

黃土高原降水量偏少，總體分布呈現(xiàn)西北低、東南高的態(tài)勢，界限明顯。降水量大于200 mm的半濕潤和半干旱面積占整個區(qū)域面積的4/5以上，因此，黃土高原以半濕潤-半干旱氣候為主。黃土高原年均溫較低，總體分布呈現(xiàn)西北低、東南高的態(tài)勢，界限也很明顯。年均溫7.0～14.0℃面積占整個區(qū)域面積的85%，黃土高原氣候適宜農(nóng)牧業(yè)發(fā)展。

數(shù)學模型與地理信息系統(tǒng)的集成將是第4代地理信息系統(tǒng)研究中的最重要問題之一，而空間插值模型是數(shù)學模型的主要方面，空間插值方法的應用與進一步研究有著廣闊的前景。目前用于插值的軟件還有很多，如 ANUSPLIN、PRISM 和 GIDS等[14]。對空間插值方法而言，沒有絕對最佳的方法，只有在特定的條件下，對于特定區(qū)域實際情況的最佳方法，在考慮插值誤差的同時，還要顧及到實用性，即是否具有簡便的插值檢驗方法，從而選出最優(yōu)插值方法。

[1] Mohamed A S.Reliabilty estimation of rainfall-runoff 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] Dirks K N，Hayl E，Stow C D，et al.High-resolution studies of rainfall on Norfolk Island.PartⅡ:interpolation of rainfall data[J].J.Hydrol，1998，208(3/4):187-193.

[4] Borga M，Vizzaccaro A.On the interpolation of hydrologic variables:formal equivalence of multiquadratic surface fitting and Kriging[J].J.Hydrol，1997，195(1/4):160-171.

[5] Hevesi J A，F(xiàn)lint A L，Isto J D.Precipitation estimation in mountainous terrain using multivariate geostatistics.part I:structural analysis[J].J.Appl.Meteor，1992，31:661-676.

[6] 中國科學院水利部水土保持研究所.黃土高原水土保持數(shù)據(jù)庫[EB/OL].[05-08-20]，http://www.loess.csdb.cn/.

[7] 黨安榮，賈海峰，易善楨，等.ArcGIS 8 Desktop地理信息系統(tǒng)應用指南[M].北京:清華大學出版社，2005.

[8] 侯景儒，黃競先.地質統(tǒng)計學的理論與方法[M].北京:地質出版社，1990.

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

[10] 任朝霞，楊達源.近50a西北干旱區(qū)氣候變化對農(nóng)業(yè)的影響[J].干旱區(qū)資源與環(huán)境，2007，21(8):48-52.

[11] 石朋，芮孝芳.降水空間插值方法的比較與改進[J].河海大學學報:自然科學版，2005，33(4):361-365.

[12] 王景雷，孫景生，張寄陽，等.基于 GIS和地統(tǒng)計學的作物需水量等值線圖[J].農(nóng)業(yè)工程學報，2004，20(5):51-54.

[13] 范敏杰，袁玉江，魏文壽，等.新疆伊犁地區(qū)夏季平均最高氣溫的重建和分析[J].干旱區(qū)研究，2008，25(1):75-81.

[14] 閻洪.薄板光順樣條插值與中國氣候空間模擬[J].地理科學，2004，24(4):163-168.