范 軍，湯渭霖，卓琳凱

（上海交通大學(xué) 海洋工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，上海200240）

1 引 言

聲吶工程中對(duì)于水中復(fù)雜形狀目標(biāo)回聲特性的預(yù)報(bào)目前主要采用兩種方法。一是基于亮點(diǎn)模型的部件法[1-2]，這種方法將復(fù)雜形狀的目標(biāo)分解為一組簡(jiǎn)單形狀的子目標(biāo)，每個(gè)子目標(biāo)的回聲用解析形式表示，計(jì)算簡(jiǎn)單且物理概念清晰。但由于子目標(biāo)的限制，對(duì)實(shí)際目標(biāo)形狀的逼近誤差較大；二是數(shù)值積分方法[3]，這種方法雖然能較精確地逼近復(fù)雜形狀的目標(biāo)，但由于要?jiǎng)澐值拿嬖獢?shù)量巨大（數(shù)萬(wàn)至數(shù)百萬(wàn)個(gè)），計(jì)算面積分的速度慢，且物理概念不清晰。

隨著聲吶技術(shù)和水下武器系統(tǒng)的發(fā)展，要求目標(biāo)回聲特性預(yù)報(bào)的精度更高、速度更快。特別是，由于水聲對(duì)抗和反對(duì)抗的發(fā)展，對(duì)于目標(biāo)回聲特性預(yù)報(bào)提出了一些新的需求：①隨著隱身技術(shù)的發(fā)展，聲吶目標(biāo)表面敷設(shè)吸聲覆蓋層，因此需要預(yù)報(bào)非剛性目標(biāo)的回聲；②隨著水中精確制導(dǎo)技術(shù)的發(fā)展，需要預(yù)報(bào)近距離目標(biāo)的回聲特性，而且要求給出二維和三維回聲圖像，計(jì)算的工作量大大增加，迫切要求提高計(jì)算速度。

根據(jù)這些需求，本文發(fā)展了水下聲吶目標(biāo)遠(yuǎn)場(chǎng)、近場(chǎng)回聲特性預(yù)報(bào)的一種數(shù)值計(jì)算方法-板塊元方法。這種方法在應(yīng)用物理聲學(xué)方法求解水中目標(biāo)散射聲場(chǎng)時(shí)，用一組平面板塊元近似目標(biāo)曲面，將所有板塊元的散射聲場(chǎng)疊加得到總散射聲場(chǎng)的近似值。由于板塊元方法把散射聲場(chǎng)的積分運(yùn)算轉(zhuǎn)化成代數(shù)運(yùn)算，使這種方法的計(jì)算速度大大提高。并把這種方法推廣到目標(biāo)近場(chǎng)回聲特性預(yù)報(bào)和敷設(shè)吸聲材料的非剛性表面的聲吶目標(biāo)回波特性預(yù)報(bào)中。

采用板塊元方法計(jì)算的關(guān)鍵是目標(biāo)幾何建模，即目標(biāo)表面板塊的劃分和獲取。建模的精度直接影響到計(jì)算結(jié)果的精度。用于目標(biāo)外形建模的方法很多，一是直接利用目標(biāo)的型值點(diǎn)采用數(shù)值曲面擬合方法得到目標(biāo)曲面，二是利用現(xiàn)成圖形軟件如3DS MAX、ANSYS、AutoCAD和I-DEAS等進(jìn)行幾何建模。本文對(duì)這兩種幾何建模方法也進(jìn)行了討論。

圖1 積分區(qū)域Fig.1 Integral domain

2 遠(yuǎn)場(chǎng)板塊元方法

首先歸納物理聲學(xué)或Kirchhoff近似方法。從散射問(wèn)題的Kirchhoff公式出發(fā)，如圖1所示目標(biāo)的散射聲場(chǎng)可以表示為：

其中s是散射體表面，n是外法線，r2是散射點(diǎn)矢徑。

為了避免解積分方程，對(duì)于高頻情況假設(shè)：

A.忽略幾何影區(qū)對(duì)聲場(chǎng)的貢獻(xiàn)：

實(shí)際積分面積是從M1，M2看去均處于亮區(qū)的那部分表面s0。

