胡　飞, 陈建国, 曾小华, 黎　蓉, 李　靖, 李　鹏

(1.中国地质大学 资源学院,湖北 武汉　430074; 2.湖北省地质调查院,湖北 武汉　430034; 3.湖北省地质灾害防治中心,湖北 武汉　430034;4.湖北省地质环境总站,湖北 武汉　430034)

0　引言

分形已被广泛应用到自然科学的众多领域,在地球化学领域中,主要应用分形技术进行化探异常下限的确定,成秋明[1-3]指出化探数据具有多重分形结构特征,其背景值往往服从正态分布或对数正态分布,而高低异常值服从多重分形分布。与统计方法相比,分形方法具有将样品空间分布和统计特征随空间尺度变化性纳入化探数据的分形过程中,本质是采用含量与频数的关系来表达,即统计大于某一含量的数据个数,主要是含量-面积法(C-A模型)和含量-频率法(C-S模型)。MAPGIS 6.7系统是中国地质调查局指定使用的地质工作软件平台,广泛应用于各地勘单位,利用MAPGIS平台实现化探数据分形处理意义重大。

1　分形C-A模型

分形C-A模型是基于分形理论,研究含量与曲面面积之间关系的方法,并将此方法用于化探异常的研究。分形是指局部与总体具有某种相似性的形状,或者是在不同尺度上看起来基本相似的形状,即具有自相似性(self-similarity),而自相似性就是局部与整体在形态、功能和信息等方面具有统计意义上的相似[4]。定量描述这种自相似性的参数称为“分维数”,用D表示。 分形分布的特点要求大于等于某一尺度的数目与物体大小之间存在幂指数关系,由此可得到含量-面积分形统计模型,即：

A(r)=Cr-D,r>0

①

式①中：r表示元素含量值；A(r)表示元素含量≥r的面积值;C>0表示比例常数;D>0表示分维数。

分形C-A模型的关键在于确定分维数,具体方法是将公式①两端同时取对数,得出线性回归方程②。

lgA(r)=-Dlgr+lgC

②

把数据[A(r1),A(r2),…,A(rn)]和(r1,r2,…,rn)分别代入公式②,采用最小二乘法进行直线拟合,计算斜率D,一般情况散点会分别在两端或多段直线上,采用分段拟合方式,为提高精度避免人为确定分段拐点所造成的误差较大,可采用最优法确定分界点,使2个或多个拟合区间的直线与原始数据之间的残差平方和Ei(i=1,2)在各区间的总和以公式③表示为最小。

③

式中：rj为分界点；D1和D2分别为相应区间斜率,即分维数,多段直线对应D1,D2,D3…即多重分维,用维数确定异常下限后,把大于异常下限的数据删除,利用所有异常点数据绘制异常图。

2　MAPGIS 6.7平台实现

MAPGIS地理信息系统是武汉中地信息工程有限公司研制的,具有自主版权的大型基础地理信息系统软件平台,是集数字制图、数据库管理以及空间分析为一体的空间地理信息系统[5-6]。MAPGIS是目前国内地勘单位广泛使用的GIS平台,在此平台上实现化探数据多维分析处理显得尤为重要。

在MAPGIS平台上,用分析C-A模型确定化探异常下限的基本方法步骤,先将野外采集的化探数据进行网格化,并绘制成平面等值线图;然后利用MAPGIS属性分析功能,统计出不同尺度的含量r值所对应的面积A(r)值;最后利用GRAPHER软件,将数据[A(r1),A(r2),…,A(rn)]和(r1,r2,…,rn)绘制在双对数坐标图上,采用最小二乘法拟合成两条或多条直线,从而求解出分维数D的估计值和对应的拐点,拐点处所对应的元素含量值即为该元素的异常下限[7]。

具体实现步骤如下:

(1) 化探数据按采样坐标X,Y,元素含量值保存为MAPGIS能导入的文本,如W.txt 。

(2) 利用MAPGIS空间分析功能DTM子模块,对W.txt文件进行网格化,选择合适的网格化方法(克里格法、最小曲率、邻近点法等),输出W.grd网格文件;根据采样数量设定合适的(r1,r2,…,rn)值,打开“平面等值线图绘制”菜单,选择W.grd文件,按r1,r2,…,rn设置参数,生成区文件W.wp 。

(3) 利用MAPGIS空间分析子模块,打开W.wp区文件,使用“属性分析”的“双属性分类统计”功能将(r1,r2,…,rn)各个值段对应的面积统计并导出为W.wb文件。

