李承幫，江鵬飛，孫軍平，衣雪娟，林建恒

(1.中國(guó)科學(xué)院聲學(xué)研究所北海研究站，山東青島 266114；2.中國(guó)科學(xué)院大學(xué)，北京 100049)

0 引 言

波數(shù)積分法是一種常用的水下聲場(chǎng)計(jì)算方法，此方法對(duì)波數(shù)域格林函數(shù)進(jìn)行直接的數(shù)值積分，因而被認(rèn)為是一種精確的聲場(chǎng)計(jì)算方法。波數(shù)積分法由Pekeris引入到水下聲場(chǎng)的計(jì)算中[1]，此后很多學(xué)者對(duì)其具體實(shí)現(xiàn)方法進(jìn)行了研究。DiNapoli 等[2]基于波數(shù)積分理論提出了聲場(chǎng)計(jì)算效率非常高的快速場(chǎng)程序(Fast Field Program)，Schmidt[3-4]提出了直接全局矩陣法實(shí)現(xiàn)聲場(chǎng)的波數(shù)積分計(jì)算，駱文于等[5-6]提出了一種可穩(wěn)定計(jì)算Pekeris 波導(dǎo)中聲場(chǎng)的波數(shù)積分方法。

波數(shù)積分法的一個(gè)重要步驟是對(duì)波數(shù)域格林函數(shù)進(jìn)行逆漢克爾變換，此處會(huì)面臨半無限區(qū)間的積分問題。一般來說，當(dāng)水平波數(shù)大于某一個(gè)臨界值時(shí)，波數(shù)域格林函數(shù)的幅值會(huì)隨著水平波數(shù)的增大而衰減，并趨向于0。因此通常的做法是，當(dāng)波數(shù)域格林函數(shù)幅值衰減到一定程度時(shí)，對(duì)波數(shù)積分區(qū)間進(jìn)行截?cái)?，忽略大于此波?shù)積分區(qū)間的聲場(chǎng)。若截?cái)嗖〝?shù)取值太大，會(huì)增加不必要的運(yùn)算量；若截?cái)嗖〝?shù)取值太小，則會(huì)給聲場(chǎng)的計(jì)算帶來誤差。積分截?cái)嘞喈?dāng)于給波數(shù)域格林函數(shù)加矩形窗，在截?cái)嗵幰肓瞬贿B續(xù)點(diǎn)，導(dǎo)致出現(xiàn)截?cái)嘈?yīng)[7]。然而頻率、介質(zhì)類型、波導(dǎo)參數(shù)、聲源與接收點(diǎn)的相對(duì)深度等參數(shù)會(huì)對(duì)波數(shù)域格林函數(shù)的衰減特性產(chǎn)生很大的影響，在計(jì)算出波數(shù)域格林函數(shù)之前很難對(duì)截?cái)嗖〝?shù)進(jìn)行合適的選取。

在不同介質(zhì)類型的海洋波導(dǎo)中，波數(shù)域格林函數(shù)會(huì)有不同的衰減特性。對(duì)于不含冰蓋層且海底為液態(tài)的海洋波導(dǎo)，波數(shù)域格林函數(shù)的衰減特性比較簡(jiǎn)單，當(dāng)水平波數(shù)大于臨界值時(shí)，波數(shù)域格林函數(shù)不存在極點(diǎn)；而對(duì)于存在彈性海底或冰蓋層的海洋波導(dǎo)，在波數(shù)域格林函數(shù)衰減較慢，在水平波數(shù)大于臨界值的較大范圍內(nèi)，很遠(yuǎn)的位置仍包含聲場(chǎng)的重要成分[7]。本文主要關(guān)注的是各層均為流體的海洋波導(dǎo)中，寬帶聲場(chǎng)波數(shù)積分計(jì)算的截?cái)嗖〝?shù)選取問題。

