0 引言

可控源电磁法是电磁勘探的一种重要手段。通常情况下是以有限长导线为发射源，观测电磁场响应数据并计算阻抗和视电阻率响应，从而进行资料的分析和解释^[1]。很多算法已应用于一维可控源数据响应的反演中，就线性迭代反演算法来说，最传统的是最小二乘法，但其收敛稳定性和分辨率较差。尚通晓等^[2]采用阻尼最小二乘法 (或称马奎特反演法) 对可控源卡尼亚视电阻率做一维反演，通过调节阻尼因子提高了解的稳定性。刘颖等^[3]采用高斯-牛顿法实现海洋可控源电磁场的一维反演，将非线性正演响应函数线性化，忽略了计算复杂的二阶导数项。何梅兴等^[4]和汤井田等^[5]先后用Occam法进行可控源一维反演，较常规的高斯-牛顿法，Occam法使得反演模型简单光滑。在上述各种模型约束反演中，如何设置最优正则化参数以平衡数据误差与模型约束，是影响反演结果的重要问题。王堃鹏等^[6]以一维可控源全资料反演为例，采用陈小斌等^[7]提出的自适应正则化方法，加快了反演速度。

但上述方法在线性化过程中忽略了正演响应的二阶导数项，损失了电磁响应函数的非线性信息；牛顿法虽然保留了二阶项，但因需要计算复杂的海森矩阵及其逆矩阵，计算量很大，所以其在可控源反演中均没有得到广泛应用。拟牛顿法对该方法进行了改进，通过构建海森矩阵逆的近似矩阵的方式减小计算量，并且较一阶迭代的梯度法收敛更快，但对内存需求很大。Nocedal等^[8-9]进一步提出了有限内存拟牛顿法，只需储存几步最新的数据信息即可构建近似矩阵，实现反演迭代，极大节省了内存，并且反演更稳定，因而适合于求解大规模优化问题^[10-11]。Newman等^[12]和Haber^[13]将该方法相继用于解决大规模的电磁反演问题。Avdeev等^[14]已成功将其运用到三维大地电磁反演中。翁爱华等^[15]也实现了三维可控源有限内存拟牛顿的反演。Avdeeva等^[16]验证了有限内存拟牛顿法在大地电磁一维反演中的可行性。

由于有限内存拟牛顿反演方法的上述优点，本文尝试将该方法应用于有限长导线源一维频率响应数据反演中。首先给出有限长导线源在均匀层状导电地表产生的电磁场虚界面法表达式，得到一维模型正演响应；然后较详细地给出有限内存拟牛顿反演理论技术方案；在此基础上，对理论模型和实际数据进行反演计算，讨论该方法的反演效果和可行性。该方法的实现为可控源资料一维反演提供新的技术手段，同时也为二维和三维反演建立较为合适的初始模型^[17-19]提供一个新的途径。

1 一维正演理论

本文讨论的反演参数为阻抗。通过位于地表水平有限长导线激发的一维均匀层状导电介质的表面阻抗由水平正交电场E_x与磁场H_y定义为

电、磁场分量E_x和H_y的计算采用虚界面法^[20-21]，具体计算模型如图 1所示。发射源Tx长度为2L，位于第z_ls层，沿x轴方向发射，其中心点在地面的投影设为坐标原点O。这里认为介电常数ε和磁导率μ与空气中的相等。

刘云鹤等^[22-25]给出了虚界面法计算接地有限长导线在任意层的电磁场的详细推导过程。当观测点在地表时，正交水平电场E_x与磁场H_y的表达式为:

