李鵬程，馮海泓，李記龍

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

0 引 言

影響水聲目標(biāo)定位、目標(biāo)探測(cè)距離的一個(gè)重要因素就是水中的聲速[1]。海水中的聲速是對(duì)聲傳播有重要影響的海洋環(huán)境參數(shù)之一，在不同深度、不同地域會(huì)出現(xiàn)不同的聲速梯度，聲傳播也會(huì)隨著聲速梯度的不同而發(fā)生變化，聲線會(huì)向聲速較小的方向發(fā)生偏轉(zhuǎn)[2]。海洋中一定范圍內(nèi)的聲速數(shù)據(jù)可以通過(guò)溫鹽深儀(Conductivity,Temperature,Depth,CTD)、聲速剖面儀等設(shè)備直接測(cè)量獲得。但是通過(guò)直接測(cè)量獲取較大范圍和較深海域的聲速需要耗費(fèi)大量的時(shí)間，并且需要投入較大的人力物力，難以滿足實(shí)際中的需求[3]。

海水中的聲速剖面反演研究為快速獲取海區(qū)的實(shí)際聲速剖面提供了新的有效方法，并且可以實(shí)現(xiàn)長(zhǎng)時(shí)間的連續(xù)監(jiān)測(cè)[4]。在聲速剖面反演方面有許多學(xué)者進(jìn)行了研究，Tolstoy等[5]將匹配場(chǎng)處理應(yīng)用到海洋聲層析中，對(duì)經(jīng)驗(yàn)聲速剖面進(jìn)行經(jīng)驗(yàn)正交函數(shù)分解，提取聲速分布的特征，通過(guò)理論計(jì)算聲場(chǎng)與實(shí)測(cè)聲場(chǎng)進(jìn)行匹配，最優(yōu)聲場(chǎng)匹配結(jié)果對(duì)應(yīng)的聲速作為聲速反演結(jié)果。張忠兵[6]針對(duì)匹配場(chǎng)處理容易受到海底參數(shù)失配影響的問(wèn)題，提出匹配波束反演方法，利用聲信號(hào)到達(dá)接收點(diǎn)的時(shí)間差進(jìn)行反演。唐俊峰等[7]采用微擾法，將聲速剖面反演由非線性優(yōu)化轉(zhuǎn)化為線性方程組的形式，利用聲信號(hào)的傳播時(shí)間，提出了基于傳播時(shí)間的聲速剖面反演方法，以測(cè)得聲信號(hào)到達(dá)接收陣的時(shí)間差構(gòu)建代價(jià)函數(shù)，利用遺傳算法全局優(yōu)化性能得到最優(yōu)結(jié)果。張維等[8-9]針對(duì)不平整海底反射時(shí)的傳播規(guī)律，提出了一種三維空間特征聲線搜索和聲傳播計(jì)算模型，研究了海底深度、聲源位置和基陣傾斜等參數(shù)失配對(duì)聲速剖面反演性能的影響。胡軍[10]提出了基于徑向基函數(shù)(Radial Basis Function,RBF)神經(jīng)網(wǎng)絡(luò)的聲速剖面反演方法，利用RBF神經(jīng)網(wǎng)絡(luò)的逼近性能和自學(xué)習(xí)能力以及遺傳算法的全局優(yōu)化性能，得到最優(yōu)參數(shù)的反演模型，實(shí)現(xiàn)了基于海表實(shí)測(cè)溫度和歷史數(shù)據(jù)的聲速剖面反演。國(guó)家海洋局第二研究所研制了應(yīng)用于近海的基于聲線傳播時(shí)間層析的海洋聲層析儀，通過(guò)反演可以觀測(cè)海域的溫度、鹽度、聲速等海洋環(huán)境參數(shù)信息[11]。

海洋大地控制網(wǎng)是由一系列布設(shè)在海洋中的測(cè)量控制點(diǎn)構(gòu)成的大地控制網(wǎng)，是國(guó)家測(cè)繪局布設(shè)的“2000國(guó)家大地控制網(wǎng)”的重要組成部分，由陸上、海面和海底三方面的控制點(diǎn)組成，其中最主要的是海底控制點(diǎn)[12]。海底控制網(wǎng)是海洋大地測(cè)量、海洋測(cè)繪和水下導(dǎo)航定位的基礎(chǔ)，是海洋活動(dòng)尤其是水下航行的重要參考基準(zhǔn)，也是海洋權(quán)益維護(hù)、海洋資源勘探、海洋環(huán)境監(jiān)測(cè)的重要基礎(chǔ)設(shè)施[13]。海底控制點(diǎn)位置標(biāo)定后，海底控制點(diǎn)、水中或水面聲信號(hào)接收點(diǎn)及其相互間的聲信號(hào)傳輸過(guò)程也為聲速剖面反演帶來(lái)了新的思路。本文將海底控制點(diǎn)作為水下固定參考點(diǎn)，通過(guò)研究利用水下固定參考點(diǎn)與水面已知位置之間的聲信號(hào)傳播時(shí)延來(lái)反演海水聲速剖面的方法，嘗試求解逆問(wèn)題，建立一種針對(duì)目標(biāo)海域進(jìn)行快速、精確的聲速剖面反演方法。