當(dāng)海洋波導(dǎo)中各層介質(zhì)都是流體時(shí)，目前最常用的截?cái)嗖〝?shù)選取方法是把波導(dǎo)內(nèi)最小聲速對(duì)應(yīng)的波數(shù)k0=ω/cmin作為一個(gè)參考波數(shù)，截?cái)嗖〝?shù)取值為kt=ξ·k0，其中ξ取一個(gè)固定的值。文獻(xiàn)[8]認(rèn)為ξ可取1.1；文獻(xiàn)[9]認(rèn)為大多數(shù)情況下ξ取1.1～1.2 時(shí)即可包含波數(shù)域格林函數(shù)中絕大部分能量。在寬帶聲場(chǎng)計(jì)算過程中，若ξ取固定值，則不能兼顧低頻段的積分精確度和高頻段的計(jì)算量。文獻(xiàn)[10]提出一種消除截?cái)嘈?yīng)的方法，在波數(shù)域格林函數(shù)截?cái)帱c(diǎn)后，添加埃爾米特(Hermite)多項(xiàng)式函數(shù)，消除截?cái)嗵幍牟贿B續(xù)性并使其迅速衰減。該方法雖然能避免出現(xiàn)截?cái)嘈?yīng)，但是會(huì)給聲場(chǎng)計(jì)算帶來一定的誤差。文獻(xiàn)[11]提出了基于預(yù)估-校正思想的最大截止波數(shù)自動(dòng)選取算法。但此方法需要在波數(shù)域格林函數(shù)計(jì)算過程中判斷是否進(jìn)行積分截?cái)?，這不利于實(shí)現(xiàn)波數(shù)域格林函數(shù)的并行計(jì)算。

針對(duì)上述問題，本文提出一種應(yīng)用于寬帶點(diǎn)聲源聲場(chǎng)計(jì)算的積分截?cái)嗖〝?shù)自適應(yīng)選取方法。首先構(gòu)造一個(gè)數(shù)學(xué)模型擬合水平波數(shù)大于參考波數(shù)時(shí)的波數(shù)域格林函數(shù)，然后利用卡爾曼濾波器對(duì)模型的參數(shù)進(jìn)行跟蹤和預(yù)測(cè)，再根據(jù)已計(jì)算頻率的模型參數(shù)估計(jì)出下一頻率的模型參數(shù)，最后將估計(jì)的模型參數(shù)代入數(shù)學(xué)模型求解出下一個(gè)頻率的截?cái)嗖〝?shù)，從而實(shí)現(xiàn)給定精確度下的積分截?cái)唷?/p>

1 截?cái)嗖〝?shù)自適應(yīng)選取方法

在利用波數(shù)積分法計(jì)算頻率為f的點(diǎn)聲源產(chǎn)生的聲場(chǎng)時(shí)，需要對(duì)波數(shù)域格林函數(shù)G(kr;f)進(jìn)行逆漢克爾變換：

式中：J0(krr)為零階貝塞爾函數(shù)，r為水平距離，kr為水平波數(shù)。當(dāng)kr大于某個(gè)臨界值時(shí)，G(kr;f)隨kr的增大而衰減，在G(kr;f)明顯衰減后進(jìn)行積分截?cái)嗖粫?huì)帶來很大誤差。本文將積分截?cái)嗖〝?shù)記為kt，kt的取值取決于波數(shù)域格林函數(shù)的衰減特性，即G(kr;f)開始衰減的位置以及衰減的速率。

考慮頻率為fn=f0+n·Δf的寬帶點(diǎn)聲源聲場(chǎng)的計(jì)算問題，其中n=1，…，N，Δf是頻率間隔。由1.1節(jié)的算例可看出，隨頻率fn的變化，G(kr;fn)呈現(xiàn)出一定的變化規(guī)律。因而可根據(jù)已計(jì)算的G(kr;f0)，G(kr;f1)，…，G(kr;fn)的衰減特性，估計(jì)G(kr;fn+1)的衰減特性。但由于G(kr;f)是數(shù)值計(jì)算的結(jié)果，為了便于跟蹤和預(yù)測(cè)它的衰減特性，需要構(gòu)造一個(gè)合適的數(shù)學(xué)模型對(duì)其進(jìn)行擬合，并利用自適應(yīng)濾波器對(duì)此數(shù)學(xué)模型的擬合參數(shù)進(jìn)行跟蹤和預(yù)測(cè)。根據(jù)模型參數(shù)的預(yù)測(cè)結(jié)果，可計(jì)算出下一個(gè)頻率滿足精度要求的截?cái)嗖〝?shù)。

