CN116090283A - 基于压缩感知和预条件随机梯度的航空电磁三维反演方法 - Google Patents

基于压缩感知和预条件随机梯度的航空电磁三维反演方法 Download PDF

Info

Publication number
CN116090283A
CN116090283A CN202211414788.1A CN202211414788A CN116090283A CN 116090283 A CN116090283 A CN 116090283A CN 202211414788 A CN202211414788 A CN 202211414788A CN 116090283 A CN116090283 A CN 116090283A
Authority
CN
China
Prior art keywords
inversion
data
dimensional
model
matrix
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.)
Pending
Application number
CN202211414788.1A
Other languages
English (en)
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.)
Jilin University
Original Assignee
Jilin 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 Jilin University filed Critical Jilin University
Priority to CN202211414788.1A priority Critical patent/CN116090283A/zh
Publication of CN116090283A publication Critical patent/CN116090283A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Algebra (AREA)
  • Software Systems (AREA)
  • Databases & Information Systems (AREA)
  • Operations Research (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • Computing Systems (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明属于地球物理电磁法反演技术领域，具体为基于压缩感知和预条件随机梯度的航空电磁三维反演方法，包括以下步骤：步骤一，筛选航空电磁实测数据中的可用数据用于反演；步骤二，使用泊松圆盘采样法对测点进行随机采样，并对采样点利用有限体积法进行正演计算得到随机测点的预测数据，并使用压缩感知方法对所有测点的预测数据进行重构；步骤三，将步骤一的实测数据和步骤二的预测数据做数据拟合，并结合模型粗糙度构建航空电磁三维正则化反演的目标函数，收敛速度可以媲美传统全批量数据三维反演，同时在保证反演精度的同时极大的提升反演计算效率。

Description

基于压缩感知和预条件随机梯度的航空电磁三维反演方法
技术领域
本发明涉及地球物理电磁法反演技术领域，具体为基于压缩感知和预条件随机梯度的航空电磁三维反演方法。
背景技术
我国国土面积近三分之二为地面人员难以接近的无人区，特别是广大西部地区，高山、沙漠和大面积森林覆盖导致这些地区地球物理勘探程度很低。航空电磁勘查技术是基于固定翼飞机、直升机或无人机等搭载平台的地球物理勘查技术。通过机载或吊挂不接地回线作为磁性发射源并采用磁传感器接收地下电磁信号以达到勘探的目的，具有采样速度快、勘探效率高、无需地面人员接近的特点。航空电磁法的这些优势完美地解决了无人区无法勘探的问题，因此发展航空电磁法是大势所趋。
航空电磁反演目的是通过在机载接收到的电磁响应来反推地下空间的真实地电模型，并以此辅助地质解释工作者得到更为准确的地质结构判断。对于航空电磁法，国内外学者做了大量的研究，由于航空电磁数据量庞大，数据反演解释方法通常采用成像或一维反演，但是成像精度低，一维反演具有极大的局限性，现如今已无法满足对地下精细结构的刻画。所各种优化算法被引用到三维航空电磁反演问题中，目前流行的包括：高斯牛顿法(Mackie and Madden，1993；Newman and Alumbaugh,2000；Sasaki et al.,2013)、拟牛顿法(Bae et al.,2012；刘云鹤等，2013)和非线性共轭梯度法(刘云鹤等，2013；Kamm andPedersen，2014)等。
传统三维反演算法需要对测区或目标区所有测点数据进行计算，同时需要大量的伴随正演获得灵敏度矩阵，耗时长、效率低。随着机器学习理论的不断崛起，随机梯度法逐渐受到关注，该方法是在机器学习训练中，随机选取某一个或者一批样本进行训练，代替传统方法全批量计算，在达到训练目的的同时大大提高训练效率。但该方法由于数据量少使得梯度方向并不是最优的，导致反演收敛速度慢，同时会引入梯度噪声，导致其可能反复跳出最优解，而这些问题可以通过预条件技术得到有效改善。因此基于压缩感知和预条件随机梯度的时间域航空电磁三维反演方法可以保证与传统三维反演方法收敛速度相当，同时保证反演精度的同时极大的提升反演计算效率。
发明内容
本部分的目的在于概述本发明的实施方式的一些方面以及简要介绍一些较佳实施方式。在本部分以及本申请的说明书摘要和发明名称中可能会做些简化或省略以避免使本部分、说明书摘要和发明名称的目的模糊，而这种简化或省略不能用于限制本发明的范围。
鉴于上述和/或现有技术中存在的问题，提出了本发明。
因此，本发明的目的是提供基于压缩感知和预条件随机梯度的航空电磁三维反演方法，收敛速度可以媲美传统全批量数据三维反演，同时在保证反演精度的同时极大的提升反演计算效率。
为解决上述技术问题，根据本发明的一个方面，本发明提供了如下技术方案：
基于压缩感知和预条件随机梯度的航空电磁三维反演方法，包括以下步骤：
步骤一，筛选航空电磁实测数据中的可用数据用于反演；
步骤二，使用泊松圆盘采样法对测点进行随机采样，并对采样点利用有限体积法进行正演计算得到随机测点的预测数据，并使用压缩感知方法对所有测点的预测数据进行重构；
步骤三，将步骤一的实测数据和步骤二的预测数据做数据拟合，并结合模型粗糙度构建航空电磁三维正则化反演的目标函数；
步骤四，对正演方程两侧求导，利用伴随正演计算雅可比矩阵；
步骤五，利用灵敏度矩阵的转置与数据拟合差的乘积，加上模型电导率参数可计算梯度(目标函数一阶导数)，灵敏度矩阵结合模型协方差和数据协方差可得到海森矩阵(目标函数二阶导数)；
步骤六，建立预条件算子，形成预条件随机梯度-高斯-牛顿方法反演方程，计算模型更新量，得到新的迭代模型；
步骤七，重复执行步骤二到步骤六，直到达到最大迭代次数或者符合终止条件，得到最终的反演结果。
作为本发明所述的基于压缩感知和预条件随机梯度的航空电磁三维反演方法的一种优选方案，其中：所述步骤一中，将航空电磁数据进行天电噪声去除、背景场去除、运动噪声去除等预处理后，筛选可用的航空电磁数据用于反演。
作为本发明所述的基于压缩感知和预条件随机梯度的航空电磁三维反演方法的一种优选方案，其中：所述步骤二中，泊松圆盘采样将所有测点的区域划分为若干相同半径且不重叠的圆圈区域，在每个圆圈里只取一个点，完成对所有测点的随机采样，避免传统高斯采样中局部样点过于密集或稀疏的问题。针对随机欠采样点，利用有限体积法进行三维正演计算，电场的双旋度方程可表示为：
其中，Eb和Es分别代表背景场和二次场的电场强度，t表示时间，σ、ε和μ分别表示电导率、介电常数和磁导率，表示哈密顿算子，时间二阶导数项由位移电流引起，此项远小于传导电流，可忽略。
对方程(1)进行规则网格有限体积法离散，随后对每个控制体积进行积分，可得到方程(1)的积分形式：
其中，表示任意控制体积V的表面积，n代表控制体积各面法向。
对方程(2)采用无条件稳定的后推欧拉法进行时间离散，考虑随机采样获得的所有发射源、所有分量和所有时间道，将其整理为统一形式：
其中
其中Nc为航空电磁空间随机采样的测点数量，Tn为程序内部计算使用的时间道数，C和D为与网格尺寸相关的参数，和表示由随机采样的测点得到的系数、电场值和右端项。
对随机测点的正演响应计算结束后，采用压缩感知对所有测点的数据进行重构，重构是一个求解最优化的过程，即求解：
其中ε代表数据噪声，S为采样矩阵获得的测量结果，ψ为稀疏变换矩阵。即通过L1范数最小化，不断寻找一个更为稀疏的再利用稀疏逆变换得到测区所有测点正演结果。
作为本发明所述的基于压缩感知和预条件随机梯度的航空电磁三维反演方法的一种优选方案，其中：所述步骤三中，根据航空电磁观测数据、正演得到的预测电磁数据构建正则化反演的目标函数为：
其中和分别是数据拟合项和模型约束项，m＝(m1,m2,…mM)T是M维模型电导率参数向量，d＝(d1,d2,…,dN)T是N维数据向量，且N＝nst×n,nst和n分别表示航空电磁测点数量和时间道个数，γ为正则化因子。
为数据拟合项，具体形式为：
其中，f(m)表示模型参数m对应正演响应，Cd表示航空电磁数据的协方差矩阵，T为转置算符。
为模型约束项，具体形式为：
其中，m0表示模型的初始参数，Cm是用来约束反演迭代过程中模型参数变化的模型协方差。
作为本发明所述的基于压缩感知和预条件随机梯度的航空电磁三维反演方法的一种优选方案，其中：所述步骤四中，计算出的灵敏度矩阵是缺失的，计算灵敏度矩阵的过程为：
首先对正演方程(3)两侧对模型参数m取微分，经过简单变换，得到：
模型响应对模型参数的偏导数作为灵敏度矩阵J：
其中，Fb和Fs分别代表背景场响应和异常体产生的异常场响应。
在背景场为一维的半空间或层状介质响应时，不随大地任意一点电导率变化而变化，故：
其中，表示随机采样数据对应的正演系数矩阵，是由系数矩阵和二次电场构成的中间参数，L包含了由正演求解的各位置所有时刻电场直接获得接收机处各实测时间道响应的插值过程，相对于传统方法计算的灵敏度矩阵J表现为列缺失。
记中间变量故：
所以通过解伴随正演方程组(13)便可得到V，再将V代入方程(12)便可求得灵敏度矩阵的转置。
作为本发明所述的基于压缩感知和预条件随机梯度的航空电磁三维反演方法的一种优选方案，其中：所述步骤五中，航空电磁正则化反演目标函数对模型参数求一次导数得到梯度：
其中，为中间变量，且
航空电磁正则化反演目标函数对模型参数的二阶导数可得到海森矩阵：
其中海森矩阵的第一项对应矩阵是满阵且非正定，求解比较复杂，故在高斯-牛顿反演方法中，该项可忽略，即海森矩阵表达式为：
作为本发明所述的基于压缩感知和预条件随机梯度的航空电磁三维反演方法的一种优选方案，其中：所述步骤六中，预条件随机梯度的高斯-牛顿方法的反演方程式为：
其中表示随机采样的测点计算出来的随机梯度，P为预条件算子，是一个对角且正定的矩阵，其形式可以表示为：
Pij＝α/β+uij                     (18)
其中β是采样率，α为权重系数，uij为噪声。
最后通过求解反演方程(17)得到模型更新量。
作为本发明所述的基于压缩感知和预条件随机梯度的航空电磁三维反演方法的一种优选方案，其中：所述步骤七中，在迭代的过程中，若没有满足迭代终止条件，则判断相邻两次迭代数据拟合的均方根误差的减小量是否小于减小量阈值，若小于减小量阈值则更新正则化因子，继续迭代计算。
作为本发明所述的基于压缩感知和预条件随机梯度的航空电磁三维反演方法的一种优选方案，其中：所述步骤七中，迭代终止条件包括：(1)数据拟合的均方根误差小于误差阈值；(2)航空电磁正则化反演目标函数收敛；(3)达到预设最大迭代次数，只要满足任何一个迭代终止条件就迭代终止。
与现有技术相比：本发明提供的基于压缩感知和预条件随机梯度的时间域航空电磁快速三维反演方法中，反演采用的数值模拟算法为基于压缩感知的三维时间域有限体积法快速正演算法，能够模拟任意电导率结构下的电磁场信号，其具有较快速度且较高精度模拟复杂地质地形结构并节省内存的优点；同时反演方法基于预条件随机梯度的高斯-牛顿法进行求解，并提出了有效的预条件算子估计方法、正则化约束方法及灵敏度求解方法，能够在保证精度和收敛速度的情况下，极大的提升三维反演计算效率。
附图说明
为了更清楚地说明本发明实施方式的技术方案，下面将结合附图和详细实施方式对本发明进行详细说明，显而易见地，下面描述中的附图仅仅是本发明的一些实施方式，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其它的附图。其中：
图1是本发明实施例提供的基于压缩感知和预条件随机梯度的时间域航空电磁快速三维反演方法的流程图；
图2是本发明实施例提供的均匀半空间模型，其中，(a)为均匀半空间模型的示意图，图中黑点为航空电磁接收点位置，距离地面40m，下面的横线表示地表分界线；(b)为正演VTEM发射的梯形波电流形状；
图3是本发明实施例提供的均匀半空间模型时间域航空电磁三维有限体积法与一维半解析解的精度验证的示意图，其中，(a)为均匀半空间模型时间域航空电磁三维有限体积法与一维半解析解的对比结果；(b)为不同时间道的单点相对误差分析图；
图4是本发明实施例提供的单个异常体的模型示意图，色标为对数分布；
图5是本发明实施例提供的早期时间道重构结构图，(a)、(b)、(c)色标为对数分布，(d)色标为线性分布，其中，(a)为早期时间道的正演响应；(b)为40％采样率随机采样获得的正演响应图；(c)为压缩感知重构获得的响应；(d)为通过压缩感知重构获得的响应与原始响应的单点相对误差；
图6是本发明实施例提供的多个异常体的模型示意图以及多种反演方法的反演结果图，色标为对数分布，其中，(a)多个异常体模型示意图；(b)为全批量数据的高斯-牛顿法反演结果；(c)为基于随机梯度的高斯-牛顿法反演结果；(d)为基于预条件随机梯度的高斯-牛顿法反演结果；
图7是本发明实施例提供的三种反演方法迭代参数图，其中，(a)为三种反演方法反演过程中RMS变化图；(b)为三种反演方法反演过程中目标函数值变化图；
表1是本发明实施例提供的不同反演方法反演迭代参数表，其中，包括三种反演方法终止时的目标函数值、RMS、迭代次数和所用时间，以及相比于全批量数据的高斯-牛顿法提升的效率；
图8采用40％采样率的不同加权因子时基于压缩感知和预条件随机梯度三维反演结果；
表2是基于压缩感知和预条件随机梯度反演与传统全数据三维反演耗时对比；
图9是加权因子为1.4和1.8的基于40％采样率的压缩感知和预条件随机梯度反演算法的不同深度反演结果切片；
图10是该实例中传统全批量数据反演40％采样预条件随机梯度反演方法的数据拟合均方误差衰减曲线。
具体实施方式
为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图对本发明的具体实施方式做详细的说明。
在下面的描述中阐述了很多具体细节以便于充分理解本发明，但是本发明还可以采用其他不同于在此描述的其它方式来实施，本领域技术人员可以在不违背本发明内涵的情况下做类似推广，因此本发明不受下面公开的具体实施方式的限制。
其次，本发明结合示意图进行详细描述，在详述本发明实施方式时，为便于说明，表示器件结构的剖面图会不依一般比例作局部放大，而且所述示意图只是示例，其在此不应限制本发明保护的范围。此外，在实际制作中应包含长度、宽度及深度的三维空间尺寸。
为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明的实施方式作进一步地详细描述。
本发明提供基于压缩感知和预条件随机梯度的航空电磁三维反演方法，收敛速度可以媲美传统全批量数据三维反演，同时在保证反演精度的同时极大的提升反演计算效率，请参阅图1-图7，包括以下步骤：
S1，筛选航空电磁实测数据中的可用数据用于反演；
将航空电磁数据进行天电噪声去除、背景场去除、运动噪声去除等预处理后，筛选可用的航空电磁数据用于反演。
航空电磁数据在采集过程中受到多种因素干扰，信号发生畸变和扭曲，从而对进一步的反演解释造成巨大困难，所以需要对航空电磁数据进行预处理，得到可用的航空电磁数据。
S2，使用泊松圆盘采样法对测点进行随机采样，并对采样点利用有限体积法进行三维正演得到随机测点的预测数据，并使用压缩感知方法对所有测点的预测数据进行重构；
泊松圆盘采样将所有测点的区域划分为若干相同半径且不重叠的圆圈区域，在每个圆圈里只取一个点，完成对所有测点的随机采样。针对随机采样点进行有限体积法正演计算时，电场的双旋度方程可表示为：
其中，Eb和Es分别代表背景场和二次场的电场强度，t表示时间，σ、ε和μ分别表示电导率、介电常数和磁导率，表示哈密顿算子，时间二阶导数项为位移电流项，此项远小于传导电流，可忽略。
图2是实施例提供的均匀半空间模型示意图。其中，(a)为均匀半空间模型的示意图，图中黑点为航空电磁接收点位置，距离地面40m，下面的横线表示地表分界线，设置空气电导率为1×108Ω·m，半空间电导率为100Ω·m；(b)为正演VTEM发射的梯形波电流形状，off-time为7.516ms。
对方程(1)进行规则网格有限体积法离散，随后对每个控制体积进行积分，可得到方程(1)的积分形式：
其中，表示任意控制体积V的表面积，n代表控制体积各面法向。
本实施例中采用梯形波作为发射波形，由于空气和地下介质存在巨大的电性差异，如果选用显式时间离散，需要有足够短的时间步长Δt以满足稳定性条件，这对于晚期道响应计算非常不利。因此，本实施例中采用隐式后推欧拉方法进行时间离散，该方法无论时间步长Δt如何选择均无条件稳定。即将一阶后推欧拉方法应用于方程式(2)，考虑随机采样获得的所有发射源、所有分量和所有时间道，将其整理为统一形式：
其中
其中Nc为航空电磁空间随机采样的测点数量，Tn为程序内部计算使用的时间道数，C和D为与网格尺寸相关的参数，和表示由随机采样的测点得到的系数、电场值和右端项。
对随机测点的正演响应计算结束后，采用压缩感知对所有测点的数据进行重构，重构是一个求解最优化的过程，即求解：
其中ε代表数据噪声，S为采样矩阵获得的测量结果，ψ为稀疏变换矩阵。即通过L1范数最小化，不断寻找一个更为稀疏的最后通过稀疏逆变换获得测区所有测点的正演响应用于计算反演目标函数。
为了检验正演过程和压缩感知重构的精度，首先设计了如图2(a)所示的半空间模型。利用本实验例所使用的有限体积法进行计算，并将计算结果与一维半解析解进行对比。模型参数如下：空气电阻率取1×108Ω·m，半空间电阻率为100Ω·m，发射线圈直径为26m，发射如图2(b)所示的梯形电流，飞行高度为40m，黑点为航空电磁接收点位置。
均匀半空间验证结果如图3所示，(a)为均匀半空间模型时间域航空电磁三维有限体积法与一维半解析解的对比结果；(b)为不同时间道的单点相对误差分析图；由图2可见，二者的单点相对误差除了最后一个时间道都小于5％，说明本实验例采用的正演算法精度较高。
为了检验本实验例采用的压缩重构的效果，我们设计了如图4所示的单异常体模型，经过随机采样后，选取一个早期道对随机采样40％的测点进行正演响应的计算，再通过压缩重构进行全区数据的重构，模型参数如下：异常体电阻率为1Ω·m，上顶面埋深32m，共厚52m。半空间电阻率为100Ω·m，发射线圈直径为26m，发射如图2(b)所示的梯形电流，飞行高度为40m，空气电阻率取1×108Ω·m。
验证结果如图5所示，(a)为早期时间道的正演响应；(b)为40％采样率随机采样获得的正演响应图；(c)为压缩感知重构获得的响应；(d)为通过压缩感知重构获得的响应与原始响应的单点相对误差。由图5可见，随机采样率为40％时可以很好地实现测区所有测点三维正演信号的计算，重构误差都在4％以内，说明本实验例采样率为40％时，采用的压缩重构算法具有很高的精度，可以用于反演研究。
S3，将S1筛选的实测数据和S2的预测数据做数据拟合，并结合模型粗糙度构建航空电磁三维正则化反演的目标函数；
根据航空电磁观测数据、正演得到的预测电磁数据构建正则化反演的目标函数为：
其中和分别是数据拟合项和模型约束项，m＝(m1,m2,…mM)T是M维模型电导率参数向量，d＝(d1,d2,…,dN)T是N维数据向量，且N＝nst×n,nst和n分别表示航空电磁测点数量和时间道个数，γ为正则化因子。
为数据拟合项，具体形式为：
其中，f(m)表示模型参数m对应正演响应，Cd表示航空电磁数据的协方差矩阵，T为转置算符；
为模型约束项，具体形式为：
其中，m0表示模型的初始参数，Cm是用来约束反演迭代过程中模型参数变化的模型协方差。
S4，对正演方程两侧求导，其中预测数据对模型参数的偏导数为灵敏度矩阵(雅可比矩阵)；
本步骤计算出的灵敏度矩阵是缺失的，计算灵敏度矩阵的过程为：
首先对正演方程(3)两侧对模型参数m取微分，经过简单变换，得到：
模型响应对模型参数的偏导数作为灵敏度矩阵J：
其中，Fb和Fs分别代表背景场响应和异常体产生的异常场响应。
在背景场为一维的半空间或层状介质响应时，不随大地任意一点电导率变化而变化，故：
其中，表示随机采样数据对应的正演系数矩阵，是有系数矩阵和二次电场构成的中间参数，L包含了由正演求解的各位置所有时刻电场直接获得接收机处各实测时间道响应的插值过程，相对于传统方法计算的灵敏度矩阵J是按列缺失的，并用零补充。
记中间变量故：
所以通过解伴随正演方程组(13)便可得到V，再将V代入方程(12)
便可求得灵敏度矩阵的转置。
S5，将灵敏度矩阵的转置与预测数据与观测数据的拟合差的乘积，加上模型电导率参数可计算目标函数的梯度，灵敏度矩阵结合模型协方差和数协方差可得到海森矩阵；
航空电磁正则化反演目标函数对模型参数求一次导数得到梯度：
其中，为中间变量，且
航空电磁正则化反演目标函数对模型参数的二阶导数可得到海森矩阵：
其中海森矩阵的第一项对应矩阵是满阵且非正定，求解比较复杂，故在高斯-牛顿反演方法中，该项可忽略，即海森矩阵表达式为：
S6，采用预条件随机梯度的高斯-牛顿方法计算模型更新量，得到新的迭代模型；
预条件随机梯度的高斯-牛顿方法的反演方程式为：
其中表示随机采样的测点计算出来的随机梯度，P为预条件算子，是一个对角且正定的矩阵，其形式可以表示为：
Pij＝α/β+uij,i＝j                  (18)
其中β是采样率，α为权重系数，uij为噪声。
最后通过求解反演方程(17)得到模型更新量。
S7，重复执行S2～S6，直到达到最大迭代次数或者符合终止条件，得到最终的反演参考模型。
在迭代的过程中，若没有满足迭代终止条件，则判断相邻两次迭代的数据拟合均方根误差的减小量是否小于减小量阈值，若小于减小量阈值则更新正则化因子，继续迭代计算。本实施例中，迭代终止条件包括：(1)数据拟合均方根误差RMS小于误差阈值，即RMS＜r0；(2)航空电磁正则化反演目标函数收敛，即||gk||≤ε；(3)达到预设最大迭代次数N，只要满足任何一个迭代终止条件就迭代终止，得到符合地下真实地电结构的反演解。其中，RMS的计算公式为：
其中，εi为数据误差，n为测点个数。
在迭代的过程中，若没有满足迭代终止条件，则判断相邻两次迭代的数据拟合均方根误差RMS的减小量是否小于减小量阈值χ，χ＞0，一般取2％，若小于减小量阈值则更新正则化因子，继续迭代计算。
本实施例提供了一个反演算例，以验证反演算法的正确性与高效性。图6是本发明实施例提供的多个异常体模型示意图及反演结果图。其中，(a)多个异常体模型示意图；(b)为全批量数据的高斯-牛顿法反演结果；(c)为基于随机梯度的高斯-牛顿法反演结果；(d)为基于预条件随机梯度的高斯-牛顿法反演结果。如图6(a)所示，为多个异常体分布不同位置，形状各不相同。异常体电阻率为1Ω·m，上顶面埋深80m，共厚120m。半空间电阻率为100Ω·m，发射线圈直径为26m，发射如图2(b)所示的梯形电流，飞行高度为40m，空气电阻率取1×108Ω·m。由图6(b)、(c)、(d)可见，三个反演方法都对目标异常体的形态、位置和电阻率值都有很好的反映。
三种反演方法的RMS和目标函数值随迭代次数变化情况如图7所示，(a)为三种反演方法反演过程中RMS变化图；(b)为三种反演方法反演过程中目标函数值变化图。同时三种反演方法参数对比如表1所示。
由图7可知，基于预条件随机梯度的高斯-牛顿法的收敛速度可以媲美全批量数据的高斯-牛顿法的收敛速度。由表1可知，全批量数据的高斯-牛顿法反演中，迭代4次时RMS下降至1.78，总耗时为70小时。基于随机梯度的高斯-牛顿法方法需要6次迭代，RMS下降至1.82，但耗时36小时，比传统方法效率有一定提升；在加入了预条件后，仅需4次迭代，RMS便下降至1.92，且耗时20h，这说明迭代次数和收敛时间均有明显改善，且迭代次数与全批量数据的高斯-牛顿法反演迭代次数相当，耗时约为20小时，与全批量数据的高斯-牛顿法反演70小时相比，效率提升约71.4％，三维反演效率得到明显改善。
实例二为基于以上实施过程开展的地下单一良导异常体模型数据反演实验。图8是采用40％采样率的不同加权因子时基于压缩感知和预条件随机梯度三维反演结果，其中黑色线框为真实模型位置。图8可知，加权因子为0.6、1.0以及1.6时预条件反演算法均能很好地恢复地下模型位置和电导率，验证了方法的有效性。表2是基于压缩感知和预条件随机梯度反演与传统全数据三维反演耗时对比，可以看出，全批量高斯牛顿反演需要55.6小时，而加权因子为1.6时仅耗时12.3小时，效率提高78％。
实例三为基于以上实施过程开展的复杂地电结构航空电磁三维反演结果。图9是加权因子为1.4和1.8的基于40％采样率的压缩感知和预条件随机梯度反演算法的不同深度反演结果切片。图中黑色线框为实际目标体位置。可以看出，在深度为50m和80m时，反演结果能够很好低恢复浅部异常，在深度为150m时，由于异常尺度较小、埋深较大导致异常信号弱，反演结果对部分大尺度异常体能够实现较好地刻画。整体反演能够体现异常的主要展布、位置和电导率。图10是该实例中传统全批量数据反演40％采样预条件随机梯度反演方法的数据拟合均方误差衰减曲线，可以看出40％的采样计算能够实现与传统全数据反演相媲美的收敛速度，但每次迭代的计算量仅是传统方法的40％左右，大幅度减少三维反演计算量，提高了反演速度。该实例中全数据反演耗时252小时，预条件随机梯度反演耗时76小时，整体效率提高约70％。
虽然在上文中已经参考实施方式对本发明进行了描述，然而在不脱离本发明的范围的情况下，可以对其进行各种改进并且可以用等效物替换其中的部件。尤其是，只要不存在结构冲突，本发明所披露的实施方式中的各项特征均可通过任意方式相互结合起来使用，在本说明书中未对这些组合的情况进行穷举性的描述仅仅是出于省略篇幅和节约资源的考虑。因此，本发明并不局限于文中公开的特定实施方式，而是包括落入权利要求的范围内的所有技术方案。

Claims (9)

1.基于压缩感知和预条件随机梯度的航空电磁三维反演方法，其特征在于，包括以下步骤：
步骤一，筛选航空电磁实测数据中的可用数据用于反演；
步骤二，使用泊松圆盘采样法对测点进行随机采样，并对采样点利用有限体积法进行正演计算得到随机测点的预测数据，并使用压缩感知方法对所有测点的预测数据进行重构；
步骤三，将步骤一的实测数据和步骤二的预测数据做数据拟合，并结合模型粗糙度构建航空电磁三维正则化反演的目标函数；
步骤四，对正演方程两侧求导，利用伴随正演计算雅可比矩阵；
步骤五，利用灵敏度矩阵的转置与数据拟合差的乘积，加上模型电导率参数可计算梯度(目标函数一阶导数)，灵敏度矩阵结合模型协方差和数据协方差可得到海森矩阵(目标函数二阶导数)；
步骤六，建立预条件算子，形成预条件随机梯度-高斯-牛顿方法反演方程，计算模型更新量，得到新的迭代模型；
步骤七，重复执行步骤二到步骤六，直到达到最大迭代次数或者符合终止条件，得到最终的反演结果。
2.根据权利要求1所述的基于压缩感知和预条件随机梯度的航空电磁三维反演方法，其特征在于，所述步骤一中，将航空电磁数据进行天电噪声去除、背景场去除、运动噪声去除等预处理后，筛选可用的航空电磁数据用于反演。
3.根据权利要求1所述的基于压缩感知和预条件随机梯度的航空电磁三维反演方法，其特征在于，所述步骤二中，泊松圆盘采样将所有测点的区域划分为若干相同半径且不重叠的圆圈区域，在每个圆圈里只取一个点，完成对所有测点的随机采样，避免传统高斯采样中局部样点过于密集或稀疏的问题。针对随机欠采样点，利用有限体积法进行三维正演计算，电场的双旋度方程可表示为：
Figure FDA0003939671050000021
其中，Eb和Es分别代表背景场和二次场的电场强度，t表示时间，σ、ε和μ分别表示电导率、介电常数和磁导率，
Figure FDA0003939671050000022
表示哈密顿算子，时间二阶导数项由位移电流引起，此项远小于传导电流，可忽略。
对方程(1)进行规则网格有限体积法离散，随后对每个控制体积进行积分，可得到方程(1)的积分形式：
Figure FDA0003939671050000023
其中，
Figure FDA0003939671050000024
表示任意控制体积V的表面积，n代表控制体积各面法向。
对方程(2)采用无条件稳定的后推欧拉法进行时间离散，考虑随机采样获得的所有发射源、所有分量和所有时间道，将其整理为统一形式：
Figure FDA0003939671050000025
其中
Figure FDA0003939671050000026
其中Nc为航空电磁空间随机采样的测点数量，Tn为程序内部计算使用的时间道数，C和D为与网格尺寸相关的参数，
Figure FDA0003939671050000027
和
Figure FDA0003939671050000028
表示由随机采样的测点得到的系数、电场值和右端项。
对随机测点的正演响应计算结束后，采用压缩感知对所有测点的数据进行重构，重构是一个求解最优化的过程，即求解：
Figure FDA0003939671050000029
其中ε代表数据噪声，S为采样矩阵获得的测量结果，ψ为稀疏变换矩阵。即通过L1范数最小化，不断寻找一个更为稀疏的x～，再利用稀疏逆变换得到测区所有测点正演结果。
4.根据权利要求1所述的基于压缩感知和预条件随机梯度的航空电磁三维反演方法，其特征在于，所述步骤三中，根据航空电磁观测数据、正演得到的预测电磁数据构建正则化反演的目标函数为：
Figure FDA0003939671050000031
其中
Figure FDA0003939671050000032
和
Figure FDA0003939671050000033
分别是数据拟合项和模型约束项，m＝(m1,m2,…mM)T是M维模型电导率参数向量，d＝(d1,d2,…,dN)T是N维数据向量，且N＝nst×n,nst和n分别表示航空电磁测点数量和时间道个数，γ为正则化因子。
Figure FDA0003939671050000034
为数据拟合项，具体形式为：
Figure FDA0003939671050000035
其中，f(m)表示模型参数m对应的正演响应，Cd表示航空电磁数据的协方差矩阵，T为转置算符。
Figure FDA0003939671050000036
为模型约束项，具体形式为：
Figure FDA0003939671050000037
其中，m0表示模型的初始参数，Cm是用来约束反演迭代过程中模型参数变化的模型协方差。
5.根据权利要求1所述的基于压缩感知和预条件随机梯度的航空电磁三维反演方法，其特征在于，所述步骤四中，计算出的灵敏度矩阵是缺失的，计算灵敏度矩阵的过程为：
首先对正演方程(3)两侧对模型参数m取微分，经过简单变换，得到：
Figure FDA0003939671050000041
模型响应对模型参数的偏导数作为灵敏度矩阵J：
Figure FDA0003939671050000042
其中，Fb和Fs分别代表背景场响应和异常体产生的异常场响应。
在背景场为一维的半空间或层状介质响应时，不随大地任意一点电导率变化而变化，故：
Figure FDA0003939671050000043
Figure FDA0003939671050000044
其中，
Figure FDA0003939671050000045
表示随机采样数据对应的正演系数矩阵，
Figure FDA0003939671050000046
是由系数矩阵和二次电场构成的中间参数，L包含了由正演求解的各位置所有时刻电场直接获得接收机处各实测时间道响应的插值过程，
Figure FDA0003939671050000047
相对于传统方法计算的灵敏度矩阵J表现列缺失。
记中间变量
Figure FDA0003939671050000048
故：
Figure FDA0003939671050000049
Figure FDA00039396710500000410
所以通过解伴随正演方程组(13)便可得到V，再将V代入方程(12)便可求得灵敏度矩阵的转置。
6.根据权利要求1所述的基于压缩感知和预条件随机梯度的航空电磁三维反演方法，其特征在于，所述步骤五中，航空电磁正则化反演目标函数对模型参数求一次导数得到梯度：
Figure FDA00039396710500000411
其中，
Figure FDA0003939671050000051
为中间变量，且
Figure FDA0003939671050000052
航空电磁正则化反演目标函数对模型参数的二阶导数可得到海森矩阵：
Figure FDA0003939671050000053
其中海森矩阵的第一项对应矩阵是满阵且非正定，求解比较复杂，故在高斯-牛顿反演方法中，该项可忽略，即海森矩阵表达式为：
Figure FDA0003939671050000054
7.根据权利要求1所述的基于压缩感知和预条件随机梯度的航空电磁三维反演方法，其特征在于，所述步骤六中，预条件随机梯度的高斯-牛顿方法的反演方程式为：
Figure FDA0003939671050000055
其中
Figure FDA0003939671050000056
表示随机采样的测点计算出来的随机梯度，P为预条件算子，是一个对角且正定的矩阵，其形式可以表示为：
Pij＝α/β+uij                     (18)
其中β是采样率，α为权重系数，uij为噪声。
最后通过求解反演方程(17)得到模型更新量。
8.根据权利要求1所述的基于压缩感知和预条件随机梯度的航空电磁三维反演方法，其特征在于，所述步骤七中，在迭代的过程中，若没有满足迭代终止条件，则判断相邻两次迭代数据拟合的均方根误差的减小量是否小于减小量阈值，若小于减小量阈值则更新正则化因子，继续迭代计算。
9.根据权利要求1所述的基于压缩感知和预条件随机梯度的航空电磁三维反演方法，其特征在于，所述步骤七中，迭代终止条件包括：(1)数据拟合的均方根误差小于误差阈值；(2)航空电磁正则化反演目标函数收敛；(3)达到预设最大迭代次数，只要满足任何一个迭代终止条件就迭代终止。
CN202211414788.1A 2022-11-11 2022-11-11 基于压缩感知和预条件随机梯度的航空电磁三维反演方法 Pending CN116090283A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211414788.1A CN116090283A (zh) 2022-11-11 2022-11-11 基于压缩感知和预条件随机梯度的航空电磁三维反演方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211414788.1A CN116090283A (zh) 2022-11-11 2022-11-11 基于压缩感知和预条件随机梯度的航空电磁三维反演方法

Publications (1)

Publication Number Publication Date
CN116090283A true CN116090283A (zh) 2023-05-09

Family

ID=86203265

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211414788.1A Pending CN116090283A (zh) 2022-11-11 2022-11-11 基于压缩感知和预条件随机梯度的航空电磁三维反演方法

Country Status (1)

Country Link
CN (1) CN116090283A (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116794735A (zh) * 2023-06-02 2023-09-22 中国自然资源航空物探遥感中心 航空磁矢量梯度数据等效源多分量联合去噪方法及装置
CN117725354A (zh) * 2024-02-18 2024-03-19 中国地质大学（北京） 一种大数据量重力与重力梯度联合的快速正反演方法及系统

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116794735A (zh) * 2023-06-02 2023-09-22 中国自然资源航空物探遥感中心 航空磁矢量梯度数据等效源多分量联合去噪方法及装置
CN116794735B (zh) * 2023-06-02 2024-03-19 中国自然资源航空物探遥感中心 航空磁矢量梯度数据等效源多分量联合去噪方法及装置
CN117725354A (zh) * 2024-02-18 2024-03-19 中国地质大学（北京） 一种大数据量重力与重力梯度联合的快速正反演方法及系统
CN117725354B (zh) * 2024-02-18 2024-04-26 中国地质大学（北京） 一种大数据量重力与重力梯度联合的快速正反演方法及系统

Similar Documents

Publication Publication Date Title
CN116090283A (zh) 基于压缩感知和预条件随机梯度的航空电磁三维反演方法
CN111551992B (zh) 岩石储层构造表征方法、装置、计算机可读存储介质及电子设备
EP2567063B1 (en) Artifact reduction in method of iterative inversion of geophysical data
Li et al. Fast inversion of large-scale magnetic data using wavelet transforms and a logarithmic barrier method
CN110568486B (zh) 基于同步稀疏低秩张量补全模型的地震信号补全方法
CN107015274A (zh) 一种缺失地震勘探数据恢复重构方法
KR20080023946A (ko) 오일러 해를 이용한 지하공동의 3차원 중력역산 방법 및이를 이용한 3차원 영상화 방법
CN110852025B (zh) 一种基于超收敛插值逼近的三维电磁慢扩散数值模拟方法
CN106483559A (zh) 一种地下速度模型的构建方法
CN115047530A (zh) 基于一维反演的地空瞬变电磁数据三维频率域解释方法
CN114970290B (zh) 一种地面瞬变电磁法并行反演方法及系统
CN111025385B (zh) 一种基于低秩和稀疏约束的地震数据重建方法
CN110119586B (zh) 轴向电导率各向异性瞬变电磁三分量三维fdtd正演方法
Wu et al. Adaptive feedback convolutional‐neural‐network‐based high‐resolution reflection‐waveform inversion
CN112444850B (zh) 地震资料速度建模方法、存储介质及计算设备
CN113267287B (zh) 冲击波超压三维时空场重建方法
CN116165722A (zh) 一种采用高斯牛顿法的回线源瞬变电磁三维快速反演方法
CN112526621B (zh) 一种基于神经网络的地空电磁数据慢扩散多参数提取方法
Zhong et al. Electrical resistivity tomography with smooth sparse regularization
CN112946760B (zh) 基于正则化方法的未爆弹三维立体成像方法、装置及系统
CN113267830A (zh) 基于非结构网格下二维重力梯度与地震数据联合反演方法
CN113970732A (zh) 一种三维频率域探地雷达双参数同步反演方法
CN113325482A (zh) 一种时间域电磁数据反演成像方法
CN110794469B (zh) 基于最小地质特征单元约束的重力反演方法
CN112379415A (zh) 基于减采样的重构低频数据多尺度全波形反演方法及装置

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