0 引言

近年来，非地震勘探在川西地区火成岩、深层结构和断裂研究方面发挥了一定的指导作用，尤其是时频电磁法在寻找火山通道、探索火山岩形成机制方面效果良好^[1-3]。

时频电磁法^[4-8]是一种可控源电磁测深方法，该方法是基于油气藏或矿藏的激发极化特性以及其与围岩之间的电性差异，采用大功率、高精度人工源多频激发，将时间域电磁方法与频率域电磁方法有机结合，利用激发极化特性和电性综合异常检测和识别油气藏。相较于天然电磁法，可控源电磁法采用人工场源激发，即便在有工业干扰的地区，也能获得较高信噪比的数据。

时频电磁法基于原始资料分析和弱异常信息提取，通过精细反演，可寻找火山通道，探索火山岩形成机制；通过精细建模、针对火山岩约束反演，可刻画火山岩电阻率的平面变化，并结合不同岩性火山岩的电阻率变化特征，优选有利火山岩储集相带，有效区分火山碎屑岩与玄武岩，为地震勘探对火山岩岩性的识别提供了依据。

时频电磁法所采用的人工场源会导致非平面波效应、阴影效应以及场源附加效应，增加数据反演的多解性，使反演成像的质量偏低、甚至造成假异常，致使地质解释出现错误^[9-10]。本文利用时频电磁法多分量联合反演解决上述问题。

国内外学者做了大量关于可控源电磁勘探的多分量正、反演研究工作。陈明生等^[11]深入研究了二维模型条件下可控线电流源的瞬变电磁响应；王若等^[12]开展了二维模型中CSAMT的正演计算；陈小斌等^[13]利用有限元直接迭代算法正演计算了二维电性模型频率域有限长线源响应；Li等^[14]采用三角形剖分、有限元离散、背景场与散射场分离等方法消除场源处的奇异性，计算了海水中有限长线源产生的电磁响应。

可控源电磁数据的反演研究也取得了一定的成果。Key^[15]在Constable等^[16]的大地电磁法Occam反演算法的基础上，提出了适用于海洋可控源电磁数据的多分量反演算法；师学明等^[17]将最优化非线性反演方法和模拟退火方法引入电磁勘探数据多分量联合反演，取得了良好的效果；王若等^[18-19]实现了CSAMT全资料反演；何梅兴等^[20]基于Occam方法实现了CSAMT一维全资料反演；邱卫忠等^[21]详细介绍了电流源瞬变电磁法利用不同电磁场分量实现地下目标体的精确探测；王志刚等^[22]利用模拟退火法，实现了时频电磁法的多分量反演，并在实际数据处理中取得了好的效果；何展翔等^[23]基于人工鱼群算法实现了TFEM的水平电场和垂直磁场的联合约束反演，有效降低了反演结果的非唯一性，提高了反演精度。

实践中，如何平衡不同分量之间(如电道与磁道振幅)数量级的差异，防止出现绝对值较大或较小的数据在反演过程中占据主导地位、而其他分量起次要作用或不起作用的情况，是多分量联合反演成败的关键。为此，本文提出了在反演过程中增加动态平衡系数的方法，校正迭代搜索方向和模型更新步长，并对理论模型和实际数据进行了测试。

1 时频电磁二维正演理论

时频电磁法正演控制表达式为

$ \nabla \times \boldsymbol{E}=-\hat{z} \boldsymbol{H} $

(1)

$ \nabla \times \boldsymbol{H}=-\hat{y} \boldsymbol{E}+\boldsymbol{J}_{\mathrm{s}} $

(2)

式中：E表示电场；H表示磁场；J_s表示源电流密度；$\hat{z}$=iωμ₀表示阻抗率，其中ω表示角频率，μ₀表示真空中的磁导率；ŷ=iεω+σ表示导纳率，其中σ表示电导率，ε表示相对介电常数。

假设地电模型沿x方向延伸，该模型可简化为一个关于y和z的函数

$ \sigma=\sigma(y, z) $

(3)

从式(2)可以看出右边存在一个源项J_s，由于在源处存在奇异性，很难对其进行数值模拟。为了消除源项引起的奇异性，可借鉴Li等^[14]提出的背景场与散射场分离方法。

2 控制方程有限差分离散

有限差分法是一种经典的数值模拟计算方法，其基本原理是用差商替代微商，将待求解的连续微分方程变换为离散的差分方程，并通过求解差分方程得到源微分方程的近似解。其优点是方法简便、易于实现。

模型的网格剖分采用矩形剖分，如图 1所示。设横向剖分节点数为N+1，即横向网格数为N；纵向剖分为M+1个节点，即纵向网格数为M。总计N×M个网格。将网格按照从上至下、从左至右的方式顺序自然编号。设待求解的电、磁场分量分别为U=E_xⁱ, i=1, 2, …, N×M, V={H_xⁱ, i=1, 2, …, N×M}，计算点位于网格的中心位置。采用九点差分格式进行离散，并对局部节点进行自然顺序编号(A₁~A₉)，如图 2所示。

