方曉剛，李宇鵬

（中國石油大學（北京），北京102249）

地質(zhì)統(tǒng)計學發(fā)源于礦業(yè)領域，并日益擴展至石油、環(huán)境、水文、農(nóng)林等領域［1－5］。理論上可以對于任何空間上即具有隨機性又具有結(jié)構(gòu)性的自然現(xiàn)象進行研究。而變差函數(shù)是經(jīng)典地質(zhì)統(tǒng)計學中一個重要工具，它能夠反映地質(zhì)變量的空間變化特征—相關性和隨機性。因為其能透過隨機性反應空間變量的結(jié)構(gòu)性，因此也稱為結(jié)構(gòu)函數(shù)，進行變差函數(shù)分析，有時也稱為空間變量的結(jié)構(gòu)分析［2］。

目前很多空間數(shù)據(jù)分析及繪圖商業(yè)軟件如Surfer，MapGIS，……等以及礦業(yè)儲量評估軟件等等均具有地質(zhì)統(tǒng)計學模塊?？梢允褂盟鼈冞M行變差函數(shù)分析以及變量空間插值等。雖然這些這些商業(yè)軟件提供了不同的分析界面，成果的可視化也不盡相同，但是如此眾多的商業(yè)軟件所使用的方法都是建立現(xiàn)有的成熟的算法之上。而往往這些算法原理都是一致的，且都是公開發(fā)表的，并且都有隨之而來的開源程序。使用這些開源的程序不僅可以節(jié)約花費在使用商業(yè)軟件的培訓和維護費用，而且基于開源程序的優(yōu)點是分析的結(jié)果過程透明，參數(shù)可控，原理明了，所以可以根據(jù)具體數(shù)據(jù)進行特殊的研究，滿足具體礦山開發(fā)的需要。

文章對地質(zhì)統(tǒng)計分析中最基本的過程：變差函數(shù)的分析過程進行了詳盡說明。這一個過程包括對數(shù)據(jù)的檢查、數(shù)據(jù)正態(tài)轉(zhuǎn)換、實驗變差函數(shù)的求取，以及理論變差函數(shù)的人工及自動擬合等。文中對每一步驟的重要性，以及所使用的開源程序作為論述的重點，而對每一步的理論基礎沒有做過多的闡述，指出一定的參考文獻，感興趣的讀者可以進一步的學習。

1 程序與數(shù)據(jù)體介紹

在眾多的公開發(fā)布的程序中，應用最廣泛的當屬美國Stanford大學的開發(fā)的Gslib地質(zhì)統(tǒng)計學程序庫（可以在http：//www.gslib.com免費下載）［6］。程序可以在任何操作系統(tǒng)下順利運行。并且有源程序代碼（Fortran 90）下載，便于感興趣進行更深入的學習和應用。其它的開源軟件則有GeoR，Gstat等［7－8］。

使用Gslib時，所有程序運行所讀入的數(shù)據(jù)文件應符合Geo－eas的格式［6］。即文本的第一行為注釋行，給出簡單的數(shù)據(jù)說明。第二行為數(shù)據(jù)列數(shù)nvar（number of variables），即共有多少列數(shù)據(jù)。由第三行開始則為這nvar列數(shù)據(jù)的名稱，每一行為一列數(shù)據(jù)的名稱。之后為nvar列數(shù)據(jù)。圖1中給出了一個數(shù)據(jù)文件的前9行。

圖1 數(shù)據(jù)文件例子

所有Gslib程序的運行都采用了參數(shù)文件的方式。即程序運行后要求輸入?yún)?shù)文本，運算所需要的設置都通過參數(shù)文件提供給程序。程序初次運行時會生成一個默認的參數(shù)文件，使用者可以根據(jù)自己的實際情況修改樣本參數(shù)文件，重新運行程序即可。程序運行時所生成的圖形文件都保存為postscript文檔格式。推薦使用GSview程序打開（Http：//pages.cs.wisc.edu/～ghost/gsview/免費下載）。

