戴洪德，柳愛利，盧建華，戴邵武，孫玉玉

（1.海軍航空工程學(xué)院 控制工程系，煙臺 264001；2.海軍航空工程學(xué)院 干部處，煙臺264001）

艦船波浪中航行時的變形分析及其IMU實時測量方法

戴洪德1，柳愛利1，盧建華1，戴邵武1，孫玉玉2

（1.海軍航空工程學(xué)院 控制工程系，煙臺 264001；2.海軍航空工程學(xué)院 干部處，煙臺264001）

大型艦船在水面上航行過程由于波浪的作用會產(chǎn)生一定的變形，該變形會降低艦船上姿態(tài)信息的精度，影響艦載武器系統(tǒng)的作戰(zhàn)效能，針對這一問題，首先對艦船的變形進行理論分析，計算艦船所承受的波浪載荷，將波浪載荷加載到艦船的有限元模型，經(jīng)計算得到的甲板各節(jié)點的角位移及線位移信息，以此來分析艦船的變形。在艦船的關(guān)鍵戰(zhàn)位點安裝由陀螺儀和加速度計組成的慣性測量單元（IMU），將IMU的輸出信息與艦船主慣導(dǎo)的輸出信息進行匹配，設(shè)計艦船變形估計的卡爾曼濾波器，實時估計IMU安裝處的甲板變形角，為艦載武器裝備提供準確的局部姿態(tài)信息，實驗結(jié)果證明了所提方法的有效性。

變形；波浪載荷；有限元；慣性測量單元；卡爾曼濾波

隨著現(xiàn)代戰(zhàn)爭對武器系統(tǒng)精確打擊能力的要求不斷提高，大型艦船上裝備的飛機、導(dǎo)彈甚至炮彈等多種武器系統(tǒng)都裝備了慣性導(dǎo)航系統(tǒng)或慣性/衛(wèi)星組合導(dǎo)航系統(tǒng)，慣導(dǎo)系統(tǒng)在進入正常工作前首先必須進行初始對準[1-2]，也就是確定慣性導(dǎo)航系統(tǒng)的初始位置、速度及姿態(tài)等信息，在艦船這樣的動基座上完成初始對準一般都采用傳遞對準的方式，當(dāng)艦船較大時，由于艦船不是一個絕對剛體，艦船海上航行時會因船體變形出現(xiàn)失調(diào)角誤差。引起這些變形誤差的因素包括日曬、老化、武器發(fā)射沖擊、轉(zhuǎn)舵操作、載荷變化等[3]，甲板縱向平面繞縱搖軸方向的變形甚至可能達幾十個角分，這樣傳遞對準的準確性會大大降低[2]，所以需要對艦船的變形進行分析并測量估計，從理論上講，艦船的撓曲變形角可以通過卡爾曼濾波進行估計[4-5]，但在實際應(yīng)用中很難建立足夠精確的模型，以消除船體變形的影響[2]。依據(jù)物理機制建立的模型受實際因素的影響，通常不具有實際應(yīng)用價值[6-7]。近年來各國的海軍、導(dǎo)航界和造船界紛紛致力于船體變形的分析、測量與校正技術(shù)的研究，以減小因船體變形引起的誤差[3，8-11]，為艦載武器裝備提供準確的局部姿態(tài)信息。

