圣安陈，杨悦，李亚彬，邹宗霖，翁爱华

吉林大学 地球探测科学与技术学院，长春 130026

0 引言

花牛山铅锌矿床位于甘、新交界部位的北山地区，是该区一处中型的多金属矿床。该区矿产资源丰富，已发现有铅、锌、金、银和钨等多处金属矿床。矿田大地构造位置位于北山塔里木板块敦煌地块北缘活动带北缘，花牛山--黑山--双鹰山早古生代裂谷裂陷带中[1]。

矿区北邻方山口--营毛沱--牛圈子早古生代裂陷带，南接红十井--柳园--大奇山晚古生代陆内裂谷带[1]。矿区内褶皱和断裂发育。褶皱常伴随近EW向断裂发育，并控制着矿区内主要矿体的分布。矿区内岩浆活动频繁且强烈，分布较为广泛，从侵入岩到喷发岩，从超基性岩到酸性岩均有出露，岩浆活动明显受EW向与NW向构造带控制[2]。

甘肃省地质局物探队、西安地质矿产研究院及廊坊物化探所先后在花牛山铅锌矿区开展了勘查研究工作，发现矿体赋存有利部位是碳酸盐岩与碎屑岩接触带。相关文献表明，花牛山铅锌矿的有利成矿地段为深大断裂带边缘，成矿主要围岩为碳酸盐岩[3]。杨建国认为花牛山铅锌矿无论是以碎屑岩和碳酸盐岩为主的容矿区，还是以中基性火山夹碎屑岩和碳酸盐岩为主的容矿区，均赋存于震旦系上统三岩组第二岩性段千枚岩、千枚状板岩和第三岩性段大理岩夹千枚岩过渡部位，矿体主要产于大理岩与千枚岩接触带附近及其附近层间破碎带[4]。

可控源音频大地电磁法能克服天然场源信号的随机性和信号弱的缺点，已成为解决地球深部资源勘探问题的一种有效手段[5]。花牛山铅锌矿找矿模型过去一直从地面地质和地球化学方向出发，较少结合物探方法。为进一步研究铅锌矿区的三维电性结构，划分成矿有利区域，廊坊物化探所在花牛山地区开展了以可控源音频大地电磁法(CSAMT)为主的地球物理勘探工作。近年来，可控源方法的一维、二维反演技术已日趋成熟。但考虑到地下电性结构的复杂性，一维、二维反演可靠性较差，采用三维反演将是一个更有效的手段。目前三维反演方法主要有Occam反演法[6]、拟牛顿法[7]、非线性共轭梯度法[8--9]和快速松弛反演法[10]等。笔者采用的有限内存拟牛顿法[11]，是由Broyden[12--14]等在1970年提出并发展起来的，是目前三维电磁法反演中最有效的拟牛顿方法之一，由于该方法不需要存储大型海森矩阵，因此克服了牛顿法计算量大、产生非下降方向的缺点。2006年Avdeeva[15--16]首次将该方法应用到大地电磁(MT)一维反演中，并进一步推广到三维反演，非常适合求解大规模的最优化问题。笔者利用该反演方法获得花牛山铅锌矿测区的三维电性模型和界面形态，并结合前人的激电法测量结果，推测出可能的铅锌成矿区。

1 矿区地质及地球物理特征

花牛山铅锌矿田由4个矿区组成，出露地层主要为震旦系洗肠井群二、三、四岩组，矿体主要赋存于震旦系洗肠井群碎屑岩--碳酸盐岩岩系中，属火山沉积碎屑岩--碳酸盐岩系有关的喷流沉积型铅锌成矿系列[4]。成矿有利地层位于千枚岩与大理岩的接触带。矿床的形成受震旦系地层、岩浆热液和构造作用影响。二岩组(Zxb)为一套浅变质的浅海相泥质岩建造，在其内分布有大量酸性斑岩和石英脉，整体岩性较为单一；三岩组(Zxc)主要为一套浅海相碳酸盐岩与泥质岩建造，岩性主要为含粒状石英的大理岩夹绢云千枚岩及粉砂质板岩，铅锌矿体产于该岩组中；四岩组(Zxd)主要为板岩和角岩等[17]。矿田南部有EW向的花西滩--花牛山压(扭)性断裂(F1)，北部有EW向的五井河隐伏断裂(F3)，F1断裂规模较大，控制着岩体的侵位，是主要的控矿构造[18]。

