CN104570101A - 一种基于粒子群算法的avo三参数反演方法 - Google Patents

Abstract

本发明提供了一种基于粒子群算法的AVO三参数反演方法，属于石油地球物理勘探领域。本方法包括：(1)输入叠前道集数据d(θ)；(2)求取纵波反射系数R_pp；以所述纵波反射系数R_pp与子波w(θ)褶积同叠前道集数据d(θ)之间的误差函数E作为目标函数；(3)使用粒子群反演的方法寻找到使所述目标函数的值为最小的纵波速度、横波速度以及密度，即为所求的三参数。本发明方法降低了对初始模型的依赖并提高了反演的精度，而且具有较强的抗噪声能力。

Description

一种基于粒子群算法的AVO三参数反演方法

技术领域

本发明属于石油地球物理勘探领域，具体涉及一种基于粒子群算法的AVO三参数反演方法，利用粒子群算法进行AVO反演得到纵波速度、横波速度以及密度值。

背景技术

利用地震资料进行储层岩性识别和流体预测一直是地球物理学家努力追求的目标，AVO分析则揭示了重复观测的反射振幅变化规律与地层弹性参数的关系。AVO反演就是利用振幅随炮检距的变化规律估算储层岩石及流体的性质并进行油气预测，是提取隐藏在地震信息中的岩性参数的重要途径，其理论基础是著名的左普利茨(Zoeppritz)方程组。

AVO反演实际上可以看做是弹性波反演的简化形式。早在1984年，Kennett基于传播矩阵的方法，建立了反射系数模型的反演方法；Tarantola(1986)采用先验约束信息对反演问题进行约束；Crase等(1990)实现了实际地震数据的非线性弹性波波形反演，并采用了四种最小化准则分别对目标函数进行最小化；Sen和Stoffa(1995)在深度域，基于传播矩阵理论，采用回归方法反演求解弹性波方程。由于计算复杂，花费巨大，这些弹性波非线性反演方法仅在理论上比较成熟，在实际中并没有广泛应用。

Zoeppritz方程的数学表达式较为复杂，物理上也不够直观，实现起来较为困难，如果采用近似表达式，弹性波反演问题可简化近似为更为容易求解的问题，即根据振幅随入射角的变化关系，由实际地震道集记录估算岩石的地震参数。Wiggins等、Lortzer和Berkhout(1993)忽略了地震数据的有限频宽，在动校正叠前偏移道集上反演带限反射系数，并采用叠后阻抗反演方法对反射系数反演，求取阻抗及密度；Smith和Gidlow(1987)采用Gardner(1974)的经验关系式对Aki和Richard(1980)近似公式进行重新推导，得到一个没有密度项的近似公式，利用最小二乘法进行迭代反演；Simmon及Backus(1996)研究了AVO波形反演，计算了波的旅行时，在反演中建立了参数协方差矩阵约束反演，反演最终得到带限的纵横波阻抗反射系数及密度反射系数剖面；Grion(1998)忽略了多次波及转换波，做了旅行时和反射系数联合反演；Conno11yP于1999年将声波阻抗和AVO反演推广到弹性波阻抗反演，信噪比得到显著改善，也使该方法得到广泛应用；马劲风(2003)研究了广义弹性阻抗的正反演方法；Boland及More(2003)基于贝叶斯的理论，利用叠前地震数据从模型参数与观测数据的联合分布推导出纵横波速度及密度的反演方程；Downton(2002)在动校前道集上做AVO稀疏脉冲波形反演，得到反射系数剖面，提高了分辨率；Hampson和Russell(2005)根据叠后波阻抗公式推导的思想，给出了叠前AVO同步反演公式；陈建江等(2007)针对参数估计的不确定性提出了一系列改进方法，运用贝叶斯理论将基于地质环境控制的鲁棒统计分析融入反演问题当中，建立了测井数据的参数协方差约束反演，提高了反演问题的稳定性。

