任秀艷,陳夢英,許偉杰
(1.中國科學院聲學研究所東海研究站,上海201815;2.中國科學院大學,北京100190)
尾流作為一類水下聲目標,具有十分重要的軍事意義,在海戰之中起著重要的作用:(1)尾流能夠影響探測設備的正常工作;(2)尾流是魚雷進行探測、攻擊目標的重要物理場;(3)尾流能作為潛艇規避的屏障[1]。因此,對各類艦船尾流的聲學特性,包括尾流的產生、持續時間和空間的分布、尾流中氣泡升浮的速度以及氣泡散射與吸收的特征、水下聲學尾流成像等方面進行研究日益重要[2]。
既準確又直觀的尾流圖像可以為下一步的尾流特征研究奠定基礎,而在尾流多波束聲吶的原始圖像中,不僅有尾流回波,還有海面回波反射以及其他噪聲干擾,因此需要采取有效的圖像分割方法把尾流從中提取出來,從而分析尾流特征。對尾流圖像處理時,可把尾流圖像看作三維表面,即x、y軸分別為圖像的橫、縱向矩陣列數,z代表強度的灰度。圖像邊緣分割已有較長的歷史,邊緣是圖像特性發生變化的位置,不同的圖像灰度不同,邊界處會有明顯的邊緣出現,利用這種特征可以進行圖像的分割。
多波束聲吶測量系統由船載分系統和水下分系統兩大部分組成,水下系統的水密艙中包含多波束聲吶。課題組于海上進行了多次實驗,采用的坐底式多波束測量系統圖如圖1所示。

圖1 測量系統示意圖Fig.1 Schematic diagram of the measuring system
原始尾流圖像如圖2所示。當主瓣收到的尾流信號與旁瓣收到的海面反射信號同時到達時,海面反射信號會干擾尾流信號。

圖2 原始尾流圖像Fig.2 Original wake image
圖2中,海平面在高度約為30 m的中間亮條處,海平面上方的虛擬強度圖像由多波束旁瓣引起,在此不做討論。根據分析,亮環即為海面回波反射干擾,中間條狀區域即為尾流圖像區域。尾流圖像邊緣提取的目的是將亮環剔除,并將條狀區域的尾流邊緣提取出來。
利用傳統的濾波方法難以去除聲吶尾流圖像中海面反射強度干擾,因此本文用背景相減法來去除海面回波干擾,構造了一個包含“亮環”噪聲并且沒有尾流圖像的背景環境,通過聲吶圖像與構造的“亮環圖像”相減得到的尾流圖像如圖3所示。
以圖3為原始圖像進行尾流圖像邊緣分割,下面分別介紹幾種邊緣分割算法,并給出這幾種邊緣分割算法提取到的尾流邊緣。

圖3 背景相減法濾波結果Fig.3 Results of background subtraction filtering
基于一階的邊緣檢測算子中,Canny算子是公認的具有良好檢測性能的一類邊緣檢測算子,因此本文利用它對尾流圖像進行邊緣提取。
首先,用高斯濾波器對尾流圖像進行平滑處理。然后,利用普通的一階邊緣檢測算子,如Sobel算子、Roberts算子等進行梯度計算,本文中選用的Sobel算子。再利用非極大值抑制的方法判定邊緣梯度值,在尾流檢測的過程中,進行的非極大值抑制是在梯度方向上的非極大值抑制,實際處理的尾流圖像的像素點是離散的二維矩陣,對于某個監測點,梯度方向兩側的點不一定存在于二維矩陣中,本文通過插值來解決這個問題。將某點梯度幅值與插值作比較,若此處梯度幅值大于兩個插值,則此點就是局部極大值,為尾流邊緣點。最后,進行閾值設置、拼接邊緣,具體思路是選取兩個閾值,認為小于低閾值的點是弱邊緣,大于高閾值的點是強邊緣,介于中間的像素點需要進一步檢查[3]。
邊緣檢測過程中尾流圖像閾值的設置很關鍵,也是難點,合適的閾值既要盡量保存尾流邊緣信息也要抑制噪聲,避免引入錯誤邊緣。Matlab自帶的Canny算子就是通過直方圖來設定閾值的。它默認選擇直方圖中幅值大小排在前 70%的最后一個閾值為高閾值,低閾值為高閾值的 0.2倍。本文在利用Canny算子處理實驗數據時借鑒了此方法,先將所有梯度值進行排序,經過多次測試,高閾值從高到低開始試驗,保證噪聲逐漸增大,因為閾值設定越高,噪聲抑制效果越好。最終選擇排在整個梯度值序列中第 98% 處的幅值作為高閾值,低閾值約為高閾值的0.2倍。這里設定的閾值是經過多次實驗后選取的。圖4、5為利用低閾值和高閾值檢測到的弱、強邊緣。