1.1 波數(shù)域格林函數(shù)的衰減特性

本節(jié)分析波數(shù)域格林函數(shù)的衰減特性，以選取合適的數(shù)學(xué)模型對(duì)其進(jìn)行擬合。在Pekeris波導(dǎo)中，水層的波數(shù)域格林函數(shù)G1(kr;f)可表示為聲源項(xiàng)與上行波、下行波的疊加：

式中：A+1(kr)、A-1(kr)分別表示上行波、下行波的波譜數(shù)；kz，1=k1是水層的波數(shù)；z是接收深度，zs是聲源深度；Sω是聲源強(qiáng)度。文獻(xiàn)[12]中指出當(dāng)水平波數(shù)kr比較大時(shí)，在包含聲源的介質(zhì)層中，式(2)中的聲源項(xiàng)具有最慢的衰減速率，因而G1(kr;f)的衰減速率主要由聲源項(xiàng)決定。

以下討論G1(kr;f)中聲源項(xiàng)的衰減特性。用?s表示式(2)中聲源項(xiàng)的對(duì)數(shù)幅值，表達(dá)式為

令ξ=kr/k1，則kr=ξ·k1。當(dāng)kr＞k1時(shí)，ξ＞1，且把上述kz，1和kr代入式(3)中并化簡(jiǎn)可得：

當(dāng)ξ＞1時(shí)，將其代入式(4)，可得：

當(dāng)ξ→1 時(shí)，?s(ξ;f) →∞；當(dāng)ξ?1 時(shí)，ξ≈將其代入式(4)中，可得：

在Pekeris 波導(dǎo)中，由式(5)和式(6)可知當(dāng)ξ＞1時(shí)，聲源項(xiàng)的對(duì)數(shù)幅值?s(ξ;f)隨ξ的變化關(guān)系近似滿足一次線性項(xiàng)與對(duì)數(shù)項(xiàng)的疊加。

在海底為液態(tài)且沒有冰蓋層的海洋波導(dǎo)中，波數(shù)域格林函數(shù)G(kr;f)也會(huì)呈現(xiàn)與Pekeris波導(dǎo)類似的衰減特性。令參考波數(shù)k0(f)為水平分層介質(zhì)中最小聲速對(duì)應(yīng)的波數(shù)：

其中：cm表示第m層介質(zhì)中的聲速。當(dāng)kr＞k0(f)時(shí)，任意一層介質(zhì)的垂直波數(shù)都為純虛數(shù)。與Pekeris 波導(dǎo)中的情況類似，當(dāng)kr＞k0(f)時(shí)，G(kr;f)會(huì)隨kr的增大而減小，并趨向于0。對(duì)G(kr;f)的幅值取對(duì)數(shù)，并且把水平波數(shù)kr與參考波數(shù)k0(f)的比ξ作為自變量，令：

其中，水平波數(shù)比ξ的表達(dá)式為

?(ξ;f)能較好地反映G(kr;f)的衰減特性。以圖1中的波導(dǎo)環(huán)境為例，首先計(jì)算出不同頻率的波數(shù)域格林函數(shù)G(kr;f)，然后計(jì)算出對(duì)應(yīng)的?(ξ;f)，計(jì)算結(jié)果如圖2 所示。可以看出當(dāng)ξ≥1，即kr≥k0(f)時(shí)，G(kr;f)開始衰減。頻率越高，G(kr;f)衰減得越快。

圖1 海洋波導(dǎo)環(huán)境特例Fig.1 A special case of marine waveguide environment

圖2 不同頻率波數(shù)域格林函數(shù)對(duì)數(shù)幅值Fig.2 Logarithmic magnitude of Green's function in wavenumber domain at different frequencies

