劉長德，顧宇翔，張進(jìn)豐

（1.中國船舶科學(xué)研究中心，江蘇無錫214082；2.中船重工奧藍(lán)托無錫軟件技術(shù)有限公司，江蘇無錫214082）

0 引 言

船舶在海上航行過程中，由于受到風(fēng)、浪、流等環(huán)境因素的干擾，運(yùn)動具有很強(qiáng)的隨機(jī)性和非線性，尤其在惡劣的海況條件下，船舶的航行和海上作業(yè)會有很大的安全隱患。如果能提前預(yù)報未來時刻的運(yùn)動狀態(tài)，獲得未來船舶運(yùn)動的暫息期，對于輔助船舶的特定作業(yè)具有重要指導(dǎo)意義。

目前，國內(nèi)外對船舶運(yùn)動極短期預(yù)報的研究一般采用離線濾波方法，利用去除高頻噪聲后的數(shù)據(jù)進(jìn)行預(yù)報結(jié)果的分析。對于船舶運(yùn)動實測數(shù)據(jù)，所采集的數(shù)據(jù)不可避免地存在高頻噪聲，從而造成信噪比嚴(yán)重降低，文獻(xiàn)[1]中以DDG-51 驅(qū)逐艦為例離線分析了測量噪聲對船舶運(yùn)動極短期預(yù)報精度的影響，尤其是在運(yùn)動峰值處會導(dǎo)致較大的預(yù)報誤差，因此進(jìn)行極短期預(yù)報時，必須對信號進(jìn)行濾波處理。傳統(tǒng)的濾波方法是將信號進(jìn)行傅里葉變換，去除掉高頻噪聲，保留有用信號，通過逆變換得到真實信號。該方法雖然能去掉噪聲，但易產(chǎn)生高頻失真，而小波變換克服了傅里葉變換中時域的瞬間變化在頻域不能反映出來的缺陷，其在低頻部分具有較高的頻率分辨率和較低的時間分辨率，在高頻部分具有較高的時間分辨率和較低的頻率分辨率，這正符合低頻信號變化緩慢，而高頻信號變化迅速的特點(diǎn)[2]。在大尺度下，可以將信號的低頻信息全局表現(xiàn)出來，在小尺度下，可以將高頻局部特性表現(xiàn)出來，其優(yōu)點(diǎn)在慣性導(dǎo)航系統(tǒng)信號分析中得到了充分體現(xiàn)[3]，這種特性使小波變換為船舶運(yùn)動信號實時濾波提供了理論基礎(chǔ)。

