基于升沉預(yù)測的鉆柱補償裝置仿真研究*

李歡 付俊 范松 陳濤 李鵬 鄭萬里

(1.寶雞石油機械有限責(zé)任公司 2.國家油氣鉆井裝備工程技術(shù)研究中心 3.四川寶石機械鉆采設(shè)備有限責(zé)任公司)

0 引 言

作為大洋鉆探及水合物取芯的關(guān)鍵設(shè)備，鉆柱補償裝置能有效減小井底鉆壓波動、減小波浪對浮式鉆井作業(yè)的影響，保障海洋鉆井作業(yè)的正常開展、提高油氣開發(fā)效率。目前，主流的鉆柱補償裝置分為游車、天車和死繩端3種，通常配備主動補償功能，并根據(jù)作業(yè)海況及補償精度要求進行純被動補償和主被動聯(lián)合補償兩種作業(yè)工況的切換[1]。

現(xiàn)階段，成熟的鉆柱補償裝置主要供應(yīng)商均為國外公司[2]。國內(nèi)相關(guān)石油裝備企業(yè)和高校也開展了一系列樣機的試制，并取得了一定的實船應(yīng)用[3-6]，但具備主動補償功能的鉆柱補償裝置還未有應(yīng)用案例。鉆柱補償裝置的核心技術(shù)在于主動補償，其控制算法的優(yōu)劣直接決定了整套鉆柱補償裝置的性能。傳統(tǒng)的主動補償控制算法通常采用基于位置(速度或加速度)反饋的多級PID閉環(huán)控制策略，該類型的反饋控制系統(tǒng)都是在船舶升沉運動發(fā)生后，才啟動控制補償機構(gòu)進行補償，屬于被動控制類型。頂驅(qū)、大鉤和鉆桿等游動部分具有較大的質(zhì)量，屬于大慣量、大滯后系統(tǒng)，傳統(tǒng)的控制策略對于消除這種滯后并不理想，對于進一步提高補償性能具有一定的局限性[7-10]。如果引入升沉預(yù)測，便可以提前獲得船舶未來的運動趨勢，并在船舶運動發(fā)生前就采取相應(yīng)的控制策略來補償船舶的運動，這是解決大慣量系統(tǒng)滯后性問題和提高系統(tǒng)補償精度的一個可行性方案[11]。

對于船舶升沉預(yù)測的方法很多，主要分為基于水動力學(xué)和基于非水動力學(xué)兩種類型。基于水動力學(xué)的預(yù)測是通過分析船舶在海上的受力，推導(dǎo)出船舶的運動狀態(tài)方程，并依據(jù)實時海況進行預(yù)測；基于非水動力學(xué)的預(yù)測則是通過采集船舶的升沉運動數(shù)據(jù)，經(jīng)過一系列處理進行的運動預(yù)測。國內(nèi)學(xué)者在該領(lǐng)域進行了大量研究，提出了基于時間序列預(yù)報、自適應(yīng)預(yù)報、ELMAN神經(jīng)網(wǎng)絡(luò)預(yù)測和卡爾曼濾波預(yù)測等多種預(yù)測方法，并且在深海起重機、波浪補償平臺和直升機平臺等領(lǐng)域進行了現(xiàn)場應(yīng)用，驗證了升沉運動預(yù)測對于改善主動波浪補償性能的可行性[12-14]。本文在前人研究的基礎(chǔ)上[15]，對死繩端鉆柱補償裝置進行仿真研究，搭建了基于卡爾曼濾波的升沉運動預(yù)測模型，研究了升沉預(yù)測對鉆柱補償裝置補償性能的影響，并進一步分析了不同海況下不同預(yù)測時刻補償性能的影響。研究結(jié)果可為后續(xù)產(chǎn)品的應(yīng)用和控制優(yōu)化提供理論依據(jù)。

1 模型概述

1.1 結(jié)構(gòu)及動力學(xué)方程

