劉益良，殷堯，朱慶杰，劉小光

(河北聯(lián)合大學(xué)河北省地震工程研究中心，河北唐山 063009)

地統(tǒng)計學(xué)運用概率論和數(shù)理統(tǒng)計方法研究地理數(shù)據(jù)的空間分布規(guī)律，已經(jīng)形成一個成熟的學(xué)科技術(shù)體系。而把地統(tǒng)計學(xué)與地理信息系統(tǒng)(GIS)整合在一起，構(gòu)建一套完整的地統(tǒng)計分析工具，用ESRI的話說，是一次革命［1］。因為只有通過計算預(yù)測模型的統(tǒng)計誤差，GIS才能夠定量描述地理空間表面模型的質(zhì)量。隨著計算機技術(shù)與地質(zhì)統(tǒng)計方法的迅速發(fā)展，一些包含地質(zhì)統(tǒng)計模塊的地理信息系統(tǒng)軟件相繼開發(fā)并得到了廣泛的應(yīng)用。有兩個地理信息系統(tǒng)軟件進行了這樣的嘗試，它們是ArcGIS和IDRISI。

隨著我國經(jīng)濟的快速發(fā)展、人口的急劇膨脹以及生態(tài)環(huán)境的變化，使得社會對水資源的開發(fā)利用等一系列的問題更加關(guān)注，我國是水資源大國，對于水資源的研究顯得尤為重要。降雨是水資源的重要來源之一，對工農(nóng)業(yè)的生產(chǎn)建設(shè)和人們的日常生活都有著重大的意義。所謂降雨量，即從天空降落到地面上未經(jīng)蒸發(fā)、滲透及流失的降水在水平面上聚集的水層深度。降雨具有間斷性和空間變異性的特點，本文應(yīng)用IDRISI軟件的地統(tǒng)計模塊對華北地區(qū)的48個氣象站1971年至2000年的年降雨量進行空間變異性分析，研究華北地區(qū)降雨的空間分布規(guī)律。

1 地統(tǒng)計與空間變異性

1.1 地統(tǒng)計

地統(tǒng)計(Geostatistics)分析作為空間統(tǒng)計分析的一個主要內(nèi)容，首先是由G.Matheron于1962提出的，是以區(qū)域化變量理論為基礎(chǔ)，以變異函數(shù)為主要工具，研究那些在空間分布上既有隨機性又有結(jié)構(gòu)性，或空間相關(guān)性和依賴性的科學(xué)。地質(zhì)統(tǒng)計學(xué)區(qū)別于經(jīng)典統(tǒng)計學(xué)的最大特點即是:地質(zhì)統(tǒng)計學(xué)既考慮到樣本值的大小，又重視樣本空間位置及樣本間的距離，彌補了經(jīng)典統(tǒng)計學(xué)忽略空間方位的缺陷。地統(tǒng)計分析理論基礎(chǔ)包括前提假設(shè)、區(qū)域化變量、變異分析和空間估值［2］。

1.2 空間變異性

在統(tǒng)計學(xué)中，協(xié)方差函數(shù)和相關(guān)系數(shù)是測量兩個隨機變量相似性和差異性的常用工具。在地理學(xué)中，一個地理要素的觀測值在兩個較近位置間的差比較遠位置間要小，地統(tǒng)計學(xué)應(yīng)用了這兩個領(lǐng)域的理論和方法，用數(shù)學(xué)工具來計算和表達地理空間的相似性和差異性，即空間依賴性?？臻g相似性即地理上的連續(xù)性(spatial continuity)，是地理現(xiàn)象或特征沿一定方向和距離，其屬性以及表達屬性的測量值沒有或很少變化的一種分布模式。差異性則相反，它在數(shù)學(xué)上的表達，特稱空間變異性(spatial variability)?？臻g變異指研究對象在空間上的變化，它是地質(zhì)統(tǒng)計學(xué)研究的基本問題。

