0 引言

可控源音频大地电磁法(CSAMT)是在音频大地电磁法(AMT)的基础上发展起来的一种人工源频率域测深方法。20世纪70年代初，加拿大多伦多大学的Strangway教授和他的学生Goldstein提出沿用AMT的测量方式观测人工供电产生的音频电磁场，称这种方法为CSAMT。由于使用人工源^[1-2]，这种方法克服了天然源观测上的一些困难，被先后应用于深部隐伏金属矿体勘探^[3-4]、油气构造勘查^[5]、地热和水文工程^[6]等方面。目前，普遍应用的均为一维、二维反演技术。近些年，随着计算机技术的飞速发展，三维反演技术开始兴起，出现了共轭梯度法^[7-8]、非线性共轭梯度法^[9]、有限内存拟牛顿法、高斯牛顿法^[10]等三维反演方法。其中，由Broyden等^[11-14]在1970年提出并发展起来的基于拟牛顿(BFGS)法的有限内存拟牛顿(L-BFGS)法^[15]继承了牛顿法收敛速度快的优点，同时又克服了牛顿法计算量大、产生非下降方向的缺点，在求解无约束优化问题时缩短了计算时间，提高了工作效率。而在求解实际问题时L-BFGS法又具有超线性收敛性和非常好的自校正能力，是最有效的拟牛顿法之一。2006年Avdeeva首先将BFGS法应用到一维大地电磁反演^[16]，并将这种方法进一步推广到三维反演中^[17]。2012年贾定宇等^[18]采用L-BFGS方法解决了可控源电磁测深一维反演问题。2013年，刘云鹤等^[19]将其应用于3D航空电磁频率域数据反演，对比显示，该方法在反演速度和反演结果的准确性上比非线性共轭梯度法(NLCG)有一定的优势。因此，本文将L-BFGS方法应用到地面可控源电磁勘探三维反演中。不同于赵宁等^[15]的研究，在反演过程中，我们采用自适应正则化因子来平衡数据拟合和模型拟合对反演结果的影响。论文首先用理论模型来验证这种方法的准确性和有效性，然后用该方法反演吉林省集安地区采集的可控源电磁数据，并用钻孔验证发现的异常，以证明本文研究方法的有效性和实用性。

1 有限差分法正演

三维频率域可控源方法正演计算基于如下异常场满足的微分方程：

式中：k_p为背景介质复波数；k为含异常后介质复波数；E_p为背景电场；E_s为异常体产生的异常电场^[20-21]。

采用虚界面法^[22-23]，利用直接积分法计算E_p^[24]。E_s的求解通常采用有限差分法^[25-26]、有限元法^[27]、积分方程法^[28-29]，本文采用的是交错网格有限差分法。将三维介质剖分成正六面体单元，即沿x、y、z 3个坐标方向分别用若干个平行的平面以不同的间距将研究区域划分成若干个小的长方体网格单元。式(1)的矢量微分方程可以分解为3个方向上的微分方程，将离散化后的矩阵方程合并，并且进行对称化处理，得到如下的矩阵方程：

式中：K是复系数大型稀疏对称矩阵；s是式(1)右端的源项。采用拟最小残差法(quasi-minimal residual, QMR)求解式(2)，求得每个边中点的电场值，再根据法拉第定律微分形式计算磁场分量H。本次研究观测数据选择阻抗Z_xy^[30]，而对于可控源数据有

式中：L_E、L_H分别为空间插值算子；E=E_s+E_p。

2 有限内存拟牛顿法反演

设正则化目标函数为

式中：m是模型参数向量；d是数据向量；模型C_d是数据方差矩阵；C_m是模型方差函数；F是正演算子；m₀是初始模型；λ为自适应正则化因子。当两次反演的拟合差变化小于0.002时，λ=λ/10^[31]。

设目标函数二次可微，取目标函数在某已知点x_k处的二阶泰勒展开式：

式中：α为搜索步长；φ_k=φ(x_k)；g_k=g(x_k)为目标函数的一阶导数；p_k表示模型搜索方向；H_k为目标函数的二阶导数，即Hessian矩阵。为了寻找极小化目标函数的搜索方向p_k，令，可得

则牛顿法的模型更新方式为

