基于声图像的海底地形边界提取算法

王媛媛，郭延恩，施国全，韦俊霞，夏顺仁

（1.浙江大学 生物医学工程教育部重点实验室，浙江 杭州310027；2.浙江省心脑血管检测技术与药效评价重点实验室，浙江 杭州310027；3.杭州中船重工第七一五声学研究所，浙江 杭州，310012）

近年来，随着海洋资源开发、海洋环境监测评价、浅海和海岸建设的日益增多，与海洋工程相关的海底声学地层剖面探测也日渐频繁.众多探测手段中，海底地形的提取尤为重要，它不仅能提供简便直观的海底样貌、岩层特征，还可以为沉积物地质属性的识别和声学参数的提取提供参考信息［1，2］.探测手段的发展促进了海底成像的多样化，数字图像处理技术的引入使得海底各反射地层显示更加清楚直观，较为成熟的图像处理手段也为海底复杂地形的提取提供了强有力的保障.

海底探测仪器（Chirp Sub-bottom profiler）发射线性调频声波，经海底各存在声阻抗差异的地层反射，由换能器接收到的能反应水体深度（垂向）信息的信号称为探测回波信号，此回波信号可以包含原始探测信号、相关信号、解析信号、包络信号以及TVG 增益下的包络信号等［3］.此前的研究多是直接基于探测回波信号的振幅、频率、相位等变化来提取相关界面.左国平等［4］采用滚动时窗来计算回波信号能量比值，从比值特征中判断海底反射层位置.罗进华等［5］在能量比值法基础上，限制时窗滚动的范围和相邻两列信号（以ping为计量单位）之间反射的时间差，提取连续边界.刘秀娟等［6］采用最大振幅法实现海底反射信号的识别，平滑后得到边界信息.丁维凤等［7-8］根据剖面的形态特征采用间隔数据跳跃搜索来完成提取过程，提高剖面解释的效率，比较能量比值法与互相关分析法在剖面提取过程中的表现，采用人工交互来修正提取结果.

总体上看，目前常见的海底界面提取方法有能量比值法、相关分析法、最大振幅法等，此类方法均直接基于探测回波信号，回波信号精度高，信息量丰富.但回波信号均是单ping存储，此类方法多是直接基于单个ping或者2 个ping，很难灵活利用多ping甚至整个空间域的关联信息，故基于此的方法多半采取交互提取，费力费时，且存在个体差异.鉴此，本文将回波信号压缩至8位灰度图像，采用基于声图像的连通域标记和动态规划寻优相结合的算法自动提取海底底形边界.

1 基于图像的海底地形边界提取分析

1.1 常见的边界提取算子

工程上提取地形边界的基本要求是准确、连续、单像素.常见的边界提取算子，如Sobel算子、Roberts算子、Canny算子等，均能通过较明显的灰度变化识别边界，但极易受到噪声影响.实际的海底声图像（如图1所示，水平轴为回波信号的个数n，单位：ping，垂直轴为垂直探测的深度δ，单位：m）对比度差，灰度较为集中且噪声类型多变［9］.因此直接应用上述算子，会产生大量粗细不一的轮廓线，识别出的边界线会存在很多断点.形态学操作可以连接断点，但同时也会带来很多伪迹，且鲁棒性差.

图1 一幅典型的海底声学图像Fig.1 Typical acoustic image of seabed

本文首先计算图像的连通域，选取满足连通域面积条件的某一边界点，进而以此点作为动态规划算法（Dynamic programming，DP）寻优的起点，与此同时引入了能量函数的概念，实施寻找能量函数最小值为目的迭代DP 算法，从而得到单像素宽度的连续边界.

1.2 边界起点的自动识别

连通域的检测和标记［10-12］是很多图像处理和机器视觉必不可少的处理步骤，以Se表示结构元素，作用原理可以表述如下：

对于二值化图像I 内任一点未被标记的点p，令X0，p＝p，采用迭代式（1）生成连通域Yp的所有元素，遍历结束后，集合Y＝｛Yp｝即为标记结果.

对图1 所示的原始声图像采用固定阈值二值化，而后进行形态学处理所得的结果如图2（a）所示，可以看出，图像上半段出现了一整条带状区域，这是声束探测中旁瓣干扰.旁瓣干扰与地形边界的信号强度相仿，且二者距离很近，甚至可能重叠.本文提出了一种灰度映射的方法，对一幅大小为M×N（M 为宽度，N 为高度）的图像，定义灰度映射函数gi如式（2），其中I （i， j） 为当前像素的灰度：

