基于爆炸聲傳播時(shí)間的聲速剖面反演

張 維，楊士莪，黃益旺，唐俊峰，宋 揚(yáng)

(哈爾濱工程大學(xué) 水聲技術(shù)重點(diǎn)實(shí)驗(yàn)室，哈爾濱 150001)

聲速剖面反演是海洋聲層析技術(shù)(Acoustic Tomography)的一種，其實(shí)質(zhì)是根據(jù)觀測(cè)信號(hào)的某些特征信息(聲線傳播時(shí)間、簡正波相位等)來反推聲速剖面的問題[1]。聲速剖面反演首先由Munk等[2]提出，他們?cè)谏詈－h(huán)境下利用聲線的相對(duì)時(shí)延進(jìn)行反演，后來，簡正波相位反演[3]、峰值匹配[4]以及匹配場(chǎng)方法[5]等多種反演方法逐漸發(fā)展起來。近年來，我國學(xué)者在利用二維特征聲線的傳播時(shí)間反演聲速剖面上做了大量工作[6－7]，即假定特征聲線在同一垂直剖面內(nèi)傳播，這在平坦海底或者海底坡度非常小的情況下是可行的。然而，在大陸架海域，由于海底坡度較大，聲線在海底反射時(shí)將在水平方向上產(chǎn)生明顯的偏轉(zhuǎn)，使得聲線在三維空間傳播。此時(shí)，若仍用二維射線聲學(xué)的方法反演聲速剖面將給特征聲線傳播時(shí)間的計(jì)算帶來不必要的誤差，從而降低反演精度。針對(duì)此問題，本文利用南海海洋環(huán)境反演實(shí)驗(yàn)數(shù)據(jù)，在三維空間搜索特征聲線，分析了聲線的水平偏轉(zhuǎn)對(duì)傳播時(shí)間的影響，并基于最快特征聲線的傳播時(shí)間進(jìn)行了聲速剖面反演。采用該方法進(jìn)行聲速剖面反演具有以下優(yōu)勢(shì):① 相比于其他特征量，例如聲壓幅度、相位、多途時(shí)延等，最快特征聲線的傳播時(shí)間易于獲得;② 聲傳播時(shí)間不受海底聲速、密度等參數(shù)的影響，在反演過程中無需知道這些參數(shù)的準(zhǔn)確值，提高了噪聲容限[8];③ 射線方法能進(jìn)行三維聲場(chǎng)的計(jì)算，進(jìn)而反演具有復(fù)雜海底的海洋環(huán)境下的聲速剖面;④ 聲傳播時(shí)間計(jì)算速度快、精度高。

聲速剖面反演是一個(gè)多維優(yōu)化的問題，研究表明，用少量的幾階經(jīng)驗(yàn)正交函數(shù)表示聲速剖面是一種有效的方法[9]。本文通過采用全局優(yōu)化效果較好且計(jì)算速度較快的量子粒子群算法[10]搜索經(jīng)驗(yàn)正交函數(shù)系數(shù)的全局最優(yōu)解來反演聲速剖面，對(duì)考慮和不考慮聲線水平偏轉(zhuǎn)兩種情況下的反演結(jié)果進(jìn)行了對(duì)比。結(jié)果表明，考慮聲線的水平偏轉(zhuǎn)能使傳播時(shí)間的計(jì)算誤差降低一個(gè)數(shù)量級(jí)，進(jìn)而有效地提高聲速剖面反演的精度。

1 實(shí)驗(yàn)描述

