0 引言

按照场源类型，电磁法可以分为人工源电磁法和天然源电磁法。时频电磁法(TFEM)是东方地球物理公司在大功率建场测深方法基础上，结合时间域瞬变电磁法(TEM)和频率域可控源音频法(CSAMT)提出的一种应用于油气勘探的人工源电磁方法。该方法主要测量和研究长导线源激发产生的电磁场中的水平电场信号和垂直磁感应信号，并根据电阻率和极化率信息对目标圈闭进行含油气评价，经过十余年的发展和应用取得了良好效果^[1-6]。大地电磁法(MT)是一种利用天然交变电磁场研究地球电性结构的地球物理勘探方法，通过观测天然电磁场的电、磁场分量，处理后可获得地下电性结构信息。经过六十多年的理论完善和野外实践，MT已经广泛应用于油气和矿产资源勘探^[7-9]、地热资源勘探^[10]、环境和工程勘查与监测等地球物理勘探领域以及地壳和上地幔电性结构研究^[11-14]、地震监测和火山电性结构调查^[15-16]等地球深部构造研究领域。

在电磁勘探中，MT采用具有丰富频谱信息的天然源电磁信号，其勘探深度可达地下几千米到上百千米，能够提供区域尺度范围内的超深层盆地结晶基底及深大断裂构造的空间展布特征。但是，由于天然源电磁信号较微弱，特别是信号中存在死频带(0.1~10 Hz)，MT信号容易受外部环境噪声的影响，影响中深部电性层的分辨^[17-20]。TFEM法采用大功率人工场源克服天然电磁场的随机性和信号微弱的缺点，具有抗干扰能力强、对中浅层电性层具有较高分辨率的优点，但目前受限于场源的最低激发频率，对深层的电性信息反映有限。由于TFEM和MT都可反映地下构造的电性结构特征，且两者频段可很好地衔接，为开展TFEM和MT数据联合反演提供了先天有利条件。通过TFEM和MT数据联合反演处理，利用两者间的相互约束与信息互补，期望既改善浅部异常的分辨率，也能探测深部地电信息，有效降低电磁数据反演的多解性，提高深层地质目标体的反演精度和分辨能力。

自上世纪70年代以来，学者们开展了这两种电磁法的联合反演，初期集中在直流电(DC)法与MT的联合反演。Vozoff等^[21]在1975年率先发表了DC与MT数据的一维联合反演研究成果；之后，Sasaki^[22]提出了非线性最小二乘法的二维MT与偶极—偶极装置的DC法联合反演；McMillan等^[23]利用边界约束方法协同反演开展了DC、CSAMT和航空电磁(AEM)数据集的联合反演，并应用于金矿实测数据，取得良好的效果；李曼等^[24]基于高斯—牛顿法开展了非结构化三角网格条件下的二维直流电阻率与AMT自适应渐进正则化联合反演，结果表明联合反演可降低单一数据反演结果的不确定性。Meju^[25]开展了TEM与MT数据的联合反演，利用TEM响应不受浅层不均匀体影响的特点校正并改善了MT的反演效果。

随着可控源电磁法尤其是海洋可控源电磁法(MCSEM)的发展，大多数电磁法数据的联合反演研究都集中在MCSEM与海洋MT。Mackie等^[26]基于非线性共轭梯度方法开展了MCSEM与海洋MT数据联合反演，并用简单模型进行了测试，结果表明联合反演能够提高电性层的分辨率；Abubakar等^[27]采用乘法目标泛函方法研究了MCSEM与海洋MT的联合反演；Gribenko等^[28]利用积分方程方法进行了MCSEM与海洋MT数据的联合反演，基于复杂模型的测试效果较好；Wiik等^[29]评估了基于积分方程的MCSEM与海洋MT数据的对比源联合反演方案(Joint contrast source inversion, CSI)，认为并非所有情况下都需要进行联合反演，但联合反演有助于提高盐底辟成像效果；赵宁等^[30]采用奥克姆反演算法对MCSEM和海洋MT数据进行联合反演，测试结果表明联合频率更低的海洋MT进行联合反演能探测到深部基底电性结构，而不影响浅部地层的反演效果；汪茂^[31]基于MT与CSAMT方法探测深度不同的特点，通过这两种方法的雅可比矩阵相互堆叠的方式实现了二维数据空间的MT与CSAMT数据的联合反演；Miller等^[32]提出了基于构造引导的三维MCSEM和MT数据的联合反演，提高了反演结果的横向分辨率和构造的保真度；何展翔等^[4]采用人工鱼群算法实现TFEM法水平电场与垂直磁场的联合约束反演，模型测试和实际剖面反演结果表明该方法能够很好地降低反演的多解性，确保勘探精度和效果。

