一種寬帶聲場干涉圖案處理的方法及其應(yīng)用

任群言，樸勝春，HERMAND Jean-Pierre

(1.哈爾濱工程大學(xué)水聲技術(shù)重點(diǎn)實(shí)驗(yàn)室，黑龍江哈爾濱150001;2.布魯塞爾自由大學(xué) 水下聲學(xué)環(huán)境實(shí)驗(yàn)室，比利時(shí) 布魯塞爾1050)

由移動(dòng)聲源(水面船只)產(chǎn)生的寬帶聲場可呈現(xiàn)出有規(guī)律的干涉條紋.俄羅斯學(xué)者最先提出了波導(dǎo)不變量理論［1］解釋和表征此干涉現(xiàn)象，其中波導(dǎo)不變量經(jīng)常被記作β，其值和波導(dǎo)特征緊密相關(guān).波導(dǎo)不變量理論現(xiàn)已被廣泛應(yīng)用于海底沉積地聲參數(shù)反演［2］，被動(dòng)聲源定位［3-5］，目標(biāo)識(shí)別［6］以及波導(dǎo)頻散效應(yīng)消除［7］.在上述應(yīng)用中，對(duì)干涉圖案的處理和特征提取(條紋斜率)是極其重要的步驟.Radon變換［4］，Hough 變換［5-8］以及二維傅里葉變換(2DFFT)［6］已被用于處理干涉圖案來估計(jì)干涉條紋的斜率.國內(nèi)也有學(xué)者用頻移補(bǔ)償方法計(jì)算寬帶聲場的波導(dǎo)不變量［9］.Brooks等［10］對(duì)一些條紋斜率提取方法的優(yōu)劣性做了比較.上述的常用處理方法一般只能用來估計(jì)干涉圖像中由波導(dǎo)決定的條紋整體斜率.由于以上處理方法的局限性，與波導(dǎo)特性密切相關(guān)的干涉條紋位置信息［11］一直沒有被利用.因此，同時(shí)從干涉圖案中提取干涉條紋的斜率和位置信息，無論對(duì)提高基于干涉圖案應(yīng)用的準(zhǔn)確度，還是對(duì)拓展更多的應(yīng)用領(lǐng)域都是很有意義的.

本文引入一種基于圖像Hessian矩陣特征值分析的多尺度線性濾波器［12］，研究了有效地檢測和分離圖案中不同寬度條紋的方法，具體條紋的斜率和位置信息可以很方便地從處理后的結(jié)果中提取出來.通過仿真數(shù)據(jù)討論了典型淺海環(huán)境下，低頻條紋隨沉積層地聲參數(shù)的變化規(guī)律及其在沉積層地聲參數(shù)反演研究的應(yīng)用.最后通過對(duì)海上實(shí)測水面船只寬帶輻射噪聲數(shù)據(jù)的處理，檢驗(yàn)了該方法對(duì)干涉聲場條紋結(jié)構(gòu)提取和特征參數(shù)估計(jì)的性能.

1 寬帶聲場的干涉現(xiàn)象

考慮水平分層波導(dǎo)中深度z0處無指向性點(diǎn)聲源的輻射聲場，在距離聲源水平距離r、深度z處水聽器接收的復(fù)聲場p可以表示為一系列簡正波的相干疊加［13］:

圖1是基于典型淺海環(huán)境模型［14］由簡正波程序Kraken［15］計(jì)算的寬帶標(biāo)量聲強(qiáng)，聲源和接收深度分別被設(shè)為3.5 m和20 m.可以認(rèn)為圖1中某一特定的條紋的強(qiáng)度在某一距離和頻帶內(nèi)不隨距離和頻率變化而變化，即

從上式中可以推導(dǎo)出波導(dǎo)不變量［13］:

式中:u和v分別表示一定頻帶范圍內(nèi)的群速度和相速度的平均值，其大小均由波導(dǎo)特性質(zhì)決定.從該式也可以看出，波導(dǎo)不變量的值取決于該條紋的具體位置(距離 r，角頻率 ω)以及條紋的斜率(dω/d r).即波導(dǎo)不變量實(shí)際上是一個(gè)隨距離和聲波頻率變化的量，其在某種程度上更應(yīng)該被當(dāng)作一個(gè)分布，而不是一個(gè)單一的值.Heaney［2］就利用波導(dǎo)不變量的分布特性，進(jìn)行了海底沉積層特征的快速分類研究.