實(shí)驗(yàn)海域如圖1所示。海底深度變化較大，1號(hào)浮標(biāo)附近海深大約為100 m，3號(hào)浮標(biāo)附近海深大于1 000 m。文獻(xiàn)[11]指出，在實(shí)驗(yàn)海域周圍布放多個(gè)聲發(fā)射點(diǎn)和接收點(diǎn)，聲信號(hào)將從不同角度不同途徑攜帶聲速信息經(jīng)到達(dá)接收點(diǎn)，途徑越多，所獲取的信息量越大，聲層析的精度越高。另外，水聽器間隔太小從聲傳播時(shí)間上不易分辨，太大則會(huì)給吊放和回收帶來困難，考慮到以上因素在實(shí)驗(yàn)中四個(gè)浮標(biāo)大致成正方形分布，每個(gè)浮標(biāo)連接一個(gè)長約50 m帶有5個(gè)水聽器的垂直陣，水聽器的深度分別為 10 m，18 m，26 m，34 m，42 m。浮標(biāo)上端都帶有GPS裝置以實(shí)時(shí)記錄該浮標(biāo)的位置，垂直陣下端系有重物使得陣始終保持垂直狀態(tài)。發(fā)射船在浮標(biāo)間按“W”形航跡運(yùn)動(dòng)，每隔一段時(shí)間投放聲彈或手榴彈。圖1中A，B，C，D是投彈點(diǎn)的四個(gè)位置，A，B，C三點(diǎn)投放聲彈，爆炸深度是100 m，D點(diǎn)投放手榴彈，爆炸深度是8 m。聲源爆炸后，各水聽器在不同距離不同深度上接收到聲信號(hào)。定義投彈點(diǎn)和浮標(biāo)之間的連線與經(jīng)度增大方向的夾角為方位角，逆時(shí)針為正。圖1中r和α分別表示聲源與浮標(biāo)之間的水平距離和方位角。

圖1 浮標(biāo)與投彈點(diǎn)位置Fig.1 The position of buoys and bomb-dropping

實(shí)驗(yàn)時(shí)，聲爆炸時(shí)刻和爆炸位置都是未知的，只有通過其他的信息估算得到。設(shè)GPS每n秒記錄一次發(fā)射船的位置;投彈前后兩次GPS記錄的時(shí)刻分別為T1和T2，對(duì)應(yīng)的發(fā)射船的水平位置分別為(x1，y1)和(x2，y2);聲源下沉平均速度為v，爆炸深度為h，發(fā)射船監(jiān)聽到爆炸聲的時(shí)刻為T0;從海面到爆炸深度的平均聲速為c0。估算聲爆炸時(shí)刻T和爆炸的水平位置(x，y)的方法如下:聲源下沉?xí)r間，該過程中發(fā)射船前行距離，爆炸聲從爆炸位置傳播到監(jiān)聽水聽器的時(shí)間則聲爆炸時(shí)刻，爆炸點(diǎn)水平位置

經(jīng)計(jì)算后，聲源與浮標(biāo)之間的水平距離和方位角如表1和表2所示。在實(shí)驗(yàn)過程中連續(xù)多天采用CTD等設(shè)備對(duì)該海域不同位置的聲速進(jìn)行了測(cè)量，測(cè)量結(jié)果如圖2所示。

表1 各聲源與浮標(biāo)之間的水平距離(km)Tab.1 The horizontal distance between sources and buoys

表2 各聲源與浮標(biāo)之間的方位角(°)Tab.2 The azimuthal angle between sources and buoys

圖2 實(shí)測(cè)聲速剖面Fig.2 The measured SSP

2 聲傳播時(shí)間計(jì)算模型

作為聲速剖面反演的代價(jià)函數(shù)，聲傳播時(shí)間是影響反演結(jié)果精確程度的最大因素，因此有必要建立合理的計(jì)算模型進(jìn)行特征聲線的搜索及傳播時(shí)間的計(jì)算。

特征聲線搜索的關(guān)鍵點(diǎn)在于出射掠射角的確定，在三維情況下還需要考慮出射方位角的問題[12－13]。為此，本文采用在垂直方向上迭代插值求掠射角，在水平方向上迭代補(bǔ)償求方位角的辦法進(jìn)行三維特征聲線的搜索。

設(shè)聲源和接收點(diǎn)的位置分別為(x0，y0，z0)和(xr，yr，zr)。由圖2可以看到，聲速在水平方向上變化很小，因此可以假設(shè)聲速只是深度z的函數(shù)。

聲線在傳播過程中滿足Snell定理和射線方程，分別如式(1)和(2)所示[14]:

式中:θ和 θ0分別是深度 z和 z0處的掠射角，(xj，yj，zj)和(xj+1，yj+1，zj+1)是第 j段聲線的起始點(diǎn)和終止點(diǎn)，uj，vj，sj分別是折射率 n在 x，y和 z方向上的分量。若方位角為 α，則 uj，vj，sj滿足:

