一种重磁数据处理方法及系统

摘要

本发明公开了一种重磁数据处理方法及系统。所述方法及系统将重磁数据的去噪、填充和扩边这三个问题统一考虑，基于凸集投影方法同时实现了重磁数据的去噪、填充和扩边处理。所述方法及系统首先对原始重磁异常数据补零，并通过计算归一化径向平均功率谱来确定迭代过程最终的截止波数，然后对重磁数据进行去噪、填充和扩边迭代直至达到预定的迭代次数，获得处理后的重磁数据。本发明提供的方法及系统原理简单、操作方便、数据处理效率高，同时具有较高的去噪、填充精度，且填充和扩边的衔接处光滑无畸变，非常适合重磁数据的同时去噪和扩充。

权利要求

1.一种重磁数据处理方法，其特征在于，所述方法包括：获取重磁异常数据；

将所述重磁异常数据中有数据缺失的位置补零，获得补零后重磁位场数据；

计算所述补零后重磁位场数据的归一化径向平均功率谱；所述计算所述补零后重磁位场数据的归一化径向平均功率谱，具体包括：采用快速傅里叶变换方法求解所述补零后重磁位场数据的重磁位场数据功率谱；

以所述重磁位场数据功率谱的中心为圆心，以基频的整数倍为半径作圆，生成多个环带；每个所述环带对应一个不同的频率；

求解各个所述环带内所有点的功率谱值的平均值，获得多个平均功率谱值；

以径向波数为横坐标，以各个所述环带内的所述平均功率谱值的对数为纵坐标，生成径向平均功率谱；

对所述径向平均功率谱作归一化处理，生成所述重磁位场数据的归一化径向平均功率谱；

根据所述归一化径向平均功率谱确定截止波数；所述根据所述归一化径向平均功率谱确定截止波数，具体包括：根据所述归一化径向平均功率谱的曲线特征，选取所述归一化径向平均功率谱下降预设比例a所对应的波数为截止波数cK；

根据所述截止波数对所述补零后重磁位场数据进行去噪、插值和扩边处理，获得处理后的重磁数据；所述根据所述截止波数对所述补零后重磁位场数据进行去噪、插值和扩边处理，获得处理后的重磁数据，具体包括：对所述补零后重磁位场数据作二维傅里叶变换，获得所述补零后重磁位场数据的频谱G(u,v)；

根据所述截止波数确定低通滤波函数 其中ck表示第k次迭代的截止波数，1≤k≤K且满足c1≤c2≤…≤cK和1＜ck≤min(M,N)，M和N是插值扩边后重磁位场数据的尺寸；D(u,v)＝[(u-M/2)2+(v-N/2)2]1/2；

采用所述低通滤波函数Tk(u,v)对所述频谱G(u,v)进行低通滤波处理，获得低通滤波后留下的谱分量；

对所述谱分量作二维傅里叶反变换，获得反变换后的重磁位场数据；

利用所述反变换后的重磁位场数据重置所述补零后重磁位场数据，获得重置后的重磁位场数据；

判断迭代次数k是否不大于预设迭代次数K，获得第一判断结果；

若所述第一判断结果为所述迭代次数k不大于所述预设迭代次数K，令k＝k+1，返回所述根据所述截止波数确定低通滤波函数的步骤；

若所述第一判断结果为所述迭代次数k大于所述预设迭代次数K，输出所述重置后的重磁位场数据作为处理后的重磁数据。

2.根据权利要求1所述的重磁数据处理方法，其特征在于，所述利用所述反变换后的重磁位场数据重置所述补零后重磁位场数据，获得重置后的重磁位场数据，具体包括：将所述补零后重磁位场数据中没有数据的位置用所述反变换后的重磁位场数据替换，有数据的位置保留，获得所述重置后的重磁位场数据。

3.一种重磁数据处理系统，其特征在于，所述系统包括：重磁异常数据获取模块，用于获取重磁异常数据；

补零模块，用于将所述重磁异常数据中有数据缺失的位置补零，获得补零后重磁位场数据；

