0 引言

时频电磁法^[1-2]是一种人工场源电磁法，广泛应用于油气勘探。该方法基于可控源音频大地电磁测深法与长偏移距瞬变电磁法，实现了同源激励、同步采集时域和频域电磁信号，提升了工作效率，提高了勘探精度^[3]，在油气勘探实践中取得良好的应用效果^[4-6]。随着计算方法和计算机硬件的发展，地球物理数据的反演逐步向精细化方向发展，时频电磁数据的反演算法也需进一步优化，以提升反演精度和速度。

从理论上来说，时间域与频率域是等效的，时域信号与频域信号可通过傅里叶变换相互转换，但实际采集信号时二者存在较大差异。陈卫营等^[7]通过对比分析瞬变电磁法与可控源音频大地电磁测深法发现，相比频域电磁法，时域电磁法对薄电性层的探测能力更强；柳建新等^[8]通过理论模型正演发现，时域及频域电磁信息对高、低阻异常体的敏感程度不同，时域磁场对低阻异常体更敏感，而频域电场对高阻异常体敏感度更高；王若等^[9]对可控源瞬变电磁法各分量进行灵敏度分析，认为不同地电条件下，电场与磁场对目标地层的分辨率存在差异。因此，综合利用电、磁分量进行联合反演是优化反演效果、提高勘探精度的有效途径。

随着时频电磁勘探方法的推广和应用，时、频电磁数据的联合反演方法也得到快速发展。柳建新等^[8]利用广义逆方法实现了时域磁场与频域电场数据的联合反演；岳明鑫等^[10]基于正则化反演实现了浅海一维时频电磁联合反演；何展翔等^[11]基于改进的人工鱼群算法实现了时频电磁数据的分步约束反演；郭鸿喜等^[12]利用改进的OCCAM算法实现了时频电磁数据多分量的联合反演；曹杨等^[13]对反演采用模型约束，反演结果精度更高。以上文献研究分析表明，电磁数据联合反演结果优于单一数据类型的反演结果。

当前，时频电磁联合反演算法大多基于最大平滑准则获取反演问题的稳定解，但反演结果常出现边界模糊的现象，不能清晰地反映真实地电信息^[13]。为获得精细化地电模型，提高反演稳定性，本文在前人研究基础上将多种约束方案引入反演算法，提出一种带约束的一维自适应时频电磁联合反演算法。通过对经典层状模型、含噪页岩气地电模型及实测数据反演测试与分析表明，带约束的时频电磁联合反演可有效提高复杂模型电阻率反演结果的分辨率，对于提高时频电磁方法油气藏识别精度具有实际意义。

1 反演方法和理论 1.1 联合反演函数

本文采用频域水平电场分量(E_x)与时域垂直磁场时间导数(dH_z/dt)进行联合反演，使用模型约束的正则化反演方法^[14-15]构建联合反演目标函数

$ \sigma(\boldsymbol{m})=\varphi_{\mathrm{F}}+\lambda \varphi_{\mathrm{T}}+\beta \varphi_{\text {model }} $

(1)

式中：σ(m)为总目标函数，其中m代表模型参数向量，本文反演算法固定模型各层厚度，仅反演模型电阻率，同时对模型参数引入自然边界约束，因此m是由各层电阻率对数值组成的向量；φ_F、φ_T分别表示频域和时域数据拟合项；λ为数据自适应平衡因子；β为正则化因子；φ_model为模型约束项。式(1)中的数据拟合项与模型约束项的矩阵形式分别为

$ \begin{aligned} \varphi_{\mathrm{F}}+\lambda \varphi_{\mathrm{T}}= & \left\|\boldsymbol{W}_{\mathrm{d}}\left[\boldsymbol{Q}_{\mathrm{F}}-F_{\mathrm{F}}(\boldsymbol{m})\right]\right\|^2+ \\ & \lambda\left\|\boldsymbol{W}_{\mathrm{d}}\left[\boldsymbol{Q}_{\mathrm{T}}-F_{\mathrm{T}}(\boldsymbol{m})\right]\right\|^2 \end{aligned} $