地質(zhì)統(tǒng)計學起源于礦業(yè)，最基本的應用仍然是礦產(chǎn)資源的品位估計和儲量計算。文章所使用的數(shù)據(jù)是由Pierre Mousset Jones提供給地質(zhì)統(tǒng)計學研究領域的一套真實的Brenda礦場3D數(shù)據(jù)體［9］。Brenda礦場是一個銅/鉬礦體，位于British Columbia（位于Central Okanagan的Peachland以西大概22km處）。數(shù)據(jù)包括各點x、y、z坐標，銅鉬的百分含量，以及取樣的長度和頂?shù)咨疃取Ｎ闹袃H選用了銅含量的數(shù)據(jù)進行說明。

2 變差函數(shù)分析

2.1 數(shù)據(jù)收集與檢查（Exploratory data analysis）

進行變差函數(shù)分析的第一步就是對原始數(shù)據(jù)的探索性分析（EDA：Exploratory data analysis）。可以說EDA本身并不能算一個方法或者技術，它主要是一種我們在進行地質(zhì)統(tǒng)計學分析時所應持有的態(tài)度，即基于數(shù)據(jù)的客觀態(tài)度。在對于礦業(yè)原始數(shù)據(jù)的EDA主要包括對原始數(shù)據(jù)空間位置分布，均值和方差統(tǒng)計，以及分布形態(tài)檢查等。通過這些簡單的統(tǒng)計可以對所取得的數(shù)據(jù)獲得一個大的了解，便于開展后期工作。

經(jīng)過簡單的直方圖繪制，可以得出數(shù)據(jù)的均值，標準差以及其他的一下分布形態(tài)參數(shù)。繪制直方圖可以使用GSLIB程序中Histplt程序。圖2是此數(shù)據(jù)集中銅鉬原始數(shù)據(jù)的直方圖。不僅可以看出其分布比較接近與對數(shù)分布，而且還可以識別出一些異常值。例如在此數(shù)據(jù)集中，銅的百分含量就出現(xiàn)了1.024這樣的一個數(shù)據(jù)點。就有必要對此點數(shù)據(jù)進行檢查。實際工作中如果有條件的話，應該重新測量。在這里，可以簡單的將其設置為全區(qū)的均值0.664。重新繪制直方圖如圖2所示。

2.2 數(shù)據(jù)正態(tài)轉(zhuǎn)換（Normal score data transformation）

經(jīng)過初步的檢查認為數(shù)據(jù)中人為的錯誤基本上消除之后就可以開始下一步的工作。由于后期的克里格估計（kriging estimation）以及隨機模擬都假設數(shù)據(jù)服從正態(tài)分布［1－2］，所以必須對數(shù)據(jù)進行正態(tài)轉(zhuǎn)換，而且正態(tài)轉(zhuǎn)換后有利于變差函數(shù)擬合步驟中基臺值的設定。

對原始數(shù)據(jù)的正態(tài)轉(zhuǎn)換，可以使用Gslib中的Nscore程序，。在此轉(zhuǎn)換中是使用了基于原始數(shù)據(jù)概率分布的分位數(shù)轉(zhuǎn)換，最大限度的保證了原始數(shù)據(jù)的分布形態(tài)［6］。使用分位數(shù)進行正態(tài)轉(zhuǎn)換一般要求樣本點的個數(shù)比較多。當點對個數(shù)較少時（少于50個）可以使用Histsmth程序?qū)υ紨?shù)據(jù)進行光滑，之后再進行正態(tài)轉(zhuǎn)換比較可靠［6］。在一般的礦業(yè)變差函數(shù)分析中，原始數(shù)據(jù)的個數(shù)應該是可以滿足使用分位數(shù)的Nscore轉(zhuǎn)換。將數(shù)據(jù)由正態(tài)恢復至原始分布的轉(zhuǎn)換則是使用Backtr程序進行。如果是使用Gslib中的Sgsim做隨機模擬的話，程序提供了一個是否自動轉(zhuǎn)換的控制參數(shù)。如果想在Sgsim模擬結(jié)果直接顯示為原始數(shù)據(jù)，則需要提供nscore程序運行時產(chǎn)生的轉(zhuǎn)換分位數(shù)表。如本例中轉(zhuǎn)換完后的數(shù)據(jù)直方圖如圖3所示。

圖2 銅含量直方圖分布

圖3 銅含量正態(tài)數(shù)據(jù)直方圖

2.3 試驗變差函數(shù)分析（Experimental Variogram Calculation）