海底反射時(shí)，設(shè)海底深度h是x，y的函數(shù)，h在x軸和y軸上的一次導(dǎo)數(shù)分別是b1和b2，令F=(1+b21)，則海底法向量可以表示為設(shè)入射向量為，反射向量為→R =，令 W=s－ bμ － bν，根據(jù)入射向j1j2j量、反射向量以及海底法向量之間的關(guān)系可以得到:

式(4)表征了各反射分量與入射分量之間的關(guān)系。當(dāng)海底梯度不為0，反射向量與入射向量在水平方向上的分量不相等，聲線發(fā)生水平偏轉(zhuǎn)。

在三維空間搜索特征聲線時(shí)，首先假設(shè)初始出射方位角α為聲源與接收點(diǎn)的連線方向，出射掠射角在設(shè)定范圍內(nèi)按一定步長分為N等分，按照式(1)～(4)所示的Snell定理，射線方程以及海底反射方程進(jìn)行單根聲線跟蹤直到聲線終點(diǎn)，然后利用:

若在二維空間搜索特征聲線，可以認(rèn)為海底反射點(diǎn)處b1≈0，b2≈0，即假設(shè)海底深度按照臺(tái)階式變化，則 uj+1≈uj，vj+1≈vj，sj+1≈ － sj，從而使得聲線只在一個(gè)垂直剖面內(nèi)傳播。此時(shí)，出射方位角為聲源與接收點(diǎn)的連線方向。只需對(duì)特征聲線的相鄰聲線進(jìn)行掠射角線性迭代插值獲得特征聲線的出射掠射角，即可確定特征聲線。

某段聲線從深度z1傳播到深度 z所需要的時(shí)間為:

在反轉(zhuǎn)點(diǎn)附近由于sinθ≈0，為了減小數(shù)值計(jì)算的誤差，假定該點(diǎn)附近的聲速梯度是恒定的，聲速軌跡為一段圓弧，式(6)可以簡化為[15]:

式中:θ0和θ'分別為弧線起始處與結(jié)束處聲線的掠射角，g是恒定的聲速梯度。

3 聲傳播時(shí)間精度分析

聲傳播時(shí)間計(jì)算的精確性是有效反演聲速剖面的關(guān)鍵，而聲線的水平偏轉(zhuǎn)導(dǎo)致二維特征聲線和三維特征聲線的傳播時(shí)間不一致，從而使得反演結(jié)果存在差異，因此有必要對(duì)考慮和不考慮聲線水平偏轉(zhuǎn)時(shí)聲傳播時(shí)間的精度進(jìn)行研究。

采用實(shí)驗(yàn)現(xiàn)場(chǎng)所測(cè)聲速剖面，分別在二維空間和三維空間搜索最快特征聲線，并計(jì)算其傳播時(shí)間。將聲源A，1號(hào)浮標(biāo)作為算例1，聲源B，2號(hào)浮標(biāo)作為算例2，聲源C，3號(hào)浮標(biāo)作為算例3，聲源D，4號(hào)浮標(biāo)作為算例4。算例中均采用34 m深度的水聽器。其中，算例2海底較平緩，算例4海底較陡峭，兩個(gè)算例的最快特征聲線傳播軌跡如圖3所示。

從圖3中可以看出，由于聲線在海底的反射，二維特征聲線與三維特征聲線差別較大，尤其是算例4，海底坡度較大，聲線的水平偏轉(zhuǎn)較嚴(yán)重。

最快特征聲線的傳播時(shí)間如表3所示。

圖3 最快特征聲線傳播軌跡Fig.3 Ray trajectory of rapidest eigenray

表3 最快特征聲線傳播時(shí)間Tab.3 Transmission time of rapidest eigenray

其中:t是從各水聽器接收信號(hào)提取的聲信號(hào)到達(dá)時(shí)間，t1是在三維空間搜索時(shí)最快特征聲線的傳播時(shí)間，t2是在二維空間搜索時(shí)最快特征聲線的傳播時(shí)間。從表3可以看出，二維特征聲線傳播時(shí)間的誤差明顯比三維特征聲線傳播時(shí)間的誤差大一個(gè)數(shù)量級(jí)。三維特征聲線的誤差主要由計(jì)算的累積誤差產(chǎn)生，一般地，距離越遠(yuǎn)，這種誤差越大。二維特征聲線的誤差主要由忽略了聲線的水平偏轉(zhuǎn)所帶來，一般地，海底地形越復(fù)雜，坡度越大，這種誤差越大。另外，二維空間搜索計(jì)算得到的聲傳播時(shí)間在四種算例中均比實(shí)測(cè)的聲傳播時(shí)間要小。

