论文专区▏多波束水体影像的归位算法研究与实现

【编者按】多波束水体影像近年来得到了广泛的应用，但由于缺少格式标准，各多波束厂商所支持的格式各不相同，给进一步的数据处理带来了一定的困难。基于Matlab平台，通过对ALL格式多波束测深数据文件与WCD格式多波束水体影像文件进行了详细的解析，并且研究了最近邻插值、改进最近邻插值、二维线性插值与三次多项式插值等多种算法，实现了多波束水体影像的回放显示与水体异常提取，并且实现了回放速度的控制，极大地提高了资料处理效率。本文发表在《海洋测绘》2016年第1期上，现编发给朋友们阅读了解。刘胜旋，1972出生，男，广西来宾人，广州海洋地质调查局，教授级高级工程师，主要从事海洋测深资料处理及方法研究。

文/刘胜旋 关永贤 宋永志

一 引言

多波束测量发展的早期，由于技术、硬件设备的原因，往往只注侧重于水深测量方面的研究与发展。近年来，随着技术的发展、硬件运算速度的提高、磁盘存储空间的快速增长，多波束测量开始了三位一体的全海底空间测量：海底测深、海底反向散射强度与水体的反向散射强度，其中第三种测量即我们通常所称的多波束水体影像（也称水柱影像，Water Column Image）测量。

近年来，利用多波束系统进行水体影像观测，并应用于沉船^[1-2]、渔群分布^[3-4]，甚至应用于冷泉^[4-5]和海洋内波^[6-7]的研究，已成为一个新兴的研究手段并取得了很好的效果。

目前，世界上各大知名的多波束厂商新一代的多波束系统均具有记录水体影像数据的功能。但各多波束系统所记录的水体影像数据的格式跟测深数据格式类似，种类众多，互补兼容。例如挪威Kongsberg多波束水体影像数据扩展名为WCD、美国Teledyne Reson多波束水体影像数据扩展名为S7K、德国Teledyne Atlas多波束水体影像数据扩展名为ASD、德国Elac-nautic多波束水体影像数据扩展名为WCI。这给数据管理、开发利用和后处理成图均带来各种困难。

国内一些学者，例如高金耀，曾将各种格式原始多波束测深数据进行了转换与统一，开发了微机版的多波束海底地形电子成图系统MBChart，实现各种数据的综合处理^[8]。张哲^[9]等进行了基于数据库的多波束测深数据处理技术研究。陆秀平^[10]等则自主研发了一套多波束测深数据处理及管理系统。基于这种思路，在国家高技术研究发展计划（863计划）的支持下，项目组开发了一套多波束水体影像数据处理系统，能够以不同的速度回放多波束水体影像，以及提取水体中的异常特征体，以ASCII文本格式输出供第三方制图软件使用，本系统基于Matlab平台开发。

二 原始数据格式解析

由于多波束数据量巨大、内容非常丰富，即使以二进制格式来存储，每小时也高达数十至数百兆比特。为节省数据存储空间，多波束系统几乎全部采用长度为1、2、4的无符号整数或有符号整数以二进制方式存储数据（除了实时参数数据包外）^[11]。

ALL格式是挪威KongsbergEM系列多波束系统采用的数据存储格式，该类型产品在多波束市场中占有很大的份额。该格式记录数据方式灵活，所有不同类型的数据均以数据包为单位，以类似积木的方式顺序存储（见图1）。每一数据包的大小随实际数据量的大小而变，不会出现无数据时以零填充浪费存储空间的现象，在节约空间上更胜一筹，因此ALL格式在EM多波束系列中的应用十分广泛^[12]。

数据包主要类型有安装参数数据包、定位数据包、姿态传感器数据包、水深数据包、时间旅程与波束指向角数据包、声速剖面数据包、水底影像数据包、水体影像数据包等等。所有数据包的前4个节均用于存储该数据包的大小。然后以STX（二进制值为02）为开始标志，依次是数据包类型、仪器型号、日期、时间、脉冲顺序号、仪器序列号、核心数据包（依数据类型的不同而不同），最后以ETX（二进制值为03）为结束标志（图1下）。

