谢 源, 李 琼, 陈 杰, 何荣胜

(成都理工大学 地球物理学院，成都 610059)

在进行储层平面非均质性研究中，变差函数是最常见的研究方法。通过计算实验变差函数，然后对其进行拟合得到相关的参数，再对其进行地质解释。然而在实验变差函数拟合的过程中，由于理论变差函数参数众多且通常为非连续可导的，所以对其拟合一直是一个难点。最初对于变差函数的拟合都采用人工的方法[1]，但是这种方法得到的变差函数参数没有一个统一的评价标准，拟合的结果往往因人而异，很难得到最佳的拟合效果。当前，变差函数拟合的方法主要有最小二乘法[2]、极大似然法[3]、线性规划法[4]、加权多项式回归法[5]和遗传算法[6]等，但是各种方法都有其相应的局限性。最小二乘方法虽然简单易行，但是由于无法确定待拟合参数的正负号，故无法得到最优解；而极大似然法则计算复杂，而且有时不能得到最优解；加权多项式以及线性规划法依然不能解决参数正负号的问题；遗传算法虽然在对于求解非线性优化时具有全局寻优的特点，但是由于其可调参数过多且结构复杂从而影响了效率。

J.Kennedy教授和R.C.Eberhar教授在1995年提出了粒子群算法(Particle Swarm Optimization, 简称为PSO)[7]。它是一种基于群智能的全局寻优算法。其基本原理是模拟鸟群在觅食的过程中其个体间共享所了解的食物源信息，从而在准确的信息下使整个群体向着食物源的方向飞行。由于PSO算法收敛速度快、设置参数少，而且不受函数是否连续可导的条件的限制，使其在非线性拟合方面得到了广泛的应用。但是传统的PSO算法具有易陷入局部极值的缺陷，有时不能很好地解决问题。因此，本文提出了一种改进的PSO算法，并利用改进的算法实现了对实验变差函数进行自动拟合。

1 变差函数原理

1.1 变差函数的定义

变差函数是J.Motheron等[8]提出的一种矩估计方法。其表示为区域化变量的增量平方的数学期望，也就是区域化变量的增量的方差。根据地质统计学理论,通过半变差函数，利用适当的数理模型(如球状模型、高斯模型、指数模型等)，可对由自然地质规律形成的储层参数的空间相关性进行定量描述。下式为实验变差函数的计算公式

(1)

式中：N(h)表示距离为h的点个数；Xi为第i个点坐标；Z(Xi)、Z(Xi+h)分别为Xi和Xi+h两点处的观测值；h为两观测点间的距离；r*(h)为实验变差函数值。

1.2 模型的选择

变差函数的模型有多种，常用的模型主要有球状模型、指数模型、高斯模型、幂函数模型、对数模型等。各种模型的储层参数在同一方向上的变化速度不同，变程的表达形式也不同，变程的大小能反映储层参数沿某一方向的变化速度的快慢。基台值是反映储层参数在某一方向变化的幅度的大小，基台值越大，参数变化幅度越大，非均质性越强；反之，基台值越小，参数变化幅度也越小，非均质性越弱[9]。球状模型能够更快速地达到基台值，当达到基台值后，变程的增加不再会影响函数值，因此实际应用中多为球状模型。球状模型的表达式

(2)

式中：C0为块金常数；C为拱高；C0+C为基台值；a为变程。

2 PSO算法及其改进

2.1 标准粒子群算法

在PSO算法中，首先假设在D维搜索空间中，初始化M个粒子，将Xi=(Xi1,Xi2, …,XiD)记为第i个粒子的所在位置，可以根据目标函数计算得到Xi的适应度值，从而通过适应度值来判断其位置的优劣；将Vi=(Vi1,Vi2, …,ViD)记为第i个粒子的飞行速度；将pbi=(pbi1,pbi2, …,pbiD)记为第i个粒子所经过的位置中最优的位置；将gbi=(gbi1,gbi2, …,gbiD)定义当前群体经过的最优位置。然后粒子通过下式进行速度和位置的迭代更新，其更新公式如下

