宋文華王寧?高大治王好忠屈科

1)(中國海洋大學(xué)海洋技術(shù)系,青島 266100)

2)(廣東海洋大學(xué)電子與信息工程學(xué)院,湛江 524088)

波導(dǎo)不變量譜值及其分離方法?

宋文華1)王寧1)?高大治1)王好忠1)屈科2)

1)(中國海洋大學(xué)海洋技術(shù)系,青島 266100)

2)(廣東海洋大學(xué)電子與信息工程學(xué)院,湛江 524088)

(2016年12月30日收到;2017年3月19日收到修改稿)

波導(dǎo)不變量有兩種定義方式,按照干涉條紋斜率的定義適合工程應(yīng)用,按照簡正波頻散的定義適合理論分析.在有躍層的波導(dǎo)中,這兩種定義方式并不完全一致.由于簡正波頻散特性差異,按照頻散特性定義的波導(dǎo)不變量β會有許多不同的取值,這些β被稱為波導(dǎo)不變量的β譜.此時(shí)聲場干涉結(jié)構(gòu)應(yīng)該用多個(gè)β來描述,而以往的β提取算法只能給出一個(gè)最佳估計(jì)值,導(dǎo)致一些信息的丟失.借鑒圖像處理中的積分投影概念,將聲強(qiáng)圖像按照角度進(jìn)行積分投影,以分離不同斜率的條紋成分;然后對投影曲線進(jìn)行傅里葉變換,以分離不同間距的條紋成分,最終實(shí)現(xiàn)各個(gè)β譜值的分離.仿真和實(shí)驗(yàn)結(jié)果表明,β譜分離算法可以從聲場干涉結(jié)構(gòu)中有效地提取不同條紋成分的β,并將其映射到二維平面內(nèi),具有更強(qiáng)的噪聲抑制能力,所以能在更低的信噪比條件下使用.

波導(dǎo)不變量β譜,條紋斜率,條紋間距,積分投影

1 引言

由于簡正波的相互干涉,淺海波導(dǎo)中的聲場會形成穩(wěn)定的條紋結(jié)構(gòu).干涉條紋斜率與聲源距離、頻率之間的關(guān)系可以用一個(gè)參數(shù)——波導(dǎo)不變量β來描述[1,2].波導(dǎo)不變量β對水體聲速的小擾動不敏感[3],具有非常好的穩(wěn)定性,在聲源定位[4?9]、陣處理[10,11]、時(shí)間反轉(zhuǎn)等[12,13]方面均有應(yīng)用.在理論分析中,人們常采用波導(dǎo)不變量的另一種定義,即聲場簡正波相慢度差與群慢度差的比值[2].波導(dǎo)不變量是簡正波頻散特性的表征,而頻散特性與水體[2,5]、海底聲學(xué)參數(shù)[14]有關(guān).這一點(diǎn)可以用來反演環(huán)境參數(shù)[15?18].需要說明的是,基于波導(dǎo)不變量的Warping變換[19]或者時(shí)頻譜分析[20]手段可以實(shí)現(xiàn)不同簡正波頻散曲線的分離,對提高反演的準(zhǔn)確度非常重要.

以上基于β的應(yīng)用都依賴于條紋斜率的獲取.從聲強(qiáng)干涉結(jié)構(gòu)中提取條紋斜率的算法主要依靠圖像處理手段,比如Hough變換[4,21]采用參數(shù)映射的方法,能夠得到最符合聲強(qiáng)干涉條紋特征的β值.Radon變換[7,22]則是通過提取條紋的傾斜角度來得到β值.Rouse ff則是在二維快速傅里葉變換(2D-FFT)的基礎(chǔ)上給出了β的取值分布[23].另外,波導(dǎo)不變量β也可以通過基于干涉條紋的波束形成得到[24].需要說明的是,上述算法都是基于干涉條紋進(jìn)行的,要求有一定的信噪比.

但是復(fù)雜波導(dǎo)中,單一數(shù)值的波導(dǎo)不變量不能有效地描述干涉條紋.在有躍層的波導(dǎo)中,一定頻段內(nèi)簡正波的頻散特性(或者說β的取值)與簡正波號數(shù)密切相關(guān)[25].在有躍層的波導(dǎo)中,簡正波可以近似分為SRBR(surface-re fl ected-and-bottomre fl ected)簡正波、RBR(refracted-and-bottomre fl ected)簡正波兩大類,每一類簡正波的頻散特性比較相近.聲場往往是由這兩類簡正波共同組成的,此時(shí)SRBR簡正波和RBR簡正波之間的相互干涉可能會使干涉結(jié)構(gòu)變得非常復(fù)雜,不存在明顯的干涉條紋.這種情況下聲場干涉結(jié)構(gòu)需要有多個(gè)β才能描述,但以往的提取算法只能給出β的最佳估計(jì)值.

為了解決這個(gè)問題,本文按照簡正波的頻散特性定義了波導(dǎo)不變量的β譜,利用β譜的概念解釋了躍層波導(dǎo)中干涉條紋的變化規(guī)律,并給出了波導(dǎo)不變量β譜的分離算法.該算法可以將不同簡正波組合對應(yīng)的β譜分離出來,與以往的波導(dǎo)不變量提取算法相比,該算法能在更低的信噪比條件下工作.

2 波導(dǎo)不變量β譜的概念及分離算法

2.1 波導(dǎo)不變量β譜

根據(jù)簡正波理論,波導(dǎo)中距離為r,深度為z,角頻率為ω的聲場聲壓p可以表示為

式中km和φm(z)分別為聲場簡正波的水平波數(shù)和模態(tài)函數(shù),αm為模態(tài)衰減系數(shù),

是與聲源深度zs有關(guān)的模態(tài)激發(fā)強(qiáng)度.干涉聲強(qiáng)表達(dá)式為

干涉條紋的斜率為

