大粗糙度起伏下一維海面聲散射建模與分析

周 江，王 斌，王文歡，范 軍

(1.上海交通大學，上海 200240；2.海洋智能裝備與系統教育部重點實驗室，上海 200240)

0 引 言

起伏海面的聲散射特性研究在水聲通信、主動聲吶探測和水聲對抗等領域有著廣泛的應用[1-5]。對于中遠距離的聲場傳播問題，大入射角甚至掠入射條件下的起伏海面聲散射特性一直是研究的重點之一。然而，對于水面目標或水面附近目標如水面艦船、水下航行器等的聲吶探測和預警問題，考慮到海面與目標之間復雜的相互作用，需要對起伏海面的聲散射進行建模與分析，尤其是小入射角下大粗糙度起伏海面的多次聲散射不容忽視。近年來國內外對粗糙海面的聲散射發展了多種精確和近似的解析方法，主要有：可以近似計算絕對軟邊界條件下粗糙海面散射聲場的積分方程法[6-7]，其優點是考慮了多次散射和影區的貢獻，但因其計算量較大無法應用于工程中；當隨機起伏方差和相關長度遠小于入射波波長時，適用于低頻近似條件的小尺度起伏粗糙面的微擾法[8-10]；當粗糙面平均曲率半徑遠大于入射波波長、隨機起伏方差和相關長度遠大于入射波波長時，適用于高頻近似條件下的大尺度起伏粗糙面的基爾霍夫近似(Kirchhoff Approximation,KA)方法[11-15]等，其優點是形式簡單，但是不能計算小入射角下大粗糙度起伏海面的多次聲散射；可以近似計算均方根斜率較小的粗糙面散射聲場的小斜率近似法[16-17]。隨著計算機技術的快速發展，一系列數值計算方法如有限元方法[18]、邊界元方法[19]等也可用于起伏海面聲散射計算，但是由于網格劃分工作量巨大，僅適用于低頻計算。針對KA方法的缺點，研究人員也做出了相應的修正，比如Mc-Cammon等[20-21]利用修正的KA方法計算了絕對軟邊界條件下正弦表面的散射場，同時也考慮影區、多次散射和表面曲率的影響。Welton[22]將KA方法與菲涅爾(Fresnel)相位近似相結合，研究了幾何影區、二階散射和發射(接收)指向函數對絕對軟邊界條件下正弦表面散射聲場的影響。盡管文獻[20-22]已經將海面多次散射引入到KA方法中，但其應用局限于確定性正弦海面聲散射情況。Richards等[23]將KA方法與射線理論相結合，根據本征聲線計算海面散射聲場，但是該方法不適用于焦散區和聲影區，難以計算大粗糙度起伏海面的多次散射聲場。

本文在傳統KA方法的基礎上考慮了海面的陰影區和亮區之間的多次散射，推導了大粗糙度起伏海面多次聲散射的迭代求解公式，從而提出了一種適用于起伏海面高頻散射聲場計算的改進KA方法。該方法拓寬了傳統KA方法的適用范圍，對海面的形式、均方根斜率和起伏高度以及聲波的入射角度沒有限制。本文以一維余弦和高斯譜海面的有限元結果為標準解，驗證了本文方法的準確性和適用性。以一維余弦海面為例，本文分析了不同入射聲波頻率和角度下海面散射強度的分布情況，并解釋了Bragg散射的產生機理。以一維高斯譜海面為例，本文討論了不同海面均方根高度和相關長度情況下高斯譜海面散射強度的變化規律。

1 改進的Kirchhoff近似方法

針對聲波入射到一維起伏海面上的情形，本文采用了如圖1所示的聲散射模型，海面S起伏高度為z=η(x)，聲波入射角和散射角分別為θi和θs，海面水平長度L取決于發射、接收指向性以及距離。利用改進的KA方法求解該模型的中高頻散射聲場，如圖2所示該方法將其聲散射分為兩部分，(1)入射聲波照射海面亮區的直接散射，忽略了海面影區對散射聲場的貢獻。(2)海面的多次聲散射，考慮了海面亮區內某點處的第一次聲散射對海面上其他點(包括亮區和影區)的作用，從而引起第二次聲散射，同樣也引起了其他高階散射。

圖1 起伏海面聲散射幾何模型Fig.1 Geometric model of acoustic scattering from undulating sea surface

