CN109307847B - 一种磁体的二阶欧拉反演法 - Google Patents

一种磁体的二阶欧拉反演法 Download PDF

Info

Publication number
CN109307847B
CN109307847B CN201811315647.8A CN201811315647A CN109307847B CN 109307847 B CN109307847 B CN 109307847B CN 201811315647 A CN201811315647 A CN 201811315647A CN 109307847 B CN109307847 B CN 109307847B
Authority
CN
China
Prior art keywords
magnet
position coordinate
value
magnetometer
coordinate
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
CN201811315647.8A
Other languages
English (en)
Other versions
CN109307847A (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.)
Harbin Engineering University
Original Assignee
Harbin Engineering 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 Harbin Engineering University filed Critical Harbin Engineering University
Priority to CN201811315647.8A priority Critical patent/CN109307847B/zh
Publication of CN109307847A publication Critical patent/CN109307847A/zh
Application granted granted Critical
Publication of CN109307847B publication Critical patent/CN109307847B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R33/00—Arrangements or instruments for measuring magnetic variables
    • G01R33/0094—Sensor arrays
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R33/00—Arrangements or instruments for measuring magnetic variables
    • G01R33/02—Measuring direction or magnitude of magnetic fields or magnetic flux
    • G01R33/022—Measuring gradient

Landscapes

  • Physics & Mathematics (AREA)
  • Condensed Matter Physics & Semiconductors (AREA)
  • General Physics & Mathematics (AREA)
  • Measuring Magnetic Variables (AREA)

Abstract

本发明提供了一种磁体的二阶欧拉反演法，将13个标量磁强计构成标量磁强计阵列，并安装GPS或高精度惯性导航系统，确定磁体的搜索区域，在每航迹段计算磁体的位置坐标和构造指数的估计值及总估计值的均值，再计算均值的最邻近整数值，得到磁体构造指数的反演值，对更新后的第K航迹段上的磁体位置坐标估计值进行解中心计算与迭代筛选，得到磁体位置坐标的反演值。本发明便于磁强计阵列装配和平面单元替换，实现了由多个测量点处的一阶和二阶磁梯度值同时反演出磁体的位置坐标与构造指数，消减了磁梯度的粗大测量误差对反演精度的影响，也免除了构造指数的主观经验选取和梯度测量数据的中心化处理。

Description

