一種寬帶聲場干涉圖案處理的方法及其應用

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

(1.哈爾濱工程大學水聲技術重點實驗室，黑龍江哈爾濱150001;2.布魯塞爾自由大學 水下聲學環境實驗室，比利時 布魯塞爾1050)

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

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

1 寬帶聲場的干涉現象

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

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

從上式中可以推導出波導不變量［13］:

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

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

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

2 多尺度線條濾波器

對一個圖像M，其在點x0處進行二階泰勒展開式可以寫為

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

二維高斯函數經常寫成如下形式:

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

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

表1 二維圖像中可能的結構，取決于λk的符號和大小Table 1 Possible patterns in 2D im age depending on the values of the eigenvaluesλk

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

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

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

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

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

由圖2可以看出，該多尺度線性濾波器可以很有效地提取干涉圖案的干涉結構，檢測和分離干涉圖案中不同寬度的條紋.該濾波器的優點就在于只要給定搜索的條紋尺度范圍，此尺度范圍內的條紋都可以被檢測出來.當干涉條紋被成功地檢測和分離開來后，蘊含在干涉條紋位置和斜率信息中的波導效應就有可能提取出來.

圖2 圖像1經多尺度線性濾波器在尺度區間［1，3］時的輸出及其在整數尺度的響應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 利用條紋位置信息的沉積層地聲參數反演

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

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

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

圖3 參考環境下的干涉圖像以及選定的觀察條紋和當沉積層參數變化時條紋位置的位置函數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 利用一組條紋位置來反演沉積層參數Table 2 Using a set of striations positions to invert for the sediment geoacoustic parameters

4 實際海上實驗數據處理結果

海上實驗在2008年夏季進行，實驗海域海深約25m.一個輻射寬帶噪聲的漁船被用作為目標船，一個垂直接收陣吊放在另外一個自由漂浮的小船上，接收目標船運動中輻射的寬帶噪聲.目標船先向接收船靠近，然后遠離.實驗約采集10 min的目標船輻射噪聲信號.

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

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

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

5 結束語

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

［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］王寧，高大治，王好忠.頻散，聲場干涉結構，波導不變量與消頻散變換［J］.哈爾濱工程大學學報，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變換提取波導不變量［J］.聲學與電子工程，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］蘇曉星，張仁和，李風華.用頻移補償方法計算聲場的波導不變量［J］.聲學技術，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.