幾乎所有的空間變量都是具有各向異性的特征。而這種空間的變異性則是由變差函數(shù)來表征。這種空間變異通過所擬合的理論模型，最終將體現(xiàn)在對空間變量的估值或者是模擬中。所以變差函數(shù)的求取是非常重要的一步。經(jīng)典地質(zhì)統(tǒng)計學中認為，在本征假設條件下，理論變差函數(shù)γ（α，｜h｜）是沿空間任一方向α，相距為｜h｜的兩個區(qū)域化變量z（x）和z（x＋h）的增量的方差，見式（1）。

式中：h稱為滯后距（lag），是一個矢量距離，決定于方向α和之間的距離增量｜h｜。在實際的計算中，把由取樣點計算而來的變差函數(shù)稱為實驗變差函數(shù)（experimental variogram），記為（h），它是理論變差函數(shù)γ（h）的估計值。計算公式見式（2）。

式中：N（h）是具有相同滯后距的點對個數(shù)。所以試驗變差函數(shù)是求取了［z（x）－z（x＋h）］的算術平均值。對于不同的滯后距分別計算相應的（h），這些點繪在h－（h）圖上則為試驗變差函數(shù)圖。

使用Gslib計算并繪制試驗變差函數(shù)圖主要用到兩個程序Gamv和Vargplt兩個程序。利用Gamv是依據(jù)式（2）來計算變差函數(shù)，使用vargplt這個程序繪圖。

在Gamv程序的運行中幾個關鍵的參數(shù)包括：變差函數(shù)計算方向、方向容限、滯后距，滯后距容限等幾個參數(shù)，見圖4。由于試驗變差函數(shù)取決于方向α，所以必須首先選定方向。而由于在選定的方向上可能的點對比較少，所以根據(jù)地質(zhì)實際情況可以設定一個容限。對于滯后距的選擇一般選擇多數(shù)樣品的間距作為滯后距，同樣由于一個滯后距不可能提供足夠多的點對計算一個較穩(wěn)定的算術平均值，所以對于一個滯后距也可以設定一個容限。同時也要考慮滯后距的個數(shù)（nlag），一般來說nlag＊lag應該覆蓋是搜索方向樣本點分布區(qū)域的1.5倍。

圖4 試驗變差函數(shù)參數(shù)設定

利用本數(shù)據(jù)集，沿x，y方向計算的試驗變差函數(shù)如圖5所示。在用Vargplt進行繪圖后，可以根據(jù)計算結(jié)果調(diào)整Gamv程序中計算試驗變差函數(shù)所需要的滯后距、滯后距個數(shù)以及角度容限等參數(shù)，然后重新計算重新繪圖。經(jīng)過這樣幾個反復的調(diào)整，使得能沿著選定的方向穩(wěn)定試驗變差函數(shù)應該能達到數(shù)據(jù)的總基臺值（Sill）值。但是注意在帶狀各向異性（zonal anisotropy）的影響可能無法達到基臺值。本例中，由于之前進行了正態(tài)轉(zhuǎn)換，所以數(shù)據(jù)的基臺值應該就是1.0。沿此方向的變程（r）則是試驗變差函數(shù)達到基臺值時的滯后距。

圖5 X、Y方向試驗變差函數(shù)圖

在繪出試驗變差函數(shù)后，應該首先檢查的是否與我們對于此區(qū)域的地質(zhì)認識符合。而另一方面試驗變差函數(shù)又可以佐證我們的地質(zhì)分析，二者應該是相輔相成的。如此例中，可以得出銅含量數(shù)據(jù)在X方向和Y方向上沒有很大的差異，基本上是具有相同的空間變異性。

2.4 理論變差函數(shù)擬合

試驗變差函數(shù)僅在所選定的方向上，在不同的滯后距位置進行了計算，而實際空間估值時是需要一個理論上滿足正定性要求，并且在空間任何方向上，任何距離都可以求取值的函數(shù)［6］。這就需要用理論的變差函數(shù)模型來擬合試驗變差函數(shù)。經(jīng)典的理論變差函數(shù)模型主要有，球狀模型，指數(shù)模型，高斯模型以及純塊金效應模型等［6］。而對于試驗變差函數(shù)的擬合需要使用Gslib中的Vmodel和Vargplt程序。