归一化径向平均功率谱计算模块，用于计算所述补零后重磁位场数据的归一化径向平均功率谱；所述归一化径向平均功率谱计算模块具体包括：重磁位场数据功率谱求解单元，用于采用快速傅里叶变换方法求解所述补零后重磁位场数据的重磁位场数据功率谱；

环带生成单元，用于以所述重磁位场数据功率谱的中心为圆心，以基频的整数倍为半径作圆，生成多个环带；每个所述环带对应一个不同的频率；

平均功率谱值求解单元，用于求解各个所述环带内所有点的功率谱值的平均值，获得多个平均功率谱值；

径向平均功率谱生成单元，用于以径向波数为横坐标，以各个所述环带内的所述平均功率谱值的对数为纵坐标，生成径向平均功率谱；

归一化径向平均功率谱生成单元，用于对所述径向平均功率谱作归一化处理，生成所述重磁位场数据的归一化径向平均功率谱；

截止波数确定模块，用于根据所述归一化径向平均功率谱确定截止波数；所述截止波数确定模块具体包括：截止波数确定单元，用于根据所述归一化径向平均功率谱的曲线特征，选取所述归一化径向平均功率谱下降预设比例a所对应的波数为截止波数cK；

重磁数据处理模块，用于根据所述截止波数对所述补零后重磁位场数据进行去噪、插值和扩边处理，获得处理后的重磁数据；所述重磁数据处理模块具体包括：傅里叶变换单元，用于对所述补零后重磁位场数据作二维傅里叶变换，获得所述补零后重磁位场数据的频谱G(u,v)；

低通滤波函数确定单元 ，用于根据所述截止波数确定低通滤波函数其中ck表示第k次迭代的截止波数，1≤k≤K且满足c1≤c2≤…≤cK和1＜ck≤min(M,N)，M和N是插值扩边后重磁位场数据的尺寸；D(u,v)＝[(u-M/2)2+(v-N/2)2]1/2；

低通滤波单元，用于采用所述低通滤波函数Tk(u,v)对所述频谱G(u,v)进行低通滤波处理，获得低通滤波后留下的谱分量；

傅里叶反变换单元，用于对所述谱分量作二维傅里叶反变换，获得反变换后的重磁位场数据；

数据重置单元，用于利用所述反变换后的重磁位场数据重置所述补零后重磁位场数据，获得重置后的重磁位场数据；

第一判断单元，用于判断迭代次数k是否不大于预设迭代次数K，获得第一判断结果；

迭代计算单元，用于当所述第一判断结果为所述迭代次数k不大于所述预设迭代次数K时，令k＝k+1，返回所述根据所述截止波数确定低通滤波函数的步骤；

重磁数据输出单元，用于当所述第一判断结果为所述迭代次数k大于所述预设迭代次数K时，输出所述重置后的重磁位场数据作为处理后的重磁数据。

4.根据权利要求3所述的重磁数据处理系统，其特征在于，所述数据重置单元具体包括：数据重置子单元，用于将所述补零后重磁位场数据中没有数据的位置用所述反变换后的重磁位场数据替换，有数据的位置保留，获得所述重置后的重磁位场数据。

说明书

一种重磁数据处理方法及系统

技术领域

[0001] 本发明涉及重力场和地磁场数据测量技术领域，特别是涉及一种重磁数据处理方法及系统。

背景技术

[0002] 实际重磁测量的测区经常是不规则的，获得的重磁数据往往存在空缺且含有明显的噪声干扰。在波数域中进行重磁数据处理和转换前必须对原始数据进行填充和扩边处理，且高频噪声干扰是造成一些转换不稳定的主要原因。当前，涉及波数域重磁数据处理过程中，大都假设输入的重磁数据是纯粹的数学矩阵或向量，既要数据无空白，又要满足快速傅里叶变换(fast Fourier transform，FFT)对数据长度的要求。显然，实测数据很难满足上述条件。因此在波数域处理前，必须对数据空白部分进行填充、对数据进行扩边，这将直接关系到数据的处理精度。常规重磁数据处理一般将数据的去噪、填充和扩边这三个步骤独立进行，各个步骤之间相互制约，费时费力，操作不便且缺乏统一性。

发明内容

