甲瑪銅多金屬礦三維變異函數(shù)計算系統(tǒng)研發(fā)及應(yīng)用

白 林,唐 攀,李 楠

(1.成都理工大學(xué) a.數(shù)學(xué)地質(zhì)四川省重點實驗室,b.地球科學(xué)學(xué)院，成都 610059；2.中國地質(zhì)科學(xué)院 礦產(chǎn)資源研究所，北京 100037)

甲瑪銅多金屬礦三維變異函數(shù)計算系統(tǒng)研發(fā)及應(yīng)用

白 林1a,唐 攀1b,李 楠2

(1.成都理工大學(xué) a.數(shù)學(xué)地質(zhì)四川省重點實驗室,b.地球科學(xué)學(xué)院，成都 610059；2.中國地質(zhì)科學(xué)院 礦產(chǎn)資源研究所，北京 100037)

采用Matlab語言編寫了三維變異函數(shù)計算系統(tǒng)，實現(xiàn)了西藏甲瑪銅多金屬礦實驗變異函數(shù)計算和理論變異函數(shù)擬合等功能。設(shè)計按方位角、傾角、滯后距的不同組合計算實驗變異函數(shù)，以供用戶選擇形態(tài)較好的實驗變異函數(shù)圖；實現(xiàn)了直接法、加權(quán)回歸法和二級套合法3種不同的理論變異函數(shù)擬合方法；針對形態(tài)較差的實驗變異函數(shù)圖，提出模擬具有擬合經(jīng)驗的地質(zhì)工作者的思維方式，設(shè)定理論變異函數(shù)曲線形態(tài)，反推塊金、變程、基臺值等參數(shù)。系統(tǒng)應(yīng)用于西藏甲瑪銅多金屬礦矽卡巖礦體Cu元素變異函數(shù)模型計算，計算結(jié)果經(jīng)交叉驗證，具有較好的均值和標(biāo)準(zhǔn)方差。

變異函數(shù)； 地質(zhì)統(tǒng)計學(xué)； Matlab

0 引言

變異函數(shù)是地質(zhì)統(tǒng)計學(xué)的重要工具。變異函數(shù)能夠表征地質(zhì)變量的空間變化特征，不同方向上的實驗變異函數(shù)可反映區(qū)域化變量的各向異性；理論變異函數(shù)的參數(shù)各有含義，變程反映了區(qū)域化變量的影響范圍，基臺值反映了區(qū)域化變量變化幅度的大?。煌瑫r，變異函數(shù)也是應(yīng)用克立格法進(jìn)行礦山儲量估算的基礎(chǔ)[1-2]。

現(xiàn)有變異函數(shù)的計算和參數(shù)擬合常借助于商業(yè)軟件如(Micromine、Surfer等)，但它們的實驗變異函數(shù)計算步驟繁瑣、人為干預(yù)多，需要人工分別確定方位角、傾角、滯后距等以選擇較好的實驗變異函數(shù)圖[3]；同時理論變異函數(shù)擬合往往需要多次調(diào)整參數(shù)才能得到較好的擬合效果，要求操作者具有豐富的地質(zhì)經(jīng)驗[4]，對于廣大地質(zhì)工作者來說一直都是具有難度的“技術(shù)活”。

對于上述實驗變異函數(shù)計算和理論變異函數(shù)擬合中的問題，作者基于Matlab語言開發(fā)了變異函數(shù)計算系統(tǒng)，以組合計算的方式方便三維實驗變異函數(shù)計算，模擬專家經(jīng)驗實現(xiàn)對實驗變異函數(shù)模型的自動擬合，并應(yīng)用于西藏甲瑪銅多金屬礦變異函數(shù)計算。

1 變異函數(shù)計算系統(tǒng)設(shè)計

變異函數(shù)計算系統(tǒng)針對西藏甲瑪銅多金屬礦研發(fā)，系統(tǒng)包括：①數(shù)據(jù)導(dǎo)入；②圖形顯示；③實驗變異函數(shù)計算；④理論變異函數(shù)擬合等模塊[5]。系統(tǒng)結(jié)構(gòu)設(shè)計各模塊功能如下(圖1)：

