CN112946760B - 基于正则化方法的未爆弹三维立体成像方法、装置及系统 - Google Patents

基于正则化方法的未爆弹三维立体成像方法、装置及系统 Download PDF

Info

Publication number
CN112946760B
CN112946760B CN202110142938.7A CN202110142938A CN112946760B CN 112946760 B CN112946760 B CN 112946760B CN 202110142938 A CN202110142938 A CN 202110142938A CN 112946760 B CN112946760 B CN 112946760B
Authority
CN
China
Prior art keywords
data
matrix
expression
magnetic
abnormal
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
Application number
CN202110142938.7A
Other languages
English (en)
Other versions
CN112946760A (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.)
Tsinghua University
Original Assignee
Tsinghua 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 Tsinghua University filed Critical Tsinghua University
Priority to CN202110142938.7A priority Critical patent/CN112946760B/zh
Publication of CN112946760A publication Critical patent/CN112946760A/zh
Application granted granted Critical
Publication of CN112946760B publication Critical patent/CN112946760B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V3/00Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation
    • G01V3/08Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation operating with magnetic or electric fields produced or modified by objects or geological structures or by detecting devices
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V3/00Electric or magnetic prospecting or detecting; Measuring magnetic field characteristics of the earth, e.g. declination, deviation
    • G01V3/38Processing data, e.g. for analysis, for interpretation, for correction

Landscapes

  • Life Sciences & Earth Sciences (AREA)
  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Electromagnetism (AREA)
  • Geophysics And Detection Of Objects (AREA)
  • Measuring Magnetic Variables (AREA)

Abstract

本发明公开了一种基于正则化方法的未爆弹三维立体成像的方法、装置以及系统,涉及探测技术领域,所述方法包括:获取目标探测区域的磁场三分量数据;根据磁场三分量数据,得到梯度张量数据和磁总场梯度数据;根据梯度张量数据,得到第一先验信息;将第一先验信息加入到正则化方法中的目标函数中;根据第二先验信息和磁总场梯度数据,对目标函数的最优解进行运算,得到正则化方法的反演结果,根据反演结果,得到未爆弹的三维立体成像。本发明将第一先验信息加入到正则化方法中的目标函数中,得到的正则化方法的反演结果的精确度更高,未爆弹的三维立体成像结果也更精确,对未爆弹的识别效果得到提升,提高了遗弃化学弹发掘工作的工作效率。

Description