(2)

$ \begin{aligned} \beta \varphi_{\text {model }}= & \beta \times\left[a_{\mathrm{s}}\left\|\boldsymbol{W}_{\mathrm{s}}\left(\boldsymbol{m}-\boldsymbol{m}_{\mathrm{ref}}\right)\right\|^2+\right. \\ & \left.a_z\left\|\boldsymbol{W}_z(\boldsymbol{m})\right\|^2\right] \end{aligned} $

(3)

式中：W_d为数据加权矩阵，是由数据协方差倒数组成的对角矩阵，其作用是基于数据协方差对观测数据进行归一化操作；W_s表示最小模型加权矩阵；W_z为模型粗糙度矩阵，定义为模型一阶差分算子，与模型划分层数及层厚度有关；Q_F、Q_T分别表示频域与时域的观测数据；F_F、F_T分别表示频域与时域的正演算子; a_s、a_z为模型约束比重；m_ref为先验模型。

构建模型约束目标函数时可同时考虑模型最小构造及平缓度^[14]。式(3) 中, $a_{\mathrm{s}}\left\|\boldsymbol{W}_{\mathrm{s}}\left(\boldsymbol{m}-\boldsymbol{m}_{\text {ref }}\right)\right\|^2$、$a_\mathrm{z}\left\|\boldsymbol{W}_\mathrm{z}(\boldsymbol{m})\right\|^2$分别对应最小模型约束与最平缓模型约束, 根据经验取$a_{\mathrm{s}}=0.01 、a_\mathrm{z}=0.99$, 因此光滑模型在反演中起主要作用; 最小模型加权矩阵$\boldsymbol{W}_{\mathrm{s}}$取单位矩阵。

反演采用高斯—牛顿最小化思想将非线性化问题转换为线性化问题，在线性化目标函数最小情况下，通过方程

$ \begin{aligned} & \left(\boldsymbol{J}_1^{\mathrm{T}} \boldsymbol{W}_{\mathrm{d}}^{\mathrm{T}} \boldsymbol{W}_{\mathrm{d}} \boldsymbol{J}_1+\lambda \boldsymbol{J}_2^{\mathrm{T}} \boldsymbol{W}_{\mathrm{d}}^{\mathrm{T}} \boldsymbol{W}_{\mathrm{d}} \boldsymbol{J}_2+\beta \boldsymbol{W}_{\mathrm{m}}^{\mathrm{T}} \boldsymbol{W}_{\mathrm{m}}\right) \delta \boldsymbol{m} \\ & \;\;\;\;\;\;\;\;\;\;\; =\boldsymbol{J}_1^{\mathrm{T}} \boldsymbol{W}_{\mathrm{d}}^{\mathrm{T}} \boldsymbol{W}_{\mathrm{d}}\left\{\boldsymbol{Q}_{\mathrm{F}}-F_{\mathrm{F}}\left[\boldsymbol{m}^{(n)}\right]\right\}+ \\ & \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\lambda \boldsymbol{J}_2^{\mathrm{T}} \boldsymbol{W}_{\mathrm{d}}^{\mathrm{T}} \boldsymbol{W}_{\mathrm{d}}\left\{\boldsymbol{Q}_{\mathrm{T}}-F_{\mathrm{T}}\left[\boldsymbol{m}^{(n)}\right]\right\}- \\ & \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\beta \boldsymbol{W}_{\mathrm{m}}^{\mathrm{T}} \boldsymbol{W}_{\mathrm{m}} \boldsymbol{m}^{(n)}+\beta a_{\mathrm{s}} \boldsymbol{W}_{\mathrm{s}}^{\mathrm{T}} \boldsymbol{W}_{\mathrm{s}}\left(\boldsymbol{m}_{\mathrm{ref}}\right) \\ & \end{aligned} $

(4)

