CN113686329A - 一种基于地磁数据的垂直高度位场测量方法 - Google Patents

一种基于地磁数据的垂直高度位场测量方法 Download PDF

Info

Publication number
CN113686329A
CN113686329A CN202110996145.1A CN202110996145A CN113686329A CN 113686329 A CN113686329 A CN 113686329A CN 202110996145 A CN202110996145 A CN 202110996145A CN 113686329 A CN113686329 A CN 113686329A
Authority
CN
China
Prior art keywords
plane
frequency
continuation
spectrum
frequency spectrum
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
Application number
CN202110996145.1A
Other languages
English (en)
Other versions
CN113686329B (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.)
National University of Defense Technology
Original Assignee
National University of Defense Technology
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 National University of Defense Technology filed Critical National University of Defense Technology
Priority to CN202110996145.1A priority Critical patent/CN113686329B/zh
Publication of CN113686329A publication Critical patent/CN113686329A/zh
Application granted granted Critical
Publication of CN113686329B publication Critical patent/CN113686329B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C21/00—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
    • G01C21/04—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by terrestrial means
    • G01C21/08—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by terrestrial means involving use of the magnetic field of the earth
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/30—Assessment of water resources

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Environmental & Geological Engineering (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • Geology (AREA)
  • Automation & Control Theory (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明公开一种基于地磁数据的垂直高度位场测量方法，步骤包括：S1.获取观测平面的地磁位场数据并进行频谱分析，得到观测平面的频谱；S2.根据观测平面的频谱确定截止频率，并根据截止频率构建自适应滤波算子；S3.以观测平面的频谱作为初值进行向下延拓的迭代更新，每次迭代时将当前延拓平面的频谱向上延拓到观测平面得到延拓值，根据观测平面的频谱与延拓值计算得到剩余谱，使用自适应滤波算子处理剩余谱并更新延拓平面的频谱；S4.将得到的延拓平面的频谱进行傅里叶反变换，得到空间域位场。本发明能够得到不同高度平面上的地磁场数据，且具有实现方法简单、效率及精度高且稳定性好数据等优点。

Description

一种基于地磁数据的垂直高度位场测量方法
技术领域
本发明涉及地磁匹配导航的地磁数据库构建中位场测量技术领域，尤其涉及一种基于地磁数据的垂直高度位场测量方法。
背景技术
地磁匹配导航是一种全天候、全地域的无源定位导航手段，隐蔽性好，且能够实现自主导航。地磁匹配导航需要完全依赖地球本身的地磁场信息，以实时测量的地磁数据为基础，通过地磁数据延拓的方法构建地磁数据库。在导航时，将载体运动轨迹上实时测量获得的地磁数据经过处理生成实时图；然后，将实时地磁图与数据库中的地磁特征信息进行匹配，从而获得载体当前的位置估计信息以达到定位导航的目的。因而完备的高精度地磁数据库是实现地磁匹配导航的前提，在地磁匹配导航技术体系中具有十分重要的地位。
目前，利用地磁场延拓技术构建地磁数据库主要针对具有垂直高度关系的位场数据，即用已知高度的地磁场信息来测算未知高度地磁场信息。一般地磁场数据的获取方式主要采用无人机搭载磁测系统进行航磁测量，但是航测磁数据仅在某一测量高度面上，未测量高度将出现数据真空带，无法满足运行在不同高度面上载体的地磁匹配导航需求。因此，研究由已知高度的地磁场测量数据来向下测算未知高度地磁场数据的方法具有十分重要的意义。
针对地磁数据库垂直高度位场的构建，目前的测量方法中不能根据磁场数据进行自适应的调整，无法充分利用数据中的高频信息，而与此同时，具有垂直高度关系的地磁数据构建的稳定性、收敛速度与磁数据的频率成分具有很大的相关性，目前的有构建方法通常都不加以区分，导致既无法根据不同频率磁数据进行调整，又无法同时保证算法稳定性与收敛速度。
为实现地磁数据库垂直高度位场的构建，需要进行地磁场向下延拓，而地磁场向下延拓属于典型的不适定问题。针对地磁场向下延拓目前几乎都是通过近似逼近的方式来求解该问题，主要有以下几种方式：
(1)积分迭代法：在众多的迭代算法中，积分迭代法表现出了较为突出的优势，尤其是向下延拓积分迭代法经模型，在无噪声的情况下向下延拓深度可达20倍点距，但是积分迭代法能够收敛到直接向下延拓的理论解，同时在观测数据含噪声时，积分迭代法会使得噪声累积，从而影响延拓精度，因此积分迭代法中的噪声累积极大地影响了延拓效果。
(2)泰勒级数迭代法：最早期的空间域位场解析向下延拓方法即是通过泰勒级数实现，但是泰勒级数的相关方法对高频的压制效果较差，容易受到噪声干扰，抗干扰能力差。
(3)最小二乘迭代法：利用梯度下降法更新迭代步长，可以加快收敛速度，但是该方法的收敛速度极慢，特别对于非“平滑解”，延拓精度极低，无法满足地磁数据库的精度标准。
综上，现有技术中地磁数据库存在在不同高度平面上的数据不完备性且精度低的问题，且基于地磁数据库测算垂直高度位场时，由于均是采用半收敛算子，会随迭代次数增加而容易发散，所以延拓结果对迭代次数选择较为敏感，且各类迭代方法的收敛速度、稳定性以及延拓精度受到数据中高频数据影响较大。
发明内容
本发明要解决的技术问题就在于：针对现有技术存在的技术问题，本发明提供一种实现方法简单、效率及精度高且稳定性好的基于地磁数据的垂直高度位场测量方法。
为解决上述技术问题，本发明提出的技术方案为：
一种基于地磁数据的垂直高度位场测量方法，步骤包括：
S1.频谱分析：获取观测平面的地磁位场数据并进行频谱分析，得到观测平面的频谱；
S2.自适应滤波算子构建：根据所述观测平面的频谱确定截止频率，并根据所述截止频率按照低通滤波器模型构建自适应滤波算子；
S3.向下延拓：以所述观测平面的频谱作为初值进行向下延拓的迭代更新，每次迭代时将当前延拓平面的频谱向上延拓到观测平面得到延拓值，根据所述观测平面的频谱与所述延拓值计算得到剩余谱，使用所述自适应滤波算子处理所述剩余谱并更新所述延拓平面的频谱，直至得到最终的延拓平面的频谱输出；
S4.测量结果输出：将步骤S3最终得到的延拓平面的频谱进行傅里叶反变换，得到空间域位场的测量结果输出。
进一步的，所述步骤S1的步骤包括：
S11.建立三维直角坐标系，其中地平面为Z＝0平面，垂直地平面向下为Z轴正方向，观测平面为T(z0)，观测平面的地磁位场数据为
Figure BDA0003233874970000021
S12.将观测平面的地磁位场数据
Figure BDA0003233874970000031
向下延拓到靠近场源的平面T(z1)，得到延拓平面位场数据
Figure BDA0003233874970000032
向下延拓高度为h＝z0-z1，z0为观测平面的高度，z1为靠近场源的平面的高度；
S13.分别将所述观测平面的地磁位场数据
Figure BDA0003233874970000033
延拓平面位场数据
Figure BDA0003233874970000034
以及空间域向上延拓核函数K(x,y)进行二维傅里叶变换，得到观测平面的频谱
Figure BDA0003233874970000035
延拓平面的频谱
Figure BDA0003233874970000036
以及向上延拓算子Hup(kx,ky)。
进一步的，所述地磁位场数据
Figure BDA0003233874970000037
延拓平面位场数据
Figure BDA0003233874970000038
以及空间域向上延拓核函数K(x,y)满足：
Figure BDA0003233874970000039
Figure BDA00032338749700000310
进一步的，所述步骤S2中确定截止频率的步骤包括：
S21.将观测平面的频谱零频移位到中心计算得到位场功率谱P(kx,ky)；
S22.以所述位场功率谱P(kx,ky)的中心作圆形成一系列环带，取所述环带内所有点的功率值作平均，得到径向平均功率谱；
S23.根据所述径向平均功率谱随频率的衰减程度变化确定出截止频率。
进一步的，所述步骤S23中，以径向频率wr为横坐标、对所述径向平均功率谱取对数为纵坐标，获取得到对数径向平均功率谱曲线，将所述对数径向平均功率谱曲线衰减程度出现差异的拐点的横坐标作为所述截止频率。
进一步的，所述自适应滤波算子通过使用低通滤波器模型，并使用所述截止频率对所述低通滤波器模型进行自适应调整构建得到，具体为：
Figure BDA00032338749700000311
其中，wc为截止频率，w为频率，n可取2或3。
进一步的，所述步骤S3的步骤包括：
S31.把观测平面的位场频谱F(kx,ky)作为向下延拓的延拓平面频谱的初值G0(kx,ky)，即G0(kx,ky)＝F(kx,ky)；
S32.将当前延拓平面频谱Gi(kx,ky)向上延拓到观测平面得到延拓值，i＝0,1,…n，并将所述观测平面的位场频谱F(kx,ky)与所述延拓值作差得到剩余谱e，即：
e＝F(kx,ky)-Hup(kx,ky)·Gi(kx,ky)
其中，Hup(kx,ky)为将空间域向上延拓核函数K(x,y)进行二维傅里叶变换得到的向上延拓算子；
S33.使用所述自适应滤波算子处理所述剩余谱e得到修正项；
S34.使用所述修正项更新当前延拓平面的频谱，即：
Gi+1(kx,ky)＝Gi(kx,ky)+ΔG
S35.判断是否达到停止迭代条件，如果是则得到最终的延拓平面频谱Gn(kx,ky)输出，否则返回步骤S32。
进一步的，所述步骤S33中按照下式得到所述修正项：
ΔG＝φ·e
其中，ΔG为所述修正项，φ为所述自适应滤波算子。
进一步的，所述步骤S35中当迭代次数达到预设最大迭代次数或迭代结果误差小于预设阈值时，判定达到停止迭代条件。
与现有技术相比，本发明的优点在于：
1、本发明通过依据地磁位场数据的频谱分析构建自适应滤波算子，在向下延拓的迭代过程中利用自适应滤波算子处理剩余谱，可以达到低通滤波效果，自适应的调整滤波算子的低通频带，并尽可能利用数据中高频有用信息而抑制噪声部分，提高向下延拓迭代的收敛速度，从而能够快速、精确的得到地磁数据空间域位场，解决传统地磁数据库中在不同高度平面上的数据不完备性、精度低的问题。
2、本发明在向下延拓的迭代过程中，利用自适应滤波算子的低通特性可以延缓发散过程，提高迭代的稳定性，解决传统采用半收敛算子的迭代方法中会随着迭代次数增加而逐渐发散的问题，保证在较长迭代区间都可以收敛到最佳延拓效果，基于数据自适应调整可以实现收敛速度快、稳定性好、延拓精度高的向下延拓迭代。
3、本发明进一步通过分析平均径向功率谱随频率的衰减程度，定性的找到一个区分位场信号与噪声的截止频率，基于该截止频率能够尽可能多的利用到高频位场信号。
附图说明
图1是本实施例向下延拓的原理示意图。
图2是本实施例基于地磁数据的垂直高度位场测量方法的实现流程示意图。
图3是在具体应用实施例中得到的不同截止频率下的自适应滤波算子与向上延拓算子滤波特性结果示意图。
图4是具体应用实施例中得到的不同埋深模型的观测平面位场数据的功率谱曲线示意图。
图5是具体应用实施例中采用本发明方法得到的磁源埋深200m模型的延拓效果示意图。
图6是具体应用实施例中采用本发明方法得到的磁源埋深1000m模型的延拓效果示意图。
图7是具体应用实施例中采用本发明方法得到的实测数据的延拓效果示意图。
具体实施方式
以下结合说明书附图和具体优选的实施例对本发明作进一步描述，但并不因此而限制本发明的保护范围。
如图1所示，本实施例基于地磁数据的垂直高度位场测量方法的步骤包括：
S1.频谱分析：获取观测平面的地磁位场数据并进行频谱分析，得到观测平面的频谱；
S2.自适应滤波算子构建：根据观测平面的频谱确定截止频率，并根据截止频率按照低通滤波器模型构建自适应滤波算子；
S3.向下延拓：以观测平面的频谱作为初值进行向下延拓的迭代更新，每次迭代时将当前延拓平面的频谱向上延拓到观测平面得到延拓值，根据观测平面的频谱与延拓值计算得到剩余谱，使用自适应滤波算子处理剩余谱并更新延拓值，直至得到最终的延拓平面的频谱输出；
S4.测量结果输出：将步骤S3最终得到的延拓平面的频谱进行傅里叶反变换，得到空间域位场的测量结果输出。
本实施例通过依据地磁位场数据的频谱分析构建自适应滤波算子，在向下延拓的迭代过程中利用自适应滤波算子处理剩余谱，可以达到低通滤波效果，自适应的调整滤波算子的低通频带，并尽可能利用数据中高频有用信息而抑制噪声部分，提高向下延拓迭代的收敛速度，能够快速、精确的得到地磁数据空间域位场，解决传统地磁数据库中在不同高度平面上的数据不完备性、精度低的问题。
上述在向下延拓的迭代过程中，还可以利用自适应滤波算子的低通特性延缓发散过程，提高迭代的稳定性，解决传统采用半收敛算子的迭代方法中会随着迭代次数增加而逐渐发散的问题，保证在较长迭代区间都可以收敛到最佳延拓效果。
本实施例中步骤S1的步骤包括：
S11.建立三维直角坐标系，其中地平面为Z＝0平面，垂直地平面向下为Z轴正方向，观测平面为T(z0)，观测平面的地磁位场数据为
Figure BDA0003233874970000061
S12.将观测平面的地磁位场数据
Figure BDA0003233874970000062
向下延拓到靠近场源的平面T(z1)，得到延拓平面位场数据
Figure BDA0003233874970000063
向下延拓高度为h＝z0-z1，z0为观测平面的高度，z1为靠近场源的平面的高度；
S13.分别将观测平面的地磁位场数据
Figure BDA0003233874970000064
延拓平面位场数据
Figure BDA0003233874970000065
以及空间域向上延拓核函数K(x,y)进行二维傅里叶变换，得到观测平面的频谱
Figure BDA0003233874970000066
延拓平面的频谱
Figure BDA0003233874970000067
以及向上延拓算子Hup(kx,ky)。
上述地磁位场数据
Figure BDA0003233874970000068
延拓平面位场数据
Figure BDA0003233874970000069
以及空间域向上延拓核函数K(x,y)之间关系可表示为：
Figure BDA00032338749700000610
Figure BDA00032338749700000611
其中，x,y为已知观测平面位场的位置坐标，ξ,γ为延拓平面位场的位置坐标。
将观测平面的地磁位场数据
Figure BDA00032338749700000612
延拓平面位场数据
Figure BDA00032338749700000613
进行二维傅里叶变换得到相应的频谱为：
Figure BDA00032338749700000614
Figure BDA00032338749700000615
应用积分变化表，得到频率域向上延拓算子Hup(kx,ky)即为：
Figure BDA00032338749700000616
通过傅里叶变换，把空间域的卷积运算转换到频率域的乘积运算，可得：
Figure BDA00032338749700000617
进一步可以简化为：
F＝Hup·G (8)
通过上述变换，可以将向下延拓问题转换为相当于已知Hup、F，求解G的解线性方程组的过程，从而使得可以与最小二乘问题联系在一起，即找到G*使得：
Figure BDA0003233874970000071
后续经过迭代求解后即可得到上述G。
迭代求解前需先确定截止频率，本实施例步骤S2中确定截止频率的步骤包括：
S21.将观测平面的频谱零频移位到中心计算得到位场功率谱P(kx,ky)；
S22.以位场功率谱P(kx,ky)的中心作圆形成一系列环带，取环带内所有点的功率值作平均，得到径向平均功率谱；
S23.根据径向平均功率谱随频率的衰减程度变化确定出截止频率wc。
由于功率谱的逆傅里叶变换是自相关谱，而观测数据中的高频噪声与位场信号是不相关的，但是由于混叠着小部分位场信号的高频成分，在功率谱的衰减应该是较为平缓，而位场信号中大部分的高低频信息会导致功率谱较大的衰减程度，即通过分析平均径向功率谱随频率的衰减程度，可以定性的找到一个区分位场信号与噪声的截止频率，同时能够尽可能多的利用到高频位场信号。
在具体应用实施例中，由式(6)计算频率域的向上延拓算子Hup(kx,ky)，以及由式(4)得到
Figure BDA0003233874970000072
后，将频谱零频移位到中心计算其功率谱P(kx,ky)，即：
Figure BDA0003233874970000073
其中M，N为网格化数据的行数、列数。
本实施例进一步以位场功率谱P(kx,ky)的中心作圆，圆的半径wr分别为基频
Figure BDA0003233874970000074
的整数倍，从而形成一系列环带；取环带内所有点的功率值作平均，即得到径向平均功率谱
Figure BDA0003233874970000075
然后分析观测平面位场的频率分布，以径向频率wr为横坐标，对径向平均功率谱取对数为纵坐标，可得到对数径向平均功率谱曲线
Figure BDA0003233874970000076
根据对数径向平均功率谱曲线
Figure BDA0003233874970000077
的衰减程度差异确定截止频率wc，具体将曲线衰减程度出现差异的拐点的横坐标作为截止频率wc。
本实施例中自适应滤波算子通过使用低通滤波器模型，并使用截止频率对低通滤波器模型进行自适应调整构建得到。自适应滤波算子φ(w)实质上是低通滤波器，且可以根据位场数据的频率成分自适应的调整通带宽度，实现对高频噪声的压制，保证算法稳定性，且调整截止频率可以改变低通频带范围，使之与位场的频率成分相对应，建立基于数据的自适应滤波算子，充分利用高频有用信号。
在具体应用实施例中，先设计低通滤波算子如下，保证迭代算法稳定性。
Figure BDA0003233874970000081
利用截止频率wc对低通滤波器φ进行自适应调整，建立基于数据-截止频率-低通滤波算子的对应关系。
Figure BDA0003233874970000082
综合截止频率wc可调性与标准低通滤波算子φ(w)，得到最终的自适应滤波算子，
Figure BDA0003233874970000083
其中，wc为截止频率，w为频率，n可取2或3。
采用上述自适应滤波算子，能够自适应的调整滤波算子的低通频带，尽可能利用数据中高频成分而压制噪声部分，提高向下延拓迭代收敛速度与实用性。
本实施例中步骤S3的步骤包括：
S31.把观测平面的位场频谱F(kx,ky)作为向下延拓的延拓平面频谱的初值G0(kx,ky)，即：
G0(kx,ky)＝F(kx,ky) (15)
S32.将当前延拓平面频谱Gi(kx,ky)向上延拓到观测平面得到延拓值，i＝0,1,…n，并将观测平面的位场频谱F(kx,ky)与延拓值作差得到剩余谱e，即：
e＝F(kx,ky)-Hup(kx,ky)·Gi(kx,ky) (16)
其中，Hup(kx,ky)为将空间域向上延拓核函数K(x,y)进行二维傅里叶变换得到的向上延拓算子；
S33.使用自适应滤波算子处理剩余谱e得到修正项，具体可按照下式计算得到：
ΔG＝φ·e (17)
S34.使用修正项更新当前延拓平面的频谱，即：
Gi+1(kx,ky)＝Gi(kx,ky)+ΔG (18)
S35.判断是否达到停止迭代条件，如果是则得到最终的延拓平面频谱Gn(kx,ky)输出，否则返回步骤S32。
剩余谱中包含位场数据中的中低频成分、高频成分与噪声部分，最佳效果是保留有用信息而滤除噪声，但是高频成分与噪声有混叠区域。本实施例利用自适应滤波算子处理剩余谱，能够达到低通滤波效果，在充分利用中低频信息的前提下，能够尽可能保留高频有用信号，抑制噪声对迭代过程的影响，提高稳定性，同时调整低通频带充分利用高频有用信息，提高收敛速度。
上述步骤S35中具体当迭代次数达到预设最大迭代次数或迭代结果误差小于预设阈值(|e|<ε(ε为很小的数)时，判定达到停止迭代条件。
本实施例中将迭代结束后得到的频谱后，具体按照下式通过傅里叶反变换得到待延拓平面位场数据：
Figure BDA0003233874970000091
磁源埋深体现在空间域为梯度变化大小，转化到频率域对应频率成分高低，磁源埋深浅的区域对应于空间中梯度较大，相当于频率域中的高频成分。本实施例通过利用径向平均功率谱分析磁数据的频率成分分布确定截止频率，基于截止频率构建自适应滤波算子，以通过对数据频率成分的分析，自适应的调整滤波算子的低通频带，使得能够尽可能利用数据中高频成分而压制噪声部分，提高迭代过程的收敛速度与稳定性，且能够保证在较长迭代区间都可以收敛到最佳延拓效果，从而基于数据自适应调整可以实现收敛速度快、稳定性好、延拓精度高的向下延拓迭代。
为验证本发明的有效性，在具体应用实施例中在向下延拓迭代过程中分别利用本发明自适应滤波算子与传统向上延拓算子并比较结果，不同截止频率下的观测平面位场数据的功率谱曲线，如图3所示，本发明自适应滤波算子与传统向上延拓算子虽然都是低通滤波器，但是本发明自适应滤波算子不仅具有低频通带，过滤带窄，且可以根据位场数据的频率成分自适应的调整通带宽度(截止频率越大，通带越长，可以使用的频率成分也就越多)，有效地利用频谱的中高频分量，提升迭代过程的性能，尤其适用于浅源位场的延拓。同时，低通特性还可以很好压制高频噪声的影响，提高算法稳定性。
为验证本发明上述方法的有效性，进一步在具体应用实例中使用本发明上述方法进行仿真实验，仿真条件具体为：
(1)采用的球体的磁异常表达式为：
Figure BDA0003233874970000101
其中μ0表示真空磁导率；D表示磁化强度的偏角；I表示地磁场倾角，也表示磁化强度倾角；m表示磁矩，它与磁化率J之间的关系：
m＝J×πR3 (18)
其中，R表示球体半径。
(2)组合球体模型的球体磁源数量为5个，参数如表1。
表1 组合球体模型参数
Figure BDA0003233874970000102
接下来，通过设置两组不同埋深的仿真试验，以验证采用本发明方法可以根据数据频率成分的差异自适应调整，得到最佳延拓效果，并且保证收敛速度与稳定性。
本实施例定义三个精度指标来定量分析本发明的性能：
1)均方根误差Δgrmse
Figure BDA0003233874970000103
2)误差矩阵Δgme
Δgme＝|ge-gt| (20)
3)相对误差变化率Δgre
Figure BDA0003233874970000104
进一步本实施例对两组不同埋深模型的观测平面进行谱分析，不同埋深模型的观测平面位场数据的功率谱曲线如图4所示。
根据对数功率谱衰减程度对应于数据中不同频率成分，按照本发明方法可以确定截止频率为：
1)埋深为200m的浅源模型截止频率为0.0041，
2)埋深为1000m的深源模型的截止频率为0.0023
从上述可看出，浅源相对于深源，频率成分中高频占比增大，对数功率谱衰减程度的拐点越靠右。
进一步本实施例将球体磁源埋深统一设置为200m，得到的结果如图5所示，其中图5(a)是Z＝200m平面的位场，图5(b)是Z＝700m平面的位场。图5(c)是本发明延拓结果，图5(d)是延拓结果图5(c)与理论平面图5(a)作差得到的误差矩阵，图5(d)是迭代过程的均方根误差曲线，分析各结果图可得到：
1)由图5(c)的延拓结果与理论平面相比，可以直观的看出本发明取得极好的延拓效果，完全拟合到理论平面，仅仅存在一些边界的吉布斯效应，可以通过扩边方法进一步提升。
2)从图5(d)的误差矩阵结果中可以直观地看到延拓误差的分布，即本发明方法在大部分梯度变化小的区域以及磁源主体区域的延拓误差几乎为0，误差仅仅分布在磁源主体的中心部分(频率成分极高)以及由于磁源靠近边界导致的褶皱部分。
3)从图5(e)的均方根误差曲线中看得出，本发明不仅具有较快收敛速度，而且保证在较长迭代区间内保持最佳延拓效果，稳定性极好，可以抑制迭代算法由于本身的半收敛性质导致在迭代后期迅速发散。
进一步本实施例将球体磁源埋深统一设置为1000m，结果如图6所示，其中图6(a)是Z＝200m平面的位场，图6(b)是Z＝700m平面的位场。图6(c)是本发明延拓结果，图5(d)是延拓结果图6(c)与理论平面图6(a)作差得到的误差矩阵，图6(d)是迭代过程的均方根误差曲线，分析各结果图可得到：
1)埋深1000m的位场数据梯度变化较小，大部分为中低频信息，可得本发明(图6(c))几乎完全与理论平面相同，延拓效果相比埋深浅模型更好。
2)通过图6(d)的误差矩阵可以看出，误差分布情况与埋深浅模型相同，但是误差范围控制在0.5nT以内，具有极小的延拓误差，在低频数据延拓中展示出极高的延拓精度。
3)通过如图6(e)的均方根误差曲线可以看出，深源模型中表现出较好的延拓精度与稳定性，本发明在此基础上还进一步提高了收敛速度。
本实施例进一步分析不同埋深下的延拓效果，表2给出了延拓精度指标：最小均方根误差、最大误差矩阵值以及最小相对误差率。
表2 本发明延拓精度指标
Figure BDA0003233874970000121
本实施例通过两组不同埋深的模型仿真试验，其中埋深200m相当于频率成分中高频分量占比大，埋深1000m相当于频率分量主要以低频成分为主，从表中看出埋深代表的频率成分对延拓效果有很大影响，高频分量越多的位场数据，向下延拓的误差越大，即根据数据频率成分来自适应调整迭代算法是十分必要的，相对于传统迭代方法，本发明在对不同频率成分的磁数据都能够表现出收敛速度快、延拓精度高以及较强稳定性等优势。
本实施例进一步实测数据采用某地区航测的磁异常数据，对实测数据的延拓效果如图7所示，其中图7(a)是航测平面的位场数据，图7(b)是将图7(a)向上延拓500m后的位场数据。图7(c)是采用本发明方法将图7(b)向下延拓500m后的延拓结果，可以看到本发明在实际使用中也可以得到较好的延拓效果，延拓结果与理论平面比较误差较小，对位场中梯度变化复杂区域也有较好的延拓结果。
上述只是本发明的较佳实施例，并非对本发明作任何形式上的限制。虽然本发明已以较佳实施例揭露如上，然而并非用以限定本发明。因此，凡是未脱离本发明技术方案的内容，依据本发明技术实质对以上实施例所做的任何简单修改、等同变化及修饰，均应落在本发明技术方案保护的范围内。

