多波束测深系统是获取海底地形数据的主要技术手段^[1]。受复杂的海洋环境、测船噪声、人工操作不合理等诸多因素影响，多波束测深数据中存在着大量异常值^[2]。对异常值加以清理从而获取可靠的海底地形是多波束水深测量数据处理的主要工作之一。目前，多波束测深异常值一般采用人工交互剔除，工作量大，可信程度与作业人员的经验有关，具有一定的主观性^[3]。自动滤波算法数据处理效率高、结果不受主观因素影响，能够快速实时显示具有较高可信度的高精度三维海底曲面，一直以来都是国内外研究的热点问题^[3]。

多波束测深异常值自动滤波方法主要有4类：第一类是使用门限阈值或基于统计的滤波方法^[4-5]，设计适用于多波束测深数据的滤波方法；第二类是基于趋势面滤波的方法，此类方法重点在于如何确定拟合多项式的阶次和拟合海底区域大小^[6-8]；第三类是使用模拟物理变换进行异常值探测的方法，例如滚动圆^[9]、模拟布料^[10]等模型；第四类是运用抗差估计理论剔除多波束测深数据异常值的方法^[11-12]。上述方法尽管在一定程度上可以有效处理平坦海底的单个离散的多波束测深异常值，然而当海底地形复杂时或测深数据质量较差、存在连续异常值情况时，仍存在异常值无法清理完全、可靠性无法保证的问题，因此并未得到广泛应用。

目前，得到广泛认可的多波束异常值自动滤波算法是CUBE(combined uncertainty and bathymetry estimator)滤波算法，该算法引入了基于“4σ”准则的多水深假设跟踪模型^[13-16]。国内学者对基于“4σ”准则的多水深假设跟踪模型研究较多^[17-22]；为应对数据质量较差的情况下出现的连续的水深异常值，CUBE算法引入贝叶斯理论对多水深假设跟踪模型进行改进^[22-24]。目前，国内少见对基于贝叶斯理论的多水深假设跟踪模型的研究^[22]。文献[25]对基于“4σ”准则和基于贝叶斯理论的多水深假设跟踪模型进行了研究和复现，证明了基于贝叶斯理论的模型处理连续异常值时，较基于“4σ”准则的模型更稳健^[25]。然而该模型参数较多，参数配置复杂，调整不便^[26]，且数据质量较差的情况下出现较多异常值时，其处理效果并不理想^[26]。

CARIS软件将CUBE滤波算法作为其核心算法之一，不断加以完善和改进^[23-24]，改进后的CUBE滤波算法已无从探究其具体实现方法。依靠现有的文献资料复现的CUBE滤波算法，无法达到CARIS内置算法相同的效果。国内长期以来使用国外软件受到诸多限制，对国产化多波束测深数据处理软件需求越发迫切。国产化多波束测深数据处理软件因缺少可靠的异常值自动清理算法而遇到阻碍。随着机器学习和计算机视觉的发展，涌现出许多成熟的机器学习聚类算法，其中基于密度的DBSCAN聚类算法凭借无须先验信息、能够识别噪声点、对样本形态不敏感等优点而被广泛应用于异常值检测和异常值识别^[27-30]，这为提出一种流程明晰、参数简单、性能可靠的多波束测深值自动清理算法、助力国产化多波束测深数据处理软件发展提供了思路。

本文在深入研究CUBE算法的基础上，充分借鉴CUBE算法框架与流程，引入DBSCAN聚类算法对异常值进行识别与聚类，对CUBE滤波算法多水深假设跟踪模型加以改进，提出一种结合不确定度与密度聚类算法的多波束异常值自动滤波算法：首先，根据CUBE滤波算法构建网格节点可吸收水深点选取模型对测深数据网格化，利用DBSCAN聚类算法对节点可吸收的水深值进行一维聚类，根据水深点之间差值构建多个水深假设；然后，使用卡尔曼滤波更新节点的水深假设和不确定度；最后，从多个水深假设中选取可信度较高的水深假设作为节点的真实水深值。实测数据和仿真试验处理结果表明，本文算法滤波结果更稳健，可以有效处理多波束测深数据单个异常值和连续异常值，适用性更广泛。

1 DBSCAN聚类算法原理与参数设置