目前，EM多波束系统记录的以ALL为后缀的数据文件，只包含测深、定位、日期、时间、姿态传感器、罗经等数据包，而水体影像数据则由于需要注册许可及数据量巨大，另外记录于一个以WCD为后缀的同名文件名中。由于该文件仅仅包含水体影像数据包而不包含定位、罗经等数据包，所以无法对水体影像进行定位，实际应用时需配合相应的ALL文件同时读取。

三 算法研究

⒈ NNIA算法

多波束系统采集的回波数据点在空间上是呈扇形分布的(可视为极坐标形式)，具有波束号、波束角、采样点、采样率、声速等信息，而显示器上的显示像素则是以点阵分布(可视为直角坐标形式)，见图2。如果不进行任何坐标转换就直接进行显示，得到的则是发生巨大畸变的原始图像。为了能够在显示器上正确还原显示所采集的回波数据，必须在采样数据和终端显示图像之间进行必要的坐标转换。NNIA为Nearest Neighbor InterpolationAlgorithm的缩写，称为最近邻插值算法。该算法是最早的一维线性变换算法，它主要运用直角坐标与极坐标之间的几何变换关系来实现^[11]。通过极坐标与直角坐标的变换关系,可得到每个采样数据点在显示器上的对应值。

该算法简单来说就是根据水体影像像素寻找最邻近的目标图像像素，并利用水体影像像素对目标图像像素赋值。

⒉ 改进NNIA算法

改进NNIA算法^[14]就是根据目标像素的位置，反过来从水体影像中寻找一个与其最邻近的像素，利用水体影像像素对目标像素赋值。算法模型见图3，对于任一目标影像像素点F，首先定位出像素点对应的四个最相邻的回波点A、B、C、D（其中A、B属于同一波束，C、D也属于同一波束，且该两个波束为相邻波束）。只需判断距F点最靠近的是哪个回波点，就可以利用该点的值对F点进行填充。

⒊ R-Theta算法

R-Theta算法^[15]是双线性插值的一种特殊实现方式，其核心思想是在两个方向分别进行一次线性插值，实现分辨率的提高。

假设f_i,j、f_i+1,j、f_i,j+1与f_i+1,j+1为待求点F周围4个最邻近的波束像素点，θ为两波束之间的角度，见图4。先用插值法找出过度点F_h1和F_h2的值：

F_h1＝f_i,j×(θ－Δθ)＋f_i,j＋1×Δθ

F_h2＝f_i+1,j×(θ－Δθ)＋f_i＋1,j+1×Δθ

然后再用同样的方法，由F_h1、F_h2的值可求取F点的值：

F＝F_h1×(1-ΔY)＋F_h2×ΔY

ΔY可由F的极半径减去f_i,j的极半径所得。

R-Theta 算法将点与点之间由于距离的舍入或截断所到来的误差考虑在内，并按照一定的误差分布关系协调点与点的值的相关程度。反映到图像上即是对图像做了一次二维的平滑滤波（低通滤波），故图像能很好的保留原始信息，图像质量有很大提高。

⒋ 三次多项式插值算法

三次多项式插值^[16]属于一种非线性插值算法。三次多项式插值即是，在给定节点的情况下，按每4个节点作为某三次函数的根值点，由此可以确定该三次函数的系数。这样对于节点之间的点可通过代入三次函数得到。该算法不需要任何边界条件，使用非常方便。

假设存在四个节点x₀、x₁、x₂、x₃（x₀＜x₁＜x₂＜x₃）。对应节点的值为y₀、y₁、y₂、y₃。

设三次多项式为：

y＝ax³＋bx²＋cx＋d       ⑴

其中a、b、c、d待定参数。将节点代入方程，整理得系数向量A：

、F_h2的值可求取F点的值：

F＝F_h1×(1-ΔY)＋F_h2×ΔY

ΔY可由F的极半径减去f_i,j的极半径所得。

R-Theta 算法将点与点之间由于距离的舍入或截断所到来的误差考虑在内，并按照一定的误差分布关系协调点与点的值的相关程度。反映到图像上即是对图像做了一次二维的平滑滤波（低通滤波），故图像能很好的保留原始信息，图像质量有很大提高。

由此可确定三次多项式的表达式，对于任意点x(x∈[x₀，x₃])可代入多项式求出相应的值。