利用充分下降条件和曲率条件来约束搜索步长α^[32]。在通过式(5)求未知变量的搜索方向p_k时，利用了目标函数二阶导数的信息，从而能快速达到局部收敛。但直接求解式(5)，需要计算和存储Hessian矩阵H并求其逆，这需要很长的计算时间和较大的存储空间。而拟牛顿法改进了这一缺点，它通过构造一个对称正定矩阵B_k来近似Hessian矩阵H_k，将式(6)写成

为了使B_k合理地逼近H_k，要求其满足

式中：s_k=x_k+1-x_k；y_k=g_k+1-g_k。通常是给定一个单位矩阵B₀，然后通过修改B_k得到B_k+1。显然，随着迭代的进行，求取矩阵B_k+1将依赖之前反演的各个B_k，从而要求较大的内存来存储B_k。L-BFGS法在上述拟牛顿法的基础上进行了改正，利用BFGS校正公式，由一个对称正定初始矩阵B₀∈R^n×n开始，利用下面的校正公式计算得出Hessian矩阵的近似：

式(9)改进了拟牛顿法^[1]，可以指定一个正整数m，使其只利用最近m次迭代的信息，即使用向量对计算B_k+1，舍弃早期迭代的信息，从而极大地减少了内存需求，且提高了计算速度^[33]。

图 1给出了有限内存拟牛顿法用于可控源资料三维反演的流程^[19]。从图 1可见首先需要求取梯度g_k。而求取g_k则需要计算雅可比矩阵，对于三维电磁法，采用伴随正演来求取雅可比矩阵。

3 理论模型试验 3.1 模型描述

为了验证反演方法的准确性，我们设计了一个两层介质模型，在埋深1 000 m处存在一个规模为4 000 m×4 000 m×500 m、电阻率为500 Ω·m的高阻异常体(图 2)。图 2a中(-4 000, 10 000)×(4 000, 18 000)为测量区域。线距400 m，点距400 m，共441个物理测点。图 2a中AB表示可控源方法观测时使用的发射源，中心点位于(0, 0)，源的长度为1 000 m，电流为1 A，发射频率从2^-3Hz开始，以倍频增加到2¹³ Hz，共17个频率。观测参数为Z_xy的实部和虚部。

3.2 响应特征

本文采用交错网格有限差分法计算了图 2模型的频率响应。模型剖分的网格区域略大于测量范围。x、y和z三个方向被剖分为40×40×30共48 000个网格，x、y方向的网格宽度均为300 m；z方向的网格第一层厚度为25 m，其下各层的厚度按系数1.10递增，总厚度为2 458 m。三维电磁反演存在着关键频率^[34-35]。在本例中，关键频率约为2^-3、2^-2、2^-1Hz，其对应的上述装置Z_xy响应平面分布模拟结果如图 3所示。从图 3可以看出，Z_xy曲线异常响应明显，实部曲线向远离源的方向偏离，虚部向靠近源的方向偏离，且随着频率的增加，异常的幅值增大。

3.3 反演结果

将理论模拟数据的Z_xy实部和虚部分别施以5%的高斯噪音后作为观测数据进行反演。初始模型为两层介质模型，第一层介质电阻率为30 Ω·m，厚300 m；第二层介质电阻率为100 Ω·m。文中用于计算背景场的背景模型与初始模型一致，且在反演过程中背景模型不进行更新。反演的网格剖分同正演一致，反演在主频3.6 GHz、内存4.0 GB、处理器为Intel Core i5的个人计算机上进行运算。反演迭代111次，每次迭代4 min，共用时6.2 h。图 4为迭代过程中拟合差(RMS)、自适应正则化因子和模型修正量的变化图。从图 4a可以看出：迭代5次后，拟合差迅速下降到3左右，说明了有限内存拟牛顿法的收敛速度较快；随着迭代次数的增加，拟合差下降速度趋于平缓，最终趋于1。正则化因子则在拟合差趋于平缓后开始降低(图 4b)，当两次迭代的拟合差的差值小于0.002时，λ减小为原来的1/10，通过降低模型在目标函数中所占的权重来尽量拟合数据。模型修正量在反演开始时增大速度较快，随着迭代次数的增加，最后也逐渐趋于平稳(图 4c)。

