基于随钻方位伽马测井的地层倾角自动识别

2018-12-26段友祥闫亚男孙歧峰李洪强

测井技术 2018年5期

段友祥,闫亚男,孙歧峰,李洪强,2

(1.中国石油大学(华东)计算机与通信工程学院,山东青岛 266580;2.中国石化集团胜利钻井工艺研究院,山东东营 257000)

0 引言

随钻方位伽马测井以常规自然伽马测井为基础,利用仪器不同方向上的多个伽马探测器进行测量,使伽马数据具有方位特性,可以根据方位特性反演得到地层构造情况,从而获得储层的边界、储层厚度和地层倾角等相关信息。在地质导向中可根据随钻方位伽马数据指导钻头准确入层,及时获得地层倾角,帮助分析地层结构,使井眼轨迹更准确地保持在目的层中,提高油气采收率。

地层倾角计算的2个必要参数分别为同一地层上不同2点之间的高度差和水平距离。目前,计算地层倾角主要有2种方案。①利用邻井信息,根据不同井钻遇同一地层时,2个钻点参数计算地层倾角。杨锦舟[1]应用邻井自然伽马建立先验模型结合人工经验进行了地质导向钻井;王谦等[2]根据邻井数据建立模型完成地层倾角计算和储层预测。但是,应用邻井信息计算倾角易受资源限制,邻井资料完整性影响着地层信息的准确性。②利用随钻方位伽马数据,根据方位数据计算地层倾角。侯爽[3]通过判断数据增量识别地层,但是不同的增量阈值会造成较大的效果差异;袁超[4]根据自然伽马图像进行了地层倾角及厚度计算,但需人工拾取计算的相关参数;McKinny K[5],Jiaxin Wang等[6]根据稳定区域均值识别边界并计算地层倾角,但是在实际钻井中数据很难稳定不变;Uzoh等[7]模拟不同地层条件下的方位伽马图像和曲线模型为人工识别服务;唐海全、肖红兵等[8]利用方位伽马曲线走势判断钻遇地层,但是无法计算地层倾角;Tyurin E等[9]在三维空间上计算地层真倾角,但仍需人工拾取地层特征。利用随钻方位伽马数据是计算地层倾角较为准确快速的方式,但是如何提取地层特征是目前难点之一。识别地层可以采用人工拾取、人机交互、计算机自动识别等方式。人工方式受经验影响较大,且速度慢;效率低;人机交互方式是人工拾取地层关键点,再由计算机计算,这种方式效率较高,但仍难以摆脱经验影响;计算机自动识别可以避免人为造成的影响和误差,而且效率更高,有利于快速作出正确的综合决策,提高开发效益,是目前的主要研究方向。

针对上述问题,本文以随钻方位伽马数据为基础,结合计算机图像处理等相关技术,对随钻方位伽马图像进行图像分割;通过地层边界曲线获取地层信息并完成地层倾角计算,以达到指导地质导向钻井的目的。该方法不仅实现了随钻方位伽马图像中的地层自动识别与倾角计算,而且满足了地质导向中的实际应用需求。

1 算法基本思想及流程

本文主要研究地层倾角自动识别,对随钻方位伽马数据进行归一化,并对归一化结果构造色阶直方图,综合考虑类内聚性和类间离散度,引入综合指标选取分割阈值;根据所选阈值对所有数据点进行阈值标记,自适应伽马数据测距,确定跟踪范围,按自定义跟踪顺序进行轮廓跟踪,获取地层边界;将哈夫变换正弦拟合与最小二乘拟合相结合,确定地层正弦公式;根据地层正弦参数计算地层倾角。算法流程见图1。

图1 算法流程图

2 伽马图像地层提取

2.1 综合指标阈值图像分割

识别地层就是要对伽马图像进行图像分割,即将伽马数据像素点分为不同的类,既要提高类内的相似度又要降低类间的相似度。本文首先对给定的随钻方位伽马数据进行归一化,再根据成像结果进行阈值分割。传统阈值分割采用图像灰度构造直方图,本文将伽马值映射为0～127级色阶,再转换为图像,但是在伽马数据→色阶→RGB值→灰度值的转化过程中,会丢失灰度值与伽马数据的正相关关系,所以本文将伽马数据映射后的色阶值作为数据输入,构造色阶直方图。

本文引用汪启伟[13]提出的综合指标,综合考虑衡量类内内聚性的熵和衡量类间离散度的类间差,选取阈值。熵值由式(1)计算,类间差由式(2)计算,综合指标由式(3)计算得到。

Hent(T)=

(1)

Hotsu(T)=pb(T)po(T)[mb(T)-mo(T)]2

(2)

(3)

式中,T为当前阈值;pb(T)、po(T)为2个类别所占总样本的比例;mb(T)、mo(T)为2个类别的类均值;went和wotsu分别为熵和类间差在计算综合指标时的权重,计算所有可选阈值下的综合指标,选取综合指标最大的值作为图像分割阈值。

