基于改进粒子群算法的阵列侧向测井快速反演

2021-06-03倪小威徐思慧

能源与环保 2021年5期

彭杰，倪小威，徐思慧，汤鹏

(1.长江大学油气资源与勘探技术教育部重点实验室，湖北武汉 430100； 2.长江大学地球物理与石油资源学院，湖北武汉 430100； 3.塔里木油田分公司油气田产能建设事业部，新疆库尔勒 841000； 4.塔里木油田分公司勘探事业部，新疆库尔勒 841000； 5.中国石油塔里木油田分公司安全环保与工程监督中心，新疆库尔勒 841000)

地球物理反演是根据观测数据、利用数学方法来估算地球物理参数的分布状态的一种有效手段。在电法测井反演中，通过对测井资料进行快速有效的一维反演，可以初步反映地层原始特征，同时也可以为测井质量好坏提供重要的参照[1]。但是这种实时反演技术对算法要求较高，要有较高的收敛速度。由于在反演过程中每一次迭代都会增大数值的计算量，这样也会造成反演效率较低[2]。对一维反演算法的优化处理可以从2个方向考虑：①重新寻找新的算法；②对传统的一维算法进行改进，使其具有符合实际生产需求的精确度和速度。对测井数据反演地层电阻率的非线性问题，国内外学者进行了许多探索。

本文以阵列侧向测井数据为研究依据，一般的阵列侧向反演算法有牛顿法、最小二乘法等[3-4]，这类算法具有较快的收敛速度，但是反演结果受初始值影响大，如何进行合理初值设定又是一个问题。随着反演算法的不断发展，部分学者将模拟退火、遗传算法及差分进化算法引入了电测井资料的反演处理中[5-7]。这类算法具备较强的全局搜索能力，同时反演结果几乎不受初值的影响，但此类算法往往收敛速度较慢，且编程复杂，不利用现场的推广使用。随着仿生学科的发展，1987年Kennedy和Eberhart基于对鸟群觅食规律的研究提出了粒子群优化算法，粒子群算法由于结构简单、易于编程实现、收敛速度快、搜索范围大等优点，已在不同领域得到了广泛应用[8]，同时也得到了油田现场工作者的青睐。并且在满足反演精度的要求下，能够替代反演过程中某些正演计算方法，目前主要以图版数字化方法来代替严格正演计算[9]。

为了解决基本粒子群算法无法平衡局部搜索(决定收敛速度)与全局搜索能力(决定寻优成功率及精度)的问题，提出了自适应非线性惯性权重，使得算法在迭代前期具有较强的全局搜索能力，利于寻找全局最优解，而在后期惯性权重较大，加强算法的收敛速度，同时非线性变化的权重也加强了种群粒子间的信息交流，更利于全局最优解的寻找。同时为进一步提高反演效率，利用预先计算的阵列侧向测井伪几何因子数据表格代替严格正演计算，在反演过程中调用线性插值算法数据，也加快了反演的计算速度。采用此反演方法进行三参数反演，反演结果与严格动态泥浆侵入正演结果符合度较好。本文方法可为现场阵列侧向测井数据的实时反演处理提供一定技术支持。

1 阵列侧向电极系结构

仪器的测井电极系结构如图1所示：主电极A0、A1与A1′分别为对应的屏蔽电极，A2与A2′、A3与A3′、A4与A4′、A5与A5′也都为相对应的屏幕电极，共5对；M1与M1′、M2与M2′、M3与M3′、M4与M4′、M5与M5′、M6与M6′、M7与M7′都是相对应的监督电极，共有7对，它们都呈主电极中心对称分布。仪器测量了4种探测深度的电阻率曲线，它们由浅到深的依次为MLR1、MLR2、MLR3、MLR4。

图1 阵列侧向测井电极系结构Fig.1 Structure of array laterolog electrode system

2 有限元正演算法理论及伪几何因子

求取侧向测井正演响应的实质在于求取一光滑的电位函数μ对侧向仪器形成的稳流场进行描述[10]，电位函数μ应满足式(1)：

(1)

式中，R为模型中不同区域的电阻率。在井眼中，R为钻井液电阻率；在地层中，R为原状地层电阻率。

给式(1)添加相对应的边界条件，规定在仪器的绝缘区域内电位函数μ的法向导数为0；距仪器源无限远的区域，电位μ也为0。这样就变成了求正演响应函数的定解问题。

对定解问题的求解往往采取等效的思想将定解问题转化为泛函φ的极值问题[11]：

(2)

