粒子群优化算法叠前非线性反演方法

摘要

随着我国经济持续稳定的增长，对能源的需求也持续快速增加，我们将面临石油天然气资源短缺的严峻考验。储层预测是石油勘探生产中的重要研究环节，优质储层的发现需要充分利用地震资料所蕴含的地下地层的构造、岩性等信息，而这些信息往往要通过地震反演得到。但常规反演基于线性简化的模型，并且有较多的人为因素干扰，往往得不到高精度的储层参数。本发明对粒子群优化算法这一课题展开研究，以粒子群优化算法的思想和原理为基础，将左普利茨(Zoeppritz)方程进行推导简化，使其适于粒子群优化算法的求解，然后将粒子群优化算法用于叠前反演。分别对二维理论模型及海上地震资料进行叠前反演取得令人满意的结果。

说明书

技术领域

本发明涉及油气田勘探技术领域，属于地震资料反演范畴，具体的说是通过粒子群优化算法解决非线性的地震反演问题。 

背景技术

利用地震资料进行储层岩性识别和流体预测一直是地球物理学家努力追求的目标。地震叠前反演即AVO反演的理论基础是著名的左普利茨(Zoeppritz)方程组。Ostrander(1984)在研究“亮点型”砂岩储层地震振幅特征过程中，发现了“含气砂岩反射振幅随偏移距增加而增大，含水砂岩反射振幅随偏移距增加而减小”的现象，这一现象极大的改善了烃类检测的能力，将人们的视线从叠后引向叠前，标志着实用AVO技术的出现。但是，左普利茨方程十分复杂，物理含义也不明确。国内外地球物理学家对其进行了很多形式的研究和简化，得到了一系列近似公式。Bortfeld(1961)利用地层厚度趋于零来逼近单界面的方法计算了平面纵波和透射波的反射系数，给出了区分流体和固体的简化公式。Aki和Richards于1980年提出的近似方程侧重描述了纵、横波速度和密度的变化对反射系数的影响。Shuey(1985)近似方程则将反射系数表示成法线入射与近、中、远不同入射的三项之和的关系式，较直观的反映了振幅与入射角的关系。Smith和Gidlow(1987)在对含气砂岩进行AVO分析时，提出了AVO分析的加权叠加处理方法，并引入了“流体因子”概念。Fatti等(1994)从Aki和Richards方程出发，在弱化地层密度项的前提下，提出了较为精简的公式用于预测含气砂岩的AVO响应。Connolly(1999)提出的弹性阻抗(EI)概念和数学表达式，极大地丰富了AVO分析技术的内涵，扩大了AVO反演技术的外延，将传统的纵波阻抗反演推广到分角度叠加数据。Ozdemir等(2001)提出了利用AVO信息进行弹性参数同步反演，称为“联合反演”。马劲风等(2003)提出了广义弹性波阻抗的概念。近几年，AVO分析和随机反演技术方法不断发展，叠前同步反演方法日益成熟，在同步反演方法中多采用非线性反演技术。AVO的提出最初仅仅是为了提高烃类检测能力，今天AVO的发展早已超出了这个范畴，渗透到地震勘探的各个领域。在裂缝检测、压力预测、油藏动态检测、油气预测、储层非均质性描述方面得到了广泛应用。 