对于多波束水体影像数据，我们可以得到：首先，对于任意一个波束，其采样间距是相等的；其次，目标图像的尺度可以像素为单位，也即采样间距的单位长度1；第三，根据插值公式，可以按照分段方法进行插值，也即每4个波束采样点求出一个三次多项式函数，然后对位于采样点间的非采样点，代入三次多项式中进行求解。

通过上述条件，令x₀＝0，x₁＝1，x₂＝2，x₃＝3，将其代入⑵式得：

对于给定的像素点F，见图5，先得到它左右两相邻波束的波束号θ_i与θ_j，然后搜索每一根波束上F点上方的两个邻近采样点及下方的两个邻近采样点的像素影像强度值：y_i1、y_i2、y_i2与y_i3，从而构成一组同波束的四个邻近采样点，代入公式⑶求得矩阵系数A，然后再利用公式⑴得到F_h1点像素强度值。使用相同方法求取另一波束上F_h2点像素强度值。最后再次运用横向的线性插值得到像素点F处的值：F＝F_h1×(θ－Δθ)＋F_h2×Δθ。

⒌ 几种算法的比较

图6各小图分别为采用以上几种算法获得的多波束水体影像图。

NNIA算法的主要缺点有：⑴由于波束与波束之间本身存在间隙，这样对目标图像上象素点的填充并不是连续的，其间会有“空洞”产生，图像会出现所谓的莫尔条纹；⑵目标图像像素位置是整数值，所以坐标点计算存在舍入或截断误差。这样，原先回波点对象素点的填充，可能会被邻近回波点的值所覆盖，图像就会丢失信息产生失真。优点是采用这种算法的图像显示速度最快。

改进NNIA算法生成的图像不会出现空隙以及莫尔纹现象，但由于是将邻近的多波束水体像素点直接拿来赋值填充，图像的噪点较多，并且显示得较为生硬。

R-Theta算法由于进行了精确插值，图像的噪点得到有效的控制，图像显得较为细腻生动。缺点是该算法的运算量相当大，求取一个象素点的值要经过多次中间变量的运算之后才能得到，这对图像的实时显示很不利。

三次多项式插值算法相对于其它三种算法，图像显得最为精细，而且噪点控制也十分有效，缺点是运算量巨大，显示速度也最慢。

由上对比分析可知，以上四种算法各具优缺点：运算速度快的，显示质量较差；显示质量越好，运算速度则越慢。因此，本软件全部集成了以上四种算法，用户可根据自己计算机的运算速度、影像数据的采样率等实际情况，选择不同的算法，以达到最佳回放显示效率。

四 软件的主要功能

多波束水体影像数据处理软件除了具有快速浏览水体影像、提取水体异常特征体的基本功能，同时还具备多种坐标转换、回放速度控制、显示分辨率控制、背景噪音消除、暂停、缩放、平移、倒退显示等其它辅助功能，见图7。

软件可以工作于两种模式：练习模式与回放模式。默认运行模式为练习模式。该模式主要是针对在没有水体影像数据文件的情况下，利用单Ping的样例数据完成各种练习操作，进而熟悉软件个各项功能。软件首先自动读取内置的一Ping水体影像数据，经过坐标转换后，将水体影像显示于程序的主界面。因此，软件启动主界面的水体影像并非光栅图像，而是矢量图像，可以进行各种变换运算。

回放模式与练习模式的功能完全相同，它们的主要区别前者是动态显示，后者是静态显示。软件目前仅支持挪威Simrad公司的EM系列多波束系统采集的多波束水体影像数据，一条测线的水体影像数据由两部分组成：*.WCD与*.ALL。一旦打开水体影像数据文件后，软件则自动进入回放模式。