[0003] 本发明的目的是提供一种重磁数据处理方法及系统，能够同时实现重磁数据的去噪、填充和扩边，从而提高重磁数据处理效率。

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

[0005] 一种重磁数据处理方法，所述方法包括：

[0006] 获取重磁异常数据；

[0007] 将所述重磁异常数据中有数据缺失的位置补零，获得补零后重磁位场数据；

[0008] 计算所述补零后重磁位场数据的归一化径向平均功率谱；

[0009] 根据所述归一化径向平均功率谱确定截止波数；

[0010] 根据所述截止波数对所述补零后重磁位场数据进行去噪、插值和扩边处理，获得处理后的重磁数据。

[0011] 可选的，所述计算所述补零后重磁位场数据的归一化径向平均功率谱，具体包括：

[0012] 采用快速傅里叶变换方法求解所述补零后重磁位场数据的重磁位场数据功率谱；

[0013] 以所述重磁位场数据功率谱的中心为圆心，以基频的整数倍为半径作圆，生成多个环带；每个所述环带对应一个不同的频率；

[0014] 求解各个所述环带内所有点的功率谱值的平均值，获得多个平均功率谱值；

[0015] 以径向波数为横坐标，以各个所述环带内的所述平均功率谱值的对数为纵坐标，生成径向平均功率谱；

[0016] 对所述径向平均功率谱作归一化处理，生成所述重磁位场数据的归一化径向平均功率谱。

[0017] 可选的，所述根据所述归一化径向平均功率谱确定截止波数，具体包括：

[0018] 根据所述归一化径向平均功率谱的曲线特征，选取所述归一化径向平均功率谱下降预设比例a所对应的波数为截止波数cK。

[0019] 可选的，所述根据所述截止波数对所述补零后重磁位场数据进行去噪、插值和扩边处理，获得处理后的重磁数据，具体包括：

[0020] 对所述补零后重磁位场数据作二维傅里叶变换，获得所述补零后重磁位场数据的频谱G(u,v)；

[0021] 根据所述截止波数确定低通滤波函数 其中ck表示第k次迭代的截止波数，1≤k≤K且满足c1≤c2≤…≤cK和1＜ck≤min(M,N)，M和N是插值扩边后重磁位场数据的尺寸；D(u,v)＝[(u-M/2)2+(v-N/2)2]1/2；

[0022] 采用所述低通滤波函数Tk(u,v)对所述频谱G(u,v)进行低通滤波处理，获得低通滤波后留下的谱分量；

[0023] 对所述谱分量作二维傅里叶反变换，获得反变换后的重磁位场数据；

[0024] 利用所述反变换后的重磁位场数据重置所述补零后重磁位场数据，获得重置后的重磁位场数据；

[0025] 判断迭代次数k是否不大于预设迭代次数K，获得第一判断结果；

[0026] 若所述第一判断结果为所述迭代次数k不大于所述预设迭代次数K，令k＝k+1，返回所述根据所述截止波数确定低通滤波函数的步骤；

[0027] 若所述第一判断结果为所述迭代次数k大于所述预设迭代次数K，输出所述重置后的重磁位场数据作为处理后的重磁数据。

[0028] 可选的，所述利用所述反变换后的重磁位场数据重置所述补零后重磁位场数据，获得重置后的重磁位场数据，具体包括：

[0029] 将所述补零后重磁位场数据中没有数据的位置用所述反变换后的重磁位场数据替换，有数据的位置保留，获得所述重置后的重磁位场数据。

[0030] 本发明还公开了一种重磁数据处理系统，所述系统包括：

[0031] 重磁异常数据获取模块，用于获取重磁异常数据；

[0032] 补零模块，用于将所述重磁异常数据中有数据缺失的位置补零，获得补零后重磁位场数据；

[0033] 归一化径向平均功率谱计算模块，用于计算所述补零后重磁位场数据的归一化径向平均功率谱；

[0034] 截止波数确定模块，用于根据所述归一化径向平均功率谱确定截止波数；

[0035] 重磁数据处理模块，用于根据所述截止波数对所述补零后重磁位场数据进行去噪、插值和扩边处理，获得处理后的重磁数据。

[0036] 可选的，所述归一化径向平均功率谱计算模块具体包括：

