基于柔性管纜運(yùn)動(dòng)監(jiān)測(cè)數(shù)據(jù)的位移計(jì)算方法

趙晨宇，羅曉蘭，張 濤，劉軍鵬

(中國(guó)石油大學(xué)(北京) 機(jī)械與儲(chǔ)運(yùn)工程學(xué)院， 北京102249)

柔性管纜是連接水下油氣開采設(shè)備與水上浮式平臺(tái)的重要設(shè)備，有順應(yīng)性強(qiáng)、安裝和回收成本低等諸多優(yōu)點(diǎn)，在深海油氣生產(chǎn)中獲得廣泛應(yīng)用。但是，由于柔性管纜長(zhǎng)期受到風(fēng)、浪、流載荷和船體往復(fù)擺動(dòng)的擠壓作用，非常容易發(fā)生過(guò)度彎曲和疲勞破壞。為此，必須通過(guò)實(shí)時(shí)監(jiān)測(cè)的柔性管纜加速度運(yùn)動(dòng)數(shù)據(jù)，獲得柔性管纜的瞬時(shí)位移，為柔性管纜彎曲強(qiáng)度和疲勞強(qiáng)度計(jì)算提供理論依據(jù)。

通過(guò)實(shí)時(shí)監(jiān)測(cè)的柔性管纜加速度運(yùn)動(dòng)數(shù)據(jù)，求取瞬時(shí)運(yùn)動(dòng)姿態(tài)角，最終獲得其瞬時(shí)位移的常用方法主要有歐拉角法、方向余弦法和四元數(shù)算法。

1) 歐拉角法。

用1組歐拉角來(lái)描述動(dòng)態(tài)坐標(biāo)系與固定坐標(biāo)系之間的關(guān)系。該方法需要求解3個(gè)姿態(tài)角的微分方程。姿態(tài)角矩陣不需要正交化，但需用三角函數(shù)計(jì)算，工作量大[1]。當(dāng)監(jiān)測(cè)對(duì)象的俯仰角達(dá)到90°時(shí)，方程會(huì)出現(xiàn)奇異現(xiàn)象，不能進(jìn)行全姿態(tài)角求解，應(yīng)用有一定的局限性。淡鵬等針對(duì)火箭大角度運(yùn)動(dòng)下的姿態(tài)角計(jì)算問(wèn)題，應(yīng)用歐拉角法對(duì)實(shí)際案例進(jìn)行了探討和分析[2]。該方法對(duì)火箭姿態(tài)角解算和衛(wèi)星姿態(tài)分析有一定的借鑒意義，但當(dāng)火箭姿態(tài)角變化較大時(shí)，其計(jì)算結(jié)果存有較大誤差。夏喜旺等研究了飛行器大姿態(tài)角的計(jì)算問(wèn)題，推薦采用四元數(shù)法來(lái)描述飛行姿態(tài)角，以避免歐拉角法存在的奇異現(xiàn)象[3]。

2) 方向余弦法。

通過(guò)定點(diǎn)旋轉(zhuǎn)的相應(yīng)2個(gè)坐標(biāo)之間的夾角余弦來(lái)表示姿態(tài)角矩陣。該方法可用于任意姿態(tài)角的計(jì)算，但計(jì)算量大，工程上并不實(shí)用。李連仲等推導(dǎo)了方向余弦矩陣的更新遞推公式，通過(guò)仿真對(duì)其遞推公式進(jìn)行驗(yàn)證，研究并提升其公式精度[4]，可滿足捷聯(lián)慣導(dǎo)系統(tǒng)應(yīng)用中的精度要求。劉恒春提出用三角函數(shù)法求解捷聯(lián)慣導(dǎo)系統(tǒng)中的方向余弦矩陣，計(jì)算結(jié)果表明此方法精度高，但其計(jì)算過(guò)程頗為復(fù)雜[5]。劉瑞華提出新的3階最小參數(shù)法來(lái)解算方向余弦矩陣微分方程，結(jié)果表明其計(jì)算精度與經(jīng)典的3階畢卡迭代數(shù)值計(jì)算方法差不多，但姿態(tài)角解算的響應(yīng)速度即實(shí)時(shí)性顯著提高[6]。

3) 四元數(shù)法(又稱四參數(shù)法)。