粒子群优化算法(Particle Swarm Optimization，PSO)，又称微粒群算法，是群智能算法(Swarm Intelligence，SI)的一种。最先是由美国学者Kennedy和Eberhart于1995年提出。粒子群算法是受鸟类群体的防御、捕食行为中的搜索策略启发而形成的。自然界中许多生物具有一定的群体行为，如鸟群、鱼群等，虽然群体中单个个体只具有简单的行为规则，但是组成的群体行为却非常复杂。很多科学家对鸟群或者鱼群的群体行为进行了研究，包括计算机仿真。粒子群优化算法从这种模型中得到启示并用于解决优化问题，将优化问题的解看作是搜索空间中的点，称之为“粒子”。每个粒子的运动根据自己和其它粒子的“飞行经验”寻优，从而达到全空间搜索最优解的目的。许多学者和研究人员在基本PSO算法的基础上在参数选择、拓扑结构以及与其他优化算法相融合方面，提出了很多改进的PSO算法。Kennedy和Eberhart于1997年提出了二进制粒子群(Binary Particle Swarm Optimization，BPSO)算法，这是基于连续空间的离散粒子群(Discrete PSO，DPSO)算法。Clerc(2000)针对旅行熵问题(Traveling Salesman Problem，TSP)提出了TSP-DPSO算法，是基于离散空间的DPSO。Jun Sun等(2003)将量子行为引入到粒子群算法中，提出了量子粒子群算法(Quantum PSO，QPSO)。2004年高鹰提出了基于模拟退火的粒子群算法(SA-PSO)。P.S.Shelokar等人于2007年提出了基于蚁群和粒子群算法的混合算法，即PSACO(Particle Swarm Ant Colony Optimization)算法。周雅兰等(2008)把分布估计算法思想引入到PSO算法中，提出基于分布估计的离散粒子群(Estimation of Distribution PSO，EDPSO)算法。林东毅等(2008)将免疫机制引入到BPSO算法中，基于决策表差别矩阵的某种属性重要性度量作为疫苗模式，提出了一种基于免疫粒子群优化的最小属性约简算法，称为免疫离散粒子群算法(IPSO算法)。粒子群算法近年来发展很快，被成功地应用于函数寻优、神经网络训练、模式识别分类、模糊系统控制以及工程等众多领域，大量实际应用证明其是有效的。它有着较好的发展前景，值得做进一步的研究。 

发明内容

本发明利用保幅处理的高分辨率地震资料进行非线性反演，在解决非线性问的过程中使用了更新的反射系数公式，并使用了粒子群算法进行寻优。这一措施可以减少人为误差，有效的提高反演的精度。 

首先做如下定义： 

在地球物理学中，连续函数x(t)的反射系数的定义为， 

$\underset{\mathrm{\Delta t}\to 0}{R_{x\left(t\right)}}=\frac{\mathrm{\Delta x}}{2\overline{x\left(t\right)}}$

以dt的采样间隔将连续函数x(t)采样成离散的函数，那么， 

$R_x=\frac{x_2-x_1}{x_2+x_1}$

其中x可以为纵波速度α、横波速度β、密度ρ。下标1表示下伏地层的参数，下标2表示上覆地层的参数。定义横纵波速度比 

$K=\frac{{\beta }_1}{{\alpha }_1}$

本专利使用保幅处理得到的反射角域共反射点道集作为输入，实现上述目的采取的技术方案如下： 

步骤1：对于地下任意位置。首先给定纵波速度反射系数、横波速度反射系数、密度反射系数、以及横纵波速度比的取值范围，设定最大迭代次数以及种群的大小。给定惯性权重和学习因子，设定容许的误差范围。 

其中惯性权重、学习因子、种群大小有如下意义： 

惯性权重：勘探能力和开发能力的平衡是影响优化算法性能的一个重要方面。对于粒子群优化算法来说，这两种能力的平衡是靠惯性权重w来实现的。较大的惯性权重使粒子在自己原来的方向上具有更大的速度，从而在原方向上飞行更远，具有更好的勘探能力；较小的惯性权重使粒子继承了较少的原方向的速度，从而飞行较近，具有更好的开发能力。 

学习因子：学习因子为非负常数，代表粒子偏好的权值，使粒子具有自我总结和向群体中优秀个体学习的能力，从而向群体内或邻域内最优点靠近。Kennedy认为，两个学习因子之和应为4.0左右，此时的搜索效果比较好。通常的做法是将他们都设为2.05。 

群体大小：当群体大小设定的很小时，陷入局部极优的可能性很高；当群体大小为一时，粒子群优化算法变为基于个体搜索的方法，一旦陷入局优，将不可能跳出；当群体大小很大时，粒子群优化算法的优化能力更好，但会导致计算时间大幅增加，并且当群体数目增长至一定水平时，继续增长将不再有显著的作用。至于究竟多少粒子参加搜索能够取得理想的效果，前人通过使用多个基准函数对数量不同的种群计算其平均适应度，认为种群数量保持在30左右时搜索效率较好。 