[0037] 重磁位场数据功率谱求解单元，用于采用快速傅里叶变换方法求解所述补零后重磁位场数据的重磁位场数据功率谱；

[0038] 环带生成单元，用于以所述重磁位场数据功率谱的中心为圆心，以基频的整数倍为半径作圆，生成多个环带；每个所述环带对应一个不同的频率；

[0039] 平均功率谱值求解单元，用于求解各个所述环带内所有点的功率谱值的平均值，获得多个平均功率谱值；

[0040] 径向平均功率谱生成单元，用于以径向波数为横坐标，以各个所述环带内的所述平均功率谱值的对数为纵坐标，生成径向平均功率谱；

[0041] 归一化径向平均功率谱生成单元，用于对所述径向平均功率谱作归一化处理，生成所述重磁位场数据的归一化径向平均功率谱。

[0042] 可选的，所述截止波数确定模块具体包括：

[0043] 截止波数确定单元，用于根据所述归一化径向平均功率谱的曲线特征，选取所述归一化径向平均功率谱下降预设比例a所对应的波数为截止波数cK。

[0044] 可选的，所述重磁数据处理模块具体包括：

[0045] 傅里叶变换单元，用于对所述补零后重磁位场数据作二维傅里叶变换，获得所述补零后重磁位场数据的频谱G(u,v)；

[0046] 低通滤波函数确定单元，用于根据所述截止波数确定低通滤波函数其中ck表示第k次迭代的截止波数，1≤k≤K且满足c1≤c2≤…≤cK和1＜ck≤min(M,N)，M和N是插值扩边后重磁位场数据的尺寸；D(u,v)＝[(u-M/2)2+(v-N/2)2]1/2；

[0047] 低通滤波单元，用于采用所述低通滤波函数Tk(u,v)对所述频谱G(u,v)进行低通滤波处理，获得低通滤波后留下的谱分量；

[0048] 傅里叶反变换单元，用于对所述谱分量作二维傅里叶反变换，获得反变换后的重磁位场数据；

[0049] 数据重置单元，用于利用所述反变换后的重磁位场数据重置所述补零后重磁位场数据，获得重置后的重磁位场数据；

[0050] 第一判断单元，用于判断迭代次数k是否不大于预设迭代次数K，获得第一判断结果；

[0051] 迭代计算单元，用于当所述第一判断结果为所述迭代次数k不大于所述预设迭代次数K时，令k＝k+1，返回所述根据所述截止波数确定低通滤波函数的步骤；

[0052] 重磁数据输出单元，用于当所述第一判断结果为所述迭代次数k大于所述预设迭代次数K时，输出所述重置后的重磁位场数据作为处理后的重磁数据。

[0053] 可选的，所述数据重置单元具体包括：

[0054] 数据重置子单元，用于将所述补零后重磁位场数据中没有数据的位置用所述反变换后的重磁位场数据替换，有数据的位置保留，获得所述重置后的重磁位场数据。

[0055] 根据本发明提供的具体实施例，本发明公开了以下技术效果：

[0056] 本发明提供一种重磁数据处理方法及系统，所述方法及系统将重磁数据的去噪、填充和扩边这三个问题统一考虑，基于凸集投影方法同时实现了重磁数据的去噪、填充和扩边处理。所述方法及系统首先对原始重磁异常数据补零，并通过计算归一化径向平均功率谱来确定迭代过程最终的截止波数，然后对重磁数据进行去噪、填充和扩边迭代直至达到预定的迭代次数，获得处理后的重磁数据。本发明提供的方法及系统原理简单、操作方便、数据处理效率高，同时具有较高的去噪、填充精度，且填充和扩边的衔接处光滑无畸变，非常适合重磁数据的同时去噪和扩充。

附图说明

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

[0058] 图1为本发明提供的重磁数据处理方法的方法流程图；

[0059] 图2为本发明提供的对补零后重磁位场数据进行去噪、插值和扩边处理的方法流程图；

[0060] 图3为本发明提供的重磁数据处理系统的系统结构图；

[0061] 图4为本发明实施例提供的1km高度理论重磁异常数据的等值线图；