具體的擬合過程是：首先通過觀察試驗變差函數(shù)確定基臺值，變程以及理論變差函數(shù)類型，比如此例中，可以設定沿X，Y方向的塊金效應值設為0.3；理論模型可以選指數(shù)模型（exponential model），基臺值取0.7。將以上參數(shù)輸入Vmodel程序的參數(shù)文件，可以計算理論變差函數(shù)值。然后用Vargplt將二者繪制在同一圖中進行觀察。如果認為需要調(diào)整，則可以重新調(diào)整Vmodel中的變程，和基臺值等，然后重新計算，繪制進行觀察。這樣進行調(diào)整直至合適。本例中擬合的結(jié)果如圖6所示。本例中，擬合的結(jié)果用數(shù)學公式表示。

圖6 理論變差函數(shù)的人工擬合結(jié)果

在人工擬合，目視檢查的過程中可以體現(xiàn)我們的主動性，而且在擬合過程中的不斷的調(diào)整理論模型的參數(shù)等可以加深我們對于研究區(qū)的認識，增加我們分析變差函數(shù)的經(jīng)驗。

2.5 理論變差函數(shù)自動擬合（varfit）

人工擬合雖然可以將我們的認識體現(xiàn)在最終的擬合結(jié)果中，但同時這也是一種主觀性，過則不及。為了避免過度的主觀性，而且為了快速的對眾多的變差函數(shù)進行擬合，不少學者發(fā)表了自動擬合的程序，此處選用了Eulogio Pardo－Iguzquiza在1999年發(fā)表的Varfit程序［10］。程序中使用了加權(quán)最小二乘法擬合各個試驗變差函數(shù)點。本例中使用自動擬合的結(jié)果如圖7所示。

圖7 使用程序varfit人工擬合結(jié)果

自動擬合的結(jié)合，用數(shù)學公式表達。

用自動擬合軟件可以減輕我們的工作量，但建模者的認識和思想是任何軟件都無法取代的。做好這項工作要求把數(shù)學與地質(zhì)相結(jié)合。研究者既要比較了解熟悉研究區(qū)的的地質(zhì)背景，同時對地質(zhì)統(tǒng)計學有一定的了解。

3 結(jié)論

利用開源程序完全可以進行地質(zhì)統(tǒng)計中的各項工作，將這套過程和方法總結(jié)起來可以很好地指導礦產(chǎn)的勘探。正確運用這一過程和方法可以最大限度的減少變差函數(shù)分析過程中的人為因素造成的失誤和錯誤，同時開源程序分析環(huán)境的搭建使得從事變差函數(shù)的過程更透明，更可控。更好的指導礦產(chǎn)的勘探工作。文章僅以地質(zhì)統(tǒng)計學分析中最基本的步驟變差函數(shù)分析為例進行了說明。實際工作中的后續(xù)工作則是盡可能的利用原始數(shù)據(jù)對擬合的理論變差函數(shù)進行交叉驗證。如果滿足誤差要求則可以將擬合的變差函數(shù)應用到空間插值或者是隨機模擬。

［1］ Matheron G.Principle of geostatistics［J］.Economic Geology，1963，58（8）：1246－1266.

［2］ Wackernagel H.Multivariate Geostatistcs［M］.Berlin Heidelberg：Springer－Verlag，2003.

［3］ Atkinson P.M.，Lioyd C.D.Geostatistics for Environmental Applications［M］.London：Springer，2010.

［4］ Webster R.，Oliver M.A.Geostatistics for environmental scientists［M］.John Wiley＆Sons，2007.

［5］ Oliver M.A.Geostatistical Applications for Precision Agriculture［M］.LondonSpringer，2010.

［6］ Deutsch C.V.，Journel A.G，GSLIB：Geostatistical Software Library and User′s Guide［D］.Oxford University Press，1997.

［7］ http：//www.leg.ufpr.br/geoR/.

［8］ http：//www.gstat.org/.

［9］ Clark I.Practical Geostatistics［M］.London，Applied Science Publishers，1979.

［10］ Pardo－Igúzquiza，Eulogio VARFIT：a fortran－77program for fitting variogram models by weighted least squares［J］.Computers ＆Geosciences，1999，25（3）：251－261.