可求得当前模型修正量$\delta \boldsymbol{m}$。式中: $\boldsymbol{J}_1 、\boldsymbol{J}_2$分别为频域与时域雅克比矩阵; $\boldsymbol{W}_{\mathrm{m}}$是模型约束矩阵, $\boldsymbol{W}_{\mathrm{m}}{ }^{\mathrm{T}} \boldsymbol{W}_{\mathrm{m}}$ $=a_{\mathrm{s}} \boldsymbol{W}_{\mathrm{s}}^{\mathrm{T}} \boldsymbol{W}_{\mathrm{s}}+a_\mathrm{z}\ \boldsymbol{W}_\mathrm{z}^{\mathrm{T}} \boldsymbol{W}_\mathrm{z}\ ; \boldsymbol{m}^{(n)}$代表第$n$次迭代时的模型向量。通过

$ \boldsymbol{m}^{(n+1)}=\boldsymbol{m}^{(n)}+d_n \delta \boldsymbol{m} $

(5)

对模型进行迭代、修改，逐步实现数据拟合，直至满足迭代终止条件。式中d_n为第n次迭代步长。

反演设置三种迭代终止阈值：反演迭代次数(20)，模型修改量(二范数1×10^－5)，拟合差(5×10^－5)。

1.2 数据加权方案

通过联合反演算法试算发现，时域及频域数据常因数据量不等或数据的量级差异过大，导致欠拟合现象。因此，本文通过对文献[8]中的平衡因子λ进行改进实现数据拟合项的均衡加权。改进的平衡因子为

$ \lambda_n=\frac{N_{\mathrm{F}}}{N_{\mathrm{T}}} \times C \times \frac{\max \left(\Delta \boldsymbol{Q}_{\mathrm{T}}\right)-\min \left(\Delta \boldsymbol{Q}_{\mathrm{T}}\right)}{\max \left(\Delta \boldsymbol{Q}_{\mathrm{F}}\right)-\min \left(\Delta \boldsymbol{Q}_{\mathrm{F}}\right)} $

(6)

式中：λ_n为第n次迭代时的平衡因子；N_F、N_T分别代表频域与时域的总数据量；C为放大因子，本文取2；ΔQ_F、ΔQ_T分别为归一化后频域与时域数据拟合差组成的数组；max(·)、min(·)分别表示取最大值和最小值。

为节约计算资源，对迭代过程进行自适应优化，其思路为：求取平衡因子初值及拟合差，若后续迭代中拟合差大于前一次迭代拟合差，则据式(6)更新平衡因子，否则沿用前次平衡因子，实现反演过程中数据平衡因子的自适应动态调整。

1.3 正则化因子的选取

反演数据与模型参数之间的巨大差异会导致反演结果不稳定，其中正则化因子会影响反演效率及收敛性，直接关系到反演成败，因此需谨慎选择。本研究选用陈小斌等^[16]提出的CMD自适应正则化方案，对时频电磁联合反演进行适配优化，增添缩放因子。改进后的正则化因子为

$ \beta_{n+1}=\frac{\varphi_{\mathrm{F}}\left[\boldsymbol{m}^{(n)}\right]+\lambda \varphi_{\mathrm{T}}\left[\boldsymbol{m}^{(n)}\right]}{\varphi_{\mathrm{F}}\left[\boldsymbol{m}^{(n)}\right]+\lambda \varphi_{\mathrm{T}}\left[\boldsymbol{m}^{(n)}\right]+\varphi_{\text {model }}\left[\boldsymbol{m}^{(n)}\right]} \times \omega $

(7)

式中：β_n+1为第n+1次迭代的正则化因子；ω为缩放因子，该参数采用阶梯式变化方案，在反演初期ω为放大因子，加强模型约束，在后期ω为缩小因子，侧重数据拟合，ω取值需考虑反演收敛性，若反演不收敛则采用较小的值。

1.4 参数搜索范围和步长约束

恰当的正则化因子可从数学角度确保迭代反演的稳定与收敛，但无法保证每次迭代结果均符合自然规律与实际情况。因此，在自然边界约束的基础上，基于可行域思想对模型电阻率参数搜索范围进行限制，具体步骤为：

(1) 根据研究区域的已知电阻率信息设定反演电阻率参数的上、下限；