鉆柱的受力和運動情況非常復(fù)雜，以450T死繩端補償裝置為研究目標(biāo)，其結(jié)構(gòu)和受力示意圖如圖1a所示。裝置的下部架體連接至鉆臺面，上部架體通過6根鋼絲繩繞過頂部滑輪與頂驅(qū)和鉆桿相連。將整個鉆柱受拉部分簡化為彈簧-質(zhì)量-阻尼系統(tǒng)，其數(shù)學(xué)模型如圖1b所示。根據(jù)瑞利法則，考慮彈簧質(zhì)量對固有頻率的影響，將彈性鉆柱總質(zhì)量的等效為Md作用在鉆柱的最底部。將鉆柱、頂驅(qū)和大鉤等游動部分的質(zhì)量用一個集中質(zhì)量MT替代，負(fù)載動力學(xué)方程如式(1)所示：

(1)

式中：Fp為被動補償液缸的拉力，N；Fa為主動補償液缸的拉力(推力)，N；Gs為上部架體的靜載荷，取147 kN；MT為集中質(zhì)量，取397 187.5 kg；Md為鉆桿等效質(zhì)量，鉆桿的質(zhì)量為114 062.5 kg；Kd為鉆柱的等效剛度，N/m；Xd為上部架體的絕對位移，m；Cd為鉆柱在鉆井液中的黏滯阻力系數(shù)，N/(m·s-1)；g為重力加速度，9.8 m/s2；Ts為提升裝置作用在死繩補償裝置上的單根鋼絲繩提升力，該力為系統(tǒng)內(nèi)力。

1—上部架體；2—下部架體；3—主動補償液缸；4—被動補償液缸；5—蓄能器；6—低壓氣瓶。

1.2 鉆柱等效剛度

等效剛度計算公式如下：

Kd=EA/L

(2)

式中：L為鉆柱長度，取9 000 m；E為鉆柱彈性模量，取206 GPa；A為鉆柱的橫截面積，按?139.7 mm鉆桿取0.004 277 m2。

1.3 鉆井液黏滯阻力系數(shù)

由于鉆井液的黏性及壓力，鉆柱相對于鉆井液運動時，會導(dǎo)致能量耗散。黏滯阻力隨著鉆井液速度和密度的增大而增大，因此當(dāng)鉆井液流速或密度較大時，該阻滯作用不能忽略。鉆井液黏滯阻力系數(shù)計算公式如下：

Cd=πCNC1ρdDvdL/8

(3)

式中：CN為摩阻系數(shù)，水基鉆井液的摩阻系數(shù)取0.1；C1為附加質(zhì)量系數(shù)，置于流體中的圓柱體取1；ρd為鉆井液密度，取1 200 kg/m3；D為鉆柱直徑，取0.139 7 m；vd為鉆井液流速，取0.3 m/s。

2 升沉預(yù)測模型研究

2.1 船舶升沉預(yù)測原理

船舶在海浪中的升沉運動可以近似地看成是一系列正弦波的合成運動。預(yù)測算法的主體思路為檢測出一段時間內(nèi)升沉運動的信息，并由此計算出未來的升沉運動。因此，對于在t0～tk時間段內(nèi)的升沉位移w(t)進行測量，經(jīng)過數(shù)據(jù)處理得到該升沉數(shù)據(jù)包含的N個正弦信號。以海平面為參考點，其升沉位移表達式如下：

(4)

式中：i=1,……,N；t0≤t≤tk；Ai、fi和φi分別為各個正弦波形所對應(yīng)幅值、頻率和相位；a(t)為直流分量，t為時間。

升沉預(yù)測模型如圖2所示。

圖2 升沉預(yù)測模型原理示意圖Fig.2 Schematic diagram of the principle of heave prediction model