1)數(shù)據(jù)導(dǎo)入模塊。將鉆孔、探槽等組合樣數(shù)據(jù)導(dǎo)入系統(tǒng)。

2)圖形顯示模塊。對組合樣數(shù)據(jù)進(jìn)行三維顯示，采用不同顏色來區(qū)分樣品的不同品位值。

3)實驗變異函數(shù)計算模塊。針對導(dǎo)入系統(tǒng)的數(shù)據(jù)，輸入滯后距、滯后距個數(shù)、方位角、傾角等，計算三維實驗變異函數(shù)。

4)理論變異函數(shù)擬合模塊。對實驗變異函數(shù)進(jìn)行擬合，計算理論變異函數(shù)變程、塊金、基臺值等參數(shù)。

1.1 數(shù)據(jù)導(dǎo)入

導(dǎo)入數(shù)據(jù)為各礦體數(shù)據(jù)，格式為：東經(jīng)、北緯、高程、元素品位值。甲瑪?shù)V山參與變異函數(shù)計算的數(shù)據(jù)為元素化學(xué)分析數(shù)據(jù)，包括:①矽卡巖陡礦體；②矽卡巖礦體；③角巖礦體三個礦體的組合樣數(shù)據(jù)。截至2012年底，參與計算的數(shù)據(jù)共包括矽卡巖陡礦體鉆孔、探槽98個，共6 830個1 m組合樣；矽卡巖礦體188個鉆孔，9 907個1 m組合樣；角巖礦體67個鉆孔，13 764個2 m組合樣。

1.2 圖形顯示

針對導(dǎo)入的礦體組合樣數(shù)據(jù),繪制其三維空間分布形態(tài)。圖2為矽卡巖礦體組合樣空間分布形態(tài)，根據(jù)礦體的幾何形態(tài)可以確定其大致的走向、傾向和垂向[6-7]。

圖2 矽卡巖礦體組合樣空間分布Fig.2 Spatial distribution of skarn deposit compositing drillhole

1.3 實驗變異函數(shù)計算

實驗變異函數(shù)計算模塊的功能是輸入方位角、傾角、滯后距、滯后距個數(shù)等，對導(dǎo)入系統(tǒng)的數(shù)據(jù)計算三維實驗變異函數(shù)，并顯示實驗變異函數(shù)圖。

實驗變異函數(shù)的計算公式[8]：

(1)

其中：γ*(h)為實驗變異函數(shù);h為步長，即在指定方向上，距離為|h|的矢量；N(h)表示步長為h的樣品對數(shù)；Z(x)、Z(x+h)為元素品位數(shù)據(jù)。

輸入數(shù)據(jù)：①原始數(shù)據(jù)。組合樣數(shù)據(jù)，包括東坐標(biāo)、北坐標(biāo)、高程、元素分析數(shù)據(jù)；②計算參數(shù)。方位角、傾角、滯后距、滯后距個數(shù)。

輸出結(jié)果：實驗變異函數(shù)(圖3)。

以矽卡巖陡礦體為例,計算其井向的實驗變異函數(shù)，如圖3所示。

圖3 實驗變異函數(shù)計算界面Fig.3 Computing of experimental variogram

為了確定礦體各主方向的方位角、傾角及滯后距，在Micromine軟件中需要分步驟畫出眾多實驗變異函數(shù)圖，對比選擇較好形態(tài)的圖形，以此先后確定礦體主方向的方位角、傾角和滯后距，計算步驟多、修改復(fù)雜。本系統(tǒng)的研發(fā)，設(shè)計出較為方便的計算方式，可根據(jù)方位角、傾角、滯后距的取值范圍進(jìn)行組合計算。如圖4中，方位角有60°、90°、120°三種情況，傾角固定為-4°，滯后距選取為50 m和70 m，則可計算出3×1×2=6個實驗變異函數(shù)圖，用戶可以直接對這6幅實驗變異函數(shù)圖進(jìn)行對比選擇。

圖4 不同方位角、傾角、滯后距組合計算界面Fig.4 Experimental variograms computing by the combination of various azimuths,dips and lags

1.4 理論變異函數(shù)模型擬合