設(shè)期望的積分截?cái)嗑_度參數(shù)為a，表示期望在?(ξ;f) =-a處進(jìn)行積分截?cái)?。?dāng)a分別取3，4，5，6 時(shí)，計(jì)算在期望截?cái)嗵帵蝨的取值隨頻率變化，結(jié)果如圖3所示。結(jié)合圖2的分析結(jié)果，可以看出在較低頻段，截?cái)嗖〝?shù)比ξt需要取遠(yuǎn)大于1的值才能獲得足夠的積分精確度。其原因是在低頻段?(ξ;f)衰減得比較慢，而在較高的頻段?(ξ;f)衰減得比較快，截?cái)嗖〝?shù)比ξt取稍大于1的值就能獲得足夠的積分精確度。

圖3 截?cái)嗖〝?shù)與參考波數(shù)的比隨頻率變化Fig.3 The ratio of truncated wavenumber to reference wavenumber versus frequency

目前常用的截?cái)嗖〝?shù)選取方法是令ξt等于某個(gè)固定的值。在寬帶聲場(chǎng)計(jì)算中，若ξt取一個(gè)較大的值，可保證在低頻聲場(chǎng)計(jì)算時(shí)獲得足夠的積分精確度，但會(huì)給高頻聲場(chǎng)計(jì)算增加不必要運(yùn)算量；若ξt取定一個(gè)較小的值，可減少高頻聲場(chǎng)的計(jì)算量，但會(huì)增大低頻聲場(chǎng)計(jì)算誤差。

1.2 波數(shù)域格林函數(shù)的模型擬合

由1.1節(jié)的分析，在海底為液態(tài)且沒有冰蓋層的海洋波導(dǎo)中，當(dāng)kr≥k0(f)時(shí)，G(kr;f)開始衰減并趨向于0。在保證精確度的條件下，為了減少計(jì)算量，截?cái)嗖〝?shù)的取值應(yīng)盡可能小，所以應(yīng)當(dāng)更關(guān)注當(dāng)kr不太大時(shí)G(kr;f)的衰減特性。

對(duì)于格林函數(shù)中的聲源項(xiàng)，根據(jù)本文1.1節(jié)中式(4)～(6)，當(dāng)ξ＞1時(shí)，?s(ξ;f)隨ξ的變化近似滿足對(duì)數(shù)項(xiàng)與線性項(xiàng)疊加的關(guān)系，且當(dāng)ξ→1 時(shí)，?s(ξ;f) →∞；對(duì)于格林函數(shù)中的齊次項(xiàng)，當(dāng)kr的值不太大時(shí)，其影響也需要考慮，但由于齊次項(xiàng)包含數(shù)值計(jì)算結(jié)果A+1(kr)和A-1(kr)，難以對(duì)其變化特性進(jìn)行理論推導(dǎo)分析。

圖2中格林函數(shù)對(duì)數(shù)幅值?(ξ;f)的數(shù)值計(jì)算結(jié)果顯示，在1 ＜ξ＜2 的范圍內(nèi)，當(dāng)ξ接近于1 時(shí)?(ξ;f)衰減速率較快；隨著ξ的增大，?(ξ;f)的衰減速率逐漸減慢。?(ξ;f)的這種變化特征與聲源項(xiàng)?s(ξ;f)所滿足的對(duì)數(shù)和線性變化特征比較相似。另外，注意到式(2)中格林函數(shù)每一項(xiàng)都包含kz，1，因而?(ξ;f)中每一項(xiàng)都包含且ξ=1 是?(ξ;f)的一個(gè)極點(diǎn)。而在數(shù)值計(jì)算中，ξ=1處格林函數(shù)的計(jì)算結(jié)果不是無窮大，且極點(diǎn)的位置可能會(huì)存在偏差。因此，構(gòu)造數(shù)學(xué)模型時(shí)，不應(yīng)完全按照格林函數(shù)表達(dá)式把ξ=1確定為極點(diǎn)，而應(yīng)該把該極點(diǎn)的位置作為一個(gè)可調(diào)整的參數(shù)，否則數(shù)學(xué)模型可能與格林函數(shù)的數(shù)值計(jì)算結(jié)果失配。