[0062] 图5为本发明实施例提供的重磁加噪数据等值线图；

[0063] 图6为本发明实施例提供的含空缺数据的重磁异常数据的等值线图；

[0064] 图7为本发明实施例提供的补零后重磁位场数据的归一化径向平均功率谱；

[0065] 图8为本发明实施例提供的经去噪和扩充处理后得到的重磁数据的等值线图；

[0066] 图9为本发明实施例提供的三种均方误差随迭代次数的变化图；

[0067] 图10为本发明实施例提供的去噪结果和理论值的残差等值线图；

[0068] 图11为本发明实施例提供的扩充结果和理论值残差等值线图。

具体实施方式

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

[0070] 本发明的目的是提供一种重磁数据处理方法及系统，能够同时实现重磁数据的去噪、填充和扩边，从而提高重磁数据处理效率。

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

[0072] 重磁数据的去噪、填充和扩边对应于信号的重构和外推。信号的重构和外推是信号处理的逆问题，具有不适定的共性特点。因此，可以考虑将重磁数据的去噪、填充和扩边统一考虑。凸集投影(Projection onto Convex Sets，POCS)方法是带限信号重构的重要理论方法之一，在图像重建和地震数据插值等领域应用广泛，且具有原理简单、易于执行和精度较高的特点。对于实际工程应用中感兴趣的重磁数据，往往是通过实际测量获得的重磁数据减去由模型计算的正常场数据而获得的重磁异常数据，因此，实测重磁异常数据可以视为带限信号，满足采用凸集投影方法处理的条件。本发明基于凸集投影原理，将重磁数据的去噪、填充和扩边这三个不适定问题统一考虑，提出一种基于凸集投影的重磁数据处理方法，能够实现重磁数据同时去噪、插值和扩边。

[0073] 图1为本发明提供的重磁数据处理方法的方法流程图。参见图1，本发明提供的一种重磁数据处理方法，具体包括：

[0074] 步骤101：获取重磁异常数据。

[0075] 所述重磁异常数据是由实际测量获得的重磁数据减去由模型计算的正常场数据而获得的重磁异常数据，可以视为带限信号，满足采用凸集投影方法处理的条件。

[0076] 步骤102：将所述重磁异常数据中有数据缺失的位置补零，获得补零后重磁位场数据。

[0077] 将含缺失数据且噪声水平为δ的重磁异常数据 有数据缺失的位置补零。 表示实际测量获得的重磁异常数据，obs是观测英文单词observation的缩写，δ代表噪声水平，x、y代表x轴和y轴上的坐标点。

[0078] 步骤103：计算所述补零后重磁位场数据的归一化径向平均功率谱。

[0079] 所述补零后重磁位场数据的归一化径向平均功率谱求解步骤如下：

[0080] 步骤(a)：采用快速傅里叶变换(FFT)方法求解所述补零后重磁位场数据的重磁位场数据功率谱。功率谱为幅度谱的平方。

[0081] 步骤(b)：以所述重磁位场数据功率谱的中心为圆心，以基频的整数倍为半径作圆，生成多个环带。

[0082] 假设基频为1，半径为基频1倍即半径为1，2倍基频即半径为2，半径为1的圆与半径为2的圆之间的部分即为环，各整数倍之间形成所述环带。每个所述环带对应一个不同的频率。

[0083] 步骤(c)：求解各个所述环带内所有点的功率谱值的平均值，获得多个平均功率谱值。

[0084] 步骤(d)：以径向波数为横坐标，以各个所述环带内的所述平均功率谱值的对数为纵坐标，生成径向平均功率谱。

[0085] 其中所述径向波数等于

[0086] 步骤(e)：对所述径向平均功率谱作归一化处理，生成所述重磁位场数据的归一化径向平均功率谱。

[0087] 此处对所述径向平均功率谱作归一化处理，是为了与后面的低通滤波函数Tk(u,v)对应起来，因为阈值函数Tk(u,v)的最大值也为1，这样方便确定截止波数。

[0088] 步骤104：根据所述归一化径向平均功率谱确定截止波数。

[0089] 根据所述归一化径向平均功率谱的曲线特征，选取所述归一化径向平均功率谱下降预设比例a所对应的波数为截止波数cK。

