海底地形测量在海洋测量中占据了极其重要的地位, 是开发、利用和保护海洋资源与环境的基础性工作。作为主要的方法之一, 多波束测深在海底地形测量中应用非常广泛。由于仪器自噪声、海况因素或声呐参数设置不合理等因素的影响, 多波束测深数据不可避免地存在异常数据, 造成虚假地形, 如何去除这些异常数据, 是多波束测深数据处理的一项重要任务。去除多波束测深异常数据, 可用手动方法, 也可用自动滤波方法, 后者近年来发展迅速。常见的自动滤波方法有常规统计检测法^[1-2]、抗差M估计法^[3]、趋势面滤波法^[4-5]、CUBE算法^[6-7]、遗传算法^[8]等。其中, 趋势面滤波方法原理简单、计算效率高, 在海量多波束测深数据滤波处理中有着广阔的应用前景。该方法利用多元回归原理, 用水深数据拟合一个趋势面, 根据拉依达准则, 将残差大于3σ的测深值看作异常值, σ为测点局部残差确定的均方根误差^[9]。趋势面滤波方法在海底地形变化简单且异常数据较少的情况下效果较好。然而, 由于多项式曲面拟合获得的趋势面包含了异常数据信息, 估计的系数向量会与最佳趋势面的系数向量产生偏差, 使构造的趋势面在异常数据较多的情况下并不准确; 此外, 异常数据也参与了σ的计算, 在异常数据较多的情况下, 3σ准则常常无法很好地检测出较小的异常。

针对传统趋势面滤波对异常数据的抗差性不足的问题, 一些学者提出了改进的趋势面滤波方法, 例如截断最小二乘估计^[10]、基于抗差估计的选权迭代趋势面拟合方法^[11]、高崩溃污染率抗差最小二乘估计^[12]等, 但是这些方法原理比较复杂、计算过程耗时较多, 失去了趋势面滤波方法在计算效率上的优势。

本文提出了一种中值滤波加权修正的方法改进了传统趋势面滤波方法, 原理简单、计算效率较高, 适用于海量多波束测深数据的滤波, 具有较强的实用价值。方法采用中值滤波加权修正的方式, 改进了趋势面滤波方法中的多项式拟合曲面的系数向量和测深点局部残差均方根误差这两个重要参数:以趋势面拟合优度的变化量作为迭代终止的原则, 不断修正由中值滤波找到的异常水深数据, 使求取的多项式拟合曲面系数向量接近真实海底地形拟合的系数向量; 此外, 采用修正后的水深数据求取测深点局部残差均方根误差σ, 极大减弱了异常数据对滤波阈值的影响。参考国内外相关文献, 目前还没有关于中值滤波加权修正的改进趋势面滤波方法的介绍, 本文提出的方法具有较强的创新性。

1 方法原理 1.1 传统趋势面滤波方法

设水深数据为$z({x_i}, {y_i})$, i=1, 2, …, n, 构造一个多项式趋势面函数$\hat z({x_i}, {y_i})$来拟合海底。以平面坐标正北方向为X轴、正东方向为Y轴, 趋势面可以表示为:

$ \begin{array}{l} \hat z({x_i}, {y_i}) = {a_p}{y^k} + {a_{p-1}}x{y^{k-1}} + \cdots + {a_{p-k}}{x^k} + \cdots \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\; + {a_5}{y^2} + {a_4}xy + {a_3}{x^2} + {a_2}y + {a_1}x + {a_0} \end{array} $

(1)

式中, $p = k(k + 3)/2$, ${a_0}, {a_1}, \cdots, {a_p}$是待定常数。将多项式回归的模型转换为线性回归的模型, 令: ${x_1} = x$, ${x_2} = y$, ${x_3} = {x^2}$, ${x_4} = xy$, ${x_5} = {y^2}$, $\cdots $, 则:

$ \hat z({x_i}, {y_i}) = {a_p}{x_p} + {a_{p-1}}{x_{p-1}} + \cdots + {a_1}{x_1} + {a_0} $

(2)

为了使整个趋势面在N个观测点上最好地逼近观测数据, 采用最小二乘法求取多项式的系数, 在N个观测点上使观测值与趋势值之差的平方和最小, 即:

$ Q = \min (\sum\limits_{i = 1}^n {{{[z({x_i}, {y_i})-\hat z({x_i}, {y_i})]}^2}} ) $

(3)

令:

$ {\boldsymbol{X}} = \left[{\begin{array}{*{20}{c}} 1&{{x_{11}}}&{{x_{21}}}&{...}&{{x_{p1}}}\\ 1&{{x_{12}}}&{{x_{22}}}&{...}&{{x_{p2}}}\\ {...}&{...}&{...}&{...}&{...}\\ 1&{{x_{1n}}}&{{x_{2n}}}&{...}&{{x_{pn}}} \end{array}} \right], {\boldsymbol{Z}} = \left[{\begin{array}{*{20}{c}} {{z_1}}\\ {{z_2}}\\ {...}\\ {{z_n}} \end{array}} \right], {\boldsymbol{A}} = \left[{\begin{array}{*{20}{c}} {{a_0}}\\ {{a_1}}\\ {...}\\ {{a_p}} \end{array}} \right] $

(4)

趋势面滤波的拟合阶次$k$通常取在10阶以内^[13], 对于变化比较简单的海底地形, 拟合阶次$k$通常取在3, 4阶以内便能满足要求, 参与计算的测深点数$n$远大于$p$($p = k(k + 3)/2$), 属于典型的超定方程组求解问题, 根据最小二乘的原理, 系数向量A的最小二乘解可表示为:

$ \mathit{\boldsymbol{A}} = {({\mathit{\boldsymbol{X}}^{\rm{T}}}\mathit{\boldsymbol{X}})^{ - 1}}{\mathit{\boldsymbol{X}}^{\rm{T}}}\mathit{\boldsymbol{Z}} $

(5)

趋势面分析通过一组水深数据拟合一个趋势函数, 这个函数从总体上反映海底地形的区域性变化趋势, 称为趋势面, 水深实测值与该点处对应的趋势面上的值之差称为残差, 记为$v({x_i}, {y_i})$, 传统趋势面滤波将$v({x_i}, {y_i})$大于3σ的测深点予以滤除(σ为测点局部残差确定的均方根误差, σ的倍数可以根据实际情况选取, 一般取为2σ或3σ)。然而, 式(5)中系数向量A以及测深点的局部残差均方根误差σ的求取都受到异常数据的影响, 异常数据较多的情况下, 传统趋势面滤波只能滤除偏离正常数据较大的异常数据, 对相对较小的异常数据不敏感。本文采用中值滤波加权修正的方式, 改进了趋势面滤波方法中的多项式拟合曲面系数向量和测深点局部残差均方根误差这两个重要参数的求取, 具有很强的抗差性。

1.2 结合中值滤波加权修正的趋势面滤波方法 1.2.1 中值滤波原理

中值滤波^[14]有一维中值滤波和二维中值滤波两种。滤波过程为:首先确定一个一维或二维窗口, 然后将其遍历序列或平面数据上的所有点, 对窗口中的数据进行排序, 以窗口数据的中值作为窗口中心位置的滤波结果。多波束测深数据的中值滤波属于二维中值滤波, 可以表示为:

$ z({u_i}, {v_i}) = {\rm{median}}\{ {z_{\rm{o}}}({u_i}, {v_i})\} $

(6)

式中, $({u_i}, {v_i})$为二维数据平面位置, ${z_{\rm{o}}}({u_i}, {v_i})$为以当前点$z({u_i}, {v_i})$为中心的窗口范围内的值, ${\rm{median\{ }}.{\rm{\} }}$表示求集合${\rm{\{ }}.{\rm{\} }}$的中值。

1.2.2 中值滤波窗口的选择

中值滤波的窗口越大, 对异常数据的滤波效果越好, 但造成的平滑作用越明显, 因此要合理选取窗口的大小。当多波束系统采用等角模式时, 测深数据在边缘波束区域分布稀疏, 在中央波束区域分布密集, 因此可以在中央波束区域设置较大的窗口, 在边缘波束区域设置较小的窗口, 中间波束区域的窗口设置介于两者之间。例如, 可以在中央波束区域设置6×6的窗口, 在边缘波束区域设置4×4的窗口, 在中间波束区域设置5×5的窗口。窗口表示沿着航迹方向同波束号的若干波束个数与垂直航迹方向的同ping的若干波束个数。

1.2.3 抗差趋势面的构造