國內(nèi)外對船舶運(yùn)動極短期預(yù)報都非常重視并展開了許多研究，文獻(xiàn)[4]對國內(nèi)外研究現(xiàn)狀進(jìn)行了歸納總結(jié)，其中利用時間序列分析法對船舶運(yùn)動姿態(tài)進(jìn)行極短期預(yù)報越來越受到重視，早期的船舶運(yùn)動極短期預(yù)報假設(shè)船舶在海浪中的運(yùn)動姿態(tài)為一平穩(wěn)的窄帶隨機(jī)過程，采用AR 模型[5]或ARMA 模型[6-8]實現(xiàn)了船舶運(yùn)動姿態(tài)的實時預(yù)報。由于船舶運(yùn)動具有較強(qiáng)非線性、非平穩(wěn)甚至混沌特性，使得隨著預(yù)報時間的增加，其預(yù)報精度下降明顯，且運(yùn)動相位誤差變大，文獻(xiàn)[9]針對船舶運(yùn)動姿態(tài)在隨機(jī)海浪作用下的非線性特性，利用Volterra 級數(shù)的非線性表征能力，給出了船舶運(yùn)動姿態(tài)時間序列的二階Volterra 自適應(yīng)預(yù)測模型。隨著機(jī)器學(xué)習(xí)及各種算法的提出，本世紀(jì)初，國內(nèi)外學(xué)者應(yīng)用人工神經(jīng)網(wǎng)絡(luò)對船舶運(yùn)動極短期預(yù)報展開了研究，Khan 等[10]建立了BP 神經(jīng)網(wǎng)絡(luò)和遺傳算法相結(jié)合的方法，并與ARMA 預(yù)報模型進(jìn)行了比較分析，驗證了神經(jīng)網(wǎng)絡(luò)非線性建模優(yōu)點(diǎn)；Khan 等[11]結(jié)合神經(jīng)網(wǎng)絡(luò)、模糊邏輯與數(shù)據(jù)融合技術(shù)實現(xiàn)了船舶運(yùn)動極短期預(yù)報的工程應(yīng)用；文獻(xiàn)[12]結(jié)合混沌相空間重構(gòu)理論和RBF神經(jīng)網(wǎng)絡(luò)對極短期預(yù)報中建模數(shù)據(jù)、輸入變量和收斂速度等進(jìn)行了系統(tǒng)分析，為極短期預(yù)報工程化提供了理論基礎(chǔ)。隨著計算能力的提升，深度神經(jīng)網(wǎng)絡(luò)開始出現(xiàn)在時間序列預(yù)報中，得益于深度學(xué)習(xí)算法和硬件技術(shù)的突破，深度神經(jīng)網(wǎng)絡(luò)在相關(guān)領(lǐng)域攻克了許多難題，涌現(xiàn)出了許多基于機(jī)器學(xué)習(xí)算法的非線性時間序列預(yù)報模型[13-15]，其中LSTM神經(jīng)網(wǎng)絡(luò)有效解決了深度網(wǎng)絡(luò)梯度爆炸、消失問題。

綜上所述，本文針對船舶運(yùn)動非線性特性和測量噪聲對極短期預(yù)報精度的影響，在多分辨率分析理論的基礎(chǔ)上，建立了船舶運(yùn)動小波分解、閾值處理和重構(gòu)算法，并通過模型試驗數(shù)據(jù)研究了不同小波基函數(shù)和分解尺度對濾波效果的影響，進(jìn)一步結(jié)合LSTM 神經(jīng)網(wǎng)絡(luò)，建立了多步直接映射船舶運(yùn)動極短期預(yù)報方法，為船舶運(yùn)動非線性逼近建模問題提供了有效的解決途徑。

1 船舶運(yùn)動信號小波濾波

1.1 小波變換

在小波變換分析中，對應(yīng)的函數(shù)空間為L2(R)，L2(R)為R 上平方可積函數(shù)構(gòu)成的空間。設(shè)函數(shù)f(t) ∈L2(R)，其關(guān)于基小波的ψa,b(t)連續(xù)小波變換的數(shù)學(xué)表達(dá)式為

連續(xù)小波變換是冗余的，計算中不可能對所有尺度因子和平移參數(shù)小波變換，且實際的處理數(shù)據(jù)都是離散的，因此本文采用離散小波變換。

對f(t) ∈L2(R)，其離散小波變換為

通常取a0=2和b0=1，則ψk,n(t) = 2-k/2ψ(2-kt - n)， k,n ∈Z。

為使得小波變換計算更加有效，所構(gòu)造的小波函數(shù)一般具有正交性，即ψk,n(t)、ψk′,n′(t)滿足

則ψj,k(t)構(gòu)成L2(R)上規(guī)范正交基，稱為二進(jìn)制正交小波基。

1.2 多分辨率分析與Mallat算法

Mallat 提出了多分辨率的概念，在多分辨率分析理論的基礎(chǔ)上，提出了一種塔式快速小波分解算法和重構(gòu)算法，Mallat算法保證了分解后數(shù)據(jù)點(diǎn)數(shù)和分解前數(shù)據(jù)點(diǎn)數(shù)相同，是無冗余的小波變換。

對于一個多分辨分析中任意子空間Vj，Wj是Vj關(guān)于Vj+1的正交補(bǔ)空間，則Vj可分解為

