0 引言

实测重力数据不仅含有大量噪声，且可能由于各种因素的限制，常常会存在部分缺失。波数域重力数据处理过程中，要求输入的重力数据无空白，且数据长度也要尽量满足快速Fourier变换对数据长度(2的幂次方)的要求。显然，实测数据很难满足上述条件。因此，在进行波数域数据处理前，必须对重力数据空白部分进行填充、对边界进行扩边，这将直接关系到数据的处理精度^[1-2]。

常用的重力数据去噪方法有Fourier域滤波法、优化滤波法、经验模态分解法、等效源法、小波分解法、滑动窗口平均法等。Fourier域滤波法的基本原理是在波数域提取并分离干扰信号，其主要问题是信号和噪声成份不易区分。由此，基于重力位场物理特性的等效层法^[3]、等效源法^[4]和基于奇异值分解、经验模态分解^[5]等的信号处理方法被引入重力数据去噪中，其目的就是有效区分有用信号和噪声。小波分解法^[6]在各类信号的去噪中都得到了广泛应用，但小波去噪的效果与小波基函数以及阈值的选择有较大关系。

对缺失数据进行填充重建，理论上所有的数据插值方法都可以使用，经典的插值方法有线性插值、样条插值、Kriging法、最小曲率法等。这些方法均有各自的理论适应性和计算特点，它们的操作共同点是：首先将测量数据与坐标结合，以网格中已知数据点为型值点，空缺位置坐标为待插值点，进行插值^[7]。

扩边实质上是信号的外推过程，常用的扩边方法有补零扩边法、对折扩边法、区域场扩边法、三方向扩边法、余弦扩边法和最小曲率法等^[1-2]。余弦扩边方法不能很好地反映实际数据的区域变化特征，并在实际数据和扩边数据的衔接处导数不连续，从而产生较严重的Gibbs效应^[1]。最小曲率法是目前重力数据插值和扩边中最受欢迎的方法，但该方法主要考虑曲面的光滑性(在没有数据点的情况下具有最小曲率表面的是球形)，因此，不能得到精确的插值结果。当网格化间距小于实际的网格间距时，稀疏数据控制区域的插值精度会快速下降^[8]。

凸集投影(Projection onto Convex Sets，POCS)方法是带限信号重构的重要理论方法之一，在图像重建^[9]和地震数据插值^[10-15]等领域应用广泛，且具有原理简单、易于执行和精度较高的特点。在实际工程应用中，往往是通过实际测量获得的重力数据减去由模型计算的正常场数据而获得重力异常数据。因此，实测重力异常数据可以视为带限信号^[16]，满足采用凸集投影方法处理的条件。在重力数据处理方面，闫浩飞等^[7]利用POCS方法实现了重力数据插值重建；曾小牛等^[17]提出一种基于凸集投影原理的重力数据填充、扩边、下延一体化方法。Ji等^[18]结合预条件共轭梯度迭代法和POCS方法实现了磁场参数的正则化反演。本文基于凸集投影原理，将重力数据的填充、去噪和扩边这三个问题统一考虑，提出基于凸集投影的规则网格重力数据同时填充、扩边和去噪方法。

1 方法原理 1.1 基于POCS的重力数据同时扩充和去噪

设含噪且存在数据缺失的观测重力网格数据为g_obs，则基于凸集投影的填充和扩边计算可表示为

$ \begin{array}{c} \boldsymbol{g}_{k}=\boldsymbol{g}_{\mathrm{obs}}+(\boldsymbol{I}-\boldsymbol{S}) \mathrm{F}^{-1} \boldsymbol{T}_{k} \mathrm{F} \boldsymbol{g}_{k-1} \\ k=1, 2, 3, \cdots, K \end{array} $

(1)

式中：g_k表示第k次迭代扩边填充后的数据；I为单位矩阵；S为采样矩阵，由0和1组成，1代表该点有数据，无需插值，反之为0，则代表该点需要进行插值；F和F^-1分别表示Fourier正、反变换；K表示最大迭代次数；T_k表示阈值矩阵，其元素满足

