最优化变尺度法在汶川科学钻探岩性识别中的应用

2014-05-10王智潘和平杜婷

测井技术 2014年3期

王智, 潘和平, 杜婷

(中国地质大学地球物理与空间信息学院, 湖北武汉 430074)

0 引言

测井岩性识别方法主要有交会图法、成像测井识别法[1]、神经网络法[2]、对应分析识别法[3]等。传统交会图法虽然在一定程度上可以区分部分岩性,但是在边界上岩性很难区分。神经网络法需要的学习样本要求尽量全面覆盖各种岩性的测井响应特征,才能保证在学习过程中学习样本的拟合误差在减小的同时待预测数据的预测误差也得到有效的控制。取心资料WFSD-2号钻孔岩心总长1 641.26 m,500 m以下取心率达到91.6%,因此将地质编录资料视为取心资料。根据岩心资料统计,WFSD-2井所钻遇地层岩性主要为岩浆岩与沉积岩。建立由不同岩石组成的体积模型不符合体积模型的概念,本质上应该建立矿物模型。由于矿物模型中矿物参数确定困难,研究表明,建立由岩石组成的体积模型识别是可行的[4-5],因此,本文采用这一方式。

本文结合汶川科钻2号井岩心资料和测井资料,依据地层组分分析原理[6],采用最优化方法[7]计算出WFSD-2号井的岩性成分含量,为该地区复杂岩性识别提供参考。

1 最优化变尺度法

牛顿法最突出的优点就是收敛速度快,但其缺点是要计算海森矩阵及逆矩阵,它的迭代公式为

Vk+1=Vk-Q-1F(Vk)

(1)

式中,V=(V1,V2,…,VM)T是用最优化方法待求的参数;k表示迭代次数;Vk+1、Vk分别为第k+1、k次迭代得到的V值;F(V)为目标函数的梯度;Q-1为海森矩阵。

为保持牛顿法收敛速度快的优点,克服它的缺点,要计算Q-1,变尺度法用HT代替Q-1,便有

(2)

为取得更大的灵活性,考虑更一般的迭代公式

(3)

式中,ak为最优步长;Hk为尺度矩阵,是随逐次迭代而变化的矩阵,即变尺度。这种迭代求解问题的方法称为变尺度法。只要在每次迭代中使目标函数值减小,并保持方向矩阵Hk+1的正定对称性,则Hk一定能逼近海森矩阵的逆矩阵。

最优化变尺度反演方法是把目标函数分解成不同尺度的分量,根据不同尺度的目标函数特征逐步搜索全局极小值点。一般在大尺度(低波数)上,目标函数极值点少且分得很开,因此很容易搜索出大尺度上的全局极小点,然后以该点为起始点,在其附近能够较容易地搜索到中等尺度上的全局极小点。如此,不断缩小尺度,再不断加入高分辨率的目标函数。最后,当目标函数的尺度降至原始尺度时,对应搜索到的全局极小点就是真正的全局极小点[8]。

变尺度方法种类较多,较为著名的有SR1、DEP和BFGS,其尺度矩阵分别为

(4)

(5)

(6)

2 建立的模型及响应方程

储集层的火山岩、碳酸盐岩、煤层、页岩或者碎屑岩都可以看成是3种基本的地层组分:孔隙流体(油、气、水)、泥质以及岩石的各种骨架矿物[4]。根据物理模型,可写出各种测井响应方程,例如声波测井响应方程为

(7)

式中,Δtmf、Δtsh、Δtmai分别为混合流体、泥质与骨架的声波时差;xmf为孔隙度;xsh为泥质含量;xmai为骨架矿物体积含量。

同理可写出其他测井响应方程,用通式表示为

(8)

式中,m为测井仪器的个数;n为组成地层组分的个数;xi为第i种组分的相对含量;Aij为第i种组分对第j种仪器的响应值;Bj为地层对第j种仪器的响应值。

式(8)中m个方程所组成的方程组,从数学角度,当mn;当m>n时,方程组为超定线性方程组,具有一个最优解,但可能出现x<0或x>1的现象,这种结果不符合实际地层的情况[5],所以需要加入约束条件

(9)

式中,xi为第i种组分的相对含量。