進(jìn)一步可得Vj中的任意函數(shù)fj(t)都存在如下多分辨率分析表示：

小波重構(gòu)算法是基于尺度函數(shù)和小波函數(shù)的二尺度關(guān)系，通過小波逆變換，使得原函數(shù)能夠根據(jù)不同尺度的分量完全恢復(fù)?？紤]第j個分辨率下的分量和

將二尺度方程代入式（5）可得

通過多分辨率分析方法，將含噪信號進(jìn)行多尺度小波變換，從時域變換到小波域，在各尺度下提取信號的小波系數(shù)，去除含噪聲的小波系數(shù)，然后用小波逆變換重構(gòu)信號。

1.3 船舶運(yùn)動小波濾波分析

采用船舶模型試驗數(shù)據(jù)進(jìn)行濾波效果驗證分析，試驗工況為艏斜浪135o、航速6 kn。為了更好地確定分解層數(shù)，采用周期圖法對船舶橫搖、縱搖運(yùn)動信號時間歷程進(jìn)行功率譜分析。圖1為測試工況下橫搖、縱搖運(yùn)動功率譜密度，由圖1可知，當(dāng)頻率大于1.5 Hz時，其對應(yīng)功率譜密度近似為零。

圖1 運(yùn)動功率譜密度（浪向135°,V=6 kn）Fig.1 Power spectrum density of roll and pitch motion(wave direction:135°,V=6 kn)

以艏斜浪135°縱搖運(yùn)動（航速V=6 kn）數(shù)據(jù)為例，采用db4小波對船舶運(yùn)動信號進(jìn)行5層分解，圖2為對應(yīng)5層分解后各層的低頻逼近信號(C1～C5)和高頻細(xì)節(jié)信號(D1～D5)。由圖2可知，隨著分解層數(shù)的增加，細(xì)節(jié)信號逐漸增強(qiáng)。

圖2 縱搖運(yùn)動db4小波5層分解（浪向135°,V=6 kn）Fig.2 Decomposition of pitch motion with db4 wavelet(wave direction:135°,V=6 kn)

圖3為各層逼近信號與細(xì)節(jié)信號的功率譜估計，由圖3可以看出，隨著分解層數(shù)的增加，細(xì)節(jié)信號功率譜密度峰值逐漸增大，且逐漸向低頻方向移動。對于船舶運(yùn)動時歷信號，由于第1～4層細(xì)節(jié)信號主要表現(xiàn)為高頻噪聲，因此可采用硬閾值函數(shù)處理小波系數(shù)，閾值選取為0，即強(qiáng)制1～4層細(xì)節(jié)信號為0。

圖3 縱搖運(yùn)動db4小波分解信號功率譜（浪向135°,V=6 kn）Fig.3 Power spectrum density of different decomposition levels(wave direction:135°,V=6 kn)

通過對細(xì)節(jié)信號的系數(shù)閾值處理后，分別對分解后的各層逼近信號進(jìn)行重構(gòu)，圖4為采用硬閾值處理后重構(gòu)信號與原始數(shù)據(jù)的對比，并與傳統(tǒng)FFT 濾波方法進(jìn)行了比較。由圖4可知，采用db4小波基函數(shù)可有效實現(xiàn)濾波，且在運(yùn)動信號時歷端點(diǎn)處，仍能逼真地反映原始信號，而傳統(tǒng)FFT 法雖然也可有效實現(xiàn)濾波，但在端點(diǎn)附近出現(xiàn)突變，使得濾波后端點(diǎn)附近的信號與真實信號偏差較大。由于在船舶運(yùn)動極短期預(yù)報中需要當(dāng)前時刻前一段歷史數(shù)據(jù)作為樣本建立模型，當(dāng)前時刻前一段歷史數(shù)據(jù)的真實性對所建立的預(yù)報模型的預(yù)報精度至關(guān)重要，因此小波濾波可有效提高船舶運(yùn)動極短期預(yù)報模型的精度。