波導(dǎo)不變量主要用來表征由環(huán)境決定的條紋斜率信息.理論和實(shí)驗(yàn)研究表明干涉條紋的位置信息也由環(huán)境特性決定，該研究結(jié)果已被用于海洋聲層析［11］.為了盡可能地挖掘干涉圖案中條紋斜率和位置所蘊(yùn)含的波導(dǎo)信息，以用于環(huán)境參數(shù)的反演研究，本文在下節(jié)研究了將一種多尺度線性濾波法用于干涉圖像處理和條紋特征的提取的方法.

圖1 基于Yellow Shark環(huán)境模型仿真計(jì)算的聲強(qiáng)的頻率－距離分布Fig.1 Predicted sound intensity frequency-space distribution for the Yellow Shark environment

2 多尺度線條濾波器

對(duì)一個(gè)圖像M，其在點(diǎn)x0處進(jìn)行二階泰勒展開式可以寫為

式中:x=x0+ δx0，▽0，σ和 H0，σ分別是該圖像在點(diǎn) x0處以尺度σ計(jì)算得到的梯度和Hessian矩陣.圖像M的微分通常由與高斯函數(shù)的微分卷積得到

二維高斯函數(shù)經(jīng)常寫成如下形式:

參數(shù)γ用來比較濾波器在各尺度時(shí)的響應(yīng)［16］，當(dāng)進(jìn)行單尺度分析時(shí)，其值就簡單設(shè)置為1.由于該卷積過程相當(dāng)于對(duì)圖像進(jìn)行一次高斯濾波，故在高斯白噪聲干擾下，本濾波器具有很強(qiáng)的噪聲抑制能力.

Hessian矩陣將2D圖像里的分析區(qū)域映射為一個(gè)橢圓，并測量和分析區(qū)域內(nèi)外的對(duì)比度.設(shè)λσ，k、u^σ，k分別表示圖像在某一分析尺度下Hessian矩陣的特征值和特征向量，并認(rèn)為|λ1|≤|λ2|.則此時(shí)u^1就表示直線的方向(聲強(qiáng)變化最小的方向)相應(yīng)的就是與直線正交的方向.表1給出了圖像Hessian矩陣的特征值與被分析圖像局部結(jié)構(gòu)之間的關(guān)系［12］.

表1 二維圖像中可能的結(jié)構(gòu)，取決于λk的符號(hào)和大小Table 1 Possible patterns in 2D im age depending on the values of the eigenvaluesλk

基于以上特征值分析，該線性濾波器在某一尺度的濾波響應(yīng)由以下2個(gè)信息來決定.其一是橢圓的偏心率R:

式中:Soval表示橢圓面積，rmax表示最大半徑.R確保該濾波器只考慮圖形中的幾何信息.另外一個(gè)信息就是信噪比S，用來測量分析區(qū)域和背景的對(duì)比度，由Frobenius矩陣的范數(shù)表示:

具體地，如果一個(gè)點(diǎn)可以被歸于某線性結(jié)構(gòu)中，那么λ1就會(huì)很小(理想情況下為0)，而λ2則具有較大的值.相應(yīng)地，R值比較小而S值就比較大.結(jié)合上面2個(gè)測量，多尺度線性濾波器在某一尺度的響應(yīng)就可以寫成:

式中:β和c是用于控制R和S的靈敏度，應(yīng)該根據(jù)不同情況選擇不同值，在對(duì)仿真數(shù)據(jù)進(jìn)行處理分析時(shí)，β、c被分別設(shè)為1.0和15.上述方程是用來測量圖案中的暗條紋，對(duì)于明條紋測量，只要將上述條件反轉(zhuǎn)即可.利用式(10)計(jì)算圖像在各個(gè)尺度σ的濾波響應(yīng)后，最后輸出最接近圖像中直線寬度時(shí)尺度的濾波結(jié)果(最大響應(yīng)):