聲強(qiáng)關(guān)于ω和r求導(dǎo)時(shí),相位是主要的變化因素,幅度項(xiàng)變化緩慢,可以忽略:

波導(dǎo)不變量根據(jù)干涉條紋的斜率定義[1,2]:

根據(jù)(4)式—(6)式,波導(dǎo)不變量的表達(dá)式為

式中Sp,mn和Sg,mn分別為第m,n號簡正波的相慢度差和群慢度差,

當(dāng)聲場由頻散特性相近的一組簡正波組成時(shí)[2],聲場簡正波群慢度與相慢度之間的函數(shù)關(guān)系與簡正波號數(shù)、頻率無關(guān),可以將(7)式中的Sp,mn和Sg,mn放到求和符號外面,(7)式可以進(jìn)一步簡化為

(6)式與(8)式的不同在于,(6)式是根據(jù)ω-r平面內(nèi)條紋斜率定義的,所以適用于工程應(yīng)用,而(8)式是根據(jù)簡正波頻散特性定義的,適用于理論分析.

對于Pekeris波導(dǎo)來說,簡正波的頻散特性非常一致,如圖1(b)所示,在Pekeris波導(dǎo)中相鄰簡正波的群慢度差Sg,m,m+1和相慢度差Sp,m,m+1在100—2500 Hz范圍內(nèi)都可以看作過原點(diǎn)的直線,并且所有簡正波都具有非常接近的斜率,這樣的情況下(6)式與(8)式定義是一致的.

但是對于有躍層的波導(dǎo)來說,簡正波的頻散特性與簡正波號數(shù)、頻率密切相關(guān)[25].圖1(c)為仿真所使用的波導(dǎo)環(huán)境,在15—30 m深度范圍內(nèi)存在躍層,圖1(d)所示為波導(dǎo)中的簡正波群慢度差Sg,m,m+1和相慢度差Sp,m,m+1的變化關(guān)系,可以看到它們不再是過原點(diǎn)的直線,其斜率不但隨頻率變化明顯,還與簡正波號數(shù)有關(guān).

在有躍層的波導(dǎo)中,根據(jù)簡正波頻散定義的波導(dǎo)不變量取值有多種可能性.每一組簡正波組合得到的βmn=?Sp,mn/Sg,mn各不相同,它們被定義為波導(dǎo)不變量的β譜.每一個(gè)譜βmn都有對應(yīng)的條紋成分,總聲場是多種條紋成分的疊加.由于β譜差異可能會非常明顯,ω-r平面內(nèi)總聲場干涉條紋會呈現(xiàn)出非常復(fù)雜的結(jié)構(gòu),甚至不再有條紋.這種情況下,很難根據(jù)條紋斜率進(jìn)行聲源定位.

圖1 (網(wǎng)刊彩色)(a)Pekeris波導(dǎo)中的水體和海底參數(shù);(b)Pekeris波導(dǎo)中相鄰簡正波相慢度差Sp,m,m+1和群慢度差Sg,m,m+1的函數(shù)關(guān)系;(c)有躍層波導(dǎo)中的水體和海底參數(shù);(d)有躍層波導(dǎo)中相鄰簡正波相慢度差Sp,m,m+1和群慢度差Sg,m,m+1的函數(shù)關(guān)系.(b)和(d)是寬帶聲場計(jì)算的結(jié)果,聲場頻率范圍為100—2500 HzFig.1.(color online)(a)Parameters of water column and sea bottom for a Pekeris waveguide;(b)functional relationship between modal phase slowness di ff erence Sp,m,m+1and modal group slowness di ff erence Sg,m,m+1for adjacent modes in the Pekeris waveguide of(a);(c)parameters of water column and sea bottom for a waveguide with thermocline;(d)functional relationship between modal phase slowness di ff erence Sp,m,m+1and modal group slowness di ff erence Sg,m,m+1for adjacent modes in the waveguide of(c).The results in(b)and(d)are calculated for a broad band,which is from 100 Hz to 2500 Hz.

β譜可以根據(jù)波導(dǎo)頻散特性直接由(8)式計(jì)算得到,但是為了利用β譜進(jìn)行聲源定位,需要給出一種能從聲強(qiáng)干涉結(jié)構(gòu)中提取分離各個(gè)β譜的算法——波導(dǎo)不變量的β譜分離算法.下面將給出β譜分離算法的原理和實(shí)現(xiàn)步驟.

2.2 波導(dǎo)不變量β譜的分離算法

對于聲場中的干涉條紋,借鑒圖像處理中的積分投影算法[26],將聲強(qiáng)圖像沿著不同角度的傾斜直線進(jìn)行積分,得到該圖像的一維投影曲線,一方面,這可以區(qū)分不同斜率的條紋成分,另一方面,不同簡正波組合的條紋間距也不同,這在投影曲線中體現(xiàn)為不同的周期,可以通過FFT進(jìn)行區(qū)分.利用干涉條紋的斜率和間距這兩個(gè)特征量,可以實(shí)現(xiàn)波導(dǎo)不變量β譜的分離.

2.2.1 聲強(qiáng)圖像的積分投影

對ω-r平面的聲強(qiáng)按照角度θ進(jìn)行積分投影.如圖2所示,取圖像中心(其頻率、距離值為(ω0,r0))為坐標(biāo)原點(diǎn),沿積分投影方向?yàn)閥軸,垂直于積分投影方向?yàn)閤軸,積分投影后得到的投影曲線是x的一維函數(shù).積分表達(dá)式為

式中γq是與簡正波號數(shù)q有關(guān)的常數(shù),稱為頻散因子,γqp是q,p號簡正波頻散因子之差,βqp是q,p號簡正波對應(yīng)的波導(dǎo)不變量β譜.(9)式所示的積分中,由于聲壓的幅度項(xiàng)Am變化很小,聲強(qiáng)中的非相干成分以及Bpq可以視為常數(shù).通過(9)式所示的積分,二維的聲強(qiáng)圖像被投影為一維的投影曲線.