一种磁体的二阶欧拉反演法
技术领域
本发明属于磁性目标反演与地质勘探领域，具体涉及一种磁体的二阶欧拉反演法。
背景技术
磁异常及其梯度的探测与定位在地质调查、矿产普查、地球科学研究、未爆爆炸物和军事目标的识别和定位中发挥着重要作用。高精度磁强计的涌现促进了磁探技术的革新，磁体定位技术也从传统的曲线匹配、匹配滤波等深度定位发展到欧拉反褶积、解析信号等空间定位。
欧拉反褶积法是一种磁测资料处理及解释的重要方法。D.T.Thompson提出用欧拉法，即通过磁测强度及其一阶导数进行磁体的深度定位，欧拉方程中的构造指数对应不同的磁体形状，表示测量点的磁场随磁体形状的衰减程度(EULDPH：A new technique formaking computer-assisted depth estimates from magnetic data,Geophysics,1982,47(1):31-37)；A.B.Reid等人将其推广到三维(Magnetic interpretation in threedimensions using Euler deconvolution,Geophysics,1990,55(1):80-91)。DanielaGerovska等人提出有限差分欧拉反褶积法，将每个窗口的背景场值视为常数重新建立方程，使用最小二乘法求解磁体的位置及构造指数(Finite-difference Eulerdeconvolution algorithm applied to the interpretation of magnetic data fromnorthern Bulgaria,Pure and applied geophysics,2005,162:591-608)。卞光浪等人提出在欧拉窗口内将背景磁场的影响视为线性变化，将构造指数作为动态变化参数代入欧拉方程，使用最小二乘法一并求解(应用改进欧拉算法解算磁性目标空间位置参数，测绘学报，2011，40(3):386-392)。
典型的欧拉反褶积法由于没有考虑地磁背景场因素影响而容易导致解的严重发散，克服这一问题的常用方法是将移动窗口内背景场视为常数，求解非线性方程组，但需要给定准确的构造指数线性化方程组。构造指数的选取多数依靠主观经验，当测区包含多种类型磁体时在各窗口内磁体的构造指数就不易事先给定；而错误的构造指数或构造指数动态变化时，也会导致解的发散，特别是深度计算结果更受其影响。
本发明提供一种同时测量一阶和二阶磁梯度的标量磁强计阵列，并给出磁体的二阶欧拉反演法；二阶欧拉反演法可看成欧拉反褶积法的梯度式，由磁体的一阶和二阶磁梯度测量值反演出磁体的位置坐标及构造指数。本发明所提供的标量磁强计阵列采用三个坐标平面等同的传感器安装布局，以便于磁强计阵列装配和平面单元替换。本发明所给出的二阶欧拉反演法通过对欧拉方程的梯度运算，去除了传统欧拉反褶积法中的地磁背景场及时变场；依据最邻近准则和均值统计确定磁体的构造指数的反演值，采用迭代法不断减小位置坐标解集的有效半径、建立筛选准则筛选出位置坐标的反演值，实现由多个测量点处的一阶和二阶磁梯度测量值同时反演磁体的位置坐标与构造指数，免除了磁体构造指数的主观经验选取，也消减了一阶和二阶磁梯度的粗大测量误差对磁体反演精度的影响。
发明内容
本发明的目的在于提供磁体的二阶欧拉反演法。
本发明的目的是这样实现的：
一种磁体的二阶欧拉反演法，具体的步骤为：
步骤1.将13个标量磁强计按照如下方式构成标量磁强计阵列，在三维空间坐标内，第一磁强计(0)位于坐标(0,0,0)点，即为测量点，第二磁强计(1)位于坐标(-0.3kl,-0.4kl,0)，第三磁强计(2)位于坐标(-0.5kl,-0.2kl,0)，第四磁强计(3)位于坐标(0.5kl,0.2kl,0)，第五磁强计(4)位于坐标(0.3kl,0.4kl,0)，第六磁强计(5)位于坐标(0,-0.3kl,-0.4kl)，第七磁强计(6)位于坐标(0,-0.5kl,-0.2kl)，第八磁强计(7)位于坐标(0,0.5kl,0.2kl)，第九磁强计(8)位于坐标(0,0.3kl,0.4kl)，第十磁强计(9)位于坐标(-0.3kl,0,-0.4kl)，第十一磁强计(10)位于坐标(-0.5kl,0,-0.2kl)，第十二磁强计(11)位于坐标(0.5kl,0,0.2kl)，第十三磁强计(12)位于坐标(0.3kl,0,0.4kl)，并将其捷联于运载体上，运载体安装GPS或高精度惯性导航系统；
步骤2.由磁体的磁异常大小和标量磁强计阵列敏感的最小梯度确定磁体的搜索区域，按“己”字形设置运载体的航迹；
步骤3.每隔NΣ个测量点将运载体的航迹划分为K段，在每一个航迹段由标量磁强计阵列对地磁总场的测量值及GPS或高精度惯性导航系统对运载体位置坐标的测量值，计算一次磁体的位置坐标和构造指数的估计值；
步骤4.计算磁体构造指数的所有估计值的均值，再计算均值的最邻近整数值，得到磁体构造指数的反演值；
步骤5.由磁体构造指数的反演值，建立关于磁体位置坐标的另一超定方程组，再使用最小二乘法重新计算第K航迹段上的磁体位置坐标的估计值；
步骤6.对更新后的第K航迹段上的磁体位置坐标估计值进行解中心计算与迭代筛选，得到磁体位置坐标的反演值。
所述的步骤3的计算磁体的位置坐标和构造指数的估计值的具体步骤为：
步骤3.1.同步采集标量磁强计阵列的各磁强计输出，由13个磁强计对地磁总场的测量值测得磁体的一阶和二阶磁梯度向量，计算NΣ个不同测量点处的磁体的一阶和二阶磁梯度值，测量点处磁体的一阶和二阶磁梯度向量为
Figure GDA0002660215480000031
其中
Figure GDA0002660215480000032
Figure GDA0002660215480000033
Tj(j＝0,1,…,12)为编号j的标量磁强计测得的地磁总场，符号'表示转置。
步骤3.2.同步采集运载体平台上的GPS或高精度惯性导航系统所提供的运载体平台在NΣ个测量点处的位置信息，获得第K航迹段上的第i(i＝1,2,…,NΣ)个测量点处的运载体位置坐标值；
步骤3.3.由NΣ个磁体的一阶和二阶磁梯度值及相应的运载体位置坐标值，根据最小二乘法计算第K航迹段上磁体的位置坐标和构造指数的估计值。
所述的步骤6磁体位置坐标反演的具体步骤为：
步骤6.1.迭代次数n＝0，设定最大迭代次数nmax、位置坐标解空间的最大半径rEuler、位置坐标解中心间距的门限值εcenter；
步骤6.2.计算位置坐标解的均值，得到位置坐标解的中心；
步骤6.3.根据rEuler设定位置坐标解的筛选准则，从位置坐标解中挑选出满足筛选准则的解组成新的位置坐标解集；
步骤6.4.迭代次数增1，如果n＞nmax，则转到步骤6；否则，对新的位置坐标解集重新计算其均值，得到新的位置坐标解中心；
步骤6.5.计算更新前与更新后的位置坐标解中心的差值的距离，如果差值的距离小于εcenter，则转到步骤6，否则，转到步骤3；
步骤6.6.迭代结束，输出磁体位置坐标的反演值。
本发明的有益效果在于：提出的一种标量磁强计阵列在三个坐标平面内采用相同的传感器布局，具有便于磁强计阵列装配和平面单元替换的优点，也是一种全新的能同时测量一阶和二阶磁梯度的标量磁强计阵列；提出的磁体的二阶欧拉反演法通过欧拉方程的梯度运算，建立了一种梯度形式的欧拉反褶积方程，去除了地磁背景场及时变场的影响；依据最邻近准则和均值统计确定磁体的构造指数的反演值，采用迭代法不断减小位置坐标解集的有效半径、设立筛选准则筛选出位置坐标的反演值，实现了由多个测量点处的一阶和二阶磁梯度值同时反演出磁体的位置坐标与构造指数，消减了磁梯度的粗大测量误差对反演精度的影响，也免除了构造指数的主观经验选取和梯度测量数据的中心化处理；与现有的欧拉反褶积反演法相比，本发明所提供的磁体的二阶欧拉反演法无需考虑地磁背景场及时变场的影响，也不必依靠主观经验选取构造指数和中心化处理梯度测量数据；能由多个测量点处的一阶和二阶磁梯度值同时反演出磁体的位置坐标与构造指数，采用迭代筛选法消减了磁梯度的粗大测量误差对反演精度的影响。
附图说明
图1为同时测量一阶和二阶磁梯度的标量磁强计阵列的传感器配置图。
图2为“己”字形运载体的航迹示意图。
图3为磁体的二阶欧拉反演法的方法流程图。
图4为球体磁体的一阶和二阶磁梯度的测量值与理论值曲线。
图5为球体磁体的位置坐标计算值曲线。
图6为无限长圆柱体磁体的位置坐标计算值曲线。
具体实施方式
下面结合附图对本发明做进一步的描述：
实施例1
一种同时测量一阶和二阶磁梯度的标量磁强计阵列，其结构由13个相同的标量磁强计组成，在三维空间坐标内，第一磁强计(0)位于坐标(0,0,0)点，即为测量点，第二磁强计(1)位于坐标(-0.3kl,-0.4kl,0)，第三磁强计(2)位于坐标(-0.5kl,-0.2kl,0)，第四磁强计(3)位于坐标(0.5kl,0.2kl,0)，第五磁强计(4)位于坐标(0.3kl,0.4kl,0)，第六磁强计(5)位于坐标(0,-0.3kl,-0.4kl)，第七磁强计(6)位于坐标(0,-0.5kl,-0.2kl)，第八磁强计(7)位于坐标(0,0.5kl,0.2kl)，第九磁强计(8)位于坐标(0,0.3kl,0.4kl)，第十磁强计(9)位于坐标(-0.3kl,0,-0.4kl)，第十一磁强计(10)位于坐标(-0.5kl,0,-0.2kl)，第十二磁强计(11)位于坐标(0.5kl,0,0.2kl)，第十三磁强计(12)位于坐标(0.3kl,0,0.4kl)。
一阶磁梯度的表达式为
Figure GDA0002660215480000051
式中，ξ＝x,y,z，δT为磁体的磁异常大小。
二阶磁梯度的表达式为
Figure GDA0002660215480000052
式中，ξ,η＝x,y,z。
由一阶和二阶磁梯度所构成的向量为
XδT＝[δTx δTy δTz δTxx δTyy δTzz δTxy δTyz δTxz]′，其中符号'表示转置。由于标量磁强计阵列所在区域内的地磁背景场及其时变场的梯度远小于磁体的磁梯度，且磁体的三阶和三阶以上的磁梯度均远小于其一阶和二阶梯度，因此有
AδTXδT＝bδT
测量点处磁体的一阶和二阶磁梯度向量为
Figure GDA0002660215480000061
其中
Figure GDA0002660215480000062
Figure GDA0002660215480000063
Tj(j＝0,1,…,12)为编号j的标量磁强计测得的地磁总场，符号'表示转置。
一种磁体的二阶欧拉反演法，具体的步骤为：
步骤1.将13个标量磁强计按照权利要求1所述的方式构成标量磁强计阵列，并将其捷联于运载体上，运载体安装GPS或高精度惯性导航系统；
步骤2.由磁体的磁异常大小和标量磁强计阵列敏感的最小梯度确定磁体的搜索区域，按“己”字形设置运载体的航迹；
步骤3.每隔NΣ个测量点将运载体的航迹划分为K段，在每一个航迹段由标量磁强计阵列对地磁总场的测量值及GPS或高精度惯性导航系统对运载体位置坐标的测量值，计算一次磁体的位置坐标和构造指数的估计值；
步骤4.计算磁体构造指数的所有估计值的均值，再计算均值的最邻近整数值，得到磁体构造指数的反演值；
Figure GDA0002660215480000071
式中，round(·)表示求最邻近的整数
步骤5.由磁体构造指数的反演值，建立关于磁体位置坐标的另一超定方程组，再使用最小二乘法重新计算第K航迹段上的磁体位置坐标的估计值；
Figure GDA0002660215480000072
式中，
Figure GDA0002660215480000073
Figure GDA0002660215480000074
和
Figure GDA0002660215480000075
为更新后的第k航迹段上磁体的位置坐标估计值；矩阵D3T(k)的表达式为
Figure GDA0002660215480000076
向量y3T的表达式为
Figure GDA0002660215480000077
步骤6.对更新后的第K航迹段上的磁体位置坐标估计值进行解中心计算与迭代筛选，得到磁体位置坐标的反演值。
所述的步骤3的计算磁体的位置坐标和构造指数的估计值的具体步骤为：
步骤3.1.同步采集标量磁强计阵列的各磁强计输出，由13个磁强计对地磁总场的测量值测得测量点处磁体的一阶和二阶磁梯度向量
Figure GDA0002660215480000081
通过下面的公式计算NΣ个不同测量点处的磁体的一阶和二阶磁梯度值；
Figure GDA0002660215480000082
步骤3.2.同步采集运载体平台上的GPS或高精度惯性导航系统所提供的运载体平台在NΣ个测量点处的位置信息，获得第K航迹段上的第i(i＝1,2,…,NΣ)个测量点处的运载体位置坐标值x(k,i)、y(k,i)和z(k,i)；
步骤3.3.由NΣ个磁体的一阶和二阶磁梯度值及相应的运载体位置坐标值，根据最小二乘法计算第K航迹段上磁体的位置坐标和构造指数的估计值
Figure GDA0002660215480000083
式中，
Figure GDA0002660215480000084
Figure GDA0002660215480000085
Figure GDA0002660215480000086
和
Figure GDA0002660215480000087
为第k航迹段上磁体位置坐标的估计值，
Figure GDA0002660215480000088
为第k航迹段上磁体构造指数的估计值；矩阵D4T(k)的表达式为
Figure GDA0002660215480000089
式中，δTx(k,i)为第k航迹段上的第i个测量点处的一阶梯度δTx，i＝1,2,…,NΣ，其他的物理量含义与之类似；向量y4T的表达式为
Figure GDA0002660215480000091
所述的步骤6磁体位置坐标反演的具体步骤为：
步骤6.1.迭代次数n＝0，设定最大迭代次数nmax、位置坐标解空间的最大半径rEuler、位置坐标解中心间距的门限值εcenter；
步骤6.2.计算位置坐标解的均值，得到位置坐标解的中心；
步骤6.3.根据rEuler设定位置坐标解的筛选准则，从位置坐标解中挑选出满足筛选准则的解组成新的位置坐标解集；
步骤6.4.迭代次数增1，如果n＞nmax，则转到步骤6；否则，对新的位置坐标解集重新计算其均值，得到新的位置坐标解中心；
步骤6.5.计算更新前与更新后的位置坐标解中心的差值的距离，如果差值的距离小于εcenter，则转到步骤6，否则，转到步骤3；
步骤6.6.迭代结束，输出磁体位置坐标的反演值。
本发明提供了磁体的二阶欧拉反演法及一种标量磁强计阵列。提供的标量磁强计阵列在三个坐标平面内采用相同的传感器安装布局，这样便于磁强计阵列装配和平面单元替换。提供的磁体的二阶欧拉反演法通过欧拉方程的梯度运算来建立梯度形式的欧拉反褶积方程，去除了传统欧拉反褶积法中的地磁背景场及时变场；依据最邻近准则和均值统计确定磁体的构造指数的反演值，采用迭代法不断减小位置坐标解集的有效半径、建立筛选准则筛选出位置坐标的反演值，实现了由多点处磁体的一阶和二阶磁梯度值同时反演出磁体的位置坐标与构造指数，消减了一阶和二阶磁梯度的粗大测量误差对反演精度的影响，也免除了磁体构造指数的主观经验选取，提高了二阶欧拉反演法对磁体位置坐标与构造指数的反演精度。
实施例2
磁体的二阶欧拉反演法及一种标量磁强计阵列
本发明属于磁性目标反演与地质勘探领域，具体涉及到磁体的二阶欧拉反演法及一种标量磁强计阵列。
磁异常及其梯度的探测与定位在地质调查、矿产普查、地球科学研究、未爆爆炸物和军事目标的识别和定位中发挥着重要作用。高精度磁强计的涌现促进了磁探技术的革新，磁体定位技术也从传统的曲线匹配、匹配滤波等深度定位发展到欧拉反褶积、解析信号等空间定位。
欧拉反褶积法是一种磁测资料处理及解释的重要方法。D.T.Thompson提出用欧拉法，即通过磁测强度及其一阶导数进行磁体的深度定位，欧拉方程中的构造指数对应不同的磁体形状，表示测量点的磁场随磁体形状的衰减程度(EULDPH：A new technique formaking computer-assisted depth estimates from magnetic data,Geophysics,1982,47(1):31-37)；A.B.Reid等人将其推广到三维(Magnetic interpretation in threedimensions using Euler deconvolution,Geophysics,1990,55(1):80-91)。DanielaGerovska等人提出有限差分欧拉反褶积法，将每个窗口的背景场值视为常数重新建立方程，使用最小二乘法求解磁体的位置及构造指数(Finite-difference Eulerdeconvolution algorithm applied to the interpretation of magnetic data fromnorthern Bulgaria,Pure and applied geophysics,2005,162:591-608)。卞光浪等人提出在欧拉窗口内将背景磁场的影响视为线性变化，将构造指数作为动态变化参数代入欧拉方程，使用最小二乘法一并求解(应用改进欧拉算法解算磁性目标空间位置参数，测绘学报，2011，40(3):386-392)。
典型的欧拉反褶积法由于没有考虑地磁背景场因素影响而容易导致解的严重发散，克服这一问题的常用方法是将移动窗口内背景场视为常数，求解非线性方程组，但需要给定准确的构造指数线性化方程组。构造指数的选取多数依靠主观经验，当测区包含多种类型磁体时在各窗口内磁体的构造指数就不易事先给定；而错误的构造指数或构造指数动态变化时，也会导致解的发散，特别是深度计算结果更受其影响。
本发明提供一种同时测量一阶和二阶磁梯度的标量磁强计阵列，并给出磁体的二阶欧拉反演法；二阶欧拉反演法可看成欧拉反褶积法的梯度式，由磁体的一阶和二阶磁梯度测量值反演出磁体的位置坐标及构造指数。本发明所提供的标量磁强计阵列采用三个坐标平面等同的传感器安装布局，以便于磁强计阵列装配和平面单元替换。本发明所给出的二阶欧拉反演法通过对欧拉方程的梯度运算，去除了传统欧拉反褶积法中的地磁背景场及时变场；依据最邻近准则和均值统计确定磁体的构造指数的反演值，采用迭代法不断减小位置坐标解集的有效半径、建立筛选准则筛选出位置坐标的反演值，实现由多个测量点处的一阶和二阶磁梯度测量值同时反演磁体的位置坐标与构造指数，免除了磁体构造指数的主观经验选取，也消减了一阶和二阶磁梯度的粗大测量误差对磁体反演精度的影响。
本发明的目的在于提供磁体的二阶欧拉反演法及一种标量磁强计阵列。
本发明的目的是这样实现的：本发明的一种同时测量一阶和二阶磁梯度的标量磁强计阵列，采用三个坐标平面全同的传感器安装布局，也就是说，三个坐标平面上的标量磁强计配置完全相同，如图1所示。编号为0的磁强计位于阵列的中心点，即为测量点；编号为1、2、3和4的磁强计位于xy平面内，相对于测量点的位置坐标分别为(-0.3kl,-0.4kl,0)、(-0.5kl,-0.2kl,0)、(0.5kl,0.2kl,0)和(0.3kl,0.4kl,0)；共有2对空间位置关于测量点对称的磁强计，分别为1-4和2-3。编号为5、6、7和8的磁强计位于yz平面内，相对于测量点的位置坐标分别为(0,-0.3kl,-0.4kl)，(0,-0.5kl,-0.2kl)，(0,0.5kl,0.2kl)和(0,0.3kl,0.4kl)；共有2对空间位置关于测量点对称的磁强计，分别为5-8和6-7。编号为9、10、11和12的磁强计位于xz平面内，相对于测量点的位置坐标分别为(-0.3kl,0,-0.4kl)，(-0.5kl,0,-0.2kl)，(0.5kl,0,0.2kl)和(0.3kl,0,0.4kl)；共有2对空间位置关于测量点对称的磁强计，分别为9-12和10-11。
一阶磁梯度的表达式为
Figure GDA0002660215480000111
式中，ξ＝x,y,z，δT为磁体的磁异常大小。
二阶磁梯度的表达式为
Figure GDA0002660215480000112
式中，ξ,η＝x,y,z。
由一阶和二阶磁梯度所构成的向量为
XδT＝[δTx δTy δTz δTxx δTyy δTzz δTxy δTyz δTxz]′，其中符号'表示转置。由于标量磁强计阵列所在区域内的地磁背景场及其时变场的梯度远小于磁体的磁梯度，且磁体的三阶和三阶以上的磁梯度均远小于其一阶和二阶梯度，因此有
AδTXδT＝bδT (3)
式中，
Figure GDA0002660215480000121
Figure GDA0002660215480000122
式中，Tj(j＝0,1,…,12)为编号j的标量磁强计测得的地磁总场。
使用最小二乘法求解式(3)所示的线性方程组得，标量磁强计阵列测得的测量点处磁体的一阶和二阶磁梯度向量为
Figure GDA0002660215480000123
标量磁强计阵列的构建及磁体的二阶欧拉反演法的具体步骤如下：
步骤1、将13个标量磁强计按图1所示的配置方式构成标量磁强计阵列，并将其捷联于运载体上，运载体同时安装有GPS或高精度惯性导航系统。
步骤2、由磁体的磁异常大小和标量磁强计阵列敏感的最小梯度确定磁体的搜索区域，按“己”字形设置运载体的航迹。
步骤3、每隔NΣ个测量点将运载体的航迹划分为K段，在每一个航迹段由标量磁强计阵列对地磁总场的测量值及GPS或高精度惯性导航系统对运载体位置坐标的测量值，计算一次磁体的位置坐标和构造指数的估计值。
步1)同步采集标量磁强计阵列的各磁强计输出，由13个磁强计对地磁总场的测量值按式(6)计算NΣ个不同测量点处的磁体的一阶和二阶磁梯度值；
步2)同步采集运载体平台上的GPS或高精度惯性导航系统所提供的运载体平台在NΣ个测量点处的位置信息，获得第k航迹段上的第i(i＝1,2,…,NΣ)个测量点处的运载体位置坐标值；
步3)由NΣ个磁体的一阶和二阶磁梯度值及相应的运载体位置坐标值，根据最小二乘法计算第k航迹段上磁体的位置坐标和构造指数的估计值。
步骤4、先计算磁体构造指数的所有估计值的均值，再计算均值的最邻近整数值，得到磁体构造指数的反演值。
步骤5、由磁体构造指数的反演值，建立关于磁体位置坐标的另一超定方程组，再使用最小二乘法重新计算第k航迹段上的磁体位置坐标的估计值。
步骤6、对更新后的第k航迹段上的磁体位置坐标估计值进行解中心计算与迭代筛选，得到磁体位置坐标的反演值。
步1)迭代次数n＝0，设定最大迭代次数nmax、位置坐标解空间的最大半径rEuler、位置坐标解中心间距的门限值εcenter。
步2)计算位置坐标解的均值，得到位置坐标解的中心。
步3)根据rEuler设定位置坐标解的筛选准则，从位置坐标解中挑选出满足筛选准则的解组成新的位置坐标解集。
步4)迭代次数增1，如果n＞nmax，则转步6)；否则，对新的位置坐标解集重新计算其均值，得到新的位置坐标解中心。
步5)计算更新前与更新后的位置坐标解中心的差值的距离，如果差值的距离小于εcenter，则转步6)；否则，转步3)。
步6)迭代结束，输出磁体位置坐标的反演值。
本发明提供了磁体的二阶欧拉反演法及一种标量磁强计阵列。提供的标量磁强计阵列在三个坐标平面内采用相同的传感器安装布局，这样便于磁强计阵列装配和平面单元替换。提供的磁体的二阶欧拉反演法通过欧拉方程的梯度运算来建立梯度形式的欧拉反褶积方程，去除了传统欧拉反褶积法中的地磁背景场及时变场；依据最邻近准则和均值统计确定磁体的构造指数的反演值，采用迭代法不断减小位置坐标解集的有效半径、建立筛选准则筛选出位置坐标的反演值，实现了由多点处磁体的一阶和二阶磁梯度值同时反演出磁体的位置坐标与构造指数，消减了一阶和二阶磁梯度的粗大测量误差对反演精度的影响，也免除了磁体构造指数的主观经验选取，提高了二阶欧拉反演法对磁体位置坐标与构造指数的反演精度。
本发明与现有技术比较具有以下优点：提出的一种标量磁强计阵列在三个坐标平面内采用相同的传感器布局，具有便于磁强计阵列装配和平面单元替换的优点，也是一种全新的能同时测量一阶和二阶磁梯度的标量磁强计阵列。提出的磁体的二阶欧拉反演法通过欧拉方程的梯度运算，建立了一种梯度形式的欧拉反褶积方程，去除了地磁背景场及时变场的影响；依据最邻近准则和均值统计确定磁体的构造指数的反演值，采用迭代法不断减小位置坐标解集的有效半径、设立筛选准则筛选出位置坐标的反演值，实现了由多个测量点处的一阶和二阶磁梯度值同时反演出磁体的位置坐标与构造指数，消减了磁梯度的粗大测量误差对反演精度的影响，也免除了构造指数的主观经验选取和梯度测量数据的中心化处理。与现有的欧拉反褶积反演法相比，本发明所提供的磁体的二阶欧拉反演法无需考虑地磁背景场及时变场的影响，也不必依靠主观经验选取构造指数和中心化处理梯度测量数据；能由多个测量点处的一阶和二阶磁梯度值同时反演出磁体的位置坐标与构造指数，采用迭代筛选法消减了磁梯度的粗大测量误差对反演精度的影响。
下面结合附图对本发明的实施方式进行详细描述：
步骤1、将编号从1到4、从5到8、从9到12的同类型标量磁强计均按图1各自集成在一个平面板上，构成三块相同的平面板块；再将这三块相同的板块与编号为0的标量磁强计装配在一起构成标量磁强计阵列，共13个标量磁强计；然后，使用无磁连接件将标量磁强计阵列捷联于运载体，运载体同时安装有GPS或高精度惯性导航系统。
步骤2、根据磁体磁异常大小和标量磁强计阵列敏感的最小梯度确定磁体的搜索区域，按“己”字形设置运载体的航迹。
步骤3、每隔NΣ个测量点将运载体的航迹划分为K段，在每一个航迹段由标量磁强计阵列对地磁总场的测量值及GPS或高精度惯性导航系统对运载体位置坐标的测量值，计算一次磁体的位置坐标和构造指数估计值。
步1)多通道数据采集器同步采集标量磁强计阵列的各磁强计输出，由13个磁强计对地磁总场的测量值按式(6)计算不同测量点处的
Figure GDA0002660215480000141
按式(7)获得不同测量点处的一阶和二阶磁梯度值。
Figure GDA0002660215480000151
步2)多通道数据采集器同步采集运载体上的GPS或高精度惯性导航系统所提供的运载体在NΣ个测量点处位置信息，得到第k航迹段上的第i个测量点处的运载体位置坐标x(k,i)、y(k,i)和z(k,i)。
步3)由NΣ个磁体的一阶和二阶磁梯度值以及相应的运载体位置坐标值，根据最小二乘法按式(8)计算第k航迹段上的磁体位置坐标及构造指数估计值。
Figure GDA0002660215480000152
式中，
Figure GDA0002660215480000153
Figure GDA0002660215480000154
Figure GDA0002660215480000155
和
Figure GDA0002660215480000156
为第k航迹段上磁体位置坐标的估计值，
Figure GDA0002660215480000157
为第k航迹段上磁体构造指数的估计值；矩阵D4T(k)的表达式为
Figure GDA0002660215480000158
式中，δTx(k,i)为第k航迹段上的第i个测量点处的一阶梯度δTx，i＝1,2,…,NΣ，其他的物理量含义与之类似；向量y4T的表达式为
Figure GDA0002660215480000161
步骤4、由步3)得到的磁体构造指数的所有估计值，根据式(11)计算磁体构造指数的估计值
Figure GDA0002660215480000162
作为磁体构造指数的反演值。
Figure GDA0002660215480000163
式中，round(·)表示求最邻近的整数。
步骤5、由磁体构造指数的反演值
Figure GDA0002660215480000164
按式(12)重新计算第k航迹段上磁体的位置坐标估计值。
Figure GDA0002660215480000165
式中，
Figure GDA0002660215480000166
Figure GDA0002660215480000167
和
Figure GDA0002660215480000168
为更新后的第k航迹段上磁体的位置坐标估计值；矩阵D3T(k)的表达式为
Figure GDA0002660215480000169
向量y3T的表达式为
Figure GDA0002660215480000171
步骤6、对更新后的第k航迹段上磁体的位置坐标估计值进行解中心计算与迭代筛选，获得磁体位置坐标的反演值。
步1)迭代次数n＝0，设定最大迭代次数nmax、位置坐标解空间最大半径rEuler、位置坐标解中心间距门限值εcenter。
步2)按式(15)计算位置坐标解中心X3c(n)，
Figure GDA0002660215480000172
式中，Pn为位置坐标解的个数。
步3)按式(16)所示的准则筛选位置坐标解，
Figure GDA0002660215480000173
即从位置坐标解挑选出满足式(16)的解组成新的位置坐标解集。
步4)迭代次数增1，如果n＞nmax，则转步6)；否则，对新的位置坐标解集按式(15)重新计算位置坐标解中心X3c(n)；
步5)计算||X3c(n)-X3c(n-1)||，如果||X3c(n)-X3c(n-1)||＜εcenter，则转步6)；否则，转步3)。
步6)结束迭代，计算位置坐标解集的均值得到磁体位置坐标的反演值。
磁体的二阶欧拉反演法的方法流程图如图2所示。
为直接地反映磁体的二阶欧拉反演法解算目标位置坐标的效果，定义磁体位置坐标的绝对反演误差，如式(17)所示。
Figure GDA0002660215480000181
式中，(x0,y0,z0)为磁体位置坐标的真值，
Figure GDA0002660215480000182
为磁体位置坐标的反演值。
仿真实验一、当地地磁场的磁倾角及磁偏角分别为45°和30°，球体磁体的半径为20米，其球心位置坐标真值为(-10米,10米,40米)，磁化强度为1.6×103安培/米，标量磁强计阵列的kl＝1米，每一个标量磁强计的测量噪声为零均值、标准差为1pT的正态分布白噪声；运载体在平面z＝0内按“己”字形路线航行，x和y方向的搜索范围分别为(-100米,100米)和(-80米,80米)，x和y方向的步进长度分别为5米和4米；最大迭代次数为1000，位置坐标解空间最大半径为10，位置坐标解中心间距门限值为1。
球体磁体的一阶和二阶磁梯度在运载体航迹上各测量点处的测量值与理论值的曲线如图4所示，其中图4(a)为δTx的测量值与理论值的曲线，图4(b)为δTy的测量值与理论值的曲线，图4(c)为δTz的测量值与理论值的曲线，图4(d)为δTxx的测量值与理论值的曲线，图4(e)为δTyy的测量值与理论值的曲线，图4(f)为δTzz的测量值与理论值的曲线，图4(g)为δTxy的测量值与理论值的曲线，图4(h)为δTyz的测量值与理论值的曲线，图4(i)为δTxz的测量值与理论值的曲线，黑色方形点线代表测量值曲线，红色圆形点线代表理论值曲线。
经迭代筛选后，得到球体磁体的位置坐标解的曲线如图5所示；对图5中的磁体位置坐标解进行均值统计，得到x、y和z方向上磁体位置坐标的反演值分别为-9.96622米、10.05941米和40.00728米，标准差分别为1.38819米、1.50928米和1.41436米；因此，位置坐标的绝对测量误差为0.06873米。反演出的球体磁体的构造指数为3，与其真值一致。
仿真实验二、当地地磁场的磁倾角及磁偏角分别为45°和30°，沿y方向无限伸长的水平圆柱体磁体的半径为1米，其中心位置坐标真值为(10米,28米)，磁化强度为3×104安培/米，标量磁强计阵列的kl＝1米，每一个标量磁强计的测量噪声为零均值、标准差为1pT的正态分布白噪声；运载体在平面z＝0内按“己”字形路线航行，x和y方向的搜索范围分别为(-100米,100米)和(-80米,80米)，x和y方向的步进长度分别为5米和4米；最大迭代次数为1000，位置坐标解空间最大半径为1，位置坐标解中心间距门限值为0.2。
经迭代筛选后，得到无限伸长的水平圆柱体磁体的位置坐标解的曲线如图6所示；对图6中的磁体位置坐标解进行均值统计，得到x和z方向上磁体位置坐标的反演值分别为9.99775米和28.02016米，标准差分别为0.22927米和0.21808米；因此，位置坐标的绝对测量误差为0.02029米。无限伸长的水平圆柱体磁体的构造指数的反演值为2，与其真值一致。
本发明有益效果说明如下：
进行了球体磁体的位置坐标与构造指数反演法的仿真实验，给出了标量磁强计阵列对一阶和二阶磁梯度的测量结果曲线、位置坐标的反演值曲线及构造指数的反演值；进行了无限伸长的水平圆柱体磁体的位置坐标与构造指数反演法的仿真实验，给出了其位置坐标的反演值曲线及构造指数的反演值；结果都表明本发明所提的标量磁强计阵列可用于测量一阶和二阶磁梯度，相应的二阶欧拉反演法能同时反演出磁体的位置坐标与构造指数。
本发明所提供的一种全新的标量磁强计阵列采用三个坐标平面相同的传感器安装布局，不仅能实现对一阶和二阶磁梯度的同时测量，而且具有便于磁强计阵列装配和平面单元替换的优点。与现有的欧拉反褶积定位法相比，本发明所提供的磁体的二阶欧拉反演法通过欧拉方程的梯度运算，去除了地磁背景场及时变场的影响；依据最邻近准则和均值统计确定磁体的构造指数的反演值，采用迭代法不断减小位置坐标解集的有效半径、建立筛选准则筛选出磁体的位置坐标的反演值，实现了由多个测量点处的一阶和二阶磁梯度值同时反演出磁体的位置坐标与构造指数；二阶欧拉反演法不必依靠主观经验选取构造指数，也不必中心化处理梯度测量数据，消减了一阶和二阶梯度的粗大测量误差对磁体位置坐标与构造指数的反演精度的影响，提高了磁体的二阶欧拉反演法的反演精度。