式中:σmin和σmax分別是最小和最大尺度，濾波器會(huì)在最小、最大尺度之間搜尋符合條件的線性結(jié)構(gòu).關(guān)于尺度的選擇問題，現(xiàn)在還沒有統(tǒng)一的結(jié)論，要根據(jù)不同情況下條紋寬度的不同來進(jìn)行選擇.在仿真數(shù)據(jù)分析中，它們分別被設(shè)為1和3，以0.5遞增.圖2給出了該濾波器應(yīng)用于圖1的結(jié)果.本文只給出了濾波器在整數(shù)尺度的輸出結(jié)果.從圖2(b)～(d)中可以看出，不同尺度下均有線性結(jié)構(gòu)出現(xiàn)，并且當(dāng)尺度增大時(shí)，濾波后會(huì)出現(xiàn)更多的線性結(jié)構(gòu).圖2(a)是最后多尺度濾波的結(jié)果，綜合了圖2(b)～(d)的結(jié)果，給出了圖像在給定尺度區(qū)間內(nèi)的線性結(jié)構(gòu)分布.

由圖2可以看出，該多尺度線性濾波器可以很有效地提取干涉圖案的干涉結(jié)構(gòu)，檢測和分離干涉圖案中不同寬度的條紋.該濾波器的優(yōu)點(diǎn)就在于只要給定搜索的條紋尺度范圍，此尺度范圍內(nèi)的條紋都可以被檢測出來.當(dāng)干涉條紋被成功地檢測和分離開來后，蘊(yùn)含在干涉條紋位置和斜率信息中的波導(dǎo)效應(yīng)就有可能提取出來.

圖2 圖像1經(jīng)多尺度線性濾波器在尺度區(qū)間［1，3］時(shí)的輸出及其在整數(shù)尺度的響應(yīng)Fig.2 The output of the multi-line filter for the scale interval from 1 to 3 with an increment of 0.5 applied to the sound intensity in Fig.1，and its selected outputs for fixed scales

3 利用條紋位置信息的沉積層地聲參數(shù)反演

仿真的環(huán)境模型以前面介紹的典型淺海環(huán)境下的海底模型為基礎(chǔ).從其他仿真研究［16］可知，在該環(huán)境模型下，沉積層的厚度和聲速對(duì)條紋位置具有決定性的作用，而沉積層吸收系數(shù)和密度對(duì)條紋位置影響不大，特別是對(duì)低頻部分的條紋.在這里我們討論利用條紋信息來獲取沉積層聲速和厚度的可能性.為了觀察當(dāng)沉積層厚度和聲速變化時(shí)，條紋位置的變化情況，研究中選取了一個(gè)具有0.5 m厚度，1 460 m/s聲速的沉積層作為參考沉積層(其他參數(shù)和文獻(xiàn)［14］保持一致).

圖3(a)是在選取參考沉積層的波導(dǎo)環(huán)境下仿真干涉圖像(濾波后)，從中選取了在低頻的3個(gè)條紋來觀察環(huán)境參數(shù)變化對(duì)它們的影響.由于3個(gè)條紋的斜率在沉積層參數(shù)變化時(shí)幾乎不發(fā)生變化，研究中只考慮沉積層參數(shù)變化對(duì)條紋水平位置(頻率)的影響.圖3(b)～(d)是這3個(gè)條紋位置與沉積層參數(shù)的關(guān)系圖，沉積層厚度從0.5 m變化到8 m，聲速從1 460 m/s變化到1 520 m/s，分別以0.5 m和10 m/s遞增.