经过近30年的发展，AVO技术已成为应用叠前资料预测油气的有效方法，在裂缝检测、油藏动态监测、油气预测、储层非均质性描述等方面得到广泛应用。

粒子群优化算法(Particle Swarm Optimization，PSO)，又称微粒群算法，是群智能算法的一种。最先是由美国学者Kennedy和Eberhart(1995)提出。算法源于对鸟群捕食行为的研究，通过对简单社会***的模拟，在多维解空间中构造“粒子群”，粒子群中每个粒子通过跟踪自己和群体所发现的最优值，修正自己的前进方向和速度，从而实现寻优。群体中每个粒子表示问题的一个可行解，并具有与目标函数相关的适应度值。粒子在搜索空间中以一定的速度飞行，并根据自身的飞行经验以及当前最优粒子的状态对速度进行动态调整，个体之间通过协作与竞争，实现对问题最优解的搜索。

基本粒子群算法中粒子向自身历史最佳位置和领域或群体历史最佳位置聚集，形成了粒子种群的快速趋同效应，容易出现陷入局部极值、早熟收敛或停滞现象。同时，它的性能也依赖于算法参数。为了克服上述不足，许多学者提出了各种改进措施，归纳总结出粒子群的有关进展主要集中在粒子群初始化、领域拓扑、参数选择以及混合策略等方面(戴朝华，2010)。如M.Richards等(2004)提出了基于centroidal voronoi tessellations(CVTs)的种群初始化方法；薛明志等(2005)采用正交设计方法对种群进行初始化；StefanJanson等(2005)提出等级PSO，Ben Niu等(2007)采用主-仆模型(master-slavermodel)，主群体在仆群体提供的最佳位置基础上开展搜索；Ratnaweera A等(2004)提出自适应时变调整策略；P.Moore等(2006)利用差分进化(DE)操作选择速度更新公式中粒子最佳位置；Fan S K等(2004)将算法与单纯形法结合，Xi-Huai Wang等(2004)将算法与模拟退化法结合；Junying Chen等(2005)将算法与爬山法结合，孙艳霞等(2008)提出一种不含随机参数、基于确定性混沌Hopfield神经网络群的粒子群模型等等。

粒子群算法近年来发展很快，被成功地应用于函数寻优、神经网络训练、模式识别分类、模糊***控制以及工程等众多领域，大量实际应用证明其是有效的。它有着较好的发展前景，值得做进一步的研究。

地震叠前反演能够保留地震反射振幅随炮检距或入射角变化的特征，可以得到高精度的、能够反映储层横向变化的多种地层弹性参数，对于研究复杂储层的空间分布及复杂油气藏的精细描述具有十分重要的意义。

发明内容

本发明的目的在于解决上述现有技术中存在的难题，提供一种基于粒子群算法的AVO三参数反演方法，将非线性反演的思想引入叠前反演过程中，形成一种新型的非线性反演方法-粒子群优化算法(Particle Swarm Optimization，PSO)。基于粒子群优化算法以及AVO三参数反演方程，本发明提出了一种叠前反演新方法，其不依赖于初始模型，不依赖于解释层位，克服了常规叠前反演强烈依赖初始模型，易陷入局部极优的缺点，是一种切实可行的叠前反演新方法。

本发明是通过以下技术方案实现的：

一种基于粒子群算法的AV0三参数反演方法，包括：

(1)输入叠前道集数据d(θ)；

(2)求取纵波反射系数R_pp;以所述纵波反射系数R_pp与子波w(θ)褶积同叠前道集数据d(θ)之间的误差函数E作为目标函数；

(3)使用粒子群反演的方法寻找到使所述目标函数的值为最小的纵波速度、横波速度以及密度，即为所求的三参数。

所述步骤(2)中是利用公式(3)求得纵波反射系数R_pp：

