甲瑪銅多金屬礦三維變異函數計算系統研發及應用

白 林,唐 攀,李 楠

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

甲瑪銅多金屬礦三維變異函數計算系統研發及應用

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

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

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

變異函數； 地質統計學； Matlab

0 引言

變異函數是地質統計學的重要工具。變異函數能夠表征地質變量的空間變化特征，不同方向上的實驗變異函數可反映區域化變量的各向異性；理論變異函數的參數各有含義，變程反映了區域化變量的影響范圍，基臺值反映了區域化變量變化幅度的大小；同時，變異函數也是應用克立格法進行礦山儲量估算的基礎[1-2]。

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

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

1 變異函數計算系統設計

變異函數計算系統針對西藏甲瑪銅多金屬礦研發，系統包括：①數據導入；②圖形顯示；③實驗變異函數計算；④理論變異函數擬合等模塊[5]。系統結構設計各模塊功能如下(圖1)：

1)數據導入模塊。將鉆孔、探槽等組合樣數據導入系統。

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

3)實驗變異函數計算模塊。針對導入系統的數據，輸入滯后距、滯后距個數、方位角、傾角等，計算三維實驗變異函數。

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

1.1 數據導入

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

1.2 圖形顯示

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

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

1.3 實驗變異函數計算

實驗變異函數計算模塊的功能是輸入方位角、傾角、滯后距、滯后距個數等，對導入系統的數據計算三維實驗變異函數，并顯示實驗變異函數圖。

實驗變異函數的計算公式[8]：

(1)

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

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

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

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

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

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

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

1.4 理論變異函數模型擬合

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

球狀模型理論變異函數形式為：

(2)

1.4.1 擬合方法選擇

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

圖5 擬合參數設定界面設計Fig.5 Setting of parameter fitting interface

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

圖6 權系數為1的擬合效果Fig.6 Fitting effect when the weight is 1

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

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

(3)

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

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

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

1.4.2 形態較差的變異函數擬合

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

圖9 形態較差的變異函數擬合方法Fig.9 Fitting of sick experimental variogram

1.4.3 擬合效果評價

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

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

表1 矽卡巖礦體Micromine計算的模型參數

Tab.1 Parameters of variogram for skarn deposit by Micromine

方向方位角傾角滯后距塊金值第一結構第二結構a1C1a2C2走向3354700.621540.114130.6傾向6614700.621850.112820.6垂向22975100.62370.11900.6

表2 矽卡巖礦體Matlab變異函數計算系統結果

Tab.2 Parameters of variogram for skarn deposit by Matlab

方向方位角傾角滯后距塊金值第一結構第二結構a1C1a2C2走向3354700.621980.254280.4傾向6614700.621850.252820.4垂向22975100.62370.25900.4

2 結論

1)作者基于Matlab語言編寫了三維變異函數計算系統，實現了實驗變異函數計算和理論變異函數擬合等功能，并應用于甲瑪礦山矽卡巖礦體Cu元素變異函數模型計算，計算結果經交叉驗證，具有較好的均值和標準方差。

2)開發的系統設計了按方位角、傾角、滯后距不同組合方式計算實驗變異函數，并分別繪制變異函數曲線，根據變異函數曲線形態選擇接近標準變異函數曲線的圖形進行擬合，回避了大多軟件繁瑣的參數調整過程。

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

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

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

[3]矯希國,劉超.變差函數的參數模擬[J].物探化探計算技術,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軟件下地質統計學在甲瑪矽卡巖型銅多金屬礦儲量計算中的應用[J].地質與勘探,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]唐菊興,王登紅,汪雄武,等.西藏甲瑪銅多金屬礦礦床地質特征及其礦床模型[J].地球學報,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].礦床地質,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]王仁鐸,胡光道.線性地質統計學[M].北京：地質出版社，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]孫洪泉.地質統計學及其應用[M].北京：中國礦業大學出版社，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

中國地質調查局項目(121201103000150003)

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

1001-1749(2016)06-0843-05

P 628+.2