空間連續(xù)性或空間變化性分布現(xiàn)象的地統(tǒng)計學(xué)描述方法有很多，半方差函數(shù)運算(semi variogram calculation)是使用最多的一類，它計算采樣數(shù)據(jù)成員間在指定的距離間隔和方向上平均的差異程度［3］。

2 克里金法

上世紀50年代初期，南非地質(zhì)學(xué)家、礦業(yè)工程師克里金(D.G.Krige)［4］在礦山工作時觀察到金屬的分布在空間上并非是純隨機的，而是在空間上具有相互聯(lián)系性。他于1952年提出克里金法，是一種無偏的最小誤差的儲量計算方法。該方法按照樣品與待估塊段的相對空間位置和相關(guān)程度來計算塊段品位及儲量，并使估計誤差為最小。

法國巴黎礦業(yè)學(xué)院馬特隆教授［5］(G.Matheron)對克里金提出的方法進行研究，認為克里金法是在考慮了礦化空間分布特征的基礎(chǔ)上，合理地改進了統(tǒng)計學(xué)，是一種傳統(tǒng)儲量計算方法與統(tǒng)計學(xué)方法高度結(jié)合起來的新方法。同時為了解決具有二重性(結(jié)構(gòu)性與隨機性)的地質(zhì)變量的條件下使用統(tǒng)計方法的問題，馬特隆教授于1962年在經(jīng)過深入探討和實踐之后，提出了區(qū)域化變量(Regionalized Variable)的概念，從而創(chuàng)立了地質(zhì)統(tǒng)計學(xué)(Geostatistics)。根據(jù)地質(zhì)統(tǒng)計學(xué)理論，地質(zhì)特征可以用區(qū)域化變量的空間分布特征來表征。而研究區(qū)域化變量的空間分布特征的主要數(shù)學(xué)工具是變差函數(shù)(Variogram)或半方差函數(shù)(Semivariogram)，二者只是在定義上相差常數(shù)系數(shù)1/2，并無本質(zhì)區(qū)別。

2.1 半方差函數(shù)

當(dāng)定義了采樣數(shù)據(jù)對(xi，xi+h)是在指定的距離間隔lag內(nèi)的數(shù)據(jù)點對集合中空間距離為h的樣本，z(x)，z(x+h)是隨機變量z在x，x+h處的采樣點值，同時假設(shè)兩點之間的方差只與距離h有關(guān)，區(qū)域性變量理論的兩個內(nèi)在假設(shè)條件是差異的穩(wěn)定性和可變性，一旦結(jié)構(gòu)性成分確定后［6］，剩余的差異變化屬于同質(zhì)變化，不同位置間的差異僅是距離的函數(shù)，定義半方差函數(shù)為

式中，n為距離為h的樣本對的數(shù)目，采樣間隔h也叫步長。半變差函數(shù)(semivariogram)是IDRISI中采用的方差函數(shù)運算(variogram calculation)的一種統(tǒng)計量，屬于空間變化性統(tǒng)計描述方法。對應(yīng)于h的γ(h)的圖被稱為半方差圖。

2.2 克里金方程

式中，λi是賦予觀測值z(xi)的權(quán)重，表示各個觀測值z(xi)對估計值的貢獻。克里金是根據(jù)無偏估計和方差最小來確定加權(quán)系數(shù)λi的，即

3 華北降雨量變異性分析

3.1 華北地區(qū)概況

我國華北地區(qū)的地理范圍有多種說法，而且現(xiàn)在還沒有對華北地區(qū)進行明確的定義。本文所指的華北地區(qū)僅是從氣象預(yù)報的角度考慮，是指位于中國北部的區(qū)域，中國蒙古邊境線以南，秦嶺淮河線以北，西鄰青藏高原，東瀕渤海與東北地區(qū)的中國廣大區(qū)域。位于北緯34°～53°，東經(jīng)97°～125°之間，包括北京市、天津市、河北省、山西省和內(nèi)蒙古自治區(qū)。