DBSCAN(density-based spatial clustering of applications with noise)聚类算法通过定义点邻域半径Eps和该邻域内包含最少的点数量MinPts来描述领域的样本分布紧密程度^[27-28]。该算法流程为：首先从聚类点集中随机选取任意一点，若该点同时满足算法所需的两个参数则该点为核心点。找到核心点后，则继续对该核心点邻域Eps内其他点依次进行判断，若某点邻域Eps内依然满足Minpts参数，则其被标记为核心点并加入类中，否则判定为边界点加入类中，直到没有符合条件的点为止。如果数据中的某个点既不满足核心点的要求，也不满足边界点的要求，则该点被归类为噪声点。算法通过不断重复上述过程，遍历数据中所有的点完成对多波束测深数据点的聚类。DBSCAN算法具有无须预先指定分类数量、可以识别噪声点、能够发现特殊形状的空间类别等优点^[28]，因此无须多波束测深数据点空间形态等先验信息，对于数据具有良好的适用性。其本质上是根据点之间的距离和邻域点数量来区分正常值与异常值，符合人工判读测深数据异常值的基本原则，使其能够适应数据质量较差的情况下测深数据异常值的自动清理工作。

合理地设置DBSCAN聚类算法的两个关键参数Eps和MinPts，可以实现多波束测深数据正常值和异常值的聚类，从而构建起多个水深假设，为剔除异常值的做准备。参数Eps描述了数据点的邻域距离阈值，用于计算核心点为圆心Eps为半径的圆区域内水深点的数量。从理论分析，Eps取值范围为(0, +∞)，从大量算例来看，参数Eps取0.8~1较为理想。考虑多波束测深数据点深度之间的差异较小，算法将参数Eps取值范围设置在(0, 1]，实际使用时取值为1。该参数设置偏小会导致聚类数量增加，类内包含的水深点数量会减少。参数MinPts用于描述核心点在Eps为半径的邻域范围内最少数据点个数的阈值，本算法中参数MinPts固定取值为1，确保吸收水深点数量较少的节点能够正常聚类而不被识别为噪声点。当数据质量较差，节点吸收的异常值过多时，MinPts的固定为1的设置可以有效地区分不同水深区间的异常值，从而确保正常值的类内包含的水深点数量多于异常值类内包含的水深点数量。

本文首先使用CUBE滤波算法构建网格节点可吸收水深点选取模型，构建推估曲面对测深数据网格化，确定网格节点可以吸收的水深点；然后使用DBSCAN聚类算法对一维水深值进行聚类。本文算法参数MinPts固定设置为1，调节Eps成为影响算法效果的重要因素。该过程将处理水深点空间三维信息转化为处理水深点一维水深信息。由于处理信息维度的压缩，此时Eps从空间欧氏距离转化为曼哈顿距离，而正常值之间水深差异较小，所以参数Eps可调整的空间也随之压缩。因此本文算法具有参数设置简单、具有较强的适用性的优点。

2 结合不确定度与密度聚类算法的多波束异常值自动滤波算法流程

结合不确定度与密度聚类算法的多波束异常值自动滤波算法主要包括推估曲面构建、水深值聚类、质心点计算、水深值和不确定度推估、选取最优估值、原始测深数据滤波6个步骤，详细流程如图 1所示。借鉴CUBE滤波算法网格节点可吸收水深点选取模型，实现网格节点与水深点关联；测深点的不确定度为衡量水深点的可信程度提供了标准，DBSCAN聚类算法提供了区分多波束测深数据正常值和异常值的工具。水深点经过卡尔曼滤波器的迭代完成水深信息的传递和不确定的收敛，得到该节点的水深假设；设置水深假设选取策略，得到推估曲面上节点的最优水深估值；最终将推估曲面覆盖到原始数据上完成多波束测深数据的自动滤波。下面详细讲解本文算法的6个步骤，选取聚类情况较为复杂的某节点作为实例数据用以必要时辅助说明。

2.1 推估曲面构建

推估曲面是独立于原始测深数据的一网格化曲面，曲面上的节点由网格化确定平面坐标，节点的水深值由其邻域的水深点推估而来。为了实现海底地形高精度高分辨率的精细化处理，海底地形网格间距设置为1 m，需先计算节点邻域可供吸收的水深点。在计算节点吸收域的同时，也要计算水深点的影响范围，文献[23]给出两种范围计算公式，节点吸收域为

式中，R_cap是节点的吸收域；μ_capdiaca为节点捕捉距离比例系数，设置为5；μ_capdiamin为最小捕捉半径，设置为0.5；G_depth是以节点为中心3倍网格为半径的圆内水深平均值。水深点的影响范围为