積分投影可以分離聲強(qiáng)中斜率不同的條紋成分.條紋的波導(dǎo)不變量譜βmn與積分投影角度的關(guān)系為

當(dāng)積分投影角度為θmn時(shí),沿著積分路徑,只有βmn對應(yīng)的條紋成分是等相位的,其他譜值的條紋成分是不等相位的.經(jīng)過(9)式所示積分之后,譜值βmn的條紋成分得到保留,其他條紋成分則被剔除.詳細(xì)的推導(dǎo)見附錄,這里只給出最終積分結(jié)果.(9)式所示的積分投影結(jié)果為

式中ω0,r0是圖像中心的頻率和距離值,L(θmn,x)是積分路徑所對應(yīng)的長度,如圖2所示,其具體表達(dá)式見附錄.L(θmn,x)一般是一個(gè)梯形的窗函數(shù).由于積分區(qū)域有限,積分投影之后,譜值不等于βmn的條紋成分會有很小的殘留項(xiàng),即R(x).

對[?π/2,π/2]范圍內(nèi)的所有角度進(jìn)行積分投影,就可以將不同譜值的條紋成分予以分離,其β譜值與積分投影角度的關(guān)系參見(11)式.積分投影不能分離β譜接近的條紋成分,需要采用其他方法.

圖2 (網(wǎng)刊彩色)聲強(qiáng)圖像的積分投影示意圖.按照一定角度積分投影后,二維聲強(qiáng)圖像變?yōu)橥队扒€.積分投影方向?yàn)閥軸,垂直于積分投影方向?yàn)閤軸.當(dāng)積分投影方向與條紋一致時(shí),投影曲線的周期即為條紋之間的間距Fig.2.(color online)The sketch map for integral projection of acoustic intensity image.The intensity image is transformed into a curve after it is projected according to an arbitrary angle.The direction of projection is labeled as y axis,while the direction perpendicular to the projection is labeled as x axis.When the projection direction is along the striations,the projected curve has the period identical with space between striations.

2.2.2 投影曲線的FFT分析

干涉條紋在r方向上的間距,就是其簡正波組合的干涉跨度.聲強(qiáng)圖像積分投影后,投影曲線的周期正好對應(yīng)干涉跨度的投影.所以,β譜接近的條紋成分,可能會有差異明顯的條紋間距.對投影曲線進(jìn)行FFT運(yùn)算可以將這些β譜映射到不同的波數(shù)位置,從而實(shí)現(xiàn)分離.

根據(jù)(12)式,將總聲強(qiáng)按照角度θmn積分投影后,得到的投影曲線是一個(gè)周期信號,其周期為

(13)式說明m,n號簡正波干涉條紋的投影曲線,經(jīng)過FFT之后映射到了kmn/(2π|sinθmn|)的波數(shù)位置處.

需要補(bǔ)充的是,根據(jù)(13)式,投影曲線進(jìn)行FFT后,β譜映射到的波數(shù)位置不僅與干涉跨度有關(guān),還與投影角度有關(guān).在實(shí)際使用中,我們往往希望同一組簡正波形成的條紋成分對應(yīng)的波數(shù)位置僅由其干涉跨度決定.這可以通過下式來實(shí)現(xiàn):

式中k是投影曲線進(jìn)行FFT之后的波數(shù)軸坐標(biāo),K的單位是m?1/(2π),它是簡正波干涉跨度的倒數(shù).

經(jīng)過(14)式,m,n號簡正波的干涉條紋無論其斜率多大,都會被映射到kmn(ω0)/(2π)的波數(shù)位置.這在實(shí)際應(yīng)用中更加方便.比如引導(dǎo)聲源與未知聲源同時(shí)到達(dá)的聲場,因距離不同,同一組簡正波的條紋斜率也不同,進(jìn)行(14)式所示的坐標(biāo)變換后,它們會被聚集在相同的波數(shù)位置,但會對應(yīng)不同的β譜值,其譜值之比等于其距離之比(由于未知聲源聲場的譜值分離結(jié)果使用了引導(dǎo)聲源的距離信息,所以其β譜是不準(zhǔn)確的).這在利用引導(dǎo)聲源進(jìn)行未知聲源定位時(shí)非常方便.

以往的波導(dǎo)不變量提取算法并不能分離各個(gè)β譜,Hough變換是通過參數(shù)映射給出最優(yōu)的β值[4,21],Radon變換只能給出條紋斜率的最優(yōu)解[7,22],Rouse ff[23]在2D-FFT方法基礎(chǔ)上,給出了波導(dǎo)不變量β的取值分布情況.它們是基于條紋的提取算法,在聲強(qiáng)干涉結(jié)構(gòu)中沒有條紋時(shí)會失效,并且它們只是利用條紋的斜率信息,沒有考慮條紋的間距,是一維變換,而β譜值分離算法能將βmn映射到(β,K)二維平面內(nèi),是二維變換,具有更強(qiáng)的噪聲抑制能力,因而能在更低的信噪比條件下工作.

3 波導(dǎo)不變量譜值分離的仿真

仿真的波導(dǎo)環(huán)境采用圖1(c)所示的剖面.分為聲源位于躍層中(zs=20 m)和躍層下(zs=30 m)兩種情況進(jìn)行分析.兩種情況下,只有聲源位置不同,其他參數(shù)完全相同.海水深度為60 m,接收深度為50 m,接收距離為20—22 km,間距為20 m,頻率為700—800 Hz,間隔為1 Hz.

3.1 聲源位于躍層下(zs=30 m)

