重磁位场正则化向下延拓方法、系统、电子设备及介质

1.本发明涉及重力场和地磁场数据测量技术领域，特别是涉及一种重磁位场正则化向下延拓方法、系统、电子设备及介质。


背景技术：

2.向下延拓能够缩短观测面与场源间的距离，是重磁位场数据解释和重磁辅助导航基准图构建的一种有效方法。众所周知，重磁位场向下延拓会放大测量数据中的高频分量而导致下延结果不稳定。为了获得更加稳定、精确的下延结果，不同学者提出了不同类型的方法，归纳起来大致可以分为三类。第一类为空间域方法，主要包括边界单元法和有限元法、样条函数法、等效源法等。一般来说，空间域方法精度较高，但计算相对较复杂。第二类方法为波数域方法，包括维纳滤波法、补偿圆滑滤波法、正则化方法、各类迭代方法等。波数域方法计算简单，但针对的一般是规则网格数据，且需考虑波数域中运算诸如栅栏效应和边界效应等问题。第三类方法为设法结合前两类方法优点的空间波数混合域方法。
3.事实上，向下延拓的不稳定性在波数域表现在下延因子对数据中高频成份的放大作用。因此，各类波数域稳定下延方法大部分都是通过附加低通滤波器或改造下延因子来压制下延因子对高频的放大作用，而低通滤波器或改造下延因子都存在截止波数或所谓正则参数的确定问题。正则参数确定方法包括l-曲线法、c-范数法、广义交叉校验法、拟最优法等。这些方法的共同点在于：首先确定一个正则参数离散变化的范围(为了提高计算精度，离散值的个数越多越好)，然后计算该范围内所有正则参数的下延结果，并最终基于一定的准则来确定最优正则参数，如l-曲线法的拐点或最大曲率、c-范数曲线、广义交叉校验和拟最优准则曲线的最小值等。由于需要遍历正则参数的变化范围，这些正则参数确定方法也被称为“序列方法”，普遍存在计算量偏大的问题。通过拟合位场径向平均功率谱径向谱曲线确定截止波数的方法虽然物理意义明确，但存在需要人工观察拐点或拟合径向谱的不足。


技术实现要素：