中值滤波作为一种非线性信号处理技术, 虽然能够有效滤除异常数据, 但是会造成细节的非线性扭曲, 使得到的海底地形失真, 为了将中值滤波与趋势面滤波的优点结合, 如果把中值滤波后的结果用来构造趋势面, 同样面临着构造的趋势面是由非线性扭曲后的海底地形所拟合的问题。事实上当多波束数据采集质量较差时, 往往存在大量簇状异常, 采用多次中值滤波才能将异常平滑, 但是多次中值滤波必然会造成更严重的海底地形非线性扭曲, 因此本文提出了一种加权水深修正的方法。水深修正的原则是对偏离中值较大的水深数据作较大的修正, 使之向中值靠近; 对偏离中值较小的水深数据作较小的修正或基本不予修正, 使修正后的数据接近真实的海底地形数据。采用多次修正的方式, 以前后两次修正后数据的拟合优度变化量为迭代是否进行的依据, 完成修正后的水深数据极大降低了异常值的影响, 以修正后的水深数据求取的多项式拟合曲面系数向量接近真实海底水深的多项式拟合曲面系数向量。

设原始多波束测深数据为$z({x_i}, {y_i})$, 其中, $({x_i}, {y_i})$是测深数据的地理位置, $z({x_i}, {y_i})$为对应位置的水深, 抗差趋势面的构造过程如下:

(1)设定窗口大小, 对数据进行中值滤波。第${\rm{Me}}{{\rm{d}}^{(k)}}({x_i}, {y_i})$次中值滤波后的水深数据表示为: $k = 1, 2, 3, \cdots $, 表示迭代次数, 以下相同。

(2) 令${D^{(k)}}({x_i}, {y_i}) = z({x_i}, {y_i}) - {\rm{Me}}{{\rm{d}}^{(k)}}({x_i}, {y_i})$, 得到水深数据与该位置第$k$次中值滤波后的数据之差。第一次迭代中, $z({x_i}, {y_i})$是原始水深数据, 之后的迭代中, 表示上一次修正后的水深值。

(3) 定义“偏离度”变量${D_{\rm{g}}}$, 表示为:

$ {D_{\rm{g}}}^{(k)}({x_i}, {y_i}){\rm{ = }}\frac{{\left| {{D^{(k)}}({x_i}, {y_i})} \right|}}{{\max \left| {({D^{(k)}}(x, y))} \right|}} $

(7)

式中, max表示求最大值。${D_{\rm{g}}}^{(k)}({x_i}, {y_i})$反映了某位置水深值相对于该处中值的偏离量与局部区域内最大偏离量之比, 显然, ${D_{\rm{g}}}^{(k)}({x_i}, {y_i})$的取值范围为0~1。

(4) 加权水深修正:对偏离该处中值较大的数据进行较大幅度的修正, 对偏离该处中值较小的数据基本不予修正, 本文引入幂函数作为水深修正量的权值。设 $ {z^{(k)}}_{{\rm{cor}}}({x_i}, {y_i})$为第$k$次迭代中修正后的水深数据。则有:

${z^{(k)}}_{{\rm{cor}}}({x_i}, {y_i}) = z({x_i}, {y_i})-{D^{(k)}}({x_i}, {y_i})S $

(8)

式中, $S = {({D_{\rm{g}}}^{(k)}({x_i}, {y_i}))^\gamma }$是以“偏离度”为变量的幂函数。如图 1所示:当${D_{\rm{g}}}$较小时, $S$的值接近为0, 即对偏离中值较小的水深值基本不予修正; 对于相同的${D_{\rm{g}}}$, γ越大, $S$越小, 即γ决定了偏离中值较大的数据向中值收敛的快慢。为了对偏离中值较小的水深值基本不予修正, 可以选择一个较大的γ值。

(5) 选取合适的阶次, 对修正后的水深数据${z^{(k)}}_{{\rm{cor}}}({x_i}, {y_i})$进行多项式最小二乘曲面拟合, 求出拟合优度${R^2}$。拟合优度是测定回归模型拟合优劣的重要指标^[15], 公式如下:

$ {R^2}(k) = 1-\frac{{\sum\limits_{j = 1}^n {{{({z_j}-{{\hat z}_j})}^2}} }}{{\sum\limits_{j = 1}^n {{{({z_j}-\bar z)}^2}} }} $

(9)

式中, ${z_j}$为第j点处第$k$次迭代修正后的水深数据, $\bar z$为平均值, ${\hat z_j}$为j点处趋势面上的值。

