劉珊珊，趙亦朋，王小秋，汪志明

(1.中國(guó)石油大學(xué)(北京)石油工程學(xué)院，北京 102249；2.中國(guó)石油集團(tuán) 工程技術(shù)研究院有限公司，北京 102200)

引 言

在鉆井軌道設(shè)計(jì)中，經(jīng)常將井眼軌道假設(shè)為某些簡(jiǎn)單的空間曲線，如空間圓弧等，這些模型形式簡(jiǎn)單，求解方便，因此應(yīng)用廣泛，能夠解決大多數(shù)的井眼軌道設(shè)計(jì)問題[1-5]。然而針對(duì)造斜同時(shí)伴隨方位變化的鉆井軌道，需要更復(fù)雜的軌道模型。此外，當(dāng)鉆井與設(shè)計(jì)軌跡不一致時(shí)，使用空間圓弧等模型不是實(shí)現(xiàn)自動(dòng)化鉆井的最佳選擇，由于其不連續(xù)性，導(dǎo)致鉆柱上產(chǎn)生顯著的應(yīng)力，扭矩和阻力增加，鉆井所需的能量也相應(yīng)增加。自然曲線井眼軌跡在鉆柱上的應(yīng)力、扭矩和阻力相對(duì)圓弧軌跡要小，能量成本明顯低于圓弧軌跡，從而降低了鉆井能耗。不同的軌跡模型在同一曲線上產(chǎn)生不同的能量消耗，有研究表明自然曲線法為最佳的軌道模型[6-9]。知名的鉆井軌道設(shè)計(jì)軟件Compass中的Build/Turn設(shè)計(jì)模型，即采用自然曲線模型。

魯港等人[10]對(duì)自然曲線模型的數(shù)值計(jì)算進(jìn)行了研究，給出了自然曲線模型在給定造斜率和轉(zhuǎn)向率的條件下，用泰勒展開求解的井段ΔX、ΔY和ΔZ的通用求解方法，未涉及約束條件下的模型反演。何樹山等人[11-12]給出了一個(gè)模型反演的基本流程，但是其初值的選擇不能保證在大斜度的情況下迭代收斂，難以保證算法的穩(wěn)定性。

本文在上述研究的基礎(chǔ)上，針對(duì)自然曲線軌道模型定點(diǎn)和垂深對(duì)齊問題的反演計(jì)算，提出了新的計(jì)算流程，并針對(duì)垂深對(duì)齊問題提出了新的初值選擇方法，通過多個(gè)算例進(jìn)行驗(yàn)算，較商業(yè)化軟件具有更高的計(jì)算精度。

1 設(shè)計(jì)模型

自然曲線法描述的是一個(gè)空間三維曲線(圖1),該曲線假設(shè)井段內(nèi)的井斜變化率和方位變化率分別為常數(shù)[6]，則

圖1 自然曲線示意圖

Δα=αB-αA，

(1)

Δφ=φB-φA，

(2)

ΔL=LA-LB，

(3)

(4)

(5)

式中：(LA,αA,φA)及(LB,αB,φB)分別為井段上、下兩個(gè)測(cè)點(diǎn)的測(cè)深、井斜角和方位角；Kα為井斜角變化率(造斜率)；Kφ為方位角變化率(轉(zhuǎn)向率)。

1.1 軌道參數(shù)計(jì)算

自然曲線模型，從A點(diǎn)以造斜率Kα和轉(zhuǎn)向率Kφ鉆進(jìn)ΔL長(zhǎng)度的井段到達(dá)B點(diǎn)，各個(gè)軌道參數(shù)的計(jì)算分別為[6]

αB=αA+KαΔL；

(6)

φB=φA+KφΔL；

(7)

(8)

(9)

(10)

(11)

式中：ΔN、ΔE、ΔZ分別為井段在北、東以及垂深方向的增量；ΔS為軌道在水平面上投影的曲線長(zhǎng)度。式(8)—(11)中，對(duì)于Kα=0和Kφ=0的特殊情況參見文獻(xiàn)[6]。