式中，R是水深点的影响半径；D_Grid是网格间距；TVU_max是IHO S-44(2020版)标准有关规定中不同等级海道测量所允许的极限垂直不确定度，计算公式具体见文献[31]，本文不再详细描述，只给出使用的参数a、b的值分别为1和0.023；TVU是水深点的垂直不确定度；μ_dist是不确定度随距离增大的比例系数，设置为2。需要说明的是，本文算法借鉴CUBE滤波算法的网格节点可吸收水深点选取模型，以上参数设置遵循CUBE滤波算法的默认设置^[23]，进行试验时，参数保持一致，以便对比相同参数下的试验效果。

只有水深点位于节点吸收域内同时节点也位于水深点的影响范围内，该水深点才可以被节点吸收，如图 2所示。这是CUBE滤波算法极具特色的“双向选取”机制，形成推估曲面完成了水深点三维信息(X, Y, Z)中平面信息(X, Y)与深度信息(Z)的分离，因此DBSCAN算法对一维水深信息进行聚类即可。

2.2 水深值聚类

该节点选取完邻域可供吸收的水深点之后，节点和水深点已经完成二维平面上的相关联。之后使用DBSCAN聚类算法处理一维水深数据，便于两个关键参数的确定。DBSCAN聚类算法良好的聚类结果是进行水深值和不确定度推估更新、最优值的选取的基础。节点通常会获取多个水深值，多波束测深数据正常值经过正确的归位计算，周围领域的测深点数较多；而异常值呈现误差的偶然性，可能散落在任意的水深区间内，且周围领域的测深点数较少。对于单个异常值，其周围没有测深点，虽然DBSCAN聚类算法将其识别噪声，但MinPts设置为1确保其也会成为一个水深假设，进入最优估值选取阶段因不确定度较大而被淘汰；对于连续异常值，其周围存在少数水深值接近的异常值测深点，被聚成一类，形成水深假设进入到最优估值选取阶段，与其他类进行比较选取最优估值。

一般情况下，聚类结果分为两类，且两类内包含的水深点数量差距较为悬殊，很容易识别正常值与异常值；当遇到较为复杂情况时，聚类数量会增加。表 1显示了较为复杂情况下的某节点水深值聚类情况。节点32个可吸收的水深值被聚成6类，其中90 m以深水深区间聚成1类有19个，50~60 m水深区间聚成2类共6个，60~70 m水深区间聚成2类共4个，40~50 m水深区间聚成1类有3个。该水深节点吸收域为3.093 m，根据多波束测深数据滤波基本原则中的水深变化区间原则，即海底地形变化较为平缓，一定范围内的水深应该处于一定的区间内，连续自然变化的海底在吸收域范围内不会出现剧烈的水深变化，结合试验数据分析，认为40~70 m水深区间的水深值是应当加以清理的异常值。DBSCAN聚类算法准确识别了水深值之间的差异并进行了聚类，聚类效果良好，达到了试验预期。

2.3 质心点计算

聚成的类中含有多个水深点数据，需要逐个进入卡尔曼滤波器进行水深值和不确定度的推估。第一个进入滤波器的水深值作为初始水深值，应能代表该类水深值的所在水深区间。在聚类算法中，通常使用类的平均值作为质心，代表该类的特点。然而质心不一定是类中的水深点，因而无法获取质心的不确定度。本文算法使用类质心的近似值作为第一个进入卡尔曼滤波器的水深值。该类水深值与平均值差异最小的水深值作为质心的近似值(表 1)，解决了质心的不确定度无法确定的问题，同时第一个进入卡尔曼滤波器的水深值又能较好地反映该类水深的平均深度。

设某一节点可吸收水深值聚类结果为M(i)，M(i)={d₁ⁱ, d₂ⁱ, …, d_jⁱ; tpu₁ⁱ, tpu₂ⁱ, …, tpu_jⁱ}，d_jⁱ表示第i个类中第j个水深值，tpu_jⁱ第i个类中第j个水深点的不确定度，则通过式(3)确定第i个类中第k个水深点为类质心近似值

类质心的近似值求取流程如图 3所示。

2.4 水深值和不确定度推估

启用卡尔曼滤波器进行水深点和不确定度的推估，推估形成的可能代表该节点水深值，称为水深假设。推估流程及卡尔曼滤波使用公式见文献[13, 16—18]。卡尔曼滤波器随着吸收水深点数量的增加，水深假设的不确定度逐渐减小。这一过程是在先验观测值(类质心的近似值)的基础上，多个后验观测值不断对先验观测值进行修正。图 4显示了某节点水深值和不确定随着水深点个数增加推估得到的水深值和不确定度，其中不确定度以误差棒表现，误差棒缩短代表的水深值不确定度越来越小，水深值的可信度越来越高，这一规律成为最优估值选取的基础。