基于正则化方法的未爆弹三维立体成像方法、装置及系统
技术领域
本发明涉及探测技术领域,特别是涉及一种基于正则化方法的未爆弹三维立体成像的方法、装置以及系统。
背景技术
目前,在我国的城市化建设中,未爆炸弹的存在会给人民群众的身体健康和人身安全造成了严重威胁,同时也会带来巨大的生命财产损失。
为了解决上述问题,需要使用技术手段对遗弃化学弹进行探测,而遗弃化学弹的探测,可以归为地球物理勘探的范畴。地球物理勘探简称“物探”,而物探主要包括以下六种方法,分别是电磁勘探、磁法探测、重力勘探、地震勘探、地热测量法、放射性测量法。
在遗弃化学弹探测工作中,磁法探测的应用成本低,速度快,而且便于实施。因此,目前日本遗弃化学弹的探测中使用的主要方法就是磁法探测。传统的磁法探测主要用于研究地质构造或探测矿产资源,其目标的空间尺度很大,而且往往不需要得到目标精确的几何形状。但是在遗弃化学弹探测中,探测目标小,而且对目标空间位置和几何形态的要求更高。因此传统的磁法探测方法往往不能直接使用。亟需一种适应小目标,且相对更精确的成像方法,从而提高遗弃化学弹发掘工作的工作效率。
发明内容
鉴于上述问题,本发明提供一种基于正则化方法的未爆弹三维立体成像的方法、装置以及系统,提出了一种适应小目标,且相对更精确的成像方法的技术方案。
第一方面,提供一种基于正则化方法的未爆弹三维立体成像的方法,所述方法包括:
获取目标探测区域的磁场三分量数据;
根据所述磁场三分量数据,得到梯度张量数据和磁总场梯度数据;
根据所述梯度张量数据,得到第一先验信息;
将所述第一先验信息加入到正则化方法中的目标函数中;
根据第二先验信息和所述磁总场梯度数据,对所述目标函数的最优解进行运算,得到所述正则化方法的反演结果,所述第二先验信息为已知先验信息;
根据所述反演结果,得到所述未爆弹的三维立体成像。
可选地,根据所述梯度张量数据,得到第一先验信息,包括:
根据所述梯度张量数据,采用欧拉反褶积法进行运算,得到所述第一先验信息。
可选地,根据所述梯度张量数据,采用欧拉反褶积法进行运算,得到所述第一先验信息,包括:
根据所述梯度张量数据,采用欧拉反褶积法进行运算所基于的欧拉反褶积公式为:
其中x0、y0、z0为异常体中心坐标,x、y、z为观测点坐标,N为构造指数,其是异常强度随深度的衰减率,Bx、By、Bz分别是磁场在x、y、z轴方向的分量,Bxx、Bxy、Bxz、Byx、Byy、Byz、Bzx、Bzy、Bzz分别是磁梯度张量数据;
基于所述欧拉反褶积公式,进行运算过程中,选取多个相邻数据点组成一个滑动窗口,使用所述滑动窗口内的数据求解所述欧拉反褶积公式组成的线性方程组,得到异常体中心坐标;
随着所述滑动窗口的滑动,计算出多组异常体中心坐标,其中包括的异常解使用主体异常距离准则和反演结果的聚散度准则进行筛选;
去除所述异常解后,剩余正常解的分布作为所述第一先验信息。
可选地,所述多组异常体中心坐标包括:Q1,Q2,…,QN,则每一组坐标都包含x,y,z三个分量,得到如下公式:
Qi=[xi yi zi]T,i=1,2,…,N
其中,T表示转置。
可选地,所述目标函数的表达式为:
Φ=||Wd(Am-d)||2+α||Wmm||2
其中,d是一个N维矢量,表示在地表观测到的某类磁异常,包括:磁总场、磁场分量或者是其它种类的数据,m是一个M维矢量,表示地下每个长方体块的物性参数,包括:磁化强度、磁化率或者其它物性参数,A是一个N×M维的矩阵,表示观测数据d与物性参数m之间的正演算子,矩阵中每个元素的值根据正演公式进行求解,Wm是模型加权函数矩阵,用来对物性参数施加权重,Wm是一个M×M维的对角矩阵;
所述模型加权函数矩阵的表达式为:
Wm=Wm1Wm2
所述模型加权函数矩阵的表达式中的Wm1,Wm2的表达式分别如下:
其中,Wm1用于提供位置约束,Wm2用于抵消深部信号的衰减,mvnpdf()代表多维正态分布函数,μ为当前点坐标,μi为第i个异常体的位置坐标,σi为协方差矩阵,z代表当前的z轴坐标,z0和β均为经验值;
可选地,将所述第一先验信息加入到正则化方法中的目标函数中,包括:
将所述第一先验信息加入到所述模型加权函数矩阵中;
令第一表达式为:
令第二表达式为:
其中,σxσyσz分别表征了欧拉反褶积结果在x,y,z方向的不确定程度,xi yi zi分别为欧拉反褶积得到的第i组解。
将所述第一表达式作为异常体中心坐标代入Wm1表达式中,将所述第二表达式作为三个方向的方差带入Wm1表达式中,即,将所述第一先验信息加入到所述模型加权函数矩阵中。
可选地,根据第二先验信息和所述磁总场梯度数据,对所述目标函数的最优解进行运算,包括:
根据第二先验信息和所述磁总场梯度数据,利用共轭梯度法对所述目标函数的最优解进行运算。
可选地,所述反演结果为一个三维矩阵,所述三维矩阵中每个元素的值代表一个长方体单元的物性参数,根据所述反演结果,进行所述未爆弹的三维立体成像,包括:
对所述三维矩阵分别沿x、y、z轴方向切片;
对切片后的三维矩阵进行二值化处理并做三维显示,得到所述未爆弹的三维立体成像。
第二方面,提供一种基于正则化方法的未爆弹三维立体成像的装置,所述装置包括:
获取三分量模块,用于获取目标探测区域的磁场三分量数据;
张量和总场梯度模块,用于根据所述磁场三分量数据,得到梯度张量数据和磁总场梯度数据;
先验信息模块,用于根据所述梯度张量数据,得到第一先验信息;
加入模块,用于将所述第一先验信息加入到正则化方法中的目标函数中;
反演模块,用于根据第二先验信息和所述磁总场梯度数据,对所述目标函数的最优解进行运算,得到所述正则化方法的反演结果,所述第二先验信息为已知先验信息;
成像模块,用于根据所述反演结果,得到所述未爆弹的三维立体成像。
可选地,所述先验信息模块具体用于:
根据所述梯度张量数据,采用欧拉反褶积法进行运算,得到所述第一先验信息;
可选地,所述反演模块具体用于:
根据第二先验信息和所述磁总场梯度数据,利用共轭梯度法对所述目标函数的最优解进行运算。
可选地,所述反演结果为一个三维矩阵,所述三维矩阵中每个元素的值代表一个长方体单元的物性参数,所述成像模块包括:
切片单元,用于对所述三维矩阵分别沿x、y、z轴方向切片;
二值化单元,用于对切片后的三维矩阵进行二值化处理并做三维显示,得到所述未爆弹的三维立体成像。
第三方面,提供一种基于正则化方法的未爆弹三维立体成像的系统,所述系统包括:探测仪,所述探测仪用于执行上述第一方面任一所述的方法。
本发明提供的基于正则化方法的未爆弹三维立体成像的方法,获取目标探测区域的磁场三分量数据;根据磁场三分量数据,得到梯度张量数据和磁总场梯度数据;根据梯度张量数据,得到第一先验信息;将第一先验信息加入到正则化方法中的目标函数中;根据第二先验信息和磁总场梯度数据,对目标函数的最优解进行运算,得到正则化方法的反演结果,根据反演结果,得到未爆弹的三维立体成像。
本发明的方法,利用梯度张量数据消除背景场的影响,使得第一先验信息对应的边界识别效果更好。而将第一先验信息加入到正则化方法中的目标函数中,再对该目标函数的最优值进行运算,得到的正则化方法的反演结果,相较于目前的探测方法的精确度更高,因此最终得到的未爆弹的三维立体成像结果也更精确,对未爆弹的识别效果自然得到了提升,从而提高了遗弃化学弹发掘工作的工作效率。
附图说明
通过阅读下文优选实施方式的详细描述,各种其他的优点和益处对于本领域普通技术人员将变得清楚明了。附图仅用于示出优选实施方式的目的,而并不认为是对本发明的限制。而且在整个附图中,用相同的参考符号表示相同的部件。在附图中:
图1是本发明实施例一种基于正则化方法的未爆弹三维立体成像的方法的流程图;
图2是本发明实施例中一种磁异常中心水平分布俯视示意图;
图3是本发明实施例中在图2基础上,去掉不好的解后的一种磁异常中心水平分布俯视示意图;
图4是本发明实施例中在图3基础上,显示竖直位置的统计学特性的统计示意图;
图5是本发明实施例中遗弃化学弹的三维成像示意图;
图6是本发明实施例一种基于正则化方法的未爆弹三维立体成像的装置的框图。
具体实施方式
为使本发明的上述目的、特征和优点能够更加明显易懂,下面结合附图和具体实施方式对本发明作进一步详细的说明。应当理解,此处所描述的具体实施例仅用以解释本发明,仅仅是本发明一部分实施例,而不是全部的实施例,并不用于限定本发明。
参照图1,示出了本发明实施例一种基于正则化方法的未爆弹三维立体成像的方法的流程图,该方法包括:
步骤101:获取目标探测区域的磁场三分量数据。
本发明实施例中,首先需要获取目标探测区域的磁场三分量数据。一般情况下,磁探测仪对目标区域进行探测时,得到的数据即为磁场三分量数据。当然也可以通过其它方式获取到目标探测区域的磁场三分量数据,本发明实施例对此不作具体限定。
步骤102:根据磁场三分量数据,得到梯度张量数据和磁总场梯度数据。
本发明实施例中,获得磁场三分量数据后,进行运算,即可得到梯度张量数据和磁总场梯度数据。一般情况下,由于实际测量得到的数据中会同时包含背景场与磁异常,所以需要滤除数据中背景场的影响。可以对磁场三分量数据沿不同方向求梯度,得到的梯度张量数据可以消除背景场的影响。
而对磁场的三分量数据求平方和运算,可以得到磁总场模量数据,再对两层不同高度的磁总场模量数据求差值,即可得到磁总场梯度数据。
步骤103:根据梯度张量数据,得到第一先验信息。
本发明实施例中,运算得到梯度张量数据之后,再对梯度张量数据采用欧拉反褶积法进行运算,可以得到第一先验信息。欧拉反褶积法属于边界识别方法,由于梯度张量数据可以更有效地屏蔽背景场的干扰,因此可以使得边界识别效果更好,更有利于提升正则化方法中对目标函数最优解求解的精确度。
步骤104:将第一先验信息加入到正则化方法中的目标函数中。
本发明实施例中,得到第一先验信息,发明人创造性的将第一先验信息加入到正则化方法中的目标函数中,这样做,是为了提高后续正则化方法中目标函数求取最优解的精确度,具体如何将第一先验信息加入到正则化方法中的目标函数中在下文说明,先不赘述。
步骤105:根据第二先验信息和磁总场梯度数据,对目标函数的最优解进行运算,得到正则化方法的反演结果,第二先验信息为已知先验信息。
本发明实施例中,在将第一先验信息加入到正则化方法中的目标函数中后,再根据第二先验信息和前面已经得到的磁总场梯度数据,对目标函数的最优解进行运算,得到正则化方法的反演结果,所谓正则化方法的反演结果就是对目标函数的最优解进行求解而得到的结果。
所谓第二先验信息为已知先验信息。所谓已知先验信息包括:磁化强度等信息,例如:根据大量测试得到的已知未爆弹的最大磁化强度:Mmax,最小磁化强度:Mmin等。
而在正则化方法的反演中,需要提前已知异常体的磁化方向(一般以地磁场的方向近似为异常体的磁化方向),而梯度张量数据对磁化方向比较敏感,因此在正则化方法的反演中,可以使用对磁化方向不敏感的磁总场梯度数据进行反演。
本发明实施例中,根据第二先验信息和磁总场梯度数据,利用共轭梯度法对目标函数的最优解进行运算,而第二先验信息是在求解最小值的迭代过程中使用,求解时需要一步一步迭代,逐渐优化得到最优解,第二先验信息就是在这个迭代过程中添加进去的。
一种优选的迭代过程为:假设已知未爆弹的最大磁化强度为Mmax,最小磁化强度为Mmin。为了将反演出的磁化强度参数约束在最大值和最小值之间,一个简单且行之有效的方法是添加强制上下限约束,即在迭代过程中令所有大于Mmax的解等于Mmax,令所有小于Mmin的解等于Mmin。对磁化强度参数采取合适强制约束可以加速收敛的速度,使反演结果与实际模型更加一致。
步骤106:根据反演结果,得到未爆弹的三维立体成像。
本发明实施例中,得到正则化方法的反演结果之后,再次进行运算即可得到未爆弹的三维立体成像。
具体的,正则化方法的反演结果为一个三维矩阵,该三维矩阵中每个元素的值代表一个长方体单元的物性参数,对三维矩阵分别沿x、y、z轴方向切片;再对切片后的三维矩阵进行二值化处理并做三维显示,即可得到未爆弹的三维立体成像。
通过上述实施例,本发明的方法利用梯度张量数据消除背景场的影响,使得第一先验信息对应的边界识别效果更好。而将第一先验信息加入到正则化方法中的目标函数中,再对该目标函数的最优值进行运算,得到的正则化方法的反演结果,相较于目前的探测方法的精确度更高,因此最终得到的未爆弹的三维立体成像结果也更精确,对未爆弹的识别效果自然得到了提升,从而提高了遗弃化学弹发掘工作的工作效率。
以下,对上述第一先验信息如何加入目标函数中进行说明。
首先,基于梯度张量数据,采用欧拉反褶积法进行运算所基于的欧拉反褶积公式为:
其中x0、y0、z0为异常体中心坐标,x、y、z为观测点坐标,N为构造指数,其是异常强度随深度的衰减率,Bx、By、Bz分别是磁场在x、y、z轴方向的分量,Bxx、Bxy、Bxz、Byx、Byy、Byz、Bzx、Bzy、Bzz分别是磁梯度张量数据;
基于上述欧拉反褶积公式,进行运算过程中,选取多个相邻数据点组成一个滑动窗口,使用滑动窗口内的数据求解欧拉反褶积公式组成的线性方程组,得到异常体中心坐标;随着滑动窗口的滑动,计算出多组异常体中心坐标,其中包括的异常解使用主体异常距离准则和反演结果的聚散度准则进行筛选;去除异常解后,剩余正常解的分布即可作为第一先验信息。
本发明实施例中,假设多组异常体中心坐标包括:Q1,Q2,…,QN,则每一组坐标都包含x,y,z三个分量,则得到如下公式:
Qi=[xi yi zi]T,i=1,2,…,N
其中,表示转置。
本发明实施例中,正则化方法反演的目的,是对目标函数求解一个最优解,目标函数的表达式为:
Φ=||Wd(Am-d)||2+α||Wmm||2
其中,d是一个N维矢量,表示在地表观测到的某类磁异常,包括:磁总场、磁场分量或者是其它种类的数据,m是一个M维矢量,表示地下每个长方体块的物性参数,包括:磁化强度、磁化率或者其它物性参数,A是一个N×M维的矩阵,表示观测数据d与物性参数m之间的正演算子,矩阵中每个元素的值根据正演公式进行求解,Wm是模型加权函数矩阵,用来对物性参数施加权重,Wm是一个M×M维的对角矩阵;
本发明实施例中,将目标函数表达式中的模型加权函数矩阵的表达式定义为:
Wm=Wm1Wm2
则模型加权函数矩阵的表达式中的Wm1,Wm2的表达式分别如下:
其中,Wm1用于提供位置约束,Wm2用于抵消深部信号的衰减,mvnpdf()代表多维正态分布函数,μ为当前点坐标,μi为第i个异常体的位置坐标,σi为协方差矩阵,z代表当前的z轴坐标,z0和β都是经验值,例如:可以取z0=1,β=3。
本发明实施例中,将第一先验信息加入到正则化方法中的目标函数中,其具体包括:将第一先验信息加入到模型加权函数矩阵中。为了实现这个目的,可以:
令第一表达式为:
令第二表达式为:
其中,σxσyσz分别表征了欧拉反褶积结果在x,y,z方向的不确定程度,xi yi zi分别为欧拉反褶积得到的第i组解。
将第一表达式作为异常体中心坐标代入Wm1表达式中,将第二表达式作为三个方向的方差带入Wm1表达式中,通过这种方式,实现了将第一先验信息加入到模型加权函数矩阵中。而模型加权函数矩阵就是目标函数中的一个参数,因此,实质上就是将第一先验信息加入到目标函数中。而最终运算得到的目标函数的最优解,相较于目前的探测方法的精确度更高,因此最终得到的未爆弹的三维立体成像结果也更精确。
以下针对上述未爆弹三维立体成像方法的正确性、实用性通过模拟仿真进行举证。
基于电磁场仿真软件CST进行仿真,用亥姆霍兹线圈模拟地磁场,对地下40cm处的遗弃化学弹进行仿真,得到地表的磁场信息。假设遗弃化学弹长度为1.994m,弹体切面的直径为0.183m,仿真得到的初始数据为磁场三分量数据。再根据磁场三分量数据运算,得到梯度张量数据和磁总场梯度数据,对梯度张量数据采用欧拉反褶积法进行运算,并且根据主体异常准则滤除距离窗口过远的点,反演得到的磁异常中心水平分布图如图2所示,图2为俯视图,其中横坐标y代表东西方向的坐标,纵坐标x代表南北方向的坐标,图2中每一个圆圈点都代表一个解。前述已知,欧拉反褶积要通过滑动窗口的方式获得很多个解,每个解都代表对异常体位置的一个估计。得到很多个解之后,用这些解的统计学特性来表示异常体的位置的。由此可以知晓,在原点附近有一系列反演得到的中心点呈柱状排列,可以在一定程度上反映异常体的水平形态和位置。
对图2得到的结果用聚散度准则进行筛选,去除孤立存在的点,得到图3所示的示意图。图3为俯视图,其中横坐标y代表东西方向的坐标,纵坐标x代表南北方向的坐标。实质上,图3所示的解就相当于是把图2所示的解中不好的解去掉,图3可以准确反映遗弃化学弹的水平位置和形态,该图反映出第一先验信息,后续可以将第一先验信息添加到模型加权函数中。
另外,对欧拉反褶积的结果进行深度统计,得到图4所示的统计示意图,图4是为了显示竖直位置的统计学特性,图4中横坐标表示深度,纵坐标表示处在这个深度区间内的解的数量,由此可以知晓,异常体中心坐标的深度主要分布在0.4m附近,与原始模型一致。
将第一先验信息添加到模型加权函数中后,根据第二先验信息和磁总场梯度数据,对目标函数的最优解进行运算,得到正则化方法的反演结果,对反演结果进行二值化之后进行三维显示,结果如图5所示,图5中,以x、y、z三轴坐标显示遗弃化学弹的三维成像,结果表明本发明的未爆弹三维立体成像方法的正确性、实用性。
基于上述未爆弹三维立体成像的方法,本发明实施例还提出一种基于正则化方法的未爆弹三维立体成像的装置,参照图6,示出了本发明实施例一种未爆弹三维立体成像的装置的框图,所述装置包括:
获取三分量模块610,用于获取目标探测区域的磁场三分量数据;
张量和总场梯度模块620,用于根据所述磁场三分量数据,得到梯度张量数据和磁总场梯度数据;
先验信息模块630,用于根据所述梯度张量数据,得到第一先验信息;
加入模块640,用于将所述第一先验信息加入到正则化方法中的目标函数中;
反演模块650,用于根据第二先验信息和所述磁总场梯度数据,对所述目标函数的最优解进行运算,得到所述正则化方法的反演结果,所述第二先验信息为已知先验信息;
成像模块660,用于根据所述反演结果,得到所述未爆弹的三维立体成像。
可选地,所述先验信息模块630具体用于:
根据所述梯度张量数据,采用欧拉反褶积法进行运算,得到所述第一先验信息;
可选地,所述反演模块650具体用于:
根据第二先验信息和所述磁总场梯度数据,利用共轭梯度法对所述目标函数的最优解进行运算。
可选地,所述反演结果为一个三维矩阵,所述三维矩阵中每个元素的值代表一个长方体单元的物性参数,所述成像模块660包括:
切片单元,用于对所述三维矩阵分别沿x、y、z轴方向切片;
二值化单元,用于对切片后的三维矩阵进行二值化处理并做三维显示,得到所述未爆弹的三维立体成像。
基于上述未爆弹三维立体成像的方法,本发明实施例还提供一种基于正则化方法的未爆弹三维立体成像的系统,所述系统包括:探测仪,所述探测仪用于执行上述步骤101~步骤106任一所述的方法。
综上所述,本发明利用梯度张量数据消除背景场的影响,使得第一先验信息对应的边界识别效果更好。而将第一先验信息加入到正则化方法中的目标函数中,再对该目标函数的最优值进行运算,得到的正则化方法的反演结果,相较于目前的探测方法的精确度更高,因此最终得到的未爆弹的三维立体成像结果也更精确,对未爆弹的识别效果自然得到了提升,从而提高了遗弃化学弹发掘工作的工作效率。
还需要说明的是,在本文中,诸如第一和第二等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来,而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。而且,术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含,从而使得包括一系列要素的过程、方法不仅包括那些要素,而且还包括没有明确列出的其他要素,或者是还包括为这种过程、方法所固有的要素。
上面结合附图对本发明的实施例进行了描述,本文中应用了具体个例对本发明的原理及实施方式进行了阐述,以上实施例的说明只是用于帮助理解本发明的方法及其核心思想;同时,对于本领域的一般技术人员,依据本发明的思想,在具体实施方式及应用范围上均会有改变之处,综上所述,本说明书内容不应理解为对本发明的限制。