1 基本原理

1.1 聲速反演模型建立

設(shè)水平方向上為x軸，垂直方向上為z軸，深度向下取正值，水面艦船的坐標(biāo)為(xF，zF)，根據(jù)海底參考點(diǎn)應(yīng)答器的定位結(jié)果可以確定水下參考點(diǎn)應(yīng)答器的坐標(biāo)為(xB，zB)。聲速在深度方向上的各層分布用C(z)表示，聲線與深度垂直向下方向的夾角為θ(x，z)，聲信號(hào)從水下參考點(diǎn)應(yīng)答器開(kāi)始發(fā)射，經(jīng)過(guò)時(shí)間τ被水面船只換能器基陣接收，如圖1所示。

根據(jù)斯涅耳(Snell)定律：

k(xF)為一條聲線所對(duì)應(yīng)的常數(shù)。對(duì)于一條聲線的傳播時(shí)間，有：

水面船只和水下參考點(diǎn)應(yīng)答器的水平距離為

根據(jù)式(1)中方程，可以得出：

代入式(2)、(3)，聲信號(hào)傳播時(shí)間為

水面船只和水下參考點(diǎn)應(yīng)答器的水平距離為

將水深等間距分為N層，水面船只接收聲信號(hào)位置記為M，根據(jù)斯涅耳定律，可得：

將式(7)、(8)離散化[14-15]，可得：

式(10)和(11)組成用來(lái)反演各層聲速的非線性方程組，方程組中聲信號(hào)傳播時(shí)延τm，每層深度hn，水下參考點(diǎn)與水面測(cè)量點(diǎn)的水平距離xB-xFm都為已知參數(shù)；未知數(shù)為各層聲速Cn和每條聲線對(duì)應(yīng)的常數(shù)值km。

1.2 牛頓迭代法求解非線性方程組

牛頓迭代法[16]是將非線性方程f(x)=0在x0點(diǎn)展開(kāi)，即：

令：

用線性方程p(x)=0 近似代替非線性方程f(x)=0，再?gòu)膒(x)=0 中解得x=x0-并令x1=x0-作為f(x)=0 的根的第一級(jí)近似值。同理，xk+1=xk-作為方程f(x)=0的根的第k+1級(jí)近似值，也作為牛頓迭代公式，xk為第k次迭代的初值。選取精度ε，判斷若|xk-x0|＜ε，則輸出滿足精度的解。

用牛頓迭代算法來(lái)求解由式(10)、(11)組成的非線性方程組，非線性方程組中方程的個(gè)數(shù)為2M，未知數(shù)為各層聲速Cn和每條聲線對(duì)應(yīng)的常數(shù)km，未知數(shù)個(gè)數(shù)為M+N+1，方程個(gè)數(shù)和未知數(shù)個(gè)數(shù)相等時(shí)方程組可解，所以M=N+1，即水面船只聲信號(hào)接收點(diǎn)的個(gè)數(shù)等于水深分層數(shù)加一。

聲速反演的流程為：

(1) 將水中聲速設(shè)為常數(shù)，即分為一層，水面選取兩個(gè)關(guān)于水下參考點(diǎn)對(duì)稱的測(cè)量點(diǎn)，通過(guò)迭代選取聲線的初始入射角，確定該測(cè)量點(diǎn)對(duì)應(yīng)聲線的常數(shù)k；

(2) 重復(fù)M次步驟(1)，確定2M個(gè)測(cè)量點(diǎn)處對(duì)應(yīng)聲線的k；

(3) 選取水下參考點(diǎn)其中一側(cè)的M個(gè)水面測(cè)量點(diǎn)，將水深分為M-1層，將各層平均聲速和每個(gè)測(cè)量點(diǎn)處的k作為初始迭代值代入，通過(guò)牛頓迭代法求解得到各層聲速。

