王 帥,王 斌,范 軍
(上海交通大學(xué) 海洋工程國(guó)家重點(diǎn)實(shí)驗(yàn)室,上海 200240)
作為主動(dòng)聲納的重要判據(jù),海洋環(huán)境中復(fù)雜目標(biāo)聲散射特性一直是海洋計(jì)算聲學(xué)的研究熱點(diǎn)與挑戰(zhàn),尤其是在淺海環(huán)境中,目標(biāo)聲散射特性會(huì)變得異常復(fù)雜。針對(duì)海洋環(huán)境中復(fù)雜目標(biāo)聲散射建模問(wèn)題,國(guó)內(nèi)外開(kāi)展了深入而廣泛的研究,涉及到了幾乎所有的海洋計(jì)算聲學(xué)方法,例如:有限元、邊界元、物理聲學(xué)、拋物方程、簡(jiǎn)正波方法、波數(shù)積分以及射線聲學(xué)等。
Ingenito[1]采用積分方程方法研究理想波導(dǎo)中目標(biāo)的聲散射,為單體散射近似奠定了理論基礎(chǔ);Sarkissian[2]在此基礎(chǔ)上用替代法改進(jìn)了表面積分方程;Jensen等[3]結(jié)合Ingenito的方法忽略了界面的多次散射,通過(guò)用物體的平面波散射函數(shù)來(lái)表示該物體,需要計(jì)算格林函數(shù),而任何涉及其中散射積分的計(jì)算都是相當(dāng)復(fù)雜的,即使是簡(jiǎn)單的球體或圓柱;Schmidt[4]和Abawi[5-6]提出了一種疊加方法,將散射場(chǎng)近似等于目標(biāo)內(nèi)表面上大量點(diǎn)源產(chǎn)生的場(chǎng),但為了獲得點(diǎn)源強(qiáng)度,表面格林函數(shù)的梯度必須精確計(jì)算,這在通常情況下是非常困難和低效的;Belibassakis等[7]提出了逐級(jí)耦合模態(tài)法,用于計(jì)算目標(biāo)必須為軸對(duì)稱(chēng)的波導(dǎo)中目標(biāo)的散射場(chǎng),而克服上述障礙的有限元方法不能直接解決低速大存儲(chǔ)的問(wèn)題;基于單散射模型和射線追蹤技術(shù),Chen 和Fan[8]提出了一種幾何聲學(xué)方法來(lái)計(jì)算淺水目標(biāo)回波,但是由于Kirchhoff 近似在高頻和聲影區(qū)的限制以及射線理論只適用于高頻,該方法不適用于涉及目標(biāo)散射聲影區(qū)域的波導(dǎo)中前向散射情況;Fan[9]結(jié)合了邊界元方法和簡(jiǎn)正波聲傳播模型建立了聯(lián)合模型,可以計(jì)算淺海波導(dǎo)在聲速剖面隨深度變化、與距離無(wú)關(guān)的情況下,復(fù)雜形狀目標(biāo)的散射聲場(chǎng)。另外結(jié)合了簡(jiǎn)正波聲傳播模型以及目標(biāo)散射的Kirchhoff近似方法建立了混合模型,一定程度上解決了大尺寸、復(fù)雜形狀目標(biāo)在淺海波導(dǎo)中散射聲場(chǎng)的工程預(yù)報(bào)問(wèn)題,僅適用于當(dāng)目標(biāo)遠(yuǎn)離波導(dǎo)界面情況,即目標(biāo)和界面兩者之間多次散射可以近似忽略。同時(shí),因?yàn)槲锢砺晫W(xué)方法不能計(jì)算目標(biāo)聲影區(qū)的散射,所以方法也不適用于淺海波導(dǎo)中的目標(biāo)前向散射的情況。Jiang[10]提出了基于格林函數(shù)-拋物方程的淺海環(huán)境中目標(biāo)聲輻射計(jì)算方法,不可避免地引入拋物方程誤差。
當(dāng)目標(biāo)與邊界的距離相對(duì)于聲波波長(zhǎng)可以比擬時(shí),不能忽略多次散射效應(yīng),基于單體散射近似的方法誤差很大。Sarkissian[2]發(fā)現(xiàn)即使在較低頻率,例如目標(biāo)離界面為10倍目標(biāo)半徑,多次散射的影響已不能忽略;Makris[12]發(fā)現(xiàn)對(duì)比自由空間目標(biāo)的散射,淺海中目標(biāo)散射最重要的區(qū)別是,當(dāng)目標(biāo)靠近界面時(shí),目標(biāo)和界面之間的多次散射將起重要作用;Fan[8]論證了淺海信道中,散射聲場(chǎng)中經(jīng)界面反射等多途散射聲的貢獻(xiàn)不可忽略,否則在散射聲壓幅值、散射方向性、頻率響應(yīng)特性等多方面的誤差都較大,這種區(qū)別在目標(biāo)靠近界面時(shí)尤其明顯;Qian[11]也研究了目標(biāo)靠近海面和海底時(shí),其耦合頻率和自由場(chǎng)差別較大,只有當(dāng)目標(biāo)達(dá)到一定深度時(shí),才逐漸和自由場(chǎng)結(jié)果基本趨于一致。
目前,考慮了目標(biāo)與界面之間多次散射的計(jì)算方法主要以分離變量或者有限元為主,但局限于理想海洋環(huán)境、簡(jiǎn)單目標(biāo)或者計(jì)算效率低下。本文提出了一種基于虛擬柱島、有限元和簡(jiǎn)正波混合數(shù)值方法。這種方法利用分離變量或者有限元的優(yōu)勢(shì),即利用有限元方法精確計(jì)算目標(biāo)附近的散射聲壓情況,并根據(jù)簡(jiǎn)正波快速計(jì)算遠(yuǎn)場(chǎng)散射聲場(chǎng)。該方法的優(yōu)勢(shì)有:(1)考慮了目標(biāo)和界面的多次散射,可以計(jì)算目標(biāo)靠近海面或者海底的情況;(2)計(jì)算速度相對(duì)于有限元方法較快,計(jì)算效率更高。同時(shí),該方法適用于低頻散射,尤其是極低頻的條件,可用于復(fù)雜目標(biāo)在淺海信道中前向散射問(wèn)題。
對(duì)于圖1 所示淺海信道中復(fù)雜目標(biāo)聲散射問(wèn)題,本文提出的虛擬柱島、有限元和簡(jiǎn)正波混合方法分為三個(gè)步驟:(1)利用簡(jiǎn)正波方法計(jì)算出無(wú)目標(biāo)情況下海洋信道中任意一點(diǎn)的入射聲場(chǎng);(2)用一個(gè)上界面是水面、下界面是海底的虛擬柱島將目標(biāo)及其所在區(qū)域從淺海信道中分離出來(lái),利用有限元方法計(jì)算該區(qū)域內(nèi)散射聲場(chǎng)空間分布,其中背景聲場(chǎng)為步驟(1)計(jì)算的入射聲場(chǎng)、虛擬柱島表面為無(wú)反射邊界;(3)利用周向諧波、深度方向簡(jiǎn)正波正交性,根據(jù)虛擬柱島表面散射聲場(chǎng)確定散射聲場(chǎng)展開(kāi)系數(shù)以及外部區(qū)域散射聲場(chǎng)。
圖1 淺海信道目標(biāo)聲散射模型Fig.1 Model of acoustic scattering from objects in shallow waters
分別對(duì)上述三個(gè)步驟展開(kāi)敘述,時(shí)間因子為e-jωt。
1.1.1 海洋信道入射聲場(chǎng)
假設(shè)水層聲速為c1、密度為ρ1,聲源深度為z0,海底聲速為c2、密度為ρ2、深度為h,在深度H處近似滿(mǎn)足聲壓自由條件。以聲源為中心建立圓柱坐標(biāo)系,其中z= 0為海面、z軸正方向指向海底。對(duì)于單極子點(diǎn)聲源,其輻射聲場(chǎng)具有圓周對(duì)稱(chēng)性,根據(jù)非齊次Helmholtz方程[13-14]、分離變量法以及邊界條件,可以得到入射聲場(chǎng)遠(yuǎn)場(chǎng)的表達(dá)式為
1.1.2 虛擬柱島內(nèi)散射聲場(chǎng)
目標(biāo)表面和外部流體相接的耦合面上法向振速連續(xù)、切向應(yīng)力連續(xù)以及切向應(yīng)力為零,由此可寫(xiě)出目標(biāo)與流體的耦合方程為
其中,K、C、M分別代表剛度矩陣、阻尼矩陣以及質(zhì)量矩陣,下標(biāo)w、s和τ分別表示聲學(xué)、力學(xué)以及耦合,u代表結(jié)構(gòu)的位移,ω代表角頻率,ps、pr分別代表彈性和剛性散射聲場(chǎng),兩者之和為總散射聲場(chǎng)。
在虛擬柱島表面上,采用完美匹配層(PML)技術(shù)[11]模擬無(wú)反射邊界,該技術(shù)通過(guò)增加吸收系數(shù)將波動(dòng)方程轉(zhuǎn)換為吸收層控制方程,即
式中,σi為吸收系數(shù),p為匹配層中的聲壓。
1.1.3 虛擬柱島外散射聲場(chǎng)
相對(duì)于單極子源而言,水下目標(biāo)散射聲場(chǎng)空間分布更為復(fù)雜,但其本質(zhì)仍是聲輻射問(wèn)題,只是聲源更為復(fù)雜而已。因此,散射聲場(chǎng)也可以表示為滿(mǎn)足遠(yuǎn)場(chǎng)輻射、邊界條件和連續(xù)條件的柱坐標(biāo)系下波動(dòng)方程形式解[17],即
根據(jù)步驟(2)計(jì)算得到的虛擬柱島表面散射聲壓,利用周向諧波、深度方向簡(jiǎn)正波正交性即可確定Anm,即
其中,r0為虛擬柱島半徑。
1.2.1 典型算例
計(jì)算參數(shù)如下:海水聲速為1500 m/s、密度為1000 kg/m3;海底深度為100 m、聲速為1700 m/s、密度為1500 kg/m3、吸收系數(shù)為0.5 dB/λ[16]、截止深度為300 m;單位幅度單極子聲源頻率為75 Hz,剛性球半徑為10 m,兩者深度均為50 m、水平相距2500 m。圖2給出了無(wú)目標(biāo)情況下聲源所在垂直平面內(nèi)入射聲場(chǎng)空間分布。
圖2 淺海信道中入射聲場(chǎng)空間分布Fig.2 Spatial distribution of incident sound field in shallow waveguide
根據(jù)仿真結(jié)果可以看到,聲波能量大部分集中在海水層中,而海底層中聲波能量隨著水平距離增加迅速衰減。圖3 給出了以剛性球?yàn)橹行?、半徑?0 m 的虛擬柱島內(nèi)剛性球球心與聲源所在垂直平面內(nèi)散射聲場(chǎng)的空間分布。其中,水平向右為x軸正方向,即前向散射,PML厚度為5 m。
圖3 全海深虛擬柱島內(nèi)散射聲壓級(jí)Fig.3 Scattering pressure level in virtual cylindrical column in double layers
可以觀察到,散射聲場(chǎng)聲能量也主要集中在海水層,透入海底層的聲波能量會(huì)隨著深度的增加而快速衰減、在PML近似無(wú)反射。為了更加直觀地觀察目標(biāo)散射聲場(chǎng)空間分布規(guī)律,圖4給出了海水層內(nèi)虛擬柱島表面散射聲場(chǎng),可直觀地看出前向散射聲場(chǎng)幅度更大,符合理論規(guī)律。
圖4 海水層虛擬柱島表面散射聲壓Fig.4 Scattering pressure on the surface of virtual cylindrical column in water layer
根據(jù)式(7)計(jì)算待定系數(shù)Anm,代入式(6)計(jì)算海洋信道中目標(biāo)的散射聲場(chǎng),圖5給出了虛擬柱島外剛性球球心與聲源所在垂直平面內(nèi)散射聲場(chǎng)的空間分布,其中周向、深度方向截?cái)嚯A次分別為10和28。
圖5 淺海信道中散射聲場(chǎng)空間分布Fig.5 Spatial distribution of scattered sound field in shallow waveguide
可以觀察到,散射聲場(chǎng)隨著水平距離增加呈現(xiàn)衰減趨勢(shì);由于高階次簡(jiǎn)正波衰減速度快,深度方向散射聲場(chǎng)起伏度逐漸削弱。散射聲波每階簡(jiǎn)正波可以分解為一對(duì)上行波和下行波,經(jīng)上、下界面多次反射向外傳播。至此,完成了三個(gè)步驟的計(jì)算,并得到了淺海信道中目標(biāo)散射聲場(chǎng)。
1.2.2 仿真驗(yàn)證
由于1.2.1 小節(jié)仿真算例沒(méi)有解析解,為此通過(guò)對(duì)比有限元計(jì)算結(jié)果驗(yàn)證本文方法的正確性,其中本文方法截?cái)嘤虬霃綖?5 m、有限元方法截?cái)嘤虬霃椒謩e為15 m、25 m 和35 m,其它計(jì)算參數(shù)和1.2.1小節(jié)相同。圖6對(duì)比了不同水平距離上前向散射聲壓級(jí)隨深度的變化規(guī)律。
圖6 不同距離上前向散射聲壓級(jí)隨深度變化規(guī)律Fig.6 Forward scattering sound pressure levels with respect to depth at different distances
可以看出,兩種方法計(jì)算結(jié)果吻合較好,驗(yàn)證了本文方法計(jì)算淺海信道中目標(biāo)散射聲場(chǎng)的準(zhǔn)確性。另外,從計(jì)算效率和計(jì)算時(shí)間上來(lái)看,三個(gè)截?cái)嘤虬霃较掠邢拊椒ǚ謩e需要15 min、70 min 和200 min,而本文方法僅需要15 min,其計(jì)算速度、效率優(yōu)勢(shì)隨著水平距離的增加而迅速增大。
考慮到界面與目標(biāo)之間多次散射,本文采用有限元方法計(jì)算虛擬柱島內(nèi)散射聲場(chǎng),為此虛擬柱島半徑與深度直接決定了本文方法的計(jì)算速度與效率。
對(duì)于有限元方法而言,計(jì)算域截?cái)嘣叫?、?jì)算速度越高,剛好覆蓋目標(biāo)的虛擬柱島最為理想。然而,虛擬柱島表面PML對(duì)目標(biāo)散射場(chǎng)瞬失波吸收效果很差,計(jì)算誤差隨著虛擬柱島半徑的減小迅速增大。為此,計(jì)算精度和速度是矛盾的,需要結(jié)合實(shí)際要求進(jìn)行平衡和取舍。
圖7 給出了不同虛擬柱島半徑情況下剛性球球心與聲源所在垂直平面散射聲場(chǎng),其半徑分別為10 m和60 m,即二分之一聲波波長(zhǎng)、三倍聲波波長(zhǎng),其它計(jì)算參數(shù)與圖3相同。
圖7 不同半徑虛擬柱島情況下目標(biāo)附近散射聲壓級(jí)Fig.7 Scattering pressure levels in virtual column with different radii
對(duì)比圖3和圖7可以明顯看出,虛擬柱島半徑為一倍或三倍聲波波長(zhǎng)時(shí)目標(biāo)附近散射聲場(chǎng)計(jì)算結(jié)果吻合較好,相比之下虛擬柱島半徑為二分之一波長(zhǎng)時(shí)計(jì)算偏差較大。
圖8 給出了根據(jù)圖3、圖7 虛擬柱島表面散射聲壓計(jì)算得到的目標(biāo)前向散射聲場(chǎng)隨深度的變化規(guī)律,其中,水平距離為35 m。
圖8 不同半徑虛擬柱島情況下目標(biāo)前向散射聲場(chǎng)隨深度變化規(guī)律Fig.8 Forward scattering fields in virtual columns of different radii with respect to depth
可以看到,隨著虛擬柱島半徑增加,淺海信道中目標(biāo)散射聲場(chǎng)計(jì)算結(jié)果逐漸收斂,例如半徑為二分之一聲波波長(zhǎng)時(shí)計(jì)算結(jié)果相對(duì)于三倍聲波波長(zhǎng)偏差接近15 dB,半徑擴(kuò)大至一倍聲波波長(zhǎng)時(shí)計(jì)算結(jié)果偏差縮小至2 dB。誤差主要來(lái)源于PML,它只能吸收向外傳播的聲波,而目標(biāo)散射聲場(chǎng)近場(chǎng)及其受上、下界面影響向內(nèi)傳播的聲波則被PML 放大,導(dǎo)致計(jì)算誤差增大,因此虛擬柱島半徑不可過(guò)小,建議在一倍聲波波長(zhǎng)以上。
在仿真的計(jì)算時(shí)間上,虛擬柱島半徑10 m、20 m 和60 m 情況下計(jì)算時(shí)間分別為5 min、15 min 和34 h,即虛擬柱島半徑越大,計(jì)算時(shí)間越長(zhǎng),但計(jì)算精度越高。因此,需要結(jié)合實(shí)際情況,權(quán)衡計(jì)算時(shí)間和計(jì)算精度。
由于聲波在海底沿深度方向傳播過(guò)程中能量快速衰減,假設(shè)特定深度方向上聲波近似衰減為0是合理的,那么,海底截?cái)嗌疃仍酱?,?jì)算精度就越高,而計(jì)算速度卻越慢。為此,海底截?cái)嗌疃仁窃摲椒ㄐ枰攸c(diǎn)考慮的另一項(xiàng)重要參數(shù)。
圖9給出了海底截取深度分別為20 m和400 m,即1倍聲波波長(zhǎng)和20倍聲波波長(zhǎng)情況下剛性球球心與聲源所在垂直平面散射聲場(chǎng),其他計(jì)算參數(shù)與圖3相同。
對(duì)比圖9和圖3,虛擬柱島海底深度為1倍聲波波長(zhǎng)時(shí)目標(biāo)附近散射聲場(chǎng)計(jì)算結(jié)果與20倍聲波波長(zhǎng)間的偏差較大,當(dāng)海底深度增加至10倍聲波波長(zhǎng)時(shí)與20倍聲波波長(zhǎng)之間的偏差已經(jīng)難以察覺(jué)。
圖10給出了根據(jù)圖3、圖9虛擬柱島表面散射聲壓計(jì)算得到的目標(biāo)前向散射聲場(chǎng)隨深度的變化規(guī)律,其中,水平距離為15 m。
圖9 不同截止深度虛擬柱島情況下目標(biāo)附近散射聲壓級(jí)Fig.9 Scattering pressure levels in virtual column of different cut-off depths
圖10 不同截止深度虛擬柱島情況下目標(biāo)前向散射聲場(chǎng)隨深度變化規(guī)律Fig.10 Forward scattering fields in virtual columns of different cut-off depths with respect to depth
可以看到,當(dāng)海底截止深度取得過(guò)小,例如一倍聲波波長(zhǎng)時(shí),計(jì)算結(jié)果未收斂、偏差較大;當(dāng)海底截止深度需要足夠大時(shí),例如10 倍聲波波長(zhǎng),即衰減幅度閾值為5 dB,計(jì)算結(jié)果已經(jīng)收斂。
在仿真的計(jì)算時(shí)間上,海底截止深度為20 m、200 m 和400 m 時(shí)分別需要2 min、15 min 和20 min,海底截止深度越大,則計(jì)算時(shí)間越長(zhǎng),計(jì)算精度越高。在實(shí)際應(yīng)用中,可以根據(jù)需要權(quán)衡計(jì)算精度和計(jì)算速度,通常情況下海底截止深度衰減幅度閾值為5 dB即可。
構(gòu)建的模型如下:海水深度為30 m,海底采用吸收系數(shù)為2 dB/λ 的泥沙;目標(biāo)是半徑為5 m 的剛性球、深度為10 m;聲源是100 Hz 的單位幅度單極子聲源,位于距離目標(biāo)水平距離400 m、深度20 m。根據(jù)第2 章分析,虛擬柱島半徑、海底截止深度分別設(shè)置為30 m、37.5 m,即兩倍聲波波長(zhǎng)和幅度衰減5 dB所對(duì)應(yīng)的深度。
圖11給出了虛擬柱島內(nèi)剛性球球心與聲源所在的垂直平面、目標(biāo)表面以及虛擬柱島表面散射聲壓級(jí)。
圖11 淺海信道中近海面剛性球散射聲場(chǎng)Fig.11 Sound scattered by a rigid sphere near the sea surface in shallow waveguide
結(jié)合示意圖和偽彩圖可見(jiàn),對(duì)于目標(biāo)而言,目標(biāo)背向散射能量集中在斜下方,可以判斷出入射聲波是從左下方斜入射的,同時(shí)目標(biāo)前向散射能量也集中在斜下方,這是由于海面反射造成的。相較于上一個(gè)案例中水平入射,目標(biāo)散射聲場(chǎng)能量明顯向海底方向偏轉(zhuǎn)。由于目標(biāo)靠近界面,受海面多次散射的影響,目標(biāo)的散射聲壓由于疊加而增強(qiáng)。
除了計(jì)算標(biāo)準(zhǔn)體如球之外,本節(jié)以簡(jiǎn)化的Benchmark模型,來(lái)說(shuō)明該計(jì)算方法對(duì)于較為復(fù)雜目標(biāo)的應(yīng)用。
目標(biāo)的結(jié)構(gòu)主體是長(zhǎng)25 m、半徑2.5 m 的圓柱體,兩端是半徑為2.5 m 的兩個(gè)半球,組合而成的圓柱組合體模型,總長(zhǎng)度為30 m,深為10 m。目標(biāo)所在的環(huán)境參數(shù)如下:海水深度為30 m,海底取吸收系數(shù)為2 dB/λ 的泥沙。100 Hz的單位幅度單極子聲源位于水下5 m,目標(biāo)和聲源水平距離為400 m,目標(biāo)頭部指向聲源。將上述結(jié)構(gòu)用計(jì)算模型構(gòu)建可以得到圖12。
圖12 淺海信道中圓柱組合體模型Fig.12 Cylindrical composite model in shallow waveguide
圖13 淺海信道中近海面組合體散射聲場(chǎng)Fig.13 Sound scattered by a cylindrical composite model near the sea surface in shallow waveguide
通過(guò)仿真計(jì)算,可以得到虛擬柱島內(nèi)目標(biāo)中心與聲源所在垂直平面、目標(biāo)表面以及虛擬柱島表面散射聲壓級(jí),其中虛擬柱島半徑、深度參數(shù)與3.1節(jié)的相同。
可以看出,目標(biāo)的前向散射聲場(chǎng)能量更強(qiáng),表面亮點(diǎn)集中在目標(biāo)主體中后段、且強(qiáng)度是非均勻的。此外,由于目標(biāo)靠近界面,受界面的多次散射影響,目標(biāo)上、下方散射空間分布是不對(duì)稱(chēng)的。
本文提出了一種基于虛擬柱島、有限元法和簡(jiǎn)正波的淺海信道中目標(biāo)低頻聲散射的混合數(shù)值方法,該方法相對(duì)于其他方法,具有的特征為:(1)考慮了目標(biāo)和界面的多次散射,可以計(jì)算的目標(biāo)和界面之間的距離較近,即目標(biāo)靠近海面或海底的散射情況;(2)有限元和簡(jiǎn)正波混合方法比傳統(tǒng)有限元方法計(jì)算速度快、效率高。另外,該混合數(shù)值方法適用于低頻散射的情況,尤其是極低頻的條件,且適用于目標(biāo)在淺海信道中前向散射的情況。為了同時(shí)滿(mǎn)足計(jì)算精度與速度的需求,虛擬柱島半徑取一倍聲波波長(zhǎng),海底截止深度對(duì)應(yīng)幅度衰減5 dB。仿真計(jì)算了近海面典型目標(biāo)在淺海信道中的聲散射特征,對(duì)散射聲場(chǎng)方位偏轉(zhuǎn)以及表面亮點(diǎn)分布進(jìn)行了分析。該方法可應(yīng)用于淺海信道中復(fù)雜目標(biāo)回聲特征的仿真和評(píng)估。