[0090] 在确定最优正则参数cK后，后续选取的阈值ck按线性方式从小递增到最大值cK。假设cK为20，迭代过程中所述阈值ck即由1，2，3，4...一直递增到20。其中阈值ck表示第k次迭代的截止波数，1≤k≤K且满足c1≤c2≤…≤cK。

[0091] 步骤105：根据所述截止波数对所述补零后重磁位场数据进行去噪、插值和扩边处理，获得处理后的重磁数据。

[0092] 图2为本发明提供的对补零后重磁位场数据进行去噪、插值和扩边处理的方法流程图。参见图2，所述步骤105具体包括：

[0093] 步骤S1：对所述补零后重磁位场数据g(x,y)作二维傅里叶(Fourier)变换，获得所述补零后重磁位场数据的频谱G(u,v)。

[0094] 步骤S 2：根据所述截止波数确定低通滤波函数其中ck表示第k次迭代的截止波数，1≤k≤K且满足c1≤
c2≤…≤cK和1＜ck≤min(M,N)。M和N是插值扩边后重磁位场数据的尺寸；例如所述补零后重磁位场数据数据是200*200大小，扩边成了256*256大小，则M和N分别为256和256。

[0095] D(u,v)表示波数域中点(u,v)与波数矩形中心的距离，即：

[0096] D(u,v)＝[(u-M/2)2+(v-N/2)2]1/2。

[0097] 步骤S3：采用所述低通滤波函数Tk(u,v)对所述频谱G(u,v)进行低通滤波处理，获得低通滤波后留下的谱分量。

[0098] 步骤S4：对低通滤波后留下的所述谱分量作二维傅里叶反变换，获得反变换后的重磁位场数据。

[0099] 步骤S5：利用所述反变换后的重磁位场数据重置所述补零后重磁位场数据，获得重置后的重磁位场数据。

[0100] 实际测量的重磁异常数据中有的位置有数据，有的位置没有数据，本步骤中重置数据即是将原来有数据的位置还采用原来的数据，没有数据的位置用变换后的数据进行重置。具体为：将所述补零后重磁位场数据中没有数据的位置用所述反变换后的重磁位场数据替换，有数据的位置保留，获得所述重置后的重磁位场数据。

[0101] 步骤S6：判断迭代次数k是否不大于预设迭代次数K，获得第一判断结果。

[0102] 步骤S7：若所述第一判断结果为所述迭代次数k不大于所述预设迭代次数K，令k＝k+1，返回所述步骤S2根据所述截止波数确定低通滤波函数的步骤。按照所述步骤S2～S6的迭代过程继续迭代，以便于恢复出缺失数据的更多细节信息。

[0103] 每次迭代过程中均更新阈值ck，阈值ck按线性方式从小递增到最大值cK，从而能够使更多的谱分量通过。

[0104] 步骤S8：若所述第一判断结果为所述迭代次数k大于所述预设迭代次数K，表示迭代次数达到预设迭代次数K，迭代过程结束，输出所述重置后的重磁位场数据作为处理后的重磁数据，即获得去噪、插值、扩边后的重磁数据。

[0105] 可见，本发明将重磁数据的去噪、填充和扩边处理这三个问题统一考虑，提出了一种基于凸集投影方法的重磁数据同时去噪、填充和扩边迭代方法。该迭代法首先对原始数据补零，并通过计算归一化径向平均功率谱来确定迭代法最终的截止波数，然后对数据进行去噪、填充和扩边迭代直至达到预定的迭代次数。该迭代方法原理简单、操作方便，具有较高的去噪、填充精度，且填充和扩边的衔接处光滑无畸变，非常适合重磁数据的同时去噪和扩充。

[0106] 本发明还提供了一种重磁数据处理系统，图3为本发明提供的重磁数据处理系统的系统结构图。参见图3，所述系统包括：

[0107] 重磁异常数据获取模块301，用于获取重磁异常数据；

[0108] 补零模块302，用于将所述重磁异常数据中有数据缺失的位置补零，获得补零后重磁位场数据；

[0109] 归一化径向平均功率谱计算模块303，用于计算所述补零后重磁位场数据的归一化径向平均功率谱；