2 仿真計(jì)算

仿真環(huán)境采用2 維坐標(biāo)，設(shè)水深為1 000 m，將水深等分為10 層，水下參考點(diǎn)坐標(biāo)為xB=0。選取常見(jiàn)的聲速剖面進(jìn)行仿真驗(yàn)證，分別為負(fù)梯度聲速剖面、正梯度聲速剖面和C形聲速剖面。

聲速剖面為負(fù)梯度聲速剖面時(shí)，設(shè)置的仿真參數(shù)如表1所示。設(shè)置好各層聲速后，隨機(jī)選取M個(gè)初始入射角，每條聲線對(duì)應(yīng)的常數(shù)k可由式(1)計(jì)算，因此仿真過(guò)程不需要選取水下參考點(diǎn)對(duì)稱的點(diǎn)，水面測(cè)量點(diǎn)的位置和時(shí)延值可以由式(10)和(11)計(jì)算獲得。建立關(guān)于聲信號(hào)傳播時(shí)間和接收位置水平距離的非線性方程組，將設(shè)置的聲速加上50%的正態(tài)分布誤差后作為迭代初始值，通過(guò)牛頓迭代法對(duì)非線性方程組迭代求解，解得各層聲速如表2所示。

表1 負(fù)梯度聲速剖面的仿真參數(shù)設(shè)置Table 1 Simulation parameter setting for negative gradient typed sound velocity profile

表2 各層反演聲速與假設(shè)聲速比較(對(duì)應(yīng)于表1)Table 2 Comparison between inverted and assumed sound velocities of each layer corresponding to Table 1

聲速剖面為正梯度聲速剖面時(shí)，設(shè)置的仿真參數(shù)如表3所示。建立關(guān)于聲信號(hào)傳播時(shí)間和接收位置水平距離的非線性方程組，將設(shè)置的聲速加上50%的正態(tài)分布誤差后作為迭代初始值，通過(guò)牛頓迭代法對(duì)非線性方程組迭代求解，解得各層聲速如表4所示。

表3 正梯度聲速剖面的仿真參數(shù)設(shè)置Table 3 Simulation parameter setting for positive gradient typed sound velocity profile

表4 各層反演聲速與假設(shè)聲速比較(對(duì)應(yīng)于表3)Table 4 Comparison between inverted and assumed sound velocities of each layer corresponding to Table 3

聲速剖面為C形聲速剖面時(shí)，設(shè)置的仿真參數(shù)如表5所示。建立關(guān)于聲信號(hào)傳播時(shí)間和接收位置水平距離的非線性方程組，將設(shè)置的聲速加上50%的正態(tài)分布誤差后作為迭代初始值，通過(guò)牛頓迭代法對(duì)非線性方程組迭代求解，解得各層聲速如表6所示。

表5 C形聲速剖面的仿真參數(shù)設(shè)置Table 5 Simulation parameter setting for C typed velocity speed profile

表6 各層反演聲速與假設(shè)聲速比較(對(duì)應(yīng)于表5)Table 6 Comparison between inverted and assumed sound velocities of each layer corresponding to Table 5

對(duì)三種聲速剖面進(jìn)行了反演仿真，仿真得到的聲速和仿真設(shè)置的聲速相比誤差較小。結(jié)果證明了反演模型的有效性。

3 海試數(shù)據(jù)驗(yàn)證

3.1 海試試驗(yàn)設(shè)計(jì)

試驗(yàn)數(shù)據(jù)來(lái)源于2021年5月某海上試驗(yàn)，在海上試驗(yàn)中，海況良好，試驗(yàn)海域的海底平坦，平均水深約為3 000 m，在海底布放1個(gè)可以發(fā)射聲信號(hào)的潛標(biāo)作為水下參考點(diǎn)，其中發(fā)射的聲信號(hào)為8～16 kHz 的寬帶復(fù)雜編碼信號(hào)，發(fā)射周期為6 s。實(shí)驗(yàn)船底部有接收聲信號(hào)的換能器基陣，船的位置可以通過(guò)GPS來(lái)精確定位，聲信號(hào)在實(shí)驗(yàn)船底換能器基陣和海底應(yīng)答器潛標(biāo)之間的信號(hào)時(shí)延也可以精確測(cè)量，其中GPS 的定位精度為0.2 m，測(cè)時(shí)精度為50 μs。潛標(biāo)標(biāo)定完畢后，實(shí)驗(yàn)船的航行軌跡如圖2所示。