圖4 Canny算子弱邊緣檢測結果Fig.4 Weak edge detection result by using Canny operator

圖5 Canny 算子強邊緣檢測結果Fig.5 Strong edge detection result by using Canny operator
最后以高閾值提取的強邊緣作為“種子”,在低閾值提取的弱邊緣中尋找其 8連通域,作為Canny算子檢測的最終邊緣[4]。Canny算子檢測的尾流最終邊緣如圖6所示。
需要說明的是,設定的高閾值不需要檢測出所有的尾流邊緣,只需要在弱邊緣的基礎上有分布,之后會利用強、弱邊緣,進一步判斷、連接尾流邊緣。從圖6中可以看出,Canny算子處理得到的尾流邊緣圖像噪聲較小。

圖6 Canny 算子最終的邊緣檢測結果Fig.6 Final edge detection result by using Canny operator
拉普拉斯(Laplace)算子采用的是二階差分計算,在數字圖像中(如尾流圖像),對于一階微分圖,認為極大值或者極小值點,是邊緣點;對于二階微分圖,認為極大值和極小值間過零的點是邊緣點[5]。
Laplace算子在點(x,y)處的定義為

用差分方程近似二階偏導,結果如下:

提取系數后得Laplace算子模板為

Laplace算子二階差分的性質對噪聲會有雙倍加強,所以在進行尾流邊緣提取時,將高斯濾波與算子結合起來使用。Laplace算子模板即卷積核,因為尾流圖像是二維圖像,需要在兩個方向求導,直接使用Laplace算子即可自動完成,而上文中使用Sobel一階算子時需要分兩個方向計算。將卷積核與尾流圖像進行卷積運算,當Laplace算子檢測出過零點時,判定有尾流邊緣存在,但不包括無意義的過零點,因為當二階導數為0時,不僅僅是在尾流邊緣,也有可能在獨立的噪聲點上。利用Laplace算子檢測尾流邊緣的結果如圖7所示。
在圖像處理中,邊緣檢測是重要的一步。尤其是在尾流測量中不僅要求較好地檢測到邊緣,還要求邊緣盡可能連續。在圖7中,利用Laplace算子檢測尾流圖像邊緣時可以檢測到較為全面的尾流邊緣信息,但是同時也有更多的噪聲。由于噪聲的影響,有些特征的邊界不清晰。同樣,經典的邊緣檢測算子由于引入了各種形式的微分運算,對噪聲非常敏感,且不能得到較好的結果,不利于后續工作,因此Laplace并不是理想的檢測尾流的算法。

圖7 拉普拉斯算子的邊緣檢測結果Fig.7 Edge detection result by using Laplace operator
利用Matlab軟件對尾流圖像進行形態學處理,可以把圖像中的尾流邊緣分割出來。腐蝕過程注重物體的邊界點,可以消除或者弱化邊界點,通過腐蝕操作可以對不同區域進行分割并且可以消除目標結構元素的點。與前文介紹過的微分算子相似,在形態學算法中,結構元素也能用模板形象定義。矩形模板中,1表示屬于結構元素,0表示不屬于結構元素。遍歷原尾流圖像中的每一個像素點,將其與模板的原點(原點需要自己定義,在本文中定義中心點為原點)對齊,若模板中心為1的位置在原圖像上都有對應的像素,則保留模板所在位置的像素點,然后取模板中所有1覆蓋下、原圖中對應像素中最小值,用這個最小值替換當前像素值。
相同圖像利用不同結構元素進行腐蝕操作的結果不同,因此結構元素的形狀和大小決定了腐蝕操作的性能。圖8給出了本文模板和它對應的結構元素形狀。

