CN111474594B - 一种三维时间域航空电磁快速反演方法 - Google Patents

一种三维时间域航空电磁快速反演方法 Download PDF

Info

Publication number
CN111474594B
CN111474594B CN202010462730.9A CN202010462730A CN111474594B CN 111474594 B CN111474594 B CN 111474594B CN 202010462730 A CN202010462730 A CN 202010462730A CN 111474594 B CN111474594 B CN 111474594B
Authority
CN
China
Prior art keywords
inversion
scale
data
time domain
wavelet
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Expired - Fee Related
Application number
CN202010462730.9A
Other languages
English (en)
Other versions
CN111474594A (zh
Inventor
齐彦福
戚志鹏
李貅
周建美
孙乃泉
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Changan University
Original Assignee
Changan University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Changan University filed Critical Changan University
Priority to CN202010462730.9A priority Critical patent/CN111474594B/zh
Publication of CN111474594A publication Critical patent/CN111474594A/zh
Application granted granted Critical
Publication of CN111474594B publication Critical patent/CN111474594B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V3/00—Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation
    • G01V3/38—Processing data, e.g. for analysis, for interpretation, for correction

Abstract

本发明公开了一种三维时间域航空电磁快速反演方法，该方法首先构造适用于时间域航空电磁数据的最佳第二代小波，对时间域航空电磁观测数据进行多尺度分解；自适应选取最优化小波阈值，提取当前尺度的特征数据；进行模型灵敏度多尺度小波分析，自适应优化时间域航空电磁反演网格；利用当前尺度上提取的特征数据和反演网格进行时间域航空电磁三维高斯‑牛顿反演；通过逐步降低数据空间和模型空间尺度，分步提取不同尺度上的时间域航空电磁特征数据参与反演，同时变尺度更新反演网格，实现时间域航空电磁数据快速精细化反演。在保证分辨率的同时可以有效提高时间域航空电磁三维反演效率。

Description

一种三维时间域航空电磁快速反演方法
技术领域
本发明涉及应用地球物理领域，具体涉及一种三维时间域航空电磁快速反演方法。
背景技术
时间域航空电磁法采用机载移动平台和瞬变电磁回线装置对地下目标实现无损探测，其具有工作速度快、分辨率高、无需地面人员接近等优点，非常适合在难以开展地面地球物理勘查工作的高山、沙漠、沼泽等地区进行快速电磁勘探，已被广泛应用于矿产资源勘查、地下水和地热资源调查以及工程环境等领域。
目前，时间域航空电磁数据解释技术主要有一维成像和三维反演两种。由于传统的一维成像方法基于水平层状模型假设，当遇到地形起伏严重、三维电性结构发育的山区时，无法揭示地下真实的电性结构，造成时间域航空电磁方法优势得不到充分发挥。而三维反演方法基于真实三维地电模型，通过对地下进行网格剖分精细刻画三维地电结构，可以有效处理地形起伏严重地区的时间域航空电磁数据，是航空电磁勘探领域前沿热点。然而，时间域航空电磁法采样密集，数据量大，造成时间域航空电磁三维计算量大、反演效率低。由于计算成本高，三维反演无法实用化，成为困扰时域航空电磁技术推广的技术难题。
对于如何有效提高时间域航空电磁三维反演效率这一技术问题，目前尚未见有成熟的技术方案。
发明内容
针对现有技术存在的缺陷或不足，本发明的目的在于，提供一种三维时间域航空电磁快速反演方法，该方法通过分步提取多尺度特征数据和变尺度更新反演网格，在保证分辨率的同时提高时间域航空电磁三维反演效率。
为了实现上述任务，本发明采取如下的技术解决方案：
步骤一，构造一种能够利用少量小波系数准确捕捉时间域航空电磁数据变化特征的第二代小波，对时间域航空电磁数据进行多尺度分解，将观测数据分解成多个包含不同尺度信息的数据子集；
步骤二，设定采样尺度，通过自适应搜索最优化小波阈值，提取当前尺度上的时间域航空电磁特征数据参与反演；
步骤三，利用第二代小波对当前反演网格模型灵敏度进行多尺度小波分解，利用小波系数识别电性突变的异常目标区域，通过设置小波阈值，对当前尺度反演网格进行局部调整，将精细网格聚焦于电性异常区域，自适应优化当前尺度反演网格；
步骤四：在当前尺度上，利用提取的特征数据和当前反演网格构建反演目标函数和反演方程，利用高斯牛顿法更新三维反演模型；具体为：
解耦正演和反演网格，利用局部网格技术建立多个小规模正演网格，降低正演问题的求解维数；
采用非结构时间域有限元算法计算正演响应，同时采用伴随正演算法显式计算灵敏度矩阵，通过结合并行技术加快正演和灵敏度矩阵计算速度；
步骤五：逐步降低数据和模型尺度，并重复步骤二至步骤四，分步提取不同尺度上的特征数据参与反演，同时变尺度更新反演网格，最终实现时间域航空电磁数据快速精细化反演。
根据本发明，步骤一中，所述的构造一种能够利用少量小波系数准确捕捉时间域航空电磁数据变化特征的第二代小波的方法为：根据观测数据随空间水平位置和时间的三维变化特征，设置最佳消失矩阶数，基于提升原理，通过求解线性方程组获得提升参数，在现有插值小波基础上构建具有高阶消失矩的第二代小波。
进一步的，在步骤二中所述自适应搜索最优化小波阈值的方法为：在初始小波阈值基础上，根据重构数据与原始数据之间的差异自动调整阈值大小，从而获得满足数据重构精度的最优化小波阈值。
进一步的，在步骤五中所述逐步降低数据和模型尺度的时机是在当前模型约束强度下，利用当前尺度抽样数据和反演网格反演收敛时，降低数据和模型尺度，同时减弱模型约束强度。
进一步的，在步骤五中所述的分步提取不同尺度上的特征数据参与反演时，首先在大尺度上抽取少量特征数据快速获得地下大尺度构造，随着数据采样尺度的降低，在小尺度上补充提取包含细节信息的特征数据，获得地下精细结构。
进一步的，步骤五中所述的变尺度更新反演参数网格的方法为：
在反演初始阶段模型约束较强的情况下，首先在大尺度上抽取少量特征数据，在粗糙的反演网格上快速获得地下电性结构的基本轮廓；
当反演收敛时，减弱模型约束强度，同时降低数据采样尺度，在小尺度上补充提取包含细节信息的特征数据；
与此同时，降低反演网格尺度，根据小波阈值局部细化反演网格，提高对小尺度精细地电结构的分辨能力。通过同步降低数据和模型尺度，逐步获得地下精细结构。
本发明的三维时间域航空电磁快速反演方法，从以下两个层面上在保证分辨率的同时提高了反演效率：
(1)在反演效率方面，通过对时间域航空电磁数据进行多尺度分解和自适应小波阈值选取，分步提取不同尺度上的特征数据参与反演，令更多的采样点聚焦到异常区域，且少量数据包含了数据的主体特征信息，减少了每次迭代的反演数据量，从而有效地提高了反演效率。
(2)在反演精度方面，通过模型灵敏度小波分析，对反演网格进行局部优化，将精细网格聚焦于电性异常区域，而对于电性均匀区域采用粗糙网格，从而更加精细地刻画了反演模型通过逐步局部细化反演网格尺度，保证了对小尺度精细地电结构的分辨能力。
附图说明
图1是本发明的三维时间域航空电磁快速反演方法实施流程图。
图2是自适应小波阈值搜索流程图。
图3是自适应反演网格优化流程图。
图4是变尺度网格更新效果图。
下面结合附图和实施例对本发明作进一步的详细说明。
具体实施方式
本发明的技术思路是，首先在大尺度上提取少量特征数据，利用粗糙反演网格快速反演地下大尺度宏观电性结构，并将当前尺度反演结果作为下一尺度的初始模型。然后，降低数据采样尺度，在小尺度上补充采样细节特征数据，同时根据模型灵敏度多尺度分析局部细化电性异常区域，更新反演网格，反演地下小尺度细节构造。最终通过逐步降低反演尺度，实现时间域航空电磁数据快速精细化反演。
如图1所示，本实施例给出一种三维时间域航空电磁快速反演方法，具体步骤为：
步骤(1)：根据航空电磁变化规律，构造一种能够利用少量小波系数准确捕捉时间域航空电磁数据变化特征的第二代小波，通过小波分析对观测数据进行多尺度分解，将观测数据分解成多个包含不同尺度信息的数据子集；
步骤(2)：根据模型约束强度，设定采样尺度，通过自适应选取小波阈值，提取当前尺度上的特征数据参与反演；
步骤(3)：利用第二代小波对当前反演网格模型灵敏度进行多尺度小波分解，通过设置小波阈值，自适应优化当前尺度反演网格；
步骤(4)：在当前尺度上，利用提取的特征数据和当前尺度反演网格进行时间域航空电磁三维高斯-牛顿反演，更新反演模型；
步骤(5)：逐步降低数据和模型尺度，并重复步骤2至步骤4，分步提取不同尺度上的特征数据参与反演，同时变尺度更新反演网格，最终实现时间域航空电磁数据快速精细化反演。
其中，步骤(1)中，第二代小波的构造基于提升理论在现有小波基础上建立提升小波：
Figure BDA0002511497350000051
Figure BDA0002511497350000052
其中，
Figure BDA0002511497350000053
和
Figure BDA0002511497350000054
分别为现有的尺度函数和小波函数，
Figure BDA0002511497350000055
和
Figure BDA0002511497350000056
分别是提升后的尺度函数和小波函数，
Figure BDA0002511497350000057
是提升参数，j表示小波层数，k表示尺度函数编号，m表示小波函数编号。
根据航空电磁数据的变化特征，设定消失矩为N阶，则提升小波系数需满足：
Figure BDA0002511497350000058
其中，x为自变量，X为自变量的值域。
通过结合(2)式和(3)式，建立关于提升系数
Figure BDA0002511497350000059
的方程组：
Figure BDA0002511497350000061
式中，所用变量含义与前述一致，求解方程组获得
Figure BDA0002511497350000062
根据(2)式即可完成提升小波构造。
步骤(2)中自适应小波阈值的选取方法如图2所示，首先设定初始小波阈值和重构误差阈值ε0，然后基于数据多尺度分析结果，将最大尺度到当前尺度所有小于阈值的小波系数置0，通过小波逆变换获得重构数据误差ε，并计算重构数据与原始数据的拟合差，如果数据重构误差ε大于重构误差阈值ε0，则说明数据采样太少，部分信息丢失，此时增大小波阈值；如果数据拟合差小于阈值，则表示利用当前数据量就可以恢复原始数据所包含的异常信息。
步骤(3)中，自适应反演网格优化过程如图3所示，首先计算所有数据点对每个反演模型单元的灵敏度累加和
Figure BDA0002511497350000063
其中n为数据点个数，Ji,j是第j个数据对第i个单元的灵敏度。然后将每个单元的灵敏度累加和
Figure BDA0002511497350000064
组装成三维矩阵
Figure BDA0002511497350000065
并利用第二代小波对累加灵敏度矩阵
Figure BDA0002511497350000066
进行多尺度分解，获得所有尺度上的小波系数和尺度系数，即：
Figure BDA0002511497350000067
其中，fd表示进行d层小波变换，S是尺度系数，φ是尺度函数，W是小波系数，ψ是小波函数，k是小波变换的层编号，q是小波块编号，m、n和l分辨为三个维度上的小波系数编号。
最后设定上限小波阈值Wu和下限小波阈值Wl，通过对小波系数大于上限小波阈值Wu的网格单元进行局部加密，而对小于下限小波阈值Wl的网格进行临近单元合并实现反演网格优化。
步骤(4)中综合考虑数据拟合和模型光滑约束建立反演目标函数如下：
Figure BDA0002511497350000071
其中，m是反演模型参数向量，mref是参考模型向量，A是正演方程系数矩阵，b是正演方程源项向量，Q是插值矩阵，dobs是观测数据向量，Wd是数据方差矩阵，Wm是模型方差矩阵，β是拉格朗日乘子。
通过求取目标函数极值可以得到反演控制方程：
Figure BDA0002511497350000072
其中，J是灵敏度矩阵，上角标T表示矩阵的转置，Δm是模型修正向量，g(m)是目标函数φ(m)的梯度向量，由下式给出：
Figure BDA0002511497350000073
采用显式方法计算灵敏度矩阵J，并由下式给出：
JT＝GTv (8)
其中，矩阵
Figure BDA0002511497350000074
列向量e是网格棱边的电场向量值，v是一个列向量，通过求解正演方程ATv＝QT获得。
反演过程中解耦正演和反演网格，并在反演网格基础上建立局部正演网格。在局部网格上利用伴随正演算法显示计算灵敏度矩阵，然后利用正演和反演网格空间映射关系，重组反演网格灵敏度矩阵。最后利用(6)式计算模型修正量，更新反演模型。
步骤(5)中，在当前模型约束强度下反演模型收敛时，降低数据和模型尺度，重复步骤(2)至步骤(4)，在小尺度上补充采样包含细节信息的特征数据，同时根据采样数据优化并更新反演网格，在当前尺度进行反演更新反演模型。
步骤(5)中，分步数据提取的方法为在反演初始阶段模型约束强度较强时，在大尺度上进行特征数据提取，获得宏观构造信息，随着模型约束强度减弱，逐步降低数据采样尺度，在小尺度上进行特征数据采样，获得细节信息。
步骤(5)中，变尺度更新反演参数网格的方法为首先利用粗糙的大尺度反演网格快速获得地下电性结构的基本轮廓，并将该尺度的反演结果作为下一尺度的反演初值，随着模型约束减弱、采样数据增加，逐步局部细化反演网格尺度，提高对小尺度精细地电结构的分辨能力，图4给出了变尺度网格更新效果图。

Claims (4)

1.一种三维时间域航空电磁快速反演方法，其特征在于，具体步骤如下：
步骤一，构造一种能够利用少量小波系数准确捕捉时间域航空电磁数据变化特征的第二代小波，对时间域航空电磁数据进行多尺度分解，将观测数据分解成多个包含不同尺度信息的数据子集；
所述的构造一种能够利用少量小波系数准确捕捉时间域航空电磁数据变化特征的第二代小波的方法为：根据观测数据随空间水平位置和时间的三维变化特征，设置最佳消失矩阶数，基于提升原理，通过求解线性方程组获得提升参数，在现有插值小波基础上构建具有高阶消失矩的第二代小波；
步骤二，设定采样尺度，通过自适应搜索最优化小波阈值，提取当前尺度上的时间域航空电磁特征数据参与反演；
步骤三，利用第二代小波对当前反演网格模型灵敏度进行多尺度小波分解，利用小波系数识别电性突变的异常目标区域，通过设置小波阈值，对当前尺度反演网格进行局部调整，将精细网格聚焦于电性异常区域，自适应优化当前尺度反演网格；
步骤四：在当前尺度上，利用提取的特征数据和当前反演网格构建反演目标函数和反演方程，利用高斯牛顿法更新三维反演模型；具体为：
解耦正演和反演网格，利用局部网格技术建立多个小规模正演网格，降低正演问题的求解维数；
采用非结构时间域有限元算法计算正演响应，同时采用伴随正演算法显式计算灵敏度矩阵，通过结合并行技术加快正演和灵敏度矩阵计算速度；
步骤五：逐步降低数据和模型尺度，并重复步骤二至步骤四，分步提取不同尺度上的特征数据参与反演，同时变尺度更新反演网格，最终实现时间域航空电磁数据快速精细化反演；所述的变尺度更新反演网格的方法为：
在反演初始阶段模型约束较强的情况下，首先在大尺度上抽取少量特征数据，在粗糙的反演网格上快速获得地下电性结构的基本轮廓；
当反演收敛时，减弱模型约束强度，同时降低数据采样尺度，在小尺度上补充提取包含细节信息的特征数据；
与此同时，降低反演网格尺度，根据小波阈值局部细化反演网格，提高对小尺度精细地电结构的分辨能力；通过同步降低数据和模型尺度，逐步获得地下精细结构。
2.如权利要求1所述的方法，其特征在于，步骤二中所述的自适应搜索最优化小波阈值的方法为：在初始小波阈值基础上，根据重构数据与原始数据之间的差异自动调整阈值大小，从而获得满足数据重构精度的最优化小波阈值。
3.如权利要求1所述的方法，其特征在于，步骤五中所述的逐步降低数据和模型尺度的时机是在当前模型约束强度下，利用当前尺度抽样数据和反演网格反演收敛时，降低数据和模型尺度，同时减弱模型约束强度。
4.如权利要求1所述的方法，其特征在于，步骤五中所述的分步提取不同尺度上的特征数据参与反演时，首先在大尺度上抽取少量特征数据快速获得地下大尺度构造，随着数据采样尺度的降低，在小尺度上补充提取包含细节信息的特征数据，获得地下精细结构。
CN202010462730.9A 2020-05-27 2020-05-27 一种三维时间域航空电磁快速反演方法 Expired - Fee Related CN111474594B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010462730.9A CN111474594B (zh) 2020-05-27 2020-05-27 一种三维时间域航空电磁快速反演方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010462730.9A CN111474594B (zh) 2020-05-27 2020-05-27 一种三维时间域航空电磁快速反演方法

Publications (2)

Publication Number Publication Date
CN111474594A CN111474594A (zh) 2020-07-31
CN111474594B true CN111474594B (zh) 2021-08-17

Family

ID=71764878

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010462730.9A Expired - Fee Related CN111474594B (zh) 2020-05-27 2020-05-27 一种三维时间域航空电磁快速反演方法

Country Status (1)

Country Link
CN (1) CN111474594B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112949134B (zh) * 2021-03-09 2022-06-14 吉林大学 基于非结构有限元方法的地-井瞬变电磁反演方法
CN113325482B (zh) * 2021-04-15 2024-01-16 成都理工大学 一种时间域电磁数据反演成像方法
CN115113286B (zh) * 2022-07-06 2023-09-15 长江大学 基于多分量频率域航空电磁数据融合三维反演方法

Family Cites Families (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CA2616390C (en) * 2005-07-28 2014-04-08 Exxonmobil Upstream Research Company Method for wavelet denoising of controlled source electromagnetic survey data
JP5596982B2 (ja) * 2010-01-08 2014-10-01 キヤノン株式会社 電磁波の測定装置及び方法
US9377548B2 (en) * 2011-11-09 2016-06-28 Chevron U.S.A. Inc. Wavelet-transform based system and method for analyzing characteristics of a geological formation
US10310138B2 (en) * 2015-09-15 2019-06-04 Exxonmobil Upstream Research Company Accelerated Occam inversion using model remapping and Jacobian matrix decomposition
CN110068816B (zh) * 2019-05-10 2023-04-28 长沙理工大学 一种基于提升格式的探地雷达信号用小波基构造方法
CN110058317B (zh) * 2019-05-10 2020-09-08 成都理工大学 航空瞬变电磁数据和航空大地电磁数据联合反演方法
CN110244275A (zh) * 2019-07-19 2019-09-17 上海交通大学 海杂波双谱重构与仿真方法

Also Published As

Publication number Publication date
CN111474594A (zh) 2020-07-31

Similar Documents

Publication Publication Date Title
CN111474594B (zh) 一种三维时间域航空电磁快速反演方法
CN110568486B (zh) 基于同步稀疏低秩张量补全模型的地震信号补全方法
Kollat et al. A computational scaling analysis of multiobjective evolutionary algorithms in long-term groundwater monitoring applications
Jiang et al. Electrical resistivity imaging inversion: An ISFLA trained kernel principal component wavelet neural network approach
CN112949134B (zh) 基于非结构有限元方法的地-井瞬变电磁反演方法
AU2014296770A1 (en) System and method for estimating a reservoir parameter using joint stochastic inversion of multisource geophysical data
Minasny et al. Pedometrics research in the vadose zone—Review and perspectives
CN114331842B (zh) 结合地形特征的dem超分辨率重建方法
Ling et al. Learning-based superresolution land cover mapping
CN108983300B (zh) 一种隧道掘进机施工条件下的瞬变电磁隧道超前预报方法
Gündoğdu et al. Three-dimensional regularized inversion of DC resistivity data with different stabilizing functionals
CN114065511A (zh) 起伏地形下大地电磁二维正演数值模拟方法、装置、设备及介质
Feng et al. Resistivity-depth imaging with the airborne transient electromagnetic method based on an artificial neural network
CN113569447B (zh) 一种基于舒尔补方法的瞬变电磁三维快速正演方法
Tao et al. Seismic surface-related multiples suppression based on SAGAN
Pessel et al. Multiscale electrical impedance tomography
Gao et al. 3D visualization monitoring and early warning of surface deformation in subsidence area based on GIS
Chen et al. Even sampling designs generation by efficient spatial simulated annealing
CN111929348A (zh) 水平多层结构大地环境下直流接地极地表电位计算方法
Peng et al. Three-dimensional spatial prediction of Zn in the soil of a former tire manufacturing plant using machine learning and readily attainable multisource auxiliary data
Wang et al. Hazard assessment of debris flows based on a PCA-GRNN model: A case study in Liaoning Province, China
Yan et al. Improved Tucker decomposition algorithm for noise suppression of 3D GPR data in road detection
CN110794469B (zh) 基于最小地质特征单元约束的重力反演方法
Jia et al. Magnetotelluric Closed-Loop Inversion
Moazam et al. The priority of microgravity focusing inversion in 3D modeling of subsurface voids

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20210817