Claims (5)

1.一种基于正则化方法的未爆弹三维立体成像的方法,其特征在于,所述方法包括:
获取目标探测区域的磁场三分量数据;
根据所述磁场三分量数据,得到梯度张量数据和磁总场梯度数据;
根据所述梯度张量数据,采用欧拉反褶积法进行运算,所基于的欧拉反褶积公式为:
其中、/>、/>为异常体中心坐标,/>、/>、/>为观测点坐标,/>为构造指数,其是异常强度随深度的衰减率,Bx、By、Bz分别是磁场在x、y、z轴方向的分量,Bxx、Bxy、Bxz、Byx、Byy、Byz、Bzx、Bzy、Bzz分别是磁梯度张量数据;
基于所述欧拉反褶积公式,进行运算过程中,选取多个相邻数据点组成一个滑动窗口,使用所述滑动窗口内的数据求解所述欧拉反褶积公式组成的线性方程组,得到异常体中心坐标;
随着所述滑动窗口的滑动,计算出多组异常体中心坐标,其中包括的异常解使用主体异常距离准则和反演结果的聚散度准则进行筛选;
去除所述异常解后,剩余正常解的分布作为第一先验信息;
其中,所述多组异常体中心坐标包括:,则每一组坐标都包含/>三个分量,得到如下公式:
其中,T表示转置;
将所述第一先验信息加入到正则化方法中的目标函数中,其包括:所述目标函数的表达式为:
其中,是一个N维矢量,表示在地表观测到的某类磁异常,包括:磁总场、磁场分量或者是其它种类的数据,/>是一个M维矢量,表示地下每个长方体块的物性参数,包括:磁化强度、磁化率或者其它物性参数,/>是一个/>维的矩阵,表示观测数据/>与物性参数/>之间的正演算子,矩阵中每个元素的值根据正演公式进行求解,/>是模型加权函数矩阵,用来对物性参数施加权重,/>是一个/>维的对角矩阵;
所述模型加权函数矩阵的表达式为:
所述模型加权函数矩阵的表达式中的的表达式分别如下:
其中,用于提供位置约束,/>用于抵消深部信号的衰减,/>代表多维正态分布函数,/>为当前点坐标,/>为第/>个异常体的位置坐标,/>为协方差矩阵,z代表当前的z轴坐标,z0和β均为经验值;
将所述第一先验信息加入到所述模型加权函数矩阵中;
令第一表达式为:
令第二表达式为:
其中,σxσyσz分别表征了欧拉反褶积结果在x,y,z方向的不确定程度,xi yi zi分别为欧拉反褶积得到的第i组解;
将所述第一表达式作为异常体中心坐标代入表达式中,将所述第二表达式作为三个方向的方差带入/>表达式中,即,将所述第一先验信息加入到所述模型加权函数矩阵中;
根据第二先验信息和所述磁总场梯度数据,对所述目标函数的最优解进行运算,得到所述正则化方法的反演结果,所述第二先验信息为已知先验信息;
根据所述反演结果,得到所述未爆弹的三维立体成像。
2.根据权利要求1所述的方法,其特征在于,根据第二先验信息和所述磁总场梯度数据,对所述目标函数的最优解进行运算,包括:
根据第二先验信息和所述磁总场梯度数据,利用共轭梯度法对所述目标函数的最优解进行运算。
3.根据权利要求1所述的方法,其特征在于,所述反演结果为一个三维矩阵,所述三维矩阵中每个元素的值代表一个长方体单元的物性参数,根据所述反演结果,进行所述未爆弹的三维立体成像,包括:
对所述三维矩阵分别沿x、y、z轴方向切片;
对切片后的三维矩阵进行二值化处理并做三维显示,得到所述未爆弹的三维立体成像。
4.一种基于正则化方法的未爆弹三维立体成像的装置,其特征在于,所述装置包括:
获取三分量模块,用于获取目标探测区域的磁场三分量数据;
张量和总场梯度模块,用于根据所述磁场三分量数据,得到梯度张量数据和磁总场梯度数据;
先验信息模块,用于根据所述梯度张量数据,采用欧拉反褶积法进行运算,所基于的欧拉反褶积公式为:
其中、/>、/>为异常体中心坐标,/>、/>、/>为观测点坐标,/>为构造指数,其是异常强度随深度的衰减率,Bx、By、Bz分别是磁场在x、y、z轴方向的分量,Bxx、Bxy、Bxz、Byx、Byy、Byz、Bzx、Bzy、Bzz分别是磁梯度张量数据;
基于所述欧拉反褶积公式,进行运算过程中,选取多个相邻数据点组成一个滑动窗口,使用所述滑动窗口内的数据求解所述欧拉反褶积公式组成的线性方程组,得到异常体中心坐标;
随着所述滑动窗口的滑动,计算出多组异常体中心坐标,其中包括的异常解使用主体异常距离准则和反演结果的聚散度准则进行筛选;
去除所述异常解后,剩余正常解的分布作为第一先验信息;
其中,所述多组异常体中心坐标包括:,则每一组坐标都包含/>三个分量,得到如下公式:
其中,T表示转置;
加入模块,用于将所述第一先验信息加入到正则化方法中的目标函数中,其包括:所述目标函数的表达式为:
其中,是一个N维矢量,表示在地表观测到的某类磁异常,包括:磁总场、磁场分量或者是其它种类的数据,/>是一个M维矢量,表示地下每个长方体块的物性参数,包括:磁化强度、磁化率或者其它物性参数,/>是一个/>维的矩阵,表示观测数据/>与物性参数/>之间的正演算子,矩阵中每个元素的值根据正演公式进行求解,/>是模型加权函数矩阵,用来对物性参数施加权重,/>是一个/>维的对角矩阵;
所述模型加权函数矩阵的表达式为:
所述模型加权函数矩阵的表达式中的的表达式分别如下:
其中,用于提供位置约束,/>用于抵消深部信号的衰减,/>代表多维正态分布函数,/>为当前点坐标,/>为第/>个异常体的位置坐标,/>为协方差矩阵,z代表当前的z轴坐标,z0和β均为经验值;
将所述第一先验信息加入到所述模型加权函数矩阵中;
令第一表达式为:
令第二表达式为:
其中,σxσyσz分别表征了欧拉反褶积结果在x,y,z方向的不确定程度,xi yi zi分别为欧拉反褶积得到的第i组解;
将所述第一表达式作为异常体中心坐标代入表达式中,将所述第二表达式作为三个方向的方差带入/>表达式中,即,将所述第一先验信息加入到所述模型加权函数矩阵中;
反演模块,用于根据第二先验信息和所述磁总场梯度数据,对所述目标函数的最优解进行运算,得到所述正则化方法的反演结果,所述第二先验信息为已知先验信息;
成像模块,用于根据所述反演结果,得到所述未爆弹的三维立体成像。
5.一种基于正则化方法的未爆弹三维立体成像的系统,其特征在于,所述系统包括:探测仪,所述探测仪用于执行以下方法:
获取目标探测区域的磁场三分量数据;
根据所述磁场三分量数据,得到梯度张量数据和磁总场梯度数据;
根据所述梯度张量数据,采用欧拉反褶积法进行运算,所基于的欧拉反褶积公式为:
其中、/>、/>为异常体中心坐标,/>、/>、/>为观测点坐标,/>为构造指数,其是异常强度随深度的衰减率,Bx、By、Bz分别是磁场在x、y、z轴方向的分量,Bxx、Bxy、Bxz、Byx、Byy、Byz、Bzx、Bzy、Bzz分别是磁梯度张量数据;
基于所述欧拉反褶积公式,进行运算过程中,选取多个相邻数据点组成一个滑动窗口,使用所述滑动窗口内的数据求解所述欧拉反褶积公式组成的线性方程组,得到异常体中心坐标;
随着所述滑动窗口的滑动,计算出多组异常体中心坐标,其中包括的异常解使用主体异常距离准则和反演结果的聚散度准则进行筛选;
去除所述异常解后,剩余正常解的分布作为第一先验信息;
其中,所述多组异常体中心坐标包括:,则每一组坐标都包含/>三个分量,得到如下公式:
其中,T表示转置;
将所述第一先验信息加入到正则化方法中的目标函数中,其包括:所述目标函数的表达式为:
其中,是一个N维矢量,表示在地表观测到的某类磁异常,包括:磁总场、磁场分量或者是其它种类的数据,/>是一个M维矢量,表示地下每个长方体块的物性参数,包括:磁化强度、磁化率或者其它物性参数,/>是一个/>维的矩阵,表示观测数据/>与物性参数/>之间的正演算子,矩阵中每个元素的值根据正演公式进行求解,/>是模型加权函数矩阵,用来对物性参数施加权重,/>是一个/>维的对角矩阵;
所述模型加权函数矩阵的表达式为:
所述模型加权函数矩阵的表达式中的的表达式分别如下:
其中,用于提供位置约束,/>用于抵消深部信号的衰减,/>代表多维正态分布函数,/>为当前点坐标,/>为第/>个异常体的位置坐标,/>为协方差矩阵,z代表当前的z轴坐标,z0和β均为经验值;
将所述第一先验信息加入到所述模型加权函数矩阵中;
令第一表达式为:
令第二表达式为:
其中,σxσyσz分别表征了欧拉反褶积结果在x,y,z方向的不确定程度,xi yi zi分别为欧拉反褶积得到的第i组解;
将所述第一表达式作为异常体中心坐标代入表达式中,将所述第二表达式作为三个方向的方差带入/>表达式中,即,将所述第一先验信息加入到所述模型加权函数矩阵中;
根据第二先验信息和所述磁总场梯度数据,对所述目标函数的最优解进行运算,得到所述正则化方法的反演结果,所述第二先验信息为已知先验信息;
根据所述反演结果,得到所述未爆弹的三维立体成像。
CN202110142938.7A 2021-02-02 2021-02-02 基于正则化方法的未爆弹三维立体成像方法、装置及系统 Active CN112946760B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110142938.7A CN112946760B (zh) 2021-02-02 2021-02-02 基于正则化方法的未爆弹三维立体成像方法、装置及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110142938.7A CN112946760B (zh) 2021-02-02 2021-02-02 基于正则化方法的未爆弹三维立体成像方法、装置及系统

