留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

一种基于改进凸集投影原理的航空重力数据插值与去噪方法

曾小牛 李夕海 刘继昊 牛超 侯维君

曾小牛, 李夕海, 刘继昊, 牛超, 侯维君. 一种基于改进凸集投影原理的航空重力数据插值与去噪方法[J]. 武汉大学学报 ● 信息科学版, 2020, 45(10): 1555-1562. doi: 10.13203/j.whugis20180470
引用本文: 曾小牛, 李夕海, 刘继昊, 牛超, 侯维君. 一种基于改进凸集投影原理的航空重力数据插值与去噪方法[J]. 武汉大学学报 ● 信息科学版, 2020, 45(10): 1555-1562. doi: 10.13203/j.whugis20180470
ZENG Xiaoniu, LI Xihai, LIU Jihao, NIU Chao, HOU Weijun. Simultaneous Interpolation and Denoising Method for Airborne Gravity Data Based on Improved Projection onto Convex Sets Theory[J]. Geomatics and Information Science of Wuhan University, 2020, 45(10): 1555-1562. doi: 10.13203/j.whugis20180470
Citation: ZENG Xiaoniu, LI Xihai, LIU Jihao, NIU Chao, HOU Weijun. Simultaneous Interpolation and Denoising Method for Airborne Gravity Data Based on Improved Projection onto Convex Sets Theory[J]. Geomatics and Information Science of Wuhan University, 2020, 45(10): 1555-1562. doi: 10.13203/j.whugis20180470

一种基于改进凸集投影原理的航空重力数据插值与去噪方法

doi: 10.13203/j.whugis20180470
基金项目: 

国家自然科学基金 41804136

国家自然科学基金 41774156

国家自然科学基金 61773389

陕西省高校科协青年人才托举计划 20180702

详细信息
    作者简介:

    曾小牛, 博士, 讲师, 主要从事重力数据处理与解释方法研究。xiaoniuzeng@163.com

    通讯作者: 李夕海, 博士, 教授。xihai_li@163.com
  • 中图分类号: P223

Simultaneous Interpolation and Denoising Method for Airborne Gravity Data Based on Improved Projection onto Convex Sets Theory

Funds: 

The National Natural Science Foundation of China 41804136

The National Natural Science Foundation of China 41774156

The National Natural Science Foundation of China 61773389

the Young Talent Fund of University Association for Science and Technology in Shanxi, China 20180702

More Information
    Author Bio:

    ZENG Xiaoniu, PhD, lecturer, specializes in gravity data processing and interpretation. E-mail: xiaoniuzeng@163.com

    Corresponding author: LI Xihai, PhD, professor. E-mail: xihai_li@163.com
图(9) / 表(3)
计量
  • 文章访问数:  567
  • HTML全文浏览量:  112
  • PDF下载量:  31
  • 被引次数: 0
出版历程
  • 收稿日期:  2019-12-27
  • 刊出日期:  2020-10-05

一种基于改进凸集投影原理的航空重力数据插值与去噪方法

doi: 10.13203/j.whugis20180470
    基金项目:

    国家自然科学基金 41804136

    国家自然科学基金 41774156

    国家自然科学基金 61773389

    陕西省高校科协青年人才托举计划 20180702

    作者简介:

    曾小牛, 博士, 讲师, 主要从事重力数据处理与解释方法研究。xiaoniuzeng@163.com

    通讯作者: 李夕海, 博士, 教授。xihai_li@163.com
  • 中图分类号: P223

摘要: 航空实测重力数据不仅噪声含量大,而且由于各种因素的限制,数据常常需要插值加密。常规航空重力数据处理一般将数据的插值和去噪这两个步骤独立进行,将这两个问题统一考虑,借鉴地震数据处理方法,提出一种基于改进凸集投影原理的规则网格重力数据同时插值和去噪迭代方法。基于EGM2008重力场模型仿真一组重力数据,验证了该迭代法的插值和去噪效果优于4种经典插值方法和部分小波阈值去噪方法。并通过某区域的实测航空重力数据进一步验证了方法的有效性。

English Abstract