(4) 利用MAPGIS库管理模块转换W.wb文件为W.dbf,在Excel软件中导入的数据文件W.dbf,累加(r1,r2,…,rn)各个值段对应的面积,求出A(r1),A(r2),…,A(rn),并利用Excel直接计算出lgri和lgA(ri)。

(5) 利用Excel或Grapher绘制lgri、lgA(ri)散点图,利用最小二乘法公式③,求出分维数D1和D2及分界点含量值j,j即为元素异常下限值。

3　应用实例

3.1　研究区地质概况

3.1.1区域地质特征

矿区位于桐柏—大别变质核杂岩隆起亚带中部。区域出露主要为TTG系列的英云闪长质片麻岩、奥长花岗质片麻岩、花岗闪长质片麻岩组合;表壳岩系为沉积其上的含铁碎屑建造、孔兹岩系及碎屑—碳酸盐岩建造夹双峰式火山建造(称大别山(岩)群);而侵入岩则主要表现为晋宁期和燕山期的花岗岩类。

造山带内岩浆活动十分频繁,发育多期岩浆岩,尤以燕山期后碰撞钙碱性花岗岩最为强烈,工区北部有灵山、新县和商城等岩体分布,规模较小的岩体主要形成于早白垩世,其中中酸性小岩体密布,岩性多为酸性富碱的花岗斑岩、似斑状花岗岩、石英斑岩及花岗闪长斑岩等。

3.1.2矿床地质特征

矿区出露地层属秦岭—大别地层区卡房—龟峰山小区,出露主要岩性有黑云二长花岗质片麻岩、斜长角闪岩、斜长角闪片岩等。广泛出露的二长花岗质片麻岩以及基性火山岩组,岩石化学性质十分稳定,岩石中节理、裂隙发育,为良好的成矿围岩,也提供了储矿空间。区内岩浆岩活动频繁,燕山期岩浆侵入尤为强烈,提供了丰富的岩浆热液和成矿物质来源。区内构造复杂,断裂构造极其发育,它们控制热液活动、地表水和地下水运移,控制成矿流体的运移,提供有利的导矿通道及成矿空间,进而控制矿体的空间分布。区内构造以断裂活动为主。按照空间展布特征可划分为北西西向、近东西向、北北东向三组,以北西西向较为发育,它们组成了区内网格状的断裂构造格架,基本与区域构造配套。其中两路口断层、狮子垴断层为工区主要断层。两路口断层区内出露长度约3 km,走向10°～30°,断层产状310°～340°∠45°～78°,形成的破碎带一般宽5～15 m,带内发育强硅化岩、构造角砾岩、碎裂岩及碎粒岩,明显与围岩面理斜交,常见断层引起的拖曳褶曲,构造岩具明显分带现象,从内到外依次为碎粉岩、角砾岩、碎裂岩,部分有早期糜棱岩残留,断裂带及其附近节理发育,沿断层常常见石英细脉沿裂隙面充填及黄铁矿化,该断层是区内规模最大的断裂构造,是商城—麻城断裂的次级断裂,且两路口钨矿床主矿体分布于其构造破碎带内。

3.2　确定异常下限

以矿区土壤化探数据为例,有样品2 600个,分析Ag、As、Bi、Mo、Pb、Sn、W、Zn元素,按C-A模型利用MAPGIS软件进行计算求出r和A(r),取对数部分统计结果如表1。

根据统计数据制作散点图,如图1。

以Ag元素为例,利用公式②③求拟合两段直线方程为:

lgA(r)= 0.272lgr+6.859 2