針對這一問題，本文首先介紹了一種基于有限元分析和切片理論相結(jié)合的艦船甲板變形計算方法，對艦船進行切片劃分，計算船體外表面切面所受的脈動壓力，根據(jù)與艦船切片相對應(yīng)的有限元模型節(jié)點的單位內(nèi)法矢，投影得到作用于有限元節(jié)點上的波浪載荷，通過有限元計算方法計算出各節(jié)點的應(yīng)力響應(yīng)和變形響應(yīng)；其次在艦船上關(guān)鍵戰(zhàn)位點安裝慣性測量單元(IMU：Inertial Measurement Unit），將其輸出與艦船主慣導(dǎo)的輸出進行匹配，設(shè)計卡爾曼濾波器估計艦面上IMU安裝處實時、精確的撓曲變形，并與有限元分析的數(shù)據(jù)進行對比，驗證了本文研究思路的正確性和所提方法的有效性。

2 基于有限元的艦船變形分析

2.1 艦船結(jié)構(gòu)有限元模型

與基于簡單梁理論的簡化公式相比，結(jié)構(gòu)有限元分析能更加準確地反映具體艦船結(jié)構(gòu)的特性和細節(jié)，而且能與波浪外載荷相匹配，從而確切地反映和模擬真實船體的結(jié)構(gòu)。要精確分析計算艦船在多種載荷作用下的整體彈性變形，首先就是要建立艦船的全船有限元分析模型[12-13]。選取艦船結(jié)構(gòu)中典型的甲板板架，采用三層甲板+底板結(jié)構(gòu)?？v桁和橫梁都是 T型材，縱桁兩端是艙壁，橫梁兩側(cè)是舷側(cè)，帶板因素計入橫梁和縱桁的慣性矩中。全船共設(shè)計有甲板船板的厚度，縱桁的腹板高度、厚度和面板寬度，橫梁的腹板高度、厚度和面板寬度共七個變量，結(jié)構(gòu)應(yīng)滿足局部強度和最大剛度條件的約束。所建艦船結(jié)構(gòu)的有限元模型主要使用 SHELL63板殼單元、BEAM188梁單元和MASS21點單元。SHELL63單元是具有彎曲剛度的板殼單元，能夠較好模擬艦船局部結(jié)構(gòu)的各種板殼元件。BEAM188單元能定義梁元的截面形狀，如工字梁、球扁鋼等，且適用于分析細長的梁。點單元MASS21用來對艦船進行質(zhì)量補償。

采用四點自由支持約束法來設(shè)定模型的邊界條件，在經(jīng)典結(jié)構(gòu)力學(xué)中，當(dāng)艦船靜止在水面上，一般將結(jié)構(gòu)模型的約束選為自由支持。自由支持的模型是一端限制平動，放開轉(zhuǎn)動，另一端限制兩個方向自由度平動，不限制轉(zhuǎn)動?；谶@個原理模型，將四點取在水線上，船體三點船艏一點，即在艦船水線中心處限制模型的三個平動自由度，限制船體模型外板在水線處的兩個交點的兩個方向平動自由度，而在船舷和水線相交的一點，僅限制沿和水線垂直方向上一個自由度，并且不限制轉(zhuǎn)動自由度。圖1建立了一個具有3層甲板的大型艦船的有限元模型。

圖1 船體部分網(wǎng)格劃分圖及3D效果圖Fig.1 Grid sketch of ship body and sketch of 3D model

2.2 海浪的數(shù)學(xué)描述及仿真模型

根據(jù)海浪隨機理論以及海浪譜的波能譜密度理論，采用多個（n個）隨機正、余弦波疊加來模擬海浪，數(shù)學(xué)模型如式(1)所示[14]：

式中：ζij、kij、iβ、ijω、iφ分別為單元規(guī)則波的波幅、波數(shù)、傳播方向、角頻率和相位；n為單元規(guī)則波的個數(shù)；M為波能譜的能量分布圖上頻率ω的劃分個數(shù)。

如果采用ITTC的雙參數(shù)譜（又稱ISSC波譜[2]）定義波能譜密度函數(shù)，在波浪周期T1和有義波高hs有統(tǒng)計資料供使用的情況下，波幅ζ定義為：

由線性波動理論，波數(shù)k與圓頻率ω滿足色散關(guān)系（Dispersion Relation），即

式中：g為重力加速度，h為水深。當(dāng)水深h大于波長的一半時，k=ω2/g。浪向角β顯示了單元規(guī)則波的傳播方向，頻率的分割程度是決定海浪復(fù)現(xiàn)精度的一個因素，隨機相位可以采用混合同余法，在0～2π內(nèi)隨機選取。

2.3 波浪中艦船載荷分析

根據(jù)海情和波浪譜資料在海浪仿真模型確定了中等海情下單一規(guī)則波的波幅ζ=3 m、波數(shù)k=0.252g-1、浪向角β=30°、ω=0.25 Hz、隨機相位iε后，并采用在此種海情下艦船六自由度運動的經(jīng)驗數(shù)據(jù)（表1），艦船航速設(shè)定為16 kn[14-15]。

表1 不同海況下艦船運動參數(shù)Tab.1 Motion parameters in different sea situations

應(yīng)用切片理論通過求解船舶在規(guī)則波中的運動方程式得到艦船各橫剖面所受的脈動壓力[16]P（表2）。將艦船有限元模型上與水面接觸的節(jié)點的坐標(biāo)向水面隨船平移坐標(biāo)系上作投影，得到這些節(jié)點上的單位內(nèi)法向矢量n?之后，可以得到作用于船體單位剖面內(nèi)的波浪載荷。

2.4 計算各節(jié)點運動數(shù)據(jù)

將2.3節(jié)分析得到的波浪載荷加到2.1節(jié)建立的有限元模型，加載方式如圖2所示。在仿真實驗中，以5 m為間隔，這樣沿水線X軸方向就確定了240×2=480個橫剖面，船體瞬時濕表面網(wǎng)格采用文獻[10]所提出的方法生成，即考慮瞬時物面變化引起的各種非線性因素，在時域內(nèi)建立起船體非線性運動方程。在每一時刻，認為船體以其瞬時平均濕表面為平均位置做簡諧運動。計算出甲板各節(jié)點的位置、速度、加速度、姿態(tài)、角速度等運動數(shù)據(jù)。

3 基于IMUs的變形實時估計模型

對于大型艦船，在甲板各戰(zhàn)位點處設(shè)置慣性測量單元IMU（陀螺+加速度計），通過網(wǎng)絡(luò)和綜合信息處理計算機構(gòu)建一個慣導(dǎo)網(wǎng)絡(luò)，各個IMU通過同步取數(shù)將測量信息送往綜合信息處理計算機，綜合信息處理計算機將各個 IMU提供的信息進行融合后再向所有戰(zhàn)位點發(fā)送姿態(tài)信息。當(dāng)某個局部基準發(fā)生故障時，綜合信息處理計算機可進行系統(tǒng)重構(gòu)，繼續(xù)向各戰(zhàn)位點發(fā)送姿態(tài)信息。通過這種方式可以為艦船上關(guān)鍵戰(zhàn)位點提供更豐富、準確的載艦狀態(tài)信息。

表2 x=0橫剖面的脈動壓力數(shù)據(jù)Tab.2 Pulsation pressure data of cross section x=0

圖2 艦船上的波浪載荷加載方式Fig.2 Wave load of ship

20世紀80年代以來，國內(nèi)外相繼發(fā)表了一些文獻，利用慣性測量匹配法來檢測變形失調(diào)角并消除誤差。俄羅斯圣彼得堡電工大學(xué)提出了采用兩套三軸激光陀螺組件測量艦船變形的方案[3]，利用兩慣性測量組件輸出的角速度，采用卡爾曼濾波來估算出艦船的變形角。國內(nèi)汪順亭等推導(dǎo)了三種慣性測量匹配法測量船體變形的實施方案[9]。

本課題組研究了在艦船上安裝慣性測量單元IMU，通過IMU輸出的速度增量以及角速度與艦船主慣導(dǎo)信息相匹配的方法，實時估計出 IMU安裝處實時、精確的變形角。在艦船中心位置的主慣性導(dǎo)航系統(tǒng)MINS（Main Inertial Navigation System）的坐標(biāo)系定義為x軸指向艦船的右側(cè)舷，y軸指向艦艏，z軸垂直于艦船向上，與x和y軸構(gòu)成右手直角坐標(biāo)系。甲板上IMU的坐標(biāo)系，沿著甲板的方向與MINS的坐標(biāo)系定義類似，他們的關(guān)系如圖3所示。

圖3 艦船主慣導(dǎo)與IMU相對位置及坐標(biāo)系示意圖Fig.3 Relative position and frame between the main INS of ship and IMU

其他用到的坐標(biāo)系作如下定義：n：導(dǎo)航坐標(biāo)系（艦船所在位置的當(dāng)?shù)氐乩碜鴺?biāo)系）；m：主慣性導(dǎo)航系統(tǒng)MINS坐標(biāo)系(oxyz坐標(biāo)系)；sr：真實的IMU坐標(biāo)系（ox′y′z′）；sc：計算的 IMU 坐標(biāo)系；：從MINS坐標(biāo)系到計算的 IMU坐標(biāo)系的坐標(biāo)變換方向余弦矩陣；：從MINS坐標(biāo)系到真實IMU坐標(biāo)系的坐標(biāo)變換方向余弦矩陣；：主慣導(dǎo)的姿態(tài)矩陣。