圖2 牛頓迭代算法流程圖Fig.2 Flow chart of Newton iterative algorithm

圖3 實(shí)驗(yàn)船直線走航軌跡Fig.3 Straight route of the trial ship

3.2 聲速反演

將潛標(biāo)作為水下參考點(diǎn)應(yīng)答器，并將潛標(biāo)設(shè)置為坐標(biāo)系的參考坐標(biāo)系原點(diǎn)，在實(shí)驗(yàn)船航線上根據(jù)分層數(shù)來(lái)選擇適量的測(cè)量點(diǎn)。

將水深等分為10層時(shí)，每層深度為300 m，水面船聲信號(hào)接收點(diǎn)需要取11 個(gè)點(diǎn)。在潛標(biāo)兩側(cè)分別選取與水下參考點(diǎn)的水平距離分別為1 000，1 500，2 000，2 500，3 000，3 500，4 000，4 500，5 000，5 500，6 000 m 處作為水面的測(cè)量點(diǎn)。將CTD 實(shí)測(cè)聲速剖面數(shù)據(jù)按分的層數(shù)每層取平均，各層的平均值作為測(cè)量聲速來(lái)和反演計(jì)算聲速做對(duì)比，各層平均聲速的計(jì)算公式為

通過(guò)前文反演流程中的步驟(2)迭代求得這11個(gè)測(cè)量點(diǎn)處的聲線初始入射角和常數(shù)k。將測(cè)量聲速加上50 %的正態(tài)分布誤差后作為迭代的初始值，通過(guò)牛頓迭代算法迭代收斂得到各層水深的聲速，結(jié)果如表7所示，計(jì)算聲速與測(cè)量聲速的平均誤差為0.79 m·s-1。

表7 各層反演聲速與測(cè)量聲速比較Table 7 Comparison between inverted and measured sound velocities of each layer

將水深等分為30層時(shí)，每層深度為100 m，水面船聲信號(hào)接收點(diǎn)取31 個(gè)點(diǎn)。反演得到的聲速剖面與實(shí)測(cè)聲速剖面如圖4所示。

圖4 反演聲速剖面與真實(shí)聲測(cè)剖面比較Fig.4 Comparisons between inversion and real sound velocity profile

在不同分層情況下，反演聲速的平均誤差如表8所示。將水深分為3～30層的過(guò)程中，反演的聲速與實(shí)測(cè)聲速之間的平均誤差逐漸減小到0.8 m·s-1左右。分層數(shù)逐漸增加后，聲速梯度包含更多層的聲速信息，聲速隨深度的變化特征會(huì)更加明顯、更加具體；在分層數(shù)較少時(shí)，聲速梯度包含的聲速信息較少，聲速剖面每層之間的特征變化較大，不能完全體現(xiàn)聲速隨深度的變化特征，所以進(jìn)行聲速反演時(shí)誤差較大。當(dāng)分層增加到一定層數(shù)時(shí)，聲速變化已經(jīng)能夠體現(xiàn)聲速剖面的變化特征，而且在等聲速分層模型中，誤差只體現(xiàn)在各層的聲速上，在牛頓迭代過(guò)程中近似值滿足精度后輸出的結(jié)果就是滿足精度的聲速，因此反演誤差趨于穩(wěn)定。

表8 不同的分層數(shù)對(duì)應(yīng)的聲速反演的平均誤差Table 8 Average errors of sound velocity inversion corresponding to different numbers of layers

4 結(jié) 論

本文通過(guò)研究海底控制點(diǎn)定位的逆過(guò)程，提出了一種新的聲速剖面反演方法。通過(guò)仿真驗(yàn)證了反演模型的有效性。再利用某海試數(shù)據(jù)對(duì)模型進(jìn)行了進(jìn)一步的驗(yàn)證，將水深分為3～30 層的過(guò)程中，反演的聲速值與實(shí)際的聲速值之間的平均誤差會(huì)逐漸減小到0.8 m·s-1左右。但是隨著分層數(shù)增加，非線性方程組變得越來(lái)越復(fù)雜，反演效率會(huì)逐漸下降。

由于解非線性方程組時(shí)要先確定解的范圍，再確定精度，所以當(dāng)初始迭代值改變時(shí)，迭代后的解會(huì)發(fā)生改變，迭代后的解是范圍內(nèi)滿足精度的解，每次的解雖然不同，但是都在提前設(shè)的初始范圍內(nèi)，滿足精度，誤差小。針對(duì)初值選取的問(wèn)題，還有待進(jìn)一步的研究。