B.物體表面滿足局部平面波邊界條件[2]：

其中φi是入射波勢(shì)函數(shù)（省略e-iωt）,φi=A/r（1）eik0r1，V（α）是表面反射系數(shù)，Zn表面聲阻抗。由表面邊界條件可以得到：

對(duì)于收—發(fā)合置的情況，有

這就是計(jì)算散射聲場(chǎng)的Kirchhoff近似公式。與文獻(xiàn)[2]不同，（3）、（4）兩式適用于任意距離，不管是近場(chǎng)還是遠(yuǎn)場(chǎng)。

首先在遠(yuǎn)場(chǎng)情況下建立板塊元方法。由（4）式取k0r?1得到遠(yuǎn)場(chǎng)條件下的聲吶目標(biāo)強(qiáng)度：

其中：

圖2 平面多邊形Fig.2 The planar polygon

圖3 目標(biāo)表面網(wǎng)格Fig.3 Target surface meshes

復(fù)雜聲吶目標(biāo)的目標(biāo)強(qiáng)度可以按下面的板塊元方法計(jì)算。先把目標(biāo)表面劃分為如圖3的許多小板塊網(wǎng)格。對(duì)亮區(qū)內(nèi)所有板塊元的散射聲場(chǎng)求和就得到目標(biāo)散射聲場(chǎng)的近似值。假如目標(biāo)表面劃分為N×M個(gè)網(wǎng)格，則可以得到一系列板塊元Si,j。則（6）式可以表示成：

其中，αij是第（i，j）個(gè)板塊元的法向與入射和反射方向的夾角。由于直接構(gòu)造的每個(gè)板塊元的取向不同，它們不可能在同一平面中。因此在應(yīng)用（7）式計(jì)算板塊元的積分前先要將不同取向的板塊元通過(guò)坐標(biāo)變換統(tǒng)一變換到某個(gè)確定的平面上。再由（7）式得到：

其中xn′，yn′是變換到二維平面后第 （i,j）個(gè)板塊的頂點(diǎn)坐標(biāo)，V（ αij′）是第 （i,j）個(gè)板塊的的局部表面反射系數(shù)，αij′是第（i,j）個(gè)板塊的法線和入射聲線夾角。在后面的具體計(jì)算中，取三角形板塊元，K=3。

3 近場(chǎng)板塊元方法

對(duì)于近場(chǎng)問(wèn)題，公式（3）、（4）中被積函數(shù)的距離因子不能直接提到積分號(hào)外，相位因子也不能按平行聲線計(jì)算。因?yàn)樵诮鼒?chǎng)中散射表面的不同部位到發(fā)射或接收點(diǎn)的距離之差可能與波長(zhǎng)相比擬，聲線不能看作平行的。這就使得（3）、（4）式的計(jì)算較復(fù)雜。這里采用如下的方法解決，雖然對(duì)于目標(biāo)整體而言，接收換能器處在近場(chǎng)，但若把目標(biāo)劃分為足夠小的小板塊，則相對(duì)于每一個(gè)小板塊，接收換能器處在遠(yuǎn)場(chǎng)?？偟纳⑸渎晥?chǎng)仍然可以表示成這些小板塊散射聲場(chǎng)的和。這樣就要求每個(gè)小板塊的尺寸滿足條件：Rmin＞D2/λ。其中Rmin是需要計(jì)算的最小距離，D是小板塊的最大尺寸，λ是入射聲波波長(zhǎng)。這樣（3）、（4）式可以進(jìn)一步簡(jiǎn)化。

圖4 板塊近場(chǎng)散射Fig.4 The near scattering from planar

考慮如圖4所示的第（i,j）個(gè)板塊，o點(diǎn)是坐標(biāo)原點(diǎn)，一般取為目標(biāo)中心。c是在板塊上選取的參考點(diǎn)，矢徑為，一般取為板塊的幾何中心處。Q點(diǎn)是板塊上任意點(diǎn)，Q點(diǎn)的矢徑為 ，發(fā)射換能器的矢徑為，接收換能器的矢徑為。令:

由于板塊較小，可以忽略因發(fā)射指向性造成的入射聲場(chǎng)在板塊上的不均勻性。設(shè)發(fā)射指向性為DT10），接收指向性為DR），接收到的散射聲波是：

對(duì)于收發(fā)合置情況則有：

其中的積分項(xiàng)與遠(yuǎn)場(chǎng)板塊元方法中的積分項(xiàng)（6）是一樣的，也可以將它轉(zhuǎn)化為平面內(nèi)關(guān)于多邊形各頂點(diǎn)坐標(biāo)的和。

其中P是目標(biāo)上能夠被入射波照射到的板塊數(shù)，它們滿足

從上面的推導(dǎo)過(guò)程可以看到，遠(yuǎn)場(chǎng)板塊元方法之所以能夠推廣到近場(chǎng)是因?yàn)檫@時(shí)雖然對(duì)整個(gè)目標(biāo)來(lái)說(shuō)測(cè)量點(diǎn)是近場(chǎng)但對(duì)于每個(gè)小板塊元卻是遠(yuǎn)場(chǎng)，因此可以應(yīng)用遠(yuǎn)場(chǎng)近似，即入射聲線可以近似為平行聲線，板塊元積分中分母上的r1，r2可以用參考距離R1，R2近似并移出積分號(hào)，最終的積分與遠(yuǎn)場(chǎng)的一樣，可以利用積分化為代數(shù)和的獨(dú)特優(yōu)點(diǎn)。

4 幾何建模

不論是板塊元計(jì)算，還是小板塊的構(gòu)造，都需要知道目標(biāo)表面的空間坐標(biāo)，這就要對(duì)目標(biāo)表面幾何進(jìn)行數(shù)值建模。目標(biāo)表面幾何建模的精度直接關(guān)系到板塊元方法預(yù)報(bào)目標(biāo)回聲特性的準(zhǔn)確性，因此是板塊元方法的關(guān)鍵所在。用于目標(biāo)外形建模的方法很多，一是直接利用目標(biāo)的型值點(diǎn)采用數(shù)值曲面擬合方法得到目標(biāo)曲面。在計(jì)算機(jī)輔助幾何外形設(shè)計(jì)（CAGD）中，適用于曲面造型的主要曲面有[6]：參數(shù)樣條曲面、Bezier曲面和B樣條曲面等，但都存在求曲面控制點(diǎn)和多曲面過(guò)渡拼接的問(wèn)題。使這種方法較為復(fù)雜和麻煩，可操作性差，二是利用現(xiàn)成圖形軟件進(jìn)行幾何建模。隨著計(jì)算機(jī)圖形技術(shù)的發(fā)展，出現(xiàn)了許多現(xiàn)成的幾何建模軟件，如3DS MAX、ANSYS、AutoCAD和I-DEAS等。利用現(xiàn)有軟件直接進(jìn)行幾何建模，可以使建模工作大大簡(jiǎn)化。特別是目前在水下目標(biāo)，如潛艇、魚(yú)雷等的設(shè)計(jì)階段大多采用了一些圖形軟件，這樣可以直接利用設(shè)計(jì)結(jié)果對(duì)目標(biāo)的回聲特性進(jìn)行理論設(shè)計(jì)和預(yù)報(bào)，可以大大節(jié)省設(shè)計(jì)費(fèi)用，提高設(shè)計(jì)效率。這種方法非常適合在實(shí)際工程中應(yīng)用。

5 方法驗(yàn)證的數(shù)值算例

5.1 遠(yuǎn)場(chǎng)板塊元方法的驗(yàn)證

為考核板塊元方法的計(jì)算精度和速度，采用解析方法、數(shù)值積分方法和遠(yuǎn)場(chǎng)板塊元方法分別計(jì)算如圖5的橢圓錐臺(tái)隨方位角φ變化的目標(biāo)強(qiáng)度。計(jì)算參數(shù)如下：a=0.1 m，b=0.02 m，β=10°，L-L1=h=0.12 m，入射平面波頻率f=1 MHz，入射角θ=90°。橢圓錐臺(tái)目標(biāo)強(qiáng)度解析表達(dá)式見(jiàn)附錄1。直接數(shù)值積分方法和板塊元方法在劃分不同網(wǎng)格數(shù)情況下計(jì)算所需時(shí)間如表1所示，計(jì)算結(jié)果如圖6和圖7所示。

