一种自动探测多波束测量异常值的组合算法

张秋明

（福建省港航管理局勘测中心（福建省港航勘察科技有限公司），福建 福州 350009）

0.引言

多波束测深系统可以有效地获得密集的三维空间数据，被广泛应用于海洋工程、海底形变监测和海洋资源开发等领域[1，2]。考虑到海洋环境的复杂性，如鱼群、仪器噪声和可变观测条件，MBS不能完全准确地获取实际的海底空间数据，在多波束测量数据中存在着异常值，导致其不能充分反映实际海底地形，给航行安全带来潜在威胁。为了满足国际海道测量规范[3]的要求，需要后处理方法以减少这些异常值的影响。

异常值的检测是多波束测深数据自动处理的关键因素。由于目前的软件不能有效地解决这一问题，常用的方法是人工编辑法。随着MBS的不断发展，采集数据集变得越来越密集，需要花费大量时间的人工编辑法变得更加笨拙[4]。研究人员已经提出了许多方法来解决这个问题。Shaw等[5]提出了一种自回归算法来检测轻微异常值。根据相邻点的情况，Du等[6]提出应用聚类技术对点进行分类。Mann等[7]应用中值滤波探测和去除浅水环境中的异常值。基于贝叶斯分析，Calder等[8]建立了CUBE算法，将稳健估计与实测数据进行比较，根据两者存在的差异大小确定异常值。Debese等[9]利用局部鲁棒表面估计和分层自适应方法对不同尺寸的海底网格进行分类。Rezvani等[2]应用M估计来衰减异常值的影响。本文结合将小波变换的定位能力与卡尔曼滤波的鲁棒动态预测相结合，提出了一种新的异常值自动检测算法，并利用数据对其性能进行检测。

1.基本原理

1.1 离散小波变换

小波变换是一种能同时进行局部分析和多分辨率分析的小波变换。母小波是各种小波变换的基础如公式（1）所示。通过尺度扩展和平移，可以得到一族小波函数如公式（2）所示：

其中，参数a是尺度因子，参数b是平移因子。离散小波变换（DWT）是对基本小波的尺度和平移的离散化。

DWT显示出非常适合于信号，因为它的时间尺度是由信号的大小决定的。当数据在不同频率下分解成一定的水平时，就会识别出隐藏的特征。

1.2 Kalman滤波

离散卡尔曼滤波的数学模型由状态方程和观测方程两部分组成。其离散形式可通过公式（3）表示：

其中，φk,k－1表示时间k-1到k的传递矩阵；Γk,k－1表示系统噪声矩阵；Hk表示k时刻的系统观测矩阵；Vk表示观测噪声；Xk表示状态变量的向量；sdk表示系统的观测向量矩阵。

2.人工编辑法

现今海洋测绘工程领域，绝大多数的工程技术人员采用的是caris软件进行内业处理。数据处理流程主要包括编辑船型配置文件、声速剖面改正、潮位改正、剔除异常值、CUBE生成、合并数据、按比例提取数据、图上展点和绘图等过程。编辑船型配置文件主要是使用各模块的相对位置信息，以及计算出来的校准参数等数据构造测船的模型；声速剖面改正主要是将以波束角和声波来回传播时间格式记录的多波束测量原始数据转换成：沿航迹、垂直航迹及深度格式的数据；剔除“飞点”（也就是异常值）是人工去除错误“跳点”的过程。每道工序也都需要一定的计算时间。特别是剔除异常值需要占据三分之二的时间，它需要一个条带一个条带的人工剔除，非常耗时。例如2km2面积的多波束扫测图，外业约一天，内业也占据了一天的时间，其中剔除异常值需要三分之二的时间，而单波束测图仅需要几分钟可出数据。

3.组合算法

3.1 构建条带数据序列

在本研究中，我们提出了一种空间模式下的异常值检测方案。在声速剖面校正和潮汐校正后，对MBS观测数据进行变换，得到在局部地理坐标系中定义的三维坐标。由于所获得的点云数据分布不均匀，不能直接作为后续处理的输入数据，为了满足算法的处理要求，将点云数据在相同距离中分成几条，这类似于将表面数据“切割”成多条数据序列（如图1所示）。条带与航行方向垂直。天然海底的形态几乎是连续的，当条带宽度较小时，单个条带中的地形点可以看成是一个窄的连续数据序列。然而，异常值会破坏系统的稳定性，这是利用条带数据序列检测异常值的理论基础。

图1 多波束测深数据划分示意图

3.2 组合算法处理