聲場的干涉條紋和簡正波幅度如圖3所示.此時(shí)聲源主要激發(fā)低號簡正波,前4號簡正波能量非常大.利用文獻(xiàn)[7]提出的Radon變換提取的條紋傾斜角度如圖4(a)所示,利用文獻(xiàn)[23]提出算法提取的β取值分布如圖4(b)所示.可以看到,傾斜角度為50?,波導(dǎo)不變量的取值為1.63,兩者結(jié)果是一致的.另外,在傾斜角度61?以及β=2.6處有個(gè)旁瓣.除此之外這兩種提取算法不能給出更多的信息.

根據(jù)β譜分離算法得到的結(jié)果如圖5所示.分離結(jié)果中出現(xiàn)非常多的峰值,分別標(biāo)記為A—F,它們的(β,K)坐標(biāo)分別為:

波數(shù)K單位為m?1/(2π).將它們的坐標(biāo)與簡正波頻散曲線得到的β(圖6)以及簡正波的水平波數(shù)差(表1)進(jìn)行比較可知,譜值A(chǔ)對應(yīng)1,2號簡正波的干涉,譜值B對應(yīng)3,4號簡正波的干涉,譜值C對應(yīng)2,4號簡正波的干涉,譜值D對應(yīng)1,4號簡正波的干涉,譜值E對應(yīng)2,5號簡正波的干涉,譜值F對應(yīng)2,6號簡正波的干涉.

圖3 (網(wǎng)刊彩色)水平均勻波導(dǎo)(圖1(c))中的仿真結(jié)果(聲源深度為30 m,聲源頻率為700—800 Hz,接收距離為20—22 km,接收深度為50 m)(a)各號簡正波接收幅度;(b)聲場干涉條紋Fig.3.(color online)The simulation result in a range-independent waveguide given in Fig.1(c)(the acoustic source is at a depth of 30 m with frequency band from 700 Hz to 800 Hz,and receivers are deployed from 20 km to 22 km at a depth of 50 m):(a)Averaged acoustic modal amplitude;(b)acoustic interference pattern.

圖4 聲源深度為30 m時(shí)從聲場干涉條紋中提取得到的條紋傾斜角度和β(a)Radon變換得到的條紋傾斜角度,算法見文獻(xiàn)[7];(b)根據(jù)2D-FFT算法得到的波導(dǎo)不變量取值分布,算法見文獻(xiàn)[23]Fig.4.The striation angle and waveguide invariant β extracted from interference pattern with source at 30 m deep:(a)Striation angle calculated by Radon transform as in Ref.[7];(b)distribution of β computed by the 2D-FFT based method which can be found in Ref.[23].

表1 圖1(c)波導(dǎo)中1—6號簡正波的水平波數(shù)差K(單位為m?1·(2π)?1,表中數(shù)據(jù)為700—800 Hz寬帶平均結(jié)果)Table 1.The horizontal wavenumber di ff erence for the fi rst six modes in the waveguide shown in Fig.1(c).Unit is m?1·(2π)?1and the result is averaged over the frequency band from 700 Hz to 800 Hz.

圖5 (網(wǎng)刊彩色)聲源深度為30 m時(shí)波導(dǎo)不變量β譜的分離結(jié)果.A—F分別表示分離得到的6個(gè)譜值,它們來源于1—6號簡正波之間的相互干涉Fig.5.(color online)The separated β spectrum when acoustic source is at the depth of 30 m.Labels A–F represent six spectrum values which can be identi fi ed,and they arise from the interference between the fi rst six acoustic modes.

圖6 (網(wǎng)刊彩色)根據(jù)簡正波頻散特性計(jì)算得到的不同簡正波組合的β譜值(簡正波的相慢度和群慢度根據(jù)Kraken計(jì)算得到,聲速剖面見圖1(c))Fig.6.(color online)Values for β spectrum calculated from acoustic modal dispersion for di ff erent combinations of modes.Modal phase slowness and group slowness are both calculated by Kraken model.The background sound speed pro fi le is given in Fig.1(c).

3.2聲源位于躍層中(zs=20 m)

聲源位于躍層中間時(shí),不但能激發(fā)高號簡正波,還會激發(fā)大量的中間模態(tài),簡正波幅度和干涉條紋結(jié)構(gòu)如圖7所示.可以看到干涉條紋結(jié)構(gòu)非常亂.條紋的傾斜角度[7]和波導(dǎo)不變量取值分布[23]如圖8所示.條紋傾斜角度有很多峰值,一部分集中在25?—31?范圍內(nèi),另一部分集中在接近垂直的傾角范圍內(nèi).波導(dǎo)不變量的取值分布最大概率取值為0.63,除此之外還有許多旁瓣.

β譜分離算法得到的結(jié)果如圖9所示,表2列出了高號簡正波的水平波數(shù)差,圖10所示為根據(jù)簡正波頻散得到的β譜值.將圖9中區(qū)域分為四部分,分別記為A1,B1,C1和D1.與聲源在30 m時(shí)明顯不同的是,圖9中在β>4的區(qū)域內(nèi)有很多峰值,比如D1區(qū)域,該部分對應(yīng)的是處于中間的簡正波,即第6—10號簡正波,其β譜是絕對值很大的正值或者負(fù)值,如圖10(a)所示.這些中間模態(tài)的β譜隨頻率劇烈變化,非常不穩(wěn)定.圖9中A1處對應(yīng)的坐標(biāo)為β=1.4,K=1.35×10?3,B1對應(yīng)坐標(biāo)為β=0.7—0.9,K=1.3×10?3—1.9×10?3,C1坐標(biāo)為β=0.7—0.9,K=2.6×10?3—3.6×10?3.波數(shù)K單位為m?1/(2π).

圖7 (網(wǎng)刊彩色)水平均勻波導(dǎo)(圖1(c))的仿真結(jié)果(聲源深度為20 m,聲源頻率為700—800 Hz,接收距離為20—22 km,接收深度為50 m)(a)各號簡正波接收幅度;(b)聲場干涉條紋Fig.7.(color online)The simulation result in a range-independent waveguide given in Fig.1(c).The acoustic source is at the depth of 20 m with frequency band from 700 Hz to 800 Hz,and receivers are deployed from 20 km to 22 km at 50 m deep:(a)Averaged acoustic modal amplitude;(b)acoustic interference pattern.