4.本发明的目的是提供一种重磁位场正则化向下延拓方法、系统、电子设备及介质，可自动、快速获取重磁位场向下延拓参数，且具有较高的向下延拓精度。
5.为实现上述目的，本发明提供了如下方案：
6.一种重磁位场正则化向下延拓方法，包括：
7.获取重磁位场数据；
8.对所述重磁位场数据进行快速傅里叶变换得到重磁位场数据功率谱；
9.以所述重磁位场数据功率谱的中心为圆心，以基频的整数倍为半径作圆，生成多个环带；
10.以径向波数为横坐标，以各环带对应的对数值作为纵坐标，得到对数分形修正径向平均功率谱；所述环带对应的对数值为对所述环带的平均功率谱值与分形修正谱的乘积
取对数得到的；
11.以所述对数分形修正径向平均功率谱的最小值对应的径向波数作为截止波数；
12.对所述重磁位场数据进行傅里叶变换得到重磁位场数据傅里叶谱；
13.根据所述截止波数、所述重磁位场数据的径向波数、所述重磁位场数据傅里叶谱以及向下延拓算子得到重磁位场的向下延拓结果。
14.可选的，所述径向波数为横坐标，以各环带对应的对数值作为纵坐标，得到对数分形修正径向平均功率谱，具体包括：
15.对于任意一个环带，计算所述环带内所有点的功率谱值的平均值，得到所述环带的平均功率谱值；
16.对所述环带的平均功率谱值与分形修正谱的乘积取对数得到所述环带对应的对数值；
17.以径向波数为横坐标，以各所述环带对应的对数值为纵坐标，得到对数分形修正径向平均功率谱。
18.可选的，所述根据所述截止波数、所述重磁位场数据的径向波数、所述重磁位场数据傅里叶谱以及向下延拓算子得到重磁位场的向下延拓结果，具体包括：
19.将所述截止波数和所述重磁位场数据的径向波数输入正则低通滤波器得到滤波值；
20.计算所述滤波值、所述重磁位场数据傅里叶谱以及向下延拓算子的乘积得到向下延拓乘积；
21.对所述向下延拓乘积进行傅里叶反变换得到重磁位场的向下延拓结果。
22.可选的，所述正则低通滤波器为：
23.其中，表示正则低通滤波器，h表示重磁位场数据向下延拓的深度，ωr表示重磁位场数据的径向波数，ωc表示截止波数。
24.一种重磁位场正则化向下延拓系统，包括：
25.获取模块，用于获取重磁位场数据；
26.重磁位场数据功率谱确定模块，用于对所述重磁位场数据进行快速傅里叶变换得到重磁位场数据功率谱；
27.环带生成模块，用于以所述重磁位场数据功率谱的中心为圆心，以基频的整数倍为半径作圆，生成多个环带；
28.对数分形修正径向平均功率谱确定模块，用于以径向波数为横坐标，以各环带对应的对数值作为纵坐标，得到对数分形修正径向平均功率谱；所述环带对应的对数值为对所述环带的平均功率谱值与分形修正谱的乘积取对数得到的；
29.截止波数确定模块，用于以所述对数分形修正径向平均功率谱的最小值对应的径向波数作为截止波数；
30.傅里叶变换模块，用于对所述重磁位场数据进行傅里叶变换得到重磁位场数据傅
里叶谱；
31.向下延拓模块，用于根据所述截止波数、所述重磁位场数据的径向波数、所述重磁位场数据傅里叶谱以及向下延拓算子得到重磁位场的向下延拓结果。
32.可选的，所述对数分形修正径向平均功率谱确定模块，具体包括：
33.平均功率谱值计算单元，用于对于任意一个环带，计算所述环带内所有点的功率谱值的平均值，得到所述环带的平均功率谱值；
34.对数值计算单元，用于对所述环带的平均功率谱值与分形修正谱的乘积取对数得到所述环带对应的对数值；
35.对数分形修正径向平均功率谱确定单元，用于以径向波数为横坐标，以各所述环带对应的对数值为纵坐标，得到对数分形修正径向平均功率谱。
36.可选的，所述向下延拓模块，具体包括：
37.滤波值计算单元，用于将所述截止波数和所述重磁位场数据的径向波数输入正则低通滤波器得到滤波值；
38.向下延拓乘积计算单元，用于计算所述滤波值、所述重磁位场数据傅里叶谱以及向下延拓算子的乘积得到向下延拓乘积；
39.向下延拓单元，用于对所述向下延拓乘积进行傅里叶反变换得到重磁位场的向下延拓结果。
40.可选的，所述正则低通滤波器为：
41.其中，表示正则低通滤波器，h表示重磁位场数据向下延拓的深度，ωr表示重磁位场数据的径向波数，ωc表示截止波数。
42.一种电子设备，包括：
43.存储器和处理器，所述存储器用于存储计算机程序，所述处理器运行所述计算机程序以使所述电子设备执行根据上述所述的重磁位场正则化向下延拓方法。
44.一种计算机可读存储介质，其存储有计算机程序，所述计算机程序被处理器执行时实现上述所述的重磁位场正则化向下延拓方法。
45.根据本发明提供的具体实施例，本发明公开了以下技术效果：
46.本发明通过以径向波数为横坐标，以对所述环带的平均功率谱值与分形修正谱的乘积取对数为纵坐标，得到对数分形修正径向平均功率谱，以对数分形修正径向平均功率谱的最小值对应的径向波数作为截止波数，可以自动确定截止波数，解决了自动、快速的问题方法，具有较高的向下延拓精度。
附图说明
47.为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图
获得其他的附图。
48.图1为本发明实施例提供的重磁位场正则化向下延拓方法流程图；
49.图2为信号谱、噪声谱及其组合随径向波数的变化示意图；
50.图3为加入高斯白噪声的加噪重力模型数据等值线图；
51.图4为图3中三种高斯白噪声水平重力数据的径向平均功率谱对比图；
52.图5为图3中三种高斯白噪声水平重力数据的分形修正径向平均功率谱对比图；
53.图6为2％高斯白噪声噪声水平重力数据在不同尺度指数β下的分形修正径向平均功率谱对比图；
54.图7为三种方法向下延拓结果与理论重力在主剖面上的对比图；
55.图8为三种方法向下延拓结果与理论重力的残差在主剖面上的对比图。
具体实施方式
56.下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
57.为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。
58.设z轴向下为正，u0(ω
x
,ωy)为观测平面z＝0上的原始数据u0(x,y)的傅里叶变换结果，uh(ω
x
,ωy)为z＝h(h＞0，场源位于z＝h平面以下)平面上的位场数据uh(x,y)的傅里叶变换结果，则它们在波数域的关系为：
[0059][0060]
其中，ω
x
和ωy分别为x和y方向的波数；为径向波数；为向下延拓算子。由于对数据高频信号及噪声的放大特性，导致下延结果uh(ω
x
,ωy)不稳定，且距离h越大，下延越不稳定。为稳定换算，通常增加一个低通滤波器l(ωr)来压制高频噪声的干扰，即
[0061][0062]
其中，u
′h(ω
x
,ωy)为稳定的下延结果。基于此，本发明实施例提供了一种重磁位场正则化向下延拓方法，如图1所示，所述方法具体包括：
[0063]
步骤101：获取重磁位场数据u0(x,y)。
[0064]
步骤102：对所述重磁位场数据进行快速傅里叶变换得到重磁位场数据功率谱。
[0065]
步骤103：以所述重磁位场数据功率谱的中心为圆心，以基频的整数倍为半径作圆，生成多个环带。
[0066]
步骤104：以径向波数为横坐标，以各环带对应的对数值作为纵坐标，得到对数分形修正径向平均功率谱；所述环带对应的对数值为对所述环带的平均功率谱值p(ωr)与分
形修正谱的乘积取对数，即得到的。
[0067]
步骤105：以所述对数分形修正径向平均功率谱的最小值对应的径向波数作为(信号和噪声大致分界的)截止波数ωc。
[0068]
步骤106：对所述重磁位场数据u0(x,y)进行傅里叶变换得到重磁位场数据傅里叶谱u0(ω
x
,ωy)。
[0069]
步骤107：根据所述截止波数、所述重磁位场数据的径向波数、所述重磁位场数据傅里叶谱以及向下延拓算子得到重磁位场的向下延拓结果。
[0070]
在实际应用中，所述以径向波数为横坐标，以各环带对应的对数值作为纵坐标，得到对数分形修正径向平均功率谱，具体包括：
[0071]
对于任意一个环带，计算所述环带内所有点的功率谱值的平均值，得到所述环带的平均功率谱值。
[0072]
对所述环带的平均功率谱值与分形修正谱的乘积取对数得到所述环带对应的对数值。
[0073]
以径向波数为横坐标，以各所述环带对应的对数值为纵坐标，得到对数分形修正径向平均功率谱。
[0074]
在实际应用中，所述根据所述截止波数、所述重磁位场数据的径向波数、所述重磁位场数据傅里叶谱以及向下延拓算子得到重磁位场的向下延拓结果，具体包括：
[0075]
将所述截止波数和所述重磁位场数据的径向波数输入正则低通滤波器得到滤波值让小于截止波数ωc的信号谱全通，仅对大于截止波数ωc的信号谱进行低通滤波。
[0076]
计算所述滤波值、所述重磁位场数据傅里叶谱以及向下延拓算子的乘积即得到向下延拓乘积。
[0077]
对所述向下延拓乘积进行傅里叶反变换即得到重磁位场的向下延拓结果。
[0078]
在实际应用中，所述正则低通滤波器为：
[0079][0080]
其中，表示正则低通滤波器，h表示重磁位场数据向下延拓的深度，ωr表示重磁位场数据的径向波数，ωc表示截止波数。
[0081]
本发明提供了一个实施例对上述方法的有效性进行验证，具体步骤如下：
[0082]
spector和grant(1970)证明空间随机不相关分布长方形棱柱体磁源的深度h与磁源径向平均功率谱p(ωr)存在简单的指数关系：
[0083][0084]
其中，a为由源强度决定的常数。对式(4)两边取对数可得：
[0085]
ln(p(ωr))＝lna-2ωrh(5)
[0086]
通过在一定波数范围内拟合对数径向平均功率谱曲线，并由拟合后直线的斜率来获得磁源的深度。因此，式(5)及其一些变形广泛应用于磁源顶部、底部和质心的估计。后续的研究表明，随机不相关磁源模型存在过渡简化的问题，磁化率分布具有典型空间自相似性和尺度不变规律，可用分形幂律模型来描述，式(4)的模型需改进为：
[0087][0088]
其中，2≤β≤4为径向谱的尺度指数。对于磁性基底尺度大、结构性强的磁性块体，通常具有β＝2.9的分形结构。
[0089]
前人的研究大都只考虑了径向谱在磁源深度估计中的应用，很少考虑信号和噪声分界点的估计。如果考虑实测数据均含有噪声(一般假设为白噪声)的情况，那么式(6)的模型显然应完善为：
[0090][0091]
其中，pn为白噪声对应的功率常数。对式(7)两边乘以并分别取对数可得：
[0092][0093]
式(8)的右边对数括号中，是随径向波数ωr单调递减的函数，而是随径向波数ωr单调递增的函数，由于a和pn均为常数，它们的和将大致呈现“v”字形，取自然对数后，单调递增的自然对数函数也不会改变曲线的大致形状；如此，式(8)的左边以分形幂律模型对数据功率谱p(ωr)进行修正后的对数功率谱也将大致呈现“v”字形。在“v”字形曲线左端对应的低、中波数区域，是由有效信号占主导的区域，而在“v”字形曲线右端对应的高波数区域，是由噪声信号占主导，由此，以“v”字形曲线极小值对应的波数作为有效信号和噪声大致分界的截止波数ωc具有一定合理性。
[0094]
令常数a＝0.0001、源深度h＝500、径向波数ωr∈[0,0.01]、噪声水平pn＝30、β＝2.9，则径向平均功率谱分形修正噪声谱两项之和及其取自然对数结果随径向波数的变化情况如图2所示。由图2所见，信号径向平均功率谱随径向波数ωr单调下降，分形修正噪声谱随径向波数ωr单调上升，它们两项之和为“凹”函数，其曲线呈现出“v”字形，同时，其对数结果同样呈现“v”字形。因此，可很自然地取对数
最小值所对应的波数作为信号和噪声大致分界的截止波数ωc。
[0095]
为验证本发明所提滤波器参数自动确定方法的有效性，采用经典的双球体重力模型来进行验证。两个球体的参数如表1所示(z坐标向下为正)。线数m和每线点数n同为512、点距δx和线距δy都为50m。
[0096]
表1双球体重力模型参数表
[0097][0098]
试验中，以h＝0m平面处的重力数据为观测数据，并分别加入零均值、均方差为该高度理论重力异常绝对均值0.2％(0.0001mgal，信噪比62.03db)、2％(0.001mgal，信噪比42.03db)和20％(0.01mgal，信噪比22.03db)的高斯白噪声以模拟实际情况，其结果如图3(a)、图3(b)和图3(c)所示，图3(a)为加入0.2％的高斯白噪声的加噪重力模型数据等值线图，图3(b)为加入2％的高斯白噪声的加噪重力模型数据等值线图，图3(c)为加入20％的高斯白噪声的加噪重力模型数据等值线图。分别计算图3(a)到图3(c)含噪重力异常数据的对数径向平均功率谱ln(p(ωr))，其结果如图4所示。由图4可知：(1)当噪声水平越小，对应的对数径向平均功率谱拐点越不明显，目测或线性拟合较难确定信号谱和噪声谱的大致分界点；(2)噪声水平越小，对数径向平均功率谱拐点对应的径向波数越大。采用分形幂律模型对图3(a)到图3(c)含噪重力异常数据径向平均功率谱p(ωr)进行修正后的对数径向平均功率谱如图5所示。由图5可知：(1)采用分形幂律模型对径向平均功率谱进行修正后得到的谱线具有明显的“v”字形特征，最小值位置很方便确定；(2)随着噪声水平的增大，修正径向平均功率谱曲线最小值对应的径向波数依次减小。在加入2％噪声情况下，分形幂律模型采用不同尺度指数(β＝2、2.9、4)的结果如图6所示。由图6可知：(1)尺度指数β越大，曲线的“v”字形特征越明显；(2)不同尺度指数β修正下曲线的最小值具有一致性，换句话说，在2≤β≤4范围内，修正径向功率谱的最小值是一致的，并不随尺度指数β的变化而变化。
[0099]
接下来，我们以h＝0m平面处的含噪重力数据为观测数据，然后采用ti khonov正则低通滤波器(式(3)、方法1)、改进tikhonov正则低通滤波器(式(4)、方法1)和改进导数迭代法(王泽et al.,2022)、方法3)将其向下延拓1000m(20倍点距)，并将得到的延拓值uc与1000m高度处的真实重力场值u
t
采用均方误差(root mean square error，rmse)
[0100][0101]
和相对误差(relative error，re)
[0102][0103]
来计算、对比延拓误差。
[0104]
通过图5，可以由添加0.2％、2％、20％不同噪声水平重力数据的修正径向平均功率谱曲线对应的最小值，确定对应的低通滤波器截止波数ωc依次为1.10
×
10-3
、9.38
×
10-4
、6.64
×
10-4
，令低通滤波器在这些截止波数ωc处降为0.5，则利用正则参数和截止波数的关系可求得对应的正则参数分别为6.57
×
10-7
、7.65
×
10-6
、2.38
×
10-4
。两种低通滤波器及改进导数迭代法相关参数及三种噪声水平下的延拓误差对比如表2、表3和表4所示。
[0105]
为了延拓效果的直观对比，将三种方法向下延拓后得到的延拓数据与h＝1000m平面上真实重力数据在两球心上方主剖面上的值放到一起，其对比如图7(a)、图7(b)和图7(c)所示，其中，图7(a)为采用三种方法对0.2％高斯白噪声重力观测数据进行处理得到的向下延拓结果与理论重力在主剖面上的对比图，图7(b)为采用三种方法对2％高斯白噪声重力观测数据进行处理得到的向下延拓结果与理论重力在主剖面上的对比图，图7(c)为采用三种方法对20％高斯白噪声重力观测数据进行处理得到的向下延拓结果与理论重力在主剖面上的对比图；三个延拓结果与真实值之差如图8(a)、图8(b)和图8(c)所示，其中，图8(a)为采用三种方法对0.2％高斯白噪声重力观测数据进行处理得到的向下延拓结果与理论重力的残差在主剖面上的对比图，图8(b)为采用三种方法对2％高斯白噪声重力观测数据进行处理得到的向下延拓结果与理论重力的残差在主剖面上的对比图，图8(c)为采用三种方法对20％高斯白噪声重力观测数据进行处理得到的向下延拓结果与理论重力的残差在主剖面上的对比图。由表2到表4及图7到图8可知：(1)正则化方法采用本发明提出的重磁位场正则化向下延拓方法获取的对应正则参数，其下延结果优于改进导数迭代法，验证了本发明所提方法的有效性；(2)在低噪声水平下，正则低通滤波器的下延误差小于改进正则低通滤波器(即式3)，但随着噪声水平的提高，改进正则低通滤波器的优势逐渐提升。
[0106]
表2理论重力模型数据加0.2％噪声向下延拓结果对比表
[0107][0108]
表3理论重力模型数据加2％噪声向下延拓结果对比表
[0109][0110]
表4理论重力模型数据加20％噪声向下延拓结果对比表
[0111][0112]
针对上述方法本发明还提供了一种重磁位场正则化向下延拓系统，包括：
[0113]
获取模块，用于获取重磁位场数据。
[0114]
重磁位场数据功率谱确定模块，用于对所述重磁位场数据进行快速傅里叶变换得到重磁位场数据功率谱。
[0115]
环带生成模块，用于以所述重磁位场数据功率谱的中心为圆心，以基频的整数倍为半径作圆，生成多个环带。
[0116]
对数分形修正径向平均功率谱确定模块，用于以径向波数为横坐标，以各环带对应的对数值作为纵坐标，得到对数分形修正径向平均功率谱；所述环带对应的对数值为对所述环带的平均功率谱值与分形修正谱的乘积取对数得到的。
[0117]
截止波数确定模块，用于以所述对数分形修正径向平均功率谱的最小值对应的径向波数作为截止波数。
[0118]
傅里叶变换模块，用于对所述重磁位场数据进行傅里叶变换得到重磁位场数据傅里叶谱。
[0119]
向下延拓模块，用于根据所述截止波数、所述重磁位场数据的径向波数、所述重磁位场数据傅里叶谱以及向下延拓算子得到重磁位场的向下延拓结果。
[0120]
作为一种可选的实施方式，所述对数分形修正径向平均功率谱确定模块，具体包括：
[0121]
平均功率谱值计算单元，用于对于任意一个环带，计算所述环带内所有点的功率谱值的平均值，得到所述环带的平均功率谱值。
[0122]
对数值计算单元，用于对所述环带的平均功率谱值与分形修正谱的乘积取对数得到所述环带对应的对数值。
[0123]
对数分形修正径向平均功率谱确定单元，用于以径向波数为横坐标，以各所述环
带对应的对数值为纵坐标，得到对数分形修正径向平均功率谱。
[0124]
作为一种可选的实施方式，所述向下延拓模块，具体包括：
[0125]
滤波值计算单元，用于将所述截止波数和所述重磁位场数据的径向波数输入正则低通滤波器得到滤波值。
[0126]
向下延拓乘积计算单元，用于计算所述滤波值、所述重磁位场数据傅里叶谱以及向下延拓算子的乘积得到向下延拓乘积。
[0127]
向下延拓单元，用于对所述向下延拓乘积进行傅里叶反变换得到重磁位场的向下延拓结果。
[0128]
作为一种可选的实施方式，所述正则低通滤波器为：
[0129]
其中，表示正则低通滤波器，h表示重磁位场数据向下延拓的深度，ωr表示重磁位场数据的径向波数，ωc表示截止波数。
[0130]
本发明实施例还提供了一种电子设备，包括：
[0131]
存储器和处理器，所述存储器用于存储计算机程序，所述处理器运行所述计算机程序以使所述电子设备执行上述实施例所述的重磁位场正则化向下延拓方法。
[0132]
本发明实施例还提供了一种计算机可读存储介质，其存储有计算机程序，所述计算机程序被处理器执行时实现上述实施例所述的重磁位场正则化向下延拓方法。
[0133]
本发明有以下技术效果：
[0134]
本发明从分形修正的位场径向平均功率谱理论模型出发，基于分形修正对数径向谱曲线的特征，提出了重磁位场正则化向下延拓方法，通过求分形修正径向谱曲线的最小值就可以自动确定向下延拓参数(截止波数)，解决了重磁位场向下延拓参数自动、快速获取的问题，且具有较高的向下延拓精度。
[0135]
本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的系统而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。
[0136]
本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。