测区位于一矿区，北部为震旦系洗肠井群三岩组大理岩夹千枚岩，南部为震旦系洗肠井群四岩组千枚岩和角岩。两者的接触带为二矿带赋存位置，位于三岩组第三岩性段大理岩夹千枚与三岩组第二岩性段千枚岩过度部位。东部、西部均有印支期花岗岩侵入。测区南端为三矿带，测点未覆盖该矿带。测区内有NW向的花黑滩--双峰山断裂(F2)与二矿带构造走向交叉，是后期次级断裂带，与矿体产状和走向有一定关系(图1)。

物性差异是地球物理工作的前提。测区内主要的岩石为千枚岩、含炭千枚岩、花岗闪长岩、大理岩和铅锌磁黄铁矿。5种主要岩石电阻率的差异较大，其中花岗闪长岩的电阻率最高，千枚岩、大理岩的电阻率相对较高，含炭千枚岩、铅锌磁黄铁矿的电阻率最低(表1)，这为可控源音频大地电磁反演异常的解释提供了基础。此外，测区内铅锌磁黄铁矿矿石、含炭千枚岩具有低阻高极化的特点。因此，含炭千枚岩是电法找矿的重要地质干扰因素。

表1 测区主要岩石电性特征[19]Table 1 Electrical characteristics of main rocks in survey area

2 数据及反演

可控源音频大地电磁法(CSAMT)是在大地电磁(MT)的基础上发展起来的一种人工源电磁测深方法，主要适用于油气、地热等资源勘探。目前，由于其勘探深度较大，逐渐成为深部找矿的重要地球物理手段。

1.中奥陶统花牛山群；2.震旦系洗肠井群四岩组；3.震旦系洗肠井群三岩组；4.震旦系洗肠井群二岩组；5.印支期正长花岗岩；6.印支期二长花岗岩；7.华力西期花岗闪长岩；8.橄榄岩；9.断裂；10.测区；11.CSAMT测线；F1.花西滩—花牛山断裂；F2.花黑滩—双峰山断裂；F3.五井河断裂。图1 花牛山矿区地质概况图Fig.1 Geological map of Huaniushan ore deposit

2.1 数据采集

本次研究采用的可控源音频大地电磁数据是廊坊物化探所独立研制的多功能电磁工作站采集的野外实测数据。测区面积约3 km2，观测方式采用赤道偶极标量方式，观测垂直测线方向的磁场分量Hy与沿测线方向的电场分量Ex，发射极距为1.5 km，接收极距为50 m，最大供电电流为 20 A，最小收发距为 8 km，共布设13条测线，分别命名为L1、L2到L13，测线L1距离发射源最远，测线L13距离发射源最近，线距为 200 m，点距为50 m，每条测线24个测点，一共312个测点。将收集到的数据进行处理，主要扫描时间序列，根据电磁干扰程度选择合适的频率段参与电场分量Ex估计。

数据质量对反演结果至关重要。图2为76号点、124号点、172号点和220号点本次反演的数据电场分量Ex。原始数据质量均较好，高频部分出现的跳点、飞点不明显，说明数据受到的噪声干扰较小，反演可以得到地下可靠的信息。

2.2 有限内存拟牛顿反演方法

根据Tikhonov正则化理论，地球物理反演问题的目标函数为：

Φ=Φd+λΦm

(1)

式中：Φd为观测数据和理论模型响应的拟合差；Φm为模型范数；λ为正则化子。

本文的正则化目标函数定义为：

(2)

式中：m是模型参数向量；d是数据向量；Cd是数据方差矩阵；Cm是模型方差函数；F是正演算子。

L--BFGS方法可以直接求取目标函数的最优化问题[20]。将模型的搜索方向定义为：

(3)

利用充分下降条件和曲率条件来约束搜索步长αk[21]。由于求解式(3)需要很长的耗时和很大的存储空间，L--BFGS在牛顿法的基础上，将Hessian矩阵的逆近似改进为：

(4)

第一步：令k=1，给定初始模型mk，初始正定矩阵Hk，允许误差范围ε>0；

第二步：计算梯度rk=-▽U(λ,mk)，当‖rk‖≤ε，输出最终解mk；

第三步：否则,利用(3)式进行步长αk的线性搜索。当步长αk小于阙值，减小正则化因子λ；如果λ小于阙值，输出最终解mk，反演终止；

第四步：更新模型mk+1=mk+αkPk，若‖rk+1‖≤ε，得到最终反演结果mk+1；

第五步：否则，利用(4)式更新正定矩阵Hk+1；

第六步：令k=k+1，返回第二步。

2.3 理论模型模拟