(3)

(4)

其中:d(1, …,D)是空间的维度；i(1, …,m)为种群的粒子数；R1、R2是介于(0,1)区间的随机数；k(1, …,n)是迭代次数；c1和c2为学习因子，其可以控制种群中个体向全体中其他优秀个体学习的能力，通常取c1与c2的值为2。式(3)中，第一项是粒子当前速度；第二项是粒子个体所经历的最好位置，这可以看作是粒子个体对自身经验的总结；第三项代表了种群的当前最优位置，可以看作是粒子间信息共享的体现。

2.2 改进粒子群算法

PSO算法在进行函数拟合时虽然简单易行，但是标准的PSO算法不能很好地提高算法的收敛性。其主要原因在于，PSO算法中粒子的收敛形式是以轨道的方式实现的；而且它在搜索过程中受最大速度的限制，使得PSO算法不是以概率1进行全局收敛。基于此，Sun等[10]提出了量子粒子群算法(Quantum-behaved Particle Swarm Optimization)。QPSO在一维空间中引入δ势阱，假设粒子在以P点为中心的δ势阱中，在粒子收敛过程中其速度和位置在量子空间中不能一起确定，其状态可以用波函数φ(Y)来表示，由式(5)确定,其中L是δ势阱的特征长度。

(5)

使用蒙特卡洛方法对粒子进行模拟，得到粒子的位置方程如下

(6)

其中:X为粒子的位置；P为吸引势的中心点；u是(0, 1)上服从均匀分布的随机数。由于L为时间t的函数，则量子粒子群算法中粒子的基本进化方程如下

(7)

然后我们将其运用于n维空间中,则吸引子Pi可以表示为Pi=(Pi1, …,Pin)，在每一维中用Pn为坐标，以Pi,j为中心建立一位δ势阱，则粒子i的波函数表示为

φ[xi,j(t+1)]=

(8)

位置方程

(9)

定义平均最好位置，其计算公式

(10)

则有

Li,j(t)=2α|mbestj(t)-xij(t)|

(11)

于是得到粒子的位置更新方程

(12)

其中:ui,j(t)～u(0,1);α为收缩扩张因子。

综上可以得出量子粒子群的更新方程如式(13)～(15)。

(13)

(14)

(15)

虽然由均匀分布产生的粒子具有较好的全局搜索能力，但是其局部搜索能力较弱。而柯西分布很容易得到距离原点较远的随机数，因此，采用柯西分布产生随机数可以使得算法很快跳出局部极值。标准QPSO算法中式(14)～(15)中参数φ1、φ2、u的值是由均匀分布产生的随机数来得到的。在改进算法中，我们采用柯西分布产生的随机数来确定φ1、φ2、u的值，并且采用线性变化的收缩扩张因子来对QPSO进行改进。收缩扩张因子α的取值如下

(16)

其中:tmax为最大迭代次数;t为当前迭代次数。

2.3 算法的实现流程

在变差函数拟合时，将所求参数一组解看作是一个粒子，将其个数看作是粒子的维度。对于一阶球状模型可以得到算法中粒子i的适应度函数如式(17)。

一阶球状模型

F(i)=

(17)

其中:hj为第j个滞后距；r*(h)为计算得到的实验变差函数值。

改进QPSO算法具体的拟合流程如图1。

3 仿真测试与实例计算

3.1 仿真测试

为了测试所改进的量子粒子群算法的性能，我们进行了一个简单的仿真测试。选择的测试函数及其性质如表1所示。

采用上述3个函数来对所改进的量子粒子群算法进行测试：模型一中设定粒子的数量为20、维度为10、迭代次数为 2 000次，分别用标准QPSO算法和改进量子粒子群算法运行30次并计算最优值的平均值；模型二中设定粒子数量为60、维度为20，迭代次数不变，运行30次并计算最优值的平均值。2种模型得到的最优值的平均值见表2。

图1 改进量子粒子群拟合流程Fig.1 Fitting process of the improved quantum particle swarm