[0110] 截止波数确定模块304，用于根据所述归一化径向平均功率谱确定截止波数；

[0111] 重磁数据处理模块305，用于根据所述截止波数对所述补零后重磁位场数据进行去噪、插值和扩边处理，获得处理后的重磁数据。

[0112] 其中，所述归一化径向平均功率谱计算模块303具体包括：

[0113] 重磁位场数据功率谱求解单元，用于采用快速傅里叶变换方法求解所述补零后重磁位场数据的重磁位场数据功率谱；

[0114] 环带生成单元，用于以所述重磁位场数据功率谱的中心为圆心，以基频的整数倍为半径作圆，生成多个环带；每个所述环带对应一个不同的频率；

[0115] 平均功率谱值求解单元，用于求解各个所述环带内所有点的功率谱值的平均值，获得多个平均功率谱值；

[0116] 径向平均功率谱生成单元，用于以径向波数为横坐标，以各个所述环带内的所述平均功率谱值的对数为纵坐标，生成径向平均功率谱；

[0117] 归一化径向平均功率谱生成单元，用于对所述径向平均功率谱作归一化处理，生成所述重磁位场数据的归一化径向平均功率谱。

[0118] 所述截止波数确定模块304具体包括：

[0119] 截止波数确定单元，用于根据所述归一化径向平均功率谱的曲线特征，选取所述归一化径向平均功率谱下降预设比例a所对应的波数为截止波数cK。

[0120] 所述重磁数据处理模块305具体包括：

[0121] 傅里叶变换单元，用于对所述补零后重磁位场数据作二维傅里叶变换，获得所述补零后重磁位场数据的频谱G(u,v)；

[0122] 低通滤波函数确定单元，用于根据所述截止波数确定低通滤波函数其中ck表示第k次迭代的截止波数，1≤k≤K且满足c1≤c2≤…≤cK和1＜ck≤min(M,N)，M和N是插值扩边后重磁位场数据的尺寸；D(u,v)＝[(u-M/2)2+(v-N/2)2]1/2；

[0123] 低通滤波单元，用于采用所述低通滤波函数Tk(u,v)对所述频谱G(u,v)进行低通滤波处理，获得低通滤波后留下的谱分量；

[0124] 傅里叶反变换单元，用于对所述谱分量作二维傅里叶反变换，获得反变换后的重磁位场数据；

[0125] 数据重置单元，用于利用所述反变换后的重磁位场数据重置所述补零后重磁位场数据，获得重置后的重磁位场数据；

[0126] 第一判断单元，用于判断迭代次数k是否不大于预设迭代次数K，获得第一判断结果；

[0127] 迭代计算单元，用于当所述第一判断结果为所述迭代次数k不大于所述预设迭代次数K时，令k＝k+1，返回所述根据所述截止波数确定低通滤波函数的步骤；

[0128] 重磁数据输出单元，用于当所述第一判断结果为所述迭代次数k大于所述预设迭代次数K时，输出所述重置后的重磁位场数据作为处理后的重磁数据。

[0129] 所述数据重置单元具体包括：

[0130] 数据重置子单元，用于将所述补零后重磁位场数据中没有数据的位置用所述反变换后的重磁位场数据替换，有数据的位置保留，获得所述重置后的重磁位场数据。

[0131] 本发明提供的重磁数据处理系统基于凸集投影方法实现了重磁数据的同时去噪、填充和扩边处理，原理简单、操作方便，极大地提高了重磁数据处理效率。

[0132] 下面采用本发明提供的重磁数据处理方法及系统对长方体磁异常模型进行仿真实验。

[0133] 理论模型采用长方体磁场无解析奇点模型。长方体的长宽高分别为8km、4km、4km，中心埋深3km，磁化强度1A·m-1，磁化方向为偏角30°倾角50°，地磁场偏角10°倾角60°。