定義變形角為主慣導(dǎo)MINS坐標(biāo)系m到真實IMU坐標(biāo)系sr之間隨時間變化的誤差角度，用ψa表示（假設(shè)為小角度）；同時用ψm(t)表示m坐標(biāo)系和 sc坐標(biāo)系之間隨時間變化的誤差角度。它們之間的關(guān)系如圖4所示。

圖4 坐標(biāo)系之間關(guān)系Fig.4 Relationship among the defined frames

3.1 系統(tǒng)方程

美國學(xué)者Kain J.E和Cloputier J.R等在上世紀末提出了基于載體坐標(biāo)系內(nèi)定義的姿態(tài)誤差的全新誤差模型，介紹了基于該誤差模型的“速度+姿態(tài)角”匹配的快速對準方法[7]，很快就引起了國內(nèi)外研究者的廣泛關(guān)注和深入研究[6，17-18]，并于1998年在美國F-16戰(zhàn)機上進行了試飛驗證，在 10 s內(nèi)達到了1 mrad（3.4377′）的對準精度[6]。作者將該思想應(yīng)用到艦船變形估計中，根據(jù)實際需要重新選擇了觀測變量，推導(dǎo)了估計變形角的 Kalman濾波器量測模型，快速準確地估計出艦船的變形角。