$R_{\mathrm{pp}}=\frac{1}{2}(1-4\frac{{\mathrm{\β}}^2}{{\mathrm{\α}}^2}{\sin }^2\mathrm{\θ})\frac{\mathrm{\Δ\ρ}}{\mathrm{\ρ}}+\frac{{\sec }^2\mathrm{\θ}}{2}\frac{\mathrm{\Δ\α}}{\mathrm{\α}}-4\frac{{\mathrm{\β}}^2}{{\mathrm{\α}}^2}{\sin }^2\mathrm{\θ}\frac{\mathrm{\Δ\β}}{\mathrm{\β}}---\left(3\right)$

其中，R_pp为纵波反射系数，α为纵波速度，β为横波速度，ρ为密度，θ为角度；

Δα=α₂-α₁，Δβ=β₂-β₁，Δρ=ρ₂-ρ₁

α=(α₁+α₂)／2，β=(β₁+β₂)／2

θ=(θ₁+θ₂)／2，ρ=(ρ₁+ρ₂)／2

下标1，2分别表示第一层介质和第二层介质。

所述步骤(2)中的误差函数E如下：

$E=\left|\right|w\left(\mathrm{\θ}\right)*R_{\mathrm{pp}}(V_p,V_s,\mathrm{\ρ},\mathrm{\θ})-d\left(\mathrm{\θ}\right)\left|\right|\⇒\min ---\left(4\right)$ 其中,R_pp(V_p,V_s，ρ,θ)·是利用公式(3)得到;其中vp为纵波速度，与公式(3)中的α含义相同，v_s为横波速度，与公式(3)中的β含义相同。

所述步骤(3)包括：

(31)将纵波速度v_p，横波速度v_s，密度ρ简记为r(t)，每个粒子的位置向量即为r(t)候选解，每个粒子的维数为3，即3个自变量：纵波速度v_p、横波速度v_s和密度ρ；粒子的位置x_i对应r(t)的值，粒子的速度v_i对应r(t)的变化值；

(32)利用公式(3)计算各个道集的纵波反射系数R_pp，所有道集的纵波反射系数R_pp形成反射系数序列；

(33)利用公式(4)计算误差函数E;

(34)判断是否满足终止条件，如果是，则停止迭代，转入步骤(35)，如果否，则将求取的误差函数E代入全局最优解Gbest，对每个粒子进行速度和位置更新，然后返回步骤(32)；

(35)输出对应的位置向量，即为所求的纵波速度v_p、横波速度v_s和密度ρ。

所述步骤(34)中的对每个粒子进行速度和位置更新是利用下面的公式实现的：

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

x_id(t+1)=x_id(t)+v_id(t+1)1≤i≤m1≤d≤n   (2)

其中，w为惯性权重因子，c₁,c₂为学习因子或加速因子，c₁为调节粒子飞向自身最好位置方向的步长，c₂为调节粒子飞向全局最好位置的步长，r₁，r₂为[0，1]内的随机数，t为当前迭代代数；第i粒子目前为止搜索到的最优位置为P_ibest，=(p_i1，p_i2，…p_im)，全局最优解G_best=(g₁，g₂,…，g_n)；

粒子群的初始位置和初始速度是随机产生的。

所述步骤(34)中的终止条件为：自适应函数E小于设定的截断误差，或者循环达到最大迭代次数。

与现有技术相比，本发明的有益效果是：

(1)能够降低对初始模型的依赖；

(2)能够提高反演的精度；

(3)具有较强的抗噪声能力。

附图说明

图1合成角度道集的地震记录。

图2-1纵波速度反演结果(圆点)与模型参数(黑色实线)对比图(不加噪声)。

图2-2横波速度反演结果(圆点)与模型参数(黑色实线)对比图(不加噪声)。

图2-3密度反演结果(圆点)与模型参数(黑色实线)对比图(不加噪声)。

图3添加20％噪声的合成角度道集的地震记录

图4-1纵波速度反演结果(圆点)与模型参数(黑色实线)对比图(加20％噪声)。