曾小牛, 李夕海, 刘继昊, 牛超, 侯维君. 一种基于改进凸集投影原理的航空重力数据插值与去噪方法[J]. 武汉大学学报 ● 信息科学版, 2020, 45(10): 1555-1562. doi: 10.13203/j.whugis20180470
引用本文: 曾小牛, 李夕海, 刘继昊, 牛超, 侯维君. 一种基于改进凸集投影原理的航空重力数据插值与去噪方法[J]. 武汉大学学报 ● 信息科学版, 2020, 45(10): 1555-1562. doi: 10.13203/j.whugis20180470
ZENG Xiaoniu, LI Xihai, LIU Jihao, NIU Chao, HOU Weijun. Simultaneous Interpolation and Denoising Method for Airborne Gravity Data Based on Improved Projection onto Convex Sets Theory[J]. Geomatics and Information Science of Wuhan University, 2020, 45(10): 1555-1562. doi: 10.13203/j.whugis20180470
Citation: ZENG Xiaoniu, LI Xihai, LIU Jihao, NIU Chao, HOU Weijun. Simultaneous Interpolation and Denoising Method for Airborne Gravity Data Based on Improved Projection onto Convex Sets Theory[J]. Geomatics and Information Science of Wuhan University, 2020, 45(10): 1555-1562. doi: 10.13203/j.whugis20180470
  • 航空实测重力数据往往存在噪声干扰, 且因测区的不规则常常导致数据需要插值和加密。经典重力数据去噪方法有Fourier频谱滤波法和小波阈值法等。Fourier频谱滤波法的基本原理是通过频谱分离的方式实现噪声干扰的去除。小波阈值法[1-2]在各类信号的去噪中都得到了广泛应用, 但小波去噪的效果也和小波基函数以及阈值的选择有较大关系。常用的重力数据插值方法有反距离加权插值法(又称Shepard法)、Kriging法、最小曲率法和径向基函数法等[3-5]。这些方法均有各自的理论适应性和计算特点, 它们的具体操作共同点是:首先将测量数据与坐标结合, 以网格中已知数据点为型值点, 以空缺位置坐标为待插值点, 计算完后再输出为网格形式以供后续计算。

    凸集投影(projection onto convex sets, POCS)方法是带限信号重构的重要理论方法之一, 在图像重建[6-7]和地震数据插值[8-11]等领域应用广泛, 具有原理简单、易于执行和精度较高的特点。对于实际应用中感兴趣的航空重力数据, 往往是通过实际测量获得的重力数据减去由模型计算的正常场数据而获得的重力异常数据, 因此, 航空实测重力异常数据可以视为带限信号[12], 满足采用凸集投影方法处理的条件。文献[13]首次提出将POCS方法应用于地形和其他地球物理场的插值;文献[14]利用POCS方法实现了重磁数据插值重建。本文将重力数据的插值和去噪问题统一考虑, 提出一种基于改进凸集投影原理的航空重力数据同时插值和去噪方法。

    • 设含缺失数据且噪声水平为δ的重力网格数据为$g_{{\rm{obs}}}^\delta \left( {x, y} \right)$, 则基于凸集投影的插值计算可表示为[14]:

      $$\begin{array}{*{20}{l}} {{\rm{}}\;\;{g_k}\left( {x, y} \right) = g_{{\rm{obs}}}^\delta \left( {x, y} \right) + }\\ {{\rm{}}\left( {\mathit{\boldsymbol{I}} - \mathit{\boldsymbol{S}}} \right){\mathit{\boldsymbol{F}}^{ - 1}}{\mathit{\boldsymbol{T}}_k}\mathit{\boldsymbol{F}}{g_{k - 1}}\left( {x, y} \right)} \end{array}$$ (1)

      式中, $k = 1, 2, 3 \cdots K$为迭代次数, K表示最大迭代次数;${g_k}\left( {x, y} \right)$表示第k次迭代插值后的数据;I为全1矩阵;S为采样矩阵, 由0和1组成, 如果某点为1, 代表该点有数据, 无需插值, 反之若为0, 则代表该点需要进行插值;FF-1分别表示Fourier正反变换;Tk表示阈值矩阵, 其元素满足:

      $${T_k}\left( {u, v} \right) = \left\{ {\begin{array}{*{20}{c}} {1, {\rm{}}\left| {{S_{k - 1}}\left( {u, v} \right)} \right| \ge {p_k}}\\ {0, {\rm{}}\left| {{S_{k - 1}}\left( {u, v} \right)} \right| < {p_k}} \end{array}} \right.$$ (2)

      式中, ${S_{k - 1}}\left( {u, v} \right)$表示第k-1次迭代得到的插值数据${g_{k - 1}}\left( {x, y} \right)$的频谱;${p_k} \in \left\{ {{p_1}, {p_2} \cdots {p_K}} \right\}$表示第k次迭代的阈值, 满足${\rm{min}}\left\{ {\left| {{S_{k - 1}}\left( {u, v} \right)} \right|} \right\} < {p_k} < {\rm{max}}\left\{ {\left| {{S_{k - 1}}\left( {u, v} \right)} \right|} \right\}$和${p_1} > {p_2} > \cdots > {p_K}$。

      显然, 公式(1)的每次迭代都将含噪声的原始数据$g_{{\rm{obs}}}^\delta \left( {x, y} \right)$带入了插值结果, 因此, 原始POCS方法对含噪声数据的插值效果不理想, 为此, 文献[9]提出了一种用于地震数据同时插值和去噪的改进凸集投影方法, 应用到重力数据, 可表示如下:

      $$\begin{array}{*{20}{l}} {{\rm{}}{g_k}\left( {x, y} \right) = {F^{ - 1}}{T_k}(F(g_{{\rm{obs}}}^\delta \left( {x, y} \right) + }\\ {{\rm{}}\;\;\;\;\;{\rm{}}\left( {I - S} \right){g_{k - 1}}\left( {x, y} \right)))} \end{array}$$ (3)

      式(1)和式(3)所示迭代法的核心是阈值pk的确定。在基于POCS法的地震和图像信号重建领域中, 一般先确定阈值pk的最大值和最小值, 再由线性、指数或数据驱动[8]模式完成从最大值到最小值的下降过程。实际应用中, 阈值pk的最大值和最小值的确定与原始数据的噪声水平有关, 对结果影响较大, 需要仔细选择[8-11]。相对地震和图像信号来说, 航空重力数据变化较平缓, 连续性较好, 频谱${S_{k - 1}}\left( {u, v} \right)$相对集中。${T_k}\left( {u, v} \right)$的实质是理想低通滤波器, 针对重力数据频谱特点及频谱的物理特征, 本文将${T_k}\left( {u, v} \right)$修改为:

      $${T_k}\left( {u, v} \right) = \left\{ {\begin{array}{*{20}{c}} {1, {\rm{}}D\left( {u, v} \right) \le {c_k}}\\ {0, {\rm{}}D\left( {u, v} \right) > {c_k}} \end{array}} \right.$$ (4)

      式中, ck表示第k次迭代的截止波数, 且满足${c_1} \le {c_2} \le \cdots \le {c_K}$和$1 < {c_k} \le {\rm{min}}\left( {M, N} \right)$, MN是插值后重力数据的尺寸;$D\left( {u, v} \right)$表示波数域中点$\left( {u, v} \right)$与波数矩形中心的距离, 即:

      $$D\left( {u, v} \right) = {[{\left( {u - M/2} \right)^2} + {\left( {v - N/2} \right)^2}]^{1/2}}$$ (5)

      显然, 式(4)和式(2)是等价的。式(4)只是将式(2)的阈值pk转化成了截止波数ck。同理, 最大截止波数cK的大小与原始数据的噪声水平是相关的。

    • 在基于POCS方法的地震和图像信号去噪、重建过程中, 阈值pk最大值和最小值的确定比较随意, 缺乏一定的准则[8-11]。对应于本文采用的公式(4)的理想低通滤波器, 截止波数cK显然最好为有用信号和噪声信号的分界点, 这样迭代过程既能保留有用信号不被削弱, 又能避免噪声对插值和去噪过程的干扰。

      文献[15]提出的径向平均功率谱能够反映重磁位场的频谱特征, 并被广泛应用于重磁数据处理[16-17]。一个假想的径向平均功率谱如图 1所示。

      图  1  径向平均功率谱示意图

      Figure 1.  Schematic Diagram for the Annular Averaged Power Spectrum

      如果假设观测噪声为白噪声, 则与重力信号不相关的白噪声功率谱应该为常数, 即噪声大致对应于径向平均功率谱的水平部分。因此, 对于重力数据的径向平均功率谱, 存在一个截止波数将重力数据信号谱和噪声谱大致分开。由此, 可将迭代法的最大截止波数cK设置为由径向平均功率谱所确定的截止波数。显然, 这样设定具有明确的物理意义, 即只采用有用信号进行插值, 同时又能消除噪声干扰。

      在通过对径向平均功率谱进行拟合获得截止波数cK以后, 运用二维快速傅立叶变换(2-dimensions fast Fourier transformation, 2D FFT), 结合迭代次数K, 迭代过程的截止波数ck依线性递增的模式从1增大到最终的截止波数cK。至此, 本文所提的迭代方法可归纳为图 2所示的过程。

      图  2  本文算法步骤示意图

      Figure 2.  Schematic Diagram of the Proposed Algorithm

    • 为验证本文提出的迭代法对重力数据同时去噪和插值的有效性, 基于EGM2008重力场模型计算某海域2 160阶的重力异常格网数据△g, 如图 3所示。试验区大小为1.05°×1.05°, 数据大小为64×64, 格网间隔为1'×1'。重力异常和地形数据的各项统计特性见表 1

      图  3  仿真的重力异常数据△g

      Figure 3.  Simulative Gravity Anomaly Data △g

      表 1  试验区域重力和地形数据统计特征

      Table 1.  Characteristics of Gravity Anomalies and Terrain Data in Test Area

      统计项 最小值 最大值 平均值 标准差
      重力/mGal -37.299 22.305 -9.099 13.360
      高程/m -3 560.430 -700.549 -2 532.543 605.167
    • 为模拟实际情况并检验方法的去噪能力, 给图 3所示的重力异常引入2种观测误差是零均值的白噪声, 标准差分别为1 mGal和3 mGal;同时, 为检验方法的插值能力, 将加噪后的数据随机抽出总数据量的1/5并在内部整体抽出两片位于异常梯度带的数据构成空白区, 如图 4(a)图 4(b)所示。

      图  4  本文方法同时插值和去噪示意图

      Figure 4.  Simultaneous Interpolation and Denoising Results for the Proposed Method

      为验证本文方法的效果, 设计以下5种插值方案:①采用反距离加权法进行插值;②采用Kriging法进行插值;③采用最小曲率法进行插值;④采用径向基函数法进行插值;⑤采用本文提出的迭代法同时插值和去噪。采用上述5种方案, 对图 4(a)4(b)所示的加噪且有缺失重力异常数据进行处理, 并以图 3所示的真实重力异常数据来检验各方案的插值精度。

      前4种方案的插值在Golden Surfer 12软件中完成, 各插值方法的参数可参考文献[18], 见表 2。根据图 2的算法步骤, 本文先计算图 4(a)图 4(b)所示重力数据的径向平均功率谱, 结果分别如图 4(c)图 4(d)所示。通过功率谱形状分析, 确定截止波数cK。选定迭代次数K次, 然后按照算法步骤进行迭代计算。插值和去噪的均方差随迭代次数k的变化见图 4(e)图 4(f)

      表 2  插值结果与原始重力异常△g的比较/mGal

      Table 2.  Comparison Results of Anomaly Gravity Differences Between Interpolated Result and Simulative Gravity Data △g/mGal

      统计项 方案 参数 最小值 最大值 平均值 标准差 均方差
      1 mGal
      误差
      权重:1
      平滑系数:0.01
      -5.628 9.820 0.205 1.368 2.722
      变异模型:线性
      Slope=1, Aniso=1, 0
      -3.070 8.408 0.143 0.838 1.672
      最大残差:1×10-6
      最大迭代次数:1×106
      -22.106 7.740 0.138 0.883 1.759
      基函数:multiquadric
      R2参数:0.05
      -2.721 8.478 0.148 0.856 1.709
      截止波数:13
      迭代次数:100
      -9.435 12.680 0.059 0.742 1.459
      3 mGal
      误差
      权重:1
      平滑系数:0.01
      -6.100 9.680 0.222 1.400 2.790
      变异模型:线性
      Slope=1, Aniso=1, 0
      -4.655 8.470 0.179 1.117 2.227
      最大残差:1×10-6
      最大迭代次数:1×106
      -27.086 9.493 0.183 1.386 2.751
      基函数:multiquadric
      R2参数:0.05
      -6.105 8.590 0.185 1.234 2.456
      截止波数:13
      迭代次数:100
      -10.895 12.586 0.118 0.937 1.849

      根据表 2, 从前4种插值计算方案中, 选择出插值结果误差最小的方案:方案②Kriging法插值结果进行小波阈值去噪试验。根据文献[1]对航空重力测量数据进行小波阈值滤波的对比研究结果, 设计如下4种去噪方案:①采用db8小波进行7层小波分解重构;②采用db10小波进行8层小波分解重构;③采用coif5小波进行8层小波分解重构;④本文方法。4种方案的去噪比较结果见表 3

      表 3  去噪结果与原始重力异常△g的比较/mGal

      Table 3.  Comparison Results of Anomaly Gravity Differences Between Denoising Result and Simulative Gravity Data △g/mGal

      统计项 方案 最小值 最大值 平均值 标准差 均方差
      1 mGal
      误差
      -6.206 8.934 0.223 1.742 1.472
      -5.926 8.794 0.072 1.700 1.442
      -4.254 7.599 0.005 1.532 1.311
      -11.257 12.456 -0.005 1.221 1.417
      3 mGal
      误差
      -8.160 10.035 0.355 2.649 2.420
      -8.637 9.595 0.077 2.640 2.416
      -6.507 8.054 0.072 2.227 2.063
      -13.245 12.177 0.030 1.592 1.848

      两种噪声条件下, 本文方法的插值和去噪最终结果见图 5(a)图 5(b), Kriging法插值后小波阈值去噪误差最小(去噪方案③)的结果分别见图 5(c)图 5(d)。4种插值和去噪结果与真实值(见图 3)的差值分布见图 6

      图  5  两种误差条件下本文方法和最优组合方法的插值和去噪结果

      Figure 5.  Interpolation and Denoising Results Given by the Proposed Method and the Optimal Combination Method Under the Conditions of the Two Error Levels

    • 1)由图 4(e)图 4(f)可知, 本文提出的迭代算法具有收敛性, 1 mGal噪声下的最终插值和去噪均方差分别为1.459 mGal和1.417 mGal, 3 mGal噪声下的最终插值和去噪均方差分别为1.849 mGal和1.848 mGal。方案①~④中插值误差最小的为方案②, 即Kriging插值法, 但本文方法(方案⑤)插值误差的平均值、标准差和均方差均小于其他4种插值方法, 验证了本文方法的插值精度较优。

      2)由表 3可知, 去噪方案③采用coif5小波进行8层小波分解阈值去噪的结果优于方案①和方案②两种小波阈值去噪结果, 而本文方法去噪误差的平均值和标准差均小于前3种方案, 只在1 mGal噪声下的误差均方差稍大于方案③, 验证了本文方法的去噪精度较优。

      3)由图 5所显示的本文方法和Kriging插值法与小波阈值去噪的组合方法最终结果对比可知, 本文方法的插值和去噪结果同理论值更一致, 进一步验证了本文方法的插值和去噪精度较优。

      4)从图 6的差值分布图可知, 本文方法差值的最大值和最小值主要分布在边界区域, 这是由于采用Fourier运算所产生的边界效应, 如果采用扩边处理, 将进一步降低误差。

      图  6  图 5的4种插值和去噪结果同理论值(图 3)的差值结果

      Figure 6.  Difference Between the Four Interpolation and Denoising Results of Fig. 5 and Fig. 3

    • 实测航空重力数据采用澳大利亚某区域的飞行测量数据, 试验区大小为0.75°×1.05°, 数据大小为24×35, 高程为655 m, 如图 7所示。

      图  7  实测重力异常数据

      Figure 7.  Measured Gravity Anomaly Data

      根据图 2的算法步骤, 本文先计算图 7所示重力数据的径向平均功率谱, 结果如图 8所示。

      图  8  实测重力异常数据的径向平均功率谱

      Figure 8.  Radial Average Power Spectrum of the Measured Gravity Anomaly Data

      通过功率谱形状分析, 确定截止波数cK。选定迭代次数K次, 然后按照算法步骤进行迭代计算。本文方法的插值和去噪最终结果见图 9(a)。采用插值精度最高的方案②和去噪效果最好的方案③组成组合方法对实测数据进行先插值后去噪, 其结果见图 9(b)

      图  9  本文方法和最优组合方法的插值和去噪结果

      Figure 9.  Interpolation and Denoising Results Given by the Proposed Method and the Optimal Combination Method

      图 9所显示的本文方法和最优组合方法插值及去噪结果与图 7的对比可知, 本文方法的插值和去噪结果同有数据区域的重力变化趋势更一致, 进一步验证了本文方法的插值和去噪精度较优。

    • 本文针对航空实测重力异常数据往往存在噪声干扰且有空缺的实际情况, 基于带限信号重建的经典算法——凸集投影算法, 提出了一种航空重力数据同时插值和去噪迭代法。从算法流程来看, 迭代法采用FFT运算, 在波数域仅涉及理想低通滤波器的使用, 只需要预先依据径向平均功率的拟合来确定截止波数和设定迭代次数即可完成迭代, 具有原理简单、实际操作方便和运算高效的特点。通过基于EGM2008重力场模型仿真的重力数据进行插值和去噪试验, 验证了本文提出的方法优于4种经典的插值方法和部分小波阈值去噪方法。并基于实测航空重力数据的对比试验进一步验证了方法的有效性。

参考文献 (18)

目录

    /

    返回文章
    返回