基于聲圖像的海底地形邊界提取算法

王媛媛，郭延恩，施國全，韋俊霞，夏順仁

（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］陳曉輝，張訓華，李鐵剛，等.渤海海峽及周邊海域淺地層結構及地層聲速的拾?。跩］.海洋地質與第四紀地質，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］劉秀娟，高抒，趙鐵虎.淺地層剖面原始數據中海底反射信號的識別及海底地形的自動提?。跩］.物探與化探，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.