從計(jì)算結(jié)果來(lái)看可以得到以下結(jié)論：

（1）當(dāng)劃分的網(wǎng)格數(shù)200×300時(shí)，兩種方法計(jì)算結(jié)果和解析解比較都較滿意，數(shù)值積分的計(jì)算速度大于板塊元方法。當(dāng)劃分的網(wǎng)格數(shù)50×100時(shí)，數(shù)值積分的計(jì)算速度大于板塊元方法，板塊元方法得到的結(jié)果精度較好，但直接數(shù)值積分的結(jié)果精度很差。

圖5 橢圓錐臺(tái)Fig.5 The truncated elliptical cone

（2）雖然網(wǎng)格劃分?jǐn)?shù)相同，計(jì)算目標(biāo)強(qiáng)度時(shí)直接積分的速度大于板塊元方法，但是板塊元方法可以用較少的網(wǎng)格劃分得到較滿意的結(jié)果，這樣所用的計(jì)算時(shí)間減少，相當(dāng)于提高了速度，大約提高了10倍。

表1 計(jì)算所需時(shí)間Tab.1 Computing time

圖6 橢圓錐臺(tái)的目標(biāo)強(qiáng)度，網(wǎng)格：200×300Fig.6 TS of the truncated elliptical cone Meshes:200×300

圖7 橢圓錐臺(tái)的目標(biāo)強(qiáng)度，網(wǎng)格：50×100Fig.7 TS of the truncated elliptical cone Meshes:50×100

5.2 近場(chǎng)板塊元方法的驗(yàn)證

為了檢驗(yàn)近場(chǎng)板塊元算法的準(zhǔn)確性，用典型問(wèn)題—半徑a=1 m的剛性球散射聲場(chǎng)的計(jì)算結(jié)果來(lái)考核。由于目標(biāo)強(qiáng)度TS是一個(gè)遠(yuǎn)場(chǎng)聲吶參數(shù)，對(duì)于近場(chǎng)問(wèn)題不適用，我們用回聲強(qiáng)度ES（Echo Strength）來(lái)描述目標(biāo)近場(chǎng)特性：

其中pi是入射波在目標(biāo)中心點(diǎn)的聲壓，ps是接收點(diǎn)的聲壓，r0是從中心點(diǎn)到接收點(diǎn)的距離。分別采用嚴(yán)格的Rayleigh簡(jiǎn)正級(jí)數(shù)解、Kirchhoff近似解析解和近場(chǎng)板塊元方法對(duì)剛性球的近場(chǎng)回聲強(qiáng)度進(jìn)行計(jì)算和對(duì)比。其中嚴(yán)格的Rayleigh簡(jiǎn)正級(jí)數(shù)解、Kirchhoff近似解析解見(jiàn)附錄B，C。計(jì)算得到回聲強(qiáng)度ES隨頻率的變化曲線，如圖8所示。從圖中可以看出：

（1）Kirchhoff近似解析方法得到的結(jié)果和用近場(chǎng)板塊元方法得到的結(jié)果完全吻合。說(shuō)明，近場(chǎng)板塊元方法作為Kirchhoff近似的數(shù)值方法有足夠高的精度。

（2）Kirchhoff近似解析解、近場(chǎng)板塊元方法得到的結(jié)果同Rayleigh簡(jiǎn)正級(jí)數(shù)解的計(jì)算結(jié)果相比，隨頻率的升高三者逐漸趨于一致。當(dāng)f>2 kHz時(shí)，三者的差在1dB以?xún)?nèi)起伏。

圖8 近場(chǎng)球回聲強(qiáng)度（ES）r=10 m，a=1 mFig.8 The near field echo strength(ES)r=10 m,a=1 m

6 板塊元方法在水下目標(biāo)回聲特性預(yù)報(bào)中的應(yīng)用

6.1 沉底水雷目標(biāo)回波特性預(yù)報(bào)[7]

