循環(huán)神經(jīng)網(wǎng)絡(luò)在淺海聲速-聲源聯(lián)合反演中的應(yīng)用

王同， 蘇林， 任群言， 王文博， 馬力

(1.中國科學(xué)院 水聲環(huán)境特性重點實驗室， 北京 100190；2.中國科學(xué)院 聲學(xué)研究所， 北京 100190；3.中國科學(xué)院大學(xué)， 北京 100049)

海水中的聲速決定了海洋中的聲傳播特性[1]，剖面受到很多海水環(huán)境參數(shù)的影響，其中有一些相對穩(wěn)定的參數(shù)，如海底的結(jié)構(gòu)、海洋的深度等；還有一些動態(tài)變化的參數(shù)，如海水的生物群、海流、溫度、鹽度等。這些動態(tài)變化的參數(shù)，會隨著時間的變化而發(fā)生變化。變化不穩(wěn)定，會導(dǎo)致海水中聲速產(chǎn)生變化，使得水層聲速剖面具有明顯的時空變化。

海水聲速剖面的顯著時間演化特性為聲速剖面的序貫反演提供了可行性。Candy[2]將擴(kuò)展卡爾曼濾波算法引入到海洋聲學(xué)反演中，將聲速剖面的時間演化過程建立為一階時間演化模型。Yardim[3]通過經(jīng)驗正交函數(shù)(empirical orthogonal function，EOF)描述聲速剖面，采用序貫濾波中的卡爾曼濾波和粒子濾波的方法對地聲參數(shù)進(jìn)行跟蹤。Carriere[4-5]在水平變化環(huán)境中將聲學(xué)反演問題改進(jìn)為狀態(tài)-空間模型，通過序貫濾波的數(shù)據(jù)同化算法進(jìn)行海洋環(huán)境狀態(tài)估計。李建龍等[6]通過海洋動力學(xué)模型與聲速的關(guān)系反演聲速場的時間與空間結(jié)構(gòu)信息并驗證了算法的有效性；金麗玲等[7]利用自回歸分析擬合方法將聲速場擾動建模為高階自回歸演化模型，通過集合卡爾曼濾波序貫的估計時間演化的海洋聲速場。這些研究將時變聲速剖面的反演問題改進(jìn)為由描述聲速剖面時變特性的狀態(tài)方程與包含測量聲場信息的測量方程組成的狀態(tài)-空間模型，應(yīng)對聲場環(huán)境的時變特性。

由于淺海聲場環(huán)境的復(fù)雜性，難以直接描述聲速剖面時變特性的狀態(tài)方程，通常為一階隨機(jī)游走過程，此時的狀態(tài)方程不能很好的預(yù)測狀態(tài)變量的走向，序貫濾波的應(yīng)用會受到限制。而深度學(xué)習(xí)的分析歸納能力強，可以將復(fù)雜的問題表示為嵌套的層次概念體系，所以嘗試?yán)蒙疃葘W(xué)習(xí)方法的循環(huán)神經(jīng)網(wǎng)絡(luò)，處理相關(guān)數(shù)據(jù)，描述聲速剖面的時變特性。學(xué)習(xí)數(shù)據(jù)構(gòu)建模型并配合序貫濾波器進(jìn)行估計，這一方法在其他學(xué)科中廣泛應(yīng)用。例如，文獻(xiàn)[8]中將長短期記憶人工神經(jīng)網(wǎng)絡(luò)(long short term memory network，LSTM)與卡爾曼濾波(Kalman filter，KF)深度融合，提出了一種KF修正的LSTM混合模型KF-LSTM，其結(jié)果與LSTM模型、淺層神經(jīng)網(wǎng)絡(luò)模型和KF模型單獨使用的結(jié)果相比，能夠有效提升性能；文獻(xiàn)[9]中將BP神經(jīng)網(wǎng)絡(luò)用于測量方程的優(yōu)化，再通過卡爾曼濾波器進(jìn)行估計；文獻(xiàn)[10]利用神經(jīng)網(wǎng)絡(luò)對非線性系統(tǒng)建立狀態(tài)空間模型，采用高階容積卡爾曼濾波對新的狀態(tài)進(jìn)行實時更新，從而達(dá)到神經(jīng)網(wǎng)絡(luò)對非線性系統(tǒng)模型的真實逼近以及對狀態(tài)值的精確估計。

