基于PSO的船體動態變形模型參數在線辨識

于東康，楊功流，謝祖輝

（北京航空航天大學儀器科學與光電工程學院，北京100191）

0 引言

現代化大型艦船通常裝備有雷達、光學觀瞄設備、跟蹤測量設備、武器發射設備等［1］。這些設備正常工作的前提，是其擁有統一的時空基準。一般而言，時間基準通常為GPS時鐘或船上當地時鐘信號，而空間基準則是船舶在陸地建造期間將各設備的空間姿態坐標進行高精度匹配而建立起來的［2］。

由于艦船并非絕對剛體，在實際航行過程中難免受日曬、海浪沖擊等外力作用而產生角形變，這使得艦上各設備失去了統一的空間基準，影響了其信息共享與協同作戰的能力。學者Mochalov等人所做的研究表明，艦船的靜態變形角最大可達1°～1.5°，其動態撓曲變形角可達數角分。為保證設備的正常工作，需對船體形變角進行精確測量并進行補償［3?4］。

近年來，各國不斷加深在船體變形測量領域的研究，并相繼出現了一系列的測量方法。其中，美國先后提出激光跟蹤法、應變測量法等方法；俄羅斯Mochalov團隊開展了基于激光陀螺組合體（LGU）的形變角測量方法研究。國內方面，秦石喬教授提出了自準直平行光管法；長春光機所采用大鋼管法研制了一套光學形變監測系統，并進行了實船驗證；學者鄭佳興針對基于姿態匹配的變形角解算方法，進行了細致的研究。

相較于傳統的光學測量法、應變測量法等，利用多套捷聯慣導系統、基于慣性量匹配的船體變形角測量方法擁有測量精度高、易于實現等優點，是時下研究的熱點，而基于角速度匹配［5］、慣性姿態匹配［6］等的多種匹配算法也陸續出現。要得到高精度的測量結果，必須解決的一項難點就是動態撓曲變形模型參數的辨識。對此，傳統的做法是通過處理一段變形測量的歷史數據來解出相應的參數，即所謂的 “先驗數據分析法”。但實際上，受環境影響，動態模型參數是實時變化的，由先驗數據解出的既定參數無法滿足需要。文獻［7］提出了基于交互式多模型濾波的變形測量方案，通過自適應估計的方法解決了參數變化的問題。文獻［8］提出，在2套INS的角增量輸出差值中包含動態撓曲角信息，并推導得到了角增量差值的自相關函數值，其為指數衰減正弦信號（EDS）形式，最終利用基于后向預測及奇異值截斷的方法（T?K法）解出了動態模型的參數，具有較高的精度。

受文獻［8］啟發，本文利用2套INS的角增量差值自相關函數為EDS形式建立了模型，采用粒子群算法進行動態撓曲變形模型參數的在線估計，并將所估參數應用于基于角速度匹配的Kalman濾波器中，進行變形角的估計。上述方法在仿真實驗中取得了良好的求參效果，且在不同程度噪聲的干擾下均獲得了較高的求參精度。結合該參數辨識方法的船體變形角估計誤差能夠控制在10″左右，滿足了工程應用的需要。

1 角速度匹配法及動態撓曲角模型

1.1 角速度匹配方程

光纖陀螺捷聯慣導系統（FOG SINS）主要由3只光纖陀螺與3只加速度計組成，它們分別安裝在由精密機械結構支撐的3個正交面上。光纖陀螺能夠實時敏感載體繞3個正交軸方向的旋轉角速度。

假定2套INS分別安裝在艦艇的相應位置，如圖1所示。其中，INS1坐標系oxyz與船體坐標系（b 系）取齊，INS2 坐標系 ox′y′z′與艦載設備坐標系（s系）取齊。

圖1 INS安裝坐標示意圖Fig.1 Installation schematic diagram of INS

總形變角φ由2部分組成，分別為靜態形變角Φ和動態形變角θ（t），即有

在認為變形角為小角度的前提下，有

在式（4）基礎上建立量測方程，以角速度差Δω作為濾波觀測量，合理設計線性Kalman濾波器，可以實現對變形角的精確估計。

1.2 動態撓曲角模型

對于短時間的觀測而言，艦艇的靜態變形角可被視作常值，而動態撓曲變形角來自于外界的動態擾動。更進一步說，其產生的主要原因為艦體在航行過程中所受到的隨機風力和流體力的作用，因此它是一種由外界隨機力驅動的隨機變量，與由白噪聲驅動的Markov過程類似。故可以用3個方向獨立的二階Markov過程來建立模型，這種簡單的二階模型適用于大多數應用場合［9］。

動態撓曲變形的自相關函數為

式（6）中，i＝x，y，z分別代表3 個軸向，μi為衰減因子，λi為支配頻率，σi為幅度標準差。自相關函數對應的濾波器方程為

