地球物理学进展

核函数阶数	对数离散	线性离散
100	7.7547×10¹⁹	1.4365×10¹⁹
200	7.7670×10²⁰	3.5176×10¹⁹
500	4.4176×10²¹	2.2730×10²⁰

核函数阶数	对数离散	线性离散
100	7.7547×10¹⁹	1.4365×10¹⁹
200	7.7670×10²⁰	3.5176×10¹⁹
500	4.4176×10²¹	2.2730×10²⁰

基于波场特征约束的瞬变电磁波场反变换优化算法研究

PDF下载

曹佳慧 ¹^,² , 戚志鹏 ¹^,²^,^* , 李文翰 ¹^,² , 李貅 ¹^,²

地球物理学进展 | 应用地球物理学（油气、金属矿、海洋、环境、工程、仪器等） 2026,41(1): 217-227

收起

地球物理学进展 |应用地球物理学（油气、金属矿、海洋、环境、工程、仪器等） 2026 , 41 (1) : 217 -227

基于波场特征约束的瞬变电磁波场反变换优化算法研究

全屏

曹佳慧¹^,², 戚志鹏¹^,²^,^*, 李文翰¹^,², 李貅¹^,²

作者信息

¹长安大学地质工程与测绘学院, 西安 710054

²长安大学地球物理场多参数综合模拟实验室(中国地球物理学会重点实验室), 西安 710054

通讯作者:

戚志鹏, 男, 1982年生, 博士生导师, 教授, 主要从事瞬变电磁场的理论与应用方面的研究. E-mail: qzhipeng@126.com

作者简介:

曹佳慧, 女, 1999年生, 硕士研究生, 主要从事瞬变电磁场的理论与应用方面的研究. E-mail: Caojiahui×d.edu.cn

Research on inverse transform optimization algorithm of transient electromagnetic wave field based on wave field feature constraint

JiaHui CAO¹^,², ZhiPeng QI¹^,²^,^*, WenHan LI¹^,², Xiu LI¹^,²

Affiliations

¹College of Geological Engineering and Geomatics, Chang'an University, Xi'an 710054, China

²Integrated Geophysical Simulation Lab of Chang'an University (Key Laboratory of Chinese Geophysical Society), Xi'an 710054, China

出版时间: 2026-02-20 doi: 10.6038/pg2026II0594

文章导航

摘要

收起

为了利用波场特征进行瞬变电磁数据解释, 提高电磁探测的解释精度, 突破传统扩散场解释能力的局限性, 本文以构建融合扩散场数据特征和波场频率特征的目标函数为核心, 提出一种基于扩散场与波动场之间的积分变换关系的多约束波场变换算法.首先, 根据波场变换理论构建目标函数框架, 在目标函数中引入考虑滤波过程的约束项和考虑波场特征的约束项; 然后, 对目标函数进行优化求解, 在计算过程中为了避免计算海森矩阵, 采用拟牛顿算法对目标函数进行优化, 获得了既满足数据约束又符合波场特征的虚拟波场.最后, 利用文中算法对层状模型和复杂模型的瞬变电磁正演数据进行波场变换, 验证算法的有效性与稳定性.数值试验结果表明, 利用文中算法得到的波场保留了波场的传播特性, 而且与模型波场正演结果吻合, 充分说明了文中反变换算法的有效性.该研究为瞬变电磁拟波动方程解释提供了理论基础, 在深部矿产精细探测与复杂地质构造识别中具有应用价值.

关键词

瞬变电磁 / 波场变换 / 波场特征约束 / 拟牛顿优化 / 滤波约束

Abstract

收起

To utilize wave field characteristics for transient electromagnetic data interpretation, enhance the interpretation accuracy of electromagnetic detection, and overcome the limitations of traditional diffusion field interpretation capabilities. This paper proposes a multi-constrained wavefield transformation algorithm based on integral transform relationships between diffusion and wavefields, with its core being the construction of an objective function that integrates diffusion field data characteristics and wavefield frequency features. First, based on the wavefield transformation theory, this study constructs an objective function framework. Within this framework, filtering process constraints and wavefield characteristic constraints are introduced. The filtering constraints mitigate the loss of wavefield resolution caused by traditional smoothing constraints, while the wavefield characteristic constraints further narrow the solution space based on wavefield propagation features. These two constraint terms are integrated into a unified objective function to enhance the accuracy of inverse transformation. Next, for optimizing the objective function, a quasi-Newton algorithm is employed for iterative optimization. During computation, this approach eliminates the need to compute the complex Hessian matrix required in traditional methods, ensuring both efficiency and robustness. Through iterative calculations, a virtual wavefield that satisfies both data constraints and wavefield characteristics is ultimately obtained. Finally, the proposed algorithm is applied to perform wavefield transformation on Transient Electromagnetic (TEM) forward modeling data from layered and complex geological models, validating its effectiveness and stability in practical applications. During the wavefield transformation process, comparative analysis with forward modeling data is conducted to investigate the algorithm's performance under diverse geological models, ensuring broad applicability and accuracy. Numerical experimental results demonstrate that the transformed wavefields retain the propagation characteristics of the wavefield and align closely with the forward modeling results, fully confirming the validity of the inverse transformation algorithm proposed in this study. This study provides a theoretical foundation for interpreting transient electromagnetic quasi-wave equations and holds application value in deep mineral exploration and complex geological structure identification

Key words

Transient electromagnetic / Wavefield transformation / The characteristic constraint of wave field / Quasi-newton optimization / Filtering constraints

引用本文

曹佳慧, 戚志鹏, 李文翰, 李貅. 基于波场特征约束的瞬变电磁波场反变换优化算法研究. 地球物理学进展, 2026 , 41 (1) : 217 -227 . DOI: 10.6038/pg2026II0594

JiaHui CAO, ZhiPeng QI, WenHan LI, Xiu LI. Research on inverse transform optimization algorithm of transient electromagnetic wave field based on wave field feature constraint[J]. Progress in Geophysics, 2026 , 41 (1) : 217 -227 . DOI: 10.6038/pg2026II0594