通過(guò)建立的4個(gè)微分方程求解四維空間向量。運(yùn)用四元數(shù)法求解姿態(tài)角，具有精度高、計(jì)算量小、可避免奇異性的優(yōu)點(diǎn)，是目前計(jì)算姿態(tài)角的常用方法之一。四元數(shù)法需通過(guò)傳感器陀螺儀輸出柔性管纜運(yùn)動(dòng)時(shí)的角速率數(shù)據(jù)，運(yùn)用四元數(shù)角速率法進(jìn)行姿態(tài)更新。黃昊等針對(duì)陀螺儀角速率的姿態(tài)算法進(jìn)行研究，理論上推導(dǎo)了幾種航姿算法的誤差表達(dá)式，計(jì)算結(jié)果表明四元數(shù)法和4階Runge-Kutta法相結(jié)合是處理該問(wèn)題的有效算法[7]。史凱等采用四元數(shù)法和4階Runge-Kutta法研究炮彈姿態(tài)，進(jìn)行仿真平臺(tái)模擬和樣機(jī)實(shí)際試驗(yàn)，證明該算法可推廣工程應(yīng)用[8]。許毛躍等解算姿態(tài)角時(shí)，發(fā)現(xiàn)計(jì)算步長(zhǎng)對(duì)使用4階Runge-Kutta法求解四元數(shù)微分方程時(shí)的計(jì)算誤差有較大影響[9]。王德春等提出了Adams預(yù)測(cè)-校正-改進(jìn)算法，提高了捷聯(lián)慣導(dǎo)系統(tǒng)的姿態(tài)角解算精度[10-11]。

吳簡(jiǎn)彤等對(duì)比分析了四元數(shù)法、歐拉角法、方向余弦法[12]，研究結(jié)果表明：運(yùn)用四元數(shù)法進(jìn)行姿態(tài)角解算性能最佳，實(shí)用性更強(qiáng)。因而可用于柔性管纜姿態(tài)角的解算，而歐拉角法與方向余弦法因其計(jì)算量大、方程會(huì)出現(xiàn)奇異現(xiàn)象、公式復(fù)雜的問(wèn)題，均不適用于求解姿態(tài)角。

本文將利用柔性管纜的實(shí)時(shí)監(jiān)測(cè)數(shù)據(jù)，運(yùn)用四元數(shù)法和4階Runge-Kutta法進(jìn)行姿態(tài)角的計(jì)算，并應(yīng)用MATLAB軟件，最終獲得柔性管纜的實(shí)時(shí)位移數(shù)據(jù)。為此，柔性管纜實(shí)時(shí)監(jiān)測(cè)系統(tǒng)是對(duì)實(shí)時(shí)位移求解精度具有重要的作用。

1 監(jiān)測(cè)系統(tǒng)簡(jiǎn)介

1.1 監(jiān)測(cè)總體方案

柔性管纜監(jiān)測(cè)系統(tǒng)由管纜運(yùn)動(dòng)記錄儀及夾具2個(gè)部分組成。此監(jiān)測(cè)裝置安裝在防彎器下方1～2 m區(qū)域內(nèi)的柔性管纜上，如圖1所示。該區(qū)域?yàn)楦邞?yīng)力或高疲勞損傷的部位，可捕獲最大的加速度監(jiān)測(cè)數(shù)據(jù)。

圖1 監(jiān)測(cè)系統(tǒng)布局示意

1.2 監(jiān)測(cè)模塊說(shuō)明及測(cè)量數(shù)據(jù)

1.2.1 檢測(cè)模塊說(shuō)明

監(jiān)測(cè)系統(tǒng)的測(cè)量單元采用EPSON公司的G354傳感器(如圖2)，此傳感器是1個(gè)6自由度慣性測(cè)量單元，包括三軸加速度計(jì)和三軸陀螺儀。其中，三軸加速度計(jì)記錄監(jiān)測(cè)點(diǎn)的三軸線性加速度數(shù)據(jù)，三軸陀螺儀記錄監(jiān)測(cè)點(diǎn)的角速率數(shù)據(jù)，監(jiān)測(cè)數(shù)據(jù)實(shí)時(shí)存儲(chǔ)在大容量SD卡中，定期提取記錄儀存儲(chǔ)的數(shù)據(jù)。監(jiān)測(cè)模塊采用高聚能鋰電池組供電，可滿足6個(gè)月連續(xù)工作的要求。

圖2 G354傳感器

1.2.2 測(cè)量數(shù)據(jù)及處理