步骤2：对各个粒子的纵波速度反射系数、横波速度反射系数、密度反射系数、以及横纵波速度比用随机数进行初始化，并设定初始速度为零，令粒子的个体最优值与初始值相同，并计算各个粒子的适应度，选择适应度最小的粒子为当前的全局最优解。 

粒子适应度的计算方法为： 

首先定义矩阵 

$M^{+}=\left(\begin{array}{cccc}D& \sqrt{1-{\left(\mathrm{KD}\right)}^2}& -\mathrm{AD}& \sqrt{1-{\left(\mathrm{KBD}\right)}^2}\\ \sqrt{1-D^2}& -\mathrm{KD}& \sqrt{1-{\left(\mathrm{AD}\right)}^2}& \mathrm{KBD}\\ 2D\sqrt{1-D^2}& \frac{1-2{\left(\mathrm{KD}\right)}^2}{K}& -2{\mathrm{CB}}^{2}D\sqrt{1-{\left(\mathrm{AD}\right)}^2}& \frac{\mathrm{BC}[1-2{\left(\mathrm{KBD}\right)}^2]}{K}\\ 1-2{\left(\mathrm{KD}\right)}^2& -2K^{2}D\sqrt{1-{\left(\mathrm{KD}\right)}^2}& -\mathrm{AC}[1-{\left(\mathrm{KD}\right)}^2]& -2{\mathrm{CB}}^2{\mathrm{DK}}^2\sqrt{1-{\left(\mathrm{KBD}\right)}^2}\end{array}\right)$

式中的A、B、C、D以及K用以下方式定义 

$A=\frac{{\alpha }_2}{{\alpha }_1}=\frac{1+R_{\alpha }}{1-R_{\alpha }}$

$B=\frac{{\beta }_2}{{\beta }_1}=\frac{1+R_{\beta }}{1-R_{\beta }}$

$C=\frac{{\rho }_2}{{\rho }_1}=\frac{1+R_{\rho }}{1-R_{\rho }}$

D＝sini₁

$K=\frac{{\beta }_1}{{\alpha }_1}$

其中i₁为入射角。再定义 

$M^{-}=\left(\begin{array}{cccc}-D& \sqrt{1-{\left(\mathrm{KD}\right)}^2}& -\mathrm{AD}& \sqrt{1-{\left(\mathrm{KBD}\right)}^2}\\ \sqrt{1-D^2}& -\mathrm{KD}& \sqrt{1-{\left(\mathrm{AD}\right)}^2}& \mathrm{KBD}\\ 2D\sqrt{{1-D}^2}& \frac{1-2{\left(\mathrm{KD}\right)}^2}{K}& -2{\mathrm{CB}}^{2}D\sqrt{1-{\left(\mathrm{AD}\right)}^2}& \frac{\mathrm{BC}[1-2{\left(\mathrm{KBD}\right)}^2]}{K}\\ 2{\left(\mathrm{KD}\right)}^2-1& -2K^{2}D\sqrt{1-{\left(\mathrm{KD}\right)}^2}& -\mathrm{AC}[1-{\left(\mathrm{KD}\right)}^2]& -2{\mathrm{CB}}^2{\mathrm{DK}}^2\sqrt{1-{\left(\mathrm{KBD}\right)}^2}\end{array}\right)$

则由左普利茨方程可以得到以下形式的解 

$R_{\mathrm{pp}}=\frac{\left|M^{-}\right|}{\left|M^{+}\right|}$

由这一表达形式可以得到粒子适应度的表达式为 

E＝||w(θ)*R_pp(R_α，R_β，R_ρ，K，θ)-d(θ)|| 

式中E为粒子适应度，w(θ)为地震子波，d(θ)为地震观测数据。 

步骤3：比较当前的全局最优的适应度是否达到了误差容许范围。如果达到误差容许范围就停止计算，输出全局最优为反演结果。否则进入步骤4。 

