近年来，频率域可控源电磁法在油气勘探、航空近地表调查和寻找隐伏金属矿藏等方面发挥着越来越重要的作用，并使得可控源电磁法反演研究取得了长足进步.近年来，国内外学者对可控源电磁法均做了大量工作，底青云等(2004，2008)，汤井田和何继善(2005)系统介绍了可控源音频大地电磁法国内外进展及针对多维问题进行了正反演研究.林昌洪等(2012)利用视电阻率信息进行了非线性共轭梯度法(NLCG)的反演.翁爱华等(2012)利用场值信息进行了非线性共轭梯度反演.刘云鹤和殷长春(2013)采用NLCG和L-BFGS对航空电磁三维反演问题进行了对比研究，Newman和Boggs(2004)用拟牛顿法进行了电磁数据的反演解释.Egbert和Kelbert(2012)对三维电磁反演框架进行了深入讨论.Alex and er等(2014)通过直接法进行了三维可控源电磁分布式计算及实际资料处理.三维反演的计算效率问题是目前研究的重点.

三维可控源电磁反演过程中，正演计算的速度是影响反演计算效率的关键因素.正演大多采用直接法或迭代法处理大型复系数稀疏方程，针对单个发射源的可控源电磁问题，三维正演采用迭代法具有速度快，占用内存少等优点.为提高三维反演速度，反演大多采用基于线搜索的方法.在线搜索中如何选择搜索方向是反演方法的核心，且不同的搜索方向形成不同的解决方案，主要包括：最速下降法、牛顿法、非线性共轭梯度法和有限内存拟牛顿法.其中，最速下降法较简单，但收敛速度慢；牛顿法收敛速度快，但需要计算二阶偏导数矩阵，计算量大.有 限内存Broyden-Fletcher-Goldfarb-Shanno(L-BFGS)(Haber，2005)是共轭梯度法的一种推广，其本质思想是由给定的初始矩阵开始，利用一阶导数和函数信息近似二阶导数，从而不需要显式地形成Hessian矩阵，相对于利用一阶导数信息的线搜索方法，具有收敛速度快的优点.

本文首先介绍了基于耦合势的三维有限体积正演算法，其次针对L-BFGS算法的反演流程、“拟正演”计算和选择搜索步长问题进行了讨论，最后根据合成数据反演结果进行了L-BFGS与NLCG对比研究，探讨了张量测量的优势，并分析了旁侧效应和阴影效应的影响.

三维电磁正演中，若控制方程中采用的频率比较低或电阻率比较高时，直接对Maxwell方程进行离散，最终形成的大型稀疏复系数线性方程组的条件数往往很大，导致方程很难收敛.

本文建立一套基于耦合势的三维有限体积算法(Weiss and Newman，2003;徐志锋和吴小平，2010). 利用库仑规范条件将麦克斯韦方程变换为与其等价的矢势与标势相耦合的亥姆霍兹方程，并将其在Yee氏交错网格上进行数值离散，最终得到一个复系数稀疏线性方程组.为减少保存方程所需的内存空间，把方程存储为CSR(行压缩)格式，采用不完全LU分解预处理的对称拟最小残差法(ILU-SQMR)求解大型复系数稀疏方程.

设σ是各向同性地层中的电导率，电磁场随时间变化关系为e^－iωt，ω=2πf.在均匀各向同性介质中频率域Maxwell方程为

其中μ₀=4π×10^－7H/m是真空中的磁导率.因为磁感应强度 B 是无散的，它可以表示为一个矢量的旋度，即：

将公式(3)代入(1)中，变形得：

如上(4)式表示 E －iω A 为无旋的，所以它可以表示为一个标量的梯度(Aruliah et al.，2001)，即：

为表述方便，引入变量Φ=(iω)^－1Ψ，得：

将(3)和(6)代入(2)中，得：

对(7)应用矢量恒等式 ▽ ×(▽ Φ)=0和▽ × ▽ × A = ▽(▽ · A)－ ▽ ² A .同时，为了保证矢势 A 的唯一性，应用库仑规范条件 ▽ · A =0，得到矢量Helmholtz方程为

式(8)具有非常稳定的离散形式.在式(7)中，▽ × ▽ ×算子隐式的保证了 ▽ · J =0，这个条件在获得式(8)的过程中丢失了.因此，为保证电流密度所满足的物理条件，对式(7)两边求散度，得到：

