張 維，楊士莪，黃益旺，唐俊峰，宋 揚

(哈爾濱工程大學 水聲技術重點實驗室，哈爾濱 150001)

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

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

1 實驗描述

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

圖1 浮標與投彈點位置Fig.1 The position of buoys and bomb-dropping

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

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

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

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

圖2 實測聲速剖面Fig.2 The measured SSP

2 聲傳播時間計算模型

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

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

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

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

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

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

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

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

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

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

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

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

3 聲傳播時間精度分析

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

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

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

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

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

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

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

4 實驗數(shù)據(jù)處理結果

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

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

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

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

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

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

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

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

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

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

圖6 現(xiàn)場實測聲速Fig.6 The measured filed SSP

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

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

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

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

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

5 結論

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

[1]沈遠海，馬遠良，屠慶平，等.淺海聲速剖面的反演方法與實驗驗證[J].西北工業(yè)大學學報，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]唐俊峰，楊士莪.由傳播時間反演海水中的聲速剖面[J].哈爾濱工程大學學報，2006，27(5):733－737.

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

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

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

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

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

[12]王恕銓.求解三維本征聲線的一種新方法[J].聲學學報，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]劉伯勝，雷家煜.水聲學原理[M].哈爾濱:哈爾濱工程大學出版社，2006:112.

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

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