函数名称表达式自变量范围最小值Sphere函数f(x)=∑Di=1x2i(-100, 100)0Rosenbrock函数f(x)=∑Di=1(100(xi+1-x2i)2+(xi-1)2)(-100, 100)0Rastrigrin函数f(x)=∑Di=1(x2i-10cos(2πxi)+10(-5.12, 5.12)0

从测试的结果我们可以看出，改进的量子粒子群算法相比较于标准的QPSO算法所求得的函数极值更加接近于0。其在求解函数的极值的精度上有明显的提高，特别是对高维函数其优越性更加明显。

表2 模型测试结果Table 2 Model test results

3.2 实例计算

对塔河油田某井的声波时差数据按照式(2)计算实验变差函数值。由于所选层段的非均质性较强，在实际计算过程中发现，当滞后距为30 m时计算得到的实验变差函数的稳健性稍好。其计算结果见表3。

首先使用传统的加权最小二乘拟合对得到的实验变差函数进行拟合，拟合模型选择球状模型。得到拟合参数值为：C0=4.810 9,C=7.923 3,a=294.533 9。

在使用改进的QPSO算法拟合时，我们选择一阶球状模型，粒子的群体规模为30，维度设置为3，迭代次数设置为 1 500次，其中参数的取值范围分别设置为：C0=[0,15],C=[0,15],a=[0,500]。独立运行程序10次，得到的待拟合的参数取值为：C0=4.325 4,C=7.342 4,a=233.416 5，此时得到F(i) =37.75。

然后将2种方法得到的参数代入球状模型得到理论变差函数值(表3)。

表3 塔河油田某井的实验变差函数及理论变差函数值Table 3 Experimental variogram and theoretical variogram value from a drilling well in Tahe area

h为滞后距；r(h)为计算得到的实验变差函数值；r*(h)1为改进的量子粒子群算法所得到的理论变差函数值；r*(h)2为最小二乘算法得到的理论变差函数值

2种拟合方法得到的变差函数如图2，从图中可以看出，蓝色的折线为计算所得的实验变差函数值，红线为利用改进的量子粒子群算法拟合得到的理论变差函数曲线，绿线为用加权最小二乘拟合得到的理论变差函数曲线。从中我们可以看出利用改进的量子粒子群算法拟合得到的理论变差函数曲线比加权最小二乘得到的曲线具有更小的变程，即其非均质性更强，拟合结果也更加贴合计算得到的实验变差函数曲线，所以本文提出的改进量子粒子群算法对变差函数的拟合取得了比较理想的效果。

图2 实验变差函数拟合结果Fig.2 Results of the experimental variogram fitting

我们引入拟合度和绝对偏差来对2种拟合方法的结果进行评价，其计算公式

(18)

(19)

分别计算改进的量子粒子群算法和加权最小二乘拟合得到的理论变差函数值的拟合度和绝对偏差见图3。从图中可以看出，采用改进的量子粒子群算法拟合得到的理论变差函数的拟合度和绝对偏差都是小于由加权最小二乘方法所得到的理论变差函数的。因此，本文的方法可以用于实验变差函数的拟合。

3.3 储层平面非均质性预测

对某储层的孔隙度数据在4个方向(0°, 45°, 90°, 135°)上分别计算实验变差函数，然后采用改进QPSO算法对实验变差函数进行拟合。在拟合时，设定粒子数量为50，维度为3，迭代次数为2 000次，得到的拟合效果如图4。

由图4可以看出，在4个方向上该储层的孔隙度变差函数具有较明显的差异，表明储层存在非均质性。图中可以看出，在0°和90°方向，变程a的值明显小于45°和135°方向的变程值，这表明在这2个方向上该储层具有较强的非均质性。在45°和135°方向，变程a值较大，表明该方向上非均质性较弱。

4 结 论

a.本文采用柯西分布产生随机数并采用线性变化的收缩扩张因子的方法来改进QPSO算法，对改进算法进行了仿真，结果表明，改进的QPSO算法比标准的QPSO算法在求解函数的极值上具有更高的精度。

b.运用改进的量子粒子群算法对实验变差函数进行拟合，相比较于传统的加权最小二乘拟合变差函数的方法，改进算法得到的结果的拟合度和绝对偏差更小，拟合的结果更接近于实验计算值，是一种更优的变差函数拟合方法。

图3 两种方法的拟合度和绝对偏差对比Fig.3 Contrast of the fitting and absolute deviations of the two methods

图4 某储层4个方向孔隙度变差函数Fig.4 Porosity variation function in four directions of a reservoir

c.运用改进QPSO方法对塔河油田某储层的平面非均质性进行研究，对储层孔隙度实验变差函数进行拟合，从拟合结果看出该储层在0°和90°方向具有较强的非均质性，在45°和135°方向非均质性较弱。

[参考文献]

[1] 王仁铎，胡光道.线性地质统计学[M].北京：地质出版社，1998.

Wang R D, Hu G D. Linear Geostatistics[M]. Beijing: Geological Publishing House, 1998. (in Chinese)

[2] 李正文，李琼，吴朝容.沉积盆地有效储集层综合识别技术[M]．成都：四川科学技术出版社，2002.

Li Z W, Li Q, Wu C R. Comprehensive Identification Technology of Effective Reservoirs in Sedimentary Basins[M]. Chengdu: Sichuan Science and Technology Press, 2002. (in Chinese)

[3] Cressie N. Fitting variogram models by weighted least squares[J]. Journal of the International Association for Mathematical Geology, 1985, 17(5): 563-586.

[4] Mardia K V, Marshall R J. Maximum likelihood parameter estimation of models for residual covariance in spatial regression[J]. Biometrika, 1984, 71(1): 135-138.

[5] 矫希国，刘超.变差函数的参数模拟[J].物探化探计算技术，1996，16(2)：157-161．

Jiao X G, Liu C. Estimation of variation parameter[J]. Computing for Geophysical and Geochemical Exploration, 1996, 16(2): 157-161. (in Chinese)

[6] 王仁铎.用加权多项式回归进行球状模型变差图的最优拟合[J].地球科学，1986，11(2)：143-151.

Wang R D. On the best fitting of the spherical model variogram with the weighted polynomial regression[J]. Earth Science, 1986, 11(2): 143-151. (in Chinese)

[7] 冯国庆，周涌沂，李允，等.利用遗传算法拟合实验变差函数[J].西南石油学院学报，2005，27(5)：23-25．

Feng G Q, Zhou Y Y, Li Y,etal. Application of genentic algorithm to match theoretical variogram[J]. Journal of Southwest Petroleum Institute, 2005, 27(5): 23-25. (in Chinese)

[8] Kennedy J, Eberhart R C. Particle swarm optimization[C]//Proceeding of IEEE International Conference on Neural Networks. Piscataway: IEEE Press, 1995: 1942-1948．

[9] Matheron G. Principles of geostatistics[J]. Economic Geology, 1963, 58(8): 1246-1266.

[10]Sun J, Xu W, Feng B. Adaptive parameter control for quantum-behaved particle swarm optimization on individual level[C]//2005 IEEE International Conference on Systems, Man and Cybernetics. IEEE, 2005, 4: 3049-3054.

敬告作者

为适应我国科技信息化建设需要，扩大作者学术交流渠道，本刊已加入《中国学术期刊(光盘版)》和《中国知网》(http://www.cnki.net)、万方数据电子出版社的《万方数据——数字化期刊群》(http://www.wanfangdata.com.cn)、教育部科技发展中心的《中国科技论文在线》、重庆维普资讯有限公司的《中文科技期刊数据库》、华艺数位艺术股份有限公司的《CEPS中文电子期刊》、北京书生网络技术有限公司的《书生数字期刊》、北京世纪超星公司的“域出版”平台、中教数据库等。作者著作权使用费与稿酬由本刊一次性给付。如果作者不同意将文章编入上述数据库，请在来稿时声明，本刊将做适当处理。

《成都理工大学学报(自然科学版)》编辑部