本文將深度學(xué)習(xí)中的循環(huán)神經(jīng)網(wǎng)絡(luò)引入到序貫濾波中，利用循環(huán)神經(jīng)網(wǎng)絡(luò)學(xué)習(xí)歷史水文數(shù)據(jù)，對淺海環(huán)境下的時變聲速剖面進(jìn)行學(xué)習(xí)建模，在此基礎(chǔ)上，利用集合卡爾曼濾波進(jìn)行了對聲速剖面和聲源位置的聯(lián)合反演，通過實測聲速剖面的仿真聲學(xué)數(shù)據(jù)驗證了方法的可行性。

1 聲速剖面時變特性模型

LSTM 神經(jīng)網(wǎng)絡(luò)由Hochreiter等[11]提出，并由Graves[12]進(jìn)行改進(jìn)，是基于RNN的一種完善，解決RNN中易出現(xiàn)的梯度消亡問題。在LSTM中，通過門結(jié)構(gòu)來對細(xì)胞狀態(tài)信息進(jìn)行增加或刪除。門結(jié)構(gòu)通常由一個sigmoid神經(jīng)網(wǎng)絡(luò)層和逐點乘積的操作組成，從而選擇性地讓信息通過。LSTM通過輸入門、遺忘門和輸出門3個門結(jié)構(gòu)來實現(xiàn)對信息的控制。

基于循環(huán)神經(jīng)網(wǎng)絡(luò)進(jìn)行模型構(gòu)建與訓(xùn)練。模型訓(xùn)練的數(shù)據(jù)集來自于某次淺海實驗中A、B2點約10 d的水文數(shù)據(jù)，數(shù)據(jù)集中包括A點聲速剖面和B點聲速剖面，每10 s記錄一次數(shù)據(jù)。利用這個數(shù)據(jù)集構(gòu)建淺海環(huán)境下的聲速剖面時變特性模型，最終達(dá)到已知A點處前1 h的水文數(shù)據(jù)，可以預(yù)測B點1 h后的聲速剖面。在Keras框架下進(jìn)行了模型的構(gòu)建與訓(xùn)練。

在模型訓(xùn)練結(jié)束后，使用測試集數(shù)據(jù)進(jìn)行模型測試。圖1為實際聲速剖面、模型預(yù)測的聲速剖面、模型預(yù)測的誤差。可以看出，模型預(yù)測值的變化趨勢，與實際聲速剖面的時間演化趨勢基本一致，只在溫躍層深度處，因聲速剖面起伏較大，模型的預(yù)測誤差變大。圖2為測試集數(shù)據(jù)進(jìn)行預(yù)測的均方根誤差，在溫躍層起伏較大時，均方根誤差較大。

圖1 聲速剖面時變特性模型的預(yù)測結(jié)果Fig.1 Predicted results based on lstm model using SSP time-varying characters

圖2 模型預(yù)測結(jié)果的均方根誤差Fig.2 Root mean squared error of predicted results based on lstm model

2 基于聲速剖面時變特性模型的聯(lián)合反演

2.1 聲速剖面參數(shù)化

EOF方法能夠有效地將聲速剖面參數(shù)化，降低表示聲速剖面所需要參數(shù)個數(shù)[13-14]。假設(shè)在一段時間內(nèi)測量得到N條聲速剖面SSP數(shù)據(jù)，在深度上進(jìn)行等間隔線性插值。離散到M個深度上：c1(zj),c2(zj),…,cN(zj),j=1,2,…,M

從而可以通過求平均值的辦法得到平均聲速剖面：

(1)

則第i條SSP在第j個深度處的聲速起伏為：

(2)

協(xié)方差矩陣R可以表示為：

(3)

(4)

對矩陣R進(jìn)行特征值分解，可以得到：

(5)

式中λ1>λ2>…>λM為對應(yīng)特征向量fn的特征值。

由于本征值快速衰減，通常前幾階數(shù)EOF系數(shù)就可以較為準(zhǔn)確地描述聲速剖面：

(6)

式中αk和fk(z)分別為第k階EOF系數(shù)和第k階EOF向量。

2.2 空間-狀態(tài)方程

本文利用集合卡爾曼濾波方法[15]對淺海環(huán)境下時變聲速剖面及聲源位置同時進(jìn)行反演，系統(tǒng)的狀態(tài)方程需要同時描述聲速剖面和聲源位置的變化。在先驗數(shù)據(jù)的基礎(chǔ)上，利用前3階經(jīng)驗正交系數(shù)表示聲速剖面:

(7)