式中：u₁²=λ²-k₁²，k₁²=-ẑ₁ŷ₁；ẑ₁=iμω，ω为角频率；ŷ₁=σ₁+iωε；；λ为波长；R=[(x-x')²+y²]^1/2，为观测点到源中心的距离；R₁、R₂分别为观测点到源两个端点的距离，R₁=[(x-L)²+y²]^1/2，R₂=[(x+L)²+y²]^1/2；I为供电电流；(x, y, z) 为观测点坐标；x'为源横坐标；F₁⁺，F₁^-，A₁⁺，A₁^-分别为TE、TM模式的正、负向传播振幅系数; J₀(x) 和J₁(x) 为别为0阶和1阶贝塞尔函数。式 (2)、(3) 中的Hankel积分采用数值滤波法^[26]或者直接积分法^[27-28]计算。对于有限长导线，在收发距离较大时，通常采用偶极子等效近似；但考虑到可控源观测可能在近区进行，式 (2)、(3) 中直接对有限长导线进行积分，采用自适应梯形积分方法进行计算^[22]。

虚线为过发射源引入的平行于层界面的虚拟界面。若源位于地表，则空气层作为一层特殊介质层。

2 有限内存拟牛顿反演 2.1 基本原理

反演直接对考虑收发距离的有限长导线产生的阻抗频率响应进行。令N个观测阻抗数据构成的数据向量d=(Z₁, Z₂，…，Z_N)^T，考虑模型光滑性约束，定义目标函数为

式中：φ(m) 是关于模型参数的正则化目标函数；m为迭代过程中产生的模型参数向量；F(m) 为模型的正演响应函数；‖Wm‖²为模型约束；W为模型加权矩阵；ξ为正则化因子，是模型约束的加权参数。

假设目标函数二次可微，在已知模型参数m_k处作二阶泰勒展开，有

式中：g_k为目标函数的一阶导数 (梯度向量)；G_k为二阶导数 (海森矩阵)；p_k为模型的搜索方向；k为迭代次数。为了找到极小化目标函数的搜索方向，令，则

在式 (6) 的迭代过程中，拟牛顿法为避免直接计算海森矩阵的逆，用一个正定对称矩阵H_k近似等于G_k^-1。给定正定对称矩阵H₀，H_k+1通过

逐步迭代计算出，并满足拟牛顿条件：

式中：s_k=m_k+1－m_k；y_k=g_k+1－g_k。采用式 (7) 计算H_k+1需要依赖前面所有的计算结果；而有限内存拟牛顿对此进行了改进，只需在迭代点m_k处给定正整数m(根据实际情况选定m值，一般取3~20) 和初始正定对称矩阵H₀，通过对H₀进行m+1次修正得到H_k+1^{[8-9, 29-30]}，避免了储存大量的矩阵来计算H_k+1。改进过程大致如下。将式 (7) 变换为

式中，I为单位矩阵。若记，V_k=(I－ρ_ky_ks_k^T)，则有

将H₀代入式 (10)，经过k+1次迭代有

当k+1 > m时，忽略旧的信息{s_i, y_i}_i=0^k－m－1，得到

由式 (12) 可见，有限内存拟牛顿法只需储存m+1个向量{s_i, y_i}_i=k－m^k和1个H₀矩阵，即可求出H_k+1。

2.2 算法步骤

1) 给定初始模型m₀、初始对称正定矩阵H₀，以及需要存储的迭代非负整数m，设定导数精度范围ε > 0。

2) 计算导数g_k，如果‖g_k‖≤ε，则算法终止，得到反演解m_k；若不满足，则令搜索方向p_k=-H_kg_k。

3) 确定搜索合适步长α_k > 0，并令其满足Wolfe-Powell条件：

式中：δ和σ为常系数，。令m_k+1=m_k+α_kp_k，通常先取α_k=1。

4) 令。根据式 (12)，利用向量对H₀迭代次，确定H_k+1。

5) 令k=k+1，重复步骤2)。

3 数值结果

理论模拟时，假设源长为1 000 m，供电电流为1 A，发射频率从2⁰Hz按的倍数逐渐增加到2¹³ Hz，源中心坐标为 (0, 0, 0)，观测坐标为 (1 000, 10 000, 0)。模型反演数据为阻抗理论数据和10%高斯噪音理论数据。反演过程中，正则化因子ξ初始值为2，采用目标函数自适应调整技术^[24]。反演时，模型剖分为17层，表层厚度为20 m，其余各层厚度按照1.1倍等比增加。反演从均匀半空间模型初始，迭代停止条件为目标函数相对变化小于10^-4。

3.1 层状模型反演

理论三层模型参数如表 1所示。反演模型为初始电阻率为100 Ω·m的均匀半空间，无噪声数据迭代38次结束，噪声数据迭代34次。反演结果如图 2(H型) 和图 3(K型) 所示。两图综合可见：1) 反演迭代前期收敛较快，目标函数能下降2个数量级，继续反演收敛变慢，而噪声的出现导致反演迭代后期目标函数基本不下降，反映目标函数受数据噪声所控制；2) 从数据本身拟合情况看，迭代结束后，无论噪声数据还是无噪声数据，反演得到的模型理论阻抗响应数据与正演的阻抗响应数据基本重合，说明有限内存拟牛顿反演方法具有较强的抗干扰能力；3) 比较反演模型与理论模型可见，反演结果恢复了地层电阻率随深度变化的真实规律，由于采用光滑模型约束，因此反演的分辨力不高；4) 虽然反演模型的初始模型与理论模型相去甚远，但反演的结果还是反映了真实模型，没有受到初始模型的影响，说明有限内存拟牛顿反演对初始模型的依赖性不强，从而为有限内存拟牛顿反演用于可控源实际数据的反演奠定基础。

