基于柔性管纜運動監測數據的位移計算方法

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

(中國石油大學(北京) 機械與儲運工程學院， 北京102249)

柔性管纜是連接水下油氣開采設備與水上浮式平臺的重要設備，有順應性強、安裝和回收成本低等諸多優點，在深海油氣生產中獲得廣泛應用。但是，由于柔性管纜長期受到風、浪、流載荷和船體往復擺動的擠壓作用，非常容易發生過度彎曲和疲勞破壞。為此，必須通過實時監測的柔性管纜加速度運動數據，獲得柔性管纜的瞬時位移，為柔性管纜彎曲強度和疲勞強度計算提供理論依據。

通過實時監測的柔性管纜加速度運動數據，求取瞬時運動姿態角，最終獲得其瞬時位移的常用方法主要有歐拉角法、方向余弦法和四元數算法。

1) 歐拉角法。

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

2) 方向余弦法。

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

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

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

吳簡彤等對比分析了四元數法、歐拉角法、方向余弦法[12]，研究結果表明：運用四元數法進行姿態角解算性能最佳，實用性更強。因而可用于柔性管纜姿態角的解算，而歐拉角法與方向余弦法因其計算量大、方程會出現奇異現象、公式復雜的問題，均不適用于求解姿態角。

本文將利用柔性管纜的實時監測數據，運用四元數法和4階Runge-Kutta法進行姿態角的計算，并應用MATLAB軟件，最終獲得柔性管纜的實時位移數據。為此，柔性管纜實時監測系統是對實時位移求解精度具有重要的作用。

1 監測系統簡介

1.1 監測總體方案

柔性管纜監測系統由管纜運動記錄儀及夾具2個部分組成。此監測裝置安裝在防彎器下方1～2 m區域內的柔性管纜上，如圖1所示。該區域為高應力或高疲勞損傷的部位，可捕獲最大的加速度監測數據。

圖1 監測系統布局示意

1.2 監測模塊說明及測量數據

1.2.1 檢測模塊說明

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

圖2 G354傳感器

1.2.2 測量數據及處理

監測模塊獲得的三軸線性加速度數據包括哥氏加速度、重力加速度和運動加速度，而柔性管纜的運動情況只與運動加速度有關，所以哥氏加速度與重力加速度作為加速度數據中的污染項，需被剔除。鑒于哥式加速度過小，可直接忽略。為剔除重力加速度，需要在加速度計水平和傾斜2種狀態下分析計算。

1) 加速度計水平狀態。

當傳感器靜止時，加速度計處于水平狀態，如圖3所示。此時，X0、Y0方向的重力加速度分量為0，而Z0軸方向的重力分量為g。

圖3 三軸加速度計處于水平狀態

2) 加速度計傾斜狀態。

當加速度計隨柔性管纜搖擺，處于傾斜狀態時。此時，加速度計各軸與Z0軸存在一定的夾角，動坐標系OX1Y1Z1各軸與定坐標系OX0Y0Z0對應各軸的夾角為φ、θ、ψ，如圖4所示。加速度計各軸在動坐標系內都有重力加速度的分量。

圖4 三軸加速度計處于傾斜狀態

要剔除動坐標系內各軸的重力加速度分量，需結合監測模塊中三軸陀螺儀輸出的角速率數據，計算任意時刻動坐標系各軸相對于定坐標系各軸的夾角，通過夾角確定各軸的重力加速度分量。為此，通過四元數微分方程及4階Runge-Kutta法對三軸角速率處理，以對柔性管纜任意時刻的姿態角進行求解。

2 柔性管纜姿態角的求解

2.1 四元數表示坐標變換

四元數的數學表達式為：

(1)

四元數能便捷地描述剛體的角運動，通過旋轉運算來實現這種轉動關系。針對圖4，OX0Y0Z0在經過旋轉變換后，得到新坐標系OX1Y1Z1，其旋轉運算表達式為：

(2)

(3)

(4)

(5)

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

(6)

坐標變換系數決定著動坐標系相對于定坐標系的位置，系數矩陣C實質上就是姿態角矩陣，經姿態角矩陣可計算運動狀態的姿態角，即，將矩陣式(6)簡化為：

(7)

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

2.2 四元數微分方程的建立

四元數微分方程一般用單位四元數來推導，四元數描述了剛體的定點轉動，單位四元數的表達式[13-14]為：

(8)

對式(8)進行求導，并將α的求導值用角速率ω表示，得到3個任意動坐標的角速率矩陣為：