式中:wk為狀態(tài)噪聲向量；L表示之前所構(gòu)建的聲速剖面時變特性模型，以A點前1 h的聲速剖面數(shù)據(jù)作為先驗數(shù)據(jù)，逐步預(yù)測B點聲速剖面，并進(jìn)行EOF分解，更新α。

(8)

(9)

式中：s為聲源的參數(shù)，Z為聲源的深度，R為水平距離；v為速度；Δt為2次測量之間的時間差；vd和va分別是聲源深度與速度的隨機(jī)變量。

則系統(tǒng)的狀態(tài)向量定義及系統(tǒng)的狀態(tài)方程和測量方程為：

(10)

xk=f(xk-1)+vk

(11)

yk=h(xk)+wk

(12)

式中：f(·)為狀態(tài)向量的狀態(tài)轉(zhuǎn)移函數(shù)；L為聲速剖面時變特性模型；F為聲源位置狀態(tài)轉(zhuǎn)移矩陣。vk和wk為過程噪聲和測量噪聲；h(·)為簡正波模型；考慮遠(yuǎn)場近似，k時刻水平距離R、深度Z處的聲源，在接受陣列處的聲壓P可以表示為：

(13)

式中：Φm(z)為本征函數(shù)，描述了第m號簡正波的模式形狀；krm為本征值，即水平波數(shù)；ρ為海水密度；M為簡正波個數(shù)；Z為聲源深度。

2.3 集合卡爾曼濾波

集合卡爾曼濾波[16]為基于蒙特卡羅采樣的卡爾曼濾波，可以解決非線性和非高斯的問題。本文在2.2節(jié)的基礎(chǔ)上，獲得由N個樣本點組成的背景集合Xk和觀測集合Zk，并計算背景集合的平均狀態(tài)分別為：

Xk={xk,i,i=1,2,…,N}

(14)

Zk={zk,i,i=1,2,…,N}

(15)

xk,i=f(xk-1)+vk

(16)

(17)

其中觀測集合在背景集合的基礎(chǔ)上構(gòu)建：

zk,i=h(xk,i)

(18)

利用最新的觀測信息，通過卡爾曼濾波更新當(dāng)前時刻的背景集合里的每個樣本點，得到分析集合：

(19)

(20)

其中Gk為卡爾曼增益：

(21)

利用分析集合的平均值對k時刻的狀態(tài)進(jìn)行更新：

(22)

3 聯(lián)合反演結(jié)果

3.1 勻速直線運動聲源仿真

假設(shè)聲源做勻速直線運動，利用KRAKEN聲場計算工具，使用實測聲速剖面數(shù)據(jù)進(jìn)行淺海波導(dǎo)下的聲場仿真。仿真波導(dǎo)的海深為84 m，海底密度為1.78 g/cm3，海底聲速為1 650 m/s，海底衰減為0.13 dB/λ。

利用LSTM聲速剖面時變特性模型進(jìn)行了聲速剖面和聲源位置的聯(lián)合反演，記為LSTM-EnKF方法；同時利用了一階時間演化模型進(jìn)行了聲速剖面和聲源位置的聯(lián)合反演，記為EnKF方法。將反演結(jié)果進(jìn)行對比。

LSTM-EnKF方法和EnKF方法下聲源參數(shù)的反演結(jié)果如圖3、圖4所示。圖3(a) 、(c) 、(e)分別為LSTM-EnKF方法下聲源深度、距離、速度的反演結(jié)果與實際的對比，圖3(b) 、(d) 、(f)為3個參數(shù)的反演誤差；圖4(a) 、(c) 、(e)分別為EnKF方法下聲源深度、距離、速度的反演結(jié)果與實際的對比，圖4(b) 、(d) 、(f)為3個參數(shù)的反演誤差?？梢园l(fā)現(xiàn)2種方法都能夠?qū)β曉磪?shù)進(jìn)行追蹤，且均展現(xiàn)了較好的追蹤性能，3個參數(shù)的誤差均在0左右。這是由于在進(jìn)行聲源追蹤時，2種方法使用了相同的狀態(tài)方程和狀態(tài)轉(zhuǎn)移矩陣。

圖3 LSTM-EnKF 方法下勻速直線運動聲源追蹤結(jié)果Fig.3 Inversion results for the source parameters using LSTM-EnKF method when the source in state of uniform motion

圖4 EnKF 方法下勻速直線運動聲源追蹤結(jié)果Fig.4 Inversion results for the source parameters using EnKF method when the source in state of uniform motion