首先，利用快速傅里葉變換(FFT)對升沉運動數(shù)據(jù)w(t)進行處理，根據(jù)采樣定理，采樣頻率應(yīng)大于升沉運動最大頻率的2倍以上，采樣信號數(shù)量應(yīng)在計算時間允許范圍內(nèi)盡可能大，以保證FFT處理的數(shù)據(jù)有足夠的精確度。其次，分析升沉數(shù)據(jù)經(jīng)過FFT處理之后的頻譜數(shù)據(jù)，對頻譜的幅值進行峰值檢測處理，確定最大的幾個幅值以及其對應(yīng)的相位，通常最大幅值的數(shù)量不宜超過20個，由此得到包含N個主要運動頻率幅值和相位的一系列正弦波數(shù)據(jù)。再次，利用峰值檢測得到的各主要運動頻率幅值和相位推導(dǎo)出狀態(tài)空間模型，搭建基于卡爾曼濾波遞推原理的觀測器。將實時升沉信號w(t)引入觀測器中，并利用卡爾曼濾波進行誤差消除，濾波完成之后得到其校正的幅值A(chǔ)obs,i和校正相位φobs,i參數(shù)。最后，利用各個頻率的校正參數(shù)重新按照式(4)進行正弦波疊加，通過設(shè)置預(yù)報時間可以實現(xiàn)對船舶升沉運動的預(yù)測。

2.2 基于卡爾曼濾波算法的觀測器設(shè)計

考慮到鉆柱補償裝置采用工控機或PLC為核心的控制系統(tǒng)，數(shù)據(jù)的采集和處理均存在一定的時間周期，故觀測器的設(shè)計應(yīng)以離散控制系統(tǒng)進行考慮。當(dāng)觀測器接收到來自峰值檢測提取到的N個正弦信號參數(shù)后，對每個信號進行離散化處理，建立基于采樣時刻的離散化狀態(tài)方程，推導(dǎo)出相應(yīng)的狀態(tài)轉(zhuǎn)換矩陣，并進行卡爾曼遞推，最終得到修正后的幅值和相位參數(shù)。

2.2.1 升沉運動狀態(tài)空間模型

對于連續(xù)系統(tǒng)，船舶的升沉運動在一定時間范圍內(nèi)可以看成一系列正弦波的疊加。對這段時間的升沉信號進行FFT和峰值檢測之后，某一時刻t0對應(yīng)的頻率為fi、幅值為Ai和相位為φi的升沉運動的表達式：

wi(t0)=Aisin(2πfit0+φi)

(5)

針對這個正弦波位移，設(shè)定狀態(tài)變量為：

(6)

單一正弦波運動狀態(tài)方程可表示為：

(7)

式中：Ai(t0)、Ci(t0)和wi(t0)分別為狀態(tài)方程的狀態(tài)轉(zhuǎn)換矩陣、輸出轉(zhuǎn)換矩陣和輸出矩陣。

根據(jù)正弦波公式，可以推導(dǎo)出其狀態(tài)轉(zhuǎn)換矩陣和輸出轉(zhuǎn)換矩陣分別為：

(8)

(9)

用Ai、Ci和分別代替Ai(t0)、Ci(t0)，對于經(jīng)過FFT變換后峰值檢測得到的N組升沉數(shù)據(jù)進行整合，得到整個升沉運動的狀態(tài)轉(zhuǎn)換矩陣A、狀態(tài)變量X、輸出矩陣C分別為：

(10)

(11)

(12)

可以得到整個升沉運動的狀態(tài)方程為：

(13)

由于鉆柱補償裝置控制系統(tǒng)為離散系統(tǒng)，對連續(xù)正弦波的狀態(tài)方程應(yīng)進行離散化處理。設(shè)控制系統(tǒng)的采樣周期為Ts，采樣初始時間為t0，任意時刻tk的表達式應(yīng)為：

tk=t0+kTs(k∈N*)

(14)

式中：N*為正整數(shù)集。

對某一頻率連續(xù)正弦波信號進行離散化處理后，得到的狀態(tài)方程的解為：

xtk+1=eATsxtk

(15)

設(shè)定狀態(tài)轉(zhuǎn)換矩陣：

Φi=eATs

(16)

單一正弦信號經(jīng)過離散化后狀態(tài)方程表示為：