監(jiān)測(cè)模塊獲得的三軸線性加速度數(shù)據(jù)包括哥氏加速度、重力加速度和運(yùn)動(dòng)加速度，而柔性管纜的運(yùn)動(dòng)情況只與運(yùn)動(dòng)加速度有關(guān)，所以哥氏加速度與重力加速度作為加速度數(shù)據(jù)中的污染項(xiàng)，需被剔除。鑒于哥式加速度過(guò)小，可直接忽略。為剔除重力加速度，需要在加速度計(jì)水平和傾斜2種狀態(tài)下分析計(jì)算。

1) 加速度計(jì)水平狀態(tài)。

當(dāng)傳感器靜止時(shí)，加速度計(jì)處于水平狀態(tài)，如圖3所示。此時(shí)，X0、Y0方向的重力加速度分量為0，而Z0軸方向的重力分量為g。

圖3 三軸加速度計(jì)處于水平狀態(tài)

2) 加速度計(jì)傾斜狀態(tài)。

當(dāng)加速度計(jì)隨柔性管纜搖擺，處于傾斜狀態(tài)時(shí)。此時(shí)，加速度計(jì)各軸與Z0軸存在一定的夾角，動(dòng)坐標(biāo)系OX1Y1Z1各軸與定坐標(biāo)系OX0Y0Z0對(duì)應(yīng)各軸的夾角為φ、θ、ψ，如圖4所示。加速度計(jì)各軸在動(dòng)坐標(biāo)系內(nèi)都有重力加速度的分量。

圖4 三軸加速度計(jì)處于傾斜狀態(tài)

要剔除動(dòng)坐標(biāo)系內(nèi)各軸的重力加速度分量，需結(jié)合監(jiān)測(cè)模塊中三軸陀螺儀輸出的角速率數(shù)據(jù)，計(jì)算任意時(shí)刻動(dòng)坐標(biāo)系各軸相對(duì)于定坐標(biāo)系各軸的夾角，通過(guò)夾角確定各軸的重力加速度分量。為此，通過(guò)四元數(shù)微分方程及4階Runge-Kutta法對(duì)三軸角速率處理，以對(duì)柔性管纜任意時(shí)刻的姿態(tài)角進(jìn)行求解。

2 柔性管纜姿態(tài)角的求解

2.1 四元數(shù)表示坐標(biāo)變換

四元數(shù)的數(shù)學(xué)表達(dá)式為：

(1)

四元數(shù)能便捷地描述剛體的角運(yùn)動(dòng)，通過(guò)旋轉(zhuǎn)運(yùn)算來(lái)實(shí)現(xiàn)這種轉(zhuǎn)動(dòng)關(guān)系。針對(duì)圖4，OX0Y0Z0在經(jīng)過(guò)旋轉(zhuǎn)變換后，得到新坐標(biāo)系OX1Y1Z1，其旋轉(zhuǎn)運(yùn)算表達(dá)式為：

(2)

(3)

(4)

(5)

將式(1)和式(3)～(5)代入式(2)，進(jìn)行變換，并用矩陣的型式表達(dá)，則有：

(6)

坐標(biāo)變換系數(shù)決定著動(dòng)坐標(biāo)系相對(duì)于定坐標(biāo)系的位置，系數(shù)矩陣C實(shí)質(zhì)上就是姿態(tài)角矩陣，經(jīng)姿態(tài)角矩陣可計(jì)算運(yùn)動(dòng)狀態(tài)的姿態(tài)角，即，將矩陣式(6)簡(jiǎn)化為：

(7)

為了求解式(6)矩陣中的四元數(shù)值q0、q1、q2、q3，需建立四元數(shù)微分方程并確定四元數(shù)的初值，再應(yīng)用4階Runge-Kutta算法對(duì)四元數(shù)進(jìn)行迭代，計(jì)算出任意時(shí)刻的四元數(shù)值。

2.2 四元數(shù)微分方程的建立

四元數(shù)微分方程一般用單位四元數(shù)來(lái)推導(dǎo)，四元數(shù)描述了剛體的定點(diǎn)轉(zhuǎn)動(dòng)，單位四元數(shù)的表達(dá)式[13-14]為：

(8)

對(duì)式(8)進(jìn)行求導(dǎo)，并將α的求導(dǎo)值用角速率ω表示，得到3個(gè)任意動(dòng)坐標(biāo)的角速率矩陣為：

(9)

式中：ωx為任意動(dòng)坐標(biāo)系中X軸的角速率；ωy為任意動(dòng)坐標(biāo)系中Y軸的角速率；ωz為任意動(dòng)坐標(biāo)系中Z軸的角速率。

