李歡 付俊 范松 陳濤 李鵬 鄭萬里
(1.寶雞石油機(jī)械有限責(zé)任公司 2.國家油氣鉆井裝備工程技術(shù)研究中心 3.四川寶石機(jī)械鉆采設(shè)備有限責(zé)任公司)
作為大洋鉆探及水合物取芯的關(guān)鍵設(shè)備,鉆柱補(bǔ)償裝置能有效減小井底鉆壓波動、減小波浪對浮式鉆井作業(yè)的影響,保障海洋鉆井作業(yè)的正常開展、提高油氣開發(fā)效率。目前,主流的鉆柱補(bǔ)償裝置分為游車、天車和死繩端3種,通常配備主動補(bǔ)償功能,并根據(jù)作業(yè)海況及補(bǔ)償精度要求進(jìn)行純被動補(bǔ)償和主被動聯(lián)合補(bǔ)償兩種作業(yè)工況的切換[1]。
現(xiàn)階段,成熟的鉆柱補(bǔ)償裝置主要供應(yīng)商均為國外公司[2]。國內(nèi)相關(guān)石油裝備企業(yè)和高校也開展了一系列樣機(jī)的試制,并取得了一定的實(shí)船應(yīng)用[3-6],但具備主動補(bǔ)償功能的鉆柱補(bǔ)償裝置還未有應(yīng)用案例。鉆柱補(bǔ)償裝置的核心技術(shù)在于主動補(bǔ)償,其控制算法的優(yōu)劣直接決定了整套鉆柱補(bǔ)償裝置的性能。傳統(tǒng)的主動補(bǔ)償控制算法通常采用基于位置(速度或加速度)反饋的多級PID閉環(huán)控制策略,該類型的反饋控制系統(tǒng)都是在船舶升沉運(yùn)動發(fā)生后,才啟動控制補(bǔ)償機(jī)構(gòu)進(jìn)行補(bǔ)償,屬于被動控制類型。頂驅(qū)、大鉤和鉆桿等游動部分具有較大的質(zhì)量,屬于大慣量、大滯后系統(tǒng),傳統(tǒng)的控制策略對于消除這種滯后并不理想,對于進(jìn)一步提高補(bǔ)償性能具有一定的局限性[7-10]。如果引入升沉預(yù)測,便可以提前獲得船舶未來的運(yùn)動趨勢,并在船舶運(yùn)動發(fā)生前就采取相應(yīng)的控制策略來補(bǔ)償船舶的運(yùn)動,這是解決大慣量系統(tǒng)滯后性問題和提高系統(tǒng)補(bǔ)償精度的一個可行性方案[11]。
對于船舶升沉預(yù)測的方法很多,主要分為基于水動力學(xué)和基于非水動力學(xué)兩種類型。基于水動力學(xué)的預(yù)測是通過分析船舶在海上的受力,推導(dǎo)出船舶的運(yùn)動狀態(tài)方程,并依據(jù)實(shí)時海況進(jìn)行預(yù)測;基于非水動力學(xué)的預(yù)測則是通過采集船舶的升沉運(yùn)動數(shù)據(jù),經(jīng)過一系列處理進(jìn)行的運(yùn)動預(yù)測。國內(nèi)學(xué)者在該領(lǐng)域進(jìn)行了大量研究,提出了基于時間序列預(yù)報(bào)、自適應(yīng)預(yù)報(bào)、ELMAN神經(jīng)網(wǎng)絡(luò)預(yù)測和卡爾曼濾波預(yù)測等多種預(yù)測方法,并且在深海起重機(jī)、波浪補(bǔ)償平臺和直升機(jī)平臺等領(lǐng)域進(jìn)行了現(xiàn)場應(yīng)用,驗(yàn)證了升沉運(yùn)動預(yù)測對于改善主動波浪補(bǔ)償性能的可行性[12-14]。本文在前人研究的基礎(chǔ)上[15],對死繩端鉆柱補(bǔ)償裝置進(jìn)行仿真研究,搭建了基于卡爾曼濾波的升沉運(yùn)動預(yù)測模型,研究了升沉預(yù)測對鉆柱補(bǔ)償裝置補(bǔ)償性能的影響,并進(jìn)一步分析了不同海況下不同預(yù)測時刻補(bǔ)償性能的影響。研究結(jié)果可為后續(xù)產(chǎn)品的應(yīng)用和控制優(yōu)化提供理論依據(jù)。
鉆柱的受力和運(yùn)動情況非常復(fù)雜,以450T死繩端補(bǔ)償裝置為研究目標(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為被動補(bǔ)償液缸的拉力,N;Fa為主動補(bǔ)償液缸的拉力(推力),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為提升裝置作用在死繩補(bǔ)償裝置上的單根鋼絲繩提升力,該力為系統(tǒng)內(nèi)力。
1—上部架體;2—下部架體;3—主動補(bǔ)償液缸;4—被動補(bǔ)償液缸;5—蓄能器;6—低壓氣瓶。
等效剛度計(jì)算公式如下:
Kd=EA/L
(2)
式中:L為鉆柱長度,取9 000 m;E為鉆柱彈性模量,取206 GPa;A為鉆柱的橫截面積,按?139.7 mm鉆桿取0.004 277 m2。
由于鉆井液的黏性及壓力,鉆柱相對于鉆井液運(yùn)動時,會導(dǎo)致能量耗散。黏滯阻力隨著鉆井液速度和密度的增大而增大,因此當(dāng)鉆井液流速或密度較大時,該阻滯作用不能忽略。鉆井液黏滯阻力系數(shù)計(jì)算公式如下:
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。
船舶在海浪中的升沉運(yùn)動可以近似地看成是一系列正弦波的合成運(yùn)動。預(yù)測算法的主體思路為檢測出一段時間內(nèi)升沉運(yùn)動的信息,并由此計(jì)算出未來的升沉運(yùn)動。因此,對于在t0~tk時間段內(nèi)的升沉位移w(t)進(jìn)行測量,經(jīng)過數(shù)據(jù)處理得到該升沉數(shù)據(jù)包含的N個正弦信號。以海平面為參考點(diǎn),其升沉位移表達(dá)式如下:
(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)對升沉運(yùn)動數(shù)據(jù)w(t)進(jìn)行處理,根據(jù)采樣定理,采樣頻率應(yīng)大于升沉運(yùn)動最大頻率的2倍以上,采樣信號數(shù)量應(yīng)在計(jì)算時間允許范圍內(nèi)盡可能大,以保證FFT處理的數(shù)據(jù)有足夠的精確度。其次,分析升沉數(shù)據(jù)經(jīng)過FFT處理之后的頻譜數(shù)據(jù),對頻譜的幅值進(jìn)行峰值檢測處理,確定最大的幾個幅值以及其對應(yīng)的相位,通常最大幅值的數(shù)量不宜超過20個,由此得到包含N個主要運(yùn)動頻率幅值和相位的一系列正弦波數(shù)據(jù)。再次,利用峰值檢測得到的各主要運(yùn)動頻率幅值和相位推導(dǎo)出狀態(tài)空間模型,搭建基于卡爾曼濾波遞推原理的觀測器。將實(shí)時升沉信號w(t)引入觀測器中,并利用卡爾曼濾波進(jìn)行誤差消除,濾波完成之后得到其校正的幅值A(chǔ)obs,i和校正相位φobs,i參數(shù)。最后,利用各個頻率的校正參數(shù)重新按照式(4)進(jìn)行正弦波疊加,通過設(shè)置預(yù)報(bào)時間可以實(shí)現(xiàn)對船舶升沉運(yùn)動的預(yù)測。
考慮到鉆柱補(bǔ)償裝置采用工控機(jī)或PLC為核心的控制系統(tǒng),數(shù)據(jù)的采集和處理均存在一定的時間周期,故觀測器的設(shè)計(jì)應(yīng)以離散控制系統(tǒng)進(jìn)行考慮。當(dāng)觀測器接收到來自峰值檢測提取到的N個正弦信號參數(shù)后,對每個信號進(jìn)行離散化處理,建立基于采樣時刻的離散化狀態(tài)方程,推導(dǎo)出相應(yīng)的狀態(tài)轉(zhuǎn)換矩陣,并進(jìn)行卡爾曼遞推,最終得到修正后的幅值和相位參數(shù)。
2.2.1 升沉運(yùn)動狀態(tài)空間模型
對于連續(xù)系統(tǒng),船舶的升沉運(yùn)動在一定時間范圍內(nèi)可以看成一系列正弦波的疊加。對這段時間的升沉信號進(jìn)行FFT和峰值檢測之后,某一時刻t0對應(yīng)的頻率為fi、幅值為Ai和相位為φi的升沉運(yùn)動的表達(dá)式:
wi(t0)=Aisin(2πfit0+φi)
(5)
針對這個正弦波位移,設(shè)定狀態(tài)變量為:
(6)
單一正弦波運(yùn)動狀態(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ù)進(jìn)行整合,得到整個升沉運(yùn)動的狀態(tài)轉(zhuǎn)換矩陣A、狀態(tài)變量X、輸出矩陣C分別為:
(10)
(11)
(12)
可以得到整個升沉運(yùn)動的狀態(tài)方程為:
(13)
由于鉆柱補(bǔ)償裝置控制系統(tǒng)為離散系統(tǒng),對連續(xù)正弦波的狀態(tài)方程應(yīng)進(jìn)行離散化處理。設(shè)控制系統(tǒng)的采樣周期為Ts,采樣初始時間為t0,任意時刻tk的表達(dá)式應(yīng)為:
tk=t0+kTs(k∈N*)
(14)
式中:N*為正整數(shù)集。
對某一頻率連續(xù)正弦波信號進(jìn)行離散化處理后,得到的狀態(tài)方程的解為:
xtk+1=eATsxtk
(15)
設(shè)定狀態(tài)轉(zhuǎn)換矩陣:
Φi=eATs
(16)
單一正弦信號經(jīng)過離散化后狀態(tài)方程表示為:
(17)
Φi為離散正弦信號狀態(tài)轉(zhuǎn)換矩陣,對其進(jìn)行矩陣變換后得到簡化計(jì)算公式為:
(18)
在FFT狀態(tài)信息未發(fā)生更新時,將式(18)帶入式(17),可以通過任意一正弦信號tk時刻狀態(tài)推算出tk+1時刻的狀態(tài)。再將式(17)根據(jù)式(11)的形式進(jìn)行組合,即可得到由N個正弦信號合并的船舶升沉運(yùn)動狀態(tài)空間模型。
2.2.2 卡爾曼濾波算法
卡爾曼濾波算法是一種能夠?qū)崿F(xiàn)快速狀態(tài)遞推濾波的方法,它結(jié)合了狀態(tài)更新和時間更新過程[16]。本文考慮到計(jì)算速度,采用線性卡爾曼濾波進(jìn)行升沉運(yùn)動預(yù)測,其遞推公式如下。
從t0開始,離散系統(tǒng)依據(jù)第tk時刻狀態(tài)計(jì)算tk+1時刻狀態(tài)變量的預(yù)測方程為:
(19)
輸出變量的估計(jì)值:
(20)
狀態(tài)變量的最優(yōu)估計(jì)方程:
(21)
卡爾曼增益方程:
(22)
預(yù)測誤差協(xié)方差更新方程:
Pk+1|k=ΦPkΦT+ΓQΓT
(23)
估計(jì)誤差協(xié)方差更新方程:
Pk+1=(1-Kk+1C)Pk+1|k
(24)
式中:Zk+1為測量的實(shí)時升沉位移,用于對狀態(tài)估計(jì)值進(jìn)行誤差校正;Γ為噪聲驅(qū)動矩陣,這里取1;Q為系統(tǒng)過程噪聲方差矩陣,這里取2N×2N的單位矩陣;R為系統(tǒng)測量噪聲方差矩陣,取-0.1;Φ為升沉運(yùn)動的狀態(tài)轉(zhuǎn)移矩陣,為N個正弦波的狀態(tài)轉(zhuǎn)移矩陣Φi根據(jù)公式(10)的形式組合而成。
卡爾曼濾波算法的具體實(shí)現(xiàn)方法如圖3所示。
圖3 卡爾曼濾波流程圖Fig.3 Kalman filtering flow chart
由圖3可知,線性卡爾曼濾波分為狀態(tài)和增益更新兩個計(jì)算回路,狀態(tài)更新回路依據(jù)每個時刻增益更新回路計(jì)算出的卡爾曼增益來計(jì)算最優(yōu)的狀態(tài)估計(jì)值。從t0開始,經(jīng)過卡爾曼k步濾波之后,得到tk時刻某一頻率的正弦波形的狀態(tài)變量為:
(25)
可以計(jì)算出各個頻率對應(yīng)的新的狀態(tài)參數(shù)為:
(26)
對所有組成升沉位移的正弦波經(jīng)過上述卡爾曼濾波校正的幅值和相位參數(shù)進(jìn)行疊加,并將校正的幅值和相位帶入式(4)進(jìn)行更新,即可得到tk的升沉運(yùn)動預(yù)報(bào)公式。根據(jù)預(yù)測時間TPred可以得到在tp=tk+TPred時刻的升沉運(yùn)動方程為:
(27)
式中:i=1,……,N;a(tp)為直流分量,且與公式(4)中的a(t)相等。
對預(yù)測效果進(jìn)行驗(yàn)證,波形方程如下:
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ù)文獻(xiàn)[15]的研究結(jié)論,升沉預(yù)測時間TPred為1 s,測量過程引入高斯白噪聲干擾,白噪聲設(shè)置為(0~1)范圍隨機(jī)數(shù)的0.05倍,根據(jù)升沉運(yùn)動實(shí)時位移得到的預(yù)測位移值如圖4所示。
圖4 升沉預(yù)測效果示意圖Fig.4 Schematic diagram of heave prediction effect
可見,由于測量噪聲的干擾,預(yù)測模型在對實(shí)時升沉運(yùn)動進(jìn)行預(yù)測時,在波浪的極值處出現(xiàn)了一定程度的波動,但從整體來看,升沉預(yù)測曲線的趨勢和幅值大小與實(shí)際運(yùn)動曲線基本保持一致,預(yù)測效果可以得到保證。
在AMESim平臺中搭建死繩鉆柱補(bǔ)償裝置仿真模型,并以式(28)對應(yīng)的波浪曲線模擬船舶升沉運(yùn)動。然后,搭建Simulink接口,采集波浪數(shù)據(jù)后在Simulink中進(jìn)行預(yù)測,并進(jìn)行聯(lián)合仿真。
以鉆井深度9 000 m、要求井底鉆壓控制在200 kN附近為仿真初始條件,控制方式采用位置PID閉環(huán)控制,以井底鉆壓波動值衡量鉆柱補(bǔ)償裝置的補(bǔ)償性能,比較未開啟預(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組最高和最低井底壓力值計(jì)算鉆壓波動,然后對5組鉆壓波動進(jìn)行均值處理,得到平均井底鉆壓波動,如表1所示。
表1 未開啟預(yù)測和開啟預(yù)測鉆壓數(shù)據(jù)Table 1 The WOB data with and without starting prediction
可見,引入1 s的升沉預(yù)測并不能減小井底鉆壓波動,無法達(dá)到改善鉆柱補(bǔ)償裝置性能的目的。這是因鉆柱補(bǔ)償裝置運(yùn)動時的固有特性導(dǎo)致,仿真中采用的液壓伺服閥0~100%階躍響應(yīng)時間為40 ms,故超前1 s的升沉信號比起主動補(bǔ)償裝置伺服液壓控制系統(tǒng)的響應(yīng)時間要長得多,主動補(bǔ)償根據(jù)提前的升沉信號進(jìn)行預(yù)補(bǔ)償,反而導(dǎo)致了井底鉆壓的紊亂,使鉆壓波動比起未開啟補(bǔ)償更高。因此,應(yīng)對預(yù)測時間對補(bǔ)償性能的影響進(jìn)行單獨(dú)研究。
在其他條件不變的前提下,對5組不同的預(yù)測時間進(jìn)行批處理計(jì)算,預(yù)測時間為0(不進(jìn)行預(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ù)計(jì)算鉆壓波動,并對每條曲線的5個鉆壓波動數(shù)據(jù)進(jìn)行均值處理,得到不同預(yù)測時間井底平均鉆壓波動情況,見表2。由圖6和表2可知,預(yù)測時間對井底鉆壓波動有顯著影響,主要體現(xiàn)在:預(yù)測時間在0.4 s時鉆壓波動最小,且隨著預(yù)測時間的延長,鉆壓波動增大;當(dāng)預(yù)測時間超過0.6 s之后,鉆壓波動比不進(jìn)行預(yù)測更大。引起該現(xiàn)象的原因是預(yù)測時間0.4 s與鉆柱補(bǔ)償裝置的響應(yīng)特性基本一致,高于或低于0.4 s將無法獲得最佳的補(bǔ)償效果。
表2 不同預(yù)測時間對鉆壓的影響Table 2 The impact of different prediction times on WOB
3.2節(jié)的研究結(jié)果是基于式(28)的海況仿真,現(xiàn)設(shè)計(jì)另外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ù)測效果明顯強(qiáng)于振幅為3.5 m之時,且預(yù)測時間0.4 s的時候補(bǔ)償效果最佳;由于海況5的升沉幅值較大,引入升沉預(yù)測對提高其補(bǔ)償效果并無明顯作用,反而會增大井底鉆壓波動??傮w來說,隨著波浪升沉幅值的增大,升沉預(yù)測對于減小井底鉆壓波動的效果逐漸減弱,且與波浪的規(guī)則與否并無太大關(guān)系。因此,鉆柱補(bǔ)償裝置在對補(bǔ)償精度要求高的場合時,建議選擇在較好的海況下開展作業(yè),且開啟升沉預(yù)測,預(yù)測時間以0.4 s為宜。
(1)提出了一種利用FFT對升沉運(yùn)動信號進(jìn)行預(yù)處理、建立狀態(tài)空間模型、利用卡爾曼濾波算法對升沉運(yùn)動進(jìn)行誤差修正,最后合成預(yù)測曲線的方法,該方法實(shí)現(xiàn)較為簡單、計(jì)算量小、實(shí)時性強(qiáng),適合鉆柱補(bǔ)償裝置這類離線平臺使用。
(2)由于鉆柱補(bǔ)償裝置和鉆桿的整體力學(xué)特性,升沉預(yù)測的效果主要受預(yù)測算法、波浪幅值、預(yù)測時間、鉆井深度、系統(tǒng)采樣頻率和邏輯處理效率等因素的影響。在只考慮波浪幅值和預(yù)測時間的前提下,升沉預(yù)測時間并不是越長越好,應(yīng)控制在0.4 s左右,此時井底鉆壓波動最小。0.4 s的預(yù)測效果只與升沉運(yùn)動的幅值有關(guān)系,幅值越高、海況越惡劣,預(yù)測對提高補(bǔ)償精度的效果越差。
(3)建議后續(xù)開展預(yù)測算法的實(shí)用化研究,依托450T死繩端補(bǔ)償裝置,編制升沉預(yù)測軟件,驗(yàn)證預(yù)測之后的補(bǔ)償效果,為下一步的現(xiàn)場應(yīng)用奠定基礎(chǔ)。