華北地區(qū)地勢西北高、東南低，山地及山前沖、洪積平原是兩大基本地貌單元。所跨的經(jīng)緯度比較廣，地貌多樣，地形地勢復(fù)雜，區(qū)域氣候空間變異性明顯。北部內(nèi)蒙古屬典型的中溫帶季風(fēng)氣候，具有降水量少而不勻、寒暑變化劇烈的顯著特點。冬季漫長而寒冷，多數(shù)地區(qū)冷季長達5個月到半年之久。其中1月份最冷，月平均氣溫從南向北由零下10℃遞減到零下32℃，夏季溫?zé)岫虝?，多?shù)地區(qū)僅有一至兩個月，部分地區(qū)無夏季。最熱月份在7月，月平均氣溫在16℃至27℃之間，最高氣溫為36℃至43℃。氣溫變化劇烈，冷暖懸殊甚大。降水量受地形和海洋遠近的影響，自東向西由500毫米遞減為50毫米左右。蒸發(fā)量則相反，自西向東由3000毫米遞減到1000毫米左右。與之相應(yīng)的氣候帶呈帶狀分布，從東向西由濕潤、半濕潤區(qū)逐步過渡到半干旱、干旱區(qū)。降水量空間變異性是根據(jù)實際觀測的降水資料，分析降水量的空間變異性。

3.2 數(shù)據(jù)處理

選取華北地區(qū)48個氣象站從1971年至2000年30年的年均降雨量進行分析［7-8］，首先將這些氣象站的名稱、經(jīng)度、緯度、高程、降雨量等數(shù)據(jù)存入Excel文檔，然后將這些數(shù)據(jù)導(dǎo)入IDRISI數(shù)據(jù)庫，如圖1所示。將數(shù)據(jù)庫中的數(shù)據(jù)生成矢量數(shù)據(jù)文件rain.vct，如圖2所示。圖中的小圓點越大，表示此氣象站的降雨量最大。使用屬性特征查詢工具可以來查詢每個降雨量觀測站的屬性數(shù)值。

2.3 半方差表面圖

半方差表面圖代表了空間統(tǒng)計［9］，它是依據(jù)方差云圖(variogram cloud)得到的。方差云圖是將每個采樣數(shù)據(jù)點與其他采樣數(shù)據(jù)點匹配并得到每個結(jié)果對的變差函數(shù)值制圖的結(jié)果。根據(jù)其分離矢量，即分離的距離和分離方向定位半方差函數(shù)值，然后顯示結(jié)果。在云圖上疊加?xùn)鸥瘢骄總€單元的云值創(chuàng)建一個柵格半方差表面圖。在IDRISI軟件地統(tǒng)計模塊中打開空間依賴性建模器，將數(shù)據(jù)處理時生成的rain作為輸入矢量變量(Variable)文件，顯示類型應(yīng)設(shè)置為表面類型，生成半方差表面圖。由于采樣點數(shù)比較少，半方差表面圖不能明顯看出降雨的連續(xù)性，改變步長參數(shù)，將步長數(shù)目設(shè)置為15，步長寬度設(shè)置為0.5，再次生成表面圖，如圖3左半部分。

柵格中央的步長距離為零，從中央沿各個方向向外移動時步長距離增加，方向讀度數(shù)從正北開始按順時針計算。當(dāng)使用IDRISI標準調(diào)色板，深藍色的顏色代表低的變異函數(shù)值，或低的變化，而綠色的顏色代表的高可變性。請注意在表面圖底部顯示的是與表面圖中心的方向以及步長的數(shù)目。從圖3可以看出，一個深藍色帶狀大約呈37度分布，即最大的連續(xù)性方向在37度左右。

2.4 方向方差圖

鑒于樣本數(shù)據(jù)比較少，為了更好地分析降雨量的變異性，我們需要從不同的角度考慮。現(xiàn)在，我們將改進我們的分析方法，使用方向性圖向不同的方向伸長。將顯示類型有表面改變?yōu)榉较蛐?，方向角度設(shè)置為37度，繪制方向方差圖。然后將方向角度改為127度，繪制127度方向方差圖。選取全方位選項，繪制全方位圖，如圖3右半部分。