Claims (9)

1.一种基于地磁数据的垂直高度位场测量方法，其特征在于，步骤包括：
S1.频谱分析：获取观测平面的地磁位场数据并进行频谱分析，得到观测平面的频谱；
S2.自适应滤波算子构建：根据所述观测平面的频谱确定截止频率，并根据所述截止频率按照低通滤波器模型构建自适应滤波算子；
S3.向下延拓：以所述观测平面的频谱作为初值进行向下延拓的迭代更新，每次迭代时将当前延拓平面的频谱向上延拓到观测平面得到延拓值，根据所述观测平面的频谱与所述延拓值计算得到剩余谱，使用所述自适应滤波算子处理所述剩余谱并更新所述延拓平面的频谱，直至得到最终的延拓平面的频谱输出；
S4.测量结果输出：将步骤S3最终得到的延拓平面的频谱进行傅里叶反变换，得到空间域位场的测量结果输出。
2.根据权利要求1所述的基于地磁数据的垂直高度位场测量方法，其特征在于，所述步骤S1的步骤包括：
S11.建立三维直角坐标系，其中地平面为Z＝0平面，垂直地平面向下为Z轴正方向，观测平面为T(z0)，观测平面的地磁位场数据为
Figure FDA0003233874960000013
S12.将观测平面的地磁位场数据
Figure FDA0003233874960000015
向下延拓到靠近场源的平面T(z1)，得到延拓平面位场数据
Figure FDA0003233874960000014
向下延拓高度为h＝z0-z1，z0为观测平面的高度，z1为靠近场源的平面的高度；
S13.分别将所述观测平面的地磁位场数据
Figure FDA0003233874960000016
延拓平面位场数据
Figure FDA0003233874960000017
以及空间域向上延拓核函数K(x,y)进行二维傅里叶变换，得到观测平面的频谱
Figure FDA0003233874960000018
延拓平面的频谱
Figure FDA0003233874960000019
以及向上延拓算子Hup(kx,ky)。
3.根据权利要求2所述的基于地磁数据的垂直高度位场测量方法，其特征在于，所述地磁位场数据
Figure FDA00032338749600000110
延拓平面位场数据
Figure FDA00032338749600000111
以及空间域向上延拓核函数K(x,y)满足：
Figure FDA0003233874960000011
Figure FDA0003233874960000012
4.根据权利要求1所述的基于地磁数据的垂直高度位场测量方法，其特征在于，所述步骤S2中确定截止频率的步骤包括：
S21.将观测平面的频谱零频移位到中心计算得到位场功率谱P(kx,ky)；
S22.以所述位场功率谱P(kx,ky)的中心作圆形成一系列环带，取所述环带内所有点的功率值作平均，得到径向平均功率谱；
S23.根据所述径向平均功率谱随频率的衰减程度变化确定出截止频率。
5.根据权利要求4所述的基于地磁数据的垂直高度位场测量方法，其特征在于，所述步骤S23中，以径向频率wr为横坐标、对所述径向平均功率谱取对数为纵坐标，获取得到对数径向平均功率谱曲线，将所述对数径向平均功率谱曲线衰减程度出现差异的拐点的横坐标作为所述截止频率。
6.根据权利要求1～5中任意一项所述的基于地磁数据的垂直高度位场测量方法，其特征在于，所述自适应滤波算子通过使用低通滤波器模型，并使用所述截止频率对所述低通滤波器模型进行自适应调整构建得到，具体为：
Figure FDA0003233874960000021
其中，wc为截止频率，w为频率，n可取2或3。
7.根据权利要求1～5中任意一项所述的基于地磁数据的垂直高度位场测量方法，其特征在于，所述步骤S3的步骤包括：
S31.把观测平面的位场频谱F(kx,ky)作为向下延拓的延拓平面频谱的初值G0(kx,ky)，即G0(kx,ky)＝F(kx,ky)；
S32.将当前延拓平面频谱Gi(kx,ky)向上延拓到观测平面得到延拓值，i＝0,1,…n，并将所述观测平面的位场频谱F(kx,ky)与所述延拓值作差得到剩余谱e，即：
e＝F(kx,ky)-Hup(kx,ky)·Gi(kx,ky)
其中，Hup(kx,ky)为将空间域向上延拓核函数K(x,y)进行二维傅里叶变换得到的向上延拓算子；
S33.使用所述自适应滤波算子处理所述剩余谱e得到修正项；
S34.使用所述修正项更新当前延拓平面的频谱，即：
Gi+1(kx,ky)＝Gi(kx,ky)+ΔG
S35.判断是否达到停止迭代条件，如果是则得到最终的延拓平面频谱Gn(kx,ky)输出，否则返回步骤S32。
8.根据权利要求7所述的基于地磁数据的垂直高度位场测量方法，其特征在于，所述步骤S33中按照下式得到所述修正项：
ΔG＝φ·e
其中，ΔG为所述修正项，φ为所述自适应滤波算子。
9.根据权利要求7所述的基于地磁数据的垂直高度位场测量方法，其特征在于，所述步骤S35中当迭代次数达到预设最大迭代次数或迭代结果误差小于预设阈值时，判定达到停止迭代条件。
CN202110996145.1A 2021-08-27 2021-08-27 一种基于地磁数据的垂直高度位场测量方法 Active CN113686329B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110996145.1A CN113686329B (zh) 2021-08-27 2021-08-27 一种基于地磁数据的垂直高度位场测量方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110996145.1A CN113686329B (zh) 2021-08-27 2021-08-27 一种基于地磁数据的垂直高度位场测量方法