正文

收起

0　引言

收起

近年来，随着深部探测与矿产资源开发逐步向复杂地形与深部构造区域扩展，工区地表形态和起伏变化愈加剧烈，地下地质条件的复杂性和不确定性也大大增加，由于地下岩层种类繁多、分布不均，且可能存在断层、褶皱等复杂地质构造，这些因素显著提升了施工和开采的难度与成本(陈本池等，1999; 薛国强等, 2008, 2011; 底青云等，2019).在这种背景下，瞬变电磁法凭借其装置灵活，施工便利等优势在复杂环境探测中得到了广泛应用(王其丁等，2024；张继锋等，2024).同时，衍生出的地空瞬变电磁法，凭借较高的工作效率与信噪比等特点，也在这些应用场景中发挥了重要作用，有效提升了勘探精度和工作效率(Xue et al., 2024; Shu et al., 2024).但是低频电磁场遵循扩散方程，地表磁场是地下感应电流的体积扩散特征表现，由于感应电流的扩散效应，导致扩散波场具有明显的体积效应，呈现的是地下介质的综合作用结果，因此在地下结构识别时很难获得精细结构，尤其在解释复杂地质结构时，难以清晰地识别地下几何结构，导致其在具体应用中的解释精度和可靠性不足(Gretskov et al., 2017; Gibert and Virieux, 1991).而波动场具有波的传播特征可以有效识别反射波和折射波，能够直观地描绘地下结构(吴琼，2012; 樊亚楠等，2019；马劼等，2023)，不仅可通过偏振分析、波速方位变换反演各向异性参数，还可以清晰刻画尖锐界面.由于扩散场方程与波动方程之间具有明确的数学转换关系，通过研究场之间的变换，可将电磁扩散场数据转换为拟波动场，结合波动场的界面识别能力和扩散场的深部探测优势，从而更好的实现解释.因此研究扩散场和波动场的变换是必要的(Gershenson, 1997; Maaø, 2007; Mittet, 2010, 2015; 王家映，1995).

Mooney(1982)将波场反变换技术用于加州某一油田折射数据反演，提升了复杂构造区低阻储层边界识别精度，为油田深部勘探提供了理论基础.李貅等(2010)将瞬变电磁数据转换为波场数据，然后使用三维曲面延拓成像技术对野外实际资料进行了处理与解释，得到的结果与实际情况吻合较好.说明瞬变电磁拟地震解释方法是一种有潜力的解释手段.获得高精度的波场就可以实现高精度的解释，因此后续学者将研究聚焦于如何获得高精度的波场之中.

关于波场变换方法的研究，李貅等(2012)提出了一种采用相关叠加技术实现多孔径数据合成的数据处理方法，提高了瞬变电磁法的分辨率.戚志鹏等(2013)将正则化共轭梯度法和预条件共轭梯度法结合起来，使用预条件正则化共轭梯度法，避免了分段叠加，降低了系数矩阵的条件数，提高了波场变换的效率.智庆全(2015)提出了一种新的算法—扫时波场变换，提高了波场变换的稳定性.鲁凯亮等(2021)提出了精细积分法，有效的提高了数据项的精度.范克睿等(2023)提出了虚拟波场降速的方法，通过对波速的控制使波场反变换更加精准高效.波场变换技术具有广阔的应用前景，随着科技的发展进步，发展更加高效、智能、精确的波场变换方法依然是一个重要的课题.薛俊杰等(2018)提出了修正式截断奇异值法，和最近提出的利用神经网络优化的方法都证明了该方法的有效性.但是，传统的波场变换仅考虑了数据约束部分，由于反演的不适定性，以及缺少先验信息，因此前人在构建正则项时采用了平滑约束，即令波场具有平滑特征，从理论上讲平滑约束必然导致波场分辨率的损失.

本文基于前人的研究基础，开展考虑波场特征约束的一种波场反变换方法，将波场的频率特征和传播特征通过设置滤波矩阵的形式作为约束项，构造反变换目标函数，该约束项充分考虑了波场的特征，已知从频域电磁场到波场的转换关系确定波场是带限的，因此在进行波场变换时可基于滤波器设计体现波场的频率分布特征，此外，论文构建了波场的传播特征约束，进一步根据波场特征进行解空间限缩，提高反变换的精度；之后，采用拟牛顿方法进行目标函数优化求解，并利用典型模型波场特征研究该方法的可行性与约束效果.

1　扩散场到波动场转换的基本原理

收起

Kunetz(1972)根据电磁场方程与弹性波动方程的类比特征，发现了扩散方程和波动方程之间存在着数学上的转换关系，Lee等推导出这一变换适用于电磁场任意分量，De Hoop(1996)通过引入尺度变换因子构建了变速度波场变换公式(Lee et al., 1989; Lee and Xie, 1993)，进一步完善了波场变换理论，具体推导如下.

已知扩散场方程

$\boldsymbol{D}(\boldsymbol{r}, t)$

和波动场方程

$\boldsymbol{W}(\boldsymbol{r}, q)$

满足式(1)和式(2)所示的方程：