为了考察有限内存拟牛顿反演方法的适应性和分辨力，设计了一个更为复杂的HKH型5层地层模型 (表 2)。该模型的最大特点是在深部2个低阻层中，存在一个相对高阻的地层。该模型的反演难点在于恢复低阻层所夹的高阻层；因为按照电磁勘探概念，电磁法对高阻反应弱，对低阻灵敏^[1]。在反演模型结构的设计上，将理论模型中的浅部3层在初始模型设定中细分为第1层厚20 m、下层层厚按1.1倍逐渐递增共17层，深部2层厚度不变，即初始模型为19层。各层反演的初始电阻率均为500 Ω·m，反演迭代35次结束。图 4给出了HKH型5层地层模型的具体反演结果。从图 4a可见，反演电阻率能较好地分辨出地下5层电阻率结构，特别是低阻层中的高阻层; 说明有限内存拟牛顿反演方法具有较高的分辨能力。图 4b显示，反演结束时，曲线拟合非常好; 这也反映在图 4c的目标函数变化曲线上，因为在迭代后期，目标函数小于10^-3，说明数据误差已经非常小。

3.2 实测数据反演

实测数据来自新疆哈密市白山泉铁矿区可控源电磁探测。工区岩性简单，中低阻的浅部第四系砂土和第三系黏土覆盖在深部高阻的花岗岩或者酸性英安岩和构造角砾岩上，低阻的磁铁矿侵入到石英片岩中，如图 5a的岩性柱状图所示。可控源采用有限长导线发射，源长为1 300 m，源中心坐标设为 (500, -5 000, 0)，发射频率为2⁵~2¹³ Hz (以的倍数递增) 共17个频率；观测点坐标为 (475, 500, 0)，收发距离约5 500 m。工区接地条件好，人文干扰小，数据质量较高，具体测量的阻抗振幅数据如图 5b中实线所示。因为已有的钻孔资料只有200 m左右，因此为了对比，在对数据反演时，将地下剖分成浅部四层5 m厚、向下每层10 m厚的模型，共32层，总厚度300 m，各层反演的初始电阻率均为50 Ω·m。反演迭代32次结束，目标函数的变化如图 5c所示。可见，随着迭代进行，目标函数逐步减小，表明一维有限内存拟牛顿反演能稳定收敛。在迭代结束时，模型阻抗响应与实际阻抗数据基本拟合 (图 5b)。反演模型 (图 5a) 电阻率曲线 (为方便对比，将反演结果截取到250 m) 显示为KH型特征：浅部电阻率较低，约100 Ω·m，随深度增加逐渐增大，在深度75 m附近，达到近500 Ω·m；然后减小，反映深部存在较厚或者良导电性异常；上扬的深部电阻率为高阻基底的反映。比较测点附近的岩性柱状图，该电阻率变化规律符合地层情况，两者具有较好的一致性。

4 结论

本文研究了有限内存拟牛顿法应用于有限长导线频率电磁测深资料反演中的可行性，给出了可控源资料反演的一种新途径。在可控源测深反演中，针对阻抗数据进行反演；有限长导线的阻抗频率响应计算由水平正交的电磁场计算；采用虚界面法计算水平均匀层状介质水平电偶源产生的电磁场并沿导线积分得到有限长导线的电磁场。利用层状理论模型和实测数据对本文提出的反演方法进行了检验。反演结果表明：

1) 一维有限内存拟牛顿法用于可控源阻抗数据的反演，反演初期目标函数收敛较快，以后逐渐趋于稳定，反演结果能恢复地下电性变化，响应曲线拟合较好；因此，一维有限内存拟牛顿法能用于可控源阻抗数据的反演。

2) 一维有限内存拟牛顿法对初始模型的要求不高，但由于该方法反演迭代次数较多，和其他反演方法一样，尽可能好的初始模型将对提高反演的计算效率有较大帮助。

3) 对噪声数据和实际数据的反演结果表明，有限内存拟牛顿反演方法对数据干扰具有较强的压制能力，从而为该方法的实际应用奠定了理论基础。

4) 有限长导线频率电磁测深资料的有限内存拟牛顿法反演的实现，为可控源资料一维反演提供一个新的途径，也为二维或者三维可控源电磁数据的反演提供一种建立合适初始模型的新手段。