無論是要了解區(qū)域化變量的變異特征，還是要進(jìn)一步進(jìn)行地質(zhì)統(tǒng)計學(xué)的計算，都必須知道變異函數(shù)的理論模型γ(h)以及模型中的參數(shù)(如變程a，基臺值C+C0和塊金常數(shù)C0等)。在實際地質(zhì)工作中，大部分的變異函數(shù)都是采用球狀模型，所以系統(tǒng)開發(fā)選擇在球狀模型的基礎(chǔ)上來討論理論變異函數(shù)的擬合問題[9]。

球狀模型理論變異函數(shù)形式為：

(2)

1.4.1 擬合方法選擇

理論變異函數(shù)擬合模塊是對實驗變異函數(shù)擬合出其理論模型參數(shù)變程、塊金、基臺值等。擬合方法包括(圖5)：①直接法；②加權(quán)多項式回歸法[10]；③二級套合法。

圖5 擬合參數(shù)設(shè)定界面設(shè)計Fig.5 Setting of parameter fitting interface

直接法較為簡單，加權(quán)回歸法需要注意的是，因為變異函數(shù)擬合時，前面3至4個點的重要性要高于后面的點，系統(tǒng)提供了給前面4個點設(shè)定權(quán)系數(shù)的功能。如圖6、圖7，針對矽卡巖陡礦體Cu元素井向?qū)嶒炞儺惡瘮?shù)，設(shè)定權(quán)系數(shù)分別為1和樣品對數(shù)N(h)時，采用加權(quán)回歸法的擬合效果。從圖6、圖7對比看得出來，權(quán)系數(shù)為樣品對數(shù)時，擬合曲線稍微對前2個點較為兼顧，擬合效果稍好。

圖6 權(quán)系數(shù)為1的擬合效果Fig.6 Fitting effect when the weight is 1

圖7 權(quán)系數(shù)為點對數(shù)的擬合效果Fig.7 Fitting effect when the weight is N(h)

二級套合球狀模型理論變異函數(shù)的數(shù)學(xué)表達(dá)式為式(3)。

(3)

二級套合的球狀模型適合對前段較徒，后段較緩的實驗變異函數(shù)曲線進(jìn)行擬合。系統(tǒng)專門設(shè)計了設(shè)定分界點功能，把全部數(shù)據(jù)點分為前，后兩部分，分界點就選在曲線由陡變緩的轉(zhuǎn)折點處。

針對矽卡巖陡礦體Cu元素井向?qū)嶒炞儺惡瘮?shù)，進(jìn)行3種擬合方法的效果對比(圖8)。圖8上方的線(深色)表示直接法，下方的線(淺色)表示加權(quán)回歸法、中間的線為二級套合法的擬合效果，從圖形上看二級套合法兼顧了前面的點和后面的點，更逼近實驗變異函數(shù)，擬合效果更好。

圖8 直接法加權(quán)回歸法二級套合法擬合效果對比Fig.8 Comparison of fitting effect of direct fitting,weighted regressing and nested-overlap

1.4.2 形態(tài)較差的變異函數(shù)擬合

計算出的實驗變異函數(shù)并非都具有標(biāo)準(zhǔn)變異函數(shù)的形態(tài)。以西藏甲瑪銅多金屬礦矽卡巖礦體Cu元素為例，計算出的實驗變異函數(shù)分布散亂、形態(tài)較差(指與標(biāo)準(zhǔn)球狀模型形態(tài)差別較大)。對于此類形態(tài)較差的實驗變異函數(shù)，往往需要地質(zhì)工作者擁有豐富的變異函數(shù)擬合經(jīng)驗，通過不斷的調(diào)整和后期的交叉驗證得到最終的擬合結(jié)果。但并非所有人都具有變異函數(shù)擬合的經(jīng)驗，作者提出了一種模擬具有變異函數(shù)擬合經(jīng)驗的地質(zhì)工作者的思維方式，由操作者根據(jù)實驗變異函數(shù)由經(jīng)驗設(shè)定理論變異函數(shù)曲線形態(tài)，由此反推理論變異函數(shù)參數(shù)。圖9為甲瑪?shù)V山矽卡巖礦體Cu元素的變異函數(shù)輔助擬合效果，黑色實心點為實驗變異函數(shù)，可以看出原始數(shù)據(jù)較為分散，直接擬合不出理論變異函數(shù)，可以通過圖中添加空心圓圈的方式，選擇參與理論變異函數(shù)擬合的數(shù)據(jù)，設(shè)定理論變異函數(shù)形態(tài)，由此擬合出模型參數(shù)。

