1 引言

磁共振技术(Magnetic Resonance Sounding，MRS)是目前世界上唯一一种能直接定性、定量分析地下水含量的地球物理探测方法(戴苗，2008；林君，2010；林君等，2012)，无需钻探即可获取地下水分部信息(Semenov，1987；Legchenko and Valla, 2002).然而地电结构是影响MRS探测过程中激发磁场传播的主要因素，不准确的电阻率分布信息会导致MRS反演结果偏离实际(Schirov et al., 1991；Mohnke and Yaramanci, 2000；Yaramanci，2000；Roy and Lubczynski, 2003；Lehmann-Horn et al., 2011).

磁共振与电磁探测方法联合反演是提高含水层解释结果准确性的有效途径之一(Braun and Yaramanci, 2003；Braun et al., 2006).国际上，Hertrich和Yaramanci (2002)最早利用垂直电测深法(Vertical Electric Sounding，VES)获取地下电阻率信息，实现了MRS与VES数据联合反演.然而VES方法探测时需要铺设很长的供电导线以获取更大的探测深度，同时受地形地貌的影响较严重.Braun和Yaramanci (2008)通过理论仿真利用MRS相位反演得到了地下电阻率的分布.然而，在实际测量时受MRS探测系统制约，难以获取准确的相位信息.瞬变电磁法(Transient Electromagnetic，TEM)利用不接地的回线发射一次磁场，在一次磁场的间歇期间接收由地下良导体受激励引起的涡流所产生的二次感应磁场，与其他电法测深技术(直流电法，垂向电测深)相比，它具有高效率、高分辨率、大探测深度等优点，能够获取地下几百米范围内的电阻率信息，因此是与MRS进行联合探测的最优选择(Raiche，1987；Nabighian，1992；薛国强等, 2007, 2008；刘桂芬，2008；齐彦福等，2015).Legchenko等(2009)率先实现了磁共振与瞬变电磁的联合解释，对比传统磁共振解释方法，解释精度显著提高.国内，林婷婷等(2013)通过磁共振与瞬变电磁联合解释，获取了内蒙正镶白旗地区的含水层准确分布(Lin et al., 2013).李貅等(2015)利用瞬变电磁方法进行快速电阻率成像，将得到的电阻率结果作为磁共振三维反演的电性模型，开展磁共振与瞬变电磁三维联合解释.然而上述方法，均是将电阻率作为先验信息，不能实现联合反演中的含水量与电阻率的实时迭代与更新.万玲等(2013)在国内提出将自适应遗传算法(AGA)应用到磁共振与瞬变电磁单一测点联合反演中.但随着研究的深入，作者发现，传统的MRS-TEM联合反演忽略了地下含水结构的横向连续性，会降低反演结果的准确度.横向约束方法最初被应用在连续垂向电测深法(CVES)的数据解释中，通过对大地模型中的电性参数施加横向约束条件，反演得到更精确的电阻率信息(Auken and Christiansen, 2004).后经不断改进，被广泛应用于航空电磁数据处理中(蔡晶等，2014；殷长春等，2016).Behroozmand等(2012)采用延伸型指数函数(Stretched-Exponential，SE)近似估计MRS信号中复杂的弛豫分布，通过横向约束方法实现了磁共振与瞬变电磁的联合成像.但联合反演中涉及到的快速正演迭代以及横向约束具体实现的核心技术文中并未阐述.

本文依据相邻测点地下电性结构具有连续性思想，基于横向约束方法，重点解决复杂地电结构探测时解释结果横向不连续的问题.首先引入外推积分法，解决了传统正演过程中基于直接积分方法导致的求解效率低的难题；其次利用LCI思想，建立了MRS与TEM横向约束联合反演目标函数，通过正演迭代实现了磁共振与瞬变电磁的横向约束联合反演.然后，利用理论模型，对比传统的无约束联合反演与LCI联合反演的数据处理效果，验证LCI联合反演方法的有效性和准确性.最后，通过处理安徽黄山野外实验数据，进一步验证LCI联合反演方法的实用性.

2 地下MRS正演计算方法

磁共振地下水探测基本原理如图 1所示，磁共振方法是在置于地表的线圈中通以I=I₀cos (2πf_Lt)的交变电流，产生的激发磁场使得地下水中的氢质子从低能级跃迁至高能级；当交变电流关断后，质子回到低能级状态，由接收线圈探测携带地下水信息的磁共振信号(孙淑琴等，2005).