为了解决模拟中场源点处的奇异性问题，我们采用将总场分为一次场和二次场之和.即选择一个参考模型σ₀，通过解析方法求出场源 J _s的响应(A _p，Φ _p)，然后利用一次场和参考模型求出二次场(A _s，Φ _s).二次场所满足的控制方程为

其中Δσ=σ－σ₀为模型电导率与参考模型电导率之差.

为验证三维数值解的正确性以及网格设计的有效性，我们利用三维程序计算了一个层状模型的视电阻率及相位响应，并与一维解析解进行比较.层状模型的第一层电阻率为10 Ωm，厚度为400 m，第二层电阻率为100 Ωm，空气的电阻率为2×10⁸ Ωm.如图 1所示.X向的电偶极子源放置于(0，4500 m，0)处，长度1000 m，观测点坐标为(0，0，0).网格剖分为46×70×34，在观测区域和源所在区域，网格是均匀的，间距为50 m.收发距之间采用先大后小的策略，首先将网格间距逐步增加到200 m，保持15个网格，之后逐步减小到50 m，然后以因子3等比增长添加边界网格.利用此网格剖分我们分别计算了从0.01 Hz到8192 Hz共24个频率的视电阻率及相位响应，结果如图 2所示.通过与解析解的对比，电阻率和相位的平均相对误差都小于2%.

图 1

Fig. 1

图 1 三维可控源电磁正演模型Fig. 1 3D model for CSEM forward modeling

图 2

Fig. 2

图 2 三维可控源电磁正演结果验证Fig. 2 Verification of forward modeling results of 3D CSAMT

本文测试计算机基本配置：CPU为Intel I7-4770K，主频3.5 GHz，内存8 GB.一次正演计算用时1小时11分钟(一次场计算耗时52 min).

按照Tikhonov正则化理论，地球物理反演问题可归结为求解以下目标函数的极小值问题(Avdeev，2005):

其中ψ(m)为总的目标函数，ψ_d(m)=‖ V(d - F(m))‖₂² 为数据目标函数项，ψ_m(m)=‖ W(m - m ₀)‖₂²为模型约束目标函数项，λ为正则化因子，使得数据目标函数项和模型约束目标函数项达到一种平衡. d 为观测数据向量，m 为待反演模型电阻率参数矢量，m ₀为初始模型矢量，F 为模型的正演响应向量，V 为数据方差矩阵， W 为模型光滑矩阵(Sasaki，2001)(本文对 W 进行了归一化处理).

目标函数确立后，反演问题就是极小化目标函数的问题.对式(12)的目标函数求梯度得：

根据式(13)中的定义，并令 A 为正演响应向量的雅可比矩阵，则目标函数的梯度可进一步表示为

式(14)中，目标函数及其梯度信息是反演的关键信 息.其主要的计算量集中在求取雅可比矩阵的转置 与一个向量的乘积上，即 A ^T q(Newman and Alumbaugh，2000; Rodi and Mackie，2001)，在附录中简述其求取过程.

针对大规模反演问题，Nocedal(1980)修正了BFGS方法并提出解决非线性最优化问题的有限内存BFGS(L-BFGS)方法.相对于BFGS方法的最大改进就是利用有限次的迭代信息修正Hessian矩阵.L-BFGS反演流程如图 3所示.

图 3

Fig. 3

图 3 三维可控源电磁L-BFGS反演流程图Fig. 3 Flowchart of 3D CSEM L-BFGS inversion

对于α的要求，不能只要求f(m _k+α_kp_k)≤f(m _k)，这样会导致收敛速率过慢.我们采用满足 Wolfe条件的线搜索方法(Jorge and Wright，2006)，公式为

其中0<c₁<c₂<1.式(15)称为Armijo条件，要求α的选值要使f(m _k+α_kp_k)至少减少c₁倍的α_k g _k^Tp_k，可设C₁=10^－4.式(16)称为曲率条件，g ^T_k+1p_k为f(m _k+α_kp_k)对α的导数，使得 m _k+α_kp_k和m _k不会距离太近.其中：

本文反演测试计算机配置同正演配置.反演区域与测区一致，网格剖分方式同2.2节网格剖分方式.