圖4 縱搖運(yùn)動濾波對比（浪向135°,V=12 kn）Fig.4 Filtering comparison of pitch motion (wave direction:135°,V=12 kn)

進(jìn)一步采用上述相同的閾值處理方法，表1 為采用dbN（N=2～8）系列、symN（N=2～8）系列小波基函數(shù)濾波效果的統(tǒng)計結(jié)果。由表1可知，dbN（N=3～8）、symN（N=3～8）系列小波對船舶運(yùn)動信號濾波效果相當(dāng)（FFT濾波均方差為0.005 0，信噪比為17.649 5）。

表1 不同小波基函數(shù)均方差統(tǒng)計結(jié)果Tab.1 Mean square deviation with different wavelet bases

2 LSTM 神經(jīng)網(wǎng)絡(luò)預(yù)報原理

2.1 LSTM 神經(jīng)網(wǎng)絡(luò)

長短期記憶（LSTM）神經(jīng)網(wǎng)絡(luò)是建立在RNN上的一種新型深度機(jī)器學(xué)習(xí)神經(jīng)網(wǎng)絡(luò)，通過在隱藏層各神經(jīng)單元中增加記憶單元，使時間序列上的記憶信息可控，在輸入、反饋與防止梯度爆炸之間建立了一個長時滯，該架構(gòu)強(qiáng)制其在特殊記憶單元的內(nèi)部狀態(tài)保持持續(xù)誤差流，從而使得梯度既不會爆炸也不會消失，其神經(jīng)網(wǎng)絡(luò)單元主要結(jié)構(gòu)如圖5所示。

圖5 LSTM神經(jīng)單元結(jié)構(gòu)圖Fig.5 LSTM unit structure

其每層神經(jīng)元設(shè)計為具有多個門的結(jié)構(gòu)，分別為輸入門、遺忘門和輸出門，具體計算步驟如下：

（1）輸入門，用來計算哪些信息保存到狀態(tài)單元中，包含兩部分：一部分作為當(dāng)前輸入信息保存到細(xì)胞狀態(tài)it，另一部分將當(dāng)前輸入產(chǎn)生的新信息添加到細(xì)胞狀態(tài)gt，從而產(chǎn)生新的記憶狀態(tài)ct。

（2）遺忘門，用于計算信息的遺忘程度，通過sigmoid 處理后為0 或1 的值，0 表示全部忘記，1 表示全部保留。

（3）輸出門，用于計算當(dāng)前時刻信息被輸出程度。

2.2 梯度計算

設(shè)y?i(t)、yi(t)為t時刻第i個單元的實際輸出與目標(biāo)輸出，定義損失函數(shù)E：

各參數(shù)梯度計算公式如下：

2.3 LSTM 的船舶運(yùn)動極短期預(yù)報建模流程

對于船舶運(yùn)動時間序列{x(t),t = 1,2,…,N}，對應(yīng)輸入向量為X(t - 1)，構(gòu)造訓(xùn)練樣本集，其中，X(t - 1) =[x(t - 1),x(t - 2),…,x(t - n)]T，期望輸出為x(t + L)（L=1,2…為預(yù)報步數(shù)）。采用反向傳播法對LSTM網(wǎng)絡(luò)結(jié)構(gòu)各節(jié)點(diǎn)和參數(shù)求解，其步驟如下：

（1）對各參數(shù)進(jìn)行初始化，前向計算每個神經(jīng)元的輸出值；

（2）反向計算每個神經(jīng)元的誤差項、誤差函數(shù)E對神經(jīng)元j的加權(quán)輸入的偏導(dǎo)數(shù)；

（3）計算每個權(quán)重的梯度，并對參數(shù)進(jìn)行更新；

（4）重復(fù)步驟（2）、（3），直到擬合回歸模型誤差滿足精度要求，利用該模型進(jìn)行未來時刻預(yù)報。