圖9 形態(tài)較差的變異函數(shù)擬合方法Fig.9 Fitting of sick experimental variogram

1.4.3 擬合效果評價

對于系統(tǒng)計算的塊金、變程、基臺值等理論變異函數(shù)模型參數(shù)，可以采用交叉驗證的方法檢驗變異函數(shù)參數(shù)的可靠性。下面列出了矽卡巖礦體在Micromine軟件中人工計算的結(jié)果和Matlab編寫的變異函數(shù)系統(tǒng)計算結(jié)果的對比(表1和表2)，表1、表2中a1、a2表示第一結(jié)構(gòu)和第二結(jié)構(gòu)的變程，C1、C2表示第一結(jié)構(gòu)和第二結(jié)構(gòu)的部分基臺值。

交叉驗證的誤差統(tǒng)計的均值越接近“0”，標(biāo)準(zhǔn)方差越接近“1”，表明擬合的變異函數(shù)參數(shù)效果越好[11]。Matlab變異函數(shù)計算系統(tǒng)計算的誤差統(tǒng)計均值、標(biāo)準(zhǔn)方差為(0.001 293,1.070 1)，相對于Micromine的(0.001 393，1.080 3)稍小，表明研發(fā)的系統(tǒng)計算結(jié)果達(dá)到了Micromine軟件的計算效果。

表1 矽卡巖礦體Micromine計算的模型參數(shù)

Tab.1 Parameters of variogram for skarn deposit by Micromine

方向方位角傾角滯后距塊金值第一結(jié)構(gòu)第二結(jié)構(gòu)a1C1a2C2走向3354700.621540.114130.6傾向6614700.621850.112820.6垂向22975100.62370.11900.6

表2 矽卡巖礦體Matlab變異函數(shù)計算系統(tǒng)結(jié)果

Tab.2 Parameters of variogram for skarn deposit by Matlab

方向方位角傾角滯后距塊金值第一結(jié)構(gòu)第二結(jié)構(gòu)a1C1a2C2走向3354700.621980.254280.4傾向6614700.621850.252820.4垂向22975100.62370.25900.4

2 結(jié)論

1)作者基于Matlab語言編寫了三維變異函數(shù)計算系統(tǒng)，實現(xiàn)了實驗變異函數(shù)計算和理論變異函數(shù)擬合等功能，并應(yīng)用于甲瑪?shù)V山矽卡巖礦體Cu元素變異函數(shù)模型計算，計算結(jié)果經(jīng)交叉驗證，具有較好的均值和標(biāo)準(zhǔn)方差。

2)開發(fā)的系統(tǒng)設(shè)計了按方位角、傾角、滯后距不同組合方式計算實驗變異函數(shù)，并分別繪制變異函數(shù)曲線，根據(jù)變異函數(shù)曲線形態(tài)選擇接近標(biāo)準(zhǔn)變異函數(shù)曲線的圖形進(jìn)行擬合，回避了大多軟件繁瑣的參數(shù)調(diào)整過程。

3)對于形態(tài)較差的變異函數(shù)曲線，提出模擬具有擬合經(jīng)驗的地質(zhì)工作者的思維方式，設(shè)定期望的理論變異函數(shù)曲線，再反推模型參數(shù)，取得較好效果，具有一定的創(chuàng)新性。

[1]CHILES J P,DELFINER P.Geostatistics:modeling spatial uncertainty[M].America:Wiley-Interscience,1999.

[2]侯景儒,尹鎮(zhèn)南,李維明,等.實用地質(zhì)統(tǒng)計學(xué)[M].北京:地質(zhì)出版社,1998.HOU J R,YIN Z N,LI W M,et al.Practical geostatistics[M].Beijing:Geological Publishing House,1998.(In Chinese)