圖8 圖像腐蝕模板Fig.8 Image corrosion template
膨脹操作與腐蝕操作相對。在遍歷圖像的像素點的過程中,保留模板所在位置的像素點,然后取模板中所有 1覆蓋下、原圖中對應像素中的最大值,用這個最大值替換當前像素值。腐蝕操作使目標圖像收縮,膨脹操作使目標圖像擴展。利用收縮或擴展后的圖像可以精確地提取原圖像的內外邊界。這樣提取的邊界理論上更精確,也能有效控制邊界線條的寬度[6]。
先對圖像進行腐蝕操作然后再進行膨脹操作的運算稱為開啟運算。對尾流圖像進行開啟運算后結果如圖9所示。
從圖9中可以看出,與尾流原始圖相比,尾流背景中的某些數據元素被消除,集中體現在-30°和30°附近,主要因為這些數據點的回波強度相對較弱。在圖9中也發現,尾流圖像內部出現了一個個閉合的小圈,形態學處理算法會把相似的鄰近的物體連接起來,這在某些圖像處理中有非常好的效果,但在尾流圖像處理中我們需要的是尾流的外邊界,至于尾流內部回波強度的大小,只需要大致范圍即可,并不需要細致的歸類。要想得到尾流完整的外邊界還需要對圖 9中的聲吶尾流數據設置閾值,進一步提取邊界,這就增加了算法的復雜性。

圖9 先腐蝕再膨脹的邊緣檢測結果Fig.9 Edge detection results of first corrosion and then expansion
先對圖像進行膨脹操作然后再進行腐蝕操作的運算稱為閉合運算。對圖像進行閉合運算后,尾流圖像進行閉運算結果如圖10所示。
由圖10可以看出,利用閉合算法可以把尾流邊緣大致勾勒出來,但是尾流區域內部也出現了更小的封閉邊緣。閉合算法與開啟算法相比,尾流元素丟失更多,在開角 -20°和20°附近也丟失了很多數據。

圖10 先膨脹再腐蝕的邊緣檢測結果Fig.10 Edge detection results of first expansion and then corrosion
分形維數一般都是通過線性擬合的方法計算。最常用的分形維數計算方法是盒計數分形算法,在對尾流邊緣提取的過程中,采用的是差分盒計數分形算法,其他的計算方法大部分都是基于盒計數的基礎上[7]。下面介紹差分盒計數分形算法的原理和實現過程。
盒計數分形維數通過式(4)計算:

其中:Nr為盒子個數;r為盒子大小,通常取r=1,2,4··2i。
把圖像放置在平面上,構成的灰度曲面如圖11所示,然后把圖 11中的平面用網格進行劃分,然后用不同大小的盒子去覆蓋圖像的灰度平面。計算覆蓋整個灰度曲面所需的盒子的個數,上述過程的示意圖11所示。

圖11 差分盒子算法示意圖Fig.11 Schematic diagram of differential box algorithm
對于不同的r可得到不同的Nr:

式中:nr(i,j)為小盒子尺度r下的像素點(i,j)處的小盒子數量。Nr這是尺度為r時覆蓋整幅圖像的盒子數,變換尺度又得到另一組數據。采用最小二乘法擬合得到的各組點(l og2r,log2Nr),并在雙對數坐標系中標出,通過對這些分布點進行最小二乘線性擬合,直線的斜率取絕對值后即為盒維數[8]。
但采用盒維數計算尾流圖像的分形維數時,首先需將圖像中所關心的尾流區域提取出來,先轉化為黑白位圖,在相應矩陣是由 1、0表示,離散的像素點很容易判斷網格是否覆蓋了圖像中的尾流區域。統計出相應網格中覆蓋尾流區域的網格數目,擬合8個點為例,即在8個尺度下,得出的尾流圖像維數為1.77,誤差e為1.5%,結果如圖12所示。