3 實例預(yù)報與分析

為了驗證預(yù)報模型的有效性，本文采用某船耐波性模型試驗數(shù)據(jù)進(jìn)行建模預(yù)報，試驗工況如表2所示。

表2 試驗工況Tab.2 Test conditions

利用小波濾波和LSTM 相結(jié)合的預(yù)報模型對表2中不同工況下的運(yùn)動進(jìn)行預(yù)報，模型訓(xùn)練數(shù)據(jù)個數(shù)為500，預(yù)報時間為1～15 s，連續(xù)預(yù)報次數(shù)為500次，并對預(yù)報精度（EPA）進(jìn)行統(tǒng)計分析。

式中：EPA表示預(yù)報精度；n為預(yù)報的數(shù)據(jù)個數(shù)；x(t)、x?(t)分別為t時刻的真實值與預(yù)報值；x?m(t)、xm(t)分別為n個預(yù)報值和實測值的平均期望；σx?、σx分別代表預(yù)報和期望的標(biāo)準(zhǔn)差。EPA反映了預(yù)報曲線與真實曲線之間的形狀相似程度。EPA越接近1，表明預(yù)報曲線與實測數(shù)據(jù)曲線相似程度越高。

圖6 為浪向180°、航速6 kn 工況下縱搖和垂蕩運(yùn)動預(yù)報精度隨預(yù)報時間的關(guān)系，由圖可以看出，隨著預(yù)報時間的增長，預(yù)報精度逐漸下降。對于該工況下不同預(yù)報時間的縱搖運(yùn)動預(yù)報精度要高于垂蕩運(yùn)動，縱搖運(yùn)動預(yù)報時間為11 s 時預(yù)報精度仍滿足0.85 以上，之后隨著預(yù)報時間增長，下降速度明顯加快。對于該工況下垂蕩運(yùn)動預(yù)報在預(yù)報精度大于0.85時，預(yù)報時間可達(dá)8 s。

圖7 為浪向180°、航速6 kn 工況下縱搖運(yùn)動超前4 s 和8 s 時的預(yù)報與試驗結(jié)果對比，由圖中可以看出，超前預(yù)報4 s 時，預(yù)報結(jié)果與試驗結(jié)果在相位及峰值上都吻合較好；超前預(yù)報8 s 時，在峰（谷）值較大的時刻，無明顯相位誤差，預(yù)報結(jié)果與試驗結(jié)果峰值誤差變大。

圖6 縱搖、垂蕩運(yùn)動不同預(yù)報時間對應(yīng)的預(yù)報精度（浪向180o，V=6 kn）Fig.6 Pitch and heave prediction precision versus prediction time(wave direction:180°,V=6 kn)

圖7 縱搖運(yùn)動預(yù)報結(jié)果（浪向180o，V=6 kn）Fig.7 Pitch motion prediction(wave direction:180°,V=6 kn)

進(jìn)一步對浪向180°、航速12 kn 工況下縱搖和垂蕩運(yùn)動，浪向135°、航速6 kn 工況下的橫搖、縱搖和垂蕩運(yùn)動開展了極短期預(yù)報研究。圖8 和圖9 分別為兩個工況下不同運(yùn)動預(yù)報精度隨預(yù)報時間的變化曲線。結(jié)果表明，隨著預(yù)報時間的增加，縱搖運(yùn)動預(yù)報精度均高于橫搖和垂蕩運(yùn)動的預(yù)報精度。浪向180°、航速12 kn工況下，在滿足預(yù)報精度EPA大于0.85時，縱搖運(yùn)動預(yù)報時間可達(dá)11 s，垂蕩運(yùn)動預(yù)報時間可達(dá)9 s；浪向135°、航速6 kn 工況下，預(yù)報精度EPA大于0.85 時，橫搖運(yùn)動預(yù)報時間可達(dá)8 s，縱搖運(yùn)動預(yù)報時間可達(dá)14 s，垂蕩運(yùn)動預(yù)報時間可達(dá)11 s。