圖8 聲源深度為20 m時(shí)從聲場干涉條紋中提取得到的條紋傾斜角度和β(a)Radon變換得到的條紋傾斜角度,算法見文獻(xiàn)[7];(b)根據(jù)2D-FFT算法得到的波導(dǎo)不變量取值分布,算法見文獻(xiàn)[23]Fig.8.The striation angle and waveguide invariant β extracted from interference pattern with source at 20 m deep:(a)Striation angle calculated by Radon transform as in Ref.[7];(b)distribution of β calculated by the 2D-FFT based method which can be found in Ref.[23].

A1區(qū)域的譜值是第10號簡正波與第11號簡正波干涉的結(jié)果,即β10,11.第10,11號簡正波的譜值β10,11和水平波數(shù)差分別如圖10(b)和表2所示.B1區(qū)域和C1區(qū)域則可以用高號簡正波來解釋.B1區(qū)域包含了11 6 m 6 15時(shí)相鄰簡正波構(gòu)成的譜值,即βm,m+1;C1區(qū)域包含了m>16的βm,m+1,以及11 6 m 6 15的βm,m+2.11號以上的簡正波形成的波導(dǎo)不變量β譜及其水平波數(shù)差如圖10(c)和表2所示.

另外,通過圖9的結(jié)果也可以看到,頻散特性一致并且干涉跨度也接近的簡正波,其β譜并不能被分離.實(shí)際上這些簡正波的干涉條紋斜率非常一致,在應(yīng)用中也沒有必要進(jìn)行區(qū)分.

圖9 (網(wǎng)刊彩色)聲源深度為20 m時(shí)波導(dǎo)不變量β譜的分離結(jié)果.A1—D1表示分離得到的4個(gè)譜值Fig.9.(color online)The separated β spectrum when acoustic source is 20 m deep.Labels A1–D1represent the four identi fi ed spectrum values.

圖10 (網(wǎng)刊彩色)根據(jù)簡正波頻散特性計(jì)算得到的相鄰簡正波譜值βm,m+1(簡正波的相慢度和群慢度根據(jù)Kraken計(jì)算得到,聲速剖面見圖1(c))(a)5 6 m 6 9;(b)m=10;(c)m>11Fig.10.(color online)The spectrum value βm,m+1calculated from acoustic modal dispersion for adjacent modes:(a)5 6 m 6 9;(b)m=10;(c)m>11.Modal phase slowness and group slowness are both calculated by Kraken model.The sound speed pro fi le is given in Fig.1(c).

3.3 信噪比的影響

定義信噪比為信號的平均幅度與高斯白噪聲標(biāo)準(zhǔn)差的比值.在聲源深度為30 m時(shí),接收深度為50 m,沒有噪聲時(shí)的結(jié)果如圖3所示,信噪比為?9 dB時(shí)的結(jié)果如圖11所示.對比圖5和圖11(d),可以看到β譜分離算法仍能分辨各個(gè)β譜,而以往的波導(dǎo)不變量提取算法都不能得到準(zhǔn)確結(jié)果.

表2 圖1(c)波導(dǎo)中高號簡正波的水平波數(shù)差(單位為m?1·(2π)?1,表中數(shù)據(jù)為700—800 Hz寬帶平均結(jié)果)Table 2.The horizontal wavenumber di ff erence for higher-order modes in the waveguide shown in Fig.1(c).Unit is m?1·(2π)?1and the result is averaged over the frequency band from 700 Hz to 800 Hz.

圖11 (網(wǎng)刊彩色)信噪比為?9 dB時(shí)的仿真結(jié)果(a)聲場干涉條紋,計(jì)算參數(shù)同圖3;(b)根據(jù)2D-FFT算法提取的波導(dǎo)不變量取值分布;(c)Radon變換得到的條紋傾斜角度;(d)β譜分離結(jié)果Fig.11.(color online)The simulation result when SNR is?9 dB:(a)Acoustic interference pattern with the same parameters as those in Fig.3;(b)the distribution for β calculated using the 2D-FFT based method;(c)the striation angle calculated using the Radon transform;(d)the separated waveguide invariant β spectrum.

3.4關(guān)于仿真結(jié)果的說明

對于β譜分離算法來說,投影角度θ通過(11)式對應(yīng)為β.根據(jù)(11)式,對于某一個(gè)投影角度θ來說,聲源距離r與β是一一對應(yīng)的;但它們耦合在一起,必須給定其中一個(gè)參數(shù)的值才能確定另一個(gè)參數(shù)的值.在上述仿真中,將聲源距離r作為已知條件,從而得到了不同簡正波組合的β譜.而在實(shí)際使用中聲源距離r是未知的,此時(shí)需要已知波導(dǎo)的參數(shù),包括背景聲速剖面和海底參數(shù),然后根據(jù)波導(dǎo)參數(shù)確定不同簡正波組合的β譜,最后根據(jù)(11)式去反推聲源的距離r.或者借助引導(dǎo)聲源(其距離rg、頻帶、深度zg等參數(shù)已知),將未知聲源的聲場干涉結(jié)構(gòu)按照引導(dǎo)聲源的距離rg進(jìn)行β譜分離,它們的(相同簡正波組合對應(yīng)的)β譜值之比就是其距離之比.

3.5 與以往方法的比較

對于ω-r平面內(nèi)的聲強(qiáng),其Radon變換結(jié)果會在條紋傾斜角度上具有最大的標(biāo)準(zhǔn)差,文獻(xiàn)[7]正是據(jù)此來提取條紋傾斜角度,但沒有進(jìn)一步分析條紋波數(shù)域中的細(xì)節(jié).文獻(xiàn)[23]將2D-FFT之后的結(jié)果轉(zhuǎn)化到極坐標(biāo)系后,進(jìn)行了徑向(即沿著波數(shù)軸)上的能量積分,最后給出了波導(dǎo)不變量的取值分布,同樣沒有利用波數(shù)域內(nèi)的細(xì)節(jié)信息.也就是說,這兩種方法都沒有使用干涉條紋的間距這一信息.

