白 林,唐 攀,李 楠

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

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

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

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

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

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

0 引言

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

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

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

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

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

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

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

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

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

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

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

1.2 圖形顯示

針對(duì)導(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í)驗(yàn)變異函數(shù)計(jì)算

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

實(shí)驗(yàn)變異函數(shù)的計(jì)算公式[8]：

(1)

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

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

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

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

圖3 實(shí)驗(yàn)變異函數(shù)計(jì)算界面Fig.3 Computing of experimental variogram

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

圖4 不同方位角、傾角、滯后距組合計(jì)算界面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)計(jì)學(xué)的計(jì)算，都必須知道變異函數(shù)的理論模型γ(h)以及模型中的參數(shù)(如變程a，基臺(tái)值C+C0和塊金常數(shù)C0等)。在實(shí)際地質(zhì)工作中，大部分的變異函數(shù)都是采用球狀模型，所以系統(tǒng)開發(fā)選擇在球狀模型的基礎(chǔ)上來討論理論變異函數(shù)的擬合問題[9]。

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

(2)

1.4.1 擬合方法選擇

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

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

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

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

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

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

(3)

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

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

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

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

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

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

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

表1 矽卡巖礦體Micromine計(jì)算的模型參數(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ù)計(jì)算系統(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ù)計(jì)算系統(tǒng)，實(shí)現(xiàn)了實(shí)驗(yàn)變異函數(shù)計(jì)算和理論變異函數(shù)擬合等功能，并應(yīng)用于甲瑪?shù)V山矽卡巖礦體Cu元素變異函數(shù)模型計(jì)算，計(jì)算結(jié)果經(jīng)交叉驗(yàn)證，具有較好的均值和標(biāo)準(zhǔn)方差。

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

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

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

[2]侯景儒,尹鎮(zhèn)南,李維明,等.實(shí)用地質(zhì)統(tǒng)計(jì)學(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].物探化探計(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)計(jì)學(xué)在甲瑪矽卡巖型銅多金屬礦儲(chǔ)量計(jì)算中的應(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報(bào),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)計(jì)學(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)計(jì)學(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)查局項(xiàng)目(121201103000150003)

白林(1980-)，男，博士，研究方向?yàn)榈厍蛱綔y(cè)與信息技術(shù)，E-mail：bailin@mail.cdut.edu.cn。

1001-1749(2016)06-0843-05

P 628+.2