CN108957571A - 一种航空重力数据插值、扩边和下延一体化方法 - Google Patents

Abstract

本发明公开一种航空重力数据插值、扩边和下延一体化的方法及系统。该方法包括：获取航空重力数据，将航空重力数据有缺失数据和边界位置补零，获得补零重力数据；采用L‑曲线法确定最大截止波数，对补零重力数据进行二维傅里叶变换处理，计算补零重力数据的频谱；对频谱进行低通滤波处理，获得频谱的谱分量；将变换插值和扩边结果作二维傅里叶反变换获得航向重力数据的向下延拓结果。考虑了实测重力数据存在数据丢失和数据长度不符合快速傅里叶变换对数据长度的要求这两个问题；利用了重力数据插值、扩边和向下延拓同属不适定反问题的共性特点，建立了这三个问题的一体化解法方案，实现了它们的同步解决。

Description

一种航空重力数据插值、扩边和下延一体化方法

技术领域

本发明涉及航空重力测量领域，特别是涉及一种航空重力数据插值、扩边和下延一体化方法。

背景技术

航空重力测量是快速高效获取重力异常的典型方法，但获取仅仅是航向高度的重力异常。对大地测量、重力导航等许多应用领域来说，经常需要通过航空重力的向下延拓来获取地面或低于航线高度的重力异常。

重力场向下延拓是经典的不适定问题，对于该类问题的求解方法分为两类：一是直接法，二是迭代法。重力下延最经典的逆Poisson积分法就有直接求解和迭代求解两种形式。代表性迭代法包括积分迭代法、Landwber迭代法、迭代Tikhonov法等。现有技术中的航空重力迭代下延方法，存在的问题包括：仅仅解决了向下延拓的问题，基本都假设输入的重力数据是纯粹的数学向量或矩阵，即数据无空白；对波数域的下延，也大都不考虑快速傅里叶变换对数据长度的要求。

显然，实测重力资料往往存在空缺，且数据长度是随意的。因此，对航空实测重力数据进行向下延拓处理前，必须对数据空白部分进行插值，且波数域向下延拓方法还有必要对数据进行扩边，这将直接关系到重力向下延拓的精度。重力数据的插值和扩边对应于信号的重构和外推。信号的重构和外推是信号处理的逆问题，也是不适定的。所以，现有技术中都是将以往航空重力插值、扩边和下延这三个问题分开解决的，精度低，亟需一种能够实现航空重力数据插值、扩边和下延的同步高精度一体化解决。

发明内容

本发明的目的是提供一种能够实现航空重力数据插值、扩边和下延的同步高精度一体化解决的航空重力数据插值、扩边和下延一体化方法。

为实现上述目的，本发明提供了如下方案：

一种航空重力数据插值、扩边和下延一体化的方法，所述方法包括：

获取航空重力数据g(x,y)，将所述航空重力数据g(x,y)有缺失数据和边界位置补零，获得补零重力数据g′(x,y)；

采用L-曲线法计算截止波数C_K；

对所述补零重力数据g′(x,y)进行二维傅里叶变换处理，获得所述补零重力数据g′(x,y)的频谱G(u,v)；

对所述频谱G(u,v)进行低通滤波处理，获得所述频谱G(u,v)的谱分量G′(u,v)；

对所述谱分量G′(u,v)作二维傅里叶反变换，获得反变换数据；

将所述航空重力数据g(x,y)中有缺失数据的位置填充反变换数据中对应位置的反变换数据，所述航空重力数据g(x,y)中无缺失数据的位置保持原重力数据；

判断迭代次数是否等于预定迭代次数K-1，如果是，获得插值和扩边结果；

将所述插值和扩边结果作二维傅里叶变换，获得所述插值和扩边结果的频谱；

将所述插值和扩边结果的频谱应用算子获得变换插值和扩边结果；其中，(u,v)表示波数域的波数；

将所述变换插值和扩边结果作二维傅里叶反变换获得所述航向重力数据的向下延拓结果。

可选的，所述采用L-曲线法计算截止波数C_K具体包括：

由向下延拓正则解的范数和残差范数共同组成的极小化泛函

来确定截止波数C_K，其中，h(x,y)＝d/[2π(x²+y²+h²)^3/2]是转换的积分核。

可选的，所述对所述频谱G(u,v)进行低通滤波处理，获得所述频谱G(u,v)的谱分量G′(u,v)具体包括：

对所述频谱G(u,v)作用阈值函数T_k(u,v)

其中，c_k表示第k次迭代的截止波数，且满足c₁≤c₂≤c₃...≤c_k和1＜c_k≤min(M,N)，其中，M和N是插值扩边后所述重力数据的尺寸；D(u,v)表示波数域波数(u,v)与波数矩形中心的距离，D(u,v)＝[(u-M/2)²+(v-N/2)²]^1/2。