经过式（2）的灰度映射，受旁瓣影响的区域灰度趋于一致，与不受影响的其他区域有较为明显的灰度差值.边界起点识别的步骤如下：

Step 1：计算原图的灰度映射函数，并提取映射函数的梯度▽gi.

Step 2：从图像上观察，旁瓣干扰的垂向范围不超过探测深度的20%，故在0.2×h 高度范围内，根据梯度差最大原则，即可去除图像上方的带状区域的干扰.h为图像垂向分辨率，取值为1 024.

Step 3：将去除干扰之后的图像二值化，而后进行连通域标定.

Step 4：扫描标定后的图像，计算各连通域的面积，满足式（3）时，记录此连通域的标号.

Step 5：找到Step4所得连通域的起点，以此点作为2.3节DP寻优的起点.

不等式（3）中C 是经验常数，表示可以认定为边界的最小连通面积值，具体应用中可根据实际情况适当微调.本文取值2 000，是经过多次实验确定的一个经验值.若C 值较小，式（3）可能找到较小的噪声区连通域；若C 较大，式（3）可能会跳过有效边界，找到较大的噪声区连通域.

按照上述步骤，得到的起点位置用白色十字标注如图2（b）所示，可以看出，本文提出的算法能够有效地避开旁瓣信号影响，找到正确的边界起点.

图2 连通域标记结果图Fig.2 Result of connectivity labeling

1.3 动态规划提取边界

DP的核心思想［13］是将问题分解为若干相似的子问题来解决，从而降低计算量和复杂度.对于一个复杂的问题，DP不直接求解，而是从易到难地解决此问题的子集，这些子集被称为“状态”.对于每个子问题，DP只保留最好的结果作为当前值，而后解决下一个子问题，直至解决该问题本身.每一个子问题解决方法类似，这种一致的方法被称为“状态方程”，因此DP算法的行为模式取决于状态和状态方程.

概念1.韦伯最小分辨值Iweibo，Coren等［14］依据韦伯感光模型和人眼视觉分辨率提出了韦伯亮度差感知函数定义，见式（4），

此函数代表了不同亮度层次上，人眼能够分辨的最小亮度差.当噪声较多时，则用8邻域内的灰度均值代替当前像素点的灰度值I （i， j）［15］.

概念2.8 邻域下三角平均灰度Iavepart，如图3阴影部分所示.观察原始的海底图像可以看出，待识别的边界位置较暗，边界之上总是有一段较亮的海水区域.不论边界如何走向（左上倾斜、右上倾斜或者水平走向），边界位置下三角内的灰度均值总是相对较小.

图3 8邻域下三角区域示意图Fig.3 Diagram of lower triangle in 8-neighborhood

概念3.像素的能量函数e.此函数综合了多种特征信息，是海底地形的识别函数，也是整个DP算法的核心函数.当能量函数值非负时，此值越小，越接近海底地形边界，能量函数定义如式（5），

式中：Iave为当前像素所在8邻域内的灰度均值，gy为垂直方向上的梯度，gx为水平方向上的梯度，dy为当前点与前一个边界的位置差值，即垂直方向上的位置差值.

能量函数的物理含义可以表述为：当前点灰度、8邻域下三角灰度均值都不应超过200，8邻域下三角灰度均值不应该超过整个8邻域的灰度均值，当前点所在位置的垂直方向梯度应该超过对应的韦伯分辨率.如不满足上述条件，将当前点的能量置成-1，表示当前点不可能成为边界点.对于满足条件的像素点，则准确地计算能量，公式中的参数为实验所得，公式中 的Ii，j、gy、gx、dy分别表 示了4个趋势：灰度最小化，垂直方向灰度差异最大化，水平方向灰度一致化，边界平稳化.由此可以明确，在非负的情况下，能量越小，越接近边界.

进一步地，定义DP的状态及状态方程如下：

定义1.状态 J （i） ，DP算法的目的是在每一列像素中寻找能量函数的最小值.i表示列数，图像左起第1列i＝0.J （i） 则表示第i 列中能量函数取得非负最小值的深度位置，即当前列的地层位置.

定义2.状态方程，状态方程表述前后2种状态之间的迭代关系，如式（6），