步骤4：更新粒子的速度，并用新的粒子速度更新粒子的位置。迭代次数增加一次，若迭代次数超过最大迭代次数，则停止计算，输出全局最优。 

粒子速度以及位置的更新方法为： 

v_id(t+1)＝wv_id(t)+c₁r₁(p_id(t)-x_id(t))+c₂r₂(g_d(t)-x_id(t)) 

x_id(t+1)＝x_id(t)+v_id(t+1) 

式中，v为速度，x为粒子位置。m为种群大小，n为粒子的维数，1≤i≤m，1≤d≤n，w称为惯性权重因子；c₁、c₂称为学习因子或加速因子，其中，c₁为调节粒子飞向自身最好位置方向的步长，c₂为调节粒子飞向全局最好位置的步长；r₁、r₂为[0，1]内的随机数；t为当前迭代代数。 

步骤5：重新计算粒子的适应度，检测各个粒子的适应度是否小于更新前的适应度。若小于则更新粒子的个体最优为粒子更新后的位置。 

步骤6：比较新的个体最优，选择适应度最小的个体最优作为全局最优。判断全局最优是否达到误差容许范围，达到误差容许范围就输出全局最优为反演结果，否则回到步骤4。直到处理完地下所有计算点。 

以上具体实施方式仅用于说明本发明，而非用于限定本发明。 

附图说明

图1是粒子群算法粒子位置更新的示意图。 

图2是一个合成的反射角域的共反射点道集。 

图3是图2所示道集使用本方法的反演得到的纵波反射系数与模型纵波反射系数的对比。 

图4(a)是图2所示道集用本方法的反演得到的纵波速度与模型纵波速度的对比； 

图4(b)是图2所示道集用本方法的反演得到的横波速度与模型横波速度的对比； 

图4(c)是图2所示道集用本方法的反演得到的密度与模型密度的对比； 

图4(d)是图2所示道集用本方法的反演得到的横纵波速度比与模型横纵波速度比的对比。 

图5是本方法在某地区实际资料中反演得到的纵横波速度比(右)与商业软件得到的纵横波速度比(左)的对比。 

具体实施方式

通过一个合成道集来说明： 

首先，利用左普利茨方程对该模型进行正演，每隔3°计算1°到30°范围内各地层的反射系数。计算出的10个反射系数序列分别与40Hz的零相位雷克子波进行褶积，得到的合成地震记录，即角道集数据如图2所示。 

将这十个角道集数据作为输入进行反演。以求解计算出的纵波反射系数R_pp与实际地震记录之间误差函数E最小作为目标函数。采用粒子群优化算法叠前反演的步骤如下： 

(1)将由模型正演得到的角道集记录作为观测值(真值)。 

(2)利用粒子群优化算法生成各角道集的反射系数序列如r(t)候选解。每个粒子的位置向量即为r(t)候选解，每个粒子的维数为4，即4个自变量纵波反射系数R_α，横波反射系数R_β，密度反射系数R_ρ，以及横纵波速度比K。 

(3)按照步骤2计算适应度函数E。将计算得到R_pp的与角道集数据比较，计算适应度函数E。 

(4)将求取的适应度函数E代入全局最优解G_best，并按照步骤4对每个粒子进行位置更新。 

(5)利用粒子群算法循环迭代，对最优解不断进行优化，直到满足终止条件(如E小于某截断误差，或达到最大迭代次数)，则此时个体最优与全局最优一致，所有粒子收敛于一点，对应的位置向量即为要求的纵波反射系数R_α，横波反射系数R_β，密度反射系数R_ρ，以及横纵波速度比K。 

在本例中利用粒子群算法反演数值模型角道集数据的反射系数时，参数取值如下：误差函数小于1e-6，种群大小m＝30，惯性因子w＝0.1，学习因子c₁＝c₂＝2.04，粒子位置取值的最大值为[0.1，0.1，0.1，0.75]，最小值为[-0.1，-0.1，-0.1，0.35]，分别对应纵波反射系数R_α，横波反射系数R_β，密度反射系数R_ρ，以及横纵波速度比K的限制范围。 