從圖3中可以看出，這3個(gè)條紋隨沉積層參數(shù)的變化大體一致.對(duì)固定的沉積層厚度，沉積層聲速增加時(shí)，條紋向高頻方向移動(dòng);對(duì)固定的沉積層聲速，沉積層厚度增加時(shí)，條紋向低頻方向移動(dòng).由于這些條紋的頻率不同，對(duì)環(huán)境的敏感程度也不大形同，所以它們的位置信息能夠提供不同的沉積層信息，可以用來聯(lián)合反演沉積層的地聲參數(shù).作為一個(gè)簡單的例子，我們假定一個(gè)沉積層的厚度為4.5 m，聲速為1 490 m/s，則相對(duì)于參考沉積層，這3個(gè)條紋位置和相應(yīng)可能的沉積層參數(shù)如表2所示.從表中可以看出，它們的交集就是所設(shè)定的沉積層參數(shù)(4.5 m，1 490 m/s).如果它們的交集不止一個(gè)解時(shí)，還可以通過比較不同解時(shí)的干涉圖像和數(shù)據(jù)的相似度進(jìn)行比較篩選.這里給出的只是一個(gè)簡單的數(shù)值仿真，但說明了利用條紋位置信息進(jìn)行沉積層參數(shù)反演的可行性.

圖3 參考環(huán)境下的干涉圖像以及選定的觀察條紋和當(dāng)沉積層參數(shù)變化時(shí)條紋位置的位置函數(shù)Fig.3 The interference structure of the reference sediment and selected striations，and their ambiguity functions of the striation position with the sediment variations

表2 利用一組條紋位置來反演沉積層參數(shù)Table 2 Using a set of striations positions to invert for the sediment geoacoustic parameters

4 實(shí)際海上實(shí)驗(yàn)數(shù)據(jù)處理結(jié)果

海上實(shí)驗(yàn)在2008年夏季進(jìn)行，實(shí)驗(yàn)海域海深約25m.一個(gè)輻射寬帶噪聲的漁船被用作為目標(biāo)船，一個(gè)垂直接收陣吊放在另外一個(gè)自由漂浮的小船上，接收目標(biāo)船運(yùn)動(dòng)中輻射的寬帶噪聲.目標(biāo)船先向接收船靠近，然后遠(yuǎn)離.實(shí)驗(yàn)約采集10 min的目標(biāo)船輻射噪聲信號(hào).

圖4是海上實(shí)驗(yàn)中深度6 m處水聽器接收目標(biāo)船寬度噪聲信號(hào)的時(shí)頻干涉圖.從圖中可以看到有規(guī)律的干涉條紋，也可以看到與水面船只的螺旋槳轉(zhuǎn)動(dòng)以及發(fā)動(dòng)機(jī)振動(dòng)有關(guān)的線譜.圖5給出了應(yīng)用多尺度濾波器處理圖4中時(shí)頻干涉圖的結(jié)果.圖5和圖4對(duì)比可以發(fā)現(xiàn)，不同寬度的干涉條紋被提取和分離開來，進(jìn)一步驗(yàn)證了該條紋提取方法在實(shí)際數(shù)據(jù)處理中的有效性.從圖5中可以看到艦船噪聲的干涉圖案相對(duì)于兩船最近位置具有極高的對(duì)稱性，說明該測試海域的海底聲學(xué)特性可以近似為和距離無關(guān).這為下一步的海底地聲模型建立和地聲參數(shù)反演提供了一定的先驗(yàn)知識(shí).

圖4 海上實(shí)驗(yàn)測量的艦船噪聲時(shí)頻干涉Fig.4 The measured acoustic pressure spectrogram of themoving surface ship

圖5 由多尺度線性濾波器得到的移動(dòng)船只寬帶噪聲場的條紋結(jié)構(gòu)Fig.5 The filtered line structure of the measured acoustic pressure spectrogram

5 結(jié)束語

本文研究了利用一種多尺度線性濾波器提取寬帶聲場干涉結(jié)構(gòu)和條紋信息的有效方法.通過對(duì)仿真數(shù)據(jù)和實(shí)驗(yàn)數(shù)據(jù)的處理結(jié)果可以看出，該方法可以有效地檢測和提取寬帶聲場中不同寬度的條紋結(jié)構(gòu)，尤其對(duì)干涉結(jié)構(gòu)的低頻部分.考慮到聲場干涉條紋的位置反映了波導(dǎo)的環(huán)境參數(shù)，本文也通過仿真計(jì)算初步研究了利用條紋位置來進(jìn)行沉積層地聲參數(shù)反演的可行性.仿真結(jié)果表明，利用提取的干涉條紋斜率和位置信息，可建立一個(gè)簡單有效的沉積層地聲參數(shù)估計(jì)方法.