基于上述對(duì)格林函數(shù)衰減特性理論和數(shù)值計(jì)算結(jié)果的分析，并參照式(5)和式(6)中?s(ξ;f)的近似結(jié)果，構(gòu)造了以下包含對(duì)數(shù)項(xiàng)和一次線性項(xiàng)的數(shù)學(xué)模型：

或

其中：?=[?1?2?3?4?5]T和θ=[θ1θ2θ3θ4]T為模型參數(shù)，ξt為截?cái)嗵幍乃讲〝?shù)比，1 ≤ξ≤ξt。根據(jù)對(duì)數(shù)函數(shù)的性質(zhì)，可證明式(10)和式(11)是等效的，而式(11)的參數(shù)更少，因此采用式(11)作為格林函數(shù)的擬合模型。在該模型中可以通過調(diào)整參數(shù)θ2去調(diào)整極點(diǎn)的位置，避免了因ξ=1附近極點(diǎn)位置的偏差而導(dǎo)致的模型失配。

利用此模型通過最小二乘法對(duì)格林函數(shù)對(duì)數(shù)幅值進(jìn)行擬合：

圖4 是利用式(11)中的數(shù)學(xué)模型對(duì)不同頻率?(ξ;f)的擬合結(jié)果，圖中粗虛線對(duì)應(yīng)圖2中各頻率的波數(shù)域格林函數(shù)對(duì)數(shù)幅值?(ξ;f)，實(shí)線是擬合結(jié)果。所有頻率的?(ξ;f)與對(duì)應(yīng)的擬合曲線幾乎重合，這表明利用式(11)對(duì)?(ξ;f)進(jìn)行擬合取得了較好的效果。

圖4 波數(shù)域格林函數(shù)對(duì)數(shù)幅值擬合結(jié)果Fig.4 Fitting results of the logarithmic magnitude of Green's function in wavenumber domain

本文構(gòu)造的式(11)中的數(shù)學(xué)模型符合格林函數(shù)的變化趨勢(shì)，具有一定的物理意義，可調(diào)整極點(diǎn)的位置，避免模型在極點(diǎn)處失配，并且復(fù)雜程度較低。因此，本文將采用式(11)中的數(shù)學(xué)模型對(duì)格林函數(shù)對(duì)數(shù)幅值?(ξ;f)進(jìn)行擬合。

1.3 模型參數(shù)自適應(yīng)濾波

卡爾曼濾波器是一種最優(yōu)線性狀態(tài)估計(jì)方法，被廣泛應(yīng)用于工程領(lǐng)域，本文利用卡爾曼濾波器對(duì)波數(shù)域格林函數(shù)的擬合參數(shù)進(jìn)行跟蹤和預(yù)測(cè)。1.2節(jié)構(gòu)造的數(shù)學(xué)模型h(ξ;θ)中包含4 個(gè)參數(shù)θj(n)，j=1，2，3，4，n是頻率序號(hào)。這4個(gè)參數(shù)之間是相互獨(dú)立的，用4個(gè)卡爾曼濾波器分別對(duì)其進(jìn)行跟蹤和預(yù)測(cè)。

式(11)中第j個(gè)參數(shù)的狀態(tài)方程為

第j個(gè)參數(shù)的觀測(cè)方程為

其中：zj(n)是狀態(tài)的觀測(cè)值；觀測(cè)矩陣H=[1 0]T；vj(n)為觀測(cè)噪聲，設(shè)觀測(cè)噪聲的協(xié)方差矩陣是Rj。