魯港等人[10]通過定義一個(gè)輔助函數(shù)λ(x)對(duì)式(8)—(11)進(jìn)行簡(jiǎn)化得

(12)

對(duì)式(12)在x=0處進(jìn)行泰勒展開得

(13)

式(13)中第5、6兩項(xiàng)為高階小量，工程上取前4項(xiàng)即可滿足計(jì)算精度要求。借助函數(shù)λ(x)，可使用

(14)

(15)

(16)

(17)

統(tǒng)一形式計(jì)算。其中，

(18)

a=λ(Δα+Δφ),b=λ(Δα-Δφ),c=λ(Δα)(其中Δα=αB-αA，Δφ=φB-φA)，

(19)

(20)

1.2 約束條件下的模型反演問題

本文將探討兩種約束條件下的反演問題：

(1)定點(diǎn)問題：給定ΔN、ΔE和ΔZ，反演求解ΔL、Kα和Kφ。

(2)垂深對(duì)齊問題：給定ΔN、ΔE和ΔZ，并要求曲線段終點(diǎn)到達(dá)指定的垂深后，穩(wěn)斜到目標(biāo)點(diǎn)，反演求解造斜段長(zhǎng)度ΔL1、穩(wěn)斜段長(zhǎng)度ΔL2、Kα和Kφ。

兩種反演問題如圖2、圖3所示。現(xiàn)有研究中，很少說明自然曲線模型的反演求解算法，本文針對(duì)這兩種問題，給出詳細(xì)穩(wěn)定的求解算法。

圖2 第一種約束條件

圖3 第二種約束條件

2 定點(diǎn)問題的求解方法

設(shè)井段起始點(diǎn)和終止點(diǎn)為A和B，已知條件為(LA,αA,φA,NA,EA,ZA)和(NB,EB,ZB)。模型求解流程如下：

(1)計(jì)算AB連線的斜率

(21)

(22)

若αA,B>0則該井段為增斜，若αA,B<0則該井段為降斜，若αA,B=0則為穩(wěn)斜。

(2)給αB一個(gè)初值，可取為[11]

(23)

根據(jù)A點(diǎn)方位角和坐標(biāo)偏移量，求解方位變化方向和方位角初值。根據(jù)φA建立變換矩陣M，對(duì)ΔN和ΔE進(jìn)行變換，確定沿方位角φA方向的坐標(biāo)增量以及方位角的變化初值Δφ0。

(24)

(25)

(26)

(4)計(jì)算井段長(zhǎng)度ΔL

(27)

(5)求解方位角φB

根據(jù)式(8)、(9)可得

(28)

式(8)中，只有φB一個(gè)未知數(shù)，但由于公式極為復(fù)雜，無法得到φB的解析表達(dá)式，通過計(jì)算分析該式在φB的臨域范圍內(nèi)函數(shù)具有單調(diào)性。因此，可以使用牛頓迭代法，快速計(jì)算φB的值。對(duì)式(28)求導(dǎo)得

(29)

(30)

需要注意的是，在式(29)、(30)中，當(dāng)出現(xiàn)兩個(gè)分子項(xiàng)|Δα±Δφ|→0的情況，由于計(jì)算機(jī)的浮點(diǎn)計(jì)算誤差，這時(shí)的計(jì)算結(jié)果不準(zhǔn)確，當(dāng)其等于零時(shí)，將會(huì)出現(xiàn)除零的錯(cuò)誤。在迭代的過程中，這一情況難以避免，因此算法必須考慮|Δα±Δφ|→0的極限問題，在此臨域內(nèi)，用極限代替計(jì)算，避免計(jì)算錯(cuò)誤的發(fā)生。

(6)求解井斜角αB