技术特征：
1.一种重磁位场正则化向下延拓方法，其特征在于，包括：获取重磁位场数据；对所述重磁位场数据进行快速傅里叶变换得到重磁位场数据功率谱；以所述重磁位场数据功率谱的中心为圆心，以基频的整数倍为半径作圆，生成多个环带；以径向波数为横坐标，以各环带对应的对数值作为纵坐标，得到对数分形修正径向平均功率谱；所述环带对应的对数值为对所述环带的平均功率谱值与分形修正谱的乘积取对数得到的；以所述对数分形修正径向平均功率谱的最小值对应的径向波数作为截止波数；对所述重磁位场数据进行傅里叶变换得到重磁位场数据傅里叶谱；根据所述截止波数、所述重磁位场数据的径向波数、所述重磁位场数据傅里叶谱以及向下延拓算子得到重磁位场的向下延拓结果。2.根据权利要求1所述的重磁位场正则化向下延拓方法，其特征在于，所述以径向波数为横坐标，以各环带对应的对数值作为纵坐标，得到对数分形修正径向平均功率谱，具体包括：对于任意一个环带，计算所述环带内所有点的功率谱值的平均值，得到所述环带的平均功率谱值；对所述环带的平均功率谱值与分形修正谱的乘积取对数得到所述环带对应的对数值；以径向波数为横坐标，以各所述环带对应的对数值为纵坐标，得到对数分形修正径向平均功率谱。3.根据权利要求1所述的重磁位场正则化向下延拓方法，其特征在于，所述根据所述截止波数、所述重磁位场数据的径向波数、所述重磁位场数据傅里叶谱以及向下延拓算子得到重磁位场的向下延拓结果，具体包括：将所述截止波数和所述重磁位场数据的径向波数输入正则低通滤波器得到滤波值；计算所述滤波值、所述重磁位场数据傅里叶谱以及向下延拓算子的乘积得到向下延拓乘积；对所述向下延拓乘积进行傅里叶反变换得到重磁位场的向下延拓结果。4.根据权利要求1所述的重磁位场正则化向下延拓方法，其特征在于，所述正则低通滤波器为：正则低通滤波器，h表示重磁位场数据向下延拓的深度，ω
r
表示重磁位场数据的径向波数，ω
c
表示截止波数。5.一种重磁位场正则化向下延拓系统，其特征在于，包括：获取模块，用于获取重磁位场数据；重磁位场数据功率谱确定模块，用于对所述重磁位场数据进行快速傅里叶变换得到重
磁位场数据功率谱；环带生成模块，用于以所述重磁位场数据功率谱的中心为圆心，以基频的整数倍为半径作圆，生成多个环带；对数分形修正径向平均功率谱确定模块，用于以径向波数为横坐标，以各环带对应的对数值作为纵坐标，得到对数分形修正径向平均功率谱；所述环带对应的对数值为对所述环带的平均功率谱值与分形修正谱的乘积取对数得到的；截止波数确定模块，用于以所述对数分形修正径向平均功率谱的最小值对应的径向波数作为截止波数；傅里叶变换模块，用于对所述重磁位场数据进行傅里叶变换得到重磁位场数据傅里叶谱；向下延拓模块，用于根据所述截止波数、所述重磁位场数据的径向波数、所述重磁位场数据傅里叶谱以及向下延拓算子得到重磁位场的向下延拓结果。6.根据权利要求5所述的重磁位场正则化向下延拓系统，其特征在于，所述对数分形修正径向平均功率谱确定模块，具体包括：平均功率谱值计算单元，用于对于任意一个环带，计算所述环带内所有点的功率谱值的平均值，得到所述环带的平均功率谱值；对数值计算单元，用于对所述环带的平均功率谱值与分形修正谱的乘积取对数得到所述环带对应的对数值；对数分形修正径向平均功率谱确定单元，用于以径向波数为横坐标，以各所述环带对应的对数值为纵坐标，得到对数分形修正径向平均功率谱。7.根据权利要求5所述的重磁位场正则化向下延拓系统，其特征在于，所述向下延拓模块，具体包括：滤波值计算单元，用于将所述截止波数和所述重磁位场数据的径向波数输入正则低通滤波器得到滤波值；向下延拓乘积计算单元，用于计算所述滤波值、所述重磁位场数据傅里叶谱以及向下延拓算子的乘积得到向下延拓乘积；向下延拓单元，用于对所述向下延拓乘积进行傅里叶反变换得到重磁位场的向下延拓结果。8.根据权利要求5所述的重磁位场正则化向下延拓系统，其特征在于，所述正则低通滤波器为：正则低通滤波器，h表示重磁位场数据向下延拓的深度，ω
r
表示重磁位场数据的径向波数，ω
c
表示截止波数。9.一种电子设备，其特征在于，包括：存储器和处理器，所述存储器用于存储计算机程序，所述处理器运行所述计算机程序以使所述电子设备执行根据权利要求1至4中任一项所述的重磁位场正则化向下延拓方法。
10.一种计算机可读存储介质，其特征在于，其存储有计算机程序，所述计算机程序被处理器执行时实现如权利要求1至4中任一项所述的重磁位场正则化向下延拓方法。

技术总结
本发明公开一种重磁位场正则化向下延拓方法、系统、电子设备及介质，涉及重力场和地磁场数据测量技术领域。所述方法包括对重磁位场数据进行快速傅里叶变换得到重磁位场数据功率谱；以重磁位场数据功率谱的中心为圆心，以基频的整数倍为半径作圆，生成多个环带；以径向波数为横坐标，以各环带对应的对数值作为纵坐标，得到对数分形修正径向平均功率谱；以对数分形修正径向平均功率谱的最小值对应的径向波数作为截止波数；对重磁位场数据进行傅里叶变换得到重磁位场数据傅里叶谱；根据截止波数、径向波数、重磁位场数据傅里叶谱以及向下延拓算子得到向下延拓结果。本发明可自动、快速获取重磁位场向下延拓参数，且具有较高的向下延拓精度。下延拓精度。下延拓精度。


技术研发人员：曾小牛 李夕海 谭笑枫 牛超 刘继昊 张云 刘志刚
受保护的技术使用者：中国人民解放军火箭军工程大学
技术研发日：2023.06.26
技术公布日：2023/9/20