設(shè)第j個(gè)參數(shù)、第n步的卡爾曼增益矩陣為Kj(n)?？柭鼮V波的每個(gè)迭代過程，都包含了預(yù)測(cè)與修正兩個(gè)環(huán)節(jié)，在這兩個(gè)環(huán)節(jié)中都會(huì)對(duì)狀態(tài)矢量sj和狀態(tài)誤差協(xié)方差矩陣Mj進(jìn)行更新。為了對(duì)兩個(gè)環(huán)節(jié)的變量加以區(qū)別，用sj(n|n-1)表示根據(jù)第n-1 步狀態(tài)矢量所預(yù)測(cè)的第n步狀態(tài)矢量，Mj(n|n-1)表示根據(jù)第n-1步狀態(tài)誤差協(xié)方差矩陣所預(yù)測(cè)的第n步狀態(tài)預(yù)測(cè)誤差協(xié)方差矩陣；用sj(n|n)表示第n步修正后的狀態(tài)矢量，Mj(n|n)表示第n步修正后的狀態(tài)誤差協(xié)方差矩陣。

對(duì)每個(gè)模型參數(shù)θj(n)分別進(jìn)行卡爾曼濾波，其具體計(jì)算為[13]

(1) 初始條件

(2) 預(yù)測(cè)

(3) 狀態(tài)預(yù)測(cè)誤差協(xié)方差矩陣

(4) 卡爾曼增益

(5) 修正

(6) 狀態(tài)誤差協(xié)方差矩陣

其中：I為單位矩陣。

1.4 截?cái)嗖〝?shù)的計(jì)算

設(shè)期望的積分截?cái)嗑_度參數(shù)為a，表示期望在波數(shù)域格林函數(shù)的幅值衰減到10-a時(shí)進(jìn)行積分截?cái)?，即?/p>

其中：kt為截?cái)嗖〝?shù)。

把自適應(yīng)濾波器輸出的第n+1個(gè)頻率的參數(shù)估計(jì)結(jié)果代入數(shù)學(xué)模型式(11)中，可估計(jì)出第n+1個(gè)頻率的波數(shù)域格林函數(shù)：

把式(21)代入式(22)中可得第n+1 個(gè)頻率的截?cái)嗖〝?shù)ξt(fn+1)即為式(23)的解。

由式(11)可知方程(23)是超越方程，需要通過數(shù)值方法求解。

1.5 截?cái)嗖〝?shù)自適應(yīng)選取流程

圖5是本文提出的截?cái)嗖〝?shù)自適應(yīng)選取方法的流程圖。本文提出的截?cái)嗖〝?shù)自適應(yīng)選取方法包含以下步驟：

圖5 截?cái)嗖〝?shù)自適應(yīng)選取流程Fig.5 The flow chart of adaptive selcetion process of truncated wavenumber

(1) 初始截?cái)嗖〝?shù)ξt(f0)，ξt(f0)的取值可以較大，以保證聲場(chǎng)計(jì)算精度。

(2) 利用波數(shù)積分法計(jì)算聲場(chǎng)和?(ξ;fn)。

(3) 利用式(11)中的數(shù)學(xué)模型對(duì)進(jìn)行?(ξ;fn)擬合，得到參數(shù)θ(n)。把上一次迭代自適應(yīng)濾波器估計(jì)的參數(shù)作為求解此最優(yōu)化問題的初值，以加快擬合速度。

(4) 把參數(shù)θ(n)輸入卡爾曼濾波器，估計(jì)下一個(gè)頻率的模型參數(shù)

(6) 返回步驟(2)計(jì)算下一個(gè)頻率的聲場(chǎng)。

2 仿真驗(yàn)證

計(jì)算圖1所示海洋波導(dǎo)的波數(shù)域格林函數(shù)，并利用本文提出的方法選取積分截?cái)嗖〝?shù)，以驗(yàn)證本文提出方法的有效性。初始截?cái)嗖〝?shù)取值是ξt(f0) =5，計(jì)算頻率范圍是50～1 000 Hz?？柭鼮V波器的參數(shù)如表1所示。

表1 卡爾曼濾波器參數(shù)的取值Table 1 Values of Kalman filter parameters