圖2 起伏海面多次聲散射機理Fig.2 The mechanism of multiple acoustic scattering from the undulating sea surface

根據惠更斯第二積分公式[24]，可得滿足自由邊界條件的起伏海面散射聲壓(略去時間因子e-iωt)為

式中：是單位表面法向量，r2是接收點的位置矢徑，rs是起伏海面的位置矢徑。為了簡化書寫形式，將向量表示為rα-rβ?Rαβ，其模|rα-rβ|?Rαβ，單位向量RαβRαβ?，下角標α和β對應不同的位置矢徑。格林函數及其法向量可表示為

式中：cosθs=·(rs)，H0(1)和H1(1)分別為第一類零階和一階漢克爾函數。文獻[25]也指出，海面散射相當于一系列偶極子聲源的輻射，此時ptot(rs)是該等效偶極子模型的總聲源密度[25]，于是有：

式中：pq，q=1，2，…表示第q階偶極子聲源密度，式(3)可以通過迭代方法進行求解。

考慮到海面的第一次聲散射是由入射聲波激勵的，且海面表面滿足自由邊界條件，根據Kirchhoff近似可得

式中：入射聲壓pinc(rs)=H0(1)(kR1s)，r1是發射點矢徑；可見性函數γ1(R1s)定義為

考慮海面的第二次聲散射，此時海面上再輻射聲源是由第一階偶極子聲源激發的，于是有：

式中：rs，r's∈S分別為場點和源點，cosθs'=·(rs')；互可見性函數γ2(Rss')定義為

依次類推，該等效模型的第q(q≥2)階偶極子聲源密度可歸納為

為了對式(1)中的積分進行數值求解，將總長度為L的海面離散成N個單元，則該起伏海面的散射聲場近似等于這N個單元各階散射聲場相干疊加之和。以如圖3所示的線單元對為例，c(m)和c(n)是它們的幾何中心點，其矢徑分別為和，A(m)和A(n)是單元上任意點，其矢徑分別為和。為了簡化書寫，將向量表示為其模單位向量上角標(m)和(n)分別對應第m和n個線單元(m≠n)。式(1)可以離散表示為

圖3 線單元對等效偶極子模型Fig.3 Equivalent dipole models of line element pairs

第n個單元的第q(q≥2)階偶極子聲源密度為

此時，對于式(1)的求解就轉化成對N個線單元的積分計算，而線單元的積分式可以進一步轉化成與線單元坐標相關的求和解析式。利用一階泰勒公式進行近似，可以得到：

于是有：

該迭代過程的終止條件為

式中：ε為收斂限定的最小正數閾值。得到起伏海面的單位長度散射強度為

在上述的推導過程中，考慮了海面的影區遮擋和多次聲散射效應，故式(9)可應用于中高頻率、大粗糙度起伏海面的單、雙站聲散射計算。而經典的物理聲學方法通常假設幾何影區對散射聲場的貢獻可以忽略，故只限于式(9)中q=1時的一階散射聲場計算，式(9)僅討論了一維情況。將式(9)中二維格林函數改為三維格林函數，即可適用于二維海面的計算，從而更加符合工程實際應用情況。另外，下文中一維海面相關結論均可推廣到二維情況。由于本文側重點在于方法的提出和驗證，以及海面一般規律總結，所以對于不同參數下任意形狀二維海面的相關散射計算不再贅述。

2 方法驗證

起伏海面可以看作是由無數個不同波長、波高和相位的諧波海面疊加形成的。基于兩種不同類型的起伏海面模型對本文所提方法進行驗證，對于規則的余弦波海面，其波面方程為

式中：-L2≤x≤L2，H、Λ和φ分別為波高、波長和相位。對于隨機起伏海面，選擇典型的高斯譜描述粗糙海面的結構和能量分布，其一維表達式為[11]

式中：δ為海面的均方根高度，決定了粗糙海面的起伏程度；l為相關長度，它決定了粗糙海面的橫向變化頻率。再利用蒙特卡洛(Monte Carlo)方法[26]，結合高斯譜進行一維海面幾何建模。圖4(a)給出了一維余弦海面高度曲線η(x)，其中H=0.4 m，Λ=1 m，φ=π，L=10 m。圖4(b)給出了基于高斯譜生成的一維海面高度曲線η(x)，其中δ=0.4 m，l=1 m，L=10 m。