图4-2横波速度反演结果(圆点)与模型参数(黑色实线)对比图(加20％噪声)。

图4-3密度反演结果(圆点)与模型参数(黑色实线)对比图(加20％噪声)。

图5是实施例中的模型示意图。

图6本发明方法的步骤框图。

具体实施方式

下面结合附图对本发明作进一步详细描述：

粒子群算法***中，可以将每个优化问题的解想象成n维搜索空间中的一个点，称之为“粒子”(Partic l e)，整个***中多个粒子共存、合作寻优。每个粒子在搜索空间中以一定的速度飞行，这个速度根据其自身的飞行经验和同伴的飞行经验来动态调整。每个粒子追随当前的最优粒子而动，并经逐代搜索最后得到全局最优解。

粒子群算法的数学描述为：假设一个n维搜索空间，空间中总的粒子数为m，每个粒子所处的位置代表目标函数的一个潜在解，它们的位置向量x_i构成一个种群X=(x₁,x₂,…，x_m)，每个粒子在空间中移动速度为v_i。第i粒子目前为止搜索到的最优位置为P_ibest，=(ρ_i1,p_i2,…ρ_im)，整个种群目前为止搜索到的最优位置为G_best=(g₁，g₂,…，g_n)，则粒子的速度和位置可以按如下公式进行变化：

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

x_id(t+1)=x_id(t)+v_id(t+1)1≤i≤m1≤d≤n   (2)

式中，w称为惯性权重因子，c₁,c₂称为学习因子或加速因子，其中，c₁为调节粒子飞向自身最好位置方向的步长，c₂为调节粒子飞向全局最好位置的步长，r₁,r₂为[0，1]内的随机数，t为当前迭代代数。粒子群初始位置和初始速度随机产生，然后按照上面两式进行迭代，直至找到满意的解。

AVO技术以左普利兹方程为理论依据，1980年Aki和Richards给出了纵波反射振幅的近似表达式：

$R_{\mathrm{pp}}=\frac{1}{2}(1-4\frac{{\mathrm{\β}}^2}{{\mathrm{\α}}^2}{\sin }^2\mathrm{\θ})\frac{\mathrm{\Δ\ρ}}{\mathrm{\ρ}}+\frac{{\sec }^2\mathrm{\θ}}{2}\frac{\mathrm{\Δ\α}}{\mathrm{\α}}-4\frac{{\mathrm{\β}}^2}{{\mathrm{\α}}^2}{\sin }^2\mathrm{\θ}\frac{\mathrm{\Δ\β}}{\mathrm{\β}}---\left(3\right)$

其中，R_pp为纵波反射系数，α为纵波速度，β为横波速度，ρ为密度，θ为角度。

Δα=α₂-α₁，Δβ=β₂-β₁，Δρ=ρ₂-ρ₁

α=(α₁+α₂)／2，β=(β₁+β₂)／2

θ=(θ₁+θ₂)／2，ρ=(ρ₁+ρ₂)／2

下标1，2分别表示第一层介质和第二层介质。

都比较小，所有的角度都是实数。

因为公式(3)中α、β、ρ及△α、Δβ、Aρ这几个参数与纵波速度、横波速度以及密度有关，对公式(3)用粒子群算法进行迭代更新，可以得到纵波速度v_p、横波速度v_s以及密度ρ。

将这叠前道集数据作为输入进行反演。以求解计算出的R_pp与实际地震记录之间误差函数E最小作为目标函数。使用粒子群反演的方法寻找满足目标函数最小的纵波速度、横波速度以及密度。

$E=\left|\right|w\left(\mathrm{\θ}\right)*R_{\mathrm{pp}}(V_p,V_s,\mathrm{\ρ},\mathrm{\θ})-d\left(\mathrm{\θ}\right)\left|\right|\⇒\min ---\left(4\right)$