式中，IE为电极系的电流；μE为电极系的电位；E为电极个数；积分区间包含仪器表面和无穷远边界包围的空间，对所有电极进行求和。

利用有限元素法求取阵列侧向测井正演响应的详细过程见文献[12]。

通常引入伪几何因子理论来描述[13]侧向类测井仪器的径向探测特性。即视电阻率可以看作是冲洗带、井眼及原状地层电阻率的加权和，而不考虑地层纵向上的电阻率变化，即一维模型。冲洗带电阻率的权值就是伪几何因子λ，λ表征了阵列侧向测井正演响应受泥浆侵入的程度大小。经过井眼校正后，排除该因素影响，那么视电阻率就等于冲洗带及原始地层的综合贡献值，伪几何因子λ可用式(3)来描述。

(3)

式中，RMLRi是随着i变化的视电阻值；Rt为地层真实电阻率；Rxo为冲洗带电阻率。

研究λ在不同Rt/Rxo的条件下随泥浆侵入深度变化的特征，并制作MLR1—MLR4对应的伪几何因子数据表格。表1、表2分别为MLR1、MLR4的伪几何因子数据，在实际反演过程中，利用线性表格插值代替严格正演计算。表中Di为泥浆侵入深度。

表1 MLR1伪几何因子数据Tab.1 Pseudo-geometric factor data of MLR1

表2 MLR4伪几何因子数据Tab.1 Pseudo-geometric factor data of MLR4

3 反演原理及算法

3.1 非线性反演原理

电阻率反演实质就是通过实际测井资料来推断模型电阻率大小。通过比较M个未知参数与N个实际数据，当它们误差在所设定的阈值以内，则可以把所求的M个参数看作实际地层参数[14]。研究基于阵列侧向测井仪器的三参数反演，地层是水平均匀层状的，模型表示如下：

(4)

式中，f(x)为所求目标函数；xi为所求模型参数；yj为对应j种方法的实际测井数据；Φj为对应第j种方法的正演响应算子，本文中，Φj满足：

Φj=λj·Rxo+(1-λj)·Rt

(5)

式中，λj为正演数据线性插值结果。其中，j=1，2，3，4时分别对应曲线为MLR1，MLR2，MLR3，MLR4。

3.2 基本粒子群算法

(6)

(7)

式中，c1、c2为加速因子，通常在[0 2]之间取值；r1、r2为[0 1]之间的随机数。

将式(7)代入式(4)中，当目标函数小于预先设定的阈值或者迭代达到最大次数，输出目标函数最小的情况下对应的位置属性作为反演结果。

3.3 改进粒子群算法

(1)采用混沌序列产生初始种群。混沌信号在全空间内具有随机性和遍历性，这一特性能在求解最优化问题时，有效地避免算法陷入局部极值[16]。常用的产生混沌序列的方法是一维logistic映射，即：

rk+1=μrk(1-rk)

(8)

式中，rk为随机数；k为迭代次数；μ为控制参数，当μ=4时，能够实现[0 1]的完全映射，从而达到完全的混沌状态，提高算法的全局寻优能力。

将式(8)引入到初始种群产生中，即：

(9)

(2)自适应非线性惯性权重。基本粒子群算法无法平衡算法的前期寻优及后期收敛能力，前人[17]为解决这个问题，在速度位置的更新公式中引入了惯性权重ω，分别提出了固值惯性权重、线性惯性权重、非线性惯性权重等方法[18-19]。通常对于全局搜索比较好的优化方法是在算法搜索前期就增强它的广度探索能力以便得到优秀的种子解，后期搜索算法应具有较好的深度捜索能力以便加快算法的收敛速度[20]，即要求惯性权重前期大、后期小。本文基于一种自适应非线性惯性权重，在满足惯性权重前期大、后期小的特点的同时，强调惯性权重的非线性变化，保证粒子间信息的充分交流。惯性权重更新公式如下[21]：

(10)

(11)

(12)

式中，ωmax=0.9；ωmin=0.4。

(3)自适应加速因子。加速因子c1决定粒子的自我学习能力，加速因子c2决定粒子的社会学习能力[22]。在反演初期，要求算法具备较强的社会学习能力，以扩大种群多样性，增加寻找全局最优解的可能；在反演后期，要求算法具备较强的自我学习能力，增强局部搜索能力，加快算法收敛速度。为满足该特点，将加速因子c1、c2作如下定义[23]：

(13)

(14)

式中，c1w=2.5；c1e=0.5；c2w=0.5；c2e=2.5。