[3]矯希國,劉超.變差函數(shù)的參數(shù)模擬[J].物探化探計算技術(shù),1996,18(02):157-161.JIAO X G,LIU C.Estimation of variation Parameter[J].Computing Techniques for Geophysical and Geochemical Exploration,1996,18(02):157-161.(In Chinese)

[4]鄭文寶,黎楓佶,唐菊興,等.基于Micromine軟件下地質(zhì)統(tǒng)計學(xué)在甲瑪矽卡巖型銅多金屬礦儲量計算中的應(yīng)用[J].地質(zhì)與勘探,2011,47( 4):726-736.ZHENG W B,LI F J,TANG J X,et al.The application of geostatistics to ore reserve calculation of the Jiama skarn type copper-polymetallic deposit based on Micromine software [J].Geology and Exploration,2011,47( 4) :726-736．(In Chinese)

[5]DEUTSCH C V,JOURNEL A G.GSLIB:Geostatistical Software Library and User's Guide[M].Oxford:Oxford University Press,1992.

[6]唐菊興,王登紅,汪雄武,等.西藏甲瑪銅多金屬礦礦床地質(zhì)特征及其礦床模型[J].地球?qū)W報,2010,31(4):495-506.TANG J X,WANG D H,WANG X W,et al.Geological features and metallogenic model of the jiama copper-polymetallic deposit in tibet[J].Acta Geoscientica Sinica,2010,31(4):495-506.(In Chinese)

[7]唐菊興,鄧世林,鄭文寶,等.西藏墨竹工卡縣甲瑪銅多金屬礦床勘查模型[J].礦床地質(zhì),2011,30(02):179-196.TANG J X,DENG S L,ZHENG W B,et al.An exploration model for jiama copper-polymetallic deposit in maizhokunggar county Tibet[J].mineral deposits,2011,30(2):179-196.(In Chinese)

[8]王仁鐸,胡光道.線性地質(zhì)統(tǒng)計學(xué)[M].北京：地質(zhì)出版社，1988.WANG R D,HU G D.Linear geostatistics[M].Beijing:Geological Publishing House,1988.(In Chinese)

[9]CRESSIE N.Fitting variogram models by weighted least squares[J].Mathematical Geology,1985,17(5):563-586.

[10]PARDO-IGUZQUIZA E.VARFIT:a Fortran-77 program for fitting variogram models by weighted least squares[J].Computers and Geosciences,1999,25(3):251-261.

[11]孫洪泉.地質(zhì)統(tǒng)計學(xué)及其應(yīng)用[M].北京：中國礦業(yè)大學(xué)出版社，1990.SUN H Q.Geostatistics and its application[M].Beijing:China University of Mining and Technology Press,1990.(In Chinese)

Development and application of three dimensional variogram computing system for Jiama copper-polymetallic deposit

BAI Lin1a,TANG Pan1b,LI Nan2

(1.Chengdu University of Technology,a.Geomathernatics key Laboratory of Chengdu University of Sichuan Province,b.College of Earth Sciences,Chengdu 610059,China;2.Institute of Mineral Resources,Chinese Academy of Geological Sciences,Beijing 100037,China)

The three dimensional variogram computing system could achieve the experimental variogram computing and theoretical variogram fitting based on Matlab for Jiama copper-polymetallic deposit.System proposed computing experimental variograms by various azimuths,dips and lags for choosing the best figure of variograms.Fitting of theoretical variogram in this system included direct fitting,weighted regressing and nested-overlap.It emphasized simulating the thinking patterns of geologists who are experienced with fitting variogram to set the fitting curve for sick experimental variogram.Furthermore,it computed the theoretical variogram parameters such as nugget,range and sill.This system was applied to computing variograms of the Cu element of Jiama copper-polymetallic skarn deposit in Tibet.After cross-validation,it obtains better mean and standard deviation.

variogram; geostatistics; Matlab

2015-09-09 改回日期：2015-10-08

中國地質(zhì)調(diào)查局項目(121201103000150003)

白林(1980-)，男，博士，研究方向為地球探測與信息技術(shù)，E-mail：bailin@mail.cdut.edu.cn。

1001-1749(2016)06-0843-05

P 628+.2