4 實(shí)驗(yàn)數(shù)據(jù)處理結(jié)果

根據(jù)17次測(cè)量的聲速剖面 c1(zi)，c2(zi)，…，c17(zi)和平均聲速剖面)求協(xié)方差矩陣R，R的每一個(gè)元素 rij可以表示為[16]:

將協(xié)方差矩陣進(jìn)行特征分解，選取前三個(gè)最大的特征值所對(duì)應(yīng)的特征向量作為經(jīng)驗(yàn)正交函數(shù)，平均聲速剖面和經(jīng)驗(yàn)正交函數(shù)分別如圖4和圖5所示。

圖4 平均聲速剖面Fig.4 The average SSP

圖5 經(jīng)驗(yàn)正交函數(shù)Fig.5 Experiential orthogonal functions

用平均聲速剖面和經(jīng)驗(yàn)正交函數(shù)表示待反演聲速剖面如式(9)所示[17]:

其中:ak(x，y)是待反演參數(shù)，由于整個(gè)實(shí)驗(yàn)過程中聲速在水平方向上變化不是很大，因此ak(x，y)可近似認(rèn)為是常數(shù)。

選取代價(jià)函數(shù)為:

其中:ri，rj分別為聲源和接收點(diǎn)的坐標(biāo)，Tj是爆炸時(shí)刻，tij是各水聽器收到爆炸聲時(shí)刻的計(jì)算值，τij是實(shí)驗(yàn)中各水聽器收到爆炸聲時(shí)刻的觀測(cè)值。將A、B、C、D四個(gè)聲源分為A和D，B和D，C和D三組，從各水聽器接收信號(hào)提取聲傳播時(shí)間，進(jìn)行三次聲速剖面反演。采用量子粒子群作為優(yōu)化算法，搜索當(dāng)代價(jià)函數(shù)最大時(shí)ak的全局最優(yōu)值并代入式(9)獲得反演的聲速剖面，并定義

為均方根誤差，其中，N為深度點(diǎn)數(shù)，c(z)為實(shí)測(cè)聲速剖面，c'(z)為反演聲速剖面。

實(shí)驗(yàn)現(xiàn)場(chǎng)所測(cè)量的聲速如圖6所示，三組反演得到的聲速剖面及誤差如圖7所示。由于聲速變化較大的是在0～200 m深度，200 m以下變化非常小，為更清楚地進(jìn)行比較，圖7中深度只截取0～200 m。圖7中用以和反演結(jié)果對(duì)比的實(shí)測(cè)聲速即為圖6中所示的現(xiàn)場(chǎng)實(shí)測(cè)聲速。

圖6 現(xiàn)場(chǎng)實(shí)測(cè)聲速Fig.6 The measured filed SSP

從圖7可以看出，三維特征聲線聲速剖面反演的結(jié)果與二維情況下相比誤差更小，更接近實(shí)測(cè)聲速剖面。

圖7 聲速剖面反演結(jié)果Fig.7 The results of SSP inversion

各組反演的均方根誤差和最大誤差如表4所示。

表4 聲速剖面反演誤差Tab.4 The error of SSP inversion

從表4可以看出，在以上三組反演結(jié)果中，第三組反演的誤差最大，忽略聲線水平偏轉(zhuǎn)時(shí)最大誤差是6.604 m/s，均方根誤差為 3.742 m/s，考慮聲線水平偏轉(zhuǎn)時(shí)最大誤差為3.587 m/s，均方根誤差為1.887m/s，也就是說考慮了聲線的水平偏轉(zhuǎn)，提高了聲傳播時(shí)間的計(jì)算精度，使得聲速剖面反演的最大誤差和均方根誤差都減小了將近一倍。

5 結(jié)論

