任秀艷,陳夢英,許偉杰
(1.中國科學(xué)院聲學(xué)研究所東海研究站,上海201815;2.中國科學(xué)院大學(xué),北京100190)
尾流作為一類水下聲目標(biāo),具有十分重要的軍事意義,在海戰(zhàn)之中起著重要的作用:(1)尾流能夠影響探測設(shè)備的正常工作;(2)尾流是魚雷進(jìn)行探測、攻擊目標(biāo)的重要物理場;(3)尾流能作為潛艇規(guī)避的屏障[1]。因此,對各類艦船尾流的聲學(xué)特性,包括尾流的產(chǎn)生、持續(xù)時(shí)間和空間的分布、尾流中氣泡升浮的速度以及氣泡散射與吸收的特征、水下聲學(xué)尾流成像等方面進(jìn)行研究日益重要[2]。
既準(zhǔn)確又直觀的尾流圖像可以為下一步的尾流特征研究奠定基礎(chǔ),而在尾流多波束聲吶的原始圖像中,不僅有尾流回波,還有海面回波反射以及其他噪聲干擾,因此需要采取有效的圖像分割方法把尾流從中提取出來,從而分析尾流特征。對尾流圖像處理時(shí),可把尾流圖像看作三維表面,即x、y軸分別為圖像的橫、縱向矩陣列數(shù),z代表強(qiáng)度的灰度。圖像邊緣分割已有較長的歷史,邊緣是圖像特性發(fā)生變化的位置,不同的圖像灰度不同,邊界處會有明顯的邊緣出現(xiàn),利用這種特征可以進(jìn)行圖像的分割。
多波束聲吶測量系統(tǒng)由船載分系統(tǒng)和水下分系統(tǒng)兩大部分組成,水下系統(tǒng)的水密艙中包含多波束聲吶。課題組于海上進(jìn)行了多次實(shí)驗(yàn),采用的坐底式多波束測量系統(tǒng)圖如圖1所示。
圖1 測量系統(tǒng)示意圖Fig.1 Schematic diagram of the measuring system
原始尾流圖像如圖2所示。當(dāng)主瓣收到的尾流信號與旁瓣收到的海面反射信號同時(shí)到達(dá)時(shí),海面反射信號會干擾尾流信號。
圖2 原始尾流圖像Fig.2 Original wake image
圖2中,海平面在高度約為30 m的中間亮條處,海平面上方的虛擬強(qiáng)度圖像由多波束旁瓣引起,在此不做討論。根據(jù)分析,亮環(huán)即為海面回波反射干擾,中間條狀區(qū)域即為尾流圖像區(qū)域。尾流圖像邊緣提取的目的是將亮環(huán)剔除,并將條狀區(qū)域的尾流邊緣提取出來。
利用傳統(tǒng)的濾波方法難以去除聲吶尾流圖像中海面反射強(qiáng)度干擾,因此本文用背景相減法來去除海面回波干擾,構(gòu)造了一個(gè)包含“亮環(huán)”噪聲并且沒有尾流圖像的背景環(huán)境,通過聲吶圖像與構(gòu)造的“亮環(huán)圖像”相減得到的尾流圖像如圖3所示。
以圖3為原始圖像進(jìn)行尾流圖像邊緣分割,下面分別介紹幾種邊緣分割算法,并給出這幾種邊緣分割算法提取到的尾流邊緣。
圖3 背景相減法濾波結(jié)果Fig.3 Results of background subtraction filtering
基于一階的邊緣檢測算子中,Canny算子是公認(rèn)的具有良好檢測性能的一類邊緣檢測算子,因此本文利用它對尾流圖像進(jìn)行邊緣提取。
首先,用高斯濾波器對尾流圖像進(jìn)行平滑處理。然后,利用普通的一階邊緣檢測算子,如Sobel算子、Roberts算子等進(jìn)行梯度計(jì)算,本文中選用的Sobel算子。再利用非極大值抑制的方法判定邊緣梯度值,在尾流檢測的過程中,進(jìn)行的非極大值抑制是在梯度方向上的非極大值抑制,實(shí)際處理的尾流圖像的像素點(diǎn)是離散的二維矩陣,對于某個(gè)監(jiān)測點(diǎn),梯度方向兩側(cè)的點(diǎn)不一定存在于二維矩陣中,本文通過插值來解決這個(gè)問題。將某點(diǎn)梯度幅值與插值作比較,若此處梯度幅值大于兩個(gè)插值,則此點(diǎn)就是局部極大值,為尾流邊緣點(diǎn)。最后,進(jìn)行閾值設(shè)置、拼接邊緣,具體思路是選取兩個(gè)閾值,認(rèn)為小于低閾值的點(diǎn)是弱邊緣,大于高閾值的點(diǎn)是強(qiáng)邊緣,介于中間的像素點(diǎn)需要進(jìn)一步檢查[3]。
邊緣檢測過程中尾流圖像閾值的設(shè)置很關(guān)鍵,也是難點(diǎn),合適的閾值既要盡量保存尾流邊緣信息也要抑制噪聲,避免引入錯(cuò)誤邊緣。Matlab自帶的Canny算子就是通過直方圖來設(shè)定閾值的。它默認(rèn)選擇直方圖中幅值大小排在前 70%的最后一個(gè)閾值為高閾值,低閾值為高閾值的 0.2倍。本文在利用Canny算子處理實(shí)驗(yàn)數(shù)據(jù)時(shí)借鑒了此方法,先將所有梯度值進(jìn)行排序,經(jīng)過多次測試,高閾值從高到低開始試驗(yàn),保證噪聲逐漸增大,因?yàn)殚撝翟O(shè)定越高,噪聲抑制效果越好。最終選擇排在整個(gè)梯度值序列中第 98% 處的幅值作為高閾值,低閾值約為高閾值的0.2倍。這里設(shè)定的閾值是經(jīng)過多次實(shí)驗(yàn)后選取的。圖4、5為利用低閾值和高閾值檢測到的弱、強(qiáng)邊緣。
圖4 Canny算子弱邊緣檢測結(jié)果Fig.4 Weak edge detection result by using Canny operator
圖5 Canny 算子強(qiáng)邊緣檢測結(jié)果Fig.5 Strong edge detection result by using Canny operator
最后以高閾值提取的強(qiáng)邊緣作為“種子”,在低閾值提取的弱邊緣中尋找其 8連通域,作為Canny算子檢測的最終邊緣[4]。Canny算子檢測的尾流最終邊緣如圖6所示。
需要說明的是,設(shè)定的高閾值不需要檢測出所有的尾流邊緣,只需要在弱邊緣的基礎(chǔ)上有分布,之后會利用強(qiáng)、弱邊緣,進(jìn)一步判斷、連接尾流邊緣。從圖6中可以看出,Canny算子處理得到的尾流邊緣圖像噪聲較小。
圖6 Canny 算子最終的邊緣檢測結(jié)果Fig.6 Final edge detection result by using Canny operator
拉普拉斯(Laplace)算子采用的是二階差分計(jì)算,在數(shù)字圖像中(如尾流圖像),對于一階微分圖,認(rèn)為極大值或者極小值點(diǎn),是邊緣點(diǎn);對于二階微分圖,認(rèn)為極大值和極小值間過零的點(diǎn)是邊緣點(diǎn)[5]。
Laplace算子在點(diǎn)(x,y)處的定義為
用差分方程近似二階偏導(dǎo),結(jié)果如下:
提取系數(shù)后得Laplace算子模板為
Laplace算子二階差分的性質(zhì)對噪聲會有雙倍加強(qiáng),所以在進(jìn)行尾流邊緣提取時(shí),將高斯濾波與算子結(jié)合起來使用。Laplace算子模板即卷積核,因?yàn)槲擦鲌D像是二維圖像,需要在兩個(gè)方向求導(dǎo),直接使用Laplace算子即可自動完成,而上文中使用Sobel一階算子時(shí)需要分兩個(gè)方向計(jì)算。將卷積核與尾流圖像進(jìn)行卷積運(yùn)算,當(dāng)Laplace算子檢測出過零點(diǎn)時(shí),判定有尾流邊緣存在,但不包括無意義的過零點(diǎn),因?yàn)楫?dāng)二階導(dǎo)數(shù)為0時(shí),不僅僅是在尾流邊緣,也有可能在獨(dú)立的噪聲點(diǎn)上。利用Laplace算子檢測尾流邊緣的結(jié)果如圖7所示。
在圖像處理中,邊緣檢測是重要的一步。尤其是在尾流測量中不僅要求較好地檢測到邊緣,還要求邊緣盡可能連續(xù)。在圖7中,利用Laplace算子檢測尾流圖像邊緣時(shí)可以檢測到較為全面的尾流邊緣信息,但是同時(shí)也有更多的噪聲。由于噪聲的影響,有些特征的邊界不清晰。同樣,經(jīng)典的邊緣檢測算子由于引入了各種形式的微分運(yùn)算,對噪聲非常敏感,且不能得到較好的結(jié)果,不利于后續(xù)工作,因此Laplace并不是理想的檢測尾流的算法。
圖7 拉普拉斯算子的邊緣檢測結(jié)果Fig.7 Edge detection result by using Laplace operator
利用Matlab軟件對尾流圖像進(jìn)行形態(tài)學(xué)處理,可以把圖像中的尾流邊緣分割出來。腐蝕過程注重物體的邊界點(diǎn),可以消除或者弱化邊界點(diǎn),通過腐蝕操作可以對不同區(qū)域進(jìn)行分割并且可以消除目標(biāo)結(jié)構(gòu)元素的點(diǎn)。與前文介紹過的微分算子相似,在形態(tài)學(xué)算法中,結(jié)構(gòu)元素也能用模板形象定義。矩形模板中,1表示屬于結(jié)構(gòu)元素,0表示不屬于結(jié)構(gòu)元素。遍歷原尾流圖像中的每一個(gè)像素點(diǎn),將其與模板的原點(diǎn)(原點(diǎn)需要自己定義,在本文中定義中心點(diǎn)為原點(diǎn))對齊,若模板中心為1的位置在原圖像上都有對應(yīng)的像素,則保留模板所在位置的像素點(diǎn),然后取模板中所有1覆蓋下、原圖中對應(yīng)像素中最小值,用這個(gè)最小值替換當(dāng)前像素值。
相同圖像利用不同結(jié)構(gòu)元素進(jìn)行腐蝕操作的結(jié)果不同,因此結(jié)構(gòu)元素的形狀和大小決定了腐蝕操作的性能。圖8給出了本文模板和它對應(yīng)的結(jié)構(gòu)元素形狀。
圖8 圖像腐蝕模板Fig.8 Image corrosion template
膨脹操作與腐蝕操作相對。在遍歷圖像的像素點(diǎn)的過程中,保留模板所在位置的像素點(diǎn),然后取模板中所有 1覆蓋下、原圖中對應(yīng)像素中的最大值,用這個(gè)最大值替換當(dāng)前像素值。腐蝕操作使目標(biāo)圖像收縮,膨脹操作使目標(biāo)圖像擴(kuò)展。利用收縮或擴(kuò)展后的圖像可以精確地提取原圖像的內(nèi)外邊界。這樣提取的邊界理論上更精確,也能有效控制邊界線條的寬度[6]。
先對圖像進(jìn)行腐蝕操作然后再進(jìn)行膨脹操作的運(yùn)算稱為開啟運(yùn)算。對尾流圖像進(jìn)行開啟運(yùn)算后結(jié)果如圖9所示。
從圖9中可以看出,與尾流原始圖相比,尾流背景中的某些數(shù)據(jù)元素被消除,集中體現(xiàn)在-30°和30°附近,主要因?yàn)檫@些數(shù)據(jù)點(diǎn)的回波強(qiáng)度相對較弱。在圖9中也發(fā)現(xiàn),尾流圖像內(nèi)部出現(xiàn)了一個(gè)個(gè)閉合的小圈,形態(tài)學(xué)處理算法會把相似的鄰近的物體連接起來,這在某些圖像處理中有非常好的效果,但在尾流圖像處理中我們需要的是尾流的外邊界,至于尾流內(nèi)部回波強(qiáng)度的大小,只需要大致范圍即可,并不需要細(xì)致的歸類。要想得到尾流完整的外邊界還需要對圖 9中的聲吶尾流數(shù)據(jù)設(shè)置閾值,進(jìn)一步提取邊界,這就增加了算法的復(fù)雜性。
圖9 先腐蝕再膨脹的邊緣檢測結(jié)果Fig.9 Edge detection results of first corrosion and then expansion
先對圖像進(jìn)行膨脹操作然后再進(jìn)行腐蝕操作的運(yùn)算稱為閉合運(yùn)算。對圖像進(jìn)行閉合運(yùn)算后,尾流圖像進(jìn)行閉運(yùn)算結(jié)果如圖10所示。
由圖10可以看出,利用閉合算法可以把尾流邊緣大致勾勒出來,但是尾流區(qū)域內(nèi)部也出現(xiàn)了更小的封閉邊緣。閉合算法與開啟算法相比,尾流元素丟失更多,在開角 -20°和20°附近也丟失了很多數(shù)據(jù)。
圖10 先膨脹再腐蝕的邊緣檢測結(jié)果Fig.10 Edge detection results of first expansion and then corrosion
分形維數(shù)一般都是通過線性擬合的方法計(jì)算。最常用的分形維數(shù)計(jì)算方法是盒計(jì)數(shù)分形算法,在對尾流邊緣提取的過程中,采用的是差分盒計(jì)數(shù)分形算法,其他的計(jì)算方法大部分都是基于盒計(jì)數(shù)的基礎(chǔ)上[7]。下面介紹差分盒計(jì)數(shù)分形算法的原理和實(shí)現(xiàn)過程。
盒計(jì)數(shù)分形維數(shù)通過式(4)計(jì)算:
其中:Nr為盒子個(gè)數(shù);r為盒子大小,通常取r=1,2,4··2i。
把圖像放置在平面上,構(gòu)成的灰度曲面如圖11所示,然后把圖 11中的平面用網(wǎng)格進(jìn)行劃分,然后用不同大小的盒子去覆蓋圖像的灰度平面。計(jì)算覆蓋整個(gè)灰度曲面所需的盒子的個(gè)數(shù),上述過程的示意圖11所示。
圖11 差分盒子算法示意圖Fig.11 Schematic diagram of differential box algorithm
對于不同的r可得到不同的Nr:
式中:nr(i,j)為小盒子尺度r下的像素點(diǎn)(i,j)處的小盒子數(shù)量。Nr這是尺度為r時(shí)覆蓋整幅圖像的盒子數(shù),變換尺度又得到另一組數(shù)據(jù)。采用最小二乘法擬合得到的各組點(diǎn)(l og2r,log2Nr),并在雙對數(shù)坐標(biāo)系中標(biāo)出,通過對這些分布點(diǎn)進(jìn)行最小二乘線性擬合,直線的斜率取絕對值后即為盒維數(shù)[8]。
但采用盒維數(shù)計(jì)算尾流圖像的分形維數(shù)時(shí),首先需將圖像中所關(guān)心的尾流區(qū)域提取出來,先轉(zhuǎn)化為黑白位圖,在相應(yīng)矩陣是由 1、0表示,離散的像素點(diǎn)很容易判斷網(wǎng)格是否覆蓋了圖像中的尾流區(qū)域。統(tǒng)計(jì)出相應(yīng)網(wǎng)格中覆蓋尾流區(qū)域的網(wǎng)格數(shù)目,擬合8個(gè)點(diǎn)為例,即在8個(gè)尺度下,得出的尾流圖像維數(shù)為1.77,誤差e為1.5%,結(jié)果如圖12所示。
圖12 最小二乘法擬合的盒維數(shù)圖(維數(shù)為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的典型關(guān)系。擬合誤差e值越低,說明擬合誤差越小,擬合效果越好??梢愿鶕?jù)維數(shù)d的數(shù)值變化和波動確定圖像的邊緣并作為分割標(biāo)準(zhǔn)來分割圖像。通過設(shè)置閾值,如果d在閾值內(nèi)則判定此區(qū)域?yàn)槲擦鲀?nèi)部,若超出閾值則判定為邊界。本文為了提取出尾流邊緣,經(jīng)過多次試驗(yàn)將d的范圍設(shè)定為1.7~1.9。
擬合 18個(gè)點(diǎn)得到的尾流圖像維數(shù)為 1.85,誤差e為3%,如圖13所示。
圖13 最小二乘法擬合的盒維數(shù)圖(維數(shù)為1.85,誤差為3%)Fig.13 Box dimension diagram of least squares fitting(Dimension is 1.85,Error is 3%)
當(dāng)輸入圖像(圖像必須為二值化圖像或者是灰度圖像)時(shí),函數(shù)的返回值為一個(gè)與輸入圖像大小一致的二值化圖像。函數(shù)會根據(jù)設(shè)置的尾流維數(shù)閾值檢測輸入圖像的邊緣,檢測到邊緣的圖像位置數(shù)值為1,其余部分用0填充。這樣在讀取尾流圖像時(shí),分形算法會對圖像邊緣做處理,勾出尾流的邊界。經(jīng)分形處理后的效果圖如圖14所示。
圖14 分形算法邊緣檢測結(jié)果Fig.14 Edge detection results of fractal algorithm
由圖 14中可以清晰地看到尾流圖像的上下邊緣,尾流邊緣被分割出來。在計(jì)算盒子的過程中,隨著小盒子大小的變化,盒子數(shù)量會發(fā)生變化,導(dǎo)致維數(shù)d也會變化,對不同圖像的影響程度也不一樣,在尾流圖像中需要對不同大小的盒子尺度進(jìn)行計(jì)算從而得到維數(shù)范圍,使用計(jì)算簡單的盒維數(shù)可以快速求出圖像的分形維數(shù)。與傳統(tǒng)的經(jīng)典圖像邊界算法相比,分形盒計(jì)數(shù)算法能夠忽略紋理內(nèi)部的變化情況,能夠很好地提取尾流圖像的邊緣輪廓,這為下一步尾流特性的研究奠定了基礎(chǔ)。該算法適應(yīng)于維數(shù)變化明顯的邊緣提取過程,但是盒維數(shù)不能很好地反映出圖像中的具體變化情況,其值僅由圖像塊中灰度的最大值和最小值決定。
通過幾種邊緣檢測算法對尾流圖像處理的結(jié)果進(jìn)行分析,得到以下結(jié)論:(1)Canny微分算子在算法實(shí)現(xiàn)過程中首先進(jìn)行高斯濾波,然后使用一階微分算子計(jì)算梯度幅值與方向,經(jīng)過非極大值抑制后,通過雙閾值判定、連接尾流邊緣,低閾值檢測的邊緣可以提取豐富的尾流信號信息,高閾值檢測的邊緣噪聲較小,不足之處是算法較為復(fù)雜過程繁瑣,也會出現(xiàn)邊緣不連續(xù)的情況,從尾流圖像的處理結(jié)果來看是較為理想的邊緣提取算法。(2)Laplace微分檢測算子在對尾流圖像邊緣提取的時(shí)候可以檢測到尾流的邊緣信息,但鑒于此種算法基于的是二階導(dǎo)數(shù)的原理,一階導(dǎo)數(shù)的局部峰值會造成Laplace算子判定為過零點(diǎn),容易被噪聲影響,所以理論上不適用于尾流圖像處理。(3)數(shù)學(xué)形態(tài)學(xué)算法不像微分算法那樣對噪聲敏感,同時(shí)也能提取到光滑直觀的尾流圖像邊緣,抑制了小的離散點(diǎn)或者尖峰,也能夠連接鄰近的物體圖像填充細(xì)小空洞,但這也造成了在尾流圖像內(nèi)部出現(xiàn)了不同程度的閉合區(qū)域,出現(xiàn)過分割的現(xiàn)象。(4)分形算法利用不同物體具有不同的分形維數(shù)這一原則,通過判斷維數(shù)奇異值來判斷尾流圖像與周圍區(qū)域的分界處,從而提取到尾流邊緣,不牽扯到微分計(jì)算過程,對噪聲敏感性弱,但該算法更適應(yīng)于維數(shù)變化明顯的邊緣提取過程,隨著分割精度的增加,要求的迭代次數(shù)也大大增加,造成運(yùn)行時(shí)間比其他算法長。