式 中：j∈ ［J （i- 1） -50，J （i- 1） ＋50］ ，表 示DP算法寻优的范围，EMax为人为设定的能量阈值，如果最小的非负能量值超过此阈值，则认定在指定范围内边界缺失.此时，采用修正项J （i- 1） 作为当前列的边界位置.

本算法以连通域标定的结果作为边界提取的起点，自起点向两侧逐列遍历，算法作用的示意图及流程图如图4所示.

2 算法评测与讨论

为客观评测文中所提出的算法，选用中船重工第715研究所于某海域采集到的海底地形数据.原始数据为32位浮点型，采用刘冬梅等［16］提出的对数域压缩方法将数据压缩至8位图像中，尽可能减少数据损失.

本文实验平台是Visual studio 2010，基于QT和VTK 的混合编程.考虑到实际的海底地形情况比较复杂，本文收集了较为典型的36例实测数据样本，其中13 例 约 为2 000ping，16 例 在2 000～3 500ping之间，其余数据较大，10 000ping左右.样本垂向分辨率均为1 024个像素.

采用文中算法对上述36例数据进行测试，结果表明：32例可以正确识别边界上的点，剩余4 例由于边界周围噪声较多识别效果欠佳.识别正确率为88.89%，在一定程度上表明了文中算法的有效性以及C 取值的合理性.

进一步地，以专业技术人员手画边界作为金标准，采用准确率A、平均欧氏距离D 和相似度S3种量化指标来评测实验结果.以C1，C2分别代表金标准曲线和本文算法提取的曲线，C1-C2表示2条曲线上每一列的纵坐标对应相减，μ，σ分别表示C1-C2的均值和标准差，length（）是求取向量长度的函数，对3种评测指标的定义如下：

提高实习生病历书写水平，科教科组织成立了实习生病历书写培训小组，按照1∶10比例，即一名教师带教10名同学，每月对学生进行病历书写的培训，同时，针对培训教师和实习同学提出的关于病历培训方面存在的问题，邀请病案室负责人进行答疑，针对存在的问题进行专项整改与训练，提高了培训效果。出院（“1+1”）考试：对将毕业的新疆医科大学实习生进行出院理论知识、大病历及各项技能操作的考试。

图4 DP算法示意图及流程图Fig.4 Diagram and flow chart of DP method

假设手画边界存在±0.3%的误差，在垂向1 024分辨率条件下，式（7）定义的准确率是落在金标准±3个像素位置内的边界点概率，表征了实验结果与金标准的整体差别；式（8）定义的平均欧式距离是在垂直方向距离归一化时，每一列实验结果与金标准距离的均值，反映的是二者在空间位置上的差距；曲线相似度定义由三角形相似定义类比得出，当一条曲线的任2点连线长度与第2条曲线上相对应的2点连线长度成比例，这2条曲线相似，对于实际问题中不完全相似的2条曲线，江浩等［17］引入置信区间的概念来衡量2条曲线相似的程度，为了保证实验结果的严谨性，本文取置信区间最小值μ-σ，μ＋［］σ ，以落在此区间之内的比例作为相似度，如式（9）.3种评测指标的数值结果如表1所示.

表1数据从编号1至36，数据文件的大小依次递减.表1对应的统计分析结果如图5所示，准确率与相似度盒图如图5（a）所示，平均欧式距离分布图如图5（b）所示.综合二者可以看出本文方法在所有实验数据上的评测结果良好，准确率分布集中性略差，但中值在95%以上；相似度中值接近95%，且比较集中.平均欧氏距离是归一化的结果，实验数据的垂向分辨率为1 024（≈103），故而在实际图像中，垂直方向平均距离最大不超过1个像素，最小则约等于0.对于36例测试样本，采用本文算法提取的海底地形边界（未经人工干预）准确率均值为92.6%，相似度均值为93.3%，与金标准的欧式距离均值为0.23个像素.准确率和相似度较高，说明2 条曲线走势及弯曲程度近似一致，平均欧式距离小，说明实验结果曲线与金标准在整体位置上距离很近，本文方法的鲁棒性较好.本文的实验结果如图6所示，可以看出即使在海底地形多变的情况下也可准确提取边界.

表1 3种评测指标数值分析Tab.1 Numerical analysis of 3evaluation indices

图5 评测指标统计分析图Fig.5 Statistical analysis of evaluation indices

图6 经典的海底地形边界提取结果Fig.6 Extracted result of typical seabed terrain boundary image