同一电流脉冲矩只能探测到特定范围内的磁共振信号，因此，只要发射一系列从小到大变化的电流脉冲矩就能获得由浅至深的地下水分布信息.电流脉冲矩记为q，则有：

(1)

其中I₀和τ_p分别是瞬变电流和脉冲长度.

磁共振信号的初始振幅为：

(2)

(3)

E₀是信号的初始振幅(林君，2010)，K(q, ρ, z)是核函数，其中w(z)是地下含水量分布，ρ(z)是地下电阻率分布，φ_L是角频率，M₀是磁矩；γ_p是氢质子的旋磁比；B_T⁺和B_R^-是垂直场B_T^⊥由于椭圆极化而产生的圆极化场顺时针和逆时针旋转分量，B_T^⊥=B_T⁺+B_T^-；ζ是椭圆极化场的相位参数(Hertrich，2008).

在层状大地模型中建立柱坐标系，电场只有E_φ这一个分量，而激发磁场在第j层的点P(r, z)处的磁场分量有B_rj和B_zj这两个分量，公式为

(4)

(5)

其中，j=0, 1, 2, …, n, u_j=, a_j和b_j为待定系数, a为发射线圈半径，J₀、J₁分别为零阶和一阶贝塞尔函数.且B_T^⊥与B_rj和B_zj有如下的关系：

(6)

本文运用外推积分法来计算上述磁共振激发场. QWE方法首先将连续积分离散化，用高斯积分法修改积分横坐标，然后预先计算并存储与变量参数无关的部分积分，从而提高整体积分的计算效率，最后发散积分部分采用Shanks变换递推计算(Kerry，2012). QWE方法具体运行过程如下.

将激发磁场的z分量，也就是公式(5)写为：

(7)

其中：

(8)

(9)

间隔断点k_i-1 < k_i是一阶贝塞尔函数的零点.公式(9)中的部分积分用高斯积分法计算，得到：

(10)

式中m是积分阶数，w是与积分横坐标x相关的权重.由于w_jJ₁(x_j)与参数a是相互独立的.因此，预先计算w_jJ₁(x_j)并存储，方便调用.

通过ε算法进行外推，ε算法基于Shanks变换的计算流程，定义如下：

(11)

(12)

(13)

式中S_n是部分积分的直接求和，公式为

(14)

将S_n^*定义为ε算法的计算结果，当n是奇数时S_n^*=ε_n⁽⁰⁾，当n是偶数时S_n^*=ε_n⁽¹⁾.当计算结果满足下面的条件时，迭代计算结束, 其中的θ是在程序开始时设置的精度值,

(15)

实际应用时为方便布线，常采用方形线圈代替圆形线圈，数据解释时再运用等效面积法近似计算大地磁场.因此本文在理论模型构建时，为接近实际，发射线圈采用方形回线.三层大地模型验证外推积分法在磁共振正演计算中的高效性和准确性，拉莫尔频率2282 Hz，发射线圈100 m×100 m单匝，地磁倾角60°，各层的具体参数如表 1所示.

在上述三层磁共振大地模型正演时，传统数值积分法计算用时为598.7 s，将近10 min；而QWE方法的计算用时仅仅为7.3 s，提速了大约80倍.在均匀半空间模型下，传统数值积分法用时大约8 min，而QWE方法仅用4 s就完成一次正演计算.图 2是该三层模型的正演结果E₀-q曲线图，为对比外推积分法和传统数值积分法的计算精度.图中结果表明，在峰值处二者趋势基本一致，而随着脉冲矩衰减过程中有细微偏差.

3 MRS与TEM横向约束联合反演方法 3.1 联合反演目标函数

用对数形式表示任意测点的磁共振与瞬变电磁观测数据为

(16)

其中，d_obsi表示第i个测点的勘测数据，V_TEMi和V_MRSi分别表示第i个测点的瞬变电磁和磁共振响应值.

正演模型参数由该测点地下各层的视电阻率、层厚度和含水量构成，对模型各参数取对数(刘斌等，2012)，则第i个测点的模型参数可表示为：

(17)

其中，m_i表示第i个测点的模型参数, ρ_{i, j}、h_{i, j}和w_{i, j}分别表示第i个测点第j层的视电阻率、层厚度以及含水量，N表示模型分层数.

