CN115563791A - 基于压缩感知重构的大地电磁数据反演方法 - Google Patents
基于压缩感知重构的大地电磁数据反演方法 Download PDFInfo
- Publication number
- CN115563791A CN115563791A CN202211258637.1A CN202211258637A CN115563791A CN 115563791 A CN115563791 A CN 115563791A CN 202211258637 A CN202211258637 A CN 202211258637A CN 115563791 A CN115563791 A CN 115563791A
- Authority
- CN
- China
- Prior art keywords
- data
- electromagnetic
- magnetotelluric
- inversion
- model
- 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.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 40
- 238000004364 calculation method Methods 0.000 claims abstract description 12
- 238000012545 processing Methods 0.000 claims abstract description 6
- 239000011159 matrix material Substances 0.000 claims description 36
- 238000005070 sampling Methods 0.000 claims description 24
- 239000013598 vector Substances 0.000 claims description 20
- 238000004422 calculation algorithm Methods 0.000 claims description 16
- 238000005259 measurement Methods 0.000 claims description 15
- 238000005457 optimization Methods 0.000 claims description 13
- 238000002939 conjugate gradient method Methods 0.000 claims description 10
- 230000009466 transformation Effects 0.000 claims description 6
- 238000009795 derivation Methods 0.000 claims description 4
- OAICVXFJPJFONN-UHFFFAOYSA-N Phosphorus Chemical compound [P] OAICVXFJPJFONN-UHFFFAOYSA-N 0.000 claims description 3
- 230000017105 transposition Effects 0.000 claims description 2
- 238000012886 linear function Methods 0.000 claims 1
- 230000006870 function Effects 0.000 description 24
- 238000010586 diagram Methods 0.000 description 15
- 238000010276 construction Methods 0.000 description 5
- 238000001514 detection method Methods 0.000 description 3
- 238000005516 engineering process Methods 0.000 description 3
- 230000005684 electric field Effects 0.000 description 2
- 230000005611 electricity Effects 0.000 description 2
- 229910052500 inorganic mineral Inorganic materials 0.000 description 2
- 239000011707 mineral Substances 0.000 description 2
- 230000008569 process Effects 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 230000002159 abnormal effect Effects 0.000 description 1
- 238000007792 addition Methods 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000005672 electromagnetic field Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000011084 recovery Methods 0.000 description 1
- 230000035945 sensitivity Effects 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 239000007787 solid Substances 0.000 description 1
- 230000001131 transforming effect Effects 0.000 description 1
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
-
- H—ELECTRICITY
- H03—ELECTRONIC CIRCUITRY
- H03M—CODING; DECODING; CODE CONVERSION IN GENERAL
- H03M7/00—Conversion of a code where information is represented by a given sequence or number of digits to a code where the same, similar or subset of information is represented by a different sequence or number of digits
- H03M7/30—Compression; Expansion; Suppression of unnecessary data, e.g. redundancy reduction
- H03M7/3059—Digital compression and data reduction techniques where the original information is represented by a subset or similar information, e.g. lossy compression
- H03M7/3062—Compressive sampling or sensing
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/04—Constraint-based CAD
Landscapes
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Physics & Mathematics (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明公开了一种基于压缩感知重构的大地电磁数据反演方法，包括以下步骤：获取随机测点分布采集的大地电磁数据，并将大地电磁数据处理成电磁阻抗数据后，对电磁阻抗数据进行分类整理；采用基于曲波变换的压缩感知重构方式对整理后的每类电磁阻抗数据进行从随机测点分布到规则测点分布的数据重构，得到重构电磁阻抗数据；以重构电磁阻抗数据作为观测电磁数据进行大地电磁数据反演计算。该大地电磁数据反演方法提升了大地电磁勘探的工作效率，并保证反演的收敛速度和结果的分辨率。
Description
技术领域
本发明属于地质探测领域，具体涉及一种基于压缩感知重构的大地电磁数据反演方法。
背景技术
鉴于全球近地表易勘探的油气和矿产资源开采有限，迫切需要在地球深部寻找新的突破。大地电磁测深是利用天然交变电磁场研究地下电性结构的重要方法，具有成本低、施工方便、探测深度大和对低阻体灵敏等技术优势，已被广泛应用于油气资源勘探和深部构造研究。
传统的一维、二维反演策略已无法满足深部资源探测和构造研究的技术需求。近年发展起来并获得广泛应用的大地电磁三维正反演算法是研究地球内部结构必不可少的技术手段，而处理技术的有效性和处理结果的可靠性直接影响地质解释人员的分析和判断。
电磁法三维正演模拟是三维反演和解释的基础和核心内容，三维正演计算的精度关系到三维反演结果的可靠性并影响最后的地质解释。目前在大地电磁中应用广泛的正演算法是基于交错网格的有限差分法，其具有思路简单，易于编程实现，计算速度快的优点。
目前三维大地电磁数据反演的主要方法为基于梯度类优化算法的确定性反演，在模型解空间搜索最优的下降方向，经有限次迭代就可以完成收敛，从而达到最优解。各种梯度类优化算法包括高斯-牛顿(Gauss-Newton，GN)法、拟牛顿(Quasi-Newton，QN)法、非线性共轭梯度(Nonlinear Conjugate Gradient，NLCG)法和有限内存的BFGS算法(LMBFGS，一种拟牛顿方法)等，这些梯度类优化算法被引用到寻找三维电阻率模型最优解问题。
通常情况下，大地电磁野外勘探施工成本高，且对于工区复杂的区域开展野外工作是十分困难的。因此，大地电磁实际观测到的数据信息量有限，但反演的未知模型参数远大于采集的数据量，因此反演方程是欠定的，导致反演问题具有非常强的多解性，这会导致反演结果的分辨率不足，进而导致反演结果不够准确。
因此，开发能够有效提高电磁数据反演分辨率的三维正反演算法是目前大地电磁领域的关键任务。
发明内容
鉴于上述，本发明的目的是提供一种基于压缩感知重构的大地电磁数据反演方法，提高大地电磁勘探的工作效率，并保证反演收敛速度和反演结果分辨率。
为实现上述发明目的，本发明实施例提供的一种基于压缩感知重构的大地电磁数据反演方法，包括以下步骤：
获取随机测点分布采集的大地电磁数据，并将大地电磁数据处理成电磁阻抗数据后，对电磁阻抗数据进行分类整理；
采用基于曲波变换的压缩感知重构方式对整理后的每类电磁阻抗数据进行从随机测点分布到规则测点分布的数据重构，得到重构电磁阻抗数据；
以重构电磁阻抗数据作为观测电磁数据进行大地电磁数据反演计算。
优选地，所述对电磁阻抗数据进行分类整理，包括：
按照采样频率、阻抗类型以及电磁阻抗数据的实部和虚部进行分类，将每个采样频率下的每类电磁阻抗数据的实部或虚部作为一类数据，当采样频率为P个、阻抗类型为Q类，那整理得到P×Q×2类电磁阻抗数据，其中，2表示电磁阻抗数据的实部或虚部这两种。
优选地，所述采用基于曲波变换的压缩感知重构算法对整理后的每类电磁阻抗数据进行从随机测点分布到规则测点分布的数据重构，包括：
对随机测点分布进行稀疏采样，得到稀疏的采样矩阵S；
采用曲波变换将重构电磁阻抗数据f以正交基Ψ的形式展开，即：
f＝Ψx
其中，x是重构电磁阻抗数据f的系数序列，称为变换域系数向量，正交基Ψ＝[Ψ1,Ψ2,…,ΨN]表示N个正交基向量，N为规则测点的个数；
考虑重构电磁阻抗数据f∈RN，采样矩阵S∈RM×N，M<<N，M为随机分布测点的个数，R为实数域，则得到线性测量值
则有：
这样根据线性测量值
和已知观测的每类电磁阻抗数据y构建重构电磁阻抗数据求解的最优化问题，表示为：
其中，ε代表数据噪声，||·||1表示L1范数，||·||2表示L2范数，通过对最优化问题求解，得到变换域系数向量x，然后，根据得到的变换域系数向量x以及f＝Ψx经过稀疏逆变换得到每类电磁阻抗数据对应的重构电磁阻抗数据f。
优选地，所述以重构电磁阻抗数据作为观测电磁数据进行大地电磁数据反演计算，包括：
(a)以重构电磁阻抗数据作为观测电磁数据，计算观测电磁数据与对地下电阻率模型进行正演的预测电磁数据之间的拟合差作为数据拟合项，并构建模型粗糙度约束项，根据数据拟合项和模型约束项构建大地电磁三维正则化反演的目标函数；
(b)对目标函数进行相对于模型参数的求导以得到目标函数梯度；
(c)根据目标函数梯度并采用非线性共轭梯度法计算模型参数更新量，进而得到新地下电阻率模型；
(d)重复执行步骤(a)-步骤(c)进行迭代反演，直到收敛条件为止，输出反演得到的最终地下电阻率模型。
优选地，采用交错网格有限差分算法对地下电阻率模型进行正演得到预测电磁数据。
优选地，构建的大地电磁三维正则化反演的目标函数Φ表示为：
其中，dobs为观测电磁数据向量，dprd为由地下电阻率模型正演计算得到的预测电磁数据向量，Wd是数据协方差矩阵，为对角阵，其对角元素是观测电磁数据中噪声标准差的倒数，m为地下电阻率模型，mref为包含模型先验信息的参考模型，Wm为模型协方差矩阵，采用差分算子或拉普拉斯算子，
表示L2范数的平方，λ表示调节权重的正则化参数。
优选地，对目标函数进行相对于模型参数的求导以得到的目标函数梯度g表示为：
其中，JT表示三维大地电磁的雅克比矩阵的转置，通过对地下电阻率模型伴随正演方法推导三维大地电磁的雅克比矩阵。
优选地，根据目标函数梯度并采用非线性共轭梯度法推导模型参数更新量的方程为：
其中，Δm表示模型参数更新量，基于该方程计算模型参数更新量。
优选地，迭代反演的收敛条件为：观测电磁数据与预测电磁数据的均方根误差小于1。
与现有技术相比，本发明具有的有益效果至少包括：
基于曲波变换的压缩感知重构方式实现随机测点分布的每类电磁阻抗数据到完整规则测点分布的大地电磁阻抗数据的高精度重构在，在此基础上，以重构电磁阻抗数据作为观测电磁数据进行大地电磁数据反演计算，这样可以显著提升大地电磁勘探的工作效率，并在一定程度上加速反演收敛和得到与规则采样的数据反演相似的分辨率。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图做简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动前提下，还可以根据这些附图获得其他附图。
图1是本发明实施例提供的基于压缩感知重构的大地电磁数据反演方法的流程图；
图2是本发明实施例提供的60％随机采样的测点与40％随机采样的测点分布示意图；
图3是本发明实施例提供的规则完整测点数据，60％随机采样测点的数据重构与40％随机采样测点的数据重构示意图；
图4是本发明实施例提供的地下三维电阻率真实模型示意图；
图5是本发明实施例提供的大地电磁规则完整测点反演结果，60％随机采样测点的数据重构反演结果和40％随机采样测点的数据重构反演结果的三维异常体对比图；
图6是本发明实施例提供的大地电磁规则完整测点反演结果，60％随机采样测点的数据重构反演结果和40％随机采样测点的数据重构反演结果的二维切片对比图；
图7是本发明实施例提供的大地电磁三维反演参数随迭代变化曲线图。
具体实施方式
为使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例对本发明进行进一步的详细说明。应当理解，此处所描述的具体实施方式仅仅用以解释本发明，并不限定本发明的保护范围。
为降低反演的非唯一性，提高反演结果的分辨率，实施例提供了提出一种基于压缩感知重构的大地电磁数据反演方法，通过稀疏变换将采样测点的的观测电磁数据进行稀疏表达，再利用压缩感知重构方式恢复出整个测区完整规则测点的重构电磁阻抗数据，可给出较高精度的测点数据重构电磁阻抗数据，进而提高反演结果的分辨率。
图1是本发明实施例提供的基于压缩感知重构的大地电磁数据反演方法的流程图。如图1所示，实施例提供的基于压缩感知重构的大地电磁数据反演方法，包括以下步骤：
步骤1，获取随机测点分布采集的大地电磁数据，并将大地电磁数据处理成电磁阻抗数据后，对电磁阻抗数据进行分类整理。
实施例中，在获得随机测点分布采集的大地电磁数据后，对大地电磁数据进行处理以得到电磁阻抗数据，然后为了实现电磁阻抗数据的压缩感知重构，还需要对电磁阻抗数据进行分类整理。具体分类整理包括：按照采样频率、阻抗类型以及电磁阻抗数据的实部和虚部进行分类，将每个采样频率下的每类电磁阻抗数据的实部或虚部作为一类数据，当采样频率为P个、阻抗类型为Q类，那整理得到P×Q×2类电磁阻抗数据，其中，2表示电磁阻抗数据的实部或虚部这两种。
步骤2，采用基于曲波变换的压缩感知重构方式对整理后的每类电磁阻抗数据进行从随机测点分布到规则测点分布的数据重构，得到重构电磁阻抗数据。
实施例中，需要采用压缩感知方法实现数据重构，其基本思想是利用信号的稀疏性质，由M个原始电磁阻抗数据重构获得N个(M<<N)重构电磁阻抗数据。压缩感知实现的两个基本条件为：(1)重构数据是稀疏的或者是可压缩的；(2)采样矩阵S是随机的，和重构数据本身是互不相干。
根据上述压缩感知实现的两个基本条件，压缩感知重构的一个关键点在于采样矩阵S的构建，因此，在重构之前，选取随机测点分布的大地电磁数据，并根据随机分布的测点构建出稀疏的采样矩阵。
图2是实施例提供的60％随机测点分布与40％随机测点分布示意图，假设完整规则的测点分布是29×29，图中实心点代表具有数据的测点位置，图2中(a)为60％随机测点的示意图，number表示观测点个数，图2中(b)为40％随机测点的示意图。
在数据重构时，采用曲波变换将重构电磁阻抗数据f以正交基Ψ的形式展开，即：
f＝Ψx (1)
其中，x是重构电磁阻抗数据f的系数序列，称为变换域系数向量，正交基Ψ＝[Ψ1,Ψ2,…,ΨN]表示N个正交基向量，N为规则测点的个数；
考虑重构电磁阻抗数据f∈RN，采样矩阵S∈RM×N，M<<N，M为随机分布测点的个数，R为实数域，则得到线性测量值
则有：
利用极不完整的观测电磁数据y恢复出完整重构电磁数据f是一个十分欠定的问题。根据压缩感知理论的两个前提条件，可以采用与S不相干的稀疏基Ψ对信号f进行展开，因此可以得到：
这样根据压缩感知理论，结合线性测量值
和已知观测的每类电磁阻抗数据y构建重构电磁阻抗数据求解的最优化问题，表示为：
其中，ε代表数据噪声，||·||1表示L1范数，||·||2表示L2范数，根据稀疏性和不相干性理论，通过对最优化问题求解，利用少量电磁阻抗数据y可得到变换域系数向量x，具体可以采用凸集投影算法求解最优化问题。
然后，根据得到的变换域系数向量x以及f＝Ψx经过稀疏逆变换得到每类电磁阻抗数据对应的重构电磁阻抗数据f。
图3是本发明实施例提供的规则完整测点数据与60％和40％随机测点的重构数据对比示意图，选取的大地电磁数据为0.1Hz频率下的ZXY阻抗数据的实部。图3中(a)为完整规则测点数据图，m表示距离；图3中(b)为60％随机采样测点数据重构示意图；图3中(c)为40％随机采样测点数据重构示意图。由图3可以看到，基于曲波变换的压缩感知数据重构结果整体的形态均具有较好的恢复，这是由于曲波变换对于二维图像具有更加优越地稀疏表示，可以更好地提取数据异常的特征。
步骤3，以重构电磁阻抗数据作为观测电磁数据，计算观测电磁数据与对地下电阻率模型进行正演的预测电磁数据之间的拟合差作为数据拟合项，并构建模型约束项，根据数据拟合项和模型约束项构建大地电磁三维正则化反演的目标函数。
实施例中，为了探测地质结构，刻画用于反演的地下电阻率模型，采用规则网格单元对地下电阻率模型进行几何剖分。同时设置几何模型信息，包括网格单元个数和尺寸信息，以及每个单元的电阻率信息等，并得到初始地下电阻率模型。
实施例中，可以采用规则六面体将计算区域进行剖分。设置几何模型信息包括电导率信息，传递数据信息包括地面测点的个数、位置以及测点频率等。控制参数包括最大迭代次数和正则化因子λ等。
图4是本发明实施例提供的地下电阻率三维理论算例模型示意图。构建均匀半空间模型作为初始电阻率模型，设置初始电阻率为100ohm-m。
实施例中，基于L2范数构建目标函数中的数据拟合项与模型粗糙度约束项，该目标函数用于正则化反演过程来优化模型参数。通过观测电磁数据与地下电阻率模型正演计算得到的预测电磁数据，并采用数据协方差矩阵Wd协助构建数据拟合项。采用差分算子定义模型协方差矩阵Wm，实现相邻单元模型参数之间的约束，进而构建模型粗糙度约束项，因此，得到的大地电磁三维正则化反演的目标函数Φ表示为：
其中，dobs为观测电磁数据向量，dprd为由地下电阻率模型正演计算得到的预测电磁数据向量，Wd是数据协方差矩阵，为对角阵，其对角元素是观测电磁数据中噪声标准差的倒数，m为地下电阻率模型，mref为包含模型先验信息的参考模型，Wm为模型协方差矩阵，
表示L2范数的平方，λ表示调节权重。
实施例中，采用交错网格有限差分算法对地下电阻率模型进行正演得到预测电磁数据。
步骤4，对目标函数进行相对于模型参数的求导以得到目标函数梯度。
实施例中，在对正则化反演的目标函数进行求解时，需要计算雅克比矩阵J和其转置JT。采用伴随正演方法推导三维大地电磁雅克比矩阵的表达式，实现雅克比矩阵的计算。
首先，根据正演数值模拟的线性方程组，将方程两边分别对模型参数m求导可得，
对上式做简单变换得到
其中，K是正演方程的系数矩阵，S为源矢量和边界条件，根据磁场Hs与电场Es的关系，定义磁场Hs可以由空间插值算子L和空间网格上的电场Es获得，则雅克比矩阵可定义为：
假设矩阵G为Nd×Nm维的矩阵，具体形式为：
则公式(8)可重写为：
进一步，可推导出雅克比矩阵的转置为：
定义矩阵：
w＝K-1LT, (12)
则有：
Kw＝LT. (13)
通过伴随正演方式求解(13)式得到矩阵w，进而通过(11)式即可得到雅克比矩阵的转置JT。因此，通过一次正演和一次伴随正演就可计算出目标函数的梯度g，表示为：
步骤5，根据目标函数梯度并采用非线性共轭梯度法计算模型参数更新量，进而得到新地下电阻率模型。
非线性共轭梯度法等反演方法根据梯度信息可以得到反演的下降方向和迭代步长，实施例采用非线性共轭梯度法求解反演问题，该方法对计算机存储要求较低，实现简单，十分适合求解大规模数据的优化问题。具体地，根据目标函数梯度并采用非线性共轭梯度法推导模型参数更新量的方程为：
其中，Δm表示模型参数更新量，基于该方程计算模型参数更新量。
实施例中，依据该预测电磁数据和正则化反演的目标函数实现大地电磁的迭代反演时，通过对正则化反演的目标函数进行求解以获得预测电磁数据和观测电磁数据的拟合差，将拟合差小于设置阈值作为反演迭代约束条件，在不满足约束条件时，依据求解正则化反演的目标函数获得的模型解向量更新方案参考模型。
步骤6，重复执行步骤3-步骤5进行迭代反演，直到收敛条件为止，输出反演得到的最终地下电阻率模型。
实施例中，依据预测电磁数据和正则化反演的目标函数实现大地电磁的迭代反演，采用非线性共轭梯度方法，不断更新模型参数，重复步骤4～步骤6直到满足反演终止条件或达到最大迭代次数，获取符合探测地质结构的最终反演参考模型。
实施例中，迭代最终条件为观测电磁数据与预测电磁数据的均方根误差(RootMean Square，RMS)小于1。其中，其中RMS的计算公式为
其中，N是观测电磁数据个数，dpre和dobs是预测电磁数据和观测电磁数据，δi表示观测电磁数据中噪声标准差，表示观测电磁数据的索引。
上述实施例提供的大地电磁数据反演方法中，通过基于曲波变换的压缩感知技术实现随机测点的数据到完整规则测点分布的大地电磁数据的高精度重构；对初始地下电阻率模型采用基于交错网格的有限差分正演计算得到预测电磁数据；并根据预测电磁数据和重构数据构建数据拟合项，进而结合模型的粗糙度项构建正则化反演的目标函数；采用伴随正演技术实现雅可比矩阵计算和存储，得到目标函数的梯度，可进一步节省计算内存；同时根据非线性共轭梯度法法实现反演模型的迭代更新，最终获得符合观测数据的地质电阻率模型。采用压缩感知方法重构数据可以显著提升大地电磁勘探方法的工作效率，并在一定程度上加速反演收敛和得到与规则采样的数据反演相似的分辨率，将在深部矿产和能源的探测、城市基础设施的建设、清洁地热能源的开发、地质灾害的预测与预警等领域发挥积极作用。
实施例提供了一个反演算例，以验证反演算法的有效性。图4是本发明实施例提供的三维真实模型示意图。如图4所示，蓝色低阻异常体的电阻率为10ohm-m，红色高阻异常体电阻率为1000ohm-m。测点位于地面，点距1000m。整理的观测频率为1.0e-3s至10s对数等间隔分布，共9个频率数据。此算例加入5％高斯随机噪声模拟真实情况。地下背景电阻率为100ohm-m。
图5中(a)-(c)分别为本发明实施例提供的大地电磁规则完整测点反演结果，60％随机采样测点的数据重构反演结果和40％随机采样测点的数据重构反演结果的三维异常体对比图；图6中(a)-(c)分别为是本发明实施例提供的大地电磁规则完整测点反演结果，60％随机采样测点的数据重构反演结果和40％随机采样测点的数据重构反演结果的二维切片对比图，从图4，图5与图6可以看出反演结果与真实模型吻合较好，说明本发明的大地电磁三维反演算法是可靠的，且采用重构后的数据进行反演可以得到与规则完整采样的数据反演相似的模型分辨率。
图7是本发明实施例提供的大地电磁三维反演参数随迭代变化曲线图。从图7可以看出，反演过程是稳定收敛的，进一步证明反演算法的正确性。
以上所述的具体实施方式对本发明的技术方案和有益效果进行了详细说明，应理解的是以上所述仅为本发明的最优选实施例，并不用于限制本发明，凡在本发明的原则范围内所做的任何修改、补充和等同替换等，均应包含在本发明的保护范围之内。
Claims (9)
1.一种基于压缩感知重构的大地电磁数据反演方法，其特征在于，包括以下步骤：
获取随机测点分布采集的大地电磁数据，并将大地电磁数据处理成电磁阻抗数据后，对电磁阻抗数据进行分类整理；
采用基于曲波变换的压缩感知重构方式对整理后的每类电磁阻抗数据进行从随机测点分布到规则测点分布的数据重构，得到重构电磁阻抗数据；
以重构电磁阻抗数据作为观测电磁数据进行大地电磁数据反演计算。
2.根据权利要求1所述的基于压缩感知重构的大地电磁数据反演方法，其特征在于，所述对电磁阻抗数据进行分类整理，包括：
按照采样频率、阻抗类型以及电磁阻抗数据的实部和虚部进行分类，将每个采样频率下的每类电磁阻抗数据的实部或虚部作为一类数据，当采样频率为P个、阻抗类型为Q类，那整理得到P×Q×2类电磁阻抗数据，其中，2表示电磁阻抗数据的实部或虚部这两种。
3.根据权利要求1所述的基于压缩感知重构的大地电磁数据反演方法，其特征在于，所述采用基于曲波变换的压缩感知重构算法对整理后的每类电磁阻抗数据进行从随机测点分布到规则测点分布的数据重构，包括：
对随机测点分布进行稀疏采样，得到稀疏的采样矩阵S；
采用曲波变换将重构电磁阻抗数据f以正交基Ψ的形式展开，即：
f＝Ψx
其中，x是重构电磁阻抗数据f的系数序列，称为变换域系数向量，正交基Ψ＝[Ψ1,Ψ2,…,ΨN]表示N个正交基向量，N为规则测点的个数；
考虑重构电磁阻抗数据f∈RN，采样矩阵S∈RM×N，M<<N，M为随机分布测点的个数，R为实数域，则得到线性测量值
则有：
这样根据线性测量值
和已知观测的每类电磁阻抗数据y构建重构电磁阻抗数据求解的最优化问题，表示为：
其中，ε代表数据噪声，||·||1表示L1范数，||·||2表示L2范数，通过对最优化问题求解，得到变换域系数向量x，然后，根据得到的变换域系数向量x以及f＝Ψx经过稀疏逆变换得到每类电磁阻抗数据对应的重构电磁阻抗数据f。
4.根据权利要求1所述的基于压缩感知重构的大地电磁数据反演方法，其特征在于，所述以重构电磁阻抗数据作为观测电磁数据进行大地电磁数据反演计算，包括：
(a)以重构电磁阻抗数据作为观测电磁数据，计算观测电磁数据与对地下电阻率模型进行正演的预测电磁数据之间的拟合差作为数据拟合项，并构建模型粗糙度约束项，根据数据拟合项和模型约束项构建大地电磁三维正则化反演的目标函数；
(b)对目标函数进行相对于模型参数的求导以得到目标函数梯度；
(c)根据目标函数梯度并采用非线性共轭梯度法计算模型参数更新量，进而得到新地下电阻率模型；
(d)重复执行步骤(a)-步骤(c)进行迭代反演，直到收敛条件为止，输出反演得到的最终地下电阻率模型。
5.根据权利要求4所述的基于压缩感知重构的大地电磁数据反演方法，其特征在于，采用交错网格有限差分算法对地下电阻率模型进行正演得到预测电磁数据。
6.根据权利要求4所述的基于压缩感知重构的大地电磁数据反演方法，其特征在于，构建的大地电磁三维正则化反演的目标函数Φ表示为：
其中，dobs为观测电磁数据向量，dprd为由地下电阻率模型正演计算得到的预测电磁数据向量，Wd是数据协方差矩阵，为对角阵，其对角元素是观测电磁数据中噪声标准差的倒数，m为地下电阻率模型，mref为包含模型先验信息的参考模型，Wm为模型协方差矩阵，采用差分算子或拉普拉斯算子，
表示L2范数的平方，λ表示调节权重的正则化参数。
7.根据权利要求6所述的基于压缩感知重构的大地电磁数据反演方法，其特征在于，对目标函数进行相对于模型参数的求导以得到的目标函数梯度g表示为：
其中，JT表示三维大地电磁的雅克比矩阵的转置，通过对地下电阻率模型伴随正演方法推导三维大地电磁的雅克比矩阵。
8.根据权利要求6所述的基于压缩感知重构的大地电磁数据反演方法，其特征在于，根据目标函数梯度并采用非线性共轭梯度法推导模型参数更新量的方程为：
其中，Δm表示模型参数更新量，基于该方程计算模型参数更新量。
9.根据权利要求6所述的基于压缩感知重构的大地电磁数据反演方法，其特征在于，迭代反演的收敛条件为：观测电磁数据与预测电磁数据的均方根误差小于1。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202211258637.1A CN115563791B (zh) | 2022-10-14 | 基于压缩感知重构的大地电磁数据反演方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202211258637.1A CN115563791B (zh) | 2022-10-14 | 基于压缩感知重构的大地电磁数据反演方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN115563791A true CN115563791A (zh) | 2023-01-03 |
| CN115563791B CN115563791B (zh) | 2024-07-05 |
Family
ID=
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117111161A (zh) * | 2023-07-25 | 2023-11-24 | 中国平煤神马控股集团有限公司 | 基于四维高密度电阻率法的矿井赋水性动态监测方法 |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102798898A (zh) * | 2012-08-20 | 2012-11-28 | 中国地质科学院矿产资源研究所 | 大地电磁场非线性共轭梯度三维反演方法 |
| US20170075030A1 (en) * | 2015-09-15 | 2017-03-16 | Brent D. Wheelock | Accelerated Occam Inversion Using Model Remapping and Jacobian Matrix Decomposition |
| CN110908001A (zh) * | 2019-12-16 | 2020-03-24 | 吉林大学 | 一种大地电磁测深信号的重构方法及系统 |
| CN112272068A (zh) * | 2020-10-23 | 2021-01-26 | 中国人民解放军空军工程大学 | 一种基于多任务压缩感知的多样化干扰估计和抑制方法 |
| CN113221393A (zh) * | 2021-01-29 | 2021-08-06 | 吉林大学 | 一种基于非结构有限元法的三维大地电磁各向异性反演方法 |
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102798898A (zh) * | 2012-08-20 | 2012-11-28 | 中国地质科学院矿产资源研究所 | 大地电磁场非线性共轭梯度三维反演方法 |
| US20170075030A1 (en) * | 2015-09-15 | 2017-03-16 | Brent D. Wheelock | Accelerated Occam Inversion Using Model Remapping and Jacobian Matrix Decomposition |
| CN110908001A (zh) * | 2019-12-16 | 2020-03-24 | 吉林大学 | 一种大地电磁测深信号的重构方法及系统 |
| CN112272068A (zh) * | 2020-10-23 | 2021-01-26 | 中国人民解放军空军工程大学 | 一种基于多任务压缩感知的多样化干扰估计和抑制方法 |
| CN113221393A (zh) * | 2021-01-29 | 2021-08-06 | 吉林大学 | 一种基于非结构有限元法的三维大地电磁各向异性反演方法 |
Non-Patent Citations (3)
| Title |
|---|
| YUNHE LIU等: "GEM3D: A 3D inversion code for geophysical electromagnetic data based on unstructured tetrahedron grid", 《GEM 2019 XI’AN: INTERNATIONAL WORKSHOP ON GRAVITY, ELECTRICAL & MAGNETIC METHODS AND T HEIR APPLICATIONS》, 5 October 2019 (2019-10-05), pages 1 - 4 * |
| 殷长春等: "地球物理三维电磁反演方法研究动态", 《中国科学:地球科学》, vol. 50, no. 03, 31 March 2020 (2020-03-31), pages 432 - 435 * |
| 苏扬: "基于稀疏正则化的电磁勘探数据反演研究", 《万方数据》, 12 August 2022 (2022-08-12), pages 1 - 161 * |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117111161A (zh) * | 2023-07-25 | 2023-11-24 | 中国平煤神马控股集团有限公司 | 基于四维高密度电阻率法的矿井赋水性动态监测方法 |
| CN117111161B (zh) * | 2023-07-25 | 2024-05-24 | 中国平煤神马控股集团有限公司 | 基于四维高密度电阻率法的矿井赋水性动态监测方法 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Jahandari et al. | 3-D minimum-structure inversion of magnetotelluric data using the finite-element method and tetrahedral grids | |
| Newman et al. | Solution accelerators for large-scale three-dimensional electromagnetic inverse problems | |
| Lochbühler et al. | Summary statistics from training images as prior information in probabilistic inversion | |
| WO2011115921A2 (en) | Uncertainty estimation for large-scale nonlinear inverse problems using geometric sampling and covariance-free model compression | |
| Wang et al. | Multiparameter full-waveform inversion of 3-D on-ground GPR with a modified total variation regularization scheme | |
| CN112949134A (zh) | 基于非结构有限元方法的地-井瞬变电磁反演方法 | |
| Kotsi et al. | Uncertainty quantification in time-lapse seismic imaging: a full-waveform approach | |
| Udosen et al. | Automated optimization of electrode locations for electrical resistivity tomography | |
| Su et al. | Sparse-promoting 3-D airborne electromagnetic inversion based on shearlet transform | |
| Cai et al. | Effective 3D-transient electromagnetic inversion using finite-element method with a parallel direct solver | |
| Li et al. | Using wavelet packet denoising and a regularized ELM algorithm based on the LOO approach for transient electromagnetic inversion | |
| Chen et al. | Three-dimensional inversion of controlled-source audio-frequency magnetotelluric data based on unstructured finite-element method | |
| Narciso et al. | A comparison between Kalman ensemble generator and geostatistical frequency-domain electromagnetic inversion: The impacts on near-surface characterization | |
| Tobely et al. | Position detection of unexploded ordnance from airborne magnetic anomaly data using 3-D self organized feature map | |
| Zhang et al. | 3D inversion of large-scale frequency-domain airborne electromagnetic data using unstructured local mesh | |
| Liao et al. | 2D inversion of magnetotelluric data using deep learning technology | |
| Ghalenoei et al. | Joint gravity and magnetic inversion with trans-dimensional alpha shapes and autoregressive noise models | |
| Dubey et al. | Joint inversion of gravity and gravity gradient and its application to mineral exploration | |
| Liu et al. | Application of sample-compressed neural network and adaptive-clustering algorithm for magnetotelluric inverse modeling | |
| Wang et al. | Geophysical electromagnetic modeling and evaluation: a review | |
| Cao et al. | 3-D Crosswell electromagnetic inversion based on IRLS norm sparse optimization algorithms | |
| CN115563791B (zh) | 基于压缩感知重构的大地电磁数据反演方法 | |
| Zhong et al. | Electrical resistivity tomography with smooth sparse regularization | |
| CN115755199A (zh) | 一种实用的非结构网格三维电磁反演平滑正则化方法 | |
| CN115563791A (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 | ||
| GR01 | Patent grant |