參考kain在文獻[7]的誤差模型，以及作者前期的研究[19-21]選取速度誤差、計算姿態(tài)誤差、真實姿態(tài)誤差為濾波器設(shè)計時的狀態(tài)方程，定義狀態(tài)向量為：

誤差方程為：

式中，aη為白噪聲過程。根據(jù)式(5)可以寫出狀態(tài)方程：

3.2 觀測方程

速度誤差在文獻[7]中定義為 IMU處的的速度與經(jīng)過桿臂誤差補償?shù)闹鲬T導(dǎo)速度之差：

則k時刻的速度誤差可以表示為：

因為加速度計輸出的是速度增量的形式，所以k時刻IMU、主慣導(dǎo)的速度及桿臂速度都可以表示為k-1時刻的速度加k-1時刻到k時刻的速度增量的和，所以：

所以：

將已知量和狀態(tài)量分開整理得：

定義速度增量量測為：

可以用狀態(tài)向量表示如下：

主慣導(dǎo)的角速度測量值為：

IMU的角速度測量值為：

量測差值定義為：

所以量測模型為：

4 計算機仿真分析

應(yīng)用第2節(jié)的艦船結(jié)構(gòu)分析理論模型，得到艦船在波浪作用下的變形數(shù)據(jù)，再根據(jù)艦船上安裝的IMU測量艦船各個部分的角運動和線運動值，利用卡爾曼濾波算法估計出艦船各戰(zhàn)位點的變形角，變形角估計框圖如圖5所示。

在仿真中，艦船主慣導(dǎo)陀螺儀常值漂移為 0.02(°)/h，隨機漂移為0.01 (°)/h，艦船主慣導(dǎo)導(dǎo)航解算周期15 ms，變形估計濾波周期75 ms，仿真總時間10 s，IMU陀螺儀常值漂移為2 (°)/h，隨機漂移為1 (°)/h，刻度系數(shù)誤差2×10-4，加速度計常值偏置200 μg，隨機偏置50 μg，刻度系數(shù)誤差2×10-4，載體所在緯度為34°，經(jīng)度為108°，桿臂長度為[100；30；5]m，根據(jù)2.4節(jié)得到的IMU安裝處的速度增量、角速度數(shù)據(jù)得到IMU的測量數(shù)據(jù)，按照第3節(jié)中推導(dǎo)的慣性匹配方法估計出IMU安裝處的甲板變形角。取2.4節(jié)計算出的 IMU安裝處對應(yīng)的節(jié)點的姿態(tài)數(shù)據(jù)為變形的真實值，與kalman濾波得到的變形角估計值比較，得到艦船甲板動態(tài)變形慣性測量方法的誤差，如圖6所示。