图 5为图 2中异常体模型的反演结果。从横切异常体的平面图(图 5a)可看出，异常分布范围与理论异常的平面位置对应。在纵断面图(图 5b、c)上，异常体的大小和在x、y方向的中心位置与理论模型也基本一致，但在z方向上异常体被拉伸，位置偏下，且异常中心电阻率比理论值低。结合模型特点，产生上述现象是由于覆盖层的电阻率较低，电磁波衰减较快，导致深部信号较弱。此外，从本质上，电磁法对高阻异常不灵敏，这造成深部高阻体的分辨能力下降。

为了测试该反演方法对初始模型的依赖性，我们将初始模型改为电阻率为100 Ω·m的均匀半空间模型，再一次进行反演。反演网格与上文一致。观测数据仍为施以5%高斯噪音的Z_xy实部和虚部。反演迭代了274次，每次迭代4 min，共用时18.2 h。反演拟合差从146.00下降到1.78，拟合效果较好。反演结果如图 6所示。从图 6中可以明显看出：地下介质为层状结构，浅部为低阻，厚度约200 m，真实模型为300 m，反演结果与真实模型较为接近；深部的高阻异常体与第一次反演结果对比，水平位置也吻合较好，但深度明显向下偏移，电阻率值也更低。将两次反演结果进行对比，可以看出，第一次反演结果要优于第二次，但后者也能与真实模型较好吻合。这说明当初始模型与真实模型不接近时，反演方法仍能反演出真实模型。

4 实际数据反演

为了验证该方法的实用性，本文将其应用到实际数据的处理中。采用的数据为吉林省集安市海关地区钼矿勘查的可控源实测数据，工区位于海关市的西南部。该区晚侏罗世在鸭绿江火山-岩浆岩带中南段形成了火山喷发、岩浆侵入相间出现成带分布的火山-岩浆岩带，这一中生界构造-岩浆岩带多金属成矿带内沿鸭绿江构造金、铅、锌、铜等矿化十分普遍，是主要的找矿远景区段。研究区域附近已经发现几个小型钼矿点^[36]，至今没有发现大规模的矿床。本次工作期望通过地面可控源电磁勘探，了解矿区的深部电性特征，寻找隐伏矿化富集地段和含钼矿岩体，为深部钻探施工提供依据。工区内布设了5条北西向测线，距离源由远到近依次为1—5线。线距为200 m，点距50 m，共140个测点，电极距离50 m。观测方式采用赤道偶极方式，发射源采用电导线源，位于工区东北部，平行于测线方向，源长为1 km，距离最近的测线约6.9 km。

4.1 区域地质背景

在构造上，矿区位于华北地台东部，在辽东台隆太子河—浑江陷褶断束中。鸭绿江断裂在区内通过。

区内出露的地层主要为古元古界集安群临江(岩)组，构成凉水—小古马岭背斜。地层总体走向310°左右，背斜北东翼地层北东倾，倾角20°~40°，背斜南西翼地层南西倾，倾角10°~40°。岩性为片麻岩、含墨黑云变粒岩、浅粒岩。

区内脉岩发育，种类繁多，多数为燕山期，尤以燕山晚期活动更为常见，主要有花岗斑岩、细粒闪长岩、闪长玢岩、石英斑岩等，多呈北东向展布。矿体主要分布在中、酸性脉岩发育区，临江(岩)组中的硅化、绿泥石化、绢云母化、碳酸盐化、黄铁矿化组合及中元古代碱长花岗岩中硅化、绢云母化、黄铁矿化组合中。

区内的鸭绿江断裂始于侏罗纪以前，具有多次继承性活动的特点，它控制了燕山期花岗岩、花岗斑岩、闪长玢岩、长石斑岩侵入体的空间展布及中生代火山喷发，控制了区域内生矿产的分布。同时，还发育北东、北西向断裂构造，它们是区内主要的容矿构造。

4.2 地球物理特征

为研究工区不同地层和岩性的电阻率特征，我们搜集了研究区岩石的电阻率信息(表 1)。从结果来看，片麻岩和花岗岩都表现为高阻，蚀变的片麻岩为次高阻，辉钼矿脉为低阻。从表 1中可看出矿脉与围岩电阻率有明显的差异，因此适合开展CSAMT工作。

4.3 数据分析和处理