2.5 最优估值选取

DBSCAN聚类算法的类数量会影响其中的水深点数量，水深假设吸收水深点的数量会最终影响其最后的不确定度大小，因此最优估值选取建立在DBSCAN聚类算法良好的聚类效果和卡尔曼滤波推估基础上，选取具有最小的不确定度的水深假设作为节点的最优估值。这种选取方式既考虑到DBSCAN聚类算法的类数量影响，同时考虑每一类中的水深点数量的影响和不确定度的影响，是三者结合下最优估值的一种选取方式，其选取公式如下

2.6 原始测深数据滤波

推估曲面构建的过程伴随着网格化过程，最优水深值的推估伴随着水深值的传递，会造成水深值在平面的位置偏移，但不影响其作为粗差剔除的参考曲面^[32]。因此需要使用推估曲面覆盖到原始测深数据上，进行异常值剔除。最优估值选取完毕，2.1节中构建的推估曲面上每个节点获取了水深值。这些水深值的不确定度很小，可以认为是现有测量条件下获取的最接近真实观测值的水深值。本文使用基于近似真值的异常值判断准则进行异常值剔除

式中，为水深推估值，本文中即为水深节点推估的最优水深值；Z_p为节点可吸收的水深点；N为该节点范围内可吸收的水深点数量；k一般取2或者3。根据误差理论，标准差用于衡量数据之间离散程度，而均方根误差衡量的是数据与真值之间的离散程度。考虑到近似真值的水深值已经获取，因此此处取为均方根差，用于衡量原始测深数据点与推估曲面最佳水深值之间的离散程度。

根据式(5)、式(6)对原始测深数据进行异常值滤波。值得注意的是，由于一个水深点可能被多个节点吸收，单个节点范围内水深点被识别为正常值并不能保证其不被剔除，需要多个节点同时识别为正常值，该水深点才能被保留。

3 试验结果与分析

为验证算法的有效性，本文通过Python语言编程实现了结合不确定度与密度聚类算法的多波束异常值自动滤波算法，并根据文献[14—16, 23]编程实现CUBE滤波算法。分别采用手工剔除法、CUBE滤波算法和结合不确定度与密度聚类算法的多波束异常值自动滤波算法对包含有异常值的实测数据和仿真数据进行滤波试验。最后使用Origin2021软件对试验结果进行了可视化显示和对比分析。

3.1 试验数据分析与预处理

试验数据为2020年浙江舟山金塘海区海底地形多波束测量数据，采用iBeam8120国产多波束测深系统，该系统工作频率400 kHz，每Ping有256个波束，实际作业开角140°，测深分辨率为6 mm。数据在采集过程中已经进行了安装偏差改正，并在后期进行了声速、潮汐及姿态改正，试验数据已排除系统误差影响。测区水深在80 m左右，水质较为浑浊，能见度较差。

数据预处理分为两步。第一步根据文献[33]，使用CARIS软件中HGM模型计算测深点的水平不确定度和垂直不确定度，并根据文献[13, 23]中计算出水深点的不确定度。第二步剔除边缘波束。图 5为剔除海底边缘波束的前后对比，图 5(a)为未剔除边缘波束的海底曲面，图 5(b)为剔除边缘波束之后的海底曲面，对比可以看到边缘波束含有大量集群的异常值，导致真实海底信息丢失，不利于节点推估真实海底情况。本文根据多波束水深测量作业经验，选取3倍实测水深值的波束开角作为试验数据。多波束Ping开角为140°，覆盖宽度约为305 m，平均水深为80 m左右，按照比例计算，即选取中央约200个测深点作为有效数据进行验证。原始数据中测点共50 944个，剔除边缘波束后共有测深点39 599个。

3.2 实测数据试验分析

使用Origin2021软件根据网格数据构建三维海底曲面，试验结果如图 6所示。对原始测深数据、手工剔除法、CUBE滤波算法、本文算法处理的多波束数据的水深均值、水深极大值、水深极小值、水深标准差、测深点数量进行了统计分析，结果见表 2。