Publications (2)

Publication Number Publication Date
CN113686329A true CN113686329A (zh) 2021-11-23
CN113686329B CN113686329B (zh) 2023-07-25

Family

ID=78583431

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110996145.1A Active CN113686329B (zh) 2021-08-27 2021-08-27 一种基于地磁数据的垂直高度位场测量方法

Country Status (1)

Country Link
CN (1) CN113686329B (zh)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114167511A (zh) * 2021-11-26 2022-03-11 兰州大学 一种基于连分式展开向下延拓的位场数据快速反演方法
CN115292660A (zh) * 2022-08-04 2022-11-04 中国自然资源航空物探遥感中心 一种自适应迭代的位场分离方法、系统、设备及计算机可读存储介质
CN116774301A (zh) * 2023-06-26 2023-09-19 中国人民解放军火箭军工程大学 重磁位场正则化向下延拓方法、系统、电子设备及介质
CN116774301B (zh) * 2023-06-26 2024-05-14 中国人民解放军火箭军工程大学 重磁位场正则化向下延拓方法、系统、电子设备及介质

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1184547A (zh) * 1995-04-07 1998-06-10 美国3M公司 采用自适应滤波和数字检测方法的电子物品监视系统
CN1877366A (zh) * 2006-07-12 2006-12-13 杨辉 重磁延拓回返垂直导数目标优化处理技术
CN107656272A (zh) * 2017-09-25 2018-02-02 厦门大学 一种在高阶时域算法下的电磁波三维逆时偏移成像方法
CN108957571A (zh) * 2018-06-01 2018-12-07 中国人民解放军火箭军工程大学 一种航空重力数据插值、扩边和下延一体化方法
CN109100816A (zh) * 2018-07-09 2018-12-28 中国人民解放军火箭军工程大学 一种重磁数据处理方法及系统
CN110365051A (zh) * 2019-07-31 2019-10-22 南京工程学院 一种自适应指令滤波反演的虚拟同步电机控制方法
CN111337992A (zh) * 2020-03-23 2020-06-26 兰州大学 一种基于位场数据向下延拓的场源深度获得方法
CN111504278A (zh) * 2020-04-20 2020-08-07 哈尔滨工程大学 基于自适应频域积分的海浪检测方法

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1184547A (zh) * 1995-04-07 1998-06-10 美国3M公司 采用自适应滤波和数字检测方法的电子物品监视系统
CN1877366A (zh) * 2006-07-12 2006-12-13 杨辉 重磁延拓回返垂直导数目标优化处理技术
CN107656272A (zh) * 2017-09-25 2018-02-02 厦门大学 一种在高阶时域算法下的电磁波三维逆时偏移成像方法
CN108957571A (zh) * 2018-06-01 2018-12-07 中国人民解放军火箭军工程大学 一种航空重力数据插值、扩边和下延一体化方法
CN109100816A (zh) * 2018-07-09 2018-12-28 中国人民解放军火箭军工程大学 一种重磁数据处理方法及系统
CN110365051A (zh) * 2019-07-31 2019-10-22 南京工程学院 一种自适应指令滤波反演的虚拟同步电机控制方法
CN111337992A (zh) * 2020-03-23 2020-06-26 兰州大学 一种基于位场数据向下延拓的场源深度获得方法
CN111504278A (zh) * 2020-04-20 2020-08-07 哈尔滨工程大学 基于自适应频域积分的海浪检测方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
J. WANG,ET AL: "A constrained scheme for high precision downward continuation of potential field data", 《PURE APPL. GEOPHYS.》 *

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114167511A (zh) * 2021-11-26 2022-03-11 兰州大学 一种基于连分式展开向下延拓的位场数据快速反演方法
CN114167511B (zh) * 2021-11-26 2023-08-11 兰州大学 一种基于连分式展开向下延拓的位场数据快速反演方法
CN115292660A (zh) * 2022-08-04 2022-11-04 中国自然资源航空物探遥感中心 一种自适应迭代的位场分离方法、系统、设备及计算机可读存储介质
CN115292660B (zh) * 2022-08-04 2023-03-24 中国自然资源航空物探遥感中心 一种自适应迭代的位场分离方法、系统、设备及计算机可读存储介质
CN116774301A (zh) * 2023-06-26 2023-09-19 中国人民解放军火箭军工程大学 重磁位场正则化向下延拓方法、系统、电子设备及介质
CN116774301B (zh) * 2023-06-26 2024-05-14 中国人民解放军火箭军工程大学 重磁位场正则化向下延拓方法、系统、电子设备及介质