Claims (3)

1.一种磁体的二阶欧拉反演法，其特征在于，具体的步骤为：
步骤1.将13个标量磁强计按照如下方式构成标量磁强计阵列，在三维空间坐标内，第一磁强计(0)位于坐标(0,0,0)点，即为测量点，第二磁强计(1)位于坐标(-0.3kl,-0.4kl,0)，第三磁强计(2)位于坐标(-0.5kl,-0.2kl,0)，第四磁强计(3)位于坐标(0.5kl,0.2kl,0)，第五磁强计(4)位于坐标(0.3kl,0.4kl,0)，第六磁强计(5)位于坐标(0,-0.3kl,-0.4kl)，第七磁强计(6)位于坐标(0,-0.5kl,-0.2kl)，第八磁强计(7)位于坐标(0,0.5kl,0.2kl)，第九磁强计(8)位于坐标(0,0.3kl,0.4kl)，第十磁强计(9)位于坐标(-0.3kl,0,-0.4kl)，第十一磁强计(10)位于坐标(-0.5kl,0,-0.2kl)，第十二磁强计(11)位于坐标(0.5kl,0,0.2kl)，第十三磁强计(12)位于坐标(0.3kl,0,0.4kl)，并将其捷联于运载体上，运载体安装GPS或高精度惯性导航系统；
步骤2.由磁体的磁异常大小和标量磁强计阵列敏感的最小梯度确定磁体的搜索区域，按“己”字形设置运载体的航迹；
步骤3.每隔NΣ个测量点将运载体的航迹划分为K段，在每一个航迹段由标量磁强计阵列对地磁总场的测量值及GPS或高精度惯性导航系统对运载体位置坐标的测量值，计算一次磁体的位置坐标和构造指数的估计值；
步骤4.计算磁体构造指数的所有估计值的均值，再计算均值的最邻近整数值，得到磁体构造指数的反演值；
步骤5.由磁体构造指数的反演值，建立关于磁体位置坐标的另一超定方程组，再使用最小二乘法重新计算第K航迹段上的磁体位置坐标的估计值；
步骤6.对更新后的第K航迹段上的磁体位置坐标估计值进行解中心计算与迭代筛选，得到磁体位置坐标的反演值。
2.根据权利要求1所述的一种磁体的二阶欧拉反演法，其特征在于，所述的步骤3的计算磁体的位置坐标和构造指数的估计值的具体步骤为：
步骤3.1.同步采集标量磁强计阵列的各磁强计输出，由13个磁强计对地磁总场的测量值测得磁体的一阶和二阶磁梯度向量，计算NΣ个不同测量点处的磁体的一阶和二阶磁梯度值，测量点处磁体的一阶和二阶磁梯度向量为
Figure FDA0002660215470000011
其中
Figure FDA0002660215470000021
Figure FDA0002660215470000022
Tj(j＝0,1,…,12)为编号j的标量磁强计测得的地磁总场，符号'表示转置；
步骤3.2.同步采集运载体平台上的GPS或高精度惯性导航系统所提供的运载体平台在NΣ个测量点处的位置信息，获得第K航迹段上的第i(i＝1,2,…,NΣ)个测量点处的运载体位置坐标值；
步骤3.3.由NΣ个磁体的一阶和二阶磁梯度值及相应的运载体位置坐标值，根据最小二乘法计算第K航迹段上磁体的位置坐标和构造指数的估计值。
3.根据权利要求1所述的一种磁体的二阶欧拉反演法，其特征在于，所述的步骤6磁体位置坐标反演的具体步骤为：
步骤6.1.迭代次数n＝0，设定最大迭代次数nmax、位置坐标解空间的最大半径rEuler、位置坐标解中心间距的门限值εcenter；
步骤6.2.计算位置坐标解的均值，得到位置坐标解的中心；
步骤6.3.根据rEuler设定位置坐标解的筛选准则，从位置坐标解中挑选出满足筛选准则的解组成新的位置坐标解集；
步骤6.4.迭代次数增1，如果n＞nmax，则转到步骤6；否则，对新的位置坐标解集重新计算其均值，得到新的位置坐标解中心；
步骤6.5.计算更新前与更新后的位置坐标解中心的差值的距离，如果差值的距离小于εcenter，则转到步骤6，否则，转到步骤3；
步骤6.6.迭代结束，输出磁体位置坐标的反演值。
CN201811315647.8A 2018-11-06 2018-11-06 一种磁体的二阶欧拉反演法 Active CN109307847B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201811315647.8A CN109307847B (zh) 2018-11-06 2018-11-06 一种磁体的二阶欧拉反演法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201811315647.8A CN109307847B (zh) 2018-11-06 2018-11-06 一种磁体的二阶欧拉反演法