(2) 记录当前模型各层电阻率原始状态，并求解模型修正量；

(3) 使用当前模型修正量对各层电阻率进行更新；

(4) 对更新后的各层电阻率依次进行范围判断，若某层电阻率超出预设范围，则将各层电阻率恢复至原始状态，再将模型修正量整体进行固定比例缩放(缩放因子取经验值0.85)，并记录缩放次数；

(5) 重复步骤(3)和步骤(4)，直至各层电阻率满足预设范围要求；

(6) 更新各层电阻率参数，同时根据缩放次数计算步长。

此方案无需修改雅克比矩阵计算式，且不影响目标函数整体下降方向^[17]，仅需通过压缩迭代步长实现反演参数搜索范围的约束，可在较大程度上保留模型修改量的同时确保反演收敛。但是，模型试算发现，经此约束后采用的迭代步长仍然较大，反演易快速陷入局部极小值，导致反演失败，因此需对约束方案再优化，即将可行域思想与WOLFE准则^[18]相结合，对迭代步长再次进行筛选。

WOLFE准则要求步长满足充分下降条件及曲率约束条件，充分下降条件下所筛选的步长能保证目标函数下降。但若步长偏小，则无法保证模型有效更新，因此补充曲率约束条件，在充分下降条件下保证步长存在，通过曲率约束确保所筛选步长有效。以参数搜索范围进行约束下求出的步长作为最大允许步长，代入线性搜索算法，可获得满足WOLFE准则的迭代步长，实现步长约束。

2 反演算例

设计经典地电模型及复杂层状模型进行算法试验，检验本文联合反演算法的准确性与可靠性，并通过与单一数据类型的约束反演结果进行对比，说明联合数据反演效果。以下反演算例中，初始模型均设定为100 Ω·m的均匀半空间。需注意的是，后文中的时频电磁联合反演拟合差是时域及频域数据拟合差之和，因此此参数仅用于判断反演算法收敛与否。

2.1 三层模型

选择典型的H型与K型地电模型进行反演试算。

H型模型(图 1a中黑色折线)从浅至深各层电阻率为100、10、100 Ω·m，对应的层厚度为300、200 m。反演模型划分为30层，每层厚度均为30 m，最后一层电阻率参数填充至最大反演深度(1 km)，电阻率搜索范围设定为1~400 Ω·m。

H型模型反演结果见图 1。由图 1a可见，单一数据反演结果及时频电磁联合数据反演较直观地反映出了低阻异常层及其深度，但后者更接近真实模型。由图 1b可见，三种反演拟合差均快速下降，反演快速收敛。真实模型、联合反演结果、频域数据反演结果的频域水平电场分量E_x拟合曲线和时域dH_z/dt拟合曲线分别见图 1c和图 1d，可见反演模型的响应曲线与实际模型均拟合较好。

建立图 2a中黑色折线所示K型地电模型，从浅至深各层电阻率为20、200、20 Ω·m，对应的厚度为300、200 m。反演模型划分为30层，每层厚度均为30 m，最后一层电阻率参数填充至最大反演深度(1 km)，电阻率搜索范围设定为1~400 Ω·m。

K型模型反演结果见图 2。由图 2a可见，三种反演结果基本能反映真实模型的基本特征，但因浅部低阻层的电流屏蔽作用，反演结果中高阻层与真实模型电阻率差异较大，相对而言，时频电磁联合反演结果更接近于真实模型。由图 2b可见，频域数据、时域数据和时频数据联合反演均能快速收敛。由图 2c和图 2d可见，反演模型在时域和频域的响应曲线与实际模型均拟合较好。

2.2 四层模型

设计相对复杂的四层HK及KH地电模型对算法进行验证。

建立图 3a黑色折线所示HK型地电模型。从浅至深各层电阻率为100、10、300、80 Ω·m，厚度为100、50、250 m。反演模型划分为30层，每层厚度均为30 m，最后一层电阻率参数填充至最大反演深度(1 km)，电阻率搜索范围设定为1~600 Ω·m。