由式(9)看出，四元數(shù)法在進(jìn)行姿態(tài)角解算時(shí)，僅需解4個(gè)一階微分方程，運(yùn)算量顯著減少。

2.3 四元數(shù)初值確定

四元數(shù)的初始值可看作是1個(gè)初始對(duì)準(zhǔn)的過(guò)程。圖4所示，設(shè)三向定坐標(biāo)分別按照角度ψ、θ、φ轉(zhuǎn)動(dòng)3次后的四元數(shù)形式[15-16]為：

(10)

(11)

(12)

將Q1，Q2，Q3連乘得到Q，并用矩陣表達(dá)，有：

(13)

當(dāng)ψ=0，θ=0，φ=0時(shí)，即初始時(shí)動(dòng)坐標(biāo)系與定坐標(biāo)系重合，把初始夾角代入式(13)，求得此時(shí)的四元數(shù)初值為：

(14)

2.4 四元數(shù)微分方程的求解

Runge-Kutta法是一種精度較高的單步長(zhǎng)算法，在工程中廣泛使用。對(duì)于四元數(shù)微分方程式(9)，其4階Runge-Kutta法表達(dá)式為：

(15)

式中：q為四元數(shù)向量；k為斜率；h為程序迭代的步長(zhǎng)。

通過(guò)4階Runge-Kutta法可解算出四元數(shù)微分方程，由四元數(shù)初始值開始逐點(diǎn)計(jì)算，陀螺儀測(cè)得的三軸角速率ωx、ωy、ωz作為輸入條件，可計(jì)算出對(duì)應(yīng)的任意時(shí)刻的四元數(shù)值，最終求解出式(7)的姿態(tài)角矩陣。

2.5 姿態(tài)角的求解

為了將式(7)姿態(tài)角矩陣的值轉(zhuǎn)化成任意時(shí)刻的姿態(tài)角，需建立動(dòng)坐標(biāo)系的旋轉(zhuǎn)矩陣，如圖5所示，初始時(shí)動(dòng)、定坐標(biāo)系相互重合，然后將z軸固定,x軸與y軸轉(zhuǎn)動(dòng)ψ角，得動(dòng)坐標(biāo)系OX1Y1Z1；再將y軸固定,x軸與z軸轉(zhuǎn)動(dòng)θ角，得動(dòng)坐標(biāo)系OX2Y2Z2；最后，將x軸固定,y軸與z軸轉(zhuǎn)動(dòng)φ角，動(dòng)坐標(biāo)系在空間中到達(dá)位置OX3Y3Z3。

圖5 三次轉(zhuǎn)動(dòng)結(jié)果

三次轉(zhuǎn)動(dòng)后的坐標(biāo)系矩陣為：

(16)

從式(16)可以看出，轉(zhuǎn)動(dòng)后的坐標(biāo)系矩陣最終形式不僅取決于轉(zhuǎn)動(dòng)次序，還取決于每次單獨(dú)轉(zhuǎn)動(dòng)角的大小。四元數(shù)旋轉(zhuǎn)矩陣與動(dòng)坐標(biāo)系旋轉(zhuǎn)矩陣數(shù)學(xué)含義相等，四元數(shù)旋轉(zhuǎn)矩陣轉(zhuǎn)動(dòng)向量是四元素，而坐標(biāo)系旋轉(zhuǎn)矩陣轉(zhuǎn)動(dòng)的向量是對(duì)應(yīng)的坐標(biāo)軸，因此，四元數(shù)旋轉(zhuǎn)矩陣與坐標(biāo)系旋轉(zhuǎn)矩陣是等效的，則：

(17)

由式(7)和(17)得：

(18)

將四元數(shù)值代入四元數(shù)旋轉(zhuǎn)矩陣，則可求得解算姿態(tài)角所需C21、C11、C31、C32、C33各元素的值，可由式(19)求得。

(19)

應(yīng)用MATLAB軟件將上述推導(dǎo)公式編程后，可完成對(duì)應(yīng)時(shí)刻的四元數(shù)值及姿態(tài)角的求解。

3 柔性管纜位移數(shù)據(jù)求解

3.1 運(yùn)動(dòng)加速度求解

監(jiān)測(cè)模塊上測(cè)得的加速度關(guān)系式為：

(20)

重力加速度在動(dòng)坐標(biāo)系中的表達(dá)式為：

(21)