为研究有限内存拟牛顿法三维反演的可靠性，设计了理论模型(图3)。模型为均匀半空间，背景电阻率为100 Ω·m，交错放置8个300 m×300 m×300 m的立方异常体，红色代表高阻体为500 Ω·m，蓝色代表低阻体为10 Ω·m。发射源AB使用1 000 m的长导线源，最小发射极距为7 500 m，发射频率为2-2、2-1、…、213Hz，共计16个频点。测区(虚线框)布设11条测线，线距设置为100 m，点距设置为50 m，每条测线21个测点，共计231个测点，观测参数为电场Ex复分量。

图3 理论模型及工作装置示意图Fig.3 Theoretica model and working device diagram

理论模型数值模拟时，测区范围内将网格剖分为20×20×25，网格长度均为50 m，由于模型的范围应略大于测区范围，因此模型以2.0倍的系数各扩边5个网格，最终生成的网格为30×30×30，共计27 000个。反演部分参数的设置，背景电阻率为100 Ω·m，正则化因子初始值为1，叠加5%的噪声到理论数据。

图4三维反演结果深度为200 m和900 m的水平切片。图中可以看出浅部异常体，高、低阻均能清晰的刻画出异常体的区域，反演的电阻率大小与理论模型异常体电阻率比较接近。深部异常体的低阻的轮廓有所展现，高阻不明显，与理论模型电阻率异常体电阻率相差很大。说明有限内存拟牛顿法浅部反演结果是比较可靠的。

a.深度200 m切片；b.深度900 m切片。图4 Ex分量复数据三维反演结果Fig.4 3D model inversed results from Ex component

2.4 实测数据反演

虽然L--BFGS反演方法相比于非线性共轭梯度法(NLCG)在反演速度和反演结果准确性方面有一定的优势[11]，但在实际应用中，三维模型模拟的计算量很大，计算时间很长。为减少模型网格数量，提高反演速度，对测线进行了16°的逆时针旋转，建立x轴指向正北，y轴指向正东，z轴指向地心的空间坐标系，其中坐标原点为发射源的中点。实测数据的频率范围为0.28～8 000 Hz,共计45个频点，本次反演使用的频率为0.2 Hz、55 Hz、22 Hz、640 Hz、1 280 Hz 5个频点。反演模型沿x方向网格30个，网距40 m，两侧以2.0倍系数各扩边5个网格。y方向网格30个，网距100 m，两侧以2.0倍系数各扩边5个网格。z方向首层厚度20 m，其下各层厚度按系数1.10递增，边界以1.5倍系数扩边5个网格，最终生成的反演网格为40(东西)×40(南北)×30(垂向，不包含空气层)，反演初始模型采用2 000 Ω·m的均匀半空间模型，扩边网格的电阻率同样设定为2 000 Ω·m，误差门限设置为5%的噪声值。

图5给出了迭代反演过程中参数λ和参数rms的变化情况，其中λ和rms在收敛时表征观察数据的拟合程度。从图中可以看出经过651次的反演迭代，λ在收敛时<10-5；拟合差降到了6.86。虽然该拟合差>1，但初始模型拟合差为700，通过反演迭代拟合差已经下降了近2个数量级，说明了反演结果的可靠性。

图5 反演参数λ(a)和rms(b)曲线图Fig.5 Curves of inverse parameter λ(a)and rms(b)

3 结果与讨论

3.1 反演结果

为了反演结果成图效果更好，将电阻率取对数后进行三维成图，图6给出了地下1 000 m±的三维立体反演结果显示。从图中可以看出，地表3个条带状低阻异常带J-D1、J-D2和J-D3，其中J-D1范围最大，位于测区西北角，下部可能赋有矿体；J-D2、J-D3围岩的电阻率较高，根据测区成矿规律，矿体产于大理岩与千枚岩接触带部位，与J-D2、J-D3边缘呈现的中、低阻过渡带较吻合，下部也可能赋有矿体。同时，李建华等[2]在该测区开展了TDIP法测量，其视极化率剖面平面图和视电阻率剖面平面图的结果圈定的异常区域与笔者的三维反演结果地表圈定的异常带吻合，证实了采用有限内存拟牛顿法三维反演的浅部信息是可靠的。

图6 可控源音频大地电磁测深反演结果三维图Fig.6 3D inversion results of CSAMT

由于大理岩和千枚岩两种岩石岩性的极化率相差不大，而电阻率相差较大，因此仅靠TDIP测量达不到找矿的可靠效果，仍需结合可控源音频大地电磁测深结果来综合探讨矿体赋存的有利部位。