利用界面附近目標(biāo)回波計(jì)算的Kirchhoff近似采用3DS MAX對(duì)水雷目標(biāo)進(jìn)行建模。目標(biāo)由半球和圓柱合成。球半徑0.266 m，柱半徑0.266 m，柱長(zhǎng)1.599 m，該目標(biāo)下方有一水平淤泥界面。計(jì)算頻率為5 000 Hz，入射角為θ，建模后共有2 702個(gè)頂點(diǎn)和5 325個(gè)板塊。建模結(jié)果如圖9所示，有界面和無(wú)界面情況目標(biāo)強(qiáng)度見(jiàn)圖10。

圖9 近似沉底水雷型目標(biāo)Fig.9 Approximate mines lying on the seabed

圖10 水雷目標(biāo)在無(wú)界面和淤泥界面情況下目標(biāo)強(qiáng)度 —淤泥海底情況 －o－無(wú)海底情況Fig.10 TS of approximate mines lying on the seabed—muddy seabed －o－no seabed

6.2 非剛性目標(biāo)回波特性預(yù)報(bào)[8]

板塊元方法不僅可以預(yù)報(bào)剛性目標(biāo)回波，而且也可以預(yù)報(bào)非剛性目標(biāo)，如單層或雙層殼體目標(biāo)的回波。如果已知敷瓦殼體的整體反射系數(shù)，也可以預(yù)報(bào)敷瓦目標(biāo)的回波[8]。作為例子，分別用修正亮點(diǎn)模型[8]和遠(yuǎn)場(chǎng)板塊元方法計(jì)算了如圖5所示橢圓錐臺(tái)殼體裸露情況和敷設(shè)消聲瓦時(shí)的目標(biāo)強(qiáng)度，如圖11所示。由圖中可見(jiàn)兩者計(jì)算結(jié)果吻合很好。

圖11 橢圓錐臺(tái)殼體目標(biāo)強(qiáng)度（TS）修正幾何亮點(diǎn)模型和板塊元方法計(jì)算結(jié)果Fig.11 Calculation results of target strength(TS)values of the truncated elliptical cone shell using the modified geometrical highlight model and calculation results using the planar elements method

6.3 水下目標(biāo)近場(chǎng)回波特性研究[9]

利用近場(chǎng)板塊元方法初步研究了簡(jiǎn)單目標(biāo)近場(chǎng)回波特性。橢圓柱體的回聲強(qiáng)度ES隨距離r的變化規(guī)律如圖12所示。橢圓的長(zhǎng)半軸為4 m，短半軸為1.5 m，橢圓柱的高度為4 m。計(jì)算的頻率為5 kHz，10 kHz。聲波垂直于橢圓柱軸線沿短半軸方向入射。從計(jì)算結(jié)果可以看出橢圓柱的ES隨距離變化特征與經(jīng)典活塞面散射聲場(chǎng)的特征[10]相似，存在明顯的近場(chǎng)區(qū)即費(fèi)涅爾（Fresnel）干涉區(qū)和遠(yuǎn)場(chǎng)區(qū)即弗郎霍夫爾（Fraunhofer）區(qū)，越靠近橢圓柱表面，干涉現(xiàn)象越明顯；隨著距離的增加，干涉峰值間隔增大；當(dāng)距橢圓柱幾何中心一定的距離時(shí)，不再出現(xiàn)干涉現(xiàn)象，ES值表現(xiàn)為單調(diào)，聲場(chǎng)逐漸由近場(chǎng)區(qū)進(jìn)入遠(yuǎn)場(chǎng)區(qū)，這與物理概念是符合的。

圖12 橢圓柱回聲強(qiáng)度（ES）隨距離變化Fig.12 Echo strength of Ellipsoid versus distance

7 結(jié) 論