本文首先提出了一種在三維空間搜索特征聲線并計(jì)算其傳播時(shí)間的方法。針對(duì)實(shí)驗(yàn)海域的實(shí)際情況計(jì)算了最快特征聲線的傳播時(shí)間，并與忽略聲線水平偏轉(zhuǎn)情況下的結(jié)果進(jìn)行了對(duì)比。通過對(duì)比發(fā)現(xiàn)考慮聲線水平偏轉(zhuǎn)能使得聲傳播時(shí)間計(jì)算的精度提高一個(gè)數(shù)量級(jí)。另外，通過對(duì)實(shí)驗(yàn)數(shù)據(jù)的處理，以最快特征聲線傳播時(shí)間為代價(jià)函數(shù)并用量子粒子群算法作為優(yōu)化算法反演了實(shí)驗(yàn)海域的聲速剖面。結(jié)果表明，考慮聲線的水平偏轉(zhuǎn)能使得反演的最大誤差和均方根誤差都減小將近一倍。由于本文所采用的聲速數(shù)據(jù)在水平方向上的變化很小(見圖2)，因此僅反演了聲速在深度方向上的變化情況，并不是三維聲速反演。在遠(yuǎn)距離聲傳播問題中，聲速在水平方向上的變化將不容忽略，如何快速有效地反演出聲速在三維空間的變化情況將是下一步的工作。

[1]沈遠(yuǎn)海，馬遠(yuǎn)良，屠慶平，等.淺海聲速剖面的反演方法與實(shí)驗(yàn)驗(yàn)證[J].西北工業(yè)大學(xué)學(xué)報(bào)，2000，18(2):212－215.

[2]Munk W H，Wunsch C.Ocean acoustic tomography A:scheme for large scale monitoring[J].Deep-Sea Res，1979，26:123－161.

[3]Shang E C.Ocean acoustic tomography based on adiabatic mode theory[J].J.Acoust Soc Am，1989，85(4):1531－1537.

[4]Skarsoulis E K，Athanassoulis G A.Ocean acoustic tomography based on peak arrivals[J].J.Acoust Soc Am，1996，100(2):797－813.

[5]Taroudakis M I，Markaki M G.On the use of matched-field processing and hybrid algorithms for vertical slice tomography[J].J.Acoust Soc Am，1997，102(2):885 －895.

[6]唐俊峰，楊士莪.由傳播時(shí)間反演海水中的聲速剖面[J].哈爾濱工程大學(xué)學(xué)報(bào)，2006，27(5):733－737.

[7]張忠兵，馬遠(yuǎn)良，倪晉平，等.基于聲線到達(dá)時(shí)差的聲速剖面反演[J].西北工業(yè)大學(xué)學(xué)報(bào)，2002，20(1):36 －39.

[8]黃益旺.淺海遠(yuǎn)距離匹配場(chǎng)聲源定位研究[D].哈爾濱:哈爾濱工程大，2005:1－2.

[9]沈遠(yuǎn)海，馬遠(yuǎn)良，屠慶平，等.淺水聲速剖面用經(jīng)驗(yàn)正交函數(shù)(EOF)表示的可行性研究[J].應(yīng)用聲學(xué)，1999，18(2):21－25.

[10]楊傳將，劉 清，黃 珍.一種量子粒子群算法的改進(jìn)方法[J].計(jì)算技術(shù)與自動(dòng)化，2009，28(1):100 －103.

[11]廖光洪，朱小華，林 巨，等.海洋聲層析觀測(cè)技術(shù)與方法[J].海洋學(xué)報(bào)，2010，32(3):14－21.

[12]王恕銓.求解三維本征聲線的一種新方法[J].聲學(xué)學(xué)報(bào)，1992，17(2):155－157.

[13]Sethian J A.3－D travel time computation using the fast matching method.[J].Geophysics，1999，64(2):516－523.

[14]Yang S E.Theory of Underwater sound propagation[M].Harbin Engineering University Press，2009.

[15]劉伯勝，雷家煜.水聲學(xué)原理[M].哈爾濱:哈爾濱工程大學(xué)出版社，2006:112.

[16]張之夢(mèng)，劉伯勝.遺傳模擬退火算法用于淺海聲速反演的仿真研究[J].哈爾濱工程大學(xué)學(xué)報(bào)，2006，27(4):505－508.

[17]張忠兵，馬遠(yuǎn)良，楊坤德，等.淺海聲速剖面的匹配波束反演方法[J].聲學(xué)學(xué)報(bào)，2005，30(2):103－107.