圖5 變形角估計框圖Fig.5 Estimation diagram of deformation angle

圖6為艦船甲板9117號節(jié)點處估計得到的變形角與實際變形角的誤差圖，從仿真結(jié)果可以看出，本文所提出的基于 IMU輸出匹配的艦船甲板動態(tài)變形估計方案，對艦船變形角的跟蹤能夠在10 s內(nèi)收斂，估計誤差的具體值為[0.367;-0.740;0.723]角分。

圖6 變形角估計誤差Fig.6 Estimation error of deformation angle

5 結(jié) 論

本文首先建立了艦船結(jié)構(gòu)的有限元模型及海浪的數(shù)學(xué)模型，分析海浪作用下艦船的受力情況，計算出艦船各個節(jié)點的變形運動，分別運用IMU的陀螺儀和加速度計測量出安裝位置的角運動和線運動，與艦船主慣導(dǎo)的輸出信息匹配，應(yīng)用卡爾曼濾波器準確估計出艦船關(guān)鍵戰(zhàn)位點的姿態(tài)信息。仿真結(jié)果證明了該方法的有效性。由于艦船上位置有限，且IMU成本較高，所以下一步將研究如何根據(jù)艦船上安裝的有限IMU得到的變形數(shù)據(jù)，估計出艦船上任意位置的變形信息。

（References）：

[1]Noureldin A，Karamat T B， Georgy J. Fundamentals of inertial navigation，satellite-based positioning and their integration[M]. Berlin:Springer-Verlag Berlin and Heidelberg GmbH &Co.K，2013.

[2]Titterton D H，Weston J L.Strapdown inertial navigation technology[M].2ed.The Institution of Electrical Engineers，2004.

[3]Mochalov A C，Kazantsev A V.Use of the ring laser units for measurement of the moving object deformation[C]//Proceedings of SPIE.2002，4680:85-92.

[4]Crisan D，Rozovskii B L.The Oxford handbook of nonlinear filtering[M].Oxford:Oxford University Press，2011.

[5]Zarchan P，Musoff H.Fundamentals of Kalman filtering:a practical approach [M].Third ed.Virginia:American Institute of Aeronautics and Astronautics，2009.

[6]Shorttelle K J，Graham W R，Rabourn L C.F-16 flight tests of a rapid transfer alignment procedure[C]//Proceedings of the IEEE Position Location and Navigation Symposium，1998:379-386.

[7]Kain J E，Cloutier J R.Rapid transfer alignment for tactical weapon applications[C]//Proceedings of the AIAA Guidance，Navigation and Control Conference，1989.

[8]萬德鈞，劉玉鋒.消減艦船變形的影響和為全艦提供高精度姿態(tài)基準[J].中國慣性技術(shù)學(xué)報，2005，13(4)：77-82.WAN De-jun，LIU Yu-feng.Summary on removing influence of ship deformation and providing accurate attitude references for warship[J].Journal of Chinese Inertial Technology，2005，13(4):77-82.

[9]汪順亭，汪湛清，朱昀炤，等.船體變形的監(jiān)測方法及其對航向姿態(tài)信息的修正[J].中國慣性技術(shù)學(xué)報，2007，15(6)：635-641.WANG Shun-ting，WANG Zhan-qing，ZHU Yun-zhao，et al.Monitoring on ship hul deformation and correction for heading and attitude information [J].Journal of Chinese Inertial Technology，2007，15(6):635-641.

[10]魏學(xué)通，馬利民，莊偉東，等.一種艦船變形測量方法的研究與應(yīng)用 [J].中國慣性技術(shù)學(xué)報，2006，14(3)：74-77.WEI Xue-tong，MA Li-min，ZHUANG Wei-dong，et al.Research on distortion surveying of vessels [J].Journal of Chinese Inertial Technology，2006，14(3):74-77.