圖8 縱搖、垂蕩運(yùn)動預(yù)報精度（浪向180o，V=12 kn）Fig.8 Pitch and heave prediction precision vs. prediction time(wave direction:180°,V=12 kn)

圖9 橫搖、縱搖和垂蕩運(yùn)動預(yù)報精度（浪向135o，V=6 kn）Fig.9 Roll,pitch and heave prediction precision vs.prediction time(wave direction:135°,V=6 kn)

圖10～14 分別為浪向180°、航速12 kn 工況下縱搖和垂蕩運(yùn)動，浪向135°、航速6 kn 工況下的橫搖、縱搖及垂蕩運(yùn)超前4 s和8 s時的預(yù)報結(jié)果與試驗結(jié)果對比。結(jié)果表明，超前預(yù)報4 s時，兩個工況下運(yùn)動預(yù)報結(jié)果與試驗結(jié)果具有較好的相似度，且相位與峰值無明顯差別；超前預(yù)報8 s時，預(yù)報結(jié)果峰值與試驗結(jié)果峰值誤差變大，但相位均吻合較好。

圖10 縱搖運(yùn)動預(yù)報結(jié)果（浪向180o，V=12 kn）Fig.10 Pitch motion prediction(wave direction:180°,V=12 kn)

圖11 垂蕩運(yùn)動預(yù)報結(jié)果（浪向180o，V=12 kn）Fig.11 Heave motion prediction(wave direction:180°,V=12 kn)

圖12 垂蕩運(yùn)動預(yù)報結(jié)果（浪向135°，V=6 kn）Fig.12 Heave motion prediction(wave direction:135°,V=6 kn)

圖13 橫搖運(yùn)動預(yù)報結(jié)果（浪向135°，V=6 kn）Fig.13 Roll motion prediction(wave direction:135°,V=6 kn)

圖14 縱搖運(yùn)動預(yù)報結(jié)果（浪向135°，V=6 kn）Fig.14 Pitch motion prediction(wave direction:135°,V=6 kn)

4 結(jié) 論

本文針對船舶運(yùn)動的非線性特性及測量噪聲對預(yù)報精度的影響，建立了基于小波多分辨率分析的濾波方法，有效實現(xiàn)了船舶運(yùn)動的濾波，并建立了基于LSTM 神經(jīng)網(wǎng)絡(luò)的多步直接映射船舶運(yùn)動極短期預(yù)報模型，通過耐波性試驗數(shù)據(jù)對濾波方法和預(yù)報模型性能進(jìn)行了驗證分析，得到了如下主要結(jié)論：

（1）通過對小波濾波算法中小波基、閾值處理以及尺度函數(shù)的分析，確定了dbN（N=3～8）系列、symN（N=3～8）系列小波基函數(shù)、分解層數(shù)5 層和1～4 層細(xì)節(jié)信號硬閾值處理方法，可有效去除船舶運(yùn)動信號中的高頻噪聲，解決了傳統(tǒng)FFT濾波方法對當(dāng)前時刻數(shù)據(jù)處理引起的端點(diǎn)突變問題，可有效提高船舶運(yùn)動極短期預(yù)報建模的魯棒性；

（2）通過不同航速、不同浪向下的預(yù)報結(jié)果與試驗數(shù)據(jù)的對比分析，驗證了所建立的基于LSTM神經(jīng)網(wǎng)絡(luò)的預(yù)報模型可有效實現(xiàn)船舶運(yùn)動歷史信息與未來時刻的非線性映射，橫搖、縱搖和垂蕩運(yùn)動預(yù)報幅值和相位均取得了令人滿意的精度，縱搖運(yùn)動預(yù)報時間可達(dá)10 s，橫搖和垂蕩運(yùn)動可達(dá)8 s以上。