反演出的反射系数与原反射系数的对比如图3所示。可以看出，反演出的反射系数(虚线)与真实的反射系数(圆圈)吻合的很好。通过反演得到的反射系数R_α、R_β、R_ρ，以及K，可以计算各层的纵波速度、横波速度、密度等弹性参数(实线)，与模型的弹性参数(虚线)对比如图4所示。由图4看出，纵波数据和密度数据与真实情况完全吻合，只是横波数据和纵横波速度比有极小的误差，产生这些误差的原因一是因为在反演过程中，横波分量的精确程度受初值设定、岩石物理关系等因素的影响，稳定性不如纵波分量；二是由于反演过程中计算机的截断误差导致。为了使误差尽量减小，在可能的情况下，不仅使用纵波数据作为反演的观测值，还可以使用转换横波的信息。甚至在拥有VSP资料的情况下，可将反射纵波、反射横波、透射纵波、透射横波的数据均作为反演的观测值，那么反演得到的各种反射系数以及由此推导出的岩石的弹性参数将会更加准确。 

参考文献 

Zoeppritz K，Erdbebenwellen VIIIB.On the reflection and propagation of seismic waves.Gottinger Naehriehten，1919，66-84 

Ostrander W J.Plane waves reflection coefficient s for gas sands at normal angles of incidence.Geophysics，1984，49(10)：1637-1648. 

Bortfeld R.Approximation to the reflection and transmission coefficients of plane longitudinal and transverse wave.Geophysics Prospecting，1961，9(4)：485-502. 

Aki K，Richards P G.Quantitative seismology：Theory and methods.W.H.Freeman and Co，1980，100-170. 

Shuey.A simplification of the Zoeppritz equations.Geophysics，1985，50(4)：609-634. 

Smith G C，Gidlow P M.Weighted stacking for rock property estimation and detection of gas.Geophysical Prospecting，1987，35：993-1014. 

Fatti J L，Vail P J，Smith G C，et al.Detection of gas in sandstone reservoirs using AVO analysis：A case seismic case history using the Geostack technique.Geophysics，1994，59(5)：1362～1376. 

Connolly P.Elastic Impedance.The Leading Edge，1999，18(4)：438～452. 

Ostrander W J.Plane-wave reflection coefficients for gas sands at non-normal angels of incidence.Exploration Geophysics，1984，49：1637-1648. 

马劲风.地震勘探中广义弹性阻抗的正反演.地球物理学报，2003，46(1)：118-124. 

杜焰.AVO技术及在苏里格气田储层含气性检测中的应用：(硕士学位论文).成都.成都理工大学，2009. 

Kennedy，Eberhart R C.Particle Swarm Optimization.International conference on Neural Networks，1995，1942-1948. 

刘波.粒子群优化算法及其工程应用.北京：电子工业出版社，2010. 

Kennedy J，Eberhart R C.A discrete binary version of the Particle Swarm Algorithm.The 1997 Conference on System，Cybernetics and Informatics，1997，4104-4109. 

Clerc M.Discrete Particle Swarm Optimization，Illustrated by the Traveling Salesman Problem.http://www.mau-ricecierc.net，2000. 

Sun J，Feng B，Xu W.Particle Swarm Optimization with Particle Having Quantum Behavior.IEEE Proc.Of Congress on Evolutionar Computation，2004，325-331. 

高鹰，谢胜利.基于模拟退火的粒子群优化算法.计算机工程与应用，2004(1)：47-50. 

Shelokar P S，Siarry P，Jayaraman V K，Kulkarni B D.Particle swarm and ant colony algorithms hybridized for improved continuous optimization.Applied Mathematics and Computation，2007，188：129-142. 

周雅兰，王甲海等.一种基于分布估计的离散粒子群优化算法.电子学报，2008，36(6)：1242-1248. 

林东毅，廖建坤.最小约简问题的一个免疫离散粒子群算法.小型微型计算机系统，2008，29(6)：1089-1092. 

贺懿.地震储层参数非线性反演与预测方法研究：(博士学位论文).青岛.中国海洋大学，2008。