Publications (2)

Publication Number Publication Date
CN112946760A CN112946760A (zh) 2021-06-11
CN112946760B true CN112946760B (zh) 2024-02-06

Family

ID=76241534

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110142938.7A Active CN112946760B (zh) 2021-02-02 2021-02-02 基于正则化方法的未爆弹三维立体成像方法、装置及系统

Country Status (1)

Country Link
CN (1) CN112946760B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114518599B (zh) * 2022-01-28 2024-05-17 北京大学 一种自适应的磁异常目标成像与检测方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108227024A (zh) * 2017-12-04 2018-06-29 中国科学院地质与地球物理研究所 一种采用全张量磁梯度数据反演地下磁化率的方法及系统
CN109856689A (zh) * 2019-02-28 2019-06-07 中国科学院遥感与数字地球研究所 一种超导航磁梯度张量数据抑噪处理方法和系统
CN111190230A (zh) * 2020-01-16 2020-05-22 哈尔滨工业大学 一种基于磁梯度张量的探测方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US10242126B2 (en) * 2012-01-06 2019-03-26 Technoimaging, Llc Method of simultaneous imaging of different physical properties using joint inversion of multiple datasets
US20140129194A1 (en) * 2012-11-02 2014-05-08 Technolmaging, Llc. Methods of three-dimensional potential field modeling and inversion for layered earth models

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108227024A (zh) * 2017-12-04 2018-06-29 中国科学院地质与地球物理研究所 一种采用全张量磁梯度数据反演地下磁化率的方法及系统
CN109856689A (zh) * 2019-02-28 2019-06-07 中国科学院遥感与数字地球研究所 一种超导航磁梯度张量数据抑噪处理方法和系统
CN111190230A (zh) * 2020-01-16 2020-05-22 哈尔滨工业大学 一种基于磁梯度张量的探测方法

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
Reducing Tensor Magnetic Gradiometer Data for Unexploded Ordnance Detection;Robert E. Bracken and Philip J. Brown;Scientific Investigations Report 2005-5046;全文 *
Three-dimensional reconstruction of a small-scale magnetic target from magnetic gradient observations;Gang Yin, Lin Zhang;Journal of Magnetism and Magnetic Materials;第482卷;全文 *
利用磁法数据的解析信号探测未爆炸弹(UXO);周帅;黄大年;王泰涵;;地球物理学进展;第31卷(第04期);全文 *
基于Extrapolation Tikhonov正则化算法的重力数据及梯度多分量数据的3D反演方法研究;刘银萍;中国博士学位论文全文数据库 基础科学辑(第8期);全文 *
基于光滑约束的磁梯度张量3D正则化反演方法;李金朋;张英堂;李志宁;范红波;尹刚;;中国有色金属学报;第26卷(第12期);全文 *
基于磁梯度张量的共轭梯度3D约束反演;李金朋;张英堂;范红波;李志宁;尹刚;;探测与控制学报;第38卷(第04期);全文 *