[11]周瑞.基于有限元的艦船推進軸系合理校中計算方法[J].中國艦船研究，2012，7(3)： 74-78.ZHOU Rui.Calculation method for ship propulsion shafting alignment based on finite element analysis[J].Chinese Journal of Ship Research，2012，7(3):74-78.

[12]何凱，唐文勇，羅凱，等.復(fù)合材料艦船全船有限元分析的建模方法研究[J].中國艦船研究，2011，6(5)：83-88.HE Kai，TANG Wen-yong，LUO Kai，et al.Finite element modeling method for composite ship hull[J].Chinese Journal of Ship Research，2011，6(5):83-88.

[13]Crisfield M A，Remmers J J，Verhoosel C V.Nonlinear finite element analysis of solids and structures [M].John Wiley &Sons，2012.

[14]Karnopp D C，Margolis D L，Rosenberg R C.System dynamics:Modeling，simulation，and control of mechatronic systems[M].John Wiley &Sons，2012.

[15]Cuomo G，Allsop W，Bruce T，et al.Breaking wave loads at vertical seawalls and breakwaters[J].Coastal Engineering，2010，57(4):424-439.

[16]Tarrant D，Roberts C，Jones D.Rapid and robust transfer alignment[C]//IEEE Proceedings of Aerospace Control System，1992.

[17]HAO Y，XIONG Z，WANG W，et al.Rapid transfer alignment based on unscented Kalman filter[C]//Proceedings of the 2006 American Control Conference.Minneapolis，Minnesota，USA，2006.

[18]戴洪德，周紹磊，陳明，李娟.基于四元數(shù)非線性誤差模型的快速傳遞對準[J].宇航學(xué)報，2010，31(10)：2328-2334.DAI Hong-de，ZHOU Shao-lei，CHEN Ming，LI Juan.Quaternion based nonlinear error model for rapid transfer alignment[J].Journal of Astronautics，2010，31(10):2328-2334.

[19]柳愛利，戴洪德.基于慣性傳感器輸出匹配的艦船變形估計方法[J].傳感技術(shù)學(xué)報，2011，24(1)：145-148.LIU Ai-li，DAI Hong-de.Deformation estimation of warship based on the match of inertial sensors’ output[J].Chinese Journal of Sensors and Actuators，2011，24(1):145-148.

Deformation analysis and IMU-based real-time measuring method for big ship sailing in wave

DAI Hong-de1，LIU Ai-li1，LU Jian-hua1，DAI Shao-wu1，SUN Yu-yu2
(1.Naval Aeronautical and Astronautical University，Department of Control Engineering，Yantai，264001，China；2.Naval Aeronautical and Astronautical University，Department of Political，Yantai，264001，China)

A large ship always occurs deformation when sailing in water because of the wave load，which will degrade the accuracy of attitude information on the ship and affect the efficiency of the shipboard weapons.Aiming at this problem，the ship’s deformation is analyzed theoretically，and the ship’s wave load is computed，then by loading the computed wave onto the finite element model，the angular displacement and linear displacement can be achieved by calculation.The IMU(Inertial Measurement Unit)，which is composed by gyros and accelerators，are installed in the key battle point of the ship，and the output of the IMU are matched with the main inertial navigation system of ship to estimate the deformations in real-time.Besides，a Kalman filter is designed for real-time estimating the deformation angle of the ship deck where the IMU is installed，and providing accurate locale attitude for shipboard weapons.The experiment results verify the availability of the presented method.

deformation;wave load;finite element;inertial measurement unit;Kalman filter

U666.1

A

1005-6734(2014)03-0327-06

10.13695/j.cnki.12-1222/o3.2014.03.010

2013-12-03；

2014-04-08

國家自然科學(xué)基金（61203168）；博士后（2013m532173）；航空基金（20135184007）

戴洪德（1981—），男，博士后，講師，從事慣性測量及卡爾曼濾波技術(shù)研究。E-mail：dihod@126.com