将局部离散方程组装并加上狄里克莱边界条件，形成一个大型稀疏线性方程

$ \bf{K} \cdot \bf{X}=\boldsymbol{S} $

(4)

式中：K表示与频率和电阻率有关的复对称矩阵；S是与频率和电阻率有关的向量；X表示电场或者磁场分量。求解该方程组，即可得到电、磁场响应。

3 多分量联合反演原理

本文采用正则化反演原理实现时频电磁多分量联合反演。首先，构造一个反演目标泛函数(简称目标泛函)。目标泛函一般由两部分组成，其中一个是反演过程稳定函数，简称稳定函数。稳定函数有多种实现方式，本文采用Occam算法所采用的模型导数作为稳定函数^[15-16]，其一般形为

$ \|\partial \mathit{\boldsymbol{m}}\|^{2}=\min $

(5)

式中：∂为正则化反演模型稳定器，一般为加权求导矩阵算子形式；m为反演模型参数；min表示取最小值。另一个是惩罚函数，用于度量反演模型的响应与观测数据之间的偏差，其一般形为

$ \|\boldsymbol{W} \boldsymbol{d}-\boldsymbol{W} f(\boldsymbol{m})\|^{2}=\min $

(6)

式中：d为参与反演的观测数据；f为正演响应计算函数；W为观测数据的权重。利用拉格朗日乘子理论可求得反演目标泛函

$ \begin{aligned} T(\boldsymbol{m}) &=\|\partial \boldsymbol{m}\|^{2}+\xi^{-1}\left[\|\boldsymbol{W} \boldsymbol{d}-\boldsymbol{W} f(\boldsymbol{m})\|^{2}-\chi^{2}\right] \\ &=\min \end{aligned} $

(7)

式中：ξ为拉格朗日乘子，用于反演过程中平衡惩罚函数与稳定函数；χ为期望拟合差。

时频电磁的响应函数一般为最小相位系统，其振幅与相位之间存在一定的联系，即振幅包含的信息等价于相位包含的信息，故时频电磁反演一般采用振幅数据。

时频电磁的电场分量E_x与磁场分量B_z的振幅量级差异较大，一般相差4个量级。如果直接将两个分量的振幅数据代入反演目标函数，可能会导致某个分量的数据占优，而另一分量的数据在反演过程中不能对反演模型的修改发挥作用，这种现象称为数据间的相互“淹没”。

为了克服不同分量数据的量级差异，笔者认为可以有两种方法解决这一问题：一是将不同分量的振幅转换为视电阻率，一般采用可控源电磁全区视电阻率，这点将另文论述，本文不做详细讨论；另一种是将两类分量数据(包括振幅和相位)分别作为两个数据集合代入目标函数

$ \begin{aligned} T(\boldsymbol{m})=&\|\partial \boldsymbol{m}\|^{2}+\xi^{-1}\left[\left\|\boldsymbol{W} \boldsymbol{d}_{E_{x}}-\boldsymbol{W} f(\boldsymbol{m})\right\|^{2}+\right.\\ &\left.\beta\left\|\boldsymbol{W} \boldsymbol{d}_{B_{z}}-\boldsymbol{W} f(\boldsymbol{m})\right\|^{2}-\chi^{2}\right]=\min \end{aligned} $

(8)

式中β为不同分量数据量级差异的平衡因子，后文简称平衡因子。

β可以采取固定赋值的方法，即处理人员根据经验设置，但是实际数据E_x与B_z的振幅量级差异是变化的，需要处理人员不断测试平衡因子的有效取值范围，工作量大，且人为因素多，故本文不采用这种方法。正则化反演是一种迭代寻优的数学过程，本文采用动态统计估算平衡算法，即在每一次迭代过程中，利用当前迭代反演模型的响应，利用下式

$ \beta=\frac{\operatorname{median}\left[\left\|\boldsymbol{W} \boldsymbol{d}_{E_{x}}-\boldsymbol{W} f(\boldsymbol{m})\right\|^{2}\right]}{\operatorname{median}\left[\left\|\boldsymbol{W} \boldsymbol{d}_{B_{z}}-\boldsymbol{W} f(\boldsymbol{m})\right\|^{2}\right]} $

(9)

估算平衡因子β。式中median是中值函数。

4 算法测试 4.1 合成数据测试

为了验证算法的正确性，设计包含两个异常体的二维理论模型(图 3a)，正演得到合成测试数据。模型背景为均匀介质，电阻率为10Ω·m。模型中有两个目标体，左边的低阻目标体电阻率为1Ω·m，右边的高阻目标体电阻率为100Ω·m。反演初始模型为10Ω·m的均匀半空间。