2.5 空間變異性分析

IDRISI探索空間變異性主要依靠觀察半方差表面圖(semivariogram surface)［3］。通過輸入不同的步長寬度(lag width)和步長數(shù)(number of lags)，半方差表面圖上清晰地顯示出沿某些方向表示方差值的顏色的均一性，即可判定數(shù)據(jù)中存在方向性差異，IDRISI的點數(shù)據(jù)圖層實現(xiàn)分級符號和分級彩色結(jié)合顯示，結(jié)合圖層的點位符號色彩分布和方差表面圖來確定數(shù)據(jù)方向相關(guān)性的長軸方向(majorrange)和短軸方向(minor range)。同時在方向方差圖上，可同時顯示4個方向的方差點系列。本文只繪制了37度，127度和全方位的這三個方向的系列。

從圖3可以注意到，37度系列隨著間隔距離的增加連續(xù)性變化最低。在與其正交的127度方向，使用相同的步長間距，變異性卻增加迅速。全方位的系列類似是所有方向的平均，因此，在這個例子中它位于兩個系列之間。表面圖和方向圖的比較是合乎邏輯的。37度和127度系列揭示了在方向上的差異程度，即各向異性程度。

3 結(jié)論

通過上述分析可以得出華北地區(qū)降雨量的分布規(guī)律，在地圖上以正北方向為0度，按順時針計算，在127度方向上降雨量變異性最大。在同一條37度斜線附近氣象站的降雨量變化很小。從圖2可以看出，從東南沿海的天津、滄州、樂亭一帶，越向西北內(nèi)陸地區(qū)降雨量越少。使用軟件分析得到的結(jié)果與實際情況相似，可見使用IDRISI軟件可以簡便、快捷、準確地對某一數(shù)據(jù)進行空間變異性分析，省去了大量的手工計算工作。我們此次應(yīng)用IDRISI軟件對華北地區(qū)1971年至2000年年均降雨量進行變異性分異，是對雨量變異規(guī)律探討的一個新的嘗試，從而使雨量變異規(guī)律的研究有了進一步的深入。本文僅使用了空間依賴性建模進行空間變異性分析，今后還可以使用模型擬合模塊可視化地模擬一模型來描述數(shù)據(jù)的空間變化特征，使用克里金插值形成一連續(xù)表面。

［1］David Hollingsworth.The Workflow Reference Mode 1 Issue 1.1［M］.The Workflow Management Coalition Specification，1995-01.

［2］靜，何必，李海濤.ArcGIS 9.3 Desktop地理信息系統(tǒng)應(yīng)用教程［M］.北京:清華大學(xué)出版社，2011，02.

［3］王茯泉.地統(tǒng)計分析在ArcGIS和IDRISI中實現(xiàn)特點的討論［J］.計算機工程與應(yīng)用，2006，42(15):210-215.

［4］Journel A G，Huijbregts C J.Mining Geostatistics［M］.Academic Press，New York.1978.

［5］孫洪泉.地質(zhì)統(tǒng)計學(xué)及其應(yīng)用［M］.北京:中國礦業(yè)大學(xué)出版社，1990.

［6］鄔倫.地理信息系統(tǒng)原理、方法和應(yīng)用［M］.北京:科學(xué)出版社，2001.

［7］中國氣象局.中國氣象科學(xué)數(shù)據(jù)共享服務(wù)網(wǎng):中國地面國際交換站氣候標準值年值數(shù)據(jù)集(1971-2000).［EB/OL］.［2012-3-20］.http://www.cma.gov.cn/2011qxfw/2011qsjgx/index.htm.

［8］天津政務(wù)網(wǎng).天津史志［EB/OL］.［2012-3-20］.http://www.tj.gov.cn/tblm/tj_sz/200709/t20 070921_22670.htm.

［9］劉光，賀小飛.地理信息系統(tǒng)實習(xí)教程［M］.北京:清華大學(xué)出版社，2003，01.