设计约束方程组的求解方法及计算机编程算法的技术细节

1前言
本文着重研究了设计约束方程组的求解方法，详细讨论了可以计算机编程实现的数值算法的技术细节。约定：除特别声明外，文中具有长度量纲的参数其单位为m，角度的单位为rad，井眼曲率和角度变化率的单位为rad/m。
2设计基础
设计井眼轨道由连续光滑的多段空间曲线构成，每个分段空间曲线称为设计井段。悬链线轨道由4个设计井段构成———直井段、圆弧过渡段、悬链线段、稳斜段。方位漂移轨道设计的方法是先在垂直剖面图上进行设计，规定每个设计井段的井斜角变化规律，然后再结合方位变化率进行空间轨道设计。
2.1井斜角函数
在每个设计井段上，井斜角随井深变化的规律是相同的，可以用下面的井斜角函数来表示：α=f(△L)(1)△L=L-Lb式中：α为设计井段上任意点处的井斜角；△L为井深增量；L为设计井段上任意点处的井深；Lb为开始点处的井深。直井段和稳斜段：α=αb(2)式中：αb为设计井段开始点处的井斜角圆弧过渡段：α=αb+kα△L(3)式中：kα为造斜率。悬链线段[3]：cotα=cotαb-△La(4)式中：a为悬链线特征参数，m。
2.2井斜单元
在设计轨道上任意取一连续曲线段，称之为井段。本文用设计井段和井段来区分设计轨道上具有不同属性的连续曲线段。设计井段上的任一井段称为一个井斜单元[2]。井斜单元小于或者等于、但不大于某个设计井段；一个设计井段可能划分为多个井斜单元。下面的情况绝不会出现：某个井斜单元的上端点在一个设计井段内部，而下端点在另一个设计井段内部。
2.3方位单元
将从井口到靶点的垂深划分为连续的m个区间，每个区间称为垂深区间，在每个垂深区间上给定方位变化率的数值。按照垂深与井深的对应关系，每个垂深区间对应于设计轨道上的一个井段，该井段上的方位角变化规律如下：准=准b+k准△L(5)式中：准b为井段开始点处的方位角；准为井段上任意点处的方位角；k准为方位变化率。具有上述性质的井段称为设计轨道上的一个方位单元，简称方位单元[2]。
2.4计算单元
对于方位单元可能会出现下面的情况：上端点在一个设计井段内部，而下端点在另一个设计井段内部。为了避免出现这种情况，可以将方位单元再划分成多个更小的井段，使得每个井段都是一个井斜单元。这种既是井斜单元又是方位单元的井段称之为计算单元，文献[2]中也称为细分单元等。但笔者认为“细分”这个词的语义不甚明晰，由于在漂移轨道设计时是以这样的井段为最小划分单元的，因此笔者认为称之为计算单元更言简意赅。
3设计约束
方程组假设整个设计轨道可以划分为n个计算单元，则可以给出下面的设计约束方程组[2]：ni=1Σ△Ni-Nt=0(6)ni=1Σ△Ei-Et=0(7)ni=1Σ△Hi-Ht=0(8)式中：△Ni、△Ei、△Hi分别为第i个计算单元的北坐标增量、东坐标增量、垂深增量；Nt、Et、Ht分别为靶点的北坐标、东坐标、垂深。已知参数包括：造斜点深度、圆弧过渡段造斜率、悬链线段初始井斜角、稳斜段井斜角、垂深区间及方位变化率。求解设计参数：定向方位角、悬链线特征参数、稳斜段的段长。
3.1坐标增量计算
坐标增量计算公式如下：△Ni=LiLi-1乙sinαcos准dL(9)△Ei=LiLi-1乙sinαsin准dL(10)△Hi=LiLi-1乙cosαdL(11)△Si=LiLi-1乙sinαdL(12)式中：Li-1、Li分别为第i个计算单元的开始井深和结束井深；△Si为水平投影长度增量。式(9)、式(10)中的定积分无法写成封闭的形式，在计算时只能使用数值积分法[4]。
3.2独立未知数与隐含未知数
方程组(6)～(8)只有3个方程组，理论上可以求出3个未知数，要求这些未知数之间是独立的。另外可以看到，在完成设计约束方程组求解之前，某些井深单元的端点井深或井斜角等参数是未知的，例如悬链线段上的井深单元。这些未知数也需要在求解过程中确定出来，不过它们都可以根据已知设计参数或者独立未知数计算出来，称之为隐含未知数。
4垂深增量公式
在坐标增量公式中，垂深增量公式(11)具有比较特殊的意义，这不仅在于积分函数只与井斜角有关、从而有可能求出积分的原函数，而且方位单元是根据已知垂深来确定的，利用垂深已知性可以确定出井深单元的其他参数来。
4.1显式公式将式(2)～(4)代入式(11)，得：直井段和稳斜段：△Hi=△Licosαi(13)圆弧过渡段：△Hi=sinαi-sinαi-1kαi(14)悬链线段：△Hi=a1sinαi-1-1sinαiii(15)
4.2推论如果知道一个井深单元的上端点的井深和井斜角，则可以计算出下端点的井深和井斜角以及段长。例如，假设造斜率或者曲线特征参数为已知数，则从式(14)、式(15)可以先求出井斜角αi，再代入式(3)、式(4)求出段长△Li和井深Li。而对于稳斜单元，可以直接从式(13)求出这些参数。
4.3用于降维处理如果能够将要求解的某个独立未知数用其他的独立未知数来表示，则可以将求解三元非线性方程组问题简化成求解二元非线性方程组问题。未知数的减少可以降低方程组的求解复杂度，从而提高计算速度。对整个设计轨道列出垂深方程，如下：Hz+sinαbkα+a1sinαb-1sinαtii+△lcosαt=Ht(16)式中：Hz为造斜点垂深；αt、△l分别为稳斜段的井斜角和段长。从式(16)解得：△l=Ht-Hz-sinαbkα+a1sinαb-1sinαtiicosαt(17)可见，稳斜段长△l可以从式(7)直接计算出来，求解设计约束方程组时只需要使用前两个方程代式(6)、式(7)即可。如果用均匀网格法来求解方程组，初始网格点个数为：200×2000×2000=8×108当稳斜段长不作为独立未知数时，初始网格点个数减少为：200×2000=4×105即求解的空间规模降低了3个数量级，求解难度降低。
5求解方程
组记方程(6)～(8)的左端分别为F1、F2和F3，无论使用哪种数值算法求解该方程组，都需要反复计算这3个值。
5.1隐含未知数的递推计算
用{hrr=0，1，…，m}表示方位单元的端点垂深序列(已知设计参数)，h0=0，hm=Hz。第r个方位单元上的方位变化率kr为已知设计参数。用Ω={Hii=0，1，…，n}表示最终得到的计算单元的端点垂深序列，H0=0，Hn=Hz。第i个计算单元上的方位变化率Ki待确定。对于设计轨道上的第j个设计井段，假设已知其上端点处的井深L0(j)、垂深H0(j)、井斜角α0(j)，则其他已知参数会出现三种情况：①已知设计井段的垂深增量△H(j)，例如直井段。在方位单元垂深序列中查找两个下标号r1和r2，使得下式成立：hr1-1≤H0(j) 5.2设计约束方程组的数值求解
经过降维处理之后的设计约束方程组(6)、(7)为二元非线性方程组，没有解析解，需要使用数值算法求数值解(近似解)。求解非线性方程组的数值算法有很多种[5，6]，大部分算法需要使用导数信息，并且迭代初始值对算法的收敛性有较大影响，如果迭代初始值选择不当，则迭代过程可能不收敛或者收敛速度很慢，在方程组有多个解的情况下，还有可能收敛到伪解。算法研制的最终目的是为钻井设计人员(用户)提供一套可靠性好的计算机软件，用户在使用该软件时，只需要给定必要的设计参数，不必设置太多的算法控制参数就能够快速求出轨道设计问题的解来。本着这一原则，下面给出一个具体的迭代算法———缩半网格法。记x、y、z为方程组(6)～(8)的3个独立未知数：定向方位角、悬链线特征参数、稳斜段的段长。方程左端分别记为F1(x、y、z)、F2(x、y、z)和F3(x、y、z)，前面已经说明从垂深增量方程可以将某个参数z表示为其他2个参数的函数[见式(17)]：z=λ(x，y)(20)再记：F(x，y)=F1(x，y，λ(x，y))2+F2(x，y，λ(x，y))2(21)假设未知数取值范围为：xmin(0)≤x≤xmax(0)(22)ymin(0)≤y≤ymax(0)(23)用步长为δx(0)和δy(0)将约束矩形划分为Mx(0)×My(0)的网格，其中：δx(0)=xmax(0)-xmin(0)Mx(0)(24)δy(0)=ymax(0)-ymin(0)My(0)(25)令：xi=xmin(0)+i×δx(0)，i=0，1，…，Mx(0)(26)yi=ymin(0)+j×δy(0)，j=0，1，…，My(0)(27)对每个网格点(xi，yj)，利用垂深方程求出对应的第3个待定参数zij=λ(xi，yj)，然后使用5.1节中的方法求出全部的隐含未知数并得到方程左端项的值，再代入式(21)求出网格点函数值Fij。求出所有的网格点函数值中最小的函数值，对应的网格点为(x(0)，y(0))，以该网格点为矩形中心将初始矩形缩小一半，得到新的约束矩形：xmin(1)=max{xmin(0)，x(0)-wx(0)}(28)xmax(1)=min{xmax(0)，x(0)-wx(0)}(29)ymin(1)=max{ymin(0)，y(0)-wy(0)}(30)ymax(1)=min{ymax(0)，y(0)-wy(0)}(31)其中：wx(0)=(xmax(0)-xmin(0))/2(32)wy(0)=(ymax(0)-ymin(0))/2(33)用新的约束矩形重复上述计算过程，直到约束矩形的边长或者最小函数值小于给定的允许误差时停止迭代过程。
6结论
①通过垂深方程可以将3个独立未知数中的一个表示为其他未知数的函数，从而使得设计约束方程组可以降维为二元非线性方程组，降低了数值求解的规模和难度。②结合隐含未知数的递推计算策略，缩半网格法能够可靠地求出降维后的设计约束方程组的数值解，特别适用于计算机编程实现。