首先，分别进行单分量(E_x和B_z)数据反演，结果见图 3b和图 3c；然后，对E_x和B_z进行联合反演，结果见图 3d。由图 3可见：①E_x数据单独反演结果(图 3b)对低阻目标的形态恢复较好，对高阻目标成像效果不佳，在高阻目标下方形成了明显的假异常；②B_z数据单独反演结果(图 3c)对高阻目标的形态恢复较好，但低阻目标下方形成了明显的假异常；③E_x和B_z联合反演结果(图 3d)较好地克服了E_x或B_z单独反演结果中的假异常现象。该理论模型反演结果证明了本文提出的联合反演算法的正确性和有效性。

4.2 反演收敛的一致性

为了验证算法收敛的一致性，抽取上节中的一道磁场振幅数据，并加入5%的高斯白噪声，如图 4所示。分别从四个不同的初始模型开始反演，以验证本文算法的收敛一致性。这四个初始模型都是均匀半空间，其电阻率分别为1、10、30、50Ω·m，E_x和B_z联合反演结果见图 5。可以发现，本文算法对不同的反演初始模型都能达到一致收敛的效果，对加噪数据能有效恢复出理论模型的形态。

4.3 反演收敛的稳定性

为了验证算法收敛的稳定性，分别在上节磁道数据中加入2%、5%、10%和20%的高斯白噪声进行反演，以研究噪声水平的变化对反演结果的影响，结果见图 6。由图可见, 随着噪声水平的增大，反演模型与真实模型的偏差越来越大, 但是反演电阻率剖面仍能清晰地反映异常体的位置及其电阻率的高、低特征，即使噪声水平超过10%时，剖面上异常体的位置也能准确识别，这说明了本文算法的收敛稳定性。

5 应用效果

近年来，川西地区YT1井二叠系获高产气流，开辟了火山岩勘探新领域。YT1井储层主要包括火山碎屑岩和玄武岩两大类。其中火山碎屑熔岩发育厚层孔隙型储层，平均孔隙度超过10%，厚度超过100m，是该区油气主力储层。火山碎屑熔岩储集空间以脱玻化溶蚀微孔、残余气孔为主，裂缝发育程度较低。玄武岩岩性致密，孔隙不发育，储集空间主要为裂缝。储层在区内分布面积较大，达6000km²，勘探潜力巨大。地震反射总体表现为丘状杂乱反射，特征不明显，储层识别难度较大。该区储层与围岩的电性特征差异明显，这是利用电磁法在川西地区开展火山岩、深层结构和断裂研究的良好物性基础，为此部署了时频电磁勘探。

此次勘探的部分测点位于三维地震区内，对过地震三维区内的时频测线进行了多分量联合约束反演，选取其中一条典型剖面进行时频电磁多分量联合反演，结果见图 7。可以看出，该剖面纵向上可划分为电阻率高低变化的多套电性层，纵向上的电性特征与电阻率测井曲线规律一致，说明时频电磁多分量联合反演结果较真实地反映了研究区电性结构。剖面中部偏移距25~30km范围的深部高阻层是火山通道断裂发育的表现，这里的断裂不仅是火山喷发通道，而且沿断裂还有侵入岩发育，这点已得到YT1井区钻探结果(图 7中电测井曲线)的证实。

为了进一步验证时频电磁多分量联合反演的效果，对YS1井和YT1井(位置见图 7)附近的点进行E_x分量单独反演及E_x与B_z联合反演，并与这两口井的电测井曲线进行对比，结果见图 8。根据钻井信息，这两口井都钻遇了二叠系火山岩。可以看出，单分量反演和联合反演的电阻率与电阻率测井曲线的变化趋势吻合较好，但多分量联合反演比E_x分量单独反演结果的纵向分辨率更高，能够更好地揭示不同岩性火山岩引起的电阻率变化细节。

6 结论

本文基于正则化反演原理，开发了适用于时频电磁E_x、B_z分量的联合反演方法。理论模型合成数据测试结果表明，相对于单分量数据反演，E_x与B_z数据联合反演对高阻、低阻目标体的形态恢复更好，可解决单分量反演对特定电性特征的目标形态恢复不佳的问题。另外，模型理论数据加噪测试结果表明该算法稳定、有效。实测数据反演结果表明，时频电磁多分量联合反演剖面较好地反映了火山岩的分布及火山通道(断裂)的存在，与该区电测井数据较吻合。

理论模型和实测数据的计算结果表明，时频电磁法基于电阻率空间分布特征，结合不同岩性火山岩电阻率变化规律，可有效刻画火山岩的分布，为优选有利火山岩储集相带提供有效信息，是一种兼具经济性、可靠性的勘探手段。

另外，激发极化特征也是研究火山岩的重要物理基础，笔者后续将在本文研究的基础上，进一步开展基于Cole-Cole模型的激发极化正、反演研究。