2.2 自适应扩充邻域轮廓跟踪

图像分割阈值点即为表现伽马射线活动对比度平面上的点,确定阈值后,进行轮廓跟踪。轮廓跟踪通过顺序找出边缘点获取地层边界。具体过程如下

(1)列优先遍历图像中的所有数据点,对数据点(x,y),如果f(x,y)=value(value为阈值),则标记该点为阈值点;否则如果满足f(x,y)与f(x+1,y)位于阈值点两侧,则选取(x,y)与(x+1,y)中与阈值差距小的点标记为阈值点。

(2)建立二维链表BoundaryList作为边界链表。

(3)按照从上到下,从左到右的顺序找到一个未标记边的阈值点A,建立一个一维链表BList。

(4)若A不是噪声点或边界端点,在A点周围一定有一个阈值点,一般来讲应该按照左、左上、上、右上、右、右下、下、左下的顺序跟踪下一个阈值点。但是以阈值点作为边界有可能出现如图2所示情况,图2中红色点为阈值点,P1与P2同属一边界,但按照上述跟踪顺序,无法跟踪到P2。当出现P2、P3、P4的情况,从P2开始,按照上述顺序,会先跟踪到P4,导致丢失P3,但是针对实际情况,期望跟踪点依次为P2、P3、P4。根据随钻方位伽马数据的特性,只要边界连续,那么在图3中灰色区域内一定有阈值点。2点在纵向上相差行数由测距决定,可设定最小值,当测距较小时,2点纵向相差行数可能较大,根据实际情况,选取L内含有的数据行数为检测范围。

图2 轮廓跟踪图示图3 伽马数据轮廓区域示例

本文分析了伽马数据中阈值点的分布属性,创建了适用于方位伽马数据的自适应轮廓跟踪算法,根据测距自动收缩或扩张跟踪范围。以先右后左、顺时针、近距离优先为原则建立了如图4所示跟踪顺序。

图4 自定义轮廓跟踪顺序

参数N由式(4)计算

N=4(L/step)+12

(4)

式中,Step为伽马数据测距;L为单边检测范围,本文中取0.2。按上述顺序跟踪下一个边界点,标记该点所属边,将该点加入到边界链表BList。

(5)以最新点为A点,重复步骤(4),直至跟踪不到阈值点,将单边链表BList加入到总边界链表BoundaryList。

(6)跟踪一条边界结束,查找下一条边界,重复步骤(3)到步骤(5)。

3 地层倾角计算

3.1 改进正弦拟合

经过轮廓跟踪获取到多条地层边界后,要对地层边界进行曲线拟合,确定每个地层平面在三维空间中的位置。钻孔与表现伽马射线活动对比度的平面相交,其形态见图5。在三维空间内投影在地层平面上显示为椭圆,在井壁的平面展开图上,观察到特征正弦响应,所以要将地层展开面上的地层边界线拟合为正弦曲线。

图5 地层与井筒相交示意图

井壁展开面上的正弦曲线表达式见式(5)

f(x)=Asin (ωx+φ)+y0

(5)

式中,A为振幅;ω为角速度;φ为初相位;y0为偏距,振幅是倾角计算中的关键参数。如图5所示,井壁展开面上是一个完整的正弦周期,即周期T是井壁展开面的宽度,角速度可由式(6)计算

(6)

在正弦公式中,只需要拟合3个参数:振幅A,初相位φ,偏距y0。

最小二乘拟合、迭代拟合[10]、基于哈夫变换的拟合[11-12]是测井数据处理中常用的3种正弦拟合算法,它们各有自己的优势和不足。最小二乘算法计算精度高,但是当范围过大时,计算量过大,计算速度慢;迭代拟合算法与最小二乘算法相比速度大幅提高,在单参数拟合时收敛效果较好,但在三参数拟合时难以保证收敛于最优值;哈夫变换算法速度最快,适合大量点中选择少数点拟合,是一个拾取和拟合相结合的过程,但是不适用于本文应用场景。为此,本文在研究分析3种拟合算法优缺点的基础上,根据伽马数据本身的特点,提出了一种适用于伽马数据的综合拟合算法。该算法过程如下:

(1)边界点排序。对轮廓跟踪中得到的二维链表BoundaryList,遍历每1条一维链表BList,对链表中的数据以x为关键字排序。

(2)确定初始偏移量。遍历链表BList中的点,对一个点(xa,ya)在BList中查询满足xb=xa+T/2的点(xb,yb),偏移y0根据式(7)计算

(7)

式(7)的原理是正弦本身的对称性。如图6所示,正弦上2点S1、S2横坐标相差半周期,2点中点Q一定在直线y=y0上。对所有计算结果,建立投票机制选择y0。