由线性最小二乘原理求解这一约束线性方程组的问题可转换成求极值问题

0≤xi≤1

(i=1,2,…,n;j=1,2,…,m)

(10)

式中,不同测井值的量纲不同,它们的测量值也存在很大差别。在实际计算中需要将式(10)中目标函数的系数A和B进行标准化处理,使各种仪器的系数A和B都成为无量纲的数,并在同一数量级上,这样可使得各种测井方法对最终结果都具有相同的贡献。本文将所有的各个组分之和等于1作为约束条件作了处理

(11)

式中,ε取0.05。以此作为约束条件时,即各个组分之总和会有5%左右的误差范围,因此,会出现总和不为1的情况。

2.1 测井目标函数

设进行岩性分析时所选用的测井方法有m种,其实测值为Li,i=1,2,…,m;所研究的岩石由n种成分组成,其体积含量Xj,j=1,2,…,n;则根据最小二乘法原理可得测井目标函数

(12)

2.2 收敛准则

收敛准则的选取要注意几个重要的因素,即目标函数达极小值的必要条件,根据目标函数的值下降以及岩性分析的约束条件和迭代次数的极限条件等因素提出收敛淮则[9]

l=‖F(x)‖≤ε1

δ=abs[F(Xj+1)-F(Xj)/F(Xj+1)]≤ε2

(13)

式中,ε1、ε2为预先给定的充分小的正数,称为终止限;l为范数;j=1,2,…,n为岩性成分的体积含量。

3 地层组分测井响应参数分析

实际资料处理中,依据岩性复杂程度选择合适的测井曲线,并遵循6个原则:①选择能够反映孔隙度变化的测井曲线;②选择能够反映岩石矿物成分的测井曲线;③选择能够反映流体性质的测井曲线;④所选测井曲线数须大于地层组分数,即保证方程组为超定线性方程组,可获得一个最优解;⑤对于最常见的砂泥岩剖面,地层组分主要包括油、水、泥岩、砂岩,可选择声波时差测井曲线、密度测井曲线、中子测井曲线、自然伽马测井曲线、深电阻率测井曲线等联立方程组进行求解;⑥对于岩性比较复杂的地层除选择⑤中提到的5条常规测井曲线外,还应选择敏感的造岩元素测井响应(ECS测井资料)或者伽马能谱测井[6]。针对以上原则,本文选取了AC、CNL、DEN、GR、U、Th、K等7条常规测井曲线。

由式(7)可知,方程系数Aij的确定最关键,流体组分的测井响应参数通过理论计算即可准确确定。岩石骨架组分的测井响应分为2部分:①岩石骨架的自然伽马、ECS测井、伽马能谱等测井响应参数的确定主要是通过录井或者岩心资料,结合电成像测井资料,挑选出典型稳定的岩石类型层段,统计出它们的平均测井响应值作为它们的测井响应参数;②确定岩石骨架的密度、中子、声波时差的测井响应参数的方法主要有4种。

(1) 应用交会图技术确定研究区各类岩石的骨架参数,即针对所划分出的岩性岩相类型,通过二元线性回归分别建立岩心分析孔隙度与中子孔隙度、密度孔隙度、声波孔隙度的关系式,则回归方程的斜率为岩石骨架与流体测井响应的差值,其截距即为相应的岩石骨架参数[5,10];

(2) 直接使用理论岩浆岩(火成岩)骨架参数值[11];

(3) 利用ECS测井资料计算骨架参数[12-13];

(4) 作某种岩浆岩岩石的中子(声波、密度)与电阻率交会图。由于岩浆岩骨架不导电,所以,当电阻率趋于无穷大,对应的测井值就为骨架值[12-13]。

由于无岩心数据和ECS测井,所以方法(1)与(3)不能使用,而岩浆岩的岩性非常复杂,并且随着地区的不同变化差异也较大,所以最终采取了方法(4)。由取心资料可知,本文研究WFSD-2井段所钻遇地层岩浆岩岩性主要有花岗岩、闪长岩、花岗闪长岩以及凝灰岩。图1至图8分别是各种岩性的声波、中子、密度测井曲线与电阻率测井曲线的交会图。根据电阻率趋于无穷大时对应的测井值基本可以确定岩石骨架参数。