根據(jù)式(8)、式(9)，使用ΔN和ΔE均可計(jì)算ΔL，考慮到數(shù)值計(jì)算有一定誤差，浮點(diǎn)運(yùn)算數(shù)越小，誤差越大，可根據(jù)ΔN和ΔE絕對(duì)值的大小，選擇絕對(duì)值較大的項(xiàng)計(jì)算當(dāng)前的ΔL。根據(jù)式(15)推導(dǎo)得

(31)

圖4 定點(diǎn)問題計(jì)算流程圖

3 垂深對(duì)齊問題的求解方法

該問題分為兩個(gè)井段(圖3)，第一個(gè)井段是自然曲線段，通過在該段造斜和轉(zhuǎn)向，在指定的垂深(TVD)對(duì)齊到目標(biāo)點(diǎn)，并由此點(diǎn)開始穩(wěn)斜到目標(biāo)點(diǎn)，即第二個(gè)井段為直線段。根據(jù)自然曲線段和直線段的坐標(biāo)偏移公式，建立坐標(biāo)增量公式

ΔL2sinαBcosφB；

(32)

ΔL2sinαBsinφB；

(33)

(34)

(35)

式中：ΔN、ΔE、ΔZ1、ΔZ2、αA、φA為已知參數(shù)；曲線段長(zhǎng)度ΔL1、直線段長(zhǎng)度ΔL2、直線段井斜角αB和直線段方位角φB為未知參數(shù)。

該問題的求解與定點(diǎn)問題的求解思路類似，但因?yàn)榇嬖谥本€段，所以稍微有些復(fù)雜。具體計(jì)算流程如圖4(為敘述的一致性，仍假設(shè)首末兩點(diǎn)為A、B)。

(1)計(jì)算AB連線的斜率

同定點(diǎn)問題(1)。

(2)給定αB初值

圖5 垂深對(duì)齊的初值問題

通過觀察，初值應(yīng)在AB連線和AM′連線之間，取∠BAM′的角平分線作為初值，即

(36)

同定點(diǎn)問題(3)。

(4)計(jì)算曲線段長(zhǎng)度ΔL1

(37)

(5)求解方位角φB

直線段的斜率

(38)

整理可得

f(φB)=sinφB(ΔN-ΔN1)-cosφB(ΔE-ΔE1)=0。

(39)

該函數(shù)在求解區(qū)間的圖形如圖6所示。雖然該區(qū)域內(nèi)只有一個(gè)解，但函數(shù)并不保持單調(diào)性，有間斷，因此不能采用迭代法求解，只能按照方位角變化方向，使用搜索算法求解。

圖6 式(39)曲線形態(tài)

(6)求解井斜角αB

根據(jù)式(32)和式(33)，得

f(αB)=(ΔN1+ΔN2-ΔN)2+(ΔE1+ΔE2-

ΔE)2=0。

(40)

圖臨域內(nèi)的函數(shù)值

圖臨域內(nèi)的導(dǎo)數(shù)值

f′(αB)=2(ΔN1+ΔN2-ΔN)(ΔN1+ΔN2)′+2(ΔE1+ΔE2-ΔE)(ΔE1+ΔE2)′=0；

(41)

(ΔN1+ΔN2)′=

(42)

(ΔE1+ΔE2)′=

(43)

4 計(jì)算實(shí)例

根據(jù)上述兩種條件下的求解方法，編制了相應(yīng)的軟件，以下給出兩種不同類型問題的計(jì)算實(shí)例，由于現(xiàn)有文獻(xiàn)中提出的初始井斜角計(jì)算方法不適用于大斜度垂深對(duì)齊問題，因此，專門提供了一個(gè)大斜度算例，用以驗(yàn)證本文算法的適用性，同時(shí)與知名的商業(yè)軟件Compass進(jìn)行對(duì)比。

(1)算例1 定點(diǎn)問題