圖4 一維余弦和高斯譜海面模型Fig.4 One-dimensional cosine and Gaussian sea surface models

下面分別利用改進的KA方法、KA方法和有限元方法計算如圖4所示的兩種一維海面模型的散射強度。圖5比較了θi=10°和f=10 kHz時三種方法計算的散射強度隨散射角變化結果，其中圖5(a)為余弦海面模型，圖5(b)為高斯譜海面模型。圖6比較了θi=10°和θs=10°下的三種方法計算的海面鏡向散射強度隨頻率變化結果，其中圖6(a)為余弦海面模型，圖6(b)為高斯譜海面模型。

圖5 一維海面模型三種方法散射強度隨散射角變化計算結果Fig.5 Variations of scattering intensity with scattering angle calculated by three different methods for the two one-dimensional sea surface models

圖6 一維海面模型三種方法散射強度隨頻率變化計算結果Fig.6 Variations of scattering intensity with frequency calculated by three different methods for the two one-dimensional sea surface models

由圖5和圖6可以看出，對于大傾角的余弦和大粗糙度的隨機起伏海面，當入射角和散射角均較小時，傳統KA方法的計算結果與有限元結果誤差較大，且誤差隨著頻率的增大呈增大趨勢，而本文所提出的改進的KA方法很好地解決了傳統KA方法存在的問題，其計算結果與有限元結果具有很好的一致性，這為改進的KA方法計算起伏海面聲散射的準確性和適用性提供了理論依據。

3 計算結果與分析

3.1 不同入射波發射參數對余弦海面散射強度的影響

為了清晰描述不同參數對海面聲散射特性的影響，在以下的計算中均利用改進的KA方法計算式(14)定義的一維余弦海面的分置散射強度。若無特殊說明，海面參數選取H=0.4 mm，Λ=1 m，L=10 m。

圖7給出了f=10 kHz下四個典型聲波入射角度下海面的總體、一階和二階的散射強度對比結果，圖7(a)～7(e)分別為θi=0°，30°，45°和60°下的對比，圖7(f)是上述4個典型入射角度下總體、一階和二階的散射強度的概率分布結果。

圖7 不同入射角下余弦海面的一階、二階和總散射強度隨散射角的變化及其概率分布Fig.7 Variations of the first order, second order and total scattering intensities with scattering angle and their probability distributions under cosine sea surface and different incidence angles

從圖7(a)和7(c)可以看出，對該最大波傾角為72°的余弦海面，當入射角為0°或30°時，在部分散射角度下海面的二階散射強度比一階散射強度高出5～15 dB(在θi=0°時補充了放大圖)，在散射強度[0，15]dB區間內海面的總散射強度概率要大于一階散射強度概率，這說明海面的二次散射對總散射聲場的貢獻不可忽略；當入射角為45°和60°時，只有少數散射角度下海面的二階散射起作用；而當入射角為60°時，在散射強度[0，15]dB區間內海面的二階散射強度概率幾乎為0，這說明海面的高階散射貢獻可以忽略。進一步分析圖7可以看出，當入射角小于45°，且散射角|θs|小于45°時，余弦海面的二階散射聲場對總散射聲場的貢獻較大，而隨著入射角和散射角|θs|的減小，二階散射的影響是逐漸增大的。這說明在小入射角和小散射角條件下，海面的多次散射效應予以考慮，而且越接近垂直入射，海面的多次散射效應越明顯。

為了更清晰地觀察聲波入射角度對散射強度的影響，圖8進一步給出了f=10 kHz下海面散射強度的入射角度-散射角度譜圖。圖8(a)～8(c)分別對應總體、一階和二階的散射強度。圖9給出了θi=0°時余弦海面散射強度的頻率-散射角度分布圖，圖9(a)～9(c)分別對應總體、一階和二階的散射強度。

圖8 余弦海面散射強度隨入射角度和散射角度變化云圖Fig.8 Variation nephograms of the first order, second order and total scattering intensities with incidence angle and scattering angle under cosine sea surface

圖9 余弦海面散射強度隨頻率和散射角度變化云圖(θi=10°)Fig.9 Variation nephograms of the first order, second order and total scattering intensities with frequency and scattering angle under cosine sea surface and θi=10°