根据本发明提供的具体实施例，本发明公开了以下技术效果：本发明公开了一种航空重力数据插值、扩边和下延一体化的方法，考虑了实测重力数据存在数据丢失和数据长度不符合快速傅里叶变换对数据长度的要求这两个问题；利用了重力数据插值、扩边和向下延拓同属不适定反问题的共性特点，建立了这三个问题的一体化解法方案，实现了它们的同步解决。

本发明公开的一种航空重力数据插值、扩边和下延一体化的方法只需预先采用L-曲线法求得正则参数和预定义迭代次数，整个迭代过程的参数较少，算法操作简便、复杂度低；提高了插值、扩边和向下延拓结果的精度。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为本发明提供的一种航空重力数据插值、扩边和下延一体化的方法的流程图；

图2为实施例2中的仿真实验中的0km高度理论重力异常数据；

图3为实施例2中的仿真实验中的1km高度理论重力异常数据；

图4为实施例2中的仿真实验中的1km高度理论加噪重力异常数据；

图5为实施例2中的仿真实验中的挖去边界和内部部分数据的加噪数据；

图6为实施例2中的仿真实验中的求正则参数截止波数的曲线图；

图7为实施例2中的仿真实验中的加噪重力异常数据扩边插值结果图；

图8为实施例2中的仿真实验中的扩边插值后下延的结果图；

图9为实施例2中的仿真实验中的三种误差随迭代次数k的变化曲线。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明的目的是提供一种能够实现航空重力数据插值、扩边和下延的同步高精度一体化解决的航空重力数据插值、扩边和下延一体化的方法及系统。

为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。

实施例1

如图1所示，一种航空重力数据插值、扩边和下延一体化的方法，所述方法包括：

步骤100：获取航空重力数据g(x,y)，将所述航空重力数据g(x,y)有缺失数据和边界位置补零，获得补零重力数据g′(x,y)；

步骤200：采用L-曲线法计算截止波数C_K；

步骤300：对所述补零重力数据g′(x,y)进行二维傅里叶变换处理，获得所述补零重力数据g′(x,y)的频谱G(u,v)；

步骤400：对所述频谱G(u,v)进行低通滤波处理，获得所述频谱G(u,v)的谱分量G′(u,v)；

步骤500：对所述谱分量G′(u,v)作二维傅里叶反变换，获得反变换数据；

步骤600：将所述航空重力数据g(x,y)中有缺失数据的位置填充反变换数据中对应位置的反变换数据，所述航空重力数据g(x,y)中无缺失数据的位置保持原重力数据；

步骤700：判断迭代次数是否等于预定迭代次数K-1，如果是，执行步骤800：获得插值和扩边结果；

步骤900：将所述插值和扩边结果作二维傅里叶变换，获得所述插值和扩边结果的频谱；

步骤1000：将所述插值和扩边结果的频谱应用算子获得所属频谱的谱分量；其中，(u,v)表示波数域的波数；

步骤1100：将所述变换插值和扩边结果作二维傅里叶反变换获得所述航向重力数据的向下延拓结果；否则，返回执行步骤300继续对所述补零重力数据g′(x,y)迭代处理，每次迭代的截止波数从1开始呈线性增大，用于恢复出缺失数据的更多细节信息。

所述步骤200：采用L-曲线法计算截止波数C_K具体包括：

由向下延拓正则解的范数和残差范数共同组成的极小化泛函

来确定截止波数C_K，其中，h(x,y)＝d/[2π(x²+y²+h²)^3/2]是转换的积分核。

所述步骤400：对所述频谱G(u,v)进行低通滤波处理，获得所述频谱G(u,v)的谱分量G′(u,v)具体包括：

对所述频谱G(u,v)作用阈值函数T_k(u,v)

其中，c_k表示第k次迭代的截止波数，且满足c₁≤c₂≤c₃...≤c_k和1＜c_k≤min(M,N)，其中，M和N是插值扩边后所述重力数据的尺寸；D(u,v)表示波数域中点(u,v)与波数矩形中心的距离，D(u,v)＝[(u-M/2)²+(v-N/2)²]^1/2。

实施例2

对长方体重力模型进行仿真实验：

