1数学模型
约定：除非特别指明，具有长度量纲的参数其单位为m，角度的单位为弧度(rad)。井眼轨道由多个设计井段“稳斜-圆弧-稳斜-圆弧-稳斜”连接而成，每个设计井段依序记为1、2等等。带有下标的物理量表示该下标所属设计井段的物理量。1.1设计方程组数学模型如下[2]：ΔU为靶点的位移矢量；ΔL1、ΔL3、ΔL5为3个稳斜段的段长；t1、t3、t5为3个稳斜段的单位方向矢量；K2、K4为两个圆弧段的井眼曲率；ε2、ε4为两个圆弧段的弯曲角；λ2、λ4为两个圆弧段的辅助参数；t1·t3等表示矢量的数积[6]。单位方向矢量ti=(li，mi，ni)由下式定义(i=1，3，5)：式中：li、mi、ni是单位方向矢量ti在3个坐标轴上的分量；αi、准i为稳斜段井斜角和方位角。2.2无量纲化式中：ΔL为靶点位移矢量的模长，│·│表示矢量求模运算[6]；ΔU为无量纲化的靶点位移矢量；p、q、r为Δu在三个坐标轴上的分量；x1、x2、x3、x4、x5、k2、k4为无量纲化参数。方程组(1)～(5)可以改写成下面的无量纲化形式：1.3未知数的组合情况待定参数个数为13，包括x1、x2、x3、x4、x5、k2、k4、ε2、ε4、α3、α5、准3、准5，而独立方程个数为7，所以必须指定其中的6个为已知设计参数。待定参数x2、x4、ε2、ε4、准3为未知数，α3为已知数；如果限定入靶方向，则α5、准5为已知数；还需要指定另外3个待定参数为已知数，共有10种情况，见表1。
2特征函数
求解多元方程组的一般策略[7]是尽量使用数学化简方法将多个未知数从方程组中化简掉，形成一个包含最少未知数的新的方程组。当k2、k4为已知设计参数时，使用正切半角公式，从式(10)～(13)可得：x2=1k21-t1·t31+t1·t3姨(14)x4=1k41-t3·t51+t3·t5姨(15)显然，x2和x4是未知数准3的单变量函数。如果设计轨迹是三维的，则t1、t3、t5为线性不相关的[6]，令：则矩阵M为非异矩阵。将式(9)写成下面的形式：(x1+x2，x2+x3+x4，x4+x5)·M=Δu(16)再用矩阵M的逆矩阵M-1右乘式(16)的两端，得：(x1+x2，x2+x3+x4，x4+x5)=Δv(17)其中：Δv=Δu·M-1=(p′，q′，r)(18)易知，Δv是未知数准3的单变量矢量函数，p′、q′、r′都是未知数准3的单变量函数。式(17)，当x5、x1、x3分别为已知数时，分别有：F1(准3)=x4+x5-r′=0(19)F2(准3)=x1+x2-p′=0(20)F3(准3)=x2+x3+x4-q′=0(21)函数Fj(准3)是关于未知数准3的单变量函数(j=1，2，3)，如果能从式(19)、式(20)或式(21)求出未知数准3，则进一步可从式(17)求出其他几个未知数，再从式(10)～(11)或式(12)～(13)求出未知数ε2和ε4。显然，函数函数Fj(准3)是非常关键的，本文称之为设计方程组的特征函数。式(19)、式(20)、式(21)称为特征方程。
3特征方程的数值解
尽管可以推导出用准3来表示的Fj(准3)的具体表达式，但是形式非常复杂。本文使用数值算法来求解特征方程Fj(准3)=0。算例1：定向井A的靶点北坐标ΔN=300m、东坐标ΔE=400m、垂深ΔH=1000m；圆弧井段井眼曲率K2=4°/30m，K4=5°/30m；稳斜井段井斜角α3=42°；入靶井段井斜角α5=50°，方位角准5=30°；入靶井段长度ΔL5=100m。4.1特征函数图像和零点区间算法根据式(14)～(19)求得特征函数F1(准3)的离散点图(见图1)。从图1可以看出，特征函数F1(准3)是分段连续曲线，有两个奇点(在该点函数值为无穷大)，有两个零点并且函数值在零点附近是单调的。因此，首先要求出零点个数及所在的区间。为此，给出下面的区间搜索算法。算法1(区间搜索算法)：给定准3的搜索步长Δ准3(例如0.1°)，零点检测值Δf(例如0.01)。令count=360/Δ准3，M=0。①对xi=i×Δ准3(i=0，1，…，count)，计算fi=F1(xi)。②对0≤i＜count，检验：fi×fi+1≤0，fi-fi+1≤Δf如果上式成立，则令M=M+1，ai=xi，bi=xi+1。③输出M、[ai，bi]，i=0，1，…，M。这里，M为零点的个数，[ai，bi]为包含零点的区间。对于算例1，求得零点个数为2，零点区间为[60.8，60.9]和[261.1，261.2]。4.2求零点的二分法当连续函数在某个区间的端点函数异号时，可以使用二分法[6]求解这个非线性方程。根据算法1，在包含零点的区间上，特征函数值在区间端点是异号的，故可以使用二分法来求解特征方程。利用二分法求得算例1的特征函数的两个零点分别为60.86444°和261.14755°。图2是使用二分法在零点区间上求解特征方程的迭代过程。从图2可以看出，只需要10次左右的迭代就可以求出精度很好的数值解。
4井眼轨道设计过程
下面以未知数组合情况8为例叙述井眼轨道设计流程。算法2(拟解析法)：①计算特征函数F1(准3)的离散点数组；②使用算法1搜索特征函数零点区间；③在每个零点区间上求解特征方程，得到特征方程的零点数组{ηii=1，2，…，M}；④对每个准3=ηi，依序逐个计算其他未知数；⑤如果x1≥0且x3≥0，则式(24)～(31)给出的是设计方程组的真解，否则是伪解(应舍弃)。从算法2的计算过程来看，除了需要使用数值算法(二分法)来求解特征方程求未知数准3之外，其他所有的未知数都是使用解析公式[式(24)～(31)]直接计算的。这一特点与基于特征多项式实数根的拟解析法[1，2]的特点相似，故称算法2为拟解析法。由于解2的x3为负数，故解2为伪解，应舍弃。经验证，解1为设计方程组的真解。
5多解情况的处理能力
从算法2可知，当特征函数有多个零点时，设计方程组可能存在多个解。下面举例说明本文算法对多解情况的处理能力。算例2：定向井B的靶点北坐标ΔN=300m，东坐标ΔE=400m，垂深ΔH=1000m；圆弧井段井眼曲率K2=4°/30m，K4=5°/30m；稳斜井段井斜角α3=42°；入靶井段井斜角α5=50°，方位角准5=30°；造斜点ΔL1=300m。根据式(14)～(20)求得特征函数F2(准3)的离散点图(见图3)。从图3可以看出，特征函数F2(准3)是分段连续曲线，有两个奇点(在该点函数值为无穷大)，有两个零点并且函数值在零点附近是单调的。图3(b)是图3(a)中零点附近的放大图形。使用算法1求得零点个数为2，零点区间为[61.0，61.1]和[83.3，83.4]。利用二分法求得特征函数的两个零点分别为61.01694°和83.30046°。可使用算法2求出两个解。解1：准3=61.01694°，x3=0.29797，x5=0.09147；解2：准3=83.30046°，x3=0.07292，x5=0.25895。经验证，解1和解2均为设计方程组的真解。最终的设计结果见表4和表5。水平投影图和垂直投影图(投影方位角为30°)见图4和图5。在大量的测试实例中发现，解最多为2个，没有发现解多于2个的情况。
6结论
①对于限定稳斜角的三维S型井眼轨道设计问题，当给定入靶方向时，存在一个以其中的一个未知数为单变量的特征函数。如果该特征函数有零点，则其他未知数可用一组解析计算公式来计算。对10种不同的未知数组合情况，给出了特征函数的构造算法，推导出了解析计算公式的具体表达式。②大量的实际算例表明，本文算法可以快速求解限定稳斜角的三维S型轨道设计问题，并且在设计问题存在多个解的情况下，可以有效地求出全部解；如果设计问题有解，则最多有2个解。③论文发表本文推导设计问题的特征函数的数学思路也可应用于其他类型的已知设计条件的设计问题。
作者：林丽娜 单位：中国石油辽河油田分公司钻采工艺研究院