1）部分信号缺失

探测频率变更或仪器故障时，原始数据可能出现部分声学信号缺失，在图像中则表现为部分纵列灰白，如图7所示.对于这种情形，本文方法也较为准确地检测出了边界，且检测出的边界位置较为连续，可信度较高.

图7 部分信号缺失的提取结果Fig.7 Result in case that part signals are missed

2）地形边界起伏较大

滑坡或崩塌是海底较为常见的灾害性地形，在图像上表现为较大的边界起伏，如图8所示.此类地形的误判会给海洋工程带来诸多困扰，故其提取也一直备受关注.

图8 地形边界起伏较大的提取结果Fig.8 Result of terrain boundary being largely fluctuated

本文提出的DP算法状态方程中涉及到一个寻优步长，表征着第i列到第i＋1列边界上下移动的最大可能.本文步长取值50，实验证明，这个步长对此类海底图像是有效的，既可以识别出起伏较大的边界，也不至于受到较远处噪声点的干扰.

3）地形边界与旁瓣信号交叠

部分地形边界会与旁瓣信号交叠，如图9所示.2.2节中灰度映射操作会将旁瓣信号所在区域全部抹去，因此交叠区内本文方法无法准确提取此类边界.但在交叠区外，本文方法可正常寻优.

图9 地形边界与旁瓣信号有重叠的提取结果Fig.9 Result in case that signals of side lobe and terrain boundary are overlapped

4）关于出错边界的人工干预

综合如上评测可得，本文所述方法可在保证较高准确程度的前提下，实现海底地形边界的自动提取.在实际应用中，用户可以干预出错边界，调整边界走向，如图10所示.

图10所示图像海底地形边界与旁瓣信号有交叠，且噪声较多.在连通域标记过程中，大面积噪声区域率先满足不等式（3），被标记为起点，如图10（a）.错误起点引发的错误边界如图10（b）右侧矩形框标注，DP寻优具有较好的容错性，如果错误起点在正确边界位置的寻优步长之内，则在小范围波动之后，DP会找到正确的边界位置；图10（b）左侧的矩形框则表示旁瓣信号与地形边界重叠而导致的出错情况.在实际应用中，所有边界错误均可人为修正.本文设计的人工干预方法如下：记录鼠标在屏幕上单击的坐标位置，以此位置作为DP 寻优的新起点，按照原有的状态方程迭代，至新的边界位置与原有边界位置重合，即停止调整.另外，按住鼠标拖拽也可精确调节边界位置，按照以上方法，调整后的边界如图10（c）所示.

图10 人工干预过程图示Fig.10 The flow of manually guided processing

3 结 语

针对自动提取海底地形的实际需求，本文提出了一种基于连通域标记与动态规划（DP）相结合的方法.采用连通域检测较为准确地定位了边界大致位置，进而以边界上的单个点作为DP 算法寻优的起点，并在DP算法中引入了能量函数的概念，同时以该能量函数最小为准则达到了提取边界的目的.所提出算法的特点是：有效利用了空间域信息，较好地运用了动态规划寻优的优势，从而自动完整地提取整条边界，鲁棒性好.目前文中所提出的算法已成功应用于作者单位所研发的海底浅地层剖面仪中，取得了良好的实测效果.当然，对于边界提取欠正确的情形，用户可利用人机交互接口视需要进行不同程度的调节，准确提取出边界.

目前的算法可以有效地提取出海底地形边界，未来工作需进一步考虑提取除海底表层以外的其他反射层的边界.

（

）：

［1］DJIKPESSE H，SOBREIRA J F F，HILL A，et al.Recent advances and trends in subsea technologies and seafloor properties characterization［J］.The Leading Edge，2013，32（10）：1214-1220.

［2］陈晓辉，张训华，李铁刚，等.渤海海峡及周边海域浅地层结构及地层声速的拾取［J］.海洋地质与第四纪地质，2012，32（1）：69-75.CHEN Xiao-hui，ZHANG Xun-hua，LI Tie-gang，et al.Shallow stratigraphic sequences and their acoustic velocity in the Bohai Strait and surrounding areas［J］.Marine Geology ＆Quaternary Geology，2012，32（1）：69-75.

［3］HENKART P.Chirp sub-bottom profiler processing-A Review explains how chirp signals may be recorded as correlates，analytic or envelope［J］.Sea Technology，2006，47（10）：35-38.