理论模型采用长方体重力场模型。长方体的长宽高分别为8km、4km、4km，中心埋深3km，剩余密度为1.0t/m³。正演计算地面(0km，200×200大小)和距离地面1km高度处(256×256大小)的重力数据分别如图2和3所示，虚框为长方体所在位置。为模拟实际情况并检验方法的抗噪声能力，给图3的重力数据增加零均值、方差为0.01的高斯白噪声，信噪比为39.63dB，其结果如图4所示。同时，为了检验方法插值和扩边能力，将图4的边界四个边界各28，共25536个数据和内部部分数据40×50共2000个数据，“挖去”构成空白区见图5所示。采用本发明提出的插值扩边下延一体化方法，将图5所示的加噪且有缺失重力数据插值扩边后，下延1km(10倍网格距)到地面，并以图5所示的1km高度处的真实重力数据检验方法的扩边和插值精度、以图2所示的地面重力数据检验方法的下延精度。

根据图1的算法步骤，先将图5所示的数据补零，然后求解截止波数c_K，如图6所示，最小值即为c_K＝12。选定迭代次数K＝100次，按照算法步骤，图5插值和扩边的结果见图7，图7向下延拓1km的结果见图8。由插值扩边结果图7和原始重力数据图4、下延结果图8和真实重力数据图2的对比可知，插值、扩边和下延的结果同理论结果一致。

为定量分析算法精度，插值、扩边和延拓精度统一由均方误差Rmse

来定量分析换算误差，其中，D_c(i)和D_t(i)分别表示计算值和理论值，Q表示计算误差的总数据量。插值、扩边和向下延拓均分误差随迭代次数k的变化如图9所示。由图9可知，算法具有收敛性，最终的插值、扩边和下延的误差分别为0.0938mGal、0.2297mGal和0.4483mGal，提高了结果的精度。

本发明的有益效果：

1)本发明考虑了实测航空重力数据往往存在数据缺失和数据长度不符合快速傅里叶变换对数据长度(2ⁿ)的要求这两个问题；

2)本发明利用了重力数据插值、扩边和向下延拓同属不适定反问题的共性特点，建立了这三个问题的一体化解法方案，实现了它们的同步解决；

3)本发明只需预先采用L-曲线法求得正则参数和预定义迭代次数，整个迭代过程的参数较少，算法操作简便、复杂度低；

4)本发明的仿真结果表明，方法具有很高的插值、扩边和向下延拓精度；

5)本发明基于凸集投影原理，将航空重力数据的插值、扩边与向下延拓这三个不适定问题统一考虑，提出航空重力数据同时插值、扩边和下延的一体化方法。

本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的系统而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。

本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。

Claims (3)

1.一种航空重力数据插值、扩边和下延一体化的方法，其特征在于，所述方法包括：

获取航空重力数据g(x,y)，将所述航空重力数据g(x,y)有缺失数据和边界位置补零，获得补零重力数据g′(x,y)；

采用L-曲线法计算截止波数C_K；

对所述补零重力数据g′(x,y)进行二维傅里叶变换处理，获得所述补零重力数据g′(x,y)的频谱G(u,v)；

对所述频谱G(u,v)进行低通滤波处理，获得所述频谱G(u,v)的谱分量G′(u,v)；

对所述谱分量G′(u,v)作二维傅里叶反变换，获得反变换数据；

将所述航空重力数据g(x,y)中有缺失数据的位置填充反变换数据中对应位置的反变换数据，所述航空重力数据g(x,y)中无缺失数据的位置保持原重力数据；

判断迭代次数是否等于预定迭代次数K-1，如果是，获得插值和扩边结果；

将所述插值和扩边结果作二维傅里叶变换，获得所述插值和扩边结果的频谱；

将所述插值和扩边结果的频谱应用算子获得向下延拓结果的频谱；其中，(u,v)表示波数域的波数；

将所述变换插值和扩边结果作二维傅里叶反变换获得所述航向重力数据的向下延拓结果。

2.根据权利要求1所述的一种航空重力数据插值、扩边和下延一体化的方法，其特征在于，所述采用L-曲线法计算截止波数C_K具体包括：

由向下延拓正则解的范数和残差范数共同组成的极小化泛函

来确定截止波数C_K，其中，h(x,y)＝d/[2π(x²+y²+h²)^3/2]是转换的积分核。

3.根据权利要求1所述的一种航空重力数据插值、扩边和下延一体化的方法，其特征在于，所述对所述频谱G(u,v)进行低通滤波处理，获得所述频谱G(u,v)的谱分量G′(u,v)具体包括：

对所述频谱G(u,v)作用阈值函数T_k(u,v)

其中，ck表示第k次迭代的截止波数，且满足c₁≤c₂≤c₃...≤c_k和1＜c_k≤min(M,N)，其中，M和N是插值扩边后所述重力数据的尺寸；D(u,v)表示波数域波数(u,v)与波数矩形中心的距离，D(u,v)＝[(u-M/2)²+(v-N/2)²]^1/2。