Publications (2)

Publication Number Publication Date
CN109307847A CN109307847A (zh) 2019-02-05
CN109307847B true CN109307847B (zh) 2021-01-05

Family

ID=65222996

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201811315647.8A Active CN109307847B (zh) 2018-11-06 2018-11-06 一种磁体的二阶欧拉反演法

Country Status (1)

Country Link
CN (1) CN109307847B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111766549B (zh) * 2020-07-07 2023-03-31 北京卫星环境工程研究所 一种可穿戴式磁场梯度探测仪及探测方法
CN112050799B (zh) * 2020-08-19 2022-11-18 哈尔滨工程大学 一种基于磁梯度张量缩并量之比的测距定位方法
CN112668146B (zh) * 2020-12-03 2021-08-24 重庆科技学院 一种基于欧拉反褶积法实用性改进的场源位置估算方法

Family Cites Families (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8253057B1 (en) * 2004-09-03 2012-08-28 Jack Hunt System and method for plasma generation
FR2942883B1 (fr) * 2009-03-06 2011-05-13 Commissariat Energie Atomique Capteur de gradient d'une composante d'un champ magnetique a aimant permanent
CN102401888B (zh) * 2011-08-24 2013-10-23 西安电子科技大学 一种电磁矢量传感器阵列耦合误差的自校正方法
UA102163C2 (ru) * 2012-02-02 2013-06-10 Владимир Николаевич Сосницкий Устройство для компенсации электромагнитных помех при измерениях биомагнитных сигналов
CN102988038B (zh) * 2012-12-14 2014-10-15 中国科学院上海微系统与信息技术研究所 一种用于无屏蔽心磁图仪的一阶梯度补偿模块及方法
CN103630139B (zh) * 2013-12-17 2015-12-02 哈尔滨工程大学 一种基于地磁梯度张量测量的水下载体全姿态确定方法
CN104091060B (zh) * 2014-06-30 2017-01-18 天津大学 一种分段式Halbach阵列永磁电机磁场计算方法
CN104330830B (zh) * 2014-10-30 2017-02-01 中国石油天然气集团公司 一种磁性体顶面埋深预测方法及装置
CN106291725B (zh) * 2015-05-13 2018-11-30 核工业北京地质研究院 一种快速反演地下地质体空间位置的方法
CN105091880B (zh) * 2015-07-17 2017-11-21 哈尔滨工程大学 一种基于标量传感器阵列的追踪定位水下远距离磁性目标的方法
CN107817457B (zh) * 2017-10-13 2020-03-17 北京工业大学 一种地磁梯度张量测量阵列的设计方法
US11953457B2 (en) * 2018-09-06 2024-04-09 Yissum Research Development Company Of The Hebrew University Of Jerusalem Ltd. Method for processing spin magnetometry data

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Magnetic dipole localization based on magnetic gradient tensor data at a single point;Gang Yin 等;《SPIE.DIGITAL LIBRARY》;20140703;全文 *
欧拉反褶积在重磁位场中应用与发展;王明 等;《物探与化探》;20121005;第36卷(第5期);全文 *

Also Published As

Publication number Publication date
CN109307847A (zh) 2019-02-05

Similar Documents

Publication Publication Date Title
CN109307847B (zh) 一种磁体的二阶欧拉反演法
CN107544042B (zh) 一种磁力计阵列校正方法
US7656160B2 (en) Determining properties of earth formations using the electromagnetic coupling tensor
US8466683B2 (en) Determining properties of earth formations using the electromagnetic coupling tensor
CN110146839A (zh) 一种移动平台磁梯度张量系统校正方法
CN111522835B (zh) 一种基于数据库特征匹配的多磁性目标位置探测方法
CN108388720A (zh) 一种基于无迹卡尔曼滤波的仿生偏振传感器多源误差标定方法
CN103353612B (zh) 一种地下目标物体的测量定位设备及测量定位方法
CN105043390B (zh) 基于泛克里金法的重力场插值方法
CN111220932B (zh) 无人机磁干扰标定方法及分布式磁异常探测系统
CN109633540B (zh) 一种磁源的实时定位系统及实时定位方法
CN109931956B (zh) 捷联式三分量磁测系统中三轴磁力仪与惯导安装误差校正方法
CN113552637A (zh) 一种航空-地面-井中磁异常数据协同三维反演方法
CN106918350B (zh) 一种应用于地磁导航中的地磁场模型误差补偿方法
Jiang et al. Scalar calibration of aeromagnetic data using BPANN and LS algorithms based on fixed-wing UAV platform
RU2572109C1 (ru) Способ калибровки электронного магнитного компаса
CN115524762A (zh) 基于三维亥姆赫兹线圈的地磁矢量测量系统补偿方法
CN110596619A (zh) 一种全张量磁梯度测量组件及其优化方法
CN115097370A (zh) 一种大平面测磁系统中自平衡矢量磁力仪转向差校准方法
RU2623192C1 (ru) Способ калибровки электронного магнитного компаса
CN104678448B (zh) 基于激光束和反光镜的磁干扰分量补偿方法
Yan et al. A compensation method in magnetic distortion through regularized inverse problems
CN104777440A (zh) 一种不需准确预知地磁倾角的岩矿石标本磁参数测量方法
CN114325511B (zh) 磁通门磁强计传感器的同点性设计计算方法及系统
CN112050804B (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
GR01 Patent grant