本文改进粒子群算法流程如图2所示。

图2 改进算法流程Fig.2 Improved algorithm flow chart

4 算法评价

算法的评价指标包括时间复杂度和空间复杂度，本文对遗传类算法性能的评价主要从寻优成功率、平均收敛代数、平均最优适应度、最优个体进化曲线及抗噪性能等这些指标进行分析[24]。

表3记录了改进粒子群算法、基本粒子群算法、马奎特算法在30次随机计算中寻优成功率、平均收敛代数、平均最优适应度等指标的对比结果。从表3可知，改进粒子群算法比后2种算法寻优成功率更高。平均收敛代数越少代表收敛速度越快，改进粒子群算法是基本粒子群算法收敛速度的2倍。本文是基于预测值与实际值的差值构建的适应度函数，最优适应度越小，表示预测值越接近真实值，算法性能越优。改进粒子群算法比基本粒子群算法性能更优，稍逊于马奎特算法。

表3 算法指标对比Tab.3 Comparison of algorithm indexes

改进粒子群算法及另外2种算法最优个体进化曲线如图3所示。由图3可知，改进粒子群算法比基本粒子群算法迭代次数少很多，收敛速度显然更快；马奎特算法比粒子群类算法有更大的初值适应度，说明初始值对马奎特算法影响严重，对粒子群类算法影响很小。

图3 三种算法最优个体进化曲线对比Fig.3 Comparison of optimal individual evolution curves of three algorithms

对抗噪性能的评价，一般在测井响应中加入高斯白噪声，信噪比用式(15)表示[25]。

(15)

构建径向阶跃三参数正演模型，模型参数设置见表4。

表4 地层模型Tab.4 Strata model

先对6个模型分别加入5%、10%、20%的高斯白噪声，然后再进行反演处理，所得结果见表5。

表5 不同信噪比下地层模型反演结果对比Tab.5 Comparison of formation model inversion results under different signal-to-noise ratios

加入0%、5%、10%、20%的高斯白噪声进行反演处理后，平均相对误差分别为2.63%、3.10%、4.53%、8.71%。可以看出，改进粒子群算法抗噪性能较好，当加入20%的高斯白噪声时，依然可以较好地反演出原状地层信息，且反演后平均相对误差小于10%。

5 反演效果验证

实际泥浆侵入地层，造成地层电阻率分布并不是阶跃突变的，而是渐变的[26]。基于渗流—对流原理，可求得在不同时间内泥浆侵入地层后地层的电阻率径向分布特征，具体求解步骤可见文献[27]。但在实际反演处理中，若考虑泥浆的动态侵入，往往造成正演影响因素繁多，正演计算复杂且耗时长，反演结果甚至往往不收敛，故提出了简易的泥浆侵入径向阶跃介质模型，以此来简化反演过程中的正演计算。但这种反演模型的适用性及精度需与利用严格渗流—对流理论正演计算的结果进行对比，分析其误差，当误差在可控范围内时，此反演方法方可投入实际运用。不同泥浆侵入条件下，阵列侧向测井三参数反演结果与基于渗流—对流理论的正演结果对比如图4、图5所示。

图4 泥浆高侵下反演效果对比Fig.4 Comparison of inversion results under high mud penetration

图5 泥浆低侵下反演效果对比Fig.5 Comparison of inversion results under low mud penetration

实现基于渗流—对流理论的正演计算的基本参数设置如下：地层孔隙度10%，地层渗透率为5×10-3μm2，井眼与地层压力差为0.1 MPa。在泥浆高侵条件下，地层水矿化度为18 g/L，泥浆矿化度为5 g/L；在泥浆低侵条件下，地层水矿化度为18 g/L，泥浆矿化度为36 g/L。

从图4、图5可知，阵列侧向测井三参数反演的侵入带电阻率基本代表了渐变侵入带电阻率的加权平均和，可以在一定程度上反映侵入带电阻率特征。最为重要的是，改进粒子群算法可以准确地反演出地层真实电阻率，且比基本粒子群算法评价效果更好；同时，改进粒子群算法、基本粒子群算法都可以比较准确地反演出侵入带半径的大小，说明此方法可以在一定程度上反演出泥浆侵入后的地层径向电阻率分布剖面，具备现场应用价值。但此方法也有一定缺陷，在泥浆高侵过程中，实际地层电阻率分布会出现低阻环带特征，但三参数反演模型并不能反映出此特征。针对这种情况，后期精细解释时，可考虑进行阵列侧向测井四参数、五参数反演，提高反演模型适用性。