如圖5所示，對于EOF前3階系數(shù)的追蹤結(jié)果，在第3階EOF系數(shù)的追蹤上，EnKF追蹤結(jié)果的誤差要略大于LSTM-EnKF追蹤結(jié)果的誤差。并且可以看出，使用一階時間演化模型的EnKF方法，在前30 min左右出現(xiàn)了較大的誤差，也正是EnKF方法在第3階EOF系數(shù)追蹤上，誤差較大的時刻。此外，圖4顯示EOF3的跟蹤比EOF1和EOF2系數(shù)差。與前2階相比，第3階EOF系數(shù)對SSP擾動的控制作用較小，導(dǎo)致聲速的微小變化。通過比較，驗證了LSTM-EnKF法在聲速剖面的追蹤上的可行性，LSTM-EnKF法能夠估計所有參數(shù)，且波動范圍較小，與真實值幾乎完全吻合。

圖5 勻速直線運動聲源下EOF反演結(jié)果Fig.5 Inversion results for the EOF parameters when the source in state of uniform motion

3.2 非勻速運動聲源仿真

在相同仿真參數(shù)下，利用KRAKEN聲場計算工具，使用實測聲速剖面數(shù)據(jù)，選取該海域某次淺海實驗移動聲源的實測運動參數(shù)，進(jìn)行非勻速運動聲源的淺海聲場仿真。使用LSTM-EnKF方法，進(jìn)行了聲速剖面和聲源位置的聯(lián)合反演。如圖7所示，與圖3中LSTM-EnKF方法對勻速直線運動聲源的追蹤結(jié)果相比較，對于實測聲源的追蹤，由于聲源做非勻速運動，追蹤結(jié)果的誤差略大。其中聲源速度的追蹤性能較差，這可能是由于在表1所示的仿真參數(shù)下，測量方程(18)所對應(yīng)的代價函數(shù)對聲源速度的敏感度較低。

圖6 勻速直線運動聲源下聲速剖面反演結(jié)果Fig.6 Inversion results for the sound speed profile when the source in state of uniform motion

圖7 LSTM-EnKF方法下非勻速直線運動聲源追蹤結(jié)果Fig.7 Inversion results for the source parameters using LSTM-EnKF method when the source in state of non-uniform motion

LSTM-EnKF方法下聲源參數(shù)的反演結(jié)果如圖8所示,上中下3個參數(shù)分別為聲源深度、水平距離及聲源速度，其中進(jìn)行了估計值與真實值的比較，并展示了聲源參數(shù)的反演誤差。圖8、9中為EOF 3階系數(shù)與聲速剖面的反演結(jié)果，第3階EOF系數(shù)追蹤效果較差，但由于第3階EOF系數(shù)對SSP擾動的控制作用較小，聲速剖面的時間演化過程仍能夠很好地描述；其誤差與圖6(b)相比略小，對比圖5(b)與圖8，這可能是由于第2階EOF系數(shù)對SSP擾動的控制要大于第3階系數(shù)。

圖8 非勻速直線運動聲源下EOF反演結(jié)果Fig.8 Inversion results for the EOF parameters when the source in state of non-uniform motion

圖9 非勻速直線運動聲源下聲速剖面反演結(jié)果Fig.9 Inversion results for the SSPs when the source in state of non-uniform motion

4 結(jié)論

1)使用循環(huán)神經(jīng)網(wǎng)絡(luò)對歷史水文數(shù)據(jù)進(jìn)行學(xué)習(xí)建模，能夠得到對聲速剖面的時變特性較為精確的描述。

2)在相同的仿真數(shù)據(jù)下，使用循環(huán)神經(jīng)網(wǎng)絡(luò)訓(xùn)練得到的模型作為狀態(tài)方程，相較于使用一階隨機(jī)游走過程，前者聯(lián)合反演結(jié)果更為精確，能夠追蹤聲速剖面波動。

在仿真條件中，聲源的移動軌跡較為平緩，應(yīng)使用復(fù)雜多變的運動軌跡來進(jìn)行方法的驗證；仿真過程和聯(lián)合反演過程中將淺海環(huán)境視為水平均勻環(huán)境，未考慮聲速剖面空間變化對聯(lián)合反演結(jié)果的影響；該方法需要測量點前一小時的聲速剖面數(shù)據(jù)作為先驗數(shù)據(jù)。在后續(xù)的工作中，要根據(jù)以上幾點，探尋利用深度學(xué)習(xí)，如何進(jìn)行水平非均勻環(huán)境下時變的聲速剖面的描述，并進(jìn)行相應(yīng)復(fù)雜環(huán)境下的目標(biāo)定位。