(17)

Φi為離散正弦信號狀態(tài)轉(zhuǎn)換矩陣，對其進行矩陣變換后得到簡化計算公式為：

(18)

在FFT狀態(tài)信息未發(fā)生更新時，將式(18)帶入式(17)，可以通過任意一正弦信號tk時刻狀態(tài)推算出tk+1時刻的狀態(tài)。再將式(17)根據(jù)式(11)的形式進行組合，即可得到由N個正弦信號合并的船舶升沉運動狀態(tài)空間模型。

2.2.2 卡爾曼濾波算法

卡爾曼濾波算法是一種能夠?qū)崿F(xiàn)快速狀態(tài)遞推濾波的方法，它結(jié)合了狀態(tài)更新和時間更新過程[16]。本文考慮到計算速度，采用線性卡爾曼濾波進行升沉運動預(yù)測，其遞推公式如下。

從t0開始，離散系統(tǒng)依據(jù)第tk時刻狀態(tài)計算tk+1時刻狀態(tài)變量的預(yù)測方程為：

(19)

輸出變量的估計值：

(20)

狀態(tài)變量的最優(yōu)估計方程：

(21)

卡爾曼增益方程：

(22)

預(yù)測誤差協(xié)方差更新方程：

Pk+1|k=ΦPkΦT+ΓQΓT

(23)

估計誤差協(xié)方差更新方程：

Pk+1=(1-Kk+1C)Pk+1|k

(24)

式中：Zk+1為測量的實時升沉位移，用于對狀態(tài)估計值進行誤差校正；Γ為噪聲驅(qū)動矩陣，這里取1；Q為系統(tǒng)過程噪聲方差矩陣，這里取2N×2N的單位矩陣；R為系統(tǒng)測量噪聲方差矩陣，取-0.1；Φ為升沉運動的狀態(tài)轉(zhuǎn)移矩陣，為N個正弦波的狀態(tài)轉(zhuǎn)移矩陣Φi根據(jù)公式(10)的形式組合而成。

卡爾曼濾波算法的具體實現(xiàn)方法如圖3所示。

圖3 卡爾曼濾波流程圖Fig.3 Kalman filtering flow chart

由圖3可知，線性卡爾曼濾波分為狀態(tài)和增益更新兩個計算回路，狀態(tài)更新回路依據(jù)每個時刻增益更新回路計算出的卡爾曼增益來計算最優(yōu)的狀態(tài)估計值。從t0開始，經(jīng)過卡爾曼k步濾波之后，得到tk時刻某一頻率的正弦波形的狀態(tài)變量為：

(25)

可以計算出各個頻率對應(yīng)的新的狀態(tài)參數(shù)為：

(26)

對所有組成升沉位移的正弦波經(jīng)過上述卡爾曼濾波校正的幅值和相位參數(shù)進行疊加，并將校正的幅值和相位帶入式(4)進行更新，即可得到tk的升沉運動預(yù)報公式。根據(jù)預(yù)測時間TPred可以得到在tp=tk+TPred時刻的升沉運動方程為：

(27)

式中：i=1,……,N；a(tp)為直流分量，且與公式(4)中的a(t)相等。

2.3 升沉預(yù)測效果

對預(yù)測效果進行驗證，波形方程如下：

w(t)=0.5sin(2π×0.1t)+0.6sin(2π×0.08t)+

0.3sin(2π×0.075t)+sin(2π×0.07t)

(28)

系統(tǒng)的采樣周期Ts取0.05 s，升沉數(shù)據(jù)采樣數(shù)量為2 000個，卡爾曼濾波的步數(shù)k為20步，根據(jù)文獻[15]的研究結(jié)論，升沉預(yù)測時間TPred為1 s，測量過程引入高斯白噪聲干擾，白噪聲設(shè)置為(0～1)范圍隨機數(shù)的0.05倍，根據(jù)升沉運動實時位移得到的預(yù)測位移值如圖4所示。