圖12 最小二乘法擬合的盒維數圖(維數為1.77,誤差為1.5%)Fig.12 Box dimension diagram of least squares fitting(Dimension is 1.77,Error is 1.5%)
圖12顯示了 lo g2(Nr)與log2r的典型關系。擬合誤差e值越低,說明擬合誤差越小,擬合效果越好??梢愿鶕S數d的數值變化和波動確定圖像的邊緣并作為分割標準來分割圖像。通過設置閾值,如果d在閾值內則判定此區域為尾流內部,若超出閾值則判定為邊界。本文為了提取出尾流邊緣,經過多次試驗將d的范圍設定為1.7~1.9。
擬合 18個點得到的尾流圖像維數為 1.85,誤差e為3%,如圖13所示。

圖13 最小二乘法擬合的盒維數圖(維數為1.85,誤差為3%)Fig.13 Box dimension diagram of least squares fitting(Dimension is 1.85,Error is 3%)
當輸入圖像(圖像必須為二值化圖像或者是灰度圖像)時,函數的返回值為一個與輸入圖像大小一致的二值化圖像。函數會根據設置的尾流維數閾值檢測輸入圖像的邊緣,檢測到邊緣的圖像位置數值為1,其余部分用0填充。這樣在讀取尾流圖像時,分形算法會對圖像邊緣做處理,勾出尾流的邊界。經分形處理后的效果圖如圖14所示。

圖14 分形算法邊緣檢測結果Fig.14 Edge detection results of fractal algorithm
由圖 14中可以清晰地看到尾流圖像的上下邊緣,尾流邊緣被分割出來。在計算盒子的過程中,隨著小盒子大小的變化,盒子數量會發生變化,導致維數d也會變化,對不同圖像的影響程度也不一樣,在尾流圖像中需要對不同大小的盒子尺度進行計算從而得到維數范圍,使用計算簡單的盒維數可以快速求出圖像的分形維數。與傳統的經典圖像邊界算法相比,分形盒計數算法能夠忽略紋理內部的變化情況,能夠很好地提取尾流圖像的邊緣輪廓,這為下一步尾流特性的研究奠定了基礎。該算法適應于維數變化明顯的邊緣提取過程,但是盒維數不能很好地反映出圖像中的具體變化情況,其值僅由圖像塊中灰度的最大值和最小值決定。
通過幾種邊緣檢測算法對尾流圖像處理的結果進行分析,得到以下結論:(1)Canny微分算子在算法實現過程中首先進行高斯濾波,然后使用一階微分算子計算梯度幅值與方向,經過非極大值抑制后,通過雙閾值判定、連接尾流邊緣,低閾值檢測的邊緣可以提取豐富的尾流信號信息,高閾值檢測的邊緣噪聲較小,不足之處是算法較為復雜過程繁瑣,也會出現邊緣不連續的情況,從尾流圖像的處理結果來看是較為理想的邊緣提取算法。(2)Laplace微分檢測算子在對尾流圖像邊緣提取的時候可以檢測到尾流的邊緣信息,但鑒于此種算法基于的是二階導數的原理,一階導數的局部峰值會造成Laplace算子判定為過零點,容易被噪聲影響,所以理論上不適用于尾流圖像處理。(3)數學形態學算法不像微分算法那樣對噪聲敏感,同時也能提取到光滑直觀的尾流圖像邊緣,抑制了小的離散點或者尖峰,也能夠連接鄰近的物體圖像填充細小空洞,但這也造成了在尾流圖像內部出現了不同程度的閉合區域,出現過分割的現象。(4)分形算法利用不同物體具有不同的分形維數這一原則,通過判斷維數奇異值來判斷尾流圖像與周圍區域的分界處,從而提取到尾流邊緣,不牽扯到微分計算過程,對噪聲敏感性弱,但該算法更適應于維數變化明顯的邊緣提取過程,隨著分割精度的增加,要求的迭代次數也大大增加,造成運行時間比其他算法長。