图6 正弦基线计算图示

(3)确定最小二乘参数范围。遍历链表BList中每一个点,得到y的最大值ymax、最小值ymin,以及y取最大值时对应的横坐标xmax,正弦公式中偏移量的遍历范围是[y0-(ymax-ymin)/6,y0+(ymax-ymin)/6],振幅的遍历范围是[0,(ymax-ymin)×2/3];由于正弦曲线本身的性质,初相位的遍历范围由式(8)确定

(8)

(4)计算残差。对3个参数的可取值自由组合,对每种组合(yi,Aj,φk)计算残差,残差Errorijk由式(9)计算

(9)

(5)确定参数。比较所有残差,残差最小的三参数取值即为该条边界的参数值。

(6)计算下一条边界。取BoundaryList中的下一条边界,重复步骤(2)到步骤(5),直至所有边界计算完毕。

3.2 倾角计算

地层倾角是地质导向中指导钻头行进的关键参数,根据正弦公式即可计算地层倾角。如图5所示,平面边界与钻孔相交的相对倾角可以用式(10)确定

(10)

式中,D为钻孔直径,计算中可认为等于钻头尺寸;H为正弦曲线在井壁展开面上的展布高度,等于正弦曲线振幅的2倍;根据Jiaxin Wang等[6]论文中的定义,DI为伽马图像的成像深度。Jiaxin Wang详细研究了不同地层倾角、地层密度、钻孔尺寸等环境对DI的影响,本文采用该思想,利用式(10),根据上文所得平面正弦参数计算地层倾角。

4 实验与分析

在随钻测井中,采用泥浆脉冲遥测方式将随钻方位伽马测井数据实时传输至地面,并对传输数据进行实时预处理与伽马成像。在实际钻井中,地质情况复杂,地层构造信息不精确。基于人工构造的地层模型,对测井响应进行正演模拟,可达到实测数据效果,而且可以提供准确的地层界面信息作对比参考。本文分析了前人研究的各种正演算法,最终采用邵才瑞[14]、Zhen Qin等[15]所用积分思想进行了随钻方位伽马数据模拟。

以随钻方位伽马模拟数据为基础,实验环境采用Visual Studio 2013和Open Scene Graph,实现了随钻方位伽马图像的地层自动识别与倾角计算。采用滑动窗口分析伽马图像,以最新更新数据为滑动窗口结束位置,使用本文算法处理滑动窗口内的图像数据。滑动窗口的窗长和滑动步长可根据实际情况设定,窗长设定过长,会导致数据量过多,降低分析效率;窗长设定过短,可能造成大倾角地层的倾角信息不完整,测得倾角存在误差;滑动步长设定过长可能会导致发现地层较晚,错过最佳控制钻头的时机;滑动步长设定过短会导致频繁分析伽马图像,降低分析效率,严重可能导致延误地层发现时机。设定更新井轨迹上一个采样点八方位数据的时间为一个单位时间,伽马数据的测距为0.1 m,八方位伽马数据横向内插为64个点数据,表1为该条件下不同窗长和滑动步长下的单位时间处理数据点数的比较。

表1 不同滑动窗口处理数据量比较

正演模型中地层模型为砂岩和泥岩,地层模型相关参数见表2。K每秒平均发射光子数3.4/gs;Th每秒平均发射光子数1×104/gs;U每秒平均发射光子数2.8×104/gs[16-17]。

表2 正演模型相关地层参数

模拟4组数据作为实验数据,图7、图8、图9、图10分别为4组数据实验结果,表3为4组数据的识别结果与实际角度的对比及误差。

表3 本文算法地层倾角识别准确性

图7 第1组数据识别结果

图8 第2组数据识别结果

图9 第3组数据识别结果

图10 第4组数据识别结果

实验组别拟合点数角度/(°)拟合算法耗时/s残差第1组7284.29迭代最小二乘拟合0.04213.582哈夫变换0.0040.1011传统最小二乘拟合25.7940.0606本文算法14.9580.0606第2组6578.69迭代最小二乘拟合0.0112.9477哈夫变换0.0040.0238传统最小二乘拟合5.8290.0101本文算法3.2810.0101第3组6563.43迭代最小二乘拟合0.0030.5174哈夫变换0.0030.0116传统最小二乘拟合1.2260.0072本文算法0.6150.0072第4组6526.57迭代最小二乘拟合0.0010.3158哈夫变换0.0030.0167传统最小二乘拟合0.1160.0027本文算法0.0620.0027

本文所用正弦拟合算法,与迭代拟合和哈夫变换拟合相比残差更小;与传统最小二乘算法相比,耗时更少,表4为4种算法耗时及残差对比。

实验证明,本文所用图像分割、轮廓跟踪、正弦拟合算法取得了较好效果,可自动识别随钻方位伽马图像中的地层并计算地层倾角,且精确度较高。