(6) 由于异常数据的存在对曲面的拟合优度有较大的影响, 当$\left| {{R^2}(k)-{R^2}(k-1)} \right| \ge \varepsilon $时, 有必要继续进行水深修正, 令$z({x_i}, {y_i}) = {z^{(k)}}_{{\rm{cor}}}({x_i}, {y_i})$, 重复步骤(1)—(5)。当$\left| {{R^2}(k)-{R^2}(k-1)} \right| < \varepsilon $时, 前后两次水深修正数据的拟合优度接近, 此时可认为修正后的水深数据已经接近真实的海底水深数据, 此时的${R^2}(k)$或${R^2}(k-1)$接近真实水深数据的拟合优度, 故执行步骤(7)。

(7) 终止迭代, 设此时的迭代次数为N, 完成水深修正的数据为${z^{(N)}}_{{\rm{cor}}}({x_i}, {y_i})$, 求出${z^{(N)}}_{{\rm{cor}}}({x_i}, {y_i})$的多项式拟合曲面系数向量, 由系数向量和测深点平面位置得到最终的趋势面${z_{{\rm{fit}}}}({x_i}, {y_i})$。

不同于传统趋势面滤波直接由原始测深数据拟合得到趋势面, 本文提出的方法由完成加权修正后的水深数据拟合得到趋势面, 构造的趋势面具有很强的抗差性。

1.2.4 抗差阈值的求取与滤波

据原始水深数据和最终的抗差趋势面得到残差$v({x_i}, {y_i})$:

$ v({x_i}, {y_i}) = z({x_i}, {y_i}) - {z_{{\rm{fit}}}}({x_i}, {y_i}) $

(10)

式中的$z({x_i}, {y_i})$为原始测深数据。根据修正后的水深数据与最终的抗差趋势面得到残差$v'({x_i}, {y_i})$:

$ v'({x_i}, {y_i}) = {z^{(N)}}_{{\rm{cor}}}({x_i}, {y_i})-{z_{{\rm{fit}}}}({x_i}, {y_i}) $

(11)

利用$v'({x_i}, {y_i})$求出局部残差均方根误差σ′, 异常数据的检测阈值取为3σ′(可根据实际情况选取σ′若干倍作为阈值, 一般取为3σ′即可), 将$v({x_i}, {y_i})$中超过3σ′的测深点视为异常数据并予以滤除。

传统趋势面滤波中, 先求出原始测深数据和以原始测深数据拟合得到的趋势面之间的残差, 再求出残差的局部均方根误差, 以一定倍数(通常为3倍)的局部残差均方根误差作为滤波的阈值; 而本文的局部残差均方根误差由完成加权修正的水深数据和根据该数据拟合的趋势面之差求取, 最终用以滤波的阈值具有很强的抗差性。

2 仿真与实测多波束测深数据滤波试验 2.1 仿真海底数据模拟试验

为了验证本文提出的滤波方法的准确性, 先利用仿真模拟的方法进行检验。采用三次多项式曲面模拟真实海底地形, 并向其中加入均值为0, 方差为0.01的高斯噪声模拟海底地形的随机变化, 模拟海底地形数据共10 000个, 向其中分别加入100, 300, 500, 1 000个偏离真实海底2.0~10.0 m的异常数据, 对比本文提出的方法和传统趋势面滤波方法的效果, 两种滤波方案都采用三阶多项式曲面进行拟合, 阈值分别取为3σ′和3σ(此处仅为了对比两种方法的抗差性)。本文提出的方法其他参数为: ε取为0.01, γ取为30。统计结果见表 1, 图 2a、图 2b、图 2c和图 2d分别为加入1000个异常数据时的原始数据、基于中值滤波加权水深修正后的数据、按本文提出的方法滤波后的数据以及传统趋势面滤波之后的数据示意图。

由表 1可知, 传统趋势面滤波只有在异常数据较少时有较好的效果, 当异常个数增多, 其抗差性不足的缺点得到明显体现, 滤波效果很差, 当异常数据在总数据中的比例超过3%时, 使用传统趋势面滤波方法有大量异常数据不能被滤除。本文提出的滤波方法的效果抗差性良好, 在异常数据占总数据10%的情况下, 依然能够准确地滤除异常数据。由图 2b可知, 加权水深修正后的测深数据已经很接近真实的海底地形, 利用修正后的水深数据构造趋势面与求取局部残差均方根误差, 大幅降低了异常数据对系数向量和滤波阈值的影响。

值得注意的是, 仿真模拟的真实的海底地形由三次多项式曲面和高斯噪声组成, 在此假设的基础上, 本文提出的方法在存在大量异常数据的情况下仍然能够取得很好的效果, 但真实的海底地形必然更加复杂, 因此, 需要进行真实的多波束测深数据的滤波检验。