(4)式中的d(θ)是实际地震记录，R_pp(V_p,V_s,ρ,θ)就是利用公式(3)得到的。其中纵波速度v_p、横波速度v_s与公式(3)中的α为纵波速度，β为横波速度表示相同的含义，只不过公式(3)中α和β为地层界面分界处平均值，如第6页第一行公式所示。(4)表示地震子波w(θ)与反射系数的褶积(即卷积)，得到理论记录，该记录与实际地震记录的差值的范数(左右两条竖杠表示范数)，即使得两者的差值最小。

本发明的一个实施例如下：

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

首先，利用Aki和Richards方程对设计的一个一维模型(该模型的结构示意如图5所示。该模型为一维模型，含泥岩、含水砂岩、低含气饱和度砂岩和高含气饱和度砂岩组成)进行正演，每隔3°计算1°到30°范围内各地层的反射系数。计算出的10个反射系数序列分别与40Hz的零相位雷克子波进行褶积，得到的合成地震记录，即角道集数据，如图1所示。

将这十个(通用本方法时，对该数字无特殊规定，该实施例中为1-30度，每个3°计算出来的即为10个，一般角度道集在30°(小角度)计算较为准确)角道集数据作为输入进行反演。适应度函数E为实际地震记录与计算出的纵波反射系数R_pp之间的差值，以求解E最小为目标函数。如图6所示，本发明采用改进的粒子群算法AVO反演的步骤如下：

(1)将由模型正演得到的角道集记录(其他实施例可用其他记录作为输入，只要是待求解问题的函数即可，如此处角度道集是待求解的纵波速度、横波速度和密度的函数)作为观测值d(θ)(真值)。

(2)利用粒子群优化算法生成各角道集的反射系数序列(是利用公式(3)计算各个角度道集的Rpp，即反射系数序列，Rpp是α、β、ρ和θ的函数，对于给定的模型，α、β、ρ和θ的值是已知的。在给定角度θ的情况下，如果α、β、ρ的值变化，则Rpp也发生变化。后述的r(t)是纵波速度v_p，横波速度v_s，密度ρ的一个集合，即将纵波速度v_p，横波速度v_s，密度ρ简记为r(t)(需注意的是v_p与公式(3)里面的α，v_s与β含义相同。))如r(t)候选解。每个粒子的位置向量即为r(t)候选解，每个粒子的维数为3，即3个自变量纵波速度v_p，横波速度v_s，密度ρ。

(3)按照步骤(2)计算适应度函数E。将计算得到R_pp的与角道集数据比较，计算适应度函数E。(参考公式(4)E计算公式，其中d(θ)由步骤(1)得到，w(θ)给定，即1-30度，每隔3°给定，Rpp由步骤(2)得到)

(4)将求取的适应度函数E代入全局最优解Gbest(即整个种群目前为止搜索到的最优位置G_best(G_best=(g₁，_g2,…，g_n)，输入的道集是由多个点组成的，每次计算时只针对1个点进行计算，每个点的最优解对应g₁到g_n间的一个，故当满足终止条件得到的E，将被代入该点的全局最优解G_best))，并按照公式(1)和公式(2)对每个粒子进行位置更新；具体更新可参考公式(1)和公式(2)，粒子群算法中粒子的位置x_i对应所求问题潜在的解，此处对应r(t)(即纵波速度v_p，横波速度v_s，密度ρ)的值，r(t)的变化会引起R_pp的变化，粒子的速度v_i对应待求变量的变化值，即r(t)的变化值；每次迭代过程中，粒子的速度v_i发生变化，产生新的r(t)的变化值，也即会产生r(t)的变化，这样回到步骤(2)重新计算，则会引起R_pp的改变；步骤(3)中的适应度函数也会发生改变，则迭代可以继续)。