当稳定的数据序列只包含孤立的异常值时，小波分析就能准确地检测出异常值的位置。然而，当存在多个连续离群点时，小波检测结果与异常值的位置并不完全一致。换句话说，检测到的异常区域可能同时包含异常值和真实探测数据。此时，我们引入一个阈值来判断高频信号是否溢出。如果高频信号连续溢出，则可确定为异常区域。同时，对这一异常区域的范围进行了标记，滤波后得到了所有条带的新的高频信号。对新的高频信号和低频信号进行重构，得到新的条带数据.它既能消除随机噪声，又能减弱孤立点的影响，大大提高了卡尔曼滤波模型的精度。根据安装在船上的不同传感器的信息，确定卡尔曼滤波的参数。联合算法检测和删除异常值的基本过程。信号对新的高频信号和低频信号进行重构，得到新的条带数据。它既能消除随机噪声，又能减弱异常值的影响，大大提高了卡尔曼滤波模型的精度。根据安装在船上的不同传感器的信息，确定卡尔曼滤波的参数。组合算法检测和删除异常值的基本过程（如图2所示）：

图2 组合算法探测和移除异常值流程图

4.实测案例

利用仿真数据和实际数据对组合算法的性能进行了测试。利用仿真数据进行测试的目的是验证该算法的有效性。实际数据用于在实际环境中测试应用程序。测量误差通常有三种：系统误差、不规则噪声和异常值。由于这些数据是在早期处理的，因此将消除MBS的系统误差。不规则噪声虽然不可避免且不可预测，但可以用正态分布来表示。因此，这些数据的测试结果是可靠的。

4.1 模拟数据测试

本文利用SeaBat T20-P宽带MBS采集的点云数据对算法进行了测试。由于边缘数据质量差，只能用中间数据进行测试。该地区地形复杂，但地形起伏不大。水深在66m-71m之间，数据集包含10952个三维点。

测量环境中的异常值可能以两种形式出现：孤立点或连续区域。为了验证该方法的可靠性和稳定性，我们还将在这两种方法中添加异常值。模拟数据（如图3（b）所示）是通过添加异常值（总数的8%）生成的，范围从1至3m。这些条带沿航行方向以一定距离（测区的测深平均距离）构造。总共产生120条带。每个条带中的水深被重新排序，以产生新的数据序列。然后，创建新的索引，并与原始数据建立关联。该区域原始数据经小波分解得到的高频信号主要是测深过程中的随机误差，可用于计算原始均方误差σ0。截止误差为σ1=2σ0，由于条带数太大，我们选取了六条有代表性的条带的处理结果为例。这些条带共有43个异常值。利用DB3小波分析了每条带的高频信号和低频信号。当高频信号大于截断误差σ1时，数据被定义为超限。（如图4所示）除了112条带中1m的异常值外，所有异常值都已被检测到。

图4为各条带卡尔曼滤波结果。滤波值与正常深度基本一致，只是与峰值有明显偏差。当比较偏差和截断误差时，可以确定超限测深点为孤立点，并对其位置和指标进行标注。当所有条带的处理完成后，标记的离群点将从计算机内存中删除。然后对所有的条带进行重构，得到新的MBS数据集。

图3 第（a）20、（b）55、（c）57、（d）68、（e）93和（f）112条带数据的小波分解和异常区域的判断结果（红线代表异常区域）

图4 第（a）20、（b）55、（c）57、（d）68、（e）93和（f）112条带组合算法的滤波结果

4.2 实测数据案例

测试案例位于福建省湄洲湾大生岛附近海域。该海域水深约为5m-35m，水流流速较大，独立渔网较多，水流流速大容易产生噪声点，渔网多也会被多波束扫测到，这些并不是真实的海底水深数据，实际数据集包含约63000个水深，是在最大开启角为160°的SeaBat T50P上采集的。由于真实海底未知，本文从检测时间和探测位置的符合度两个方面对组合算法与人工编辑过程进行了比较。结果（如表1和图5所示）：

图5 实测数据异常值探测和移除的结果比较：（a）组合算法和（b）手工编辑

表1 组合算法和手工编辑在实际应用中比较

5.结论

该方法不需要迭代运算，适用于多核处理。根据模拟数据和实际数据的测试结果，该方法具有较强的鲁棒性和较高的去除效率。但是，由于测区边界点分布不均匀，难以构造连续的数据条带，无法对异常点进行完全检测和剔除。未发现的异常值主要集中在MBS海底的边界地区。此外，但该方法将一些少量的自然尖峰标记为异常值。该方法适用于天然水域，地势平缓，未施工区域，若有礁石等区域会被误删，因此，需要进一步优化各种参数，如小波阈值、卡尔曼递推方程的参数和条带宽度。