(-3.91

lgA(r)=-2.627lgr-0.090 2

(-2.95

因此分维数为D1=0.272,D2=2.627,异常下限r0=0.052 3,分界点r0的地质意思是含量<0.052 3为背景数据,>0.052 3为异常数据。经分析表明:

(1) 在双对数(lgr-lgA)图中,元素含量与面积的投影点呈现曲线分布,可用两条或多条直线去拟合,显示出一种连续多重分型特征。其中第一段直线非常平坦,D1数值接近零,大致反应了元素含量的低值背景,第二段直线斜率较陡,反应元素含量的高值区及异常区,两条直线的拐点对应元素含量可视为异常下限。

(2) 图中分维值D定量反应了元素含量在空间上的丛集集度和不均匀程度,D值越大,则元素空间分布均匀,且主要集中在低含量区,在地质上成矿性能极小,如As元素D2=8.48;D值越小,反应元素低值含量点到高值含量点的频率下降得越慢,元素在空间上

分布越不均匀,存在着较多高含量点,在地质上有富集成矿的可能,如W元素D2=1.04、Mo元素D2=1.43,在两路口地区主攻矿种就是钨钼矿,分形分维值真实反映了地质事实,经长期分形数据处理笔者认为D2值<2,尤其是<1.8有极大可能性富集成矿[8]。

(3) 图中有些元素多重分维,分维数D3明显小于D2,在图像上曲线明显偏离第二条直线段,如Mo、Pb、Sn、Zn元素,说明这些元素在高含量点基础上有二次局部富集[9],在地质上极有可能是成矿元素的伴生元素,且本区Ag与Pb、Zn富集没有相关性,Mo与Pb、Zn富集有一定相关性,W、Bi富集具有较强相关性,Mo与W、Bi富集也存在相关性。

表1　lgr(含量)与lgA(面积)统计结果Table 1　Statistical results of lgr (content) and lgA (area)

图1　两路口矿区土壤化探元素C-A分形双对数图Fig.1　Chart of double logarithm of C-A fractal dimension of soil elements from geochemical survey in Liangluokou Mining Area

3.3　绘制异常图

依据分维确定的拐点处对应的元素值,可直接绘制元素异常图,以W元素为例,依据多重分维D1、D2、D3,确定区域异常下限为7.39 μg/g,局部异常(矿致异常)下限为40 μg/g绘制异常图,如图2。

图2　W元素分形异常图Fig.2　Fractal anomaly of W element

图3　W元素异常图Fig.3　Anomaly map of W element

对比传统方法绘制的异常图,如图3。分形C-A模型处理方法保留了所有异常数据,其区域异常能较好反应矿区构造特征,体现W元素富集与东西向构造的密切关系;其局部异常面积小于传统方法,更为有效地体现了最小面积、最大含矿率原则。

4　结论

基于分形C-A模型,确定化探异常下限主要有以下特点：

(1) 反映化探数据元素特征。传统方法确定地球化学元素异常下限值的前提,是元素含量呈正态分布或者对数正态分布,所以传统方法通常需要多次剔除特高值,不断检验是否符合正态分布,即“检验—计算—再检验—再计算”的过程,分形C-A模型法确定异常下限值时,无需数据服从正态分布,也不用对数据进行特高值剔除处理,大大降低了人为因素的干扰,能较客观地反映研究区成矿元素的分布形态。

(2) 异常下限确定有理论依据。对服从多维分形的化探数据运用最小二乘法拟合曲线,斜率不同的直线段反映的是不同期次的地球化学场,所对应的分界值往往可作为区分背景和异常的临界值[10-11],这为研究地球化学元素的分布规律及其异常下限确定提供了理论支持,且元素分维值D2越小表示该元素越有富集成矿的可能,可利用分维值快速确定主攻矿种。

(3)C-A模型分析实现容易。仅利用MAPGIS和Excel就能实现化探数据分形C-A模型分析,操作简单易行。结合传统方法确定异常下限进行对比分析,能取得较好的效果,提高矿产预测的可靠性和客观性,对指导矿区下一步工作具有参考意义。

参考文献：

[1]Cheng Q,Agterberg F P,Bonham-Carter G F.A spatial analysis method for geochemical anomaly separation[J].J Geochem Explor,1996，56:185-195.

[2]Cheng Q.Multifractality and spatial statistics[J].Computer & Geosciences,1999,5(9):949-962.

[3]Cheng Q,Agterberg F P,Bonhaml,Ballantyne S B.The separation of geochemical anomalies from background by fractal methods[J].J Geochem Explor,1994，51:109-130.

[4]成秋明.空间模式的广义自相似分析与矿产资源评价[J].地球科学,2004,29(6):733-743.

[5]张钊,韦龙明,陈三明.MAPGIS在地质填图及化探数据处理中的应用[J].物探化探计算技术,2010,32(2):119-122.

[6]吴信才.MAPGIS地理信息系统[M].北京:电子工业出版社,2007.

[7]杨震,梅红波.基于MAPGIS的含量—面积法确定异常下限：以个旧西区为例[J].物探化探计算技术,2013,35(4):477-480.

[8]李堃,胡光道,刘才泽.个旧地区地球化学元素的多重分型特征与找矿预测[J].矿产地质,2006，20(4/5):498-502.

[9]张庆敏,岳云娟.基于小波模极大值求取多重分形谱多重分形克里格算法探究[J].鲁东大学学报：自然科学版,2013(1):22-25.

[10]陈建国,王仁铎,陈永清.利用分形统计学提取化探数据中的隐蔽信息并圈定地球化学异常[J].地球科学,1998,25(3):311-318.

[11]成秋明,赵鹏大,陈建国,等.奇异性理论在个旧锡铜矿产资源预测中的应用:成矿弱信息提取和复合信息分解[J].地球科学,2009,34(2):232-242.