$\left\{\begin{array}{l}\nabla^{2} \boldsymbol{D}(r, t)-\mu \sigma(r) \frac{\partial}{\partial t} \boldsymbol{D}(r, t)=S(r, t) \\ \boldsymbol{D}(\boldsymbol{r}, 0)=0, \boldsymbol{D}_{\Gamma}=\boldsymbol{D}\left(\boldsymbol{r}_{\mathrm{b}}, t\right), t>0\end{array}\right., $

(1)

$\left\{\begin{array}{l}\nabla^{2} \boldsymbol{W}(\boldsymbol{r}, q)-\mu \boldsymbol{\sigma}(\boldsymbol{r}) \frac{\partial^{2}}{\partial q^{2}} \boldsymbol{W}(\boldsymbol{r}, q)=\boldsymbol{F}(\boldsymbol{r}, q) \\ \boldsymbol{W}(\boldsymbol{r}, 0)=0, \boldsymbol{W}_{\mathrm{r}}=\boldsymbol{W}\left(\boldsymbol{r}_{\mathrm{b}}, q\right), q>0\end{array}\right., $

(2)

其中，

$\boldsymbol{S}$

是扩散场激励，

$\boldsymbol{F}$

为波场激励，

$\varepsilon$

是介电常数，μ是磁导率、σ是电导率，Γ为

$\boldsymbol{r}=\boldsymbol{r}_{\mathrm{b}}$

的边界，t为时间变量，q为虚拟时间.对式(2)进行拉普拉斯变换，令s=p²可得：

$\left\{\begin{array}{l}\nabla^2 \hat{\boldsymbol{D}}\left(\boldsymbol{r}, p^2\right)-\mu \sigma(\boldsymbol{r}) p^2 \hat{\boldsymbol{D}}\left(r, p^2\right)=\hat{\boldsymbol{S}}\left(\boldsymbol{r}, p^2\right) \\\hat{\boldsymbol{D}}_{\Gamma}=\hat{\boldsymbol{D}}\left(\boldsymbol{r}_{\mathrm{b}}, p^2\right), -\frac{\pi}{4}<arg (p)<\frac{\pi}{4} \\\nabla^2 \hat{\boldsymbol{W}}(\boldsymbol{r}, p)-\varepsilon \mu(\boldsymbol{r}) p^2 \hat{\boldsymbol{W}}(\boldsymbol{r}, p)=\hat{\boldsymbol{F}}(\boldsymbol{r}, p) \\\hat{\boldsymbol{W}}_{\Gamma}=\hat{\boldsymbol{W}}\left(\boldsymbol{r}_{\mathrm{b}}, p\right), -\frac{\pi}{2}<\arg (p)<\frac{\pi}{2}\end{array} .\right.$

(3)

令

$\hat{\boldsymbol{D}}\left(\boldsymbol{r}_{\mathrm{b}}, p^{2}\right)=\hat{\boldsymbol{W}}\left(\boldsymbol{r}_{\mathrm{b}}, p\right) ; \hat{\boldsymbol{F}}(\boldsymbol{r}, p)=\hat{\boldsymbol{S}}\left(r, p^{2}\right) ; \varepsilon \mu(\boldsymbol{r})=\mu \sigma(\boldsymbol{r})$

，将式(3)的两个等式左右两边求差，得到：

$\left\{\begin{array}{l}\nabla^{2} \hat{\boldsymbol{G}}(\boldsymbol{r}, p)-\mu \sigma(\boldsymbol{r}) p^{2} \boldsymbol{G}(\boldsymbol{r}, p)=0 \\\hat{\boldsymbol{G}}_{\Gamma}=0, -\frac{\pi}{4}<\arg (p)<\frac{\pi}{4}\end{array}, \right.$

(4)

其中

$\hat{\boldsymbol{G}}(\boldsymbol{r}, p)=\hat{\boldsymbol{D}}\left(\boldsymbol{r}, p^{2}\right)-\hat{W}(\boldsymbol{r}, p)$

.式(4)左右两侧同时乘以

$\hat{\boldsymbol{G}}^{*}(\boldsymbol{r}, p)$

, 并对其积分得到：

$\left\{\begin{array}{l}\int_{\mathrm{V}}|\nabla \times \hat{\boldsymbol{G}}(\boldsymbol{r}, p)|^{2} \mathrm{~d} v+\mu p^{2} \int_{\mathrm{V}} \boldsymbol{\sigma}(r)|\hat{\boldsymbol{G}}(\boldsymbol{r}, p)|^{2} \mathrm{~d} v=0 \\ -\frac{\pi}{4}<\arg (p)<\frac{\pi}{4}\end{array}\right., $

(5)

其中，当

$\hat{\boldsymbol{D}}\left(\boldsymbol{r}, p^{2}\right)=\hat{\boldsymbol{W}}(\boldsymbol{r}, p)$

时式(5)才能成立，忽略空间向量

$\boldsymbol{r}$

，可得：

$\begin{equation*}\int_{0}^{\infty} \boldsymbol{D}(t) \mathrm{e}^{-p^{2} t} \mathrm{~d} t=\int_{0}^{\infty} \boldsymbol{W}(q) \mathrm{e}^{-p q} \mathrm{~d} q, \end{equation*}$

(6)

将p恢复到s的形式可得：

$\begin{equation*}\int_{0}^{\infty} \boldsymbol{D}(t) \mathrm{e}^{-s t} d t=\int_{0}^{\infty} \boldsymbol{W}(q) \mathrm{e}^{-\sqrt{s q}} \mathrm{~d} q . \end{equation*}$

(7)

对式(7)进行拉普拉斯逆变换，可得：

$\begin{equation*}\boldsymbol{D}(t)=\int_{0}^{\infty} \boldsymbol{W}(q) L^{-1}\left(\mathrm{e}^{-\sqrt{s q}}\right) \mathrm{d} q, \end{equation*}$

(8)

已知特殊拉普拉斯变换对

$L\left(\frac{k}{2 \sqrt{\pi t^{3}}} \mathrm{e}^{\frac{-k^{2}}{4 t}}\right)=\mathrm{e}^{-k \sqrt{s}}$

（

$k>0$

），因此式(9)可以写成式(10)的形式，这就是波场变换的积分公式：

$\begin{equation*}\boldsymbol{D}(t)=\frac{1}{2 \sqrt{\pi t^{3}}} \int_{0}^{\infty} q \mathrm{e}^{-q^{2} / 4 t} \boldsymbol{W}(q) \mathrm{d} q .\end{equation*}$

(9)

只要求出该积分公式，便能实现波动场向扩散场的转换，这是一个适定的过程.但是本文研究的是波场反变换，只需要采用数值方法，借助线性插值去对被积函数进行近似处理，把积分问题转变成求和问题，并对核函数的特征展开分析，为后续的反变换做好准备.使用数值方法将式(9)写成式(10)的形式：

$\begin{equation*}\boldsymbol{D}\left(t_{i}\right)=\sum_{j=1}^{n} k_{i, j} \boldsymbol{W}\left(q_{j}\right), \end{equation*}$

(10)

其中：

$\begin{equation*}k_{i, j}=-1 / \sqrt{\pi t_{i}} \mathrm{e}^{-q^{2} / 4 t} \mid{ }_{q_{j}}^{q_{j+1}} .\end{equation*}$

(11)

由

$k_{i, j}$

组成的矩阵

$\boldsymbol{K}$

是离散后的核函数矩阵，为了确定q域的离散方式，分别对比了核函数在不同离散方式下的条件数和形态，表 1分别对比了将核函数在线性离散和对数离散条件下时100阶，200阶和500阶的条件数.

由于野外工作时，瞬变电磁法时间域一般取[10^-4~10^-2 s]，图 1绘制了该时间域范围内的不同离散方式下的核函数展布图来分析该积分核的特征.

由图 1和表 1可以看出，不管什么样的离散方式积分核条件数都很大，说明波场反变换是一个不适定的过程，而在线性离散条件下，核函数条件数较小；随着时间变大，核函数的值变化较大，这有助于后续的计算，因此本文后续讨论基于q域的线性离散.

2　基于拟牛顿算法的波场优化

收起

2.1　优化目标函数构造

已知波场反变换是一个不适定的过程，构造合适的优化目标函数是提高求解质量的重要步骤(彭国民等，2018).由于系数矩阵K通常缺乏良好的条件数，没有办法直接求解，因此需要引入正则化约束构建稳定泛函将反变换问题转换为优化问题实现波场的反变换.传统优化过程中在不知道函数特征情况下往往采用光滑约束构造光滑泛函，该方法先验假设依赖性强，默认假设解具有全局光滑性.本节通过引入波场的频率特征将滤波约束融入到目标函数中，对波场特征得以保留，代替光滑约束.目标函数定义为：

$\begin{equation*}\varPhi(\boldsymbol{W})=\|\boldsymbol{K} \boldsymbol{W}-\boldsymbol{D}\|^{2}+\lambda\|(\boldsymbol{F}-\boldsymbol{I}) \boldsymbol{W}\|^{2}, \end{equation*}$

(12)

$\boldsymbol{F}$

是滤波约束矩阵，控制了所添加约束的意义；

$\boldsymbol{I}$

是单位矩阵；λ是正则化参数，平衡约束项在求解过程中的权重占比.

该正则化项就是滤波约束项，假设需要一个的n阶滤波器h对波场值进行滤波.为了能够将其应用在目标函数并构造滤波约束项，必须将其转化为滤波矩阵 F，使公式(13)成立：

$\begin{align*}\boldsymbol{W}(q) * h(n) & =\sum\nolimits_{m=-\infty}^{+\infty} \boldsymbol{W}(m) h(n-m) \\& =\boldsymbol{F} \boldsymbol{W}(q) , \end{align*}$

(13)

推导可知矩阵 F是一个toeplitz矩阵，可以写成式(14)的形式：

$\boldsymbol{F}=\left[\begin{array}{ccccc}h(1) & 0 & \cdots & 0 & 0 \\h(2) & h(1) & \cdots & \vdots & \vdots \\h(3) & h(2) & \cdots & 0 & 0 \\\vdots & h(3) & \cdots & h(1) & 0 \\h(n-1) & \vdots & \cdots & h(2) & h(1) \\h(n) & h(n-1) & \vdots & \vdots & h(2) \\0 & h(n) & \cdots & h(n-2) & \vdots \\0 & 0 & \cdots & h(n-1) & h(n-2) \\\vdots & \vdots & \vdots & h(n) & h(n-1) \\0 & 0 & 0 & \cdots & h(n)\end{array}\right] .$

(14)

该优化目标函数由两部分组成.第一部分是拟合误差项，用于衡量预测值与实际观测值之间的差异，构成最小二乘法的核心，最小化这一项的目的是尽可能找到一个接近实际数据的解.第二部分是滤波约束项，该项的作用在于筛选目标函数的解，从而在众多可能的解中选出最优解.智庆全提出的扫时波场变换有效提升了波场反变换的稳定性，扫时波场变换通过设定一个时间窗口在时间域内移动，对每个窗口内的场值进行波场变换提高了数据精度，然后将结果相关叠加起来，该过程也具有低通滤波特征，但是叠加可能对数据精度造成破坏.本文的目标函数在此基础上改进而来，但是不需要建立时间窗只需计算出合理的滤波方程，而且将数据约束和滤波约束放在一个方程式中，不仅简化了计算过程、节省了计算时间还保证结果即符合波场的频率特征又满足精度要求.

2.2　约束项参数设计原理

上节指出了约束项的类型，本节着重讨论滤波约束矩阵 F和正则化参数λ如何选取.滤波约束矩阵是约束项的核心，主体滤波器的阶数的大小影响计算效率和滤波效率，截至频率的大小影响滤波效果，因此滤波器参数的选取十分重要.正则化参数起到均衡最小二乘项和滤波约束项的残差的作用，该参数的选取也十分关键，过大会引起波场特征消失，背离了计算的意义；过小会导致约束不起作用，无法得到更优的解(吴国忱和田军，2014; 郭荣文等，2023).本节采用了一对相互对应的波场与扩散场问题的解析解来证明约束项的合理性，具体形式见式(15)和式(16)：

$\boldsymbol{D}=\frac{\mu}{2 \pi} \frac{\mathrm{e}^{-\mu \sigma \rho^{2} / 4 t}}{t},$

(15)

$\boldsymbol{W}=\frac{\mu}{4 \pi} \frac{H(q-\sqrt{\mu \sigma \rho})}{\sqrt{q^{2}-\mu \sigma \rho^{2}}}, $

(16)

其中，

$\rho=\sqrt{x^{2}+z^{2}}, \boldsymbol{H}$

为单位阶跃函数.

由于地震波通常具有较多的低频成分，本文设计了一种FIR(有限脉冲响应)低通滤波器，以有效去除信号中的高频成分，进而使数据更加平滑且清晰.在滤波器设计中，低阶滤波器的滤波效率较低；高阶滤波器则需要更多的计算资源，增加了实时处理的难度和计算时间，因此在选择滤波器阶数时，需要根据具体应用的需求进行权衡.为确保满足设计规格并保持高效性，本文选择了100阶的滤波器并且经过试验确定了正则化参数的最佳范围为[0.1~0.2].滤波器的截止频率是其频率响应的关键特性，通常用于区分通过频率和被衰减频率的范围.如果截止频率设置过大，会导致失真，可能使更多的高频信号通过，尤其是在设计滤波器的目标是去除高频噪声或干扰信号时.这样会让信号中不需要的高频成分得以保留，进而影响最终输出的质量；而截止频率过小，则会滤除过多的信号部分，导致有用的低频成分被削弱或丢失，甚至可能完全丧失一些关键频段，从而影响信号的完整性和可用性.图 2是不同的截止频率下的滤波效果.

由图 2可知，当截至频率为5 Hz时会滤掉有用的信号使结果失真，当截至频率为35 Hz干扰信号过多导致滤波效果较差，通过调试得到最优截止频率为17.5 Hz，该截至频率下反变换后的波场值和解析解基本吻合，可以有效的平滑数据.

值得注意的是，图 2的滤波结果是对波场反变换结果直接卷积滤波器得到的，该方法只适合验证滤波约束的可行性，无法直接用于实地测量中，因此需要对目标函数优化求解.

2.3　基于拟牛顿算法的目标函数优化

梯度法、牛顿法、拟牛顿法等方法都可以对目标函数进行优化，但是他们的工作方式和原理都各有不同，梯度法计算相对简单，对内存需求较低，但是收敛速度较慢，计算效率较低；牛顿法理论上收敛很快，但是由于要计算二阶导数和海森(Hesse)矩阵导致计算成本变大；由于拟牛顿法通常比梯度法更有效并且还可以避免计算海森矩阵，因此本节选用拟牛顿法对目标函数优化.

拟牛顿法是对牛顿法的一种改进，它不直接计算Hesse矩阵的逆，而是通过迭代逐步构造一个近似的矩阵.这使得算法在计算上更加高效，尤其是在大规模问题中(李貅等，2005).构造一个矩阵B来逼近Hesse矩阵，已知拟牛顿法的核心公式为：

$\begin{equation*}\boldsymbol{B}_{k+1}=\boldsymbol{B}_{k}+\frac{y_{k} y_{k}^{\mathrm{T}}}{y_{k}^{\mathrm{T}} s_{k}}-\frac{\boldsymbol{B}_{k} s_{k} s_{k}^{\mathrm{T}} \boldsymbol{B}_{k}}{s_{k}^{\mathrm{T}} \boldsymbol{B}_{k} s_{k}} \end{equation*}.$

(17)

初始矩阵可以设置为单位矩阵或对角矩阵，避免了直接求逆操作，算法的步骤如下所示：

(1) 初始化

$\boldsymbol{W}^{(0)}$

和允许误差ε，计算目标函数

$\varPhi(\boldsymbol{W})$

的梯度，初始化拟牛顿矩阵

$\boldsymbol{B}_{0}$

，设a=0；

(2) 在每次迭代中计算新的搜索方向

$p_{k}=-\boldsymbol{B}_{k} \nabla \varPhi\left(\boldsymbol{W}_{k}\right)$

;

(3) 通过线性搜索确定步长，使目标函数沿着方向p_k下降最快；

(4) 使用公式(17)更新矩阵

$\boldsymbol{B}_{k}, s_{k}=\boldsymbol{W}_{k+1}-\boldsymbol{W}_{k}$

是当前步长，

$y_{k}$

是梯度差；

(5) 判断精度，当满足收敛条件时停止迭代，否者a=a+1，转步骤(2).

通过使用上述拟牛顿法对目标函数进行全局优化，得到在滤波约束和波场特征约束下的虚拟波场结果.图 3展示了优化结果与解析解的对比.

从图 3可以看出，该方法进行变换得到的波场值和理论值吻合较好，相比较不添加滤波项直接优化，无意义的起伏较少、变换结果连续且平滑.

由于波场具有传播过程中的能量为零的特征，因此可以设计一个波场传播特征约束，目标函数可以改写成式(18)的形式：

$\begin{align*}\varPhi(\boldsymbol{W})= & \|\boldsymbol{K} \boldsymbol{W}-\boldsymbol{D}\|^{2}+\lambda_{1}\|(\boldsymbol{F}-\boldsymbol{I}) \boldsymbol{W}\|^{2}+ \\& \lambda_{2}\left\|\int_{0}^{q} W(q) \mathrm{d} q\right\|^{2} .\end{align*}$

(18)

这一改进使得波场的物理特性能得到更好地保留，从而提升了波场反变换在复杂环境下的适应性.为验证该算法的准确性，设计了如图 4所示的简单二层模型.该模型的第一层厚度为200 m，电阻率为100 Ω · m；第二层电阻率为10 Ω · m.使用电性源作为发射源，发射源端点坐标分别为(50，0)和(-50，0)，使用阶跃波作为发射波形，测点范围为[200 m~500 m].

对上述模型正演，并将正演得到的场值进行波场反变换，变换结果如图 5所示.从图 5可以观察到，该二层模型中，地层分界面清晰可辨，存在直达波和下层界面的反射波.由于波场的传播特性，导致随着x的变大，直达波的同相轴随着一条斜率恒定的曲线变化，而反射波的同相轴随着一条斜率逐渐变大的曲线变化.

3　复杂模型波场变换数值模拟

收起

在使用电法或电磁法等技术对地下进行探测，以了解地下结构和储层情况时.当地下存在一些导电性较好的岩石或矿体等区域时，会出现电阻率较低的异常现象.这种异常现象表示存在一些特殊的地质构造或矿化体，如含水层、矿石矿床、破碎带等.低阻异常区的存在可以提供相关地质信息，对于勘探和开发具有一定的指导意义.因此，在地球物理勘探中，对低阻异常区的探测和解释是非常重要的.将上述情形简化为如图 6所示的水平低阻板状体模型，来验证滤波约束算法的准确性，该模型电阻率为1000 Ω · m，在地面以下100 m的地方有一个电阻率为100 Ω · m的异常，该异常长度为400 m宽度为400 m，高度为100 m.使用阶跃波作为发射波形，使用长电性源作为发射源，发射源端点坐标为(600，-600，0)(-600，-600，0).

对上述模型正演，并绘制对应的单点衰减曲线和多测道图，结果如图 7所示.从图 7a可以看出，在时间域[10^-5~10^-2 s]区间内能明显看到异常的特征，绘制该时间域范围内的多测道图，由图 7b可以看出由于中间区域有一个低阻异常体，场值衰减速度较左右两边慢，呈现一个拱形.对该正演结果做波场反变换，变换结果如图 8所示.

由图 8所示的剖面图可以看出，在[-300~-200 m]和[200~300 m]的范围内只有直达波传播，而在[-200~200 m]区域内同时包含了直达波和来自下层界面的反射波，可以清晰地观察到，该区域内存在一个明显的异常体.该异常体的上下分界面十分清晰，且界面过渡明显.

地陷通常指的是地面或地下由于自然或人为原因发生下沉或塌陷的现象，地震、过度开采地下水等行为都有可能导致地面塌陷，该区域通常与矿产资源的富集相关，尤其是一些重要的能源和金属矿床.地陷褶皱也常常是油气藏的重要形成机制之一，褶皱带中的岩层可能作为油气的储层，顶部可能形成良好的盖层，能够聚集油气.因此，对地陷褶皱区域的研究具有重要的科学价值和实际意义.将上述地质条件简化为如图 9所示的地陷褶皱模型来验证滤波约束算法的准确性，该模型第一层层厚为50 m，电阻率为1000 Ω · m；第二层厚度为100 m，电阻率为100 Ω · m；第三层电阻率为1000 Ω · m，在横坐标[-250~250 m]处有一个V型的地陷区域，使用阶跃波作为发射波形，使用电性源作为发射源，发射源端点坐标为(600，-600，0)(-600，-600，0).

对该模型正演，并绘制对应的单点衰减曲线和多测道图，结果如图 10所示.由图 10a可知，在[10^-6~10^-3 s]的时间区间内能看到明显异常，绘制该区间的多测道图，可以明显看到，由于电场在高阻层衰减快，在低阻层衰减慢，所以图 10b呈现一个向下凹陷的形状，对该正演结果做波场反变换，变换结果如图 11所示.

该剖面同时包含了直达波和来自下层界面的反射波，反射波的波形呈现两侧较高、中间较低的特点.由于高频成分会随着地下的深度增加而衰减，因此波场在中间区域的分辨率低于两侧区域.该剖面清晰地反映了地陷层界面的趋势和位置，同时也较好地刻画了凹陷底部与下层界面的变化，整体上更接近真实模型.

4　结论

收起

文中通过构建波场变换目标函数，将波场变换转化为优化问题，通过引入考虑波场频率特征的滤波约束项，并使用拟牛顿法对目标函数实现优化，得到了较好的反变换结果，之后引入了考虑波场传播特征的约束，并使用简单模型证明这一算法的有效性，接着对复杂模型试算表明，该方法能明显找到地层或异常的边界，进一步证明该方法的可行性.

通过目标函数求解得到的虚拟波场，既能保证给定的数据精度，又能准确反映波场的特征；滤波约束和波场传播特征约束的引入能够很好地保留波场的物理特性，得到有效、稳定的波场，从而为后续工作提供了更为可靠的基础.

这个变换算法不仅可以用于波场变换，也可以用于子波变换，子波是激发波形特征，而且瞬变电磁发射波形是可以记录的，这样可以得到确定的子波波形，对于应用波动方程的处理可以等效为已知子波信息，为拟地震处理提出一种新的思路.

参考文献

收起

文献

收起

参考文献引证文献

排序方式：

Chen B C, Zhou F T, Li J M. 1999. The wavefield transformation study of transient electromagnetic field. Geophysical and Geochemical Exploration (in Chinese), 23(3): 195-201.

De Hoop A T. 1996. A general correspondence principle for time-domain electromagnetic wave and diffusion fields. Geophysical Journal International, 127(3): 757-761, doi: 10.1111/j.1365-246X.1996.tb04054.X.

Di Q Y, Zhu R X, Xue G Q, et al. 2019. New development of the Electromagnetic (EM) methods for deep exploration. Chinese J. Geophys. (in Chinese), 62(6): 2128-2138, doi: 10.6038/cjg2019M0633.

Fan K R, Wang G, Lu K L, et al. 2023. Inverse wavefield transform of vertical magnetic signal of galvanic source TEM when slowing down pseudo-wavefield velocity. Chinese J. Geophys. (in Chinese), 66(4): 1743-1757, doi: 10.6038/cjg2022P0796.

Fan Y N, Li X, Qi Z P, et al. 2019. Study on Born approximation algorithm of TEM pseudo wave-field. Progress in Geophysics (in Chinese), 34(2): 529-536, doi: 10.6038/pg2019CC0358.

Gershenson M. 1997. Simple interpretation of time-domain electromagnetic sounding using similarities between wave and diffusion propagation. Geophysics, 62(3): 763-774.

Gibert D, Virieux J. 1991. Electromagnetic imaging and simulated annealing. Journal of Geophysical Research: Solid Earth, 96(B5): 8057-8067, doi: 10.1029/91JB00278.

Gretskov G A, Epov M I, Antonov E Y. 2017. Wave transforms of transient electromagnetic field in conductive earth. Russian Geology and Geophysics, 58(6): 744-751, doi: 10.1016/j.rgg.2016.09.031.

Guo R W, Hou S, Liu J X, et al. 2023. Efficient magnetotelluric forward modeling in the spherical coordinate system based on a regularization constraint. Chinese J. Geophys. (in Chinese), 66(9): 3947-3958, doi: 10.6038/cjg2022Q0651.

Kunetz G. 1972. Processing and interpretation of magnetotelluric soundings. Geophysics, 37(6): 1005-1021, doi: 10.1190/1.1440310.

Lee K H, Liu G, Morrison H F. 1989. A new approach to modeling the electromagnetic response of conductive media. Geophysics, 54(9): 1180-1192.

Lee K H, Xie G Q. 1993. A new approach to imaging with low-frequency electromagnetic fields. Geophysics, 58(6): 780-796.

Li X, Xue G Q, Song J P, et al. 2005. An optimize method for transient electromagnetic field-wave field conversion. Chinese J. Geophys. (in Chinese), 48(5): 1185-1190.

Li X, Qi Z P, Xue G Q, et al. 2010. Three dimensional curved surface continuation image based on TEM pseudo wave-field. Chinese J. Geophys. (in Chinese), 53(12): 3005-3011, doi: 10.3969/j.issn.0001-5733.2010.12.025.

Li X, Xue G Q, Liu Y A, et al. 2012. A research on TEM imaging method based on synthetic-aperture technology. Chinese J. Geophys. (in Chinese), 55(1): 333-340, doi: 10.6038/j.issn.0001-5733.2012.01.034.

Lu K L, Li X, Qi Z P, et al. 2021. A precise integration transform algorithm for transformation from the transient electromagnetic diffusion field into the pseudo wave field. Chinese J. Geophys. (in Chinese), 64(9): 3379-3390, doi: 10.6038/cjg2021O0454.

Ma J, Li X, Qi Z P, et al. 2023. Layered Green function Born approximation imaging method for transient electromagnetic virtual wave fields. Chinese J. Geophys. (in Chinese), 66(6): 2631-2645, doi: 10.6038/cjg2022Q0356.

Maaø F A. 2007. Fast finite-difference time-domain modeling for marine-subsurface electromagnetic problems. Geophysics, 72(2): A19-A23.

Mittet R. 2010. High-order finite-difference simulations of marine CSEM surveys using a correspondence principle for wave and diffusion fields. Geophysics, 75(1): F33-F50.

Mittet R. 2015. Seismic wave propagation concepts applied to the interpretation of marine controlled-source electromagnetics. Geophysics, 80(2): E63-E81.

Mooney W D. 1982. Application of Wave Field Continuation to the Inversion of Refraction Data. Geophysical Journal International, 70(1), 199-216, doi: 10.1111/j.1365-246X.1982.tb06590.x.

Peng G M, Liu Z, Yu H Z, et al. 2018.3D gravity inversion based on Cauchy distribution constraint and fast proximal objective function optimization. Chinese J. Geophys. (in Chinese), 61(12): 4934-4941, doi: 10.6038/cjg2018L0776.

Qi Z P, Li X, Wu Q, et al. 2013. A new algorithm for full-time-domain wave-field transformation based on transient electromagnetic method. Chinese J. Geophys. (in Chinese), 56(10): 3581-3595, doi: 10.6038/cjg20131033.

Shu C C, Qi Y F, Yang S Y, et al. 2024.3-D forward modeling of semi-airborne TEM method with IP effect based on goal-oriented adaptive finite element algorithm. IEEE Transactions on Geoscience and Remote Sensing, 62: 5926211, doi: 10.1109/TGRS.2024.3454113.

Wang J Y. 1995. Magnetotelluric Sounding Interpretation Method (in Chinese). Beijing: Petroleum Industry Press.

Wang Q D, Li W H, Liu W T, et al. 2024. Research on multi-source transient electromagnetic imaging method in TBM tunnel environment. Progress in Geophysics (in Chinese), 39(6): 2463-2470, doi: 10.6038/pg2024HH0478.

Wu G C, Tian J. 2014. AVO inversion based on scale constraints. Oil Geophysical Prospecting (in Chinese), 49(6): 1157-1164, 1190.

Wu Q. 2012. Forward modeling of large loop transient electromagnetic field and wave-field transform theory's study [Master's thesis](in Chinese). Xi'an: Chang'an University.

Xue G Q, Li X, Di Q Y. 2008. Research progress in TEM forward modeling and inversion calculation. Progress in Geophysics (in Chinese), 23(4): 1165-1172.

Xue G Q, Li X, Qi Z P, et al. 2011. Study of sharpening the TEM pseudo-seismic wave-form. Chinese J. Geophys. (in Chinese), 54(5): 1384-1390, doi: 10.3969/j.issn.0001-5733.2011.05.028.

Xue G Q, Zhou N N, Su B X, et al. 2024. Geophysical exploration strategy for Cu-Ni-Co deposits in China: A review. Geophysics, 89(1): WB25-WB34.

Xue J J, Zhong H S, Li H, et al. 2018. An improvement to the transformation algorithm for the transient electromagnetic field. Chinese J. Geophys. (in Chinese), 61(12): 5077-5083, doi: 10.6038/cjg2018M0271.

Zhang J F, Wang S, Sun N Q, et al. 2024. Transient electromagnetic 3D numerical simulation and application of underground pipeline leakage. Progress in Geophysics (in Chinese), 39(3): 1206-1217, doi: 10.6038/pg2024HH0354.

Zhi Q Q. 2015. Study on wave field transformation and migration imaging of MTEM data [Master's thesis] (in Chinese). Xi'an: Chang'an University.

陈本池, 周凤桐, 李金铭. 1999. 瞬变电磁场的波场变换研究. 物探与化探, 23(3): 195-201.

底青云, 朱日祥, 薛国强, 等. 2019. 我国深地资源电磁探测新技术研究进展. 地球物理学报, 62(6): 2128-2138, doi: 10.6038/cjg2019M0633.

范克睿, 王刚, 鲁凯亮, 等. 2023. 虚拟波场降速条件下电性源瞬变电磁垂直磁场分量的波场反变换方法. 地球物理学报, 66(4): 1743-1757, doi: 10.6038/cjg2022P0796.

樊亚楠, 李貅, 戚志鹏, 等. 2019. 瞬变电磁虚拟波场Born近似算法研究. 地球物理学进展, 34(2): 529-536, doi: 10.6038/pg2019CC0358.

郭荣文, 侯生, 柳建新, 等. 2023. 基于正则化约束项的球坐标系大地电磁快速正演. 地球物理学报, 66(9): 3947-3958, doi: 10.6038/cjg2022Q0651.

李貅, 薛国强, 宋建平, 等. 2005. 从瞬变电磁场到波场的优化算法. 地球物理学报, 48(5): 1185-1190.

李貅, 戚志鹏, 薛国强, 等. 2010. 瞬变电磁虚拟波场的三维曲面延拓成像. 地球物理学报, 53(12): 3005-3011, doi: 10.3969/j.issn.0001-5733.2010.12.025.

李貅, 薛国强, 刘银爱, 等. 2012. 瞬变电磁合成孔径成像方法研究. 地球物理学报, 55(1): 333-340, doi: 10.6038/j.issn.0001-5733.2012.01.034.

鲁凯亮, 李貅, 戚志鹏, 等. 2021. 瞬变电磁扩散场到虚拟波场的精细积分变换算法. 地球物理学报, 64(9): 3379-3390, doi: 10.6038/cjg2021O0454.

马劼, 李貅, 戚志鹏, 等. 2023. 瞬变电磁虚拟波场层状Green函数Born近似成像方法. 地球物理学报, 66(6): 2631-2645, doi: 10.6038/cjg2022Q0356.

彭国民, 刘展, 于会臻, 等. 2018. 基于柯西分布约束和快速近端目标函数优化的三维重力反演方法. 地球物理学报, 61(12): 4934-4941, doi: 10.6038/cjg2018L0776.

戚志鹏, 李貅, 吴琼, 等. 2013. 从瞬变电磁扩散场到拟地震波场的全时域反变换算法. 地球物理学报, 56(10): 3581-3595, doi: 10.6038/cjg20131033.

王家映. 1995. 大地电磁拟地震解释法. 北京: 石油工业出版社.

王其丁, 李文翰, 刘文韬, 等. 2024. TBM隧道环境下的多源瞬变电磁成像方法研究. 地球物理学进展, 39(6): 2463-2470, doi: 10.6038/pg2024HH0478.

吴国忱, 田军. 2014. 基于尺度约束项的AVO反演方法. 石油地球物理勘探, 49(6): 1157-1164, 1190.

吴琼. 2012. 大回线源电磁场正演与波场变换理论研究[硕士论文]. 西安: 长安大学.

薛国强, 李貅, 底青云. 2008. 瞬变电磁法正反演问题研究进展. 地球物理学进展, 23(4): 1165-1172.

薛国强, 李貅, 戚志鹏, 等. 2011. 瞬变电磁拟地震子波宽度压缩研究. 地球物理学报, 54(5): 1384-1390, doi: 10.3969/j.issn.0001-5733.2011.05.028.

薛俊杰, 钟华森, 李海, 等. 2018. 瞬变电磁波场转换算法的改进. 地球物理学报, 61(12): 5077-5083, doi: 10.6038/cjg2018M0271.

张继锋, 王尚, 孙乃泉, 等. 2024. 地下管道渗漏的瞬变电磁三维数值模拟及应用. 地球物理学进展, 39(3): 1206-1217, doi: 10.6038/pg2024HH0354.

智庆全. 2015. MTEM波场变换与偏移成像方法研究[硕士论文]. 西安: 长安大学.

2026年第41卷第1期

PDF下载

210

引用本文

BibTeX

文章信息

doi: 10.6038/pg2026II0594

接收时间：2025-04-22
首发时间：2026-09-23
出版时间：2026-02-20

补充材料

相关文章

文章信息

作者

出版历史

收稿日期：2025-04-22
修回日期：2025-11-08

基金

作者信息

¹长安大学地质工程与测绘学院, 西安 710054

²长安大学地球物理场多参数综合模拟实验室(中国地球物理学会重点实验室), 西安 710054

通讯作者:

戚志鹏, 男, 1982年生, 博士生导师, 教授, 主要从事瞬变电磁场的理论与应用方面的研究. E-mail: qzhipeng@126.com

参考文献

分享链接

https://castjournals.cast.org.cn/joweb/dqwlxjz/CN/10.6038/pg2026II0594

分享至

全文二维码

扫描看全文

引用本文

BibTeX

本文的引用情况

2种不同金属材料的力学参数

科 Family	属数 Number of genus	种数 Number of species	占总种数比例 Percentage of total species (%)	属 Genus	种数 Number of species	占总种数比例 Percentage of total species (%)
鹅膏菌科Amanitaceae	2	11	5.26	鹅膏菌属 Amanita	10	4.78
小菇科 Mycenaceae	2	12	5.74	丝盖伞属 Inocybe	5	2.39
多孔菌科 Polyporaceae	8	14	6.70	蜡蘑属 Laccaria	5	2.39
红菇科 Russulaceae	3	23	11.00	小皮伞属 Marasmius	6	2.87
				小菇属 Mycena	11	5.26
				光柄菇属 Pluteus	5	2.39
				红菇属 Russula	17	8.13
				栓菌属 Trametes	5	2.39

关闭全屏

BibTeX
EndNote
RefWorks
TxT