β譜分離算法的步驟和原理是對圖像進(jìn)行不同角度的積分投影,然后對投影曲線進(jìn)行FFT運(yùn)算,投影角度θ(即條紋斜率)通過(11)式對應(yīng)為β,而FFT變換的波數(shù)坐標(biāo)則根據(jù)(14)式映射為簡正波波數(shù)差.這樣,就將ω-r平面內(nèi)的聲強(qiáng)干涉條紋斜率映射到β譜,條紋間距則映射到波數(shù)域內(nèi).文獻(xiàn)[7]和文獻(xiàn)[23]都是一維映射,而β譜分離算法是二維映射,將噪聲散布到更大的范圍內(nèi),因而β譜分離算法可以在更低的信噪比條件下工作.

4 實(shí)驗(yàn)數(shù)據(jù)的處理

實(shí)驗(yàn)數(shù)據(jù)來源于2012年青島附近海域的漁船實(shí)驗(yàn).實(shí)驗(yàn)海區(qū)海深為30 m,聲源為小漁船,漁船螺旋槳深度約為1.5 m,接收水聽器靠近海底的深度.實(shí)驗(yàn)海區(qū)測得的聲速剖面如圖12所示,躍層位于2—9 m的深度范圍,聲速從1530 m/s降到1492 m/s.實(shí)驗(yàn)過程中,漁船向遠(yuǎn)離水聽器的方向行駛,使用全球定位系統(tǒng)(GPS)記錄漁船的位置.

圖12(a)實(shí)驗(yàn)海區(qū)的聲速剖面;(b)350 Hz時(shí)前3號簡正波的模態(tài)函數(shù)Fig.12.(a)Sound speed pro fi le in the experiment area;(b)acoustic modal functions at the frequency of 350 Hz.

圖13 為漁船噪聲的聲強(qiáng)譜,聲強(qiáng)譜是接收信號FFT幅度譜的平方,圖中橫軸代表漁船與水聽器之間的距離,根據(jù)GPS數(shù)據(jù)計(jì)算得到.實(shí)驗(yàn)中,水聽器記錄的漁船噪聲中有明顯的線譜成分,而漁船噪聲連續(xù)譜的信噪比非常低,圖中看不到任何條紋成分.

對于圖13所示的聲強(qiáng)譜,利用文獻(xiàn)[23]中的方法得到的(歸一化的)波導(dǎo)不變量取值分布如圖14(a)所示.可以看到由于線譜成分的存在,波導(dǎo)不變量的取值分布在β=0處有很大的峰值,而在除了β=0之外的區(qū)域,β取值分布曲線存在很多極大值,但是這些極大值并不突出,不能排除這些極大值是環(huán)境噪聲干擾的結(jié)果.

采用β譜分離算法得到的分離結(jié)果如圖14(b)所示.β=0附近的區(qū)域是線譜成分導(dǎo)致的,除此之外,還有A,B兩個(gè)區(qū)域存在明顯的峰值.實(shí)際上,圖14(a)中的分布曲線就是圖14(b)沿著橫軸方向積分疊加的結(jié)果.由于圖14(a)中分布曲線沒有利用波數(shù)域內(nèi)的細(xì)節(jié)信息,因而得不到有效的β值.

圖13 (網(wǎng)刊彩色)水聽器記錄的漁船聲強(qiáng)譜(橫軸的距離根據(jù)GPS數(shù)據(jù)計(jì)算得到)Fig.13.(color online)The sound intensity spectrum recorded by a hydrophone.The range axis is computed from GPS data.

圖14 (網(wǎng)刊彩色)(a)根據(jù)文獻(xiàn)[23]得到的波導(dǎo)不變量取值分布;(b)波導(dǎo)不變量的β譜分離結(jié)果Fig.14.(color online)(a)The distribution of β calculated by the method in Ref.[23];(b)the separated β spectrum.

根據(jù)文獻(xiàn)[27,28]可知,青島附近海域的海底地質(zhì)類型為砂-粉砂類型,所以選定海底聲速值為1700 m/s,密度為2.0 g/cm3,海底衰減系數(shù)為0.2 dB/λ.聲速剖面如圖12所示,根據(jù)Kraken可以計(jì)算出各號簡正波的水平波數(shù),確定各號簡正波的波導(dǎo)不變量β譜以及水平波數(shù)差,結(jié)果分別如圖15和表3所示.

表31 —5號簡正波的水平波數(shù)差K(單位為m?1·(2π)?1,表中數(shù)據(jù)為300—400 Hz寬帶平均結(jié)果)Table 3.The horizontal wavenumber di ff erence for modes 1–5.Unit is m?1·(2π)?1and the result is averaged over the frequency band from 300 Hz to 400 Hz.

圖14(b)中的A,B兩個(gè)區(qū)域可以用2,3,4號簡正波的β譜來解釋.B區(qū)域的峰值坐標(biāo)為(1.35,7×10?3);A區(qū)域范圍較大,峰值在A的右下側(cè),為(1.2,4.2×10?3).根據(jù)圖15和表3,B區(qū)域代表的是2,4號簡正波形成的譜值β2,4,而A區(qū)域中的峰值則是β3,4.另外,2,3號簡正波形成的譜值β2,3為1.45,對應(yīng)的波數(shù)為3×10?3m?1/(2π),對應(yīng)A區(qū)域中的左上部分.由于2,3號簡正波干涉跨度與3,4號簡正波干涉跨度差別不大,所以兩者一起構(gòu)成了A區(qū)域.由于聲源很淺,激發(fā)的1號簡正波強(qiáng)度比其他簡正波小得多,所以在分析中可以不考慮1號簡正波.