式中：Ax、Ay、Az分別為加速度計(jì)x、y、z三軸的讀數(shù)；ax、ay、az分別為動(dòng)坐標(biāo)x、y、z三軸運(yùn)動(dòng)加速度的分量；gx、gy、gz分別為動(dòng)坐標(biāo)x、y、z三軸重力加速度的分量。

通過(guò)求解出的式(18)姿態(tài)角數(shù)據(jù)，可以求解重力加速度在三軸上的分量，并在監(jiān)測(cè)模塊獲得的線性加速度數(shù)據(jù)中將其剔除。

3.2 監(jiān)測(cè)點(diǎn)位移求解

利用剔除后的動(dòng)坐標(biāo)系的加速度ax、ay、az，應(yīng)用FFT頻域計(jì)算位移，先將運(yùn)動(dòng)加速度的時(shí)域信號(hào)轉(zhuǎn)換到頻域。然后根據(jù)傅立葉變換的性質(zhì)，經(jīng)過(guò)頻域的積分和濾波后，通過(guò)傅里葉逆變換返回時(shí)域，此時(shí)得到的數(shù)組即為時(shí)域數(shù)組積分的結(jié)果。時(shí)域和頻域中信號(hào)之間的對(duì)應(yīng)關(guān)系可用離散傅里葉變換實(shí)現(xiàn)，對(duì)應(yīng)關(guān)系為：

(22)

(23)

式中：a(n)為加速度信號(hào)；X(k)為a(n)的頻譜；N為采樣數(shù)據(jù)量；n、k均為非負(fù)整數(shù)。

對(duì)信號(hào)做離散化處理，并通過(guò)傅里葉變換，則有：

(24)

(25)

式中：v(n)、d(n)分別為速度信號(hào)和位移信號(hào)的離散化表示;ωk=2πkΔf，Δf=Fs/N,為頻率分辨率。

根據(jù)上述算法，并參考王濟(jì)[17]在頻域中對(duì)振動(dòng)加速度信號(hào)的處理，設(shè)計(jì)具體計(jì)算程序，應(yīng)用MATLAB軟件計(jì)算得出柔性管纜位移結(jié)果，為柔性管纜彎曲強(qiáng)度和疲勞強(qiáng)度計(jì)算提供一定理論依據(jù)。

4 仿真結(jié)果

通過(guò)MATLAB軟件對(duì)上述公式進(jìn)行編程計(jì)算，可得以下計(jì)算結(jié)果。圖6為姿態(tài)角隨時(shí)間的變化曲線。

圖6 姿態(tài)角信號(hào)

圖7為重力加速度分量隨時(shí)間的變化曲線。

圖7 重力加速度分量信號(hào)

圖8為剔除加速度計(jì)中的重力加速度分量后，運(yùn)動(dòng)加速度隨時(shí)間的變化曲線。

圖8 加速度信號(hào)

圖9為通過(guò)FFT變換和頻域計(jì)算的方法求解出的位移隨時(shí)間變化曲線。

圖9 位移信號(hào)

由計(jì)算結(jié)果可知，重力加速度總量在單一象限內(nèi)隨時(shí)間而周期變化。其中，重力加速度在x軸方向上的分量變化最大，在z軸方向上的分量變化最小。同時(shí)，對(duì)應(yīng)的位移信號(hào)均隨時(shí)間而周期變化，其中位移信號(hào)在z軸方向上變化最大，在x軸方向上變化最小，說(shuō)明柔性管纜在垂直方向受到外載荷的影響最大。因此，在關(guān)于柔性管纜的防護(hù)工作中，應(yīng)著重考慮垂直方向上的影響。

5 結(jié)語(yǔ)

本文應(yīng)用四元數(shù)法和4階Runge-Kutta法對(duì)柔性管纜進(jìn)行姿態(tài)角更新求解，應(yīng)用FFT變換和頻域計(jì)算的方法推導(dǎo)出柔性管纜3個(gè)方向位移的計(jì)算公式，并運(yùn)用MATLAB計(jì)算其位移值。該方法具有計(jì)算量小，運(yùn)算快，程序簡(jiǎn)單等優(yōu)點(diǎn)。通過(guò)對(duì)位移數(shù)據(jù)的分析，柔性管纜在垂直方向受到外載荷的影響最大，為柔性管纜彎曲強(qiáng)度和疲勞強(qiáng)度計(jì)算提供一定理論依據(jù)，也為柔性管纜的設(shè)計(jì)與防護(hù)提供數(shù)據(jù)支持。