文献调研表明，截至目前鲜有陆上TFEM与MT数据联合反演相关研究成果发表，为此，本文基于奥克姆反演方法提出TFEM与MT数据的联合反演。首先，推导了一维模型TFEM响应正演公式，并将正演结果与三维积分方程正演结果进行对比，验证公式的正确性；然后，采用奥克姆反演算法，研究了一维TFEM与MT的联合反演问题，利用高/低阻模型对联合反演程序进行试算；最后，对实测数据进行处理。

1 TFEM正演模拟

MT一维正演原理参见文献[18, 33]，本文不再赘述，下面主要介绍TFEM的一维正演。

对于图 1所示包含N层地层的一维地电模型，地层的电导率和厚度分别为σ_i、h_i，i=1，…，N，图中IdL表示水平电偶极子，I表示电流，θ表示夹角∠xOr，r表示极距。取电磁场时谐因子为e^－iωt，其中ω为角频率。假设磁导率μ恒定不变，在TFEM研究频率范围(10^－2~10⁴ Hz)内位移电流可以忽略，麦克斯韦方程可写为

$ \left\{\begin{array}{l} \nabla \times \boldsymbol{E}=\mathrm{i} \omega \mu \boldsymbol{H} \\ \nabla \times \boldsymbol{H}=\sigma \boldsymbol{E}+\boldsymbol{J}_s \end{array}\right. $

(1)

式中：E、H分别表示电场和磁场；σ表示电导率；J_s表示源电流密度^[34-35]。

引用磁类型的矢量位A，磁场H可写为

$ \boldsymbol{H}=\frac{1}{\mu} \nabla \times \boldsymbol{A} $

(2)

基于式(1)和式(2)，电场E可写做

$ \boldsymbol{E}=\mathrm{i} \omega \boldsymbol{A}+\frac{1}{\mu \sigma} \nabla(\nabla \cdot \boldsymbol{A}) $

(3)

磁矢量势A可以基于汉克尔变换方程求得

$ \boldsymbol{A}(r)=\frac{1}{2 \pi} \int_0^{\infty} \hat{\boldsymbol{A}}(\lambda, z) J_0(\lambda r) \lambda \mathrm{d} \lambda $

(4)

式中：J₀(λr)是第一类零阶贝塞尔函数；λ是积分变量，为r方向的波矢量分量；$ \mathit{\boldsymbol{\widehat A}}(\lambda , z)$为傅里叶变换磁矢量势，其中z表示垂向距离。关于如何求取A以及电磁场分量，何继善^[17]、汤井田等^[20]、Key ^[34]、薛国强等^[35]给出了完整的推导公式，在此不再赘述。由于TFEM勘探野外测量的是与电场分量E_x和磁场分量H_z相关的参数^[1-3]，因而这里只需给出层状模型这两个分量的计算公式^{[17, 35]}

$ \left\{\begin{array}{l} E_x=\frac{\mathrm{i} \alpha \mu_0 I \mathrm{~d} L}{4 \pi}\left[\int_0^{\infty}\left(\frac{1}{K_0+K_1 G_1}-\frac{1}{\xi_0+\xi_1 G_2}\right) J_0(\lambda r) \lambda \mathrm{d} \lambda+\cos 2 \theta \int_0^{\infty}\left(\frac{1}{K_0+K_1 G_1}+\frac{1}{\xi_0+\xi_1 G_2}\right) J_2(\lambda r) \lambda \mathrm{d} \lambda\right] \\ H_z=\frac{I \mathrm{~d} L}{2 \pi} \sin \theta \int_0^{\infty} \frac{1}{K_0+K_1 G_1} J_1(\lambda r) \lambda^2 \mathrm{~d} \lambda \end{array}\right. $

(5)

式中：J₁(λr)、J₂(λr)分别为第一类1阶和2阶贝塞尔函数；$ K_i=\sqrt{\lambda^2-\kappa_i^2}, \xi_i=\frac{\kappa_i^2}{K_i}, \kappa_i$表示第i层的波数；G₁、G₂均为地层因子。

TFEM的激发源通常是长接地导线源，对电磁场沿导线方向积分即可求得TFEM的一维正演响应。为了验证一维TFEM正演算法的正确性，将其一维正演结果与三维积分方程正演结果进行对比^[36]。

设计一个5层地电模型，参数如表 1所示。发射源AB长度为7.5 km，收发距为8 km，电流强度为1 A，激发频率范围为10~1/30 Hz，共31个频点。TFEM一维正演频率按照实测数据提取部分频点；三维积分方程正演的最高频和最低频与一维正演频率一致，并采用对数等间隔取其间的29个频点。

图 2为一维TFEM模拟与三维积分方程正演的E_x和H_z实部和虚部分量对比。可以看出，一维与三维的正演结果几乎重合，验证了一维TFEM正演公式的正确性，为后续反演提供了可靠的基础。

2 TFEM与MT数据联合反演算法

根据经典的奥克姆反演算法^[37]，TFEM与MT数据联合反演问题的目标泛函可以表示为

$ U = ||\partial \mathit{\boldsymbol{m|}}{\mathit{\boldsymbol{|}}^2} + {\alpha ^{ - 1}}\left\{ {||\mathit{\boldsymbol{W}}[\mathit{\boldsymbol{d}} - F(\mathit{\boldsymbol{m}})]|{|^2} - X_*^2} \right\} $

(6)

式中：d表示观测数据向量；m表示模型向量；W表示数据加权矩阵；F表示正演算子；α为拉格朗日乘子，用以平衡粗糙度和数据拟合误差，当α取较大值时，反演以搜索光滑模型为主，反之则以搜索最小拟合误差为主；X_*²为数据拟合误差期望值。等式右边第一项为模型粗糙度因子，第二项为数据拟合误差。因为TFEM和MT反演的模型参数都是电阻率，所以这里的模型参数m为电阻率向量。假设有N_p个模型参数、M组数据d, 则

$ \|\partial \boldsymbol{m}\|^2=\sum\limits_{i=2}^{N_{\mathrm{p}}}\left(m_i-m_{i-1}\right)^2 $

(7)

$ \boldsymbol{W}=\left[\begin{array}{c} \boldsymbol{W}_{\mathrm{MT}} \\ \boldsymbol{W}_{\mathrm{TFEM}} \end{array}\right] $

(8)

式中：m_i表示第i个模型参数的值；W_MT和W_TFEM分别表示MT数据和TFEM数据的加权矩阵。数据加权矩阵W一般以每个观测数据标准差的导数作为该矩阵的对角元素，对每个观测数据在反演过程中起到加权作用^{[30, 34]}。由于TFEM与MT数据的量级差别比较大，在联合反演时需采用合适的加权因子去平衡TFEM和MT的数据权重。本文在上述加权方法基础上进行改进，提出数据加权矩阵

$ \left\{\begin{array}{l} \boldsymbol{W}=\left[\begin{array}{c} \boldsymbol{W}_{\mathrm{MTnew}} \\ \boldsymbol{W}_{\mathrm{TFEMnew}} \end{array}\right] \\ \boldsymbol{W}_{\mathrm{MTnew}}=\beta \sqrt{\operatorname{mean}\left(\phi_0^{\mathrm{MT}}\right)} \boldsymbol{W}_{\mathrm{MT}} \\ \boldsymbol{W}_{\mathrm{TFEMnew}}=(1-\beta) \sqrt{\operatorname{mean}\left(\phi_0^{\mathrm{TFEM}}\right)} \boldsymbol{W}_{\mathrm{TFEM}} \end{array}\right. $

(9)

式中：$ \phi_0^{\mathrm{MT}}$和$ \phi_0^{\text {TFEM }}$分别表示初始模型的MT和TFEM正演数据与相应观测数据拟合差；mean(·)表示求平均值；β为平衡因子，取值范围为0~1。利用式(9)求取满足泛函极值的模型参数即为反演结果，即

$ \frac{\partial U}{\partial \boldsymbol{m}}=\alpha^{-1}(\boldsymbol{W} \boldsymbol{J})^{\mathrm{T}} \boldsymbol{W}(\boldsymbol{J} \boldsymbol{m}-\boldsymbol{d})+\boldsymbol{R}^{\mathrm{T}} \boldsymbol{R} \boldsymbol{m}=0 $

(10)

式中：R为N_p×N_p阶粗糙度矩阵，形式如下

$ \boldsymbol{R}=\left[\begin{array}{ccccc} 0 & & & & \\ -1 & 1 & & & \\ & -1 & 1 & & \\ & & \ddots & \ddots & \\ & & & -1 & 1 \end{array}\right] $

(11)

J为M×N_p阶雅克比矩阵，其矩阵元素的表达式为

$ J_{i j}=\frac{\partial F_i(\boldsymbol{m})}{\partial m_j} $

(12)

计算雅克比矩阵是反演中重要的一步。由于采用TFEM与MT数据联合的方式，偏导数矩阵由两部分组成：一部分为TFEM的偏导数矩阵，采用电磁场对每层电导率的汉克尔变换公式解析求导获得^[38]；另外一部分为MT的偏导数矩阵，可以很容易地用差分方法求取。例如，第k+1次迭代的参数模型表达式为

$ \begin{aligned} \boldsymbol{m}_{k+1}=[ & \left.\alpha \boldsymbol{R}^{\mathrm{T}} \boldsymbol{R}+\left(\boldsymbol{W} \boldsymbol{J}_k\right)^{\mathrm{T}} \boldsymbol{W} \boldsymbol{J}_k\right]^{-1}\left(\boldsymbol{W} \boldsymbol{J}_k\right)^{\mathrm{T}} \cdot \\ & \boldsymbol{W}\left[\boldsymbol{d}_k-F(\boldsymbol{m})+\boldsymbol{J}_k \boldsymbol{m}_k\right] \end{aligned} $

(13)

基于上式，以上一次的计算结果作为此次计算的初始模型。α在第1次迭代时给出初始值，之后每次迭代都按照线性搜索方式计算得出。如果这种迭代算法最后收敛，就得到了泛函公式(6)的解，并且与初始模型不相关，最后的结果即是满足数据拟合误差前提下的光滑模型。反演数据拟合差为

$ \mathrm{RMS}=\sqrt{\frac{1}{N} \sum\limits_{i=1}^{N_{\mathrm{p}}} \frac{d_i-F_i(\boldsymbol{m})}{\delta_i{ }^2}} $

(14)

式中δ_i为第i个参数观测数据的标准误差。

3 模型算例

为验证联合反演算法的有效性，设计一层状模型进行测试。模型电阻率—深度曲线见图 3，从浅至深模型包括5个高、低阻互层，模型背景电阻率为100 Ω·m，具体参数见表 1。TFEM长导线发射源长7.5 km，收发距为8 km，发射31个频点，频率范围为9.6~0.032 Hz；模拟计算两个分量E_x和H_z；MT电磁场模拟38个频点，频率为320~0.001 Hz。

分别对TFEM、MT以及TFEM+MT三种数据进行反演，结果见图 3，可以看出：TFEM数据反演结果对浅部高阻层的恢复比MT数据更好，而MT数据反演结果对深部电阻率的恢复效果比TFEM数据反演结果更好；TFEM+MT数据的联合反演结果结合了两者的优点，对浅部高阻和深部基底的电阻率反演结果都优于TFEM数据或者MT数据单独反演的结果。图 4为图 3所对应的不同电磁数据反演迭代拟合差曲线，可见经50次迭代后，三种数据的反演拟合差均低于1.0。

一维模型算例结果表明，TFEM数据与MT数据的联合反演能够在一定程度上提高深、浅部电性层的分辨率。

4 实测数据反演

实测数据来源于中国东部M盆地。TFEM测线长度为15 km，点距为200 m，共76个测点，采集了E_x和H_z两个分量，数据采集频率范围为31.250~0.025 Hz，共36个频点，数据处理频率为9.600~0.032 Hz，共31个频点；MT数据与TFEM数据同点位采集，频率范围为320~0.0011 Hz，共38个频点。两种电磁数据质量较好，一级品率均大于85%。

图 5为不同电磁数据拟二维反演电阻率剖面。TFEM和MT数据都加入了5%误差。图 5a为TFEM数据单独反演的电阻率剖面，图 5b为MT数据单独反演电阻率剖面，图 5c为TFEM+MT数据联合反演电阻率剖面。对比图 5a与图 5b可以看出，TFEM和MT两种电磁数据单独反演的地电结构在深度小于6 km范围内基本一致，地层电阻率呈现低—高—低—高的四大套电阻层，1 km深处的低阻层电性特征一致性较好，TFEM数据单独反演剖面对深度400 m左右的高阻层反映较连续，MT数据单独反演剖面对该套高阻层的揭示不连续；深度大于6 km时，MT数据单独反演剖面揭示了一套低阻层，而TFEM数据反演剖面因为低频频率不够低，对深部地层的信息揭示不足；在测线水平方向13~15 km、埋深4 km区域，TFEM数据反演剖面揭示了一套较厚的低阻层，该低阻层在MT数据反演剖面上电阻率等值线出现轻微的扭曲。这种不吻合现象有可能是地层电阻率各向异性导致。TFEM与MT数据联合反演既保留了浅部连续高阻层，也保留了埋深大于8 km的深部低阻层。因此，TFEM与MT数据联合反演综合了两种电磁数据单独反演的优势，能够提供更丰富、更全面的地层电性信息。

5 结论

基于麦克斯韦方程，给出了一维TFEM电磁场分量正演的计算公式，通过对层状模型进行模拟正演，并与三维积分方程方法计算结果进行对比，验证了一维TFEM正演算法的正确性。

采用新的加权方法，基于奥克姆反演算法将TFEM与MT数据进行了联合反演。模型和实测数据试算结果表明，这两种电磁数据联合反演既可以提高浅部电性异常的分辨率，也能优化深部电性层的探测精度，提高电磁勘探解决地质问题的能力。

本文研究针对的是一维情况下人工源与天然源电磁数据的联合反演，实际地层电阻率更多情况下呈现二维、三维特性，甚至是各向异性^[39-41]，因此有必要开展二维、三维甚至是各向异性电阻率模型的联合反演研究，这也是我们下一步的研究方向。