圖15 (網(wǎng)刊彩色)實(shí)驗(yàn)海域波導(dǎo)環(huán)境中根據(jù)頻散特性得到的波導(dǎo)不變量β譜Fig.15.(color online)The spectrum values for β calculated by modal dispersion in a waveguide whose parameters are set according to the experiment area.

5 結(jié)論

波導(dǎo)不變量既可以根據(jù)聲場ω-r平面內(nèi)的干涉條紋斜率來定義,也可以根據(jù)簡正波的頻散來定義.在Pekeris波導(dǎo)中,簡正波的相慢度差與群慢度差存在一致的比例關(guān)系,所以這兩種定義方式具有等價(jià)性.但在復(fù)雜波導(dǎo)中,簡正波相慢度差和群慢度差之比不但與頻率有關(guān),還與簡正波號數(shù)有關(guān).由于簡正波的頻散特性差異明顯,聲場中的干涉結(jié)構(gòu)會比較復(fù)雜,甚至不再存在條紋結(jié)構(gòu).此時(shí)按照條紋斜率定義的波導(dǎo)不變量不再是一個(gè)單一的數(shù)值,不能直接在工程上使用.

為了解決這一問題,本文將不同簡正波組合的相慢度差與群慢度差之比稱為β譜,并給出了β譜的分離算法.波導(dǎo)不變量的各個(gè)β譜可以直接根據(jù)簡正波頻散給出,但是需要給出從ω-r平面內(nèi)分離出各個(gè)β譜的算法以滿足工程需求.β譜分離算法的實(shí)現(xiàn)步驟是先對聲強(qiáng)圖像進(jìn)行積分投影,然后對投影曲線進(jìn)行傅里葉變換;投影角度(即條紋斜率)按照特定公式映射到β譜值,條紋間距映射到簡正波波數(shù)差,從而實(shí)現(xiàn)不同β譜的分離.

文獻(xiàn)中的波導(dǎo)不變量提取算法都將β視為一個(gè)單一數(shù)值的量,只能給出β的最佳解.這些算法大多為基于條紋的算法,只利用了干涉條紋的斜率信息.不同于以往的提取方法,β譜分離算法沒有丟失ω-r平面內(nèi)的任何信息,而且是二維平面上的映射,能更加有效地抑制噪聲的影響.仿真和實(shí)驗(yàn)分析結(jié)果表明,該方法能有效地實(shí)現(xiàn)不同β譜的分離,同時(shí)能在更低的信噪比條件下工作.

附錄聲強(qiáng)線積分結(jié)果的相關(guān)推導(dǎo)

ω-r平面內(nèi)的聲強(qiáng)按照圖2所示y方向進(jìn)行積分投影,積分公式參見(9)式.附錄將給出投影曲線的表達(dá)式,投影角度為θmn.

為了方便討論,將(9)式分為兩部分:第一部分記為G0(x),對應(yīng)譜值為βmn的干涉條紋的積分投影;第二部分記為G1(x),對應(yīng)其他譜值干涉條紋的積分投影.

首先證明沿著y方向G0(x)中的被積函數(shù)是同相疊加,而G1(x)中的被積函數(shù)不是同相疊加.以過原點(diǎn)的積分路徑為例說明,即x=0,此時(shí)有

一般來說,聲強(qiáng)平面的頻率、距離寬度Bω和Br遠(yuǎn)小于其頻率、距離中心值ω0和r0,根據(jù)泰勒級數(shù),保留一階小量,(A4)式可以簡化為

將(A3)式代入可得

下面推導(dǎo)G0(x)的積分表達(dá)式.G0(x)的積分表達(dá)式實(shí)際上是聲強(qiáng)在x軸上的值與一個(gè)窗函數(shù)L(θmn,x)的乘積,該窗函數(shù)來源于積分路徑的長短,如圖2所示.根據(jù)圖2,x軸上頻率和距離坐標(biāo)可以表示為

聯(lián)合(A3)式,得到x軸上的被積函數(shù)相位為

將(A7)代入,可得

因此,G0(x)的表達(dá)式為

根據(jù)θmn與a的大小關(guān)系(圖A1),積分路徑長度L(θmn,x)表達(dá)式可以寫為幾種形式.

圖A1積分路徑長度L(θmn,x)Fig.A1.The illustration map for integral length L(θmn,x).

需要說明的是,由于積分區(qū)域是有限的,(A1)式中G1(x)積分結(jié)果并不為零.實(shí)際上,傾斜角度與θmn差距較大的干涉條紋成分,其被積函數(shù)的相位在積分路徑上變化非?？?其積分結(jié)果可以看作是零,而對于那些傾斜角度與θmn差距不大的干涉條紋成分,其被積函數(shù)的相位在積分路徑上的變化可能并不大,因而其積分結(jié)果并不等于零.所以G1(x)在積分區(qū)域有限的情況下積分結(jié)果記為R(x),通常情況下,R(x)值很小,可以忽略.

[1]Chuprov S D(Ed.By Brekhovskikh L M and Andreevoi I B)1982 Ocean Acoustics:Current State(Moscow:Nauka)p71

[2]Brekhovskikh L M,Lysanov Yu P 2002 Fundamentals of Ocean Acoustics(Moscow:AIP Press)pp143–148

[3]Cockrell K L,Schmidt H 2011 J.Acoust.Soc.Am.130 72

[4]Turgut A,Orr M,Rouse ffD 2010 J.Acoust.Soc.Am.127 73

[5]Cockrell K L,Schmidt H 2010 J.Acoust.Soc.Am.127 2780

[6]Thode A M,Kuperman W,Dspain G L,Hodgkiss W 2000 J.Acoust.Soc.Am.107 278

[7]Thode A M 2000 J.Acoust.Soc.Am.108 1582

