2. 武汉大学海洋研究院, 湖北 武汉 430079;
3. 武汉大电气与自动化学院, 湖北 武汉 430072
2. Institude of Marine Science and Technology, Wuhan University, Wuhan 430079, China;
3. School of Electrical Engineering and Automation, Wuhan Universality, Wuhan 430072, China
多波束测深系统(multibeam bathymetric system, MBS)因具有高密度、全覆盖测深特点[1-2],在港口、航道等工程中有着广泛的应用。实际应用中,特别在海岸防波堤抛石工程中需要边施工边测量,对抛石间较大缝隙进行填补,这就需要比较准确和清晰的水下地形信息。然而,施工区域水下地形测量受悬浮物、石块遮挡、系统噪声和测量船噪声等影响,MBS获取的测深点云存在大量粗差,严重影响了测深数据对海底测量对象的准确描述,不能给予抛石施工准确的指导信息,必须给予有效滤除。据此,国内外学者开展了大量研究。目前,常见滤波方法主要有COP法[1]、Ware法[1]、Knight & Wells法[1]、Eag(RDANH)法[3]、趋势面法[1, 4-5]、抗差估计法[6-7]、Bayes估计[8]、中值滤波[9]、均值滤波、局部方差检测和小波分析相结合的滤波方法[10]以及选权迭代加权平均[11]等方法。上述方法对于海量MBS点云数据处理存在速度慢,无法满足边施工边测量,近实时获取水下地形的需求,且普适性不高,乱石区等复杂海床测深数据滤波性能欠佳等不足,导致无法准确获取水下地形信息,无法正确指导抛石施工。在文献[12]提出MBS测深点的水平和垂直不确定度理论计算方法的基础上,文献[13—15]提出了CUBE滤波方法。由于该方法具有滤波高效、可靠、抗差、稳定等特点[16-19],目前在CARIS和HYPACK等国际商用MBS测深数据处理软件中广泛使用[20]。当测深点在格网节点周围均匀分布,存在少量粗差时,CUBE算法可准确估计格网节点水深值。然而,在乱石区,现有CUBE算法会遇到如下问题:其一,MBS测量会因乱石遮挡产生测量盲区,导致盲区测深数据空白,波束可达区密集,测深数据分布(特别是在边坡乱石区)不均匀性问题突出,导致现有CUBE算法对边坡乱石区测深数据处理时估计节点的深度不准确;其二,现有CUBE算法仅开展深度估计,未顾及平面位置估计,乱石区的上述数据特征会引起估计的水深点位置不准确,从而引起地形特征混乱;其三,边坡乱石区地形变化复杂,波束在乱石间会产生多次回波后到达换能器,或因环境复杂,水中存在大量遮挡物时导致测深数据中存在大量粗差,基于现有CUBE算法的最优水深估值选择原则[13]将很难得到准确的水深,还需要借助手工编辑剔除粗差,也因此显著降低了数据滤波的效率,严重影响施工的进度。
为此,针对CUBE算法不足,在深入研究的基础上,提出了顾及平面位置估计和开展二次CUBE估计的滤波算法,以期解决上述问题,实现MBS测深数据对乱石区地形的真实反映。
1 CUBE算法原理及其不足 1.1 CUBE算法CUBE是基于规则格网估计各格网节点最优水深的滤波方法,主要包括如下3个部分。
1.1.1 测深点深度不确定度传递在根据水深确定格网节点搜索范围内,将测点水平和垂直不确定度σH, j、σV, j传递到待估格网点上[13]
式中,δij为测深点到待估格网点的平面距离;SH为水平不确定度的影响因子,默认为1.96;α为地形起伏影响因子,默认为2;Δmin为格网大小。
1.1.2 多重估计以Kalman滤波方法建立估计方程[13],式(2)和式(3)分别给出了状态方程与观测方程
式中,w[n]、v[n]为系统噪声和观测噪声,均服从正态分布。由于模型是为了获取格网节点的唯一水深,因此有W[n]=0 m2(文献[13])。
式(4)—式(9)给出了滤波估计过程[13]
式中,
定义图 1(b)中σf为预测误差,图 1(c)中d1和d2分别为输入水深值与估值1、估值2差值,则CUBE算法多重估计更新过程为:
(1) 初始化一个估值后或者只存在一个估值时,进行如图 1(a)或图 1(b)所示的过程。
(2) 存在多个估值时(如图 1(c))选择与输入水深值差异最小估值,重复图 1(a)或1(b)。
(3) 重复步骤(1)和(2),对所有数据格网节点进行水深估计。
1.1.3 最优水深估值选择原则[13]原则1:选取估值纳入点数最多的一个估值作为最优估值。
原则2:在待估格网点周围一定范围内,寻找单个估计值的格网点水深为参考水深,寻找待估格网节点估值中与该参考水深最接近的估值作为最优估值。
原则3:综合原则1和2。
1.2 边坡乱石区CUBE算法的不足(1) 忽略了边坡乱石区测深点的分布特点。
边坡乱石区MBS测量时,因大坡度地形、石块(图 2(c))遮挡,导致大量测量盲区(图 2(a)),测深点在垂直方向上间断分布(图 2(a)、(c))。现有CUBE算法基于测点均匀分布(图 2(b))和多重估计原理开展估值图 2(c)中的测点分布不均匀问题会导致对节点Node的估计产生两个水深估值A和B。无论是选择A估值还是B估值,最终会赋值给Node作为该节点的估值,进而造成节点平面位置和水深估值不准,引起地形特征错位甚至错误。
(2) 无法消除大量粗差影响。
传统CUBE滤波给出了3个水深估值优选原则,其优缺点如下:
原则1:简单,易于实施,当地形平缓和测深数据中存在少量粗差时,可有效剔除粗差,但当粗差较多且集中分布时,易获得错误水深估值。
原则2:当地形平缓、粗差较少或者离散分布时,可有效消除粗差,具有抗差性,但当地形复杂、数据质量较差时,多重估计后将有大量的格网节点出现多个估值情况,此时很难找到只有单个估值情况。如取舍不当,会导致待估格网点水深错误。
原则3:情况同原则1和2。
实际数据处理时,如图 3(a)中红色区域为边坡乱石区,采用CUBE算法对该区水深数据进行滤波,每个格网节点将获取多个估值,因此原则2、3无法实施,只能采用原则1选择最终估值。然而,针对图 3(b)矩形区测深数据中存在大量粗差,借助原则1获取的估值如图 3(c),可看出传统CUBE并未完全消除粗差的影响。
2 顾及数据特征的二次CUBE滤波算法
在边坡乱石区,为解决因原始波束数据中存在大量粗差产生的由传统CUBE滤波带来的节点估值不准问题(如图 3(c)所示),提出一种快速二次CUBE滤波方法。
2.1 顾及位置估计的一次CUBE滤波针对1.2节中问题(1),本文提出在传统CUBE算法的基础上,在进行水深估计的同时对水深估值的平面位置估计。平面位置估计与水深估计类似,将水深值z用平面坐标(x、y)代替,用每个波束的水平不确定度代替水深不确定度。采用类似式(5)、式(7)和式(8)的形式,则有
在传统CUBE滤波方法的基础上,增加位置估计,对原始的测深数据开展一次滤波。顾及平面估计的一次CUBE滤波所得散点(图 4中所示的三角网节点)不再是每个格网节点。
2.2 二次CUBE滤波
提取一次CUBE滤波所得散点水深值、不确定度及对应的每个点上多重估计水深。如图 4所示,图中中心点(xm, yn)为待检测点,以该点为中心,直径为a的圆范围内测深点对该点估值进行检测,a=2x×Δmin(x≥1, x∈N)。节点上估值的垂直不确定度为
式中,k表示节点(xi, yj)的第k个估计值;nk表示该估值的吸收的测深点数;σv, j表示该估值纳入的第j个测深点的垂直不确定度。
将周围点(xi, yj)的不确定度σ(i, j)传递到待估点(xm, yn)上。本文顾及平面位置影响和乱石区地形梯度变化显著特点,给出传递模型如下
式中,μ为地形坡度;σi, j2为不确定度;Δmin为第1次CUBE滤波时设置的格网大小。
当水深小于50 m,波束大小为0.5°×1°,仪器量程分辨率为1.25 cm,条带覆盖角120°,假定节点不确定度σ(i, j)为0.15 m,则根据式(14),表 1列出了当周围点距待估点距离为1 Δmin时不同μ将对应不同大小的σi, j值。由表 1可知,假设相邻节点水深值差异为0.5 m,为容许该水深差,σi, j应该大于0.5 m,此时μ应大于3,取值4。μ不能过大,否则可能会导致粗差滤除失败。μ值参数需根据当前测量水深的不确定度和地形起伏来综合确定。
然后,采用多重估计对水深点值与周围一定范围内的水深值进行一致性检测(如图 1)。若发现该水深值带有粗差,利用周围水深值采用多重估计和最优参考水深选择方法获取最优水深参考值dr。dr的选取过程如下:
(1) 将一次CUBE滤波获得的待估点水深作为初始估值d0,利用直径为a区域内的所有点的平面和水深一次估值结果,对待估点进行多重估计,并统计每个估值纳入的点数ξi。
(2) 计算每个估值吸纳区域内的点数占所有有效点数的百分比
(3) 根据各估值的百分比εi,寻找最大的εmax。若εmax≥ε',则εmax对应的估值即最优水深参考值dr;否则,表明未能找到任何可靠的水深参考值,删除或标记待检测节点为不可靠点。
最后,给出如下方法,开展二次水深估计:
(1) 从一次CUBE滤波多重估计结果中获取与dr差异最小的水深估值d(m, n, index)
式中,index为对应索引数; count为待估点上的一次CUBE滤波获取的多重估值个数。
(2) 然后计算地形梯度变化量ddiff
(3) 顾及地形梯度ddiff,基于式(17)的原则确定最终水深估值df
式中,df为最终水深优选估值。
获取可靠水深值的同时,更新该点平面位置。对所有一次CUBE滤波后点进行上述二次滤波处理,最终完成乱石区测深数据的滤波。
3 阈值参数的确定 3.1 搜索范围直径a的确定如图 4所示,极端情况下,假设黑色点及中心点全为粗差点,随着直径a增大,区域粗差率降低。理论上,粗差率越低,越易获得无粗差影响的参考值,但对于复杂地形,a过大,最终获得的参考值与待估点的地形相关性将大大降低,使得待估点水深估值优选结果不可靠。所以参数a取值时既要考虑粗差率,又要顾及地形起伏。
当直径a取2、4、6倍Δmin,图 4中圆形区包含的粗差格网点数占该区域总点数(包含区域边缘点)的百分比分别为100%、42.9%和20.5%。粗差率一定的情况下,地形变化较大,a应取较小值,这样获取参考值时纳入的节点与待估点越靠近,参考值与待估点水深值在地形上相关性更大,最终优选得到的水深值越可靠。粗差率通常是未知的,但与地形坡度、水深和测量设备等相关。试验区大量试验表明,水深小于50 m,地形坡度大于30°,Δmin=0.5 m,粗差大量聚集时,a取6 Δmin较为合适;地形坡度小于30°,a取6~8倍Δmin较为合适。
3.2 有效点数阈值η给定η是确保以待估点为中心,直径为a范围内有足够多的有效水深值点数。如果有效点数太少,获取的参考水深可靠性将大大降低。地形边缘或数据稀疏地方,有效点数较少,这种点的可靠性较差(往往表现为破碎地形),可以将其删除或者标记。有效点数取值可以根据直径为a的搜索区域的外接正方形包含的点数得到,即η=ceil[(a/Δmin+1)2/2],其中ceil函数为向上取整。
3.3 阈值ε'的确定阈值ε'用于判断待检测节点水深值和选择可靠的参考水深估值。如图 4所示,若图中较大的黑色格网圆节点及中心点皆为异常点,由3.1节可得,当a=6 Δmin时,粗差率为20.5%。如果ε'小于粗差率,可能将待估点水深值判断为可靠值,从而导致错误检测,因此,ε'必须大于粗差率,才能保证获得正确的检测结果。同时,为了获得可靠的水深参考值,阈值ε'不宜过大。因此,ε'应满足20.5%<ε'<50%。
4 试验及分析为验证本文方法对乱石区MBS测深点滤波的有效性,试验采用以色列海岸防波堤抛石工程施工区测量数据。该数据采用Sonic2024多波束测量仪器[21]采集,表 2为Sonic2024各项误差分量[13]。测区水深1~25 m,最大浪高小于0.25 m,海洋潮汐影响0.1 m左右。边坡抛石区石块大小不一,呈现不规则分布。根据施工要求,每次抛石后都要对水下抛石区域进行测量,准确掌握每次抛石情况,为下一次抛石提供指导。受乱石区多次反射波等影响,测量数据中存在大量粗差,如图 5中椭圆形虚线区域所示。
误差项 | 数值 |
设备位置偏移误差 | 5 mm |
波束开角误差 | 0.05° |
GPS-RTK定位误差 | <10 cm |
GPS时延误差 | 0 ms |
Roll/pitch测量误差 | 0.01° |
Yaw测量误差 | 0.01° |
IMU时延误差 | 0 ms |
声速剖面测量误差 | 0.5 m/s |
表层声速测量误差 | 0.25 m/s |
Heave固定测量误差 | 0.05 m |
Heave测量误差比例因子 | 5% |
静吃水测量误差 | 0.02 m |
动吃水(squat)测量误差 | 0.02 m |
吃水量变化量(loading)测量误差 | 0.01 m |
船速SOG误差 | 0.01 m/s |
潮位测量误差 | 0.02 m |
潮位空间变化误差 | 0.02 m |
以下采用3种数据处理方案开展试验,分析本文方法的性能。
方案1:采用Caris中传统CUBE算法滤波。基于CARIS软件(内置CUBE滤波)对测深数据处理,结果如图 6(a)和图 7(a)所示。方案1基于格网数据估值,估值结果仍存在大量粗差,错误地反映了边坡乱石区地形。
方案2:二次CUBE算法滤波。试验中Δmin=0.5 m、a=3 m、η=25、ε'=25%、μ=4。试验结果如图 6(b)、(d)和图 7(b)所示。因同时开展了水深和位置估计,因此估值结果不再是格网节点,而是一系列散点,采用三角网构建地形模型更为恰当(图 6(d))。利用该方案检测出了大量的异常测深点如图 7(d)所示(图 7(a)、(b)、(c)均为(d)中红色矩形区域侧视图),由图 7(a)、(d)可知,检测出的异常点分布与实际吻合,最终滤波结果粗差剔除彻底。抽取图 5中边坡乱石区矩形区域内(图 5(b))的原始测深数据,采用方案1和方案2处理,滤波结果分别如图 6(a)、(c)和(b)、(d)所示。可以看出:
(1) 方案2块石地形特征更加明显,石头间隙更加分明,粗差被彻底消除。
(2) 比较图 6(c)、(d),方案2得到的估计点在平坦地区几乎与方案1得到的格网节点分布一致,在复杂地形区则存在明显差异,但对乱石区地形描述更真实。
(3) 方案2顾及了乱石区测深数据特点,同时开展了位置和水深估计,并基于二次滤波消除了粗差影响,因此为测深点赋予了准确的三维估值,从而也实现了对地形的真实表达。
方案3:CARIS手工编辑。方案3基于CARIS滤波后出现粗差时常采用的一种补救措施,因边坡乱石区地形复杂,有时难以判断准确的地形信息,需要有丰富的水深数据处理经验。如图 7(c)所示,异常测深点被剔除,基本反映了实际地形变化,但过程非常耗时。以方案3结果为参考,方案1和方案2滤波获得的水深值与之比较,差值分布如图 8所示,统计结果如表 3所示。
(1) 方案2偏差最小,且偏差分布符合正态分布。
(2) 方案1偏差均值较大,主要因其机理不足所致,也进一步表明了本文工作的必要性。
(3) 方案1中的问题未在方案2结果中出现(如图 7(b)和图 8(b)、(c));方案2实现了自动处理,且将方案1的成果精度提高了近5倍,从而也验证了本文方法的有效性。
5 结论本文基于数据特征提出的二次CUBE滤波算法,给出了位置和水深的同步估计方法、顾及地形梯度的平面和水深不确定度传递模型、参考水深的多重估计和优选算法并顾及地形梯度的二次水深估计模型,解决了传统CUBE滤波算法在乱石区出现的水深估计不准确、地形特征模糊和粗差剔除能力不足等问题,实现了边坡乱石区水深的准确估计。试验表明,二次CUBE滤波算法实现了粗差的自动检测和测深点的准确估值;在乱石区,大大提高了传统CUBE算法水深估值的精度,形成的地形特征更加准确和明晰;在平坦海床,取得了与传统CUBE滤波算法近似相同的地形结果。相对传统CUBE滤波算法,二次CUBE算法的抗差性更强,数据处理效率更高。
[1] |
赵建虎.多波束深度及图像数据处理方法研究[D].武汉: 武汉大学, 2002. ZHAO Jianhu. Research on Multibeam Depth and Image Data Processing[D]. Wuhan: Wuhan University, 2002. |
[2] |
赵建虎.
现代海洋测绘[M]. 武汉: 武汉大学出版社, 2008.
ZHAO Jianhu. Modern Marine Surveying and Charting[M]. Wuhan: Wuhan University Press, 2008. |
[3] | WILSON M F J, CONNELL B O, BROWN C, et al. Multiscale Terrain Analysis of Multibeam Bathymetry Data for Habitat Mapping on the Continental Slope[J]. Marine Geodesy, 2007, 30(1-2): 3–35. DOI:10.1080/01490410701295962 |
[4] |
董江, 任立生.
基于趋势面的多波束测深数据滤波方法[J]. 海洋测绘, 2007, 27(6): 25–28.
DONG Jiang, REN Lisheng. Filter of MBS Sounding Data Based on Trend Surface[J]. Hydrographic Surveying and Charting, 2007, 27(6): 25–28. DOI:10.3969/j.issn.1671-3044.2007.06.007 |
[5] |
王海栋, 柴洪洲, 翟天增, 等.
多波束测深异常的两种趋势面检测算法比较[J]. 海洋通报, 2010, 29(2): 182–186.
WANG Haidong, CHAI Hongzhou, ZHAI Tianzeng, et al. Comparison of Two Trend Surface Detection Algorithms of Multibeam Bathymetry Outlier[J]. Marine Science Bulletin, 2010, 29(2): 182–186. DOI:10.3969/j.issn.1001-6392.2010.02.011 |
[6] |
张志伟, 暴景阳, 肖付民.
抗差估计的多波束测深数据内插方法[J]. 测绘科学, 2016, 41(10): 14–18.
ZHANG Zhiwei, BAO Jingyang, XIAO Fumin. An Interpolation Method of Multibeam Data Based on Robust Estimation[J]. Science of Surveying and Mapping, 2016, 41(10): 14–18. DOI:10.16251/j.cnki.1009-2307.2016.10.004 |
[7] |
赵东明, 柴洪洲, 王海栋.
基于抗差最小二乘配置的海底地形生成研究[J]. 系统仿真学报, 2010, 22(9): 2091–2094, 2099.
ZHAO Dongming, CHAI Hongzhou, WANG Haidong. Research on Seabed Terrain Generation Based on Robust Least-squares Collocation[J]. Journal of System Simulation, 2010, 22(9): 2091–2094, 2099. |
[8] |
黄贤源, 隋立芬, 翟国君, 等.
利用Bayes估计进行多波束测深异常数据探测[J]. 武汉大学学报(信息科学版), 2010, 35(2): 168–171.
HUANG Xianyuan, SUI Lifen, ZHAI Guojun, et al. Outliers Detection of Multi-beam Data Based on Bayes Estimation[J]. Geomatics and Information Science of Wuhan University, 2010, 35(2): 168–171. |
[9] | MANN M, AGATHOKLIS P, ANTONIOU A. Automatic Outlier Detection in Multibeam Data Using Median Filtering[C]//2001 IEEE Pacific Rim Conference on Communications, Computers and Signal Processing, Victoria, BC, Canada: IEEE, 2001. |
[10] |
阳凡林, 刘经南, 赵建虎.
多波束测深数据的异常检测和滤波[J]. 武汉大学学报信息科学版, 2004, 29(1): 80–83.
YANG Fanlin, LIU Jingnan, ZHAO Jianhu. Detecting Outliers and Filtering Noises in Multi-beam Data[J]. Geomatics and Information Science of Wuhan University, 2004, 29(1): 80–83. DOI:10.13203/j.whugis2004.01.018 |
[11] |
何义斌, 吴书帮, 谢洪燕, 等.
多波束异常测深数据检测方法实践[J]. 测绘科学, 2004, 29(1): 50–52.
HE Yibin, WU Shubang, XIE Hongyan, et al. Study of Detection and Filter of Outlier on the Sounding Data of MES[J]. Science of Surveying and Mapping, 2004, 29(1): 50–52. DOI:10.3771/j.issn.1009-2307.2004.01.017 |
[12] | HARE R. Depth and Position Error Budgets for Multibeam Echosounding[J]. International Hydrographic Review, 1995, 72(2): 37–69. |
[13] | CALDER B R, MAYER L A. Automatic Processing of High-rate, High-density Multibeam Echosounder Data[J]. Geochemistry, Geophysics, Geosystems, 2013, 4(6): 253–278. DOI:10.1029/2002GC000486 |
[14] | CALDER B. Automatic Statistical Processing of Multibeam Echosounder Data[J]. International Hydrographic Review, 2003, 4(1): 53–68. |
[15] | CALDER B R, MAYER L A. Robust Automatic Multi-beam Bathymetric Processing[C]//U.S. Hydrographic Conference (US HYDRO). Norfolk, VA, USA: [s.n.], 2001. |
[16] |
王德刚, 叶银灿.
CUBE算法及其在多波束数据处理中的应用[J]. 海洋学研究, 2008, 26(2): 82–88.
WANG Degang, YE Yincan. The Theory of CUBE Algorithm and Its Application in the Processing of Multi-beam Data[J]. Journal of Marine Sciences, 2008, 26(2): 82–88. DOI:10.3969/j.issn.1001-909X.2008.02.012 |
[17] | QI Xiaodi, ZHOU Xinghua. Automatic Processing of the Multibeam Echosounding Data Based on CUBE Algorithm[C]//IEEE International Conference on Information and Automation. Yinchuan, China: IEEE, 2013. |
[18] |
贾帅东, 张立华, 曹鸿博.
基于CUBE算法的多波束水深异常值剔除[J]. 测绘科学, 2010, 35(s1): 57–59, 94.
JIA Shuaidong, ZHANG Lihua, CAO Hongbo. A Method for Eliminating Outliers of Multibeam Echosounder Data Based on CUBE[J]. Science of Surveying and Mapping, 2010, 35(s1): 57–59, 94. DOI:10.16251/j.cnki.1009-2307.2010.s1.022 |
[19] |
黄辰虎, 陆秀平, 侯世喜, 等.
利用CUBE算法剔除多波束测深粗差研究[J]. 海洋测绘, 2010, 30(3): 1–5.
HUANG Chenhu, LU Xiuping, HOU Shixi, et al. Study on Detecting Outliher of Multibeam Sounding Based on CUBE Algorithm[J]. Hydrographic Surveying and Charting, 2010, 30(3): 1–5. DOI:10.3969/j.issn.1671-3044.2010.03.001 |
[20] | VÁSQUEZ M E. Tuning the CARIS Implementation of CUBE for Patagonian Waters[D]. Fredericton, New Brunswick: University of New Brunswick, 2007. |
[21] |
郭军, 刘胜旋, 关永贤, 等.
浅水多波束系统SONIC 2024在码头测深中的应用[J]. 测绘工程, 2016, 25(7): 46–50, 56.
GUO Jun, LIU Shengxuan, GUAN Yongxian, et al. Application of SONIC 2024 Multibeam System to the Terminal Survey[J]. Engineering of Surveying and Mapping, 2016, 25(7): 46–50, 56. DOI:10.19349/j.cnki.issn1006-7949.2016.07.010 |