對(duì)精確度參數(shù)a取不同值、頻率間隔Δf取不同值的情況進(jìn)行試驗(yàn)。圖6 是Δf分別為1、5 和10 Hz 時(shí)，a分別為3，4，5，6 的情況下，?(ξt(f);f)隨頻率變化，其中ξt(f)為本文方法選取的截?cái)嗖〝?shù)。

圖6 采用本文方法計(jì)算聲場(chǎng)時(shí)截?cái)嗵幉〝?shù)域格林函數(shù)的對(duì)數(shù)幅值隨頻率變化圖Fig.6 Logarithmic magnitude of Green's function in wavenumber domain at the truncation versus frequency when the method proposed in this paper is used for sound field calculation

圖6中的結(jié)果表明在上述試驗(yàn)條件下，本文提出的方法都能自適應(yīng)地選取積分截?cái)嗖〝?shù)，實(shí)現(xiàn)在期望的精確度附近進(jìn)行積分截?cái)?。在較低頻段，實(shí)際截?cái)嗵幏?[ξt(f);f]會(huì)出現(xiàn)波動(dòng)，隨著頻率的增大，波動(dòng)逐漸減小。這是因?yàn)榇_定截?cái)嗑_度a后，當(dāng)頻率較低時(shí)，?(ξ;f)隨頻率f的變化比較劇烈，因而模型參數(shù)的預(yù)測(cè)誤差比較大；當(dāng)頻率較高時(shí)，?(ξ;f)隨頻率f的變化比較緩慢，因而模型參數(shù)的預(yù)測(cè)誤差較小，截?cái)嗑_度誤差變小。

圖6說明了本文提出的自適應(yīng)截?cái)嗖〝?shù)選取方法能在寬帶聲場(chǎng)波數(shù)積分計(jì)算過程中，把截?cái)嗵幍母窳趾瘮?shù)幅值穩(wěn)定地保持在某一個(gè)精度附近。作為對(duì)比，圖7給出了當(dāng)ξt取固定值時(shí)，聲場(chǎng)計(jì)算過程中截?cái)嗵幉〝?shù)域格林函數(shù)的對(duì)數(shù)幅值隨頻率變化?？梢钥闯霎?dāng)ξt取固定值時(shí)，聲場(chǎng)計(jì)算過程中截?cái)嗑_度隨頻率的增大而不斷提高。

圖7 采用ξt取固定值法計(jì)算聲場(chǎng)時(shí)截?cái)嗵幉〝?shù)域格林函數(shù)的對(duì)數(shù)幅值隨頻率變化圖Fig.7 Logarithmic magnitude of Green's function in wavenumber domain at the truncation versus frequency when the fixed ξt method is used for sound field calculation

為了評(píng)估本文方法對(duì)聲場(chǎng)計(jì)算速度的影響，在其他條件都相同的情況下，分別利用常用的ξt取固定值的方法和本文提出的自適應(yīng)截?cái)嗖〝?shù)選取方法，計(jì)算50～1 000 Hz 寬帶聲場(chǎng)100 次，對(duì)比兩種方法的平均用時(shí)。為了量化對(duì)比兩種方法的積分截?cái)嗑_度，計(jì)算寬帶平均精確度amean：

其中：n為頻率序號(hào)，N為頻點(diǎn)個(gè)數(shù)，?[ξt(fn);fn]為截?cái)嗵幐窳趾瘮?shù)對(duì)數(shù)幅值。