[8]Zhao Z D,Wang N,Gao D Z,Wang H Z 2010 Chin.Phys.Lett.27 064301

[9]Yang T C 2003 J.Acoust.Soc.Am.113 1342

[10]Tao H,Krolik J L 2008 J.Acoust.Soc.Am.123 1338

[11]Su X X,Zhang R H,Li F H 2006 Acta Acustica 31 305(in Chinese)[蘇曉星,張仁和，李風(fēng)華2006聲學(xué)學(xué)報(bào)31 305]

[12]Hodgkiss W S,Song H C,Kuperman W A 1999 J.Acoust.Soc.Am.105 1597

[13]Kim S,Kuperman W A,Hodgkiss W S,Song H C,Edelmann G F,Akal T 2003 J.Acoust.Soc.Am.114 145

[14]Shang E C,Wu J R,Zhao Z D 2012 J.Acoust.Soc.Am.131 3691

[15]Potty G,Miller J,Lynch J,Smith K 2000 J.Acoust.Soc.Am.108 973

[16]Bonnel J,Gervaise C,Nicolas B,Mars J 2012 J.Acoust.Soc.Am.131 119

[17]Heaney K D 2004 IEEE J.Ocean.Eng.29 88

[18]Ren Q Y,Hermand J P 2013 J.Acoust.Soc.Am.133 82

[19]Lu L C,Ma L 2015 Acta Phys.Sin.64 024305(in Chinese)[鹿力成,馬力2015物理學(xué)報(bào)64 024305]

[20]Bonnel J,Le T G,Nicolas B,Jéróme I M 2013 IEEE Signal Proc.Mag.30 120

[21]Hough P V C 1962 US Patent No.3069654

[22]Deans S R 1983 The Radon Transform and Some of Its Applications(New York:Wiley)

[23]Rouse ffD 2001 Waves in Random Media 11 377

[24]Rouse ffD,Zurk L M 2011 J.Acoust.Soc.Am.130 76

[25]Zhao Z D,Wu J R,Shang E C 2015 J.Acoust.Soc.Am.138 223.

[26]Zhang W Z,Chen Q,Du D,Sun Z G 2005 J.Tsinghua Univ.(Sci.&Tech.)45 1446(in Chinese)[張文增,陳強(qiáng),都東,孫振國2005清華大學(xué)學(xué)報(bào)45 1446]

[27]Liu J,Saito Y,Kong X H,Wang H,Wen C,Yang Z G,Nakashima R 2010 Marine Geology 278 54

[28]Zhang J Q 2012 Ph.D.Dissertation(Qingdao:Ocean University of China)(in Chinese)[張軍強(qiáng)2012博士學(xué)位論文(青島:中國海洋大學(xué))]

PACS:43.30.+m,43.60.+dDOI:10.7498/aps.66.114301

Concept of waveguide invariant spectrum and algorithm for its extraction?

Song Wen-Hua1)Wang Ning1)?Gao Da-Zhi1)Wang Hao-Zhong1)Qu Ke2)
1)(Department of Marine Technology,Ocean University of China,Qingdao 266100,China)
2)(School of Electronic and Information Engineering,Guangdong Ocean University,Zhanjiang 524088,China)

30 December 2016;revised manuscript

19 March 2017)

There are two kinds of de fi nitions for waveguide invariant β.One is de fi ned according to the striation slope of acoustic interference patterns,and the other is de fi ned on the basis of dispersion characteristics of acoustic modes.The fi rst de fi nition is appropriate for engineering applications,while the second is suitable for theoretical analysis.However,the two de fi nitions are not consistent with each other for a waveguide with thermoclines,because modal dispersion in such a waveguide can be very di ff erent for di ff erent modes and di ff erent frequencies.In such cases,the waveguide invariant de fi ned according to modal dispersion can take many di ff erent values,which are referred to as the spectrum of waveguide invariant(β spectrum for short)in the paper.Each β spectrum can be related to some interference striation patterns with corresponding striation slopes.The sound fi eld is composed of many modes,so the interference pattern is the summation of many components of di ff erent striation slopes and may be very complicated as a result of the diversity of β spectrum.In such a case one single β is not able to describe the complicated interference pattern adequately;instead multiple values of β spectrum are required.From the point of view of engineering application,however,the present β-extracting methods can only give one optimal value,and thus a lot of information is lost.In this paper an algorithm for doing so,called β spectrum separation technique,is proposed.By adopting the concept of integral projection used in digital image processing,the image of acoustic intensity is projected at di ff erent angles to separate out the striations of di ff erent slopes;and then fast Fourier transform(FFT)is applied to the projected curve in order to isolate striations of di ff erent spacing from each other.The values for β spectrum can be computed according to striation slopes,which are also mapped into the positions of their corresponding acoustic modal horizontal wavenumber di ff erences in the wavenumber domain.The applicability of this algorithm for the extraction of β spectrum is tested and veri fi ed by simulation results and experiment data.It is shown that the algorithm can separate out each β spectrum of di ff erent intensity components from acoustic interference structure.The algorithm maps β spectrum into a two-dimensional plane thereby being able to suppress noise more e ff ectively and work in the condition of low signal-to-noise ratio compared with the already-existing β-extracting algorithms.

waveguide invariant spectrum,striation slope,striation separation distance,integral projection

10.7498/aps.66.114301

?國家自然科學(xué)基金(批準(zhǔn)號:11674294,11374271)和廣東省自然科學(xué)基金(批準(zhǔn)號:2014A030310256)資助的課題.

?通信作者.E-mail:wangyu@public.qd.sd.cn

?2017中國物理學(xué)會Chinese Physical Society

http://wulixb.iphy.ac.cn

*Project supported by the National Natural Science Foundation of China(Grant Nos.11674294,11374271)and the Natural Science Foundation of Guangdong Province,China(Grant No.2014A030310256).

?Corresponding author.E-mail:wangyu@public.qd.sd.cn