甲玛铜多金属矿三维变异函数计算系统研发及应用

白 林,唐 攀,李 楠

(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