［1］CHUPROV S.Interference structure of a sound field in a layered ocean［J］.Ocean Acoustic，Current State，1982:71-91.

［2］HEANEY K.Rapid geoacoustic characterization using a surface ship of opportunity［J］.IEEE Journal of Oceanic Engineering，2004，29(1):88-99.

［3］WESTWOOD E.Broadband matched-field source localization［J］.The Journal of the Acoustical Society of America，1992，91:2777-2789.

［4］TURGUT A，ORR M，ROUSEFF D.Broadband source localization using horizontal-beam acoustic intensity striations［J］.The Journal of the Acoustical Society of America，2010，127:73-83.

［5］D＇SPAIN G，KUPERMANW.Application of waveguide invariants to analysis of spectrograms from shallow water environment that vary in range and azimuth［J］.The Journal of the Acoustical Society of America，1999，106:2454-2468.

［6］AN L，WANG Z，LU J.Calculating the waveguide invariant by passive sonar LOFAR gram image［C］ //Proceedings of 14th International Conference on Mechatronics and Machine Vision in Practice.Xiamen，China，2007.

［6］王寧，高大治，王好忠.頻散，聲場干涉結(jié)構(gòu)，波導(dǎo)不變量與消頻散變換［J］.哈爾濱工程大學(xué)學(xué)報(bào)，2010，31(7):825-831.WANG Ning，GAO Dazhi，WANG Haozhong.Mode dispersion，interference，dispersion，waveguide invariance and the dispersion transform［J］.Journal of Harbin Engineering University，2010，31(7):825-831.

［7］田玲愛，劉福臣，周士弘.利用Hough變換提取波導(dǎo)不變量［J］.聲學(xué)與電子工程，2009，4:22-24.TIAN Ling＇ai，LIU Fuchen，ZHOU Shihong.Extraction of waveguide invariant using Hough transform［J］.Acoustics and Electronics Engineering，2009，4:22-24.

［8］蘇曉星，張仁和，李風(fēng)華.用頻移補(bǔ)償方法計(jì)算聲場的波導(dǎo)不變量［J］.聲學(xué)技術(shù)，2007，26(6):1073-1076.SU Xiaoxing，ZHANG Renhe，LI Fenghua.Calculation of the waveguide invariant of acoustic fields by using frequency-shift compensation method［J］.Technical Acoustics，2007，26(6):1073-1076.

［9］BROOKS L，KIDNERM，ZANDER A，etal.Striation processing of spectrogram data［C］ //Proceedings of the 13th International Congress on Sound and Vibration(ICSV13).Vienna，Austria，2006.

［10］KUZ＇K V M，PERESELKOV SA，PETNIKOV E A.The possibility of reconstruction of two-dimensional random inhomogeneities in a shallow sea by frequency shifts of the spatial interference structure of the sound field［J］.Physics of Wave Phenomena，2008，16(1):42-51.

［11］FRANGI A F，NIESSEN W J，VINCKEN K L，et al.Multiscale vessel enhancement filtering［C］//Medical Image Computing and Computer-Assisted Intervention(MICCAT＇98).Cambridge MA，USA，1998.

［12］BREKHOVSKIKH L，LYSANOV Y.Fundamentals of ocean acoustic［M］.3rd ed.New York:Springer Verlag，2003:143-148.

［13］HERMAND JP.Broad-band geoacoustic inversion in shallow water from waveguide impulse response measurements on a single hydrophone:theory and experiment result［J］.IEEE Journal of Oceanic Engineering，1999，24(1):41-66.

［14］PORTERM.The Kraken normal mode program［EB/OL］.［2011-01-15］.Available:http://oalib.hlsresearch.com/Modes/kraken.pdf.1997.

［15］LINDEBERG T.Edge detection and ridge detection with automatic scale selection［C］ //Proc Conf on Comp Vis and Pat Recog.San Francisco，1996.

［16］REN Q Y，PIAO SC，HERMAND JP.Space-frequency distribution of the vector field of broad-band sound in shallow water［C］ //Oceans＇10 MTS/IEEE Seattle Conference-Innerspace:A Global Responsibility.Seattle，2010.