在回放模式下，通过按键盘的空格键即可暂停水体影像的回放,然后点击“拾取异常”即可通过多边形圈选的方式开展水体异常体拾取操作。输出数据的位置与输入水体影像文件相同。在拾取水体异常时，由于异常周围存在大量的信号水平较强的背景噪声，因此必须仔细沿着异常的轮廓圈选信号，否则强背景噪声的存在不仅会夸大异常体的规模，也大大降低了异常体信号提取的效率。因此，通过设置适合的背景噪声强度，将所有强度弱于指定值的信号均过滤消除，仅突出异常主体。经过背景噪声消除以后再进行异常提取，只需大致沿着异常体的边缘采集即可，可以提高提取效率。在实时回放时，如果归位后的水体影像分辨率大大高于显示器的分辨率，这将消耗大量的运算资源，从而影响回放速度，但显示质量却没有相应地提高。因此，软件提供了一个调节分辨率的工具，可以设置适当的分辨率，从而实现水体影像回放速度的控制。

其它的缩放、平移、倒退显示等辅助功能，由于篇幅显示，不再详细介绍。

五 结语

多波束水体影像在海洋测绘领域将会得到越来越广泛的应用,但由于各多波束厂商所支持的数据格式各不相同，这给进一步的处理与研究带来了困难。目前，一般只有厂商配套软件才能够完成影像回放以及异常提取等操作，使用范围受到很大的限制。本文通过解析ALL格式的多波束测深数据与WCD格式的多波束水体影像数据，并研究了多种插值算法，实现了水体影像数据的回放与异常提取,为进一步的研究奠定基础。

在实际应用时，可根据计算机的运算速度、水体影像数据的质量等情况选择不同的算法，可实现显示速度与显示质量之间获得一个平衡点。或者将水体影像调整到适合的分辨率，从而实现了显示速度的调整，能够很好地提高了资料处理效率。

参考文献

[1]John E. Hughes Clarke,MikeLamplugh,Kalman Czotter.Multibeam Water Column Imaging :Improved WreckLeast-Depth Determination[C]//Canadian Hydrographic Conference.May,2006.

[2]John E.HughesClarke.Applications of multibeam water column imaging for hydrographicsurvey[C]//The Hydrographic Journal.April,2006.

[3]Keir Colbo, TetjanaRoss,Craig Brown,et al.A review of oceanographic applications of water columndata from multibeam echosounders[J].Estuarine, Coastal and ShelfScience.2014(145)41-56.

[4]Gary D,Melvin.NormanA,Cochrane.Multibeam Acoustic Detection of Fish andWater Column Targets at High-FlowSites[J].Estuaries and Coasts 2015(38)(Suppl 1):227–240.

[5]A.T.Jones,J.Greinert,D.A.Bowden,etal．Acousticand visual characterisation of methane-rich seabed seeps at Omakere Ridge onthe Hikurangi Margin,New Zealand[J].Marine Geology;2010(272) 154-169.

[6]刘胜旋,崔秀云,海洋内波对多波束测深的影响[J].海洋测绘,2012,32(6):27-29.

[7]阳凡林,韩李涛,王瑞富,等.多波束声纳水柱影像探测中底层水域目标的研究进展[J].山东科技大学学报,2013,32(6):75-83.

[8]高金耀,方银霞,徐赛英.原始多波束数据的格式转换与统一[J].海洋通报,2002,21(6):68-74.

[9]张哲,顾锦华,谢丛涛,等.多波束数据处理系统技术研究[J].海洋测绘,2010,30(02):52-54.

[10]陆秀平,黄谟涛,翟国君,等.多波束测深数据处理及管理系统设计与开发[J].海洋测绘，2006,26(06):1-5.

[11]Kongsberg Simrad OperatorManual.EM series datagram formats[ S/OL].[850-160692/M January2010].http://www.kongsbergsimrad.com/.

[12],,,等．EM300X系列多波束ALL数据格式解析[J].:16-22.

[13]朱玉颖,冯佳梅,姚远程.图像几何变换的DSP算法研究与实现[J].微型机与应用,2012,29(4):43-48.

[14]周健,钱进．B超图像的计算机实时成像研究[J].声学技术,2003(3):195-198.

[15]尚兆梅,陈波,彭勇.R-theta算法的一种快速实现方式[J].科学技术与工程,2010,10(31):7803-7806.

[16]周健,钱进.一维非线性插值在重建B型超声图像中的应用[J].苏州大学学报(自然科学版),2003,19(2):55-58.

■论文专区的文章均为在《海洋测绘》上刊发的论文，若其他公众平台转载，请备注论文作者，并说明文章来源，版权归《海洋测绘》所有。