變形角的估計精度與動態撓曲模型的3個參數（σ2、μ、λ）的準確性有著直接的關系。由于艦船的航行環境具有多變性和不可預測性的特點，受其影響的動態模型參數也是實時變化的。如何精確辨識出模型參數并將其應用于Kalman濾波，是一項極具挑戰性的工作。

1.3 Kalman濾波設計

考慮光纖陀螺漂移的影響，式（4）變化為

引入陀螺漂移誤差模型

式（10）中，κi為隨機漂移不規則系數，σ′i為表征漂移不規則程度的均方差。由式（7）、式（10）可給出系統狀態方程的矩陣形式

2 基于粒子群算法的參數辨識

2.1 角增量差值自相關函數

將角速度輸出在采樣時間內進行積分，推導2套INS的角增量匹配關系。記采樣周期為t0，考慮第k次采樣，對式（4）兩邊同時積分得到

整理后，得到2套INS在第k個采樣周期內的角增量匹配方程

式（14）中，δΘ1k和 δΘ2k分別為INS1和 INS2在單采樣周期內的角增量。由φk＝Φ ＋θk，可將其進一步寫作

通過所估靜態變形角 Φ，可對式（15）中的-δΘ1k×Φ項進行在線補償，從而將其消掉，即化為

在式（17）的基礎上，文獻［8］詳細推導了 ΔΘ的自相關函數形式，這里直接給出

2.2 粒子群算法的參數辨識原理

現在面臨的工作是，如何利用式（18）辨識出其中的參數。學者吳偉提出，式（18）的形式為指數衰減正弦波（EDS）疊加的形式，即

式（23）中，信號y（n）由M個EDS信號構成，as為復振幅，w（n）為Gauss白噪聲。

基于這點，吳偉選用了一種解EDS信號的方法——T?K法，來進行參數的求解。其做法是利用RΔΘ（n）的復共軛數據建立后向預測方程，通過求解預測誤差濾波多項式的根得到未知參數。在確定預測誤差濾波多項式系數的過程中，還需用到奇異值截斷的方法。作為一種純數學手段，T?K法能夠得到較高的求參精度，但它的短板是計算量大，對信噪比的門限要求高［10?11］。

作為一種智能算法，粒子群算法（PSO）因其全局搜索能力強，被廣泛用于參數估計，并表現出了良好的估計性能［12］。PSO算法的基本思路為：首先在n維空間內初始化s個隨機粒子（隨機解），定義為X＝（X1，X2，…，Xs），并定義其中第i個粒子的位置為Xi，速度為Vi，其個體極值為pi，粒子的群體極值為BestSi。在每次迭代過程中，粒子通過式（24）、式（25）完成位置和速度的更新

式（24）中，w（t）為時變慣性權重，kg＝1，2，…，G，為當前進化代數，r1、r2為分布于區間［0，1］的隨機數，c1為局部學習因子，c2為全局學習因子。

利用PSO算法進行動態變形角模型參數的辨識，首先需要明確算法模型。回顧式（18），已知它的形式為指數衰減正弦信號的疊加，借助式（23），可將其進一步簡化為

這里，i＝1，2，…，N為數據樣本的個數，a、b、c為待辨識參數。以式（27）作為算法模型進行編碼，并將適應度函數設為如下形式

式（28）中，y（i）為數據樣本真值，J 值越小，越接近真實值，模型參數越準確。PSO算法的具體實現流程如圖2所示。

圖2 PSO算法流程圖Fig.2 Flow chart of PSO

2.3 動態變形模型參數辨識的實現流程

對利用PSO算法進行動態變形模型參數辨識的思路進行整理：首先推導了2套INS輸出角增量差值的自相關函數形式，即式（18）。其中，在簡化角增量差值ΔΘk的過程中，設定了對 -δΘ1k×Φ項進行在線補償。接著，通過對比式（18）和式（21），得出可以利用角增量差值的自相關函數對動態變形模型的參數進行辨識。選用粒子群算法，利用指數衰減正弦信號通式形式輔助建立算法模型，并合理設置了PSO參數以完成求參工作。

依照上述思路，動態變形模型參數的在線辨識流程設置如下，如圖3所示。

圖3 基于PSO算法的動態變形模型參數在線辨識流程Fig.3 Online estimation procedure of dynamic flexure model parameters based on PSO

1）求2套INS在每個采樣周期內的角增量輸出差值ΔΘk，并補償 -δΘ1k×Φ項。其中，第1次靜態變形角Φ可由粗對準大致給出，之后由變形算法估計結果給出；

2）求解 ΔΘk自相關函數值 RΔΘ（n）；

3）將 RΔΘ（n）代入粒子群算法，利用式（28）進行參數辨識；