(9)

式中：ωx為任意動坐標系中X軸的角速率；ωy為任意動坐標系中Y軸的角速率；ωz為任意動坐標系中Z軸的角速率。

由式(9)看出，四元數法在進行姿態角解算時，僅需解4個一階微分方程，運算量顯著減少。

2.3 四元數初值確定

四元數的初始值可看作是1個初始對準的過程。圖4所示，設三向定坐標分別按照角度ψ、θ、φ轉動3次后的四元數形式[15-16]為：

(10)

(11)

(12)

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

(13)

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

(14)

2.4 四元數微分方程的求解

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

(15)

式中：q為四元數向量；k為斜率；h為程序迭代的步長。

通過4階Runge-Kutta法可解算出四元數微分方程，由四元數初始值開始逐點計算，陀螺儀測得的三軸角速率ωx、ωy、ωz作為輸入條件，可計算出對應的任意時刻的四元數值，最終求解出式(7)的姿態角矩陣。

2.5 姿態角的求解

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

圖5 三次轉動結果

三次轉動后的坐標系矩陣為：

(16)

從式(16)可以看出，轉動后的坐標系矩陣最終形式不僅取決于轉動次序，還取決于每次單獨轉動角的大小。四元數旋轉矩陣與動坐標系旋轉矩陣數學含義相等，四元數旋轉矩陣轉動向量是四元素，而坐標系旋轉矩陣轉動的向量是對應的坐標軸，因此，四元數旋轉矩陣與坐標系旋轉矩陣是等效的，則：

(17)

由式(7)和(17)得：

(18)

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

(19)

應用MATLAB軟件將上述推導公式編程后，可完成對應時刻的四元數值及姿態角的求解。

3 柔性管纜位移數據求解

3.1 運動加速度求解

監測模塊上測得的加速度關系式為：

(20)

重力加速度在動坐標系中的表達式為：

(21)

式中：Ax、Ay、Az分別為加速度計x、y、z三軸的讀數；ax、ay、az分別為動坐標x、y、z三軸運動加速度的分量；gx、gy、gz分別為動坐標x、y、z三軸重力加速度的分量。

通過求解出的式(18)姿態角數據，可以求解重力加速度在三軸上的分量，并在監測模塊獲得的線性加速度數據中將其剔除。

3.2 監測點位移求解

利用剔除后的動坐標系的加速度ax、ay、az，應用FFT頻域計算位移，先將運動加速度的時域信號轉換到頻域。然后根據傅立葉變換的性質，經過頻域的積分和濾波后，通過傅里葉逆變換返回時域，此時得到的數組即為時域數組積分的結果。時域和頻域中信號之間的對應關系可用離散傅里葉變換實現，對應關系為：

(22)

(23)

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

對信號做離散化處理，并通過傅里葉變換，則有：

(24)

(25)

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

根據上述算法，并參考王濟[17]在頻域中對振動加速度信號的處理，設計具體計算程序，應用MATLAB軟件計算得出柔性管纜位移結果，為柔性管纜彎曲強度和疲勞強度計算提供一定理論依據。

4 仿真結果

通過MATLAB軟件對上述公式進行編程計算，可得以下計算結果。圖6為姿態角隨時間的變化曲線。

圖6 姿態角信號

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

圖7 重力加速度分量信號

圖8為剔除加速度計中的重力加速度分量后，運動加速度隨時間的變化曲線。

圖8 加速度信號

圖9為通過FFT變換和頻域計算的方法求解出的位移隨時間變化曲線。

圖9 位移信號

由計算結果可知，重力加速度總量在單一象限內隨時間而周期變化。其中，重力加速度在x軸方向上的分量變化最大，在z軸方向上的分量變化最小。同時，對應的位移信號均隨時間而周期變化，其中位移信號在z軸方向上變化最大，在x軸方向上變化最小，說明柔性管纜在垂直方向受到外載荷的影響最大。因此，在關于柔性管纜的防護工作中，應著重考慮垂直方向上的影響。

5 結語

本文應用四元數法和4階Runge-Kutta法對柔性管纜進行姿態角更新求解，應用FFT變換和頻域計算的方法推導出柔性管纜3個方向位移的計算公式，并運用MATLAB計算其位移值。該方法具有計算量小，運算快，程序簡單等優點。通過對位移數據的分析，柔性管纜在垂直方向受到外載荷的影響最大，為柔性管纜彎曲強度和疲勞強度計算提供一定理論依據，也為柔性管纜的設計與防護提供數據支持。