Also Published As

Publication number Publication date
CN112946760A (zh) 2021-06-11

Similar Documents

Publication Publication Date Title
Mousavi et al. STanford EArthquake Dataset (STEAD): A global data set of seismic signals for AI
Geng et al. 3D inversion of airborne gravity-gradiometry data using cokriging
Dannemiller et al. A new method for determination of magnetization direction
Uieda et al. Robust 3D gravity gradient inversion by planting anomalous densities
Fedi et al. Understanding imaging methods for potential field data
Davis et al. Automatic detection of UXO magnetic anomalies using extended Euler deconvolution
CA2793504C (en) Windowed statistical analysis for anomaly detection in geophysical datasets
Meng et al. Microseismic monitoring of stimulating shale gas reservoir in SW China: 1. An improved matching and locating technique for downhole monitoring
CN109425906B (zh) 一种磁异常探测矢量磁目标识别方法
Pugh et al. A Bayesian method for microseismic source inversion
KR100831932B1 (ko) 오일러 해를 이용한 지하공동의 3차원 중력역산 방법 및이를 이용한 3차원 영상화 방법
Lois et al. A new automatic S-onset detection technique: Application in local earthquake data
Lagos et al. Microseismic event location using global optimization algorithms: An integrated and automated workflow
CN109407149B (zh) 基于Hessian矩阵的地震相干数据裂缝检测方法
CN110361792B (zh) 一种地球物理数据融合及成像方法、介质与设备
CN106680869A (zh) 微地震事件的检测和定位方法与装置
CN112305591B (zh) 隧道超前地质预报方法、计算机可读存储介质
Lindenbaum et al. Multiview kernels for low-dimensional modeling of seismic events
Yin et al. Automatic detection of multiple UXO-like targets using magnetic anomaly inversion and self-adaptive fuzzy c-means clustering
CN112946760B (zh) 基于正则化方法的未爆弹三维立体成像方法、装置及系统
Chambers et al. Investigation of induced microseismicity at Valhall using the Life of Field Seismic array
Lee et al. GPU‐accelerated automatic microseismic monitoring algorithm (GAMMA) and its application to the 2019 Ridgecrest earthquake sequence
Kirkendall et al. Imaging cargo containers using gravity gradiometry
CN110687598A (zh) 一种加速微震数值模拟的方法及装置
Chen et al. Joint inversion of gravity gradient tensor data based on L1 and L2 norms

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