$ T_{k}(u, v)=\left\{\begin{array}{ll} 1 & \left|S_{k-1}(u, v)\right| \geqslant p_{k} \\ 0 & \left|S_{k-1}(u, v)\right| <p_{k} \end{array}\right. $

(2)

其中：S_k-1表示第k-1次迭代得到的扩充数据g_k-1的频谱；u、v分别是x、y方向的波数；p_k∈{p₁, p₂, …, p_K}表示第k次迭代的阈值，满足p₁＞p₂＞…＞p_K和min{|S_k-1(u, v)|}＜p_k＜max {|S_k-1(u, v)|}。从式(1)迭代过程可知，POCS方法是利用整个数据的谱(能量)进行扩充。

显然，式(1)的每次迭代都将含噪声的原始数据g_obs带入了插值结果，因此，传统的POCS方法对含噪数据的插值效果不理想。为此，曹静杰等^[12]提出了一种用于地震数据重建的改进凸集投影方法，本文将其应用到重力数据波数域的同时去噪和扩充，可表示为

$ \begin{array}{c} \boldsymbol{g}_{k}=\mathrm{F}^{-1} \boldsymbol{T}_{k} \mathrm{F}\left[\boldsymbol{g}_{\text {obs }}+(\boldsymbol{I}-\boldsymbol{S}) \boldsymbol{g}_{k-1}\right] \\ k=1, 2, 3, \cdots, K \end{array} $

(3)

由式(1)和式(3)可知，该类迭代方法的核心是确定波数域阈值p_k。在基于POCS方法的地震和图像信号去噪重建中，一般先确定阈值p_k的最大值和最小值，再由线性^[10]、指数或数据驱动^[11]模式完成从最大值到最小值的下降过程。实际应用中，p_k的最大值和最小值的确定与原始数据的噪声水平有关，对结果影响较大，需要谨慎选择^[10-12]。相对地震和图像信号来说，重力数据变化较平缓，连续性较好，频谱S_k-1(u, v)相对集中。T_k(u, v)的实质是理想低通滤波器，针对重力数据频谱特点及频谱的物理特征，本文将T_k(u, v)修改为

$ T_{k}(u, v)=\left\{\begin{array}{ll} 1 & D(u, v) \leqslant c_{k} \\ 0 & D(u, v)>c_{k} \end{array}\right. $

(4)

式中：c_k表示第k次迭代的截止波数，满足c₁≤c₂≤…≤c_K和1＜c_k≤min(M, N)，M和N是扩充后重力数据的大小；D(u, v)表示波数域中点(u, v)与波数矩形中心的距离

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

(5)

显然，式(4)与式(2)是等价的。式(4)只是将式(2)的阈值p_k转化成截止波数c_k。同理，最大截止波数c_K的大小与原始数据噪声水平相关。

1.2 基于径向平均功率谱的最大截止波数c_K确定方法

在基于POCS方法的地震和图像信号去噪、重建过程中，阈值最大值和最小值的确定比较随意，缺乏一定的准则^[10-12]。对应于式(4)的理想低通滤波器，截止波数c_K显然最好为有效信号与噪声的分界点，这样迭代过程既能保留有效信号应用于扩充和填充，又能避免噪声的干扰。

Spector等^[19]提出的径向平均功率谱能够反映重磁位场的频谱特征，并被广泛应用于重磁位场数据的处理^{[3, 20-22]}。一个假想的径向平均功率谱如图 1所示，则垂直叠加的深源和浅源产生重力异常的径向平均功率谱可以采用下式进行拟合

$ S(u, v)=A^{2} \mathrm{e}^{-2 h \sqrt{u^{2}+v^{2}}} $

(6)

式中：A为能谱幅值；h为源深度。如果观测噪声为白噪，则与重力信号不相关的白噪声功率谱应该为常数，即噪声大致对应于径向平均功率谱的水平部分。这样，对于重力异常的径向平均功率谱，存在一个截止波数将重力异常的信号谱与噪声谱大致分开。由此，可将迭代法的最大截止波数c_K设置为由径向平均功率谱拟合所确定的划分信号谱与噪声谱的截止波数。显然，这样的设定具有明确的物理意义，即只采用有效信号的频谱进行扩充，同时又能消除噪声谱。

通过径向平均功率谱可获得截止波数c_K，结合迭代次数K，迭代过程的截止波数c_k依线性递增的模式从1增大到最终的截止波数c_K。本文所提的迭代方法可归纳为图 2。

2 理论模型实验

理论模型由位于不同深度、不同大小的4个长方体组合而成(图 3a)，各长方体的参数见表 1。

点、线距均设为50m，平面剖分网格大小为256×256，正演计算地面(z=0)的重力响应如图 3b所示。为检验方法的去噪能力，给图 3b所示的重力数据增加零均值、标准差为1.00mGal的高斯白噪声(图 3c)。为检验方法的填充和扩边能力，首先，随机去掉图 3c所示数据总量的1/4，并在内部和边界“挖去”部分数据，得到缺失数据(图 3d)。图 3c内部“挖去”的数据选取在异常梯度带及边部异常平缓处，这样的选择具有一定的典型性。采用本文提出的扩边、填充及去噪法，对图 3d所示的加噪且有缺失的重力数据进行处理，并与图 3b所示的真实重力数据进行对比。采用均方根误差(Root Mean Square Error，RMSE)进行定量分析

$ \mathrm{RMSE}=\sqrt{\frac{1}{Q} \sum\limits_{i=1}^{Q}\left[D_{\mathrm{c}}(i)-D_{t}(i)\right]^{2}} $

(7)

式中：D_c(i)和D_t(i)分别表示计算值和理论值；Q表示参与误差计算的总数据量。

根据图 2的算法步骤，选定最大迭代次数K，先计算图 3d重力数据的径向平均功率谱(图 3e)。根据形状特征将功率谱分为三段：1~5(频段1)、5~18(频段2)和18~129(频段3)。频段1主要对应低频异常信息，频段2主要对应中高频信息，而频段3的功率谱主要对应高斯白噪声。因此，确定截止波数c_K为频段2和频段3的分界点17。最终的填充、去噪和扩边结果见图 4a中的上图。对比图 4a中的上图与原始重力数据(图 3b)可知，去噪、填充的结果同理论值非常接近，扩边的结果稍差，但形态特征非常一致，且已知数据与缺失数据的连接处光滑无畸变。

为了对比分析，分别采用经典的最小曲率法和Kriging插值法分别与小波分解去噪、余弦函数扩边法组成两组组合方法对图 3d所示数据分步进行插值、去噪和扩边，并与本文方法的处理结果进行对比。最小曲率和Kriging插值法均可在Golden Surfer 12软件中实现：最小曲率插值法的最大残差和最大迭代次数参数分别选定为1.0×10^-5和1.0×10⁵；Kri-ging插值的变异函数模型选定为Gaussian函数。小波分解采用Matlab提供的小波分析工具，分别试验了常用的Symlets小波系中的sym1~sym8小波和Daubechies小波系中db1~db8小波的去噪效果，选定去噪结果与理论值均方根误差最小的sym4小波进行最终的去噪实验；余弦函数扩边法采用Cooper等^[23]提供的“taper2d”程序。

采用这两组组合方法对图 3d所示数据按填充、去噪、扩边的顺序进行试验。图 4展示了本文方法填充、去噪和扩边残差。最终结果与理论值(图 3b)相比较的填充、去噪和扩边的均方根误差分别为：组合方法1为0.86、0.03、4.80mGal；组合方法2为0.22、0.03、4.80mGal。图 5是填充、去噪和扩边均方根误差随迭代次数k的变化曲线。由图可知，算法具有收敛性，最终的填充、去噪和扩边均方根误差分别为0.08、0.01、0.82mGal。由图 4和图 5可知，本文方法的填充、去噪和扩边结果与理论值更一致，残差和均方根误差也更小。

3 实测重力数据对比实验

为检验本文方法在实际数据处理中的实用性，采用美国地质调查局(United States Geological Survey，USGS)公布的阿富汗均衡重力异常数据进行对比实验。该重力数据融合了2006年的航空实测数据以及地面已有的重力数据，经各项处理并最终归算至离地7km的高度、数据网格距为1km的重力异常(图 6a)。由图可见，融合后的重力数据存在较大面积的空白区。

3.1 局部完整数据对比实验

首先从原始数据中截取256km×256km大小的完整数据(图 6a黑色方框区域)进行方法验证，图 6b所示即是所选区域的数据。同时，因原始数据的高频噪声经过低通滤波，为模拟实际情况，对图 6a数据添加零均值、标准差为1.00mGal的高斯白噪声，结果如图 6c。与理论模型实验一样，为了检验本文所提方法的填充和扩边能力，将图 6c所示方框内数据随机去掉数据总量的1/4，并在边界和内部“挖去”部分数据构成空白区(图 6d)。

采用本文提出的扩充去噪迭代法对图 6d数据进行处理，并对比图 6b所示的真实重力数据以检验方法的填充、去噪和扩边精度。根据图 2所示算法步骤，首先计算图 6d所示数据的径向平均功率谱，如图 7所示。根据功率谱形状，将其分为1~6(频段1)、6~18(频段2)、18~49(频段3)和49~129(频段4)四个频段。频段1主要对应低频异常信息，频段2主要对应中高频信息，而频段3和频段4主要对应高斯白噪声。因此，确定截止波数c_K为频段2和频段3的分界点18。选定迭代次数K=200，然后按照图 2所示的算法步骤进行迭代计算。

与理论模型实验过程一样，采用经典的最小曲率法和Kriging插值法分别与小波分解去噪、余弦函数扩边法组成两组组合方法，对图 6d所示的数据分步进行填充、去噪和扩边，并与本文方法的处理结果(图 8a上图)进行对比。经对比Symlets小波系和Daubechies小波系试验，小波分解采用使去噪均方根误差最小的sym8小波进行去噪；两种插值法和余弦函数扩边的实现和参数均与理论模型一致。组合方法1和2的填充、扩边和去噪结果分别见图 8b上图和图 8c上图。组合方法1处理结果(图 8b上图)相比于图 6b的填充、去噪和扩边均方根误差分别为3.36、0.20、11.75mGal；组合方法2处理结果(图 8c上图)相比于图 6b的填充、去噪和扩边均方根误差分别为1.97、0.19、11.56mGal。图 8中图展示了图 8上图与图 6b相比的填充和去噪残差，图 8下图展示了图 8上图与图 6b相比的扩边残差。填充、去噪和扩边的均方根误差随迭代次数k的变化见图 9。由图可知，本文方法具有收敛性，最终的填充、去噪和扩边的均方根误差分别为1.22、0.15、5.37mGal。对比图 8所显示的本文方法与经典组合方法的填充、去噪和扩边结果，以及残差和均方根误差，可以看出，同理论数据计算结果一样，本文方法的填充、去噪和扩边结果优于组合方法，其残差和均方根误差也均小于组合方法。

3.2 完整数据对比实验

为更进一步对比方法在数据缺失较多情况下的效果，采用提出的迭代法和两组组合方法对图 6c所示的重力数据进行对比实验。首先计算图 6c所示数据的径向平均功率谱(图 10)。基于对功率谱形状的分析，将其分为1~12(频段1)、12~82(频段2)、82~270(频段3)和270~512(频段4)四段。频段1主要对应低频异常信息，频段2和频段3主要对应中高频信息，而频段4主要对应高斯白噪声。因此，确定截止波数c_K为频段3与频段4的分界点270。选定迭代次数K=200。本文方法的填充、去噪和扩边结果见图 11a。

与前两次实验过程类似，采用经典的最小曲率法和Kriging插值法与小波分解去噪对图 6c所示的数据分步进行填充、扩边和去噪。为保证收敛，采用最小曲率法插值时，插值和扩边的参数分别确定为1.0×10^-2和1.0×10⁵；Kriging插值法的变异函数依然为Gaussian函数。去噪方面，经对比Symlets小波系和Daubechies小波系试验结果，小波分解均采用使去噪均方根误差最小的sym8小波进行去噪。组合方法1填充、扩边和去噪的结果见图 11b，小波去噪的均方根误差为0.30mGal。组合方法2填充、扩边和去噪的结果见图 11c，小波去噪的均方根误差为0.28mGal。

对比图 11a与图 11b、图 11c可知，虽然采用最优小波去噪的均方根误差略小于本文方法，但本文方法的填充和扩边结果基本保持了真实数据的形态特征，已知数据和空白数据的连接处光滑无畸变，而最小曲率法的填充和扩边结果仅在靠近已知数据的边缘形态较好，在离已知数据较远的区域则效果很差。Kriging插值法的填充结果较好，但在离已知数据较远的区域扩边效果同样不甚理想。

图 12为去噪均方根误差随迭代次数k的变化，最终的去噪均方根误差为0.32mGal，可见算法具有收敛性。

4 结论

本文针对实际重力测量获得的重力异常数据往往存在噪声且有空缺的实际情况，基于带限信号重建的经典算法——凸集投影算法，提出了一种重力数据同时填充、扩边和去噪的迭代法，并利用理论模型和实测数据验证了方法的收敛性和有效性。从算法流程来看，迭代法采用FFT运算、在波数域仅涉及理想低通滤波器，因此只需要预先设定迭代次数以及依据拟合径向平均功率确定截止波数即可完成迭代，具有原理简单、实际操作方便和运算高效的特点。从理论和实测数据处理的实验效果来看，在扩边和补空的衔接处，数据光滑无畸变，方法取得了较好的填充、扩边和去噪效果。本文方法的处理结果优于由最小曲率法和Kriging插值与小波分解去噪、余弦扩边方法组成的组合方法。