4）將辨識得到的參數應用于角速度匹配算法中進行變形角估計，將所估靜態變形角Φ代回至步驟1，進行下一步解算。

對于PSO算法的參數設置，通過嘗試和經驗，將種群規模s設置為120；學習因子c1設為1.3，c2設為1.7；為保證收斂速度，同時避免陷入局部最優，將最大進化代數設為500；使用時變權重w（t），以確保PSO算法在全程中都擁有良好的搜索性能。

3 仿真實驗

3.1 仿真條件

設置載體以5m/s的速度做定速直航行駛，其在3個方向的角運行模型如下

式（29）中，幅值 A＝5°、B＝4°、C＝4°，角頻率ω1＝2π/7、ω2＝2π/6、ω3＝2π/8，初始相位 ψ0＝θ0＝γ0＝0。

表1 動態撓曲變形模型參數設置Table 1 Parameter setting of dynamic flexure deformation model

在仿真實驗中，其他參數設置如下。

圖4 動態撓曲變形角曲線Fig.4 Diagram of dynamic deformation angle curve

采樣頻率設為100Hz，采樣時間為10min。

3.2 仿真結果

1）考慮無噪聲情況，即動態撓曲變形信號的SNR→∞。選用20s的仿真數據進行角增量差值自相關函數RΔΘ（n）的計算，結果如圖5所示。顯然，2套INS在俯仰、橫滾、航向3個方向的RΔΘ（n）均為指數衰減正弦信號形式。截取前15s的結果，代入PSO算法進行參數辨識。對50次獨立實驗的結果求取均值，結果如表2所示。

圖5 無噪聲條件下2套INS在3個軸向上角增量差值的自相關函數曲線Fig.5 Autocorrelation function curve of angular increment output difference in three axes without noise

表2 不同信噪比情況下50次獨立求參的結果均值Table 2 Mean value of the result of 50 independent parameter estimation experiments under different SNR conditions

2）在實際情況下，受艦船發動機振動等因素的影響，動態變形信號會帶有高頻噪聲。首先考慮動態撓曲變形信號的SNR＝20dB的情況，人為地給動態撓曲信號加入噪聲，并設置信噪比為20dB，考察在受噪聲干擾情況下PSO算法的參數辨識能力。將角增量差值自相關函數曲線截取前15s結果，同樣利用PSO算法進行50次獨立求參實驗并取均值，得到的結果如表2所示。

3）考慮動態撓曲變形信號的SNR＝10dB的情況。設置信噪比為10dB，進一步考察在高噪聲干擾下的求參效果，最終結果由表2給出。

將辨識結果與真實參數進行比對，不難看出在不同程度的噪聲干擾的情況下，利用PSO算法所辨識出的動態變形模型參數均具有較高的精度。將3組參數均應用于角速度匹配Kalman濾波器中進行變形角估計，得到的估計精度相當，這里給出了在SNR＝10dB條件下得到的參數測量結果，如圖6所示。

結合PSO辨識算法之后的船體變形角估計結果與真實值十分接近，且由圖6（c）可以清楚地看到，動態變形角估計曲線能夠很好地跟隨真實值的變化情況。圖7給出了總變形角的估計誤差曲線，對3個軸向的變形角估計誤差求均方根（RMS），得到的結果為：10.11″、9.74″、10.25″。

圖6 船體變形角估計結果Fig.6 Result of ship hull deformation angle estimation

圖7 總變形角估計誤差曲線Fig.7 Estimation error curves of total deformation angle

3.3 結果分析

仿真結果證實，本文提出的基于PSO算法的辨識方法能夠有效辨識出船體動態變形模型參數，并且在不同程度的噪聲干擾下均具有較高的辨識精度。

從結合該參數辨識算法的船體變形角估計結果來看，靜態變形與動態變形曲線在3s以內即可收斂，且動態變形角曲線能夠很好地跟隨實際值，總變形角估計誤差能夠保證在10″左右，可滿足工程應用的需要。

4 結論

本文通過將2套INS的輸出角增量進行匹配，提出可以利用角增量差值進行船體動態變形模型參數的在線辨識。結合前人工作，推導了角增量差值自相關函數的具體形式，并借助指數衰減正弦信號通式建立了辨識算法模型。選擇PSO算法作為辨識工具，合理設計了算法流程，設置了相應參數，最終細化為1套完整的基于角速度匹配并結合動態變形模型參數在線辨識算法的船體變形測量流程。通過仿真結果，得到如下結論：

1）本文提出的基于PSO算法的動態變形模型參數辨識方案，在不同信噪比條件下均具有良好的辨識性能；

2）結合該在線辨識方法的船體變形測量算法，所估靜態變形角與動態變形角均能夠很好地跟隨真實值的變化，證明了參數辨識的可靠性；

3）最終的總變形角估計誤差在10″左右，且算法簡單有效，具備較高的工程應用價值。