可控源数据采集频率为2^-3~2¹³ Hz，以2的指数倍递增，共17个频率。图 7为典型测点的频率-视电阻率曲线。从图 7可见，工区的电阻率很高，基本在1 000 Ω·m以上。因此，视电阻率从大致512 Hz开始就进入了近区，呈现45°上升。由于干扰和观测误差的存在，在频率-视电阻率曲线上会出现非正常的突跳，如202号测点，频率大于100 Hz时曲线突跳剧烈，只能观察出曲线的大致规律。114号测点数据质量较好，频率-视电阻率曲线平滑。

反演针对阻抗Z_xy进行，单位采用实际观测中的常用单位^[37]，对应的复阻抗频率响应如图 8a所示。由于阻抗的振幅和相位有一定的内在联系^[38]，我们提出借助Nyquist图对误差数据进行必要的处理^[39]。图 8b为对应的Nyquist图，从图中看出曲线从低频到高频逆时针旋转，不符合我们经过大量数值模拟和实际数据得出的其沿顺时针旋转的规律，需要对其进行编辑。图 9为阻抗编辑完成的阻抗数据。

4.4 反演参数

实际数据在反演之前，需要进行必要的圆滑。由于高频数据质量较差(图 7a)，我们只截取0.125~256. 000 Hz频率范围内10个频率的数据。反演数据为阻抗分量Z_xy。由于实际数据误差较大，我们将阻抗振幅的5%作为数据误差，让反演更加有效地进行。

反演初始模型参数包括网格大小和初始背景电阻率。网格大小决定了计算的精度，网格越小，计算的精度越高，但计算时间越长。让两个测点中间有两个网格节点，既保证了计算精度，又节约计算时间。如果测点较多，两个测点中间有一个网格节点也可。本次研究点距为50 m，线距200 m，所以将网格大小定为25 m×100 m。z方向的网格大小根据勘探深度的需要进行调整，本文目标深度为1 km，所以第一层厚定为10 m，层厚按照系数1.1依次递增。

初始模型的背景电阻率对反演影响很大。这里根据研究区的地质资料和岩石的地球物理特征，将初始模型的电阻率选为5 000 Ω·m。反演在CPU为3.60 GHz主频的电脑上运算，内存4 Gb。迭代一次7 min，迭代85次，用时650 min。反演拟合差从12.9降到1.53。

4.5 结果及讨论

我们实际感兴趣的是高阻岩体及与之相关的低阻矿化蚀变。图 10为CSAMT三维反演获得的高阻和低阻分布立体图。由图 10可见，电阻率整体以中高阻为主。工区的东部有一条沿北东向延伸的低阻条带，倾向为北西，推测为鸭绿江断裂(图 10中F₁所示)。在该断裂的深部有一处高阻带向上隆起，电阻率大约为10 000 Ω·m，推测为花岗斑岩，是与钼矿床成矿有关的岩石。推测成矿过程是经过热液的上涌或侵入，将矿质携带到一定的部位，并导致围岩发生蚀变。蚀变带是找矿的重要标志。

研究区内蚀变主要有硅化、绿泥石化、绢云母化、碳酸盐化、黄铁矿化等。由表 1岩石的电性特征可知，矿化蚀变后的岩石表现为中低阻，因此中低阻带可作为找矿标志。为了突出异常的特征，将高阻的岩浆岩和中低阻带单独画出(图 10b)。从图 10b中可以看出，低阻带为一条由下向上延伸的斜条带状，延伸至-500 m左右，倾角约45°。

根据花岗岩和低阻带的位置，推测该低阻为矿化蚀变带，并在(-100, 0, 0)处进行钻探验证，图 10b中黑色线为ZK0806钻孔位置及井身空间形态。钻探揭示，在178 m深度上存在厚度90 m的矿化蚀变，蚀变倾向北西，倾角约50°。将钻探剖面与过井三维反演模型的电阻率剖面进行对比(图 11)发现，钻探结果与可控源异常的位置具有非常好的一致性。

5 结论与建议

本文将有限内存拟牛顿法应用到地面三维可控源电磁勘探中，得出如下认识：

1) 数值模拟结果与理论模型一致，证明了研究方法的准确性和可行性，并且有限内存拟牛顿法收敛速度快，为三维地面可控源电磁反演的实际应用提供了基础。

2) 三维反演结果中推测了一个潜在的钼矿矿脉，并与钻探结果基本一致，证明了文中反演方法的可行性和实用性。

3) 采用可控源方法来探测蚀变钼矿是可行的，可以将其应用到其他类似类型的钼矿探测工作中。