图1 花岗岩声波时差—电阻率交会图*非法定计量单位,1 ft=12 in=0.304 8 m,下同

图2 花岗岩中子值—电阻率交会图

图3 花岗岩密度值—电阻率交会图

图4 闪长岩声波时差—电阻率交会图

图5 闪长岩中子值—电阻率交会图

图6 花岗闪长岩声波时差—电阻率交会图

图7 凝灰岩声波时差—电阻率交会图

图8 凝灰岩中子值—电阻率交会图

部分岩性的某些骨架参数无法用这种方法确定,如闪长岩、花岗闪长岩与凝灰岩骨架密度值以及花岗闪长岩的中子骨架参数值,其原因可能是由于井眼条件差,造成密度测井不能反映真实的地层密度;或是由于电阻率受岩性、井眼、泥浆侵入等影响,视电阻率不能反映真实的电阻率而导致电阻率测井值与孔隙度的相关性太差。此时可以充分利用部分已知骨架参数确定同一岩性的未知骨架参数[13],例如声波、中子和密度测井值都与孔隙度呈线性关系,若已知岩性的声波时差与中子值,根据线性关系得到中子和声波时差值为已知骨架值时对应的密度值即为相应岩性的密度值。选取岩性稳定、井眼条件好的井段分别作其交会图。

图9为闪长岩声波时差与密度交会图。由前所述方法(见图4)可知,闪长岩的声波时差为52 μs/ft,计算得闪长岩的密度为2.68 g/cm3。

图9 闪长岩声波时差与密度值交会图

图10、图11分别为花岗闪长岩声波时差与中子、密度交会图。由前所述方法(见图6)可知,花岗闪长岩的声波时差为51 μs/ft,计算得花岗闪长岩的中子值为13.67 p.u.,计算得花岗闪长岩的密度值为2.73 g/cm3。

图10 花岗闪长岩声波时差与中子值交会图

图11 花岗闪长岩声波时差与密度值交会图

图12为凝灰岩声波时差与密度值交会图。由前所述方法(见图7)可知,凝灰岩的声波时差为55 μs/ft,计算得凝灰岩的密度值为2.87 g/cm3。

图12 凝灰岩声波时差与密度值交会图

WFSD-2井岩浆岩各种岩性骨架的测井响应值见表1。

表1 WFSD-2井各种岩性骨架的测井响应值

4 应用效果及实例

为验证上述方法对WFSD-2井岩浆岩的岩性识别能力,结合地质取心资料分别选取WFSD-2井各个典型稳定段的花岗岩、闪长岩、花岗闪长岩、凝灰岩层段作为样本,计算其岩性成分含量(见表2至表5)。

表2 花岗岩样本识别结果

表3 闪长岩样本识别结果

表4 花岗闪长岩样本识别结果

表5 凝灰岩样本识别结果

花岗岩99个样本,识别率为100%;闪长岩145个样本,与岩心编录相符合的为110个样本,识别率约为78%;花岗闪长岩158个样本,与岩心编录相符合的为110个样本,识别率约为70%;凝灰岩283个样本,与岩心编录相符合的为233个样本,识别率约为82%。对比可见,该方法对花岗岩的识别度是最高的,其次是凝灰岩,闪长岩与花岗闪长岩识别度误差较大,闪长岩多被识别为花岗闪长岩和凝灰岩。从测井曲线看,由于闪长岩被绿泥石化或含有杂质或破碎,使得闪长岩的测井响应特征与花岗闪长岩和凝灰岩相似,导致部分深度处识别误差较大。花岗闪长岩多被识别为凝灰岩,在测井曲线响应特征上较为相似,容易被误判。另外,这些岩性在岩心编录和测井人工识别的过程当中本来就比较难以区分,如闪长岩与花岗闪长岩的矿物含量百分比问题,会影响到岩性判别。因此,本文方法产生的这些误差部分原因是由于岩心编录和测井人工识别引起,在程序识别当中是可以被允许的。

总体上,该方法对WFSD-2井岩浆岩的岩性识别是可行的,岩性成分含量的计算结果较可靠,证实本文方法的可行性。