[0134] 正演计算距离地面1km高空处的256×256大小的磁异常数据，如图4所示。图4为本发明实施例提供的1km高度理论重磁异常数据的等值线图。图5为本发明实施例提供的重磁加噪数据等值线图。图6为本发明实施例提供的含空缺数据的重磁异常数据的等值线图。图7为本发明实施例提供的补零后重磁位场数据的归一化径向平均功率谱。图8为本发明实施例提供的经去噪和扩充处理后得到的重磁数据的等值线图。图9为本发明实施例提供的三种均方误差随迭代次数的变化图。图10为本发明实施例提供的去噪结果和理论值的残差等值线图。图11为本发明实施例提供的扩充结果和理论值残差等值线图。图4-图11的横坐标和纵坐标均表示距离，单位为km。图中nT为磁异常单位。

[0135] 为模拟实际情况并检验本发明重磁数据处理方法的抗噪能力，给图4所示的重磁异常数据增加零均值、方差为0.5nT的高斯白噪声(信噪比为32.74dB)，其结果如图5所示。为检验本发明方法的填充和扩边能力，将图5所示的重磁异常数据内部和边界“挖去”部分数据构成空白区见图6所示。图6中四个边界各去除了28个数据，内部去掉的数据分布选取在正异常梯度带、正负异常交接处及边部异常平缓处，具有一定的典型性。采用本发明提供的重磁数据处理方法，对图6所示的加噪且有缺失的重磁异常数据进行处理，并以图4所示的真实重磁异常数据来检验本发明方法的去噪、填充和扩边精度。

[0136] 根据图1所示的算法步骤，先将图6所示的数据空白处补零，再计算补零后重磁异常数据的归一化径向平均功率谱，如图7所示。根据归一化径向平均功率谱的变化特征，选择其下降0.5所对应的波数为截止波数cK＝12。选定迭代次数K＝50次，按照图1所示的方法对补零后重磁异常数据进行迭代去噪、填充和扩边处理，获得处理后的重磁数据，图6所示数据的去噪、填充和扩边的结果见图8。由图8和原始重磁异常数据图4的对比可知，采用本发明方法得到的去噪、填充的结果同理论值非常一致，扩边的结果稍差，但也保持了形态特征，且已知数据和空白数据连接处光滑无畸变。

[0137] 为定量分析本发明方法的数据处理精度，统一由均方误差(Rmse)来定量分析去噪、填充和扩边结果的误差。所述均方误差(Rmse)计算公式为：

[0138]

[0139] 其中，Dc(i)和Dt(i)分别表示计算值和理论值，Q表示计算误差的总数据量。去噪、填充和扩边的均方误差随迭代次数k的变化如图9所示。去噪和填充的结果与理论值的残差分别见图10和图11。由图9可知，本发明提供的方法具有收敛性，最终的去噪、填充和扩边结果和理论值之间的误差统计见表1：

[0140] 表1去噪、填充和扩边的误差统计表(nT)

[0141]  误差最小值 误差最大值 均方误差
去噪 -0.70 0.62 0.09
填充 -2.22 3.68 0.77
扩边 -10.88 11.22 3.68

[0142] 由表1可见，本发明提供的重磁数据处理方法的去噪和填充具有较高的精度。

[0143] 综上可见，本发明提供的重磁数据处理方法及系统与常规重磁数据处理方法相比，至少具有以下优点：

[0144] 1、本发明提供的重磁数据处理方法及系统考虑了实测重磁位场数据存在的数据n缺失和数据长度不符合快速傅里叶变换对数据长度(2)的要求这两个问题，并采用插值方法解决了数据缺失的问题，采用扩边处理解决了数据长度问题。比如将200*200扩边至256*
256，满足了快速傅里叶变换对数据长度最好为2n的要求。

[0145] 2、本发明提供的方法及系统还利用了重磁位场数据去噪、插值和扩边同属不适定反问题的共性特点，建立了这三个问题的一体化解法方案，采用迭代法实现了它们的同步解决。

[0146] 3、本发明提供的方法及系统只需预先采用归一化径向平均功率谱确定截止波数和预定义迭代次数，采用理想低通滤波器进行迭代，整个迭代过程的参数较少、算法操作简便、复杂度低，能够大大提高数据处理效率。

[0147] 4、仿真结果表明，本发明提供的方法及系统具有很高的去噪和插值精度，且插值和扩边的衔接处光滑无畸变。

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

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