3.2 控矿断裂分析

断裂是花牛山铅锌矿区主要控矿因素,断裂带的发育部位，常伴随矿体的产出。为了分析地下断裂F2(花黑滩--双峰山断裂)走向，将电阻率模型沿测线方向切片。根据测区地质资料，断裂F2位于测区西部，只选取了测线L2--L8的切片图。图7中L2(y=10 200 m)、L3(y=10 000 m)、L4(y=9 800 m)、L6(y=9 400 m)、L7(y=9 200 m)、L8(y=9 000 m)测线切片平面图可以明显看出地下地层为二层电性结构，浅部地层电阻率较低，深部地层电阻率较高，说明测区的铅锌矿带埋深较浅。图中断裂F2均有所呈现，L2断面图中断裂F2地表露出位置为x=200 m处，L4断面图中断裂F2地表露出位置为x=0 m处，L6断面图中断裂F2地表露出位置为x=-300 m处，L8断面图中断裂F2地表露出位置为x=-400 m处，综合判断分析，断裂F2走向为北西方向，倾角45°，向南西倾斜，延伸较浅(<500 m)，为次级断裂带。

图7 三维反演电阻率模型沿测线的垂直切片图Fig.7 Vertical slice of 3D inversion resistivity model along survey line

3.3 成矿有利区圈定

由于断裂F2位于测区西部，矿体的产出常在断裂带边缘，因此取出L1、L8、L13测线切片三维图来划分成矿区域。图8中L1测线的北端有明显的低阻异常，异常区延伸至地下400 m±，根据测区地质资料，该区域地层为低电阻率的含炭千枚岩地层，是主要的找矿干扰地层。L8测线中部的地下低阻异常区明显，异常带南北方向长约500 m，东西方向长约800 m，延伸至地下300 m±，与图6地表的低阻异常带J-D2，J-D3吻合。测区地质资料显示该区域地层为中电阻率大理岩夹千枚岩地层，说明浅部异常可能是低阻铅锌矿引起的，为找矿有利区域。L13测线下方地层均为高阻地层，非找矿有利区域。

图8 三维反演模型电阻率切片图Fig.8 Resistivity slice of 3D inversion model

为了更直观地展示铅锌矿带的成矿区域和成矿地层接触带，将三维反演结果进行电阻率划分，图9显示了测区不同岩性的三维立体图像。结合图8找矿的预测结果和测区地面地质资料，图中显示的黄色区域为中电阻率岩层，可能表征千枚岩，绿色区域为中电阻率岩层，可能表征大理岩，测区左上方有明显的低阻层覆盖，推测为含炭千枚岩地层，与图6的三维反演结果吻合。同时，可以推测出测区的西部地区地层浅部表征低阻含炭千枚岩，深部表征中阻震旦系洗肠井群大理岩；中部地区地层表征中阻震旦系洗肠井群大理岩；东部地区地层表征中阻震旦系洗肠井群千枚岩，浅部可能有高阻印支期花岗岩侵入。图中白色区域代表两种岩石的岩性接触带，紫色的低阻异常位于大理岩与千枚岩接触带附近，L8测线之下，地下约200 m±。同时李健华等[2]的结果地表圈定该部位的异常也表现为高极化率，与图6中J-D2条带状低阻异常区对应，推测该低阻异常为铅锌矿带，走向为近东西方向，与已知二矿带对应。断裂F2与铅锌矿带走向相近，位于矿体边缘，证实了该断裂是测区铅锌矿带的控矿构造。

图9 三维地层分布和成矿有利区域电性结构图Fig.9 3D stratigraphic distribution and electrical structure map of favorable regional mineralization

4 结论

(1)使用L--BFGS法进行反演得到的三维电性模型，模型清晰地展示了测区的电性特征，显示测区西北角大范围的低阻异常主要是由含炭千枚岩地层引起的，铅锌矿找矿地层电性特征标志为低、中阻过渡带部位。

(2)根据三维反演结果，探明了断裂F2的情况，圈定出测区内铅锌矿带的成矿有利区域，认为该矿体产于大理岩与千枚岩接触带附近，断裂F2为该矿带的控矿构造。

(3)采用可控源方法探测铅锌矿是可行的，可以将其应用到其他类似的铅锌矿勘探工作中。

致谢本文采用的可控源数据由科研任务(2011YQ05006010)所属国家重大科学仪器研发计划开发的仪器采集。资料处理和建模基于三维电磁数据处理平台EMDesk。感谢吉林金太地球探测技术有限公司提供的三维反演计算资源。