2.2 实测数据试验及分析

试验选取的数据来自舟山某海域的R2 Sonic2024多波束测深系统实测水深数据, 范围约为190 m× 160 m, 共158 720个离散水深数据。经过安装偏差校正、姿态改正、声速改正、潮位改正之后, 如图 3a所示, 发现其中含有大量异常数据, 造成的虚假地形严重影响了真实的海底DEM和等深线成图, 需要予以滤除。

该区域真实海底地形并不复杂, 但为了保证每个子区的真实地形能够被多项式曲面很好的拟合, 按照地理位置, 将测区划分成了8个子区。采用传统的趋势面滤波的方法对其进行滤波处理, 调节趋势面拟合的阶次和异常检测的阈值, 发现无论如何都不能将异常数据很好滤除, 得到的最佳滤波情况如图 3b所示。可以看到, 传统趋势面滤波之后, 测区边缘部分的小异常和测区中部的粗大异常被滤除, 但不能将测区中部较小的异常滤掉, 这是因为传统趋势面滤波方法抗差性不足, 导致趋势面系数向量与滤波阈值的求取均受到异常数据的影响, 在异常数据较多的情况下, 对相对较小的异常值不敏感。

采用本文提出的基于中值滤波加权修正的趋势面滤波方法, 基于MATLAB编程语言编写了相关程序, 其中的参数选取为:趋势面拟合的阶次取为3, γ取为30, ε取为0.01, 滤波门限取为3σ′。在CPU型号为Intel(R) Core(TM) i5-6300U、主频为2.4 GHz的环境下, 运行时间为6.738 s。滤波后的DEM和等深线如图 3c所示。从图 3c可以看出, 经本文提出的方法滤波后的DEM和等深线图基本没有因异常数据引起的虚假地形和虚假等深线, 能够准确反映真实的海底地形。

采用多波束数据处理软件Caris Hips手动剔除粗差, 作为实际海底地形的一种参考方案(实际勘测工作中, 多波束测深异常数据滤波仍以手动方式为主), 结果如图 3d所示。对比图 3c与图 3d, 两者无论是DEM还是等深线都很接近。然而, 在异常数据较多且包含许多较小的异常时, 手动滤波往往是一个很耗费时间的过程, 而且很容易受到操作人员自身判断的影响, 与手动滤波相比, 本文提出的方法在滤波效率上有明显的优势, 而且获得的滤波结果更加客观。

将趋势面滤波、本文方法和Caris Hips软件手动滤波后的结果展开统计分析。由表 2可知, 与手动滤波结果相比, 传统趋势面滤波之后的水深平均值偏小、方差偏大, 与图 3b中DEM显示的含有部分偏浅的异常数据的结果相符合。本文提出的方法滤波后的水深均值和方差都接近手动滤波的结果, 均值和方差略微偏小可能是因为滤除了更多手动滤波时难以发现的较小且偏深的异常数据引起的。

3 结论

本文提出的方法改进了传统的趋势面滤波方法, 虽然引入了中值滤波加权修正的步骤, 但是仍然能保证较高的滤波效率, 对158 720个离散数据进行滤波处理仅用6.738 s, 远高于手动滤波的效率。

基于中值滤波加权修正的方法构造的趋势面受异常数据的影响大幅降低, 相比传统的趋势面滤波, 该趋势面能较真实地反映正常水深数据的最小二乘拟合。

利用修正后的水深数据与趋势面的残差求取滤波阈值受异常数据的影响大幅降低, 相比传统的趋势面滤波, 滤波阈值的选取更加合理。

当海底真实地形变化比较简单且含有较多的异常数据时, 使用本文提出的方法能够很好地将异常数据滤除, 由于其强大的抗差性, 对较小的异常数据也有很好的滤波效果。但是, 同其他的自动滤波方法一样, 本文提出的方法也存在一定的局限:在局部起伏剧烈的复杂地区, 本文提出的方法滤波效果不佳, 也可能会平滑真实的地形, 因为这种情况下, 真实的海底很难被多项式曲面模拟; 当存在海底真实小目标时, 如探测海底沉船桅杆时, 本文提出的方法可能会将其检测为异常数据并滤除。因此, 在多波束测深数据滤波工作中, 当海底地形变化复杂或需要保留真实小目标时, 应结合手动滤波的方法进行处理。