表2 和表3 分別為采用ξt取固定值的方法和本文提出方法進(jìn)行聲場(chǎng)計(jì)算的平均用時(shí)和帶寬平均精確度。在平均用時(shí)比較接近的情況下(如表2中ξt=1.1，1.2，表3中a=3,4,5,6的情況)，采用本文方法選取截?cái)嗖〝?shù)能獲得更高的寬帶平均精確度；在寬帶平均精確度比較接近的情況下(如表2中ξt=2.0和表3 中a=3 的情況)，采用本文方法選取截?cái)嗖〝?shù)所需平均時(shí)間更少。這是因?yàn)棣蝨取固定值的方法無法兼顧低頻精確度和高頻計(jì)算量。由圖7可知，當(dāng)ξt取固定值1.1或1.2時(shí)，雖然計(jì)算量不大，但低頻聲場(chǎng)計(jì)算的積分截?cái)嗑_度太低，導(dǎo)致寬帶平均精確度較低；而當(dāng)ξt取固定值2.0或3.0時(shí)，雖然聲場(chǎng)計(jì)算整體的積分截?cái)嗑茸銐蚋?，但高頻聲場(chǎng)計(jì)算量較大，導(dǎo)致計(jì)算的平均用時(shí)較大。本文提出的方法在寬帶聲場(chǎng)計(jì)算過程中，積分截?cái)嗑_度可穩(wěn)定地維持在某一個(gè)值上，在保證了低頻積分截?cái)嗑_度的同時(shí)不會(huì)給高頻聲場(chǎng)計(jì)算增加不必要的計(jì)算量。

表2 水平波數(shù)比取固定值方法聲場(chǎng)計(jì)算平均用時(shí)和寬帶平均精確度Table 2 Average time and broadband average accuracy of sound field calculation with the fixed horizontal wavenumber ratio method

表3 本文提出的方法聲場(chǎng)計(jì)算平均用時(shí)和寬帶平均精確度Table 3 Average time and broadband average accuracy of sound field calculation with the method proposed in this paper

3 結(jié) 論

由于頻率、介質(zhì)類型、波導(dǎo)參數(shù)、聲源與接收點(diǎn)的相對(duì)深度等參數(shù)會(huì)對(duì)波數(shù)域格林函數(shù)的衰減特性產(chǎn)生很大的影響，在采用波數(shù)積分法計(jì)算不同海洋波導(dǎo)的聲場(chǎng)時(shí)需要選取不同的積分截?cái)嗖〝?shù)。針對(duì)波數(shù)積分寬帶聲場(chǎng)計(jì)算過程中，現(xiàn)有積分截?cái)嗖〝?shù)選取方法所存在的問題，本文提出一種應(yīng)用于流體介質(zhì)聲場(chǎng)計(jì)算的積分截?cái)嗖〝?shù)的自適應(yīng)選取方法。利用數(shù)學(xué)模型擬合波數(shù)域格林函數(shù)的衰減特性，根據(jù)已計(jì)算的波數(shù)域格林函數(shù)預(yù)測(cè)下一頻率的波數(shù)域格林函數(shù)的衰減特性，然后估計(jì)出積分截?cái)嗖〝?shù)。仿真試驗(yàn)結(jié)果表明，本文提出的方法實(shí)現(xiàn)了在給定的精確度附近進(jìn)行積分截?cái)啵诒ＷC積分精確度的同時(shí)不會(huì)引入太多額外的計(jì)算量。

本文提出的波數(shù)積分法截?cái)嗖〝?shù)自適應(yīng)選取方法除了可以用于寬帶多頻點(diǎn)的聲場(chǎng)計(jì)算外，還可以拓展到其他參數(shù)漸變的聲場(chǎng)計(jì)算場(chǎng)景。

本文關(guān)注的是各層都是流體介質(zhì)的海洋波導(dǎo)的情況，當(dāng)海底是彈性介質(zhì)或波導(dǎo)中包含冰蓋層時(shí)，波數(shù)域格林函數(shù)的衰減特性會(huì)有所變化。能否把本文提出的截?cái)嗖〝?shù)自適應(yīng)選取方法，應(yīng)用到包含彈性介質(zhì)的海洋波導(dǎo)的波數(shù)積分聲場(chǎng)計(jì)算，有待進(jìn)一步研究。

構(gòu)造函數(shù)模型擬合圖3中的曲線，可以得到給定精確度時(shí)截?cái)嗖〝?shù)隨頻率變化的經(jīng)驗(yàn)公式，這可能是解決此問題的另一種思路。但圖3中的曲線會(huì)隨波導(dǎo)參數(shù)的改變而改變，如何建立曲線參數(shù)與波導(dǎo)參數(shù)之間的關(guān)系有待深入研究。