联合反演的目标函数由瞬变电磁与磁共振两部分数据共同构成，建立加权联合目标函数为

(18)

其中，Φ是联合反演目标函数值；V_TEM^obs是TEM观测数据；A_TEM是TEM正演得到的理论数据；V_MRS^obs是MRS观测数据；E_MRS是MRS正演得到的理论数据；D_α和D_β分别是TEM和MRS数据或噪声的不确定度；φ_TEM和φ_MRS分别是根据噪声水平引入的TEM和MRS数据反演加权系数；C是平滑度矩阵；m是联合反演模型；λ是正则化参数.为兼顾反演结果的稳定性、计算精度与计算效率，λ可根据公式(19)调整，公式为

(19)

(20)

(21)

其中，λ^k+1是第k+1次迭代的正则化参数(万玲，2013).

3.2 光滑度约束矩阵

为了减少相邻测点之间的电阻率、层厚度与含水量参数的差异，在反演过程中引入模型光滑约束.根据两点光滑约束方程可得到：

(22)

等式两边减去R_pm₀，得到：

(23)

其中，Δr_p=-R_pm₀；e_rp为相邻两测点之间模型参数的变化；R_p为光滑度约束矩阵，公式为

(24)

利用三点光滑约束方程减少相邻三个测点之间地电参数的差异，能够增强光滑约束强度，公式为

(25)

3.3 层厚度约束矩阵

对纵向各层进行视深度约束，能够达到光滑层界面以及增加各层之间连续性的目的，公式为

(26)

等式两边减去R_tm₀，得到：

(27)

其中，Δr_t=-R_tm₀；e_rt为相邻两个测点之间层界面深度的变化；h_{i, j}和t_{i, j}分别为第i个测点第j层的厚度和层底深度；R_t为层厚度约束矩阵，表示为：

(28)

3.4 加权约束矩阵

为避免约束过度导致反演边界模糊，反演结果偏离实际(翁爱华等，2008)，引入加权约束矩阵，公式为

(29)

(30)

其中，R′ _p=αR_p；Δr′_p=-αR_pm₀；R′_t=βR_t；Δr′_t=-βR_tm₀；α和β分别为两个约束方程的加权系数，通常α和β的取值范围为[0.1, 10].

3.5 整体反演方程建立及求解

根据联合反演数据d_obs与模型参数m之间存在的复杂非线性关系，对数据拟合方程进行一阶泰勒展开，公式为

(31)

其中，d_obs是联合反演数据；F是联合响应函数；m₀是模型初值；m是联合反演迭代模型；G是雅克比矩阵；e_obs是截断误差.移项后可得：

(32)

其中，Δd_obs=d_obs-F(m₀)；Δm=m-m₀；G=diag (G₁, G₂, …, G_M)为各测点的雅克比矩阵按对角线排列后的新雅克比矩阵，公式为

(33)

则任意测点的雅克比矩阵可以表示为

(34)

其中，G_{i, j}为第i个测点第j个模型参数的雅克比矩阵；V_iz为第i个测点的垂直磁场模型正演值；m_{i, j}为第i个测点第j个模型参数；G_i为第i个测点的雅克比矩阵；N为模型分层数.

联合数据拟合方程与约束方程，得到总体反演方程为

(35)

简化后可得：

(36)

为了选取最合适的Δm，引入残差平方和函数为

(37)

当Δm=Δm^*时，S(Δm)取最小值，记作：

(38)

通过对S(Δm)进行微分求最值，可以得到：

(39)

其中，C^-1为协方差矩阵C的逆矩阵，C表示误差矩阵中各项的协方差矩阵C_obs、C_prior、C_rp和C_rt组成的对角矩阵，公式为

(40)

则公式(39)的解即为最终需要求得的模型修正值，公式为

(41)

在系数矩阵中加入阻尼因子λ²，将系数矩阵变为非奇异矩阵，使反演计算的解唯一确定，则有：

(42)

其中，I是单位矩阵.

将矩阵A进行奇异值分解(SVD分解)，即A=UΛV^T，代入公式(42)后可得：

(43)

模型修正迭代式可表示为

(44)

其中，m^k+1是第k+1次的迭代模型.

4 算法有效性验证 4.1 理论模型验证