图 4模型包含两个高阻异常体和两个低阻异常 体，埋深100 m，异常体大小为300 m×300 m×300 m，其中，低阻异常体电阻率为10 Ωm，高阻异常体电阻 率为1000 Ωm，背景为均匀半空间，电阻率为100 Ωm. 线电流源长度1000 m，坐标为(-500 m，4500 m，0)到(500 m，4500m，0)，发射频率从0.1～4096 Hz，共15个频率.测网范围为[-800 m，-800 m，0]×[800 m，800 m，0]，共1024个物理测点，观测方式如图 9a所示，反演采用的合成数据为理论数据加入5%的高斯噪声，初始模型为电阻率为100 Ωm的均匀半空间，λ为1.

图 4

Fig. 4

图 4 4个异常体模型Fig. 4 Models with 4 anomalous targets

图 5和图 6分别为NLCG和L-BFGS反演结果在XZ、YZ、XY三个平面上的切片图.如图所示，两种方法的反演结果都较好的确定了异常体的位置，其中，低阻异常体的电阻率比较接近真实电阻率值.需要指出的是，两种优化策略具有相似的内存需求(NLCG为916MB，L-BFGS为952 MB)，但是，从图 7给出的目标函数值随迭代次数的变化可以看出，NLCG方法过早地陷入了局部极小，导致NLCG法在20次迭代以后，目标函数下降很慢，而L-BFGS法具有相当的稳定性.

图 5

Fig. 5

图 5 NLCG三维反演结果(30次迭代)Fig. 5 3D inversion results of NLCG (30 iterations)

图 6

Fig. 6

图 6 L-BFGS三维反演结果(30次迭代)Fig. 6 3D inversion results of L-BFGS (30 iterations)

图 7

Fig. 7

图 7 L-BFGS与NLCG目标函数对比Fig. 7 Comparison of L-BFGS objective function and NLCG objective function

L-BFGS方法采用1作为迭代步长就可以满足线搜索的条件，如图 8a所示.相比之下，NLCG方法需要迭代步长特别小才能满足线搜索的条件.如果从求取目标函数梯度的次数来考虑反演执行效率，如图 8b所示，NLCG法大多需要计算两次目标函数梯度，用时23小时31分，L-BFGS法大多只需要计算一次目标函数梯度，用时12小时17分.因此，L-BFGS方法的执行效率高于NLCG方法.

图 8

Fig. 8

图 8 反演中步长和目标函数梯度随迭代次数变化(a)步长随迭代次数变化；(b)目标函数梯度随迭代次数变化.Fig. 8 Step size (a) and objective function (b) varying with iteration number in inversion

我们设计三组模型，通过两个标量测量算例和一个张量测量算例的对比，在一定程度上说明张量测量的优势.模型包含两个高阻异常体和两个低阻 异常体，埋深100 m，异常体大小为300 m×300 m×300 m，其中，低阻异常体电阻率为10 Ωm，高阻异常体电阻率为1000 Ωm，背景为均匀半空间，电阻率为100 Ωm.反演采用的合成数据为理论数据加入5%的高斯噪声，初始模型为电阻率为100 Ωm的均匀半空间，λ为1.

(1)观测区域如图 9a所示.综合图 6结果，在z方向上，异常体的位置和电阻率控制的较好.然而，由于该观测方式，E_x、H_y分量强，E_y、H_x分量存在弱区，因此，在Y方向上，低阻体被拉长.

图 9

Fig. 9

图 9 观测装置平面图Fig. 9 Plane view of survey apparatus

(2)观测区域如图 9b所示.采用发射源坐标为 平行于Y轴(4500 m，-500 m，0)到(4500 m，500 m，0). 综合图 10结果，在z方向上，异常体的位置和电阻率控制的较好.然而，由于该观测方式，E_y、H_x分量强，E_x、H_y分量存在弱区，在X方向上，低阻体被拉长.

图 10

Fig. 10

图 10 观测装置2的反演结果Fig. 10 Inversion results of the second survey apparatus

(3)观测区域如图 9c所示.采用发射源坐标分别为平行于X轴(-500m，4500m，0)到(500m，4500m，0)；平行于Y轴(4500 m，-500 m，0)到(4500 m，500 m，0).综合图 11结果，发射源采用两个不同的且相距很远的极化.通过张量数据的三维反演结果与标量数据的三维结果对比，张量测量对于恢复真实异常体模型(特别是圈定异常体边界)要好于标量测量.