起始點(diǎn)A：測(cè)深=1 000 m,井斜角=30°，方位角=0°，垂深=1 000 m，北坐標(biāo)=0 m，東坐標(biāo)=0 m；靶點(diǎn)B點(diǎn)：垂深=1 100 m，北坐標(biāo)=100 m，東坐標(biāo)=50 m，本文與Compass軟件計(jì)算結(jié)果對(duì)比見表1。

表1 算例1計(jì)算結(jié)果

(2)算例2 垂深對(duì)齊問題

起始點(diǎn)A：測(cè)深=1 000 m,井斜角=30°，方位角=0°，垂深=1 000 m，北坐標(biāo)=0 m，東坐標(biāo)=0 m；造斜段末點(diǎn)垂深=1 080 m；靶點(diǎn)B點(diǎn)，垂深=1 120 m，北坐標(biāo)=100 m，東坐標(biāo)=50 m，本文與Compass軟件計(jì)算結(jié)果見表2。

表2 算例2計(jì)算結(jié)果

(3)算例3 垂深對(duì)齊問題(大斜度井)

起始點(diǎn)A：測(cè)深=1 000 m,井斜角=30°，方位角=0°，垂深=1 000 m，北坐標(biāo)=0 m，東坐標(biāo)=0 m；造斜段末點(diǎn)垂深=1 100 m；靶點(diǎn)B點(diǎn)，垂深=1 120 m，北坐標(biāo)=400 m，東坐標(biāo)=200 m，本文與Compass軟件計(jì)算結(jié)果見表3。

表3 算例3計(jì)算結(jié)果

表1—表3中，下劃線的數(shù)據(jù)為計(jì)算出的關(guān)鍵參數(shù)，分別為段長(zhǎng)、造斜率和轉(zhuǎn)向率，根據(jù)這些關(guān)鍵參數(shù)和垂深約束條件，計(jì)算出每個(gè)井段的北坐標(biāo)和東坐標(biāo)等參數(shù)。

由表1可見，本文算法的最大誤差為1.0×10-5，Compass計(jì)算的最大誤差為1.13×10-3，本文算法精度高于Compass。由表2、表3可見，本文算法與Compass計(jì)算的結(jié)果均較為精確，在1.0×10-5數(shù)量級(jí)沒有誤差，可認(rèn)為精度一致。綜上所述，針對(duì)定點(diǎn)問題本文算法計(jì)算精度高于Compass計(jì)算精度，針對(duì)垂深對(duì)齊問題精度一致。當(dāng)然，從工程角度來看，兩者計(jì)算結(jié)果均能滿足需求。

需要指出的是，前文給出了嚴(yán)格的數(shù)學(xué)公式，有利于理解問題的本質(zhì)。但是在實(shí)際的軟件中，完全采用這些公式，特別是導(dǎo)數(shù)公式，往往使程序過于繁瑣，因此在具體實(shí)現(xiàn)過程中，可以使用式(14)—(17)計(jì)算函數(shù)及差商實(shí)現(xiàn)求解，同樣可以達(dá)到收斂的效果，簡(jiǎn)化程序的編寫。

5 結(jié) 論

造斜段伴隨方位變化的自然曲線模型，是一種非常規(guī)的三維鉆井軌道設(shè)計(jì)模型，可減少鉆井摩阻。本文詳細(xì)討論了兩種約束條件下的模型反演求解問題，分別給出了兩種問題的求解流程。計(jì)算實(shí)例證明，針對(duì)定點(diǎn)問題本文算法計(jì)算精度高于Compass計(jì)算精度，針對(duì)垂深對(duì)齊問題精度一致。由于商業(yè)軟件價(jià)格昂貴，并且每年需繳納大量維護(hù)費(fèi)用，因此研究國(guó)產(chǎn)替代軟件非常必要，該算法應(yīng)用于鉆井工程設(shè)計(jì)軟件，對(duì)于實(shí)現(xiàn)國(guó)產(chǎn)軟件替代進(jìn)口軟件具有重要的意義。