本文建立了水下目標(biāo)遠(yuǎn)場(chǎng)、近場(chǎng)回聲特性預(yù)報(bào)的一種數(shù)值計(jì)算方法—板塊元方法。這種方法在應(yīng)用物理聲學(xué)方法計(jì)算水中目標(biāo)散射聲場(chǎng)時(shí)，用一組平面板塊元近似目標(biāo)曲面，求出板塊元的散射聲場(chǎng)矢量和作為散射聲場(chǎng)的近似值，得到目標(biāo)的回聲特性。板塊元方法可以計(jì)算剛性、表面敷設(shè)粘彈性材料的目標(biāo)遠(yuǎn)場(chǎng)和近場(chǎng)的目標(biāo)回聲特性。與直接數(shù)值積分相比，這種方法具有運(yùn)算速度快，精度高的特點(diǎn)。在水下目標(biāo)回波特性預(yù)報(bào)的許多方面具有廣闊的應(yīng)用前景。

[1]湯渭霖.聲納目標(biāo)回波的亮點(diǎn)模型[J].聲學(xué)學(xué)報(bào),1994,19(2):92-100.

[2]湯渭霖.用物理聲學(xué)方法計(jì)算非硬表面的聲散射[J].聲學(xué)學(xué)報(bào),1993,18(1):45－53.

[3]范 軍.復(fù)雜目標(biāo)回波特性研究[D].上海:上海交通大學(xué),2002.

[4]周建江,舒永澤.用板塊法近似計(jì)算任意復(fù)雜形體目標(biāo)的雷達(dá)散射截面[J].電子科學(xué)學(xué)刊,1992,14(1):71-75.

[5]蘇東林,王寶法,張采來(lái).復(fù)雜目標(biāo)外形擬合及散射截面計(jì)算[J].系統(tǒng)工程與電子技術(shù),1990(12):18-24.

[6]李振格.AUTO CAD用戶參考手冊(cè)[M].北京:海洋出版社,1991.

[7]萬(wàn) 琳,范 軍,湯渭霖.沉底水雷目標(biāo)強(qiáng)度與回聲信混比[J].聲學(xué)學(xué)報(bào),2003,28(5):429-433.

[8]范 軍,湯渭霖.非剛性表面聲吶目標(biāo)回波的修正幾何亮點(diǎn)模型[J].聲學(xué)學(xué)報(bào),2001,26(6):545-550.

[9]李建魯,范 軍,湯渭霖.水下簡(jiǎn)單形狀目標(biāo)回聲的近、遠(yuǎn)場(chǎng)過(guò)渡特性[J].上海交通大學(xué)學(xué)報(bào),2001,35(12):1846-1850.

[10]何祚鏞.聲學(xué)理論基礎(chǔ)[M].北京:國(guó)防工業(yè)出版社,1981.

附錄A 橢圓錐臺(tái)解析表達(dá)式

附錄B 近場(chǎng)剛性球體嚴(yán)格Rayleigh簡(jiǎn)正級(jí)數(shù)解

球的散射波勢(shì)函數(shù)可以設(shè)為

利用剛性邊界條件：

圖B.1剛性球近場(chǎng)散射

Fig.B.1 Near scattering from rigid sphere

得到：

收發(fā)合置情況下，r=r0，θ=π，有：

按定義

附錄C 近場(chǎng)剛性球體Kirchhoff近似解

取坐標(biāo)系如圖C.1，設(shè)觀察點(diǎn)M與球心的距離是r0，面元與觀察點(diǎn)的距離是r。球坐標(biāo)系的極角是θ，球面上面元的法線與入射聲線的夾角是Θ。由三角關(guān)系式得到：從中解出，

圖C.1剛性球的近場(chǎng)散射Fig.C.1 Near scattering from rigid sphere

將以上結(jié)果代入（4）式得到：面積分變成了對(duì)θ的一維積分，積分限為0～θm，θm=cos-1（a/ r0）。此式難以導(dǎo)出解析形式。這里討論一種有實(shí)際意義的近似，設(shè)目標(biāo)處在不十分近的近場(chǎng)中，也就是目標(biāo)上任何一個(gè)面元到觀察點(diǎn)的距離r滿足k0r>>1。在算例中，球的半徑a=1 m，聲波頻率大于1 kHz，只要r0＞3.5 m就有k0r＞10。在這種情況下（C.2）式可簡(jiǎn)化成

將（C.1）式代入并令ξ=r0/a為無(wú)因次的距離變量，得到：