應(yīng)用開源程序進(jìn)行礦業(yè)地質(zhì)統(tǒng)計(jì)學(xué)中的三維變差函數(shù)分析

方曉剛，李宇鵬

（中國石油大學(xué)（北京），北京102249）

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

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

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

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

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

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

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

2 變差函數(shù)分析

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

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

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

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

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

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

圖2 銅含量直方圖分布

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

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

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

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

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

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

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

圖4 試驗(yàn)變差函數(shù)參數(shù)設(shè)定

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

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

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

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

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

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

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

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

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

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

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

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

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

3 結(jié)論

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

［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.