为了验证磁共振与瞬变电磁横向约束联合反演计算方法的有效性，设计了三层和四层大地模型，分别用传统一维无约束联合反演和基于横向约束的联合反演两种方式对模型进行反演，获取电阻率和含水量分布信息，并与理论模型进行对比分析.

三层大地模型如图 3a和图 3b所示，低阻含水体从中间向两端逐渐变大，层边界清晰.具体参数设计如下：TEM采用100 m×100 m单匝发射线圈，发射电流和接收线圈等效面积全部归一化；MRS采用150 m×50 m单匝发射线圈和50 m×50 m单匝接收线圈，地磁倾角60°，地磁偏角0°，拉莫尔频率设为2282 Hz，激发脉冲矩设为20组，分别从0.1-20 As按指数分布排列；排布10个测点，相邻测点距离为10 m；测区第一层是电阻率为400 Ωm的地表高阻层，第二层是电阻率为100 Ωm的低阻含水层，第三层是电阻率为500 Ωm的高阻基底，含水量分别为10%、30%和5%.

初始模型随机设置，进行10次迭代反演得到反演结果如图 3所示：传统一维无约束联合反演结果如图 3c和图 3d所示，在第9号点和第10号点处出现错误异常，使得测线反演结果横向不连续，其他测点的层边界位置清晰；适当调整电阻率、含水量以及深度约束权重，反演结果如图 3e和图 3f所示，层边界位置清晰，有效地抑制了参数突变.

四层大地模型如图 4a和图 4b所示，低阻含水体沿测线方向呈阶梯状变化，层边界清晰，各层的层厚度不变，层深逐渐增加.硬件参数与环境条件与三层模型一致，排布20个测点，测区第一层是电阻率为500 Ωm的地表高阻层，第二层是电阻率为300 Ωm的浅层低阻含水层，第三层是电阻率为200 Ωm的深层低阻含水层，底层是电阻率为100 Ωm的低阻基底，含水量分别为5%、30%、20%和10%.

初始模型随机设置，进行10次迭代反演得到反演结果如图 4所示：无约束联合反演结果如图 4c和图 4d所示，在1-4号测点处出现错误异常，破坏了测线整体反演结果的连续性，18-20号测点的第三层深度明显与理论模型不符，其他测点的层边界和断层位置清晰；适当调整约束权重的反演结果如图 4e和图 4f所示，有效地抑制了参数突变，电阻率和含水层阶梯状变化明显.

4.2 实测数据验证

应用本文所述磁共振与瞬变电磁横向约束联合反演方法对安徽黄山滑坡实验数据进行联合反演，检验算法的实用性.

通过调查发现，滑坡后缘高程690 m左右；前缘剪出口高程约625 m；滑坡顺主滑方向长150 m左右，滑坡体宽约40 m，平均厚度10 m左右.

野外工作中，瞬变电磁探测采用25 m×25 m单匝发射线圈，用等效面积为1760 m²的多匝小回线接收；磁共振探测采用75 m×25 m单匝发射线圈，25 m×25 m双匝接收线圈；测线方向自西向东排布17个测点，点距10 m；测试点地磁场强度48403 nT，地磁倾角为68.25°，拉莫尔频率为2061 Hz.

设置四层初始模型进行联合反演，高程解释结果如图 5所示：6-15号测点位置的地表至地下10 m (含水量大于30%)以及地下20~30 m (含水量大于20%)为低阻层含水区，地下10~20 m含水量为18%-8%，此处含水量值较小.根据反演结果推测，滑体位置位于6-15号测点地下10~20 m处，6号点为滑坡剪出口，分界面清晰.

5 结论

本文详细介绍了磁共振与瞬变电磁横向约束联合反演的具体实施方法.依次建立了光滑度矩阵、加权横向约束矩阵和纵向层厚度约束矩阵，与数据目标函数共同构成联合反演目标函数.与传统的一维无约束联合反演方法相比，横向约束联合反演可增强解释结果剖面参数的连续性.此外，权系数的选择也在实际数据处理时发挥了重要的作用.为了将横向约束联合反演解释结果与传统的一维无约束联合反演解释结果进行对比，设置了不同的多层大地电阻率及含水量模型，证明了横向约束联合反演方法在MRS与TEM联合反演中的有效性.最后，通过实测数据验证了本文所述磁共振与瞬变电磁横向约束联合反演计算方法的实用性.