圖4 升沉預(yù)測效果示意圖Fig.4 Schematic diagram of heave prediction effect

可見，由于測量噪聲的干擾，預(yù)測模型在對實時升沉運動進行預(yù)測時，在波浪的極值處出現(xiàn)了一定程度的波動，但從整體來看，升沉預(yù)測曲線的趨勢和幅值大小與實際運動曲線基本保持一致，預(yù)測效果可以得到保證。

3 仿真結(jié)果分析

在AMESim平臺中搭建死繩鉆柱補償裝置仿真模型，并以式(28)對應(yīng)的波浪曲線模擬船舶升沉運動。然后，搭建Simulink接口，采集波浪數(shù)據(jù)后在Simulink中進行預(yù)測，并進行聯(lián)合仿真。

3.1 預(yù)測效果分析

以鉆井深度9 000 m、要求井底鉆壓控制在200 kN附近為仿真初始條件，控制方式采用位置PID閉環(huán)控制，以井底鉆壓波動值衡量鉆柱補償裝置的補償性能，比較未開啟預(yù)測(預(yù)測時間0 s)和開啟預(yù)測(預(yù)測時間1 s)的升沉預(yù)測鉆壓波動情況，結(jié)果如圖5所示。

圖5 開啟預(yù)測和未開啟預(yù)測井底鉆壓波動情況Fig.5 Bottomhole WOB fluctuations with and without starting prediction

考慮到波形的不規(guī)則性，在橫坐標(biāo)的20、35、50、65和80 s附近取5組最高和最低井底壓力值計算鉆壓波動，然后對5組鉆壓波動進行均值處理，得到平均井底鉆壓波動，如表1所示。

表1 未開啟預(yù)測和開啟預(yù)測鉆壓數(shù)據(jù)Table 1 The WOB data with and without starting prediction

可見，引入1 s的升沉預(yù)測并不能減小井底鉆壓波動，無法達到改善鉆柱補償裝置性能的目的。這是因鉆柱補償裝置運動時的固有特性導(dǎo)致，仿真中采用的液壓伺服閥0～100%階躍響應(yīng)時間為40 ms，故超前1 s的升沉信號比起主動補償裝置伺服液壓控制系統(tǒng)的響應(yīng)時間要長得多，主動補償根據(jù)提前的升沉信號進行預(yù)補償，反而導(dǎo)致了井底鉆壓的紊亂，使鉆壓波動比起未開啟補償更高。因此，應(yīng)對預(yù)測時間對補償性能的影響進行單獨研究。

3.2 預(yù)測時間影響

在其他條件不變的前提下，對5組不同的預(yù)測時間進行批處理計算，預(yù)測時間為0(不進行預(yù)測)、0.2、0.4、0.6和0.8 s，仿真得到的井底鉆壓波動情況如圖6所示。

圖6 不同預(yù)測時間的鉆壓波動情況Fig.6 The effect of prediction times on the WOB fluctuation

按照表1的方法對每條曲線取5組峰值數(shù)據(jù)計算鉆壓波動，并對每條曲線的5個鉆壓波動數(shù)據(jù)進行均值處理，得到不同預(yù)測時間井底平均鉆壓波動情況，見表2。由圖6和表2可知，預(yù)測時間對井底鉆壓波動有顯著影響，主要體現(xiàn)在：預(yù)測時間在0.4 s時鉆壓波動最小，且隨著預(yù)測時間的延長，鉆壓波動增大；當(dāng)預(yù)測時間超過0.6 s之后，鉆壓波動比不進行預(yù)測更大。引起該現(xiàn)象的原因是預(yù)測時間0.4 s與鉆柱補償裝置的響應(yīng)特性基本一致，高于或低于0.4 s將無法獲得最佳的補償效果。

表2 不同預(yù)測時間對鉆壓的影響Table 2 The impact of different prediction times on WOB

3.3 海況的影響

3.2節(jié)的研究結(jié)果是基于式(28)的海況仿真，現(xiàn)設(shè)計另外5種海況，分別研究海況對預(yù)測效果的影響。海況波形曲線表見表3。