HK型模型反演结果见图 3。由图可见，三种反演结果均能反映真实模型的基本特征，但频域数据反演结果中高阻异常层对应的电阻率偏高，对基底电阻率刻画不准确；时域数据反演得到的低阻异常与真实模型接近，高阻异常与真实模型偏差较大，且在异常边界处反演结果出现明显波动。对比单一数据类型与联合数据反演结果可知，时频电磁联合反演能更准确识别模型特征，反演得到的高阻异常电性特征更接近实际模型。

建立图 4a中黑色折线所示KH型地电模型，从浅至深各层电阻率为20、200、40、200 Ω·m，对应各层厚度为50、200、100 m。反演模型划分为30层，每层厚度均为30 m，最后一层电阻率参数填充至最大反演深度(1 km)，电阻率搜索范围为1~600 Ω·m。

KH型模型反演结果见图 4，可见时域数据、频域数据、时频数据联合反演结果对不同电阻率的地层识别效果有所差异。基于频域数据反演曲线可识别浅部高阻异常层，但其厚度大于真实值，且存在一定幅度的波动，对基底的电阻率反演结果偏差也较大。时域数据反演的电阻率异常较实际模型的真实数值偏小，对低阻层的反映不明显。时频电磁数据联合反演结果对KH型模型的拟合效果相对最佳，能较准确反映高、低阻异常层的深度，更接近真实模型。

2.3 实测页岩气模型

根据川南扬子地区下古生界龙马溪组气层页岩气勘探成果^[19-20]，结合测井曲线建立图 5a中包含7个地层的地电模型(黑色折线)。从浅至深各层电阻率分别90、60、500、150、100、15、210 Ω·m，对应各层厚度为400、600、300、200、500、300 m，其中第6层为目标层。对该模型的正演电磁场添加3%随机噪声。反演模型划分为30层，每层厚度均为100 m，反演最大深度设为3 km，电阻率搜索范围为1~700 Ω·m。

页岩气模型及反演、拟合结果见图 5。由图可见，时域数据、频域数据、时频数据联合反演过程均稳定收敛，实现数据拟合。由图 5a可见：频域数据单独反演由于受浅部高阻层屏蔽作用的影响，下伏低阻目标层无异常显示；时域数据反演的电阻率异常幅值太小，无法反映低阻目标层的电性特征，对上覆的高阻层电性特征反映也不明显；相对于单一的时域数据和频域数据反演结果，时频电磁联合反演结果最接近于真实模型，可有效识别低阻目标层及上覆高阻层。

2.4 实测数据反演

为验证时频电磁数据联合反演的效果，选用中国N油田一单点实测数据进行反演测试。反演模型划分为34层，每层厚度均为100 m，反演最大深度为3.4 km，电阻率搜索范围为1~500 Ω·m。

图 6a为该单点数据分别基于时域数据、频域数据、时频电磁数据反演的电阻率曲线。拟合误差曲线见图 6b。图 7所示为测点附近的电阻率测井曲线(油水同层顶面深度约2900 m)。由图 6可知，随着迭代次数增加，三种反演方案的拟合差逐步下降，数据拟合效果较好。对比图 6与图 7，时频电磁联合反演的结果最接近于相邻电阻率测井曲线，对油水同层的低阻异常还原效果最好，证明时频电磁数据联合反演对油藏的分辨能力较强。

3 结论

为提高时频电磁反演的稳定性，提高反演结果的分辨率，本文将多种约束方案引入反演算法。基于正则化反演方法，添加数据约束、模型约束及反演步长约束，提出一种一维自适应时频电磁联合反演方法。该方法对时域和频域数据权重因子及正则化因子采用自适应调整方案，同时将模型参数搜索范围约束与WOLFE准则相结合，在限制参数搜索范围的同时约束迭代步长，保证反演过程的稳定性。通过经典层状模型、实际页岩气模型及实测电磁数据的反演，证明了本文提出的联合反演方案效果明显优于单一数据类型反演，可有效提高复杂模型电阻率反演结果的分辨率，对于提高时频电磁法油气藏识别精度具有实际意义。