圖8(b)和圖9(b)中的亮線可以用Bragg散射來解釋。Bragg散射的根本原因是海面波浪間散射回波的相干疊加。Bragg散射示意圖如圖10所示。α和β是相鄰的兩個波谷，入射聲波通過點α的路徑長度為A+D，通過點β的路徑長度為D+E。只有通過兩點α和β的聲程差是聲波波長λ的整數倍，即滿足條件(D-E)-(B-A)=±nλ時，Bragg散射才會發生。根據幾何關系：(D-E)=Λsinθs和(B-A)=Λsinθi，Bragg散射的條件[27-28]為

圖10 Bragg散射示意圖Fig.10 Schematic diagram of Bragg scattering

圖8、9中黑色虛線均是根據式(16)畫出的Bragg散射線。由圖8(b)和圖9(b)可以看出，黑色線的位置與亮線吻合得非常好，且亮線的數量隨著頻率的提高而增多。這說明在中高頻段，當聲波波長遠小于海面波長時，Bragg散射對海面一階散射聲場的影響很大，而且隨著頻率的提高，Bragg散射效應更加明顯。在圖8(c)和圖9(c)中，另一個發現是，海面的二階聲散射亮線和一階聲散射亮線出現在相同的位置。

3.2 不同海面參數對高斯譜海面散射強度的影響

3.2.1 不同均方根高度下的高斯譜海面散射強度

圖11給出了不同波高下海面的散射強度隨散射角變化結果，其中θi=10°和f=10 kHz。保持高斯譜海面相關長度不變，改變均方根高度δ=0.4和1.2。隨著高斯譜海面均方根高度的增大，海面粗糙度增大，海面二階散射對總散射聲場的貢獻越來越大。當聲波以10°入射到均方根高度δ=1.2的高斯譜海面時，海面的多次散射無法忽略。在小入射角條件下，海面的均方根高度越大，海面的多次散射效應越明顯。

圖11 不同海面均方根高度對高斯譜海面散射強度的影響Fig.11 Influence of different root mean square heights of sea surface on the scattering intensity from Gaussian sea surface

3.2.2 不同相關長度下的高斯譜海面散射強度

圖12給出了不同波高下海面的散射強度隨散射角變化結果，其中θi=10°和f=10 kHz。保持海面均方根高度不變，改變相關長度l分別為0.5和1.5。隨著高斯譜海面相關長度的增大，海面二階散射對總散射聲場的貢獻越來越小，海面的一階散射對總散射聲場起主要貢獻。當聲波以10°入射到相關長度l=1.5的高斯譜海面時，海面的多次散射可以忽略。在小入射角條件下，海面的相關長度越小，海面的多次散射效應越明顯。

圖12 不同海面相關長度對高斯譜海面散射強度的影響Fig.12 Influence of different correlation lengths of sea surface on the scattering intensity from Gaussian sea surface

4 結 論

本文針對起伏海面高頻聲散射計算問題，提出了一種改進的物理聲學方法。與經典KA方法的不同之處在于，這種方法考慮了海面的陰影區和亮區之間的多次聲散射，較大地提高了計算精度。本文以一維余弦和一維高斯譜海面為例，該方法與有限元方法計算結果的比較，表明該方法適應于中高頻、小入射角下大傾角起伏海面的聲散射計算。根據一維余弦海面散射強度的計算結果，分析了不同入射聲波角度和頻率對海面散射強度的影響規律，解釋了Bragg散射的產生機理，同時討論了不同海面均方根高度和相關長度情況下高斯譜海面散射強度的變化規律。結果表明：

(1) 在小入射角和小散射角條件下，均方根高度較大的高斯起伏海面的多次散射效應要予以考慮。隨著入射角的減小和海面均方根高度的增大，海面的多次散射對總散射強度的貢獻越來越大；

(2) 在小入射角和小散射角條件下，相關長度較小的高斯起伏海面的多次散射效應要予以考慮。隨著入射角的減小和海面相關長度的減小，海面的多次散射對總散射強度的貢獻越來越大；

(3) 當聲波波長遠小于余弦海面波長時，Bragg散射對海面聲散射的影響較大。隨著余弦海面波長的增大和聲波頻率的提高，Bragg散射條紋增多，海面的Bragg散射效應更加明顯。