由实测数据试验结果对比分析可知：①手工剔除法、CUBE滤波算法和本文算法的结果水深值范围均小于原始实测数据，说明了3种方法均可以滤掉原始实测数据中的异常值点；由图 6可知，经过处理后测深数据的异常值均减少(或完全剔除); ②手工剔除法和本文算法结果水深极小值均为88.62 m，极大值均在94~95 m，见图 6(b)、(d)，比较真实地还原了连续的自然海底变化，而CUBE滤波算法结果极小值为47.88 m，极大值为110.64 m，说明存在一些异常值未能剔除，见图 6(c)；③手工剔除法和本文算法滤波效果较为理想，CUBE滤波算法滤波处理后的测深数据仍存在着个别异常值点，经与原始实测数据对比分析，发现CUBE滤波算法处理结果残留的异常值点主要出现在原始测深数据中连续异常值出现的区域。

3.3 仿真数据试验分析

为进一步说明本文算法的普适性，下面进行仿真试验，对比分析CUBE滤波算法和本文算法处理单个异常值和连续异常值的能力。设计两个试验对比方案，在上述海区人为设置单个异常值和连续异常值两种形态的异常值进行仿真滤波试验。从手工剔除异常值的多波束数据中选取一个近似正方形的区域作为仿真算例的数据，选取区域内包含15 Ping数据，每Ping包含18个测深点，共270个测深点，该区域边长约为15 m，面积约为225 m²。

方案一：随机生成10个均值为80 m、标准差为3 m的高斯噪声点，随机加入该区域内，得到其三维海底曲面图和局部放大图如图 7(a)所示，图 7(c)为添加单个噪声局部放大区域，添加的单个噪声点明显凸出于周围的地形趋势，可认为是该区域的异常值，应该加以剔除。

方案二：随机生成均值为80 m的高斯噪声点20个，添加到该区域内连续的两Ping数据内，每Ping数据各添加10个噪声点，得到其三维海底曲面图如图 7(b)所示，图 7(d)为添加噪声区域局部放大图，添加入的连续噪声点明显凸出于周围的地形趋势，可认为是该区域的异常值，应该加以剔除。使用CUBE滤波算法和本文算法分别对方案一、方案二进行滤波处理，结果见图 7(e)—(h)。统计该区域添加单个噪声点、添加连续噪声点、CUBE单个噪声点滤波、本文算法单个噪声点滤波、CUBE连续噪声点滤波、本文算法连续噪声点滤波处理后的水深均值、水深极大值、水深极小值、水深标准差、测深点数量，统计分析见表 3。

由仿真试验结果分析可知：①方案一中，CUBE滤波算法和本文算法与添加单个噪声的数据对比，极小值均为88.62 m，与表 2手工剔除法的极小值数据一致，说明处理单个异常值时，CUBE滤波算法和本文算法据能够识别并加以剔除；图 7(e)、(f)与图 7(a)相比，凸出周围地形趋势的异常值被剔除，得到了连续变化的真实自然海底，与图 6(b)手工剔除法一致。②方案二中，表 3中CUBE滤波算法处理连续噪声点结果极小值为76.14 m，而本文算法处理结果极小值为88.62 m，与表 2手工剔除法的极小值数据一致，证明CUBE滤波算法未能将添加的连续异常值剔除完全，本文算法可以较好地剔除连续异常值；图 7(g)显示了CUBE滤波算法无法完全剔除异常值，推估曲面构建受到异常值的影响，从而影响对测深数据异常值的剔除，导致海底地形不真实；图 7(h)显示本文算法剔除了连续异常值，得到了同图 6(b)相近的连续变化的真实自然海底。

通过上述实测数据与仿真数据对比分析，可以得出以下结论：处理单个异常值时，CUBE滤波算法与本文算法均可以识别异常值并加以剔除；处理连续异常值时，CUBE滤波算法无法有效剔除异常值，而本文算法可以对连续异常值有效清理，显示出了清理不同形态多波束测深数据异常值的良好特性。

4 结语

本文借鉴CUBE滤波算法网格节点可吸收水深点选取模型，使用DBSCAN聚类算法处理一维水深数据，通过卡尔曼滤波推估最优估值并给出了最优估值的选取策略，提出了一种结合不确定度与密度聚类算法的多波束异常值自动滤波算法，并通过实测数据和仿真算例验证：本文算法和CUBE滤波算法清理单个异常值时滤波效果较为理想，处理具有连续异常值的多波束测深数据时，CUBE滤波算法无法将异常值完全清理，本文算法能够准确识别连续异常值并对异常值加以剔除，显示了本文算法对不同形态的异常值具有良好的适用性。本文算法流程明晰、参数简单、性能可靠，具有实际的工程应用价值。然而为适应海量多波束测深数据处理，本文算法在程序编写与运行效率方面仍有待提高。