［4］左国平，王彦春，隋荣亮.利用能量比法拾取地震初至的一种改进方法［J］.石油物探，2003，43（4）：345-347.ZUO Guo-ping，WANG Yan-chun，SUI Rong-liang.An improved method for first arrival pickup using energy ratio［J］.Geophysical Prospecting for Petroleum，2003，43（4）：345-347.

［5］罗进华，丁维凤，潘国富.改进的滚动时窗法实现海底浅地层剖面反射层位自动拾取的研究［J］.物探化探计算技术，2008，30（5）：363-367.LUO Jin-hua，DING Wei-feng，PAN Guo-fu.Research on automatic picking of the reflection horizons of subbottom profile based on the improved moving time-window method［J］.Computing Techniques for Geophysical and Geochemical Exploration，2008，30（5）：363-367.

［6］刘秀娟，高抒，赵铁虎.浅地层剖面原始数据中海底反射信号的识别及海底地形的自动提取［J］.物探与化探，2009，33（5）：576-579.LIU Xiu-juan，GAO Shu，ZHAO Tie-hu.The recognition of the seabed reflection signal and the automatic pickup of seabed topography from the original data of sub-bottom profile ［J］.Geophysical and Geochemical Exploration，2009，33（5）：576-579.

［7］丁维凤，罗进华，来向华，等.浅地层剖面交互拾取解释技术研究［J］.海洋科学，2008，32（9）：1-6.DING Wei-feng，LUO Jin-hua，LAI Xiang-hua，et al.The research of interactive interpretation picking for subbottom profile［J］.Marine Sciences，2008，32（9）：1-6.

［8］丁维凤，潘国富，苟铮慷，等.基于能量比与互相关法的地震剖面反射同相轴交互自动拾取研究［J］.海洋学报，2012，34（5）：87-91.DING Wei-feng，PAN Guo-fu，GOU Zheng-kang，et al.The research of interactive auto pickup of seismic enents based on energy ratio and cross-correlation ［J］.Acta Oceanologica Sinica，2012，34（5）：87-91.

［9］宋坤坡，夏顺仁，徐清.考虑小波系数相关性的超声图像降噪算法［J］.浙江大学学报：工学版，2010，44（11）：2203-2208.SONG Kun-po，XIA Shun-ren，XU Qing.Algorithm considering correlation of wavelet coefficients for ultrasound image denoising［J］.Journal of Zhejiang University：Engineering Science，2010，44（11）：2203-2208.

［10］HE L，CHAO Y，SUZUKI K，et al.Fast connectedcomponent labeling ［J］.Pattern Recognition，2009，42：1977-1987.

［11］CAO Y，HAO X，XIA S.An improved region-growing algorithm for mammographic mass segmentation［C］∥Sixth International Symposium on Multispectral Image Processing and Pattern Recognition.Yichang：International Society for Optics and Photonics，2009：74971O-1-74971O-7.

［12］GONZALEZ R，WOODS R.Digital image processing［M］.Beijing：Publishing House of Electronics Industry，2002：434-437.

［13］WAGNER D B.Dynamic programming ［J］.The Mathematica Journal，1995，5（4）：42-51.

［14］COREN S，WARD L M，JENNS J T，Sensation and Perception（4th ed）［M］.Fort Worth，USA：Cold Spring Harcourt Brace College Publishers，1994.

［15］郭礼华，李建华，杨树堂.综合区域和边界信息的图像自适应分割技术［J］.上海交通大学学报，2005，39（4）：522-526.GUO Li-hua，LI Jian-hua，YANG Shu-tang.Adaptive image segmentation combining region and boundary information［J］.Journal of Shanghai Jiaotong University，2005，39（4）：522-526.

［16］刘冬梅.高动态范围图像显示算法的研究［D］，上海：上海交通大学，2009.LIU Dong-mei.Study on high dynamic range image displaying algorithm［D］，Shanghai：Shanghai Jiaotong University，2009.

［17］江浩，褚衍东，郭丽峰.曲线形态相似性的定义与度量［J］.云南 民 族 大 学 学 报：自 然 科 学 版，2009，18（4）：316-318.JIANG Hao，CHU Yan-dong，GUO Li-feng.Definition and measurement of shape similarity for curves［J］.Journal of Yunnan University Nationalities：Natural Sciences Edition，2009，18（4）：316-318.