(5)利用改进的粒子群算法循环迭代，对最优解不断进行优化(主要的迭代步骤如前所述步骤(4)，)，直到满足终止条件(如自适应函数E小于某截断误差，或循环达到最大迭代次数(两者满足其一即可终止，细究的话，条件1为自然终止，条件2为强制终止。在假设不知道实际地层的纵波速度v_p，横波速度v_s，密度ρ的情况下，观测值d(θ)是从野外采集到的，假设地层的r(t)为某个值，通过比较其产生的记录(w(θ)与R_pp的褶积)与实际野外采集记录的差值，如果差值达到某个条件，则假设的r(t)是正确的，如果差值没满足条件，则需要对r(t)进行修改直到满足收敛条件。修改的规则即为步骤4所述。)，则此时个体最优与全局最优一致，所有粒子收敛于一点，对应的位置向量即为要求的纵波速度v_p，横波速度v_s，密度ρ。

在本例中利用改进的粒子群算法反演数值模型角道集数据的反射系数时，参数取值如下：误差函数小于1e-6，种群大小M=30，最小惯性因子w₁=0.1，最大惯性因子w₂=0.6，学习因子c₁=1.8，c₂=1.8，粒子位置取值的最大值为[3000,2000,2500]，最小值为[2000,1000,1500]，分别对应纵波速度v_p，横波速度v_s，密度ρ的限制范围。使用用本方法时，对这些数字的取值范围无具体规定，但须与实际情况符合，如纵波速度横波速度的取值范围，数量级在千。

反演出的弹性参数与模型的弹性参数的对比如图2-1至图2-3所示。可以看出，反演出的弹性参数(圆点)与真实的弹性参数(黑色实线)每一层都完全吻合。另外，为了试验本反演算法的抗噪性，对输入的角道集资料加入20％的噪声，如图3所示，然后按照步骤(1)到步骤(5)进行计算，反演得到的纵波速度v_p，横波速度v_s，密度ρ(圆点)与模型的弹性参数(黑色实线)对比如图4-1至图4-3所示。可以看出，在有噪声影响的情况下，反演结果与模型参数基本吻合，但数值上略有小的波动，但差异不大。在有噪声的情况下，反演得到的三个弹性参数，纵波速度受噪音的影响最小，横波速度次之，密度再次。认为在反演过程中，横波速度和密度的精确程度受到初值设定、岩石物理关系等因素的影响，稳定性不如纵波速度。为了使误差尽量减小，在可能的情况下，不仅使用纵波数据作为反演的观测值，还可以使用转换横波的信息。甚至在拥有VSP资料的情况下，可将反射纵波和反射横波的数据均作为反演的观测值，那么反演得到的弹性参数将会更加准确。

本发明属于石油勘探领域，涉及一种叠前资料的AVO(振幅相对偏移距变化)反演技术。AVO反演技术可充分利用地震资料得到其蕴含的地下地层的构造、岩性、物性等信息，但常规反演基于线性简化的模型，并且有较多的人为因素干扰，往往得不到高精度的储层岩性物性参数。本发明以粒子群优化算法的思想和原理为基础，将左普利茨(Zoeppritz)方程进行推导简化，使其适于粒子群优化算法的求解，然后将粒子群算法用于叠前反演，对理论模型进行叠前反演，取得令人满意的结果，并且具有较强的抗噪声能力。

上述技术方案只是本发明的一种实施方式，对于本领域内的技术人员而言，在本发明公开了应用方法和原理的基础上，很容易做出各种类型的改进或变形，而不仅限于本发明上述具体实施方式所描述的方法，因此前面描述的方式只是优选的，而并不具有限制性的意义。

Claims (6)

1.一种基于粒子群算法的AVO三参数反演方法，其特征在于：所述方法包括：

(1)输入叠前道集数据d(θ)；

(2)求取纵波反射系数R_PP;以所述纵波反射系数R_PP与子波w(θ)褶积同叠前道集数据d(θ)之间的误差函数E作为目标函数；

(3)使用粒子群反演的方法寻找到使所述目标函数的值为最小的纵波速度、横波速度以及密度，即为所求的三参数。

2.根据权利要求1所述的基于粒子群算法的AVO三参数反演方法，其特征在于：所述步骤(2)中是利用公式(3)求得纵波反射系数R_pp：

$R_{\mathrm{pp}}=\frac{1}{2}(1-4\frac{{\mathrm{\β}}^2}{{\mathrm{\α}}^2}{\sin }^2\mathrm{\θ})\frac{\mathrm{\Δ\ρ}}{\mathrm{\ρ}}+\frac{{\sec }^2\mathrm{\θ}}{2}\frac{\mathrm{\Δ\α}}{\mathrm{\α}}-4\frac{{\mathrm{\β}}^2}{{\mathrm{\α}}^2}{\sin }^2\mathrm{\θ}\frac{\mathrm{\Δ\β}}{\mathrm{\β}}---\left(3\right)$

其中，R_pp为纵波反射系数，α为纵波速度，β为横波速度，ρ为密度，θ为角度；

Δα=α₂-α₁，Δβ=β₂-β₁，Δρ=ρ₂-ρ₁

α=(α₁+α₂)／2，β=(β₁+β₂)／2

θ=(θ₁+θ₂)／2，ρ=(ρ₁+ρ₂)／2

下标1，2分别表示第一层介质和第二层介质。

3.根据权利要求2所述的基于粒子群算法的AVO三参数反演方法,其特征在于:所述步骤(2)中的误差函数E如下: $E=\left|\right|w\left(\mathrm{\θ}\right)*R_{\mathrm{pp}}(V_p,V_s,\mathrm{\ρ},\mathrm{\θ})-d\left(\mathrm{\θ}\right)\left|\right|\⇒\min ---\left(4\right)$

其中，R_pp(V_p，V_s，ρ，θ)·是利用公式(3)得到的；其中V_p为纵波速度，与公式(3)中的α含义相同，V_s为横波速度，与公式(3)中的β含义相同。

4.根据权利要求3所述的基于粒子群算法的AVO三参数反演方法，其特征在于：所述步骤(3)包括：

(31)将纵波速度v_p，横波速度v_s，密度ρ简记为r(t)，每个粒子的位置向量即为r(t)候选解，每个粒子的维数为3，即3个自变量：纵波速度vp、横波速度v_s和密度ρ；粒子的位置x_i对应r(t)的值，粒子的速度v_i对应r(t)的变化值；

(32)利用公式(3)计算各个道集的纵波反射系数R_pp，所有道集的纵波反射系数R_pp形成反射系数序列；

(33)利用公式(4)计算误差函数E;

(34)判断是否满足终止条件，如果是，则停止迭代，转入步骤(35)，如果否，则将求取的误差函数E代入全局最优解Gbest，对每个粒子进行速度和位置更新，然后返回步骤(32)；

(35)输出对应的位置向量，即为所求的纵波速度v_p、横波速度v_s和密度ρ。

5.根据权利要求4所述的基于粒子群算法的AVO三参数反演方法，其特征在于：所述步骤(34)中的对每个粒子进行速度和位置更新是利用下面的公式实现的：

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

x_id(t+1)=x_id(t)+v_id(+1)1≤i≤m1≤d≤n   (2)

其中，w为惯性权重因子，c₁,c₂为学习因子或加速因子，c₁为调节粒子飞向自身最好位置方向的步长，c₂为调节粒子飞向全局最好位置的步长，r₁，r₂为[0，1]内的随机数，t为当前迭代代数；第i粒子目前为止搜索到的最优位置为P_ibest，=(p_i1，p_i2,…p_im)，全局最优解G_best=(g₁，g₂,…，g_n)；

粒子群的初始位置和初始速度是随机产生的。

6.根据权利要求5所述的基于粒子群算法的AVO三参数反演方法，其特征在于：所述步骤(34)中的终止条件为：自适应函数E小于设定的截断误差，或者循环达到最大迭代次数。