Also Published As

Publication number Publication date
CN113686329B (zh) 2023-07-25

Similar Documents

Publication Publication Date Title
CN113686329B (zh) 一种基于地磁数据的垂直高度位场测量方法
CN104614718B (zh) 基于粒子群算法的激光雷达波形数据分解的方法
CN111007571B (zh) 一种基于三维构造张量的航磁数据地质体边界识别方法
CN108680966B (zh) 海洋可控源电磁勘探噪声降噪效果评估方法
CN105929380B (zh) 一种卫星激光高度计全波形激光雷达数据去噪方法
CN107870355B (zh) 一种复杂地形条件下的克希霍夫型波束偏移方法
CN110989002B (zh) 地空时域电磁系统浅部低阻异常体数据解释方法
CN109100816B (zh) 一种重磁数据处理方法及系统
CN109459745B (zh) 一种利用辐射噪声估计运动声源速度的方法
CN107632964A (zh) 一种平面地磁异常场向下延拓递归余弦变换法
CN105445798B (zh) 一种基于梯度处理的全波形反演方法和系统
CN110361788B (zh) 一种空-地联合三维重力数据特征分析及密度反演方法
CN112882115A (zh) 基于gwo优化小波阈值的大地电磁信号去噪方法及系统
CN111427096A (zh) 全张量重力梯度仪数据质量评价和滤波处理方法
CN111665556B (zh) 地层声波传播速度模型构建方法
CN113124882B (zh) 一种背景磁场未知情形下的多磁偶极子磁源反演定位方法
Luo et al. Accurate tree roots positioning and sizing over undulated ground surfaces by common offset GPR measurements
CN113866836A (zh) 一种基于归一化磁异常导数标准差的多目标边界识别方法
CN115979245A (zh) 一种自校准估计的磁感应网络定位方法
CN114488217B (zh) 基于深度学习的高轨卫星cei信号频率估计方法
CN113504568B (zh) 一种基于小生境差分进化算法的中值滤波方法
CN111142134B (zh) 一种坐标时间序列处理方法及装置
CN113852433A (zh) 一种基于计算机视觉的无线信道阴影衰落模型预测方法
Wang et al. A New Potential-Field Downward Continuation Iteration Method Based on Adaptive Filtering
CN113985494A (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