表3 海況波形曲線表Table 3 Sea state waveform curve table

表3中，海況1為比式(28)所示更為惡劣的海況，海況2和海況3比起海況1振幅更加均勻，海況4和海況5為相同頻率不同幅值的正弦波。分別將這5種海況帶入仿真模型，得到不同海況、不同預(yù)測時間對井底鉆壓波動的影響，如圖7～圖11及表4所示。

圖7 海況1預(yù)測時間影響情況Fig.7 The effect of prediction times on the WOB fluctuation under sea state 1

圖8 海況2預(yù)測時間影響情況Fig.8 The effect of prediction times on the WOB fluctuation under sea state 2

圖9 海況3預(yù)測時間影響情況Fig.9 The effect of prediction times on the WOB fluctuation under sea state 3

圖10 海況4預(yù)測時間影響情況Fig.10 The effect of prediction times on the WOB fluctuation under sea state 4

圖11 海況5預(yù)測時間影響情況Fig.11 The effect of prediction times on the WOB fluctuation under sea state 5

表4 不同海況和預(yù)測時間下平均鉆壓波動表Table 4 Average WOB fluctuation table under different sea conditions and prediction times

對比表2和表4的數(shù)據(jù)可以看出：由于海況1比式(28)對應(yīng)的海況更加惡劣，所以其井底鉆壓波動也更大；海況1在預(yù)測時間為0.2和0.4 s時，井底鉆壓波動比未開啟預(yù)測有小幅下降，但其余預(yù)測時刻對于改善井底鉆壓波動意義并不大；對于海況2和海況3，其波浪曲線比海況1振幅更加均勻，且未出現(xiàn)大幅度的突變，海況3的井底鉆壓波動明顯優(yōu)于海況2，出現(xiàn)該現(xiàn)象的原因在于海況2的平均振幅要小于海況3；對于海況4和海況5，其波形為理想的正弦波信號，振幅1.5 m時的預(yù)測效果明顯強于振幅為3.5 m之時，且預(yù)測時間0.4 s的時候補償效果最佳；由于海況5的升沉幅值較大，引入升沉預(yù)測對提高其補償效果并無明顯作用，反而會增大井底鉆壓波動。總體來說，隨著波浪升沉幅值的增大，升沉預(yù)測對于減小井底鉆壓波動的效果逐漸減弱，且與波浪的規(guī)則與否并無太大關(guān)系。因此，鉆柱補償裝置在對補償精度要求高的場合時，建議選擇在較好的海況下開展作業(yè)，且開啟升沉預(yù)測，預(yù)測時間以0.4 s為宜。

4 結(jié)論及建議

(1)提出了一種利用FFT對升沉運動信號進行預(yù)處理、建立狀態(tài)空間模型、利用卡爾曼濾波算法對升沉運動進行誤差修正，最后合成預(yù)測曲線的方法，該方法實現(xiàn)較為簡單、計算量小、實時性強，適合鉆柱補償裝置這類離線平臺使用。

(2)由于鉆柱補償裝置和鉆桿的整體力學(xué)特性，升沉預(yù)測的效果主要受預(yù)測算法、波浪幅值、預(yù)測時間、鉆井深度、系統(tǒng)采樣頻率和邏輯處理效率等因素的影響。在只考慮波浪幅值和預(yù)測時間的前提下，升沉預(yù)測時間并不是越長越好，應(yīng)控制在0.4 s左右，此時井底鉆壓波動最小。0.4 s的預(yù)測效果只與升沉運動的幅值有關(guān)系，幅值越高、海況越惡劣，預(yù)測對提高補償精度的效果越差。

(3)建議后續(xù)開展預(yù)測算法的實用化研究，依托450T死繩端補償裝置，編制升沉預(yù)測軟件，驗證預(yù)測之后的補償效果，為下一步的現(xiàn)場應(yīng)用奠定基礎(chǔ)。