图 11

Fig. 11

图 11 张量测量反演结果Fig. 11 Inversion results of tensor measurement

在实际工作中，某个区域可能无法布置测网，设计一组模型进行三维反演用以讨论旁侧效应的效果.测试模型为均匀半空间背景中有一个低阻异常 体，如图 12所示.异常体埋深100 m，大小为200 m×200 m×200 m，其中，低阻异常体电阻率为10 Ωm，背景电阻率为100 Ωm.线电流源长度1000 m，坐标为(-500 m，4500m，0)到(500m，4500 m，0)，发射频率从0.1 Hz～4096 Hz，共15个频率.反演采用的合成数据为理论数据加入5%的高斯噪声，初始模型为电阻率为100 Ωm的均匀半空间.正则化因子为1.测网分为3种：

图 12

Fig. 12

图 12 异常体模型Fig. 12 Models of anomalous targets

(1)X从-800 m到800 m，Y从325 m到800 m，如图 13a.

图 13

Fig. 13

图 13 观测方式2平面图Fig. 13 Plane view of survey scheme 2

(2)X从-800 m到800 m，Y从-800 m到-325 m，如图 13b.

(3)X从-800 m到800 m，Y从-800 m到 -325 m；X从-800 m到800 m，Y从325 m到800 m，如图 13c.

首先，在异常体两边分别布置测网，利用旁侧效应进行三维反演，反演结果如图 14和15所示，反演出的异常位置偏离异常体真实位置，并偏向所布置测网的下方.然后，在异常两边同时布置测网，反演结果如图 16所示，反演出的异常位置相对准确，电阻率值也相对于图 14和15有较好的结果.

图 14

Fig. 14

图 14 观测方式一反演结果Fig. 14 Inversion results of the survey scheme 1

同时，我们利用三维反演处理实际数据过程中，阴影效应会影响反演结果，如图 15所示.为了避免阴影效应的影响，从图 16中所得结果，我们应该在测网外增加可控源电磁控制点，使得三维反演具有足够完备的数据，从而减弱阴影效应的影响.

图 15

Fig. 15

图 15 观测方式二反演结果Fig. 15 Inversion results of the survey scheme 2

图 16

Fig. 16

图 16 观测方式三反演结果Fig. 16 Inversion results of the survey scheme 3

本文正演采用基于耦合势的有限体积方法，反演采用L-BFGS方法.首先，通过L-BFGS算法与NLCG算法对比分析，发现三维反演中L-BFGS比NLCG算法的执行效率更高，更适合处理三维大数据量的反演问题.其次，针对测网数据，张量观测相对于标量观测，前者的三维反演结果相对于后者，对圈定异常体边界和反映真实电阻率值具有更好的效果.最后，针对无法布置测点区域，可利用旁侧效应进行三维反演.同时，实际生产中为避免阴影效应影响，应该采用三维反演处理实际数据，并在测网外增加CSEM的控制点.在本文工作的基础上，为提高计算精度，可采用加密网格的方式，然而，这样会带来更多的未知量，导致反演会更加奇异，下一步可考虑将正演模型网格和反演网格分离(Commer and Newman，2008)，进而提高三维反演的稳定性.

致谢 向两位审稿人对本文提出的修改意见表示衷心感谢.

附录A A ^T q 计算

在正演过程中，令

表示由Maxwell方程组离散得到的线性方程组，其中，e 为电场分量构成的向量，s 为源项.为了使反演更加稳定，我们定义正演响应为视电阻率的自然对数，即：

向量 a _i和b _i表示将 e 组合为观测点处的E_x和H_y的插值向量.在此我们认为插值向量和模型无关，则：

其中向量 c _i 定义为

从式(A1)中我们可以得到

将上式代入式(A3)中，得：

根据式(A8)，我们得到：

令向量 r 满足：

则：

在上述的计算过程中，和是可以通过解析方法求出的，且 K 是稀疏的，所以这部分计算量不是很大.除此之外，需要做两次正演，第一次求取 e，第二次通过式(A10)求取向量 r，这是整个过程中计算量最大的部分.这里所谓的正演，指的是求解线性方程组.