CN113516754B - 一种基于磁异常模量数据的三维可视化成像方法 - Google Patents
一种基于磁异常模量数据的三维可视化成像方法 Download PDFInfo
- Publication number
- CN113516754B CN113516754B CN202110282040.XA CN202110282040A CN113516754B CN 113516754 B CN113516754 B CN 113516754B CN 202110282040 A CN202110282040 A CN 202110282040A CN 113516754 B CN113516754 B CN 113516754B
- Authority
- CN
- China
- Prior art keywords
- matrix
- level set
- data
- set function
- batch
- 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.)
- Active
Links
- 238000003384 imaging method Methods 0.000 title claims abstract description 32
- 230000002159 abnormal effect Effects 0.000 title claims abstract description 15
- 230000000007 visual effect Effects 0.000 title claims abstract description 7
- 239000011159 matrix material Substances 0.000 claims abstract description 108
- 238000000354 decomposition reaction Methods 0.000 claims abstract description 37
- 230000035945 sensitivity Effects 0.000 claims abstract description 24
- 238000005457 optimization Methods 0.000 claims abstract description 18
- 238000005516 engineering process Methods 0.000 claims abstract description 16
- 230000000903 blocking effect Effects 0.000 claims abstract description 12
- 238000007781 pre-processing Methods 0.000 claims abstract description 11
- 238000000034 method Methods 0.000 claims description 27
- 238000012800 visualization Methods 0.000 claims 1
- 238000004364 calculation method Methods 0.000 abstract description 9
- 230000000694 effects Effects 0.000 abstract description 8
- 230000006870 function Effects 0.000 description 102
- 238000005259 measurement Methods 0.000 description 11
- 238000012217 deletion Methods 0.000 description 4
- 230000037430 deletion Effects 0.000 description 4
- 238000001514 detection method Methods 0.000 description 4
- 230000010354 integration Effects 0.000 description 4
- 238000010586 diagram Methods 0.000 description 3
- 230000015572 biosynthetic process Effects 0.000 description 2
- 238000005755 formation reaction Methods 0.000 description 2
- 239000000463 material Substances 0.000 description 2
- 230000007935 neutral effect Effects 0.000 description 2
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 description 1
- 102100029469 WD repeat and HMG-box DNA-binding protein 1 Human genes 0.000 description 1
- 101710097421 WD repeat and HMG-box DNA-binding protein 1 Proteins 0.000 description 1
- 230000002547 anomalous effect Effects 0.000 description 1
- 238000005520 cutting process Methods 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 238000000605 extraction Methods 0.000 description 1
- 238000005286 illumination Methods 0.000 description 1
- 238000009434 installation Methods 0.000 description 1
- 239000000696 magnetic material Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 239000002245 particle Substances 0.000 description 1
- 230000008569 process Effects 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 238000004441 surface measurement Methods 0.000 description 1
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
-
- 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/14—Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms
-
- 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
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/70—Denoising; Smoothing
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Mathematical Physics (AREA)
- Mathematical Optimization (AREA)
- Data Mining & Analysis (AREA)
- Pure & Applied Mathematics (AREA)
- Computational Mathematics (AREA)
- Mathematical Analysis (AREA)
- Software Systems (AREA)
- Algebra (AREA)
- Databases & Information Systems (AREA)
- General Engineering & Computer Science (AREA)
- Computing Systems (AREA)
- Geometry (AREA)
- Computer Graphics (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明公开了一种基于磁异常模量数据的三维可视化成像方法，包括：获取磁异常模量数据和先验信息；构建正演模型积分核矩阵；采用小批量随机梯度下降下的分块截断奇异值分解技术对正演模型积分核矩阵进行预处理；构建优化目标函数；初始化水平集函数；采用小批量随机梯度下降下的分块截断奇异值分解技术计算预测模量数据并代入计算优化目标函数，判断是否满足迭代终止条件；若不满足，则采用小批量随机梯度下降下的分块截断奇异值分解技术计算梯度流；根据梯度流更新水平集函数；迭代直至满足终止条件，根据输出的水平集函数的零水平集面确定磁性体的三维成像结构。本发明引入了水平集函数的敏感度加权项，提高了深层区域的水平集反演效果。
Description
技术领域
本发明涉及地球物理成像技术领域，尤其是涉及一种基于磁异常模量数据的三维可视化成像方法。
背景技术
基于磁异常数据对磁性体成像在地球物理探测中有广泛的应用，相应的探测方法称为磁法勘探。磁法勘探利用埋藏物的磁性差异所引起的磁异常来查明地下构造，其广泛地运用于研究大地构造和寻找磁性矿体。近年来，磁法勘探还在水文、工程、环境勘察方面崭露头角，在探测地下热源、含水破碎带、地下管道、地下电缆等方面均取得了良好的效果。此外，基于磁异常数据的成像方法还被用于军事探测(例如舰艇、飞行器的探测)，以及大型集装箱货物检测中。
目前，磁异常水平集反演方法采用的是线性的磁异常数据，包括磁异常投影数据与磁异常梯度数据，其数据类型对噪声敏感，极易受到噪声干扰，相应的成像结果也很受影响。此外，现有的磁异常水平集反演方法对深层区域的成像效果不佳，在海量观测数据的应用场景中计算效率低。
发明内容
本发明旨在至少解决现有技术中存在的技术问题之一。为此，本发明提出一种基于磁异常模量数据的三维可视化成像方法。本发明提出的方法运用磁异常模量数据反演水平集函数，采用的数据类型可以消除随机噪声的影响，相应的成像结果更为稳定。本发明引入了水平集函数的敏感度加权项，使得刻画磁性体边界面的水平集函数在数据敏感度低的区域更易分布，提高了深层区域的磁性体成像效果。此外，本发明提出了小批量随机梯度下降下的分块截断奇异值分解技术，采用该技术计算预测模量数据以及目标梯度流，有效提高了成像方法的计算效率。
根据本发明的第一方面实施例的一种基于磁异常模量数据的三维可视化成像方法，包括：
获取磁异常模量数据和先验信息；
构建正演模型积分核矩阵；
采用小批量随机梯度下降下的分块截断奇异值分解技术对正演模型积分核矩阵进行预处理操作；
引入水平集函数敏感度加权项，构建优化目标函数；
初始化水平集函数；
采用小批量随机梯度下降下的分块截断奇异值分解技术计算预测模量数据；
根据预测模量数据计算优化目标函数值，判断是否满足迭代终止条件；
若不满足迭代终止条件，则采用小批量随机梯度下降下的分块截断奇异值分解技术计算优化目标函数对水平集函数的目标梯度流；
根据目标梯度流更新水平集函数；
回到“采用小批量随机梯度下降下的分块截断奇异值分解技术计算预测模量数据”的步骤。
若满足迭代终止条件，则输出水平集函数；
根据水平集函数确定磁性体对应的三维成像结构。
根据本发明的一些实施例，引入水平集函数敏感度加权项，构建优化目标函数，包括：引入水平集函数敏感度加权项、批数据偏差项、水平集函数正则项；通过所述水平集函数敏感度加权项、所述批数据偏差项、所述水平集函数正则项的线性组合构建优化目标函数。
根据本发明的一些实施例，采用小批量随机梯度下降下的分块截断奇异值分解技术对正演模型积分核矩阵进行预处理操作，包括：对正演模型积分核矩阵按列分块，使得一个分块子矩阵对应地下区域的同一层(深度)；对每一个分块子矩阵进行奇异值分解操作，使得每个分块子矩阵对应三个矩阵的乘积；对每个分块子矩阵分解后的三个矩阵进行截断操作，并对截断后的矩阵作乘积合并，使得每个分块子矩阵对应两个矩阵；将预处理操作后的矩阵块存储备用：对每一个分块子矩阵，存储其对应的两个矩阵。
根据本发明的一些实施例，采用小批量随机梯度下降下的分块截断奇异值分解技术计算预测模量数据，包括：按层读入分块子矩阵对应的两个矩阵；对每一个分块子矩阵对应的两个矩阵，取第一个矩阵的特定行构成新矩阵，并与第二个矩阵相乘，得到这一分块子矩阵产生的对应本批次的预测分量数据；将所有分块子矩阵产生的预测分量数据求和，得到本批次对应的预测分量数据；计算预测分量数据的模量，即得预测模量数据。
根据本发明的一些实施例，当计算优化目标函数对水平集函数的目标梯度流时，采用小批量随机梯度下降下的分块截断奇异值分解技术计算批数据偏差项对水平集函数的梯度流，包括：按层读入分块子矩阵对应的两个矩阵；对每一个分块子矩阵对应的两个矩阵，取第一个矩阵的特定行构成新矩阵，并与第二个矩阵参与运算，得到该层区域上批数据偏差项对水平集函数的梯度流；将每层区域上的梯度流矩阵合并成为一个大矩阵，得到整个区域上批数据偏差项对水平集函数的梯度流。
本发明的附加方面和优点将在下面的描述中部分给出，部分将从下面的描述中变得明显，或通过本发明的实践了解到。
附图说明
本发明的上述和/或附加的方面和优点从结合下面附图对实施例的描述中将变得明显和容易理解，其中：
图1为本发明实施例的一种基于磁异常模量数据的三维可视化成像方法的流程示意图；
图2a为本发明实施例的地表测量点分布的示意图；
图2b为本发明实施例的磁异常模量数据的示意图；
图3a为本发明实施例的水平集函数初始值的零水平集面的示意图；
图3b为本发明实施例的水平集函数初始值的示意图：x＝0.5(km)截面上的水平集函数初始值。；
图4a为本发明实施例的三维成像结果；
图4b为本发明实施例的三维成像结果在x＝0.25(km)截面上的磁性体形状；
图4c为本发明实施例的三维成像结果在x＝0.75(km)截面上的磁性体形状。
具体实施方式
下面详细描述本发明的实施例，实施例的示例在附图中示出，其中自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附图描述的实施例是示例性的，仅用于解释本发明，而不能理解为对本发明的限制。
本发明的描述中，除非另有明确的限定，设置、安装、连接等词语应做广义理解，所属技术领域技术人员可以结合技术方案的具体内容合理确定上述词语在本发明中的具体含义。
首先，对本申请中涉及的若干名词进行解析：
1、Heaviside函数：H(x)＝1,x≥0；H(x)＝0,x<0.
2、delta函数：一般指狄拉克δ函数，狄拉克δ函数是一个广义函数，在物理学中常用其表示质点、点电荷等理想模型的密度分布，该函数在除了零以外的点取值都等于零，而其在整个定义域上的积分等于1。计算上我们采用数值delta函数，具体形式见实施例。
3、拉普拉斯运算：对一个函数f做拉普拉斯运算，即求f所有的非混合二阶偏导数之和。例如在三维空间中，
参照图1，根据本发明第一方面实施例的一种基于磁异常模量数据的三维可视化成像方法，包括：
步骤S1000，获取磁异常模量数据和先验信息。
其中，磁异常模量数据是磁性体在诱导磁场下产生的磁异常分量数据的模量；先验信息包括：诱导磁场的方位角(即诱导磁场方向)、磁性体的平均磁化率数值κ0。如图2a所示，在本发明实施例中，磁性体的测量点分布在地表上方100米，共有10000个数据测量点。图2b展示了本发明实施例获取的磁异常模量数据，图2b中绘出了磁异常模量数据在20nT、40nT、60nT、80nT和100nT处的等高线，其中nT表示单位纳特斯拉(nanotesla)。磁异常模量数据可以由如下公式描述：
其中，G表示测量点所在区域，d(r)表示测量点r处的磁异常模量，g1(r)、g2(r)、g3(r)分别表示磁异常分量数据在x(东)、y(北)、z(深度)三个方向的分量。由于磁异常模量数据是非线性数据，非线性叠加在统计上可以消除随机噪声的影响，使得相应的成像结果更为稳定。
步骤S1010，构建正演模型积分核矩阵。
具体的，先构建磁异常数据正演模型，然后分割地下区域得到网格节点坐标，并录入数据测量点的位置坐标，最后基于这些信息构建积分核矩阵，得到三个积分核矩阵：Kt(t＝1,2,3)。具体如下所示：
(i)构建正演模型。
本发明实施例使用以下公式模拟磁异常分量数据gt(t＝1,2,3)：
其中，B0表示诱导磁场强度，Ω表示地下区域，κ表示磁化率函数，Kt(t＝1,2,3)表示积分核函数。磁异常模量数据d运用分量数据通过公式(1)得到。
本发明实施例运用水平集函数表达式对磁化率函数κ建模：
其中为引入的水平集函数，H(·)表示标准的Heaviside函数，κ0是步骤S1000中输入先验信息中的平均磁化率数值。
积分核函数Kt(t＝1,2,3)的形式如下：
其中，表示诱导磁场的单位方向向量，(·)表示标准的向量内积。
(ii)分割地下区域Ω。
例如在矩形计算区域，x方向取nx个网格点，y方向取ny个网格点，z方向取nz个网格点，则地下区域被分割成N＝nx·ny·nz个网格点。
(iii)构建积分核矩阵。
记数据测量点位置坐标分别为r1,r2,…,rM，地下区域网格点坐标分别为在这些离散点上计算公式(4)中的积分核函数/>得到三个积分核矩阵：Kt(t＝1,2,3)。其中矩阵Kt为M行N列，它的第i行j列元素为/>
步骤S1020，采用小批量随机梯度下降下的分块截断奇异值分解技术对正演模型积分核矩阵进行预处理操作。
首先，对正演模型积分核矩阵Kt(t＝1,2,3)进行分块操作，使得一个分块子矩阵对应地下区域的同一层(深度)。随后，对每一个分块子矩阵进行奇异值分解操作，使得每个分块子矩阵对应三个矩阵的乘积。接下来，对每个分块子矩阵分解后的三个矩阵进行截断操作，并对截断后的矩阵作乘积合并，使得每个分块子矩阵对应两个矩阵。最后，将预处理操作后的矩阵块存储备用：对每一个分块子矩阵，存储其对应的两个矩阵。具体如下所示：
(i)积分核矩阵分块操作。
将正演模型积分核矩阵Kt按列分块，使得一个分块子矩阵对应地下区域的同一层(深度)，即同一个分块子矩阵对应的网格点的深度坐标(z坐标)相同。
例如记Kt的分块结果为一共分为nz个子矩阵，每个子矩阵为M行nx·ny列；/>的第i行j列元素为/>对应同一个/>的所有网格点的z坐标(深度坐标)相同。
(ii)子矩阵奇异值分解操作。
对每一个子矩阵进行奇异值分解：
其中()T表示矩阵转置，为对角阵。
(iii)截断操作。
选择适当的阈值将对角阵/>中小于阈值/>的元素去掉，得到低阶的近似矩阵保留/>和/>中与近似矩阵/>相对应的列，得到/>和/>计算/>
(iv)矩阵预处理结果的存储。
将矩阵(t＝1,2,3；h＝1,2,...nz)存储备用，每一个矩阵存为一个文件，一共存为6nz个文件。
步骤S1030，引入水平集函数敏感度加权项，构建优化目标函数。
由于磁异常模量数据通常在探测区域的一侧获得(例如地表)，对于远离数据测量面的深层区域，磁性体产生的信号较弱，极易被浅层区域磁性体产生的磁异常掩盖，使得成像方法对深层区域的成像效果不佳。我们引入水平集函数的敏感度加权项，使得刻画磁性体边界面的水平集函数在数据敏感度低的区域更易分布，提高了深层区域的磁性体成像效果。具体如下所示：
(i)引入水平集函数敏感度加权项Em，具体形式为：
其中权函数表达式如下，
权函数随积分核函数/>的值减小而减小，而积分核函数/>反映了测量数据的敏感度，因此权函数/>在数据敏感度低的区域(例如深层区域)有更小的数值。这使得在优化(极小化)Em时，水平集函数/>在深层区域更易分布，有利于提高深层区域的成像效果。
(ii)引入批数据偏差项Ed。
本发明实施例使用小批量随机梯度下降的优化方法。该方法逐批地使用测量数据，每使用过一遍所有测量数据便称完成了一期更新，当一期更新结束时便开始下一期更新。在每一期更新的开始，先将所有测量数据进行一次随机排列，然后将数据分批，使得每批数据具有相同的个数。随后依次逐批地使用这些数据对水平集函数进行迭代更新。
例如在第R期的第k次迭代，记该批次用到的测量数据下标集为Sk＝{i1,k,i2,k,…,ib,k}，其中b为每批数据的容量，引入针对批数据的数据偏差项Ed：
其中，表示第in,k个预测模量数据，/>表示第in,k个测量模量数据。
(iii)引入水平集函数正则项Er，形式如下：
其中，表示梯度运算；|·|表示取模长运算，例如在三维情况下，
(iv)构建优化目标函数E。
在本发明实施例中，优化目标函数E为批数据偏差项Ed、水平集函数敏感度加权项Em与水平集函数正则项Er的线性组合：
E(φ)＝Ed(φ)+μ1Em(φ)+μ2Er(φ). (10)
其中，μ1和μ2两个参数用来控制Em和Er在优化目标函数中的权重。
步骤S1040，初始化水平集函数。
设定水平集函数的初始值，具体操作如下。
(i)设定未知磁性体的初始边界面。在对未知磁性体分布完全没有先验信息的情况下，通常将初始边界面设置为球面或者椭球面。
(ii)将水平集函数的初始值设定为到初始边界面的有向距离函数。例如可以先将水平集函数/>的值设置为在初始边界面内部为1、在初始边界面外部为-1，然后做标准的水平集重初始化操作，使得/>的值成为一个有向距离。图3a和图3b展示了本发明实施例中的水平集函数的初始值情况。图3a示出了水平集函数初始值的零水平集面，其为一个球面，其中点A为亮点(即光照点)，光朝向点A照射到图3a所示的球体上；图3b展示了x＝0.5(km)截面上的水平集函数初始值情况。
步骤S1050，采用小批量随机梯度下降下的分块截断奇异值分解技术计算预测模量数据。
以第R期的第k次迭代为例，记本批次用到的测量数据下标集为Sk＝{i1,k,i2,k,…,ib,k}。
(I)计算本批次对应的预测数据分量
①原理解释
根据正演模型(公式(2))，磁异常分量数据的矩阵乘法表达式如下：
其中表示网格宽度，/>是由正演模型积分核矩阵Kt的第i1,k,…,ib,k行组成的矩阵，κ是由各个网格点上的磁化率组成的N维行向量。本发明实施例基于小批量随机梯度下降下的分块截断奇异值分解来计算(11)式。对照步骤S1020中矩阵Kt的按列分块方法，同样对/>和κ按列分块：
则有
其中由公式(5)中/>的第i1,k,…,ib,k行构成。在步骤S1020中的截断奇异值分解下，有/>取出/>的第i1,k,…,ib,k行构成矩阵/>则有
②实施方式
(i)取为b维零向量，即/>设置h＝1。
(ii)将矩阵读入内存，用/>的第i1,k,…,ib,k行构成/>删除/>
(iii)赋值然后令h＝h+1。
(iv)h≤nz时，回到第(ii)步；h>nz时，终止循环，令
(II)计算本批次对应的预测模量数据
步骤S1060，根据预测模量数据计算优化目标函数值，判断是否满足迭代终止条件。具体操作如下。
(i)根据当前的水平集函数值和预测模量数据，按照上述公式(10)计算优化目标函数E(φ)的值。
(ii)当目标函数值小于某个预设值时，判定为满足迭代终止条件；当目标函数值大于等于预设值时，判定为不满足迭代终止条件。
步骤S1070，若不满足迭代终止条件，则采用小批量随机梯度下降下的分块截断奇异值分解技术计算优化目标函数对水平集函数的目标梯度流。
优化目标函数E对水平集函数φ的梯度由三项求和得到：
以第R期的第k次迭代为例，记本批次用到的测量数据下标集为Sk＝{i1,k,i2,k,…,ib,k}。
(I)的计算：小批量随机梯度下降下的分块截断奇异值分解技术
①原理解释
其中是引入的数值delta函数，具体形式如下
表示零水平集附近的小区域，参数τ用以控制迭代区域宽度，建议选取为0.5倍网格宽度；/>表示指示函数：/>
记
则有
为了节省内存并且提高算法效率，本发明实施例基于小批量随机梯度下降下的分块截断奇异值分解来计算(20)式中的(20)式的矩阵乘法形式如下：
其中Dt是b维行向量：
是由正演模型积分核矩阵Kt的第i1,k,…,ib,k行组成的矩阵；/>是由指示函数/>组成的N维对角矩阵。考虑到(12)式中/>的分块方式，将对角阵Λ分块：
使得
因此向量F可以按照分块的方式计算：
考虑到(15)式中对的截断奇异值分解近似，有
②实施方式
(1)计算行向量F
(i)令h＝1。
(ii)取Fh为nx·ny维的零向量。
(iii)令t＝1。
(iv)将矩阵读入内存，用/>的第i1,k,…,ib,k行构成/>删除/>
(v)赋值
(vi)令t＝t+1。若t≤3，回到第(iv)步，否则进入第(vii)步。
(vii)令h＝h+1。若h≤nz，回到第(ii)步；否则终止循环，得到向量F：
(2)根据行向量F计算其中Fj表示行向量F的第j个元素。
(II)的计算：
(III)的计算：
其中Δ表示拉普拉斯运算。
步骤S1080，根据目标梯度流更新水平集函数。
(i)沿目标梯度流负方向迭代水平集函数φ：
其中Δt表示选取的迭代步长。
(ii)对水平集函数φ做标准的水平集重初始化操作，得到更新后的水平集函数。
回到步骤S1050。
更新水平集函数后，回到步骤S1050重新计算预测模量数据，直至步骤S1060中判定满足迭代终止条件。
步骤S1090，若满足迭代终止条件，则输出水平集函数。
步骤S1110，根据水平集函数确定磁性体对应的三维成像结构。
水平集函数的零水平集面，即对应反演磁性体的边界面。通过展示水平集函数的零水平集面，即得到磁性体的三维成像结构。图4展示了本发明实施例的反演结果。图4a展示了输出的水平集函数的零水平集面，即反演磁性体的边界面，给出了磁性体的三维直观结构；图4b和图4c展示了两个切面上的磁性体形状，图4b为纵切面x＝0.25(km)上的磁性体边界面，图4c为纵切面x＝0.75(km)上的磁性体边界面。
在本发明的一些实施例中，引入水平集函数敏感度加权项，构建优化目标函数，包括：
引入水平集函数敏感度加权项、批数据偏差项和水平集函数正则项。引入的水平集函数敏感度加权项Em，具体形式如上述公式(6)、(7)所示；引入的批数据偏差项Ed，具体形式如上述公式(8)所示；引入的水平集函数正则项Er，具体形式如上述公式(9)所示。
构建的优化目标函数为水平集函数敏感度加权项、批数据偏差项和水平集函数正则项的线性组合，具体形式如上述公式(10)所示。
通过引入水平集函数的敏感度加权项，使得刻画磁性体边界面的水平集函数在数据敏感度低的区域更易分布，提高了深层区域的磁性体成像效果。
在本发明的一些实施例中，采用小批量随机梯度下降下的分块截断奇异值分解技术对正演模型积分核矩阵进行预处理操作，包括：
将正演模型积分核矩阵Kt按列分块，使得一个分块子矩阵对应地下区域的同一层(深度)，即同一个分块子矩阵对应的网格点的深度坐标(z坐标)相同。例如记Kt的分块结果为/>一共分为nz个子矩阵，每个子矩阵为M行nx·ny列；其中/>的第i行j列元素为/>对应同一个/>的所有网格点/>的z坐标(深度坐标)相同。
对每一个分块子矩阵进行奇异值分解操作，使得每个分块子矩阵对应三个矩阵的乘积。具体的，通过上述公式(5)对每一个子矩阵进行奇异值分解，得到分解矩阵和/>
对每个分块子矩阵分解后的三个矩阵进行截断操作，并对截断后的矩阵作乘积合并，使得每个分块子矩阵对应两个矩阵。具体的，选择适当的阈值将分解结果的对角阵中小于阈值/>的元素去掉，得到低阶的近似矩阵/>保留/>和/>中与近似矩阵/>相对应的列，得到/>和/>计算/>
对矩阵预处理结果进行存储。具体的，将矩阵(t＝1,2,3；h＝1,2,...nz)存储备用，每一个矩阵存为一个文件，一共存为6nz个文件。
在本发明的一些实施例中，采用小批量随机梯度下降下的分块截断奇异值分解技术计算预测模量数据，包括：按层读入分块子矩阵对应的两个矩阵；对每一个分块子矩阵对应的两个矩阵，取第一个矩阵的特定行构成新矩阵，并与第二个矩阵相乘，得到这一分块子矩阵产生的对应本批次的预测分量数据；将所有分块子矩阵产生的预测分量数据求和，得到本批次对应的预测分量数据；计算预测分量数据的模量，即得预测模量数据。具体实施方式如下。
(I)计算本批次对应的预测数据分量Sk＝{i1,k,i2,k,…,ib,k},t＝1,2,3.
(i)取为b维零向量，即/>设置h＝1。
(ii)将矩阵读入内存，用/>的第i1,k,…,ib,k行构成/>删除/>
(iii)赋值然后令h＝h+1。
(iv)h≤nz时，回到第(ii)步；h>nz时，终止循环，令
(II)计算本批次对应的预测模量数据
在本发明的一些实施例中，当计算优化目标函数E对水平集函数φ的目标梯度流时，采用小批量随机梯度下降下的分块截断奇异值分解技术计算批数据偏差项Ed对水平集函数φ的梯度流/>包括：按层读入分块子矩阵对应的两个矩阵；对每一个分块子矩阵对应的两个矩阵，取第一个矩阵的特定行构成新矩阵，并与第二个矩阵参与运算，得到该层区域上批数据偏差项对水平集函数的梯度流；将每层区域上的梯度流矩阵合并成为一个大矩阵，得到整个区域上批数据偏差项对水平集函数的梯度流。具体实施方式如下。
(1)计算行向量F
(i)令h＝1。
(ii)取Fh为nx·ny维的零向量。
(iii)令t＝1。
(iv)将矩阵读入内存，用/>的第i1,k,…,ib,k行构成/>删除/>
(v)赋值
(vi)令t＝t+1。若t≤3，回到第(iv)步，否则进入第(vii)步。
(vii)令h＝h+1。若h≤nz，回到第(ii)步；否则终止循环，得到向量F：
(2)根据行向量F计算批数据偏差项Ed对水平集函数φ的梯度流 其中Fj表示行向量F的第j个元素。
上面结合附图对本发明实施例作了详细说明，但是本发明不限于上述实施例，在所属技术领域普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下作出各种变化。
以上所描述的装置实施例仅仅是示意性的，其中作为分离部件说明的单元是或者也可以不是物理上分开的，即可以位于一个地方，或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或者全部模块来实现本实施例方案的目的。
在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示意性实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不一定指的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。
尽管已经示出和描述了本发明的实施例，本领域的普通技术人员可以理解：在不脱离本发明的原理和宗旨的情况下可以对这些实施例进行多种变化、修改、替换和变型，本发明的范围由权利要求及其等同物限定。
Claims (3)
1.一种基于磁异常模量数据的三维可视化成像方法，其特征在于，包括：
获取磁异常模量数据和先验信息；
构建正演模型积分核矩阵；
采用小批量随机梯度下降下的分块截断奇异值分解技术对正演模型积分核矩阵进行预处理操作，具体包括：对正演模型积分核矩阵按列分块，使得一个分块子矩阵对应地下区域的同一层；对每一个分块子矩阵进行奇异值分解操作，使得每个分块子矩阵对应三个矩阵的乘积；对每个分块子矩阵分解后的三个矩阵进行截断操作，并对截断后的矩阵作乘积合并，使得每个分块子矩阵对应两个矩阵；将预处理操作后的矩阵块存储备用：对每一个分块子矩阵，存储其对应的两个矩阵；
引入水平集函数敏感度加权项、批数据偏差项，构建优化目标函数；
初始化水平集函数；
采用小批量随机梯度下降下的分块截断奇异值分解技术计算预测模量数据；
根据预测模量数据计算优化目标函数值，判断是否满足迭代终止条件；
若不满足迭代终止条件，则采用小批量随机梯度下降下的分块截断奇异值分解技术计算优化目标函数对水平集函数的目标梯度流，具体包括：按层读入分块子矩阵对应的两个矩阵；对每一个分块子矩阵对应的两个矩阵，取第一个矩阵的特定行构成新矩阵，并与第二个矩阵参与运算，得到该层区域上批数据偏差项对水平集函数的梯度流；将每层区域上的梯度流矩阵合并成为一个大矩阵，得到整个区域上批数据偏差项对水平集函数的梯度流；
根据目标梯度流更新水平集函数；
回到采用小批量随机梯度下降下的分块截断奇异值分解技术计算预测模量数据的步骤；
若满足迭代终止条件，则输出水平集函数；
根据水平集函数确定磁性体对应的三维成像结构。
2.根据权利要求1所述的一种基于磁异常模量数据的三维可视化成像方法，其特征在于，引入水平集函数敏感度加权项、批数据偏差项，构建优化目标函数，包括：
引入水平集函数敏感度加权项、批数据偏差项、水平集函数正则项；
通过所述水平集函数敏感度加权项、所述批数据偏差项、所述水平集函数正则项的线性组合构建优化目标函数。
3.根据权利要求1所述的一种基于磁异常模量数据的三维可视化成像方法，其特征在于，采用小批量随机梯度下降下的分块截断奇异值分解技术计算预测模量数据，包括：
按层读入分块子矩阵对应的两个矩阵；
对每一个分块子矩阵对应的两个矩阵，取第一个矩阵的特定行构成新矩阵，并与第二个矩阵相乘，得到这一分块子矩阵产生的对应本批次的预测分量数据；
将所有分块子矩阵产生的预测分量数据求和，得到本批次对应的预测分量数据；
计算预测分量数据的模量，即得预测模量数据。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110282040.XA CN113516754B (zh) | 2021-03-16 | 2021-03-16 | 一种基于磁异常模量数据的三维可视化成像方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110282040.XA CN113516754B (zh) | 2021-03-16 | 2021-03-16 | 一种基于磁异常模量数据的三维可视化成像方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113516754A CN113516754A (zh) | 2021-10-19 |
| CN113516754B true CN113516754B (zh) | 2024-05-03 |
Family
ID=78061751
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110282040.XA Active CN113516754B (zh) | 2021-03-16 | 2021-03-16 | 一种基于磁异常模量数据的三维可视化成像方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113516754B (zh) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114518599B (zh) * | 2022-01-28 | 2024-05-17 | 北京大学 | 一种自适应的磁异常目标成像与检测方法 |
| CN116068903B (zh) * | 2023-04-06 | 2023-06-20 | 中国人民解放军国防科技大学 | 一种闭环系统鲁棒性能的实时优化方法、装置及设备 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104008249A (zh) * | 2014-06-09 | 2014-08-27 | 桂林电子科技大学 | 基于动态轮廓模型的地面核磁共振反演方法 |
| CN110097550A (zh) * | 2019-05-05 | 2019-08-06 | 电子科技大学 | 一种基于深度学习的医学图像分割方法及系统 |
| CN110338840A (zh) * | 2015-02-16 | 2019-10-18 | 深圳迈瑞生物医疗电子股份有限公司 | 三维成像数据的显示处理方法和三维超声成像方法及系统 |
| WO2020260671A1 (en) * | 2019-06-26 | 2020-12-30 | Cerebriu A/S | An improved medical scan protocol for in-scanner patient data acquisition analysis |
Family Cites Families (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| AU2010347706B2 (en) * | 2010-03-03 | 2015-04-23 | Brain Research Institute Foundation Pty Ltd | Image processing system |
| US11055908B2 (en) * | 2018-08-27 | 2021-07-06 | Exxonmobil Research And Engineering Company | Inversion of large, nearly-homogeneous geobodies via ultra low-dimensional shape representation using orthogonal basis functions |
| US11315221B2 (en) * | 2019-04-01 | 2022-04-26 | Canon Medical Systems Corporation | Apparatus and method for image reconstruction using feature-aware deep learning |
-
2021
- 2021-03-16 CN CN202110282040.XA patent/CN113516754B/zh active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104008249A (zh) * | 2014-06-09 | 2014-08-27 | 桂林电子科技大学 | 基于动态轮廓模型的地面核磁共振反演方法 |
| CN110338840A (zh) * | 2015-02-16 | 2019-10-18 | 深圳迈瑞生物医疗电子股份有限公司 | 三维成像数据的显示处理方法和三维超声成像方法及系统 |
| CN110097550A (zh) * | 2019-05-05 | 2019-08-06 | 电子科技大学 | 一种基于深度学习的医学图像分割方法及系统 |
| WO2020260671A1 (en) * | 2019-06-26 | 2020-12-30 | Cerebriu A/S | An improved medical scan protocol for in-scanner patient data acquisition analysis |
Non-Patent Citations (5)
| Title |
|---|
| A level-set algorithm for the inverse problem of full magnetic gradient tensor data;Wenbin Li et al.;Applied Mathematics Letters;20200422;第107卷;第1-6页 * |
| A multiple level set method for three-dimensional inversion of magnetic data;Wenbin Li, Wangtao Lu, Jianliang Qian, and Yaoguo Li;GEOPHYSICS;20170930;第5卷(第82期);第150-239页 * |
| Joint inversion of surface and borehole magnetic data: A level-set approach;Wenbin Li, Jianliang Qian, and Yaoguo Li;GEOPHYSICS;20200131;第5卷(第80期);第1-58页 * |
| SIMULTANEOUSLY RECOVERING BOTH DOMAIN AND VARYING DENSITY IN INVERSE GRAVIMETRY BY EFFICIENT LEVEL-SET METHODS;Wenbin Li，Jianliang Qian;Inverse Problems and Imaging;20201130;第3卷(第15期);第387-413页 * |
| 基于χ~2准则的磁梯度张量3D聚焦反演方法;李金朋;张英堂;范红波;李志宁;尹刚;刘敏;;武汉大学学报(信息科学版);20180205(第02期);第255-261页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113516754A (zh) | 2021-10-19 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Gumerov et al. | Fast radial basis function interpolation via preconditioned Krylov iteration | |
| US9619590B2 (en) | Uncertainty estimation for large-scale nonlinear inverse problems using geometric sampling and covariance-free model compression | |
| Martínez et al. | PSO: A powerful algorithm to solve geophysical inverse problems: Application to a 1D-DC resistivity case | |
| US8571842B2 (en) | Method of determining parameter from sparse measurement data | |
| CN113516754B (zh) | 一种基于磁异常模量数据的三维可视化成像方法 | |
| Laloy et al. | Mass conservative three‐dimensional water tracer distribution from Markov chain Monte Carlo inversion of time‐lapse ground‐penetrating radar data | |
| KR100831932B1 (ko) | 오일러 해를 이용한 지하공동의 3차원 중력역산 방법 및이를 이용한 3차원 영상화 방법 | |
| CN109636912A (zh) | 应用于三维声呐图像重构的四面体剖分有限元插值方法 | |
| Wang et al. | Full waveform inversion based on the ensemble Kalman filter method using uniform sampling without replacement | |
| CN115238550A (zh) | 自适应非结构网格的滑坡降雨的地电场数值模拟计算方法 | |
| Ghalenoei et al. | Trans-dimensional gravity and magnetic joint inversion for 3-D earth models | |
| CN114896564A (zh) | 采用自适应泰森多边形参数化的瞬变电磁二维贝叶斯反演方法 | |
| CN113109883B (zh) | 基于等参变换全球离散网格球坐标下卫星重力场正演方法 | |
| Zou et al. | Deep Neural Helmholtz Operators for 3D elastic wave propagation and inversion | |
| bin Waheed et al. | A holistic approach to computing first-arrival traveltimes using neural networks | |
| CN116466402B (zh) | 一种基于地质信息和电磁数据联合驱动的电磁反演方法 | |
| CN116859478A (zh) | 一种基于瞬变电磁法成像的地下水模拟方法及系统 | |
| CN114200541B (zh) | 一种基于余弦点积梯度约束的三维重磁联合反演方法 | |
| Zhang et al. | Bayesian variational time-lapse full waveform inversion | |
| CN115563791A (zh) | 基于压缩感知重构的大地电磁数据反演方法 | |
| Yin et al. | A fast 3D gravity forward algorithm based on circular convolution | |
| CN114779356A (zh) | 一种基于阵列电阻率的地层电性剖面快速成像方法 | |
| Vatankhah et al. | Generalized L $ _p $-norm joint inversion of gravity and magnetic data using cross-gradient constraint | |
| Liu et al. | Permittivity inversion of GPR images based on DeepLab | |
| US11143788B2 (en) | Efficient solutions of inverse problems |
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 |