叠前混合非线性反演方法及计算机可读存储介质与流程
本发明涉及油气地球物理技术领域，更具体地，涉及一种叠前混合非线性反演方法。
背景技术：
自1919年描述各向同性介质界面下纵波入射时反射和透射系数的zoeppritz方程诞生以来，上世纪60年代以后，许多学者依据实际勘探地震需要将其纵波入射时的纵波反射(p-p)和横波反射(p-sv)公式近似成多种解析形式，目前大约已有三十多种，形成了p-p波反演、p-sv波反演和p-p与p-sv波联合反演等三类反演方法，被广泛地应用于实际储层和流体预测之中。
叠前反演通常是指avo反演，其狭义理论基础为zoeppritz方程p-p部分或其近似式和褶积模型。时域褶积模型理论基本观点是地下某点反射振幅是其反射系数与地震子波的褶积，而反射系数是联系地下弹性参数的桥梁，反演即据振幅等参数反演地下的弹性参数信息。依据平面简谐波传播理论，当纵波(p波)传至弹性界面时根据应力与应变的关系和连续条件，在满足snell定律的条件下可以获得不同弹性介质的反演和透射系数。尽管叠前反演是利用多次覆盖资料求取弹性信息的必经之路，但实际资料的叠前反演却是严重的不适定问题，具有很强的多解性。不同方法或流程对于同一资料反演的结果常差别甚远，而同一方法与流程反演的不同弹性参数精度也很不相同，难于满足勘探、开发与生产的需求。
由于zoeppritz方程完整形式过于复杂且物理意义不太明确，诸多学者从区分固液相态与突出泊松比、体现速度和密度相对变化、幂级数或射线参数、突出弹性模量等四大类进行了约十三种形式的近似，比如最为典型的aki&richards近似等，从而实现了弹性参数在实际资料情况下的反演。依据经典的近似关系，叠前的反演方程为非线性方程组，主要包括纵波反射系数、横波反射系数、密度反射梯度和横纵波速度比等四个未知量。忽略地震子波、地质模型、噪音干扰等因素，依据对弹性未知量的求解策略不同，现有的叠前反演方法大致可以分为叠前线性反演和叠前非线性反演两种(zhangetal.,2013)。
叠前线性反演方法常常假设横纵波速度比为已知的常数，从而将非线性方程组线性化求解，直接求解包含纵波反射系数、横波反射系数和密度梯度的线性近似方程，这种方法操作简便，运算量小，是工业界最为常用的方法，然而由于其直接线性化的近似导致非线性量的反演精度较低，且其结果强烈依赖于初始模型和子波估算质量(sun,1999)。非线性反演方法则直接采用非线性寻优方案对纵波反射系数、横波反射系数、密度梯度和横纵波速度比进行直接的四参量非线性求解。叠前非线性的求解算法大致可以分为从目标函数出发(比如最速下降法，牛顿法，共轭梯度法，非线性规划法，非线性最小二乘法等)和从自变量出发(比如模拟退火法，原子跃迁法，量子退火法，量子遗传法，人工神经元法，遗传算法，蚁群算法，粒子群算法，免疫算法等)两类(张远银，2015)。得益于计算机的快速进步，非线性算法可以直接被用于求解地球物理问题，然而，叠前非线性反演方法的运算量巨大，且实际地质问题常常存在多个目标函数极值从而导致运算陷入局部寻优(zhangetal.,2013)。因此，有必要开发一种具有较高的反演精度和反演速度的叠前混合非线性反演方法。
公开于本发明背景技术部分的信息仅仅旨在加深对本发明的一般背景技术的理解，而不应当被视为承认或以任何形式暗示该信息构成已为本领域技术人员所公知的现有技术。
技术实现要素：
本发明提出了一种叠前混合非线性反演方法，其能够将传统四参数非线性问题降解为两参数非线性问题，实现高效率、高精度地求取各种弹性参数，提高了运算效率与精度。
根据本发明的一方面，提出了一种叠前混合非线性反演方法。所述方法可以包括：对于叠前地震数据进行叠前线性反演，获得横波反射系数rs的初始值与横纵波速度比rk的初始值；基于所述横波反射系数rs、所述横纵波速度比rk、密度梯度rd，构造组合弹性参量ks；将叠前纵波反射系数公式表示为关于所述密度梯度rd、纵波反射系数rp和组合弹性参量ks的线性方程，对所述叠前纵波反射系数公式进行线性求解，获得所述密度梯度rd、纵波反射系数rp和所述组合弹性参量ks；构造关于所述横波反射系数rs和横纵波速度比rk的目标函数，基于所述横波反射系数rs的初始值与横纵波速度比rk的初始值对所述目标函数进行非线性求解，进而获得最终的横波反射系数与最终的横纵波速度比。
优选地，所述叠前纵波反射系数公式为：
其中，rpp(θ)表示叠前纵波反射系数，θ表示反射角与入射角的平均值。
优选地，所述组合弹性参量ks为：
优选地，所述目标函数为：
其中，f(rs,rk)表示目标函数。
优选地，所述目标函数为：
其中，f(rs,rk)表示目标函数。
优选地，通过标准粒子群算法对所述目标函数进行非线性求解。
根据本发明的另一方面，提出了一种计算机可读存储介质，其上存储有计算机程序，其中，所述程序被处理器执行时实现以下步骤：对于叠前地震数据进行叠前线性反演，获得横波反射系数rs的初始值与横纵波速度比rk的初始值；基于所述横波反射系数rs、所述横纵波速度比rk、密度梯度rd，构造组合弹性参量ks；将叠前纵波反射系数公式表示为关于所述密度梯度rd、纵波反射系数rp和组合弹性参量ks的线性方程，对所述叠前纵波反射系数公式进行线性求解，获得所述密度梯度rd、纵波反射系数rp和所述组合弹性参量ks；构造关于所述横波反射系数rs和横纵波速度比rk的目标函数，基于所述横波反射系数rs的初始值与横纵波速度比rk的初始值对所述目标函数进行非线性求解，进而获得最终的横波反射系数与最终的横纵波速度比。
优选地，所述叠前纵波反射系数公式为：
其中，rpp(θ)表示叠前纵波反射系数，θ表示反射角与入射角的平均值。
优选地，所述组合弹性参量ks为：
优选地，所述目标函数为：
其中，f(rs,rk)表示目标函数。
本发明的有益效果在于：将反演方程中的线性量与非线性量分离，采用线性反演的方法快速准确反演线性量，降低方程的非线性化程度；采用粒子群非线性反演方法反演非线性量，利用线性反演的结果约束或指导非线性反演过程，提高反演的精度与速度，提高弹性参数反演精度与效率。
本发明的方法具有其它的特性和优点，这些特性和优点从并入本文中的附图和随后的具体实施方式中将是显而易见的，或者将在并入本文中的附图和随后的具体实施方式中进行详细陈述，这些附图和具体实施方式共同用于解释本发明的特定原理。
附图说明
通过结合附图对本发明示例性实施例进行更详细的描述，本发明的上述以及其它目的、特征和优势将变得更加明显，其中，在本发明示例性实施例中，相同的参考标号通常代表相同部件。
图1示出了根据本发明的叠前混合非线性反演方法的步骤的流程图。
图2a、图2b分别示出了根据本发明的一个实施例的目标函数的误差值与部分放大的目标函数的误差值的示意图。
图3a、图3b分别示出了根据本发明的一个实施例的横波反射系数的误差值与横纵波速度比的误差值的示意图。
图4a、图4b、图4c、图4d分别示出了根据本发明的一个实施例的纵波反射系数、横波反射系数、密度梯度、横纵波速度比与理论值的对比示意图。
图5a、图5b、图5c分别示出了根据叠前线性反演、叠前非线性反演与本发明的一个实施例的纵波反射系数的示意图，图5d示出了三者的误差对比示意图。
图6a、图6b、图6c分别示出了根据叠前线性反演、叠前非线性反演与本发明的一个实施例的横波反射系数的示意图，图6d示出了三者的误差对比示意图。
图7a、图7b、图7c分别示出了根据叠前线性反演、叠前非线性反演与本发明的一个实施例的密度梯度的示意图，图7d示出了三者的误差对比示意图。
图8a、图8b、图8c分别示出了根据叠前线性反演、叠前非线性反演与本发明的一个实施例的横纵波速度比的示意图，图8d示出了三者的误差对比示意图。
图9示出了根据叠前线性反演、叠前非线性反演与本发明的一个实施例的计算时间的对比示意图。
具体实施方式
下面将参照附图更详细地描述本发明。虽然附图中显示了本发明的优选实施例，然而应该理解，可以以各种形式实现本发明而不应被这里阐述的实施例所限制。相反，提供这些实施例是为了使本发明更加透彻和完整，并且能够将本发明的范围完整地传达给本领域的技术人员。
图1示出了根据本发明的叠前混合非线性反演方法的步骤的流程图。
根据本发明的叠前混合非线性反演方法可以包括：
步骤101，对于叠前地震数据进行叠前线性反演，获得横波反射系数rs的初始值与横纵波速度比rk的初始值。
具体地，叠前反演的基本原理如下，其基于gildow(1987)对aki&richards公式(1980)的重新整理与简化：
其中，rpp(θ)表示某一界面的叠前纵波反射系数，θ表示反射角与入射角的平均值，
在非线性方程组的实际计算过程中，公式中的rp，rs，rk和rd均为未知数据，rpp(θ)和θ为已知数据。针对公式(1)进行线性反演则将rk与rs直接以常数化近似，从而可以把非线性方程组化简为线性方程组求解，即可根据公式(1)通过线性反演方法求取横波反射系数rs的初始值与横纵波速度比rk的初始值。
步骤102，基于横波反射系数rs、横纵波速度比rk、密度梯度rd，构造组合弹性参量ks；在一个示例中，组合弹性参量ks为：
步骤103，将叠前纵波反射系数公式表示为关于密度梯度rd、纵波反射系数rp和组合弹性参量ks的线性方程，对叠前纵波反射系数公式进行线性求解，获得密度梯度rd、纵波反射系数rp和组合弹性参量ks。
实际上，非线性方程组的求解过程并不一定绝对地线性或者非线性，其往往可以是多个线性量与分线性量的组合，这也意味着叠前反演可以依照未知量的属性特征分解为线性和非线性求解的组合，从而提高运算效率。将公式(1)的线性量与非线性量进行分离，可得到叠前纵波反射系数公式为：
其中，rpp(θ)表示叠前纵波反射系数。
具体地，将横波反射系数rs、横纵波速度比rk、密度梯度rd代入公式(2)，构造组合弹性参量ks。将公式(1)表示为关于密度梯度rd、纵波反射系数rp和组合弹性参量ks的线性方程即为公式(3)，并进行线性求解，获得密度梯度rd、纵波反射系数rp和组合弹性参量ks。
步骤104，构造关于横波反射系数rs和横纵波速度比rk的目标函数，基于横波反射系数rs的初始值与横纵波速度比rk的初始值对目标函数进行非线性求解，进而获得最终的横波反射系数与最终的横纵波速度比。
在一个示例中，目标函数为：
其中，f(rs,rk)表示目标函数。
在一个示例中，目标函数为：
其中，f(rs,rk)表示目标函数。
具体地，构造关于横波反射系数rs和横纵波速度比rk的目标函数，其中目标函数可以为公式(4)或公式(5)，基于横波反射系数rs的初始值与横纵波速度比rk的初始值对目标函数进行非线性求解，获得最终的横波反射系数与最终的横纵波速度比。
在一个示例中，通过标准粒子群算法对目标函数进行非线性求解。
具体地，采用的非线性算法为标准粒子群算法，该算法是对鸟类群体的防御、捕食等行为仿真，将问题的搜索空间类比于鸟类的飞行空间，算法中每个优化问题的可行解都是搜索空间中的一只鸟，称之为“粒子”。每个粒子的运动可以用上述几条规则来描述，根据自己和其它粒子的“飞行经验”寻优，从而达到全空间搜索最优解的目的。假设一个n维搜索空间，每一维空间中总的粒子数为m，每个粒子所处的位置代表目标函数的一个潜在解，它们的位置量xi构成一个种群x＝(x1,x2,…xm)，每个粒子在空间中移动速度为vi。第i粒子目前为止搜索到的最优位置为pibest＝(pi1,pi2,…pim)，整个种群目前为止搜索到的最优位置为gbest＝(g1,g2,…gn)，则粒子的速度和位置可以按如下公式进行变化：
vid(t+1)＝wvid(t)+c1r1(pid(t)-xid(t))+c2r2(gd(t)-xid(t))(6)
xid(t+1)＝xid(t)+vid(t+1)1≤i≤m1≤d≤n(7)
其中，w表示惯性权重因子，cl为调节粒子飞向自身最好位置方向的步长，c2为调节粒子飞向全局最好位置的步长，r1,r2为[0，1]内的随机数，t为当前迭代代数。为避免由于速度过大而使粒子飞过了最优解的位置，有必要对速度的最大值加以限制，限定速度上限值为vmax，当vid>vmax时，令vid＝vmax；当vid<-vmax时，令vid＝-vmax。粒子群初始位置和初始速度随机产生，然后按照公式(6)与公式(7)进行迭代，直至找到满意的解或者达到最大迭代次数。
本方法将传统四参数非线性问题降解为两参数非线性问题，实现高效率、高精度地求取各种弹性参数，提高了运算效率与精度。
应用示例
为便于理解本发明实施例的方案及其效果，以下给出一个具体应用示例。本领域技术人员应理解，该示例仅为了便于理解本发明，其任何具体细节并非意在以任何方式限制本发明。
根据本发明的叠前混合非线性反演方法包括以下步骤：
对于叠前地震数据进行叠前线性反演，获得横波反射系数rs的初始值与横纵波速度比rk的初始值；
将横波反射系数rs、横纵波速度比rk、密度梯度rd代入公式(2)，构造组合弹性参量ks；
将公式(1)表示为关于密度梯度rd、纵波反射系数rp和组合弹性参量ks的线性方程即为公式(3)，并进行线性求解，获得密度梯度rd、纵波反射系数rp和组合弹性参量ks；
构造关于横波反射系数rs和横纵波速度比rk的目标函数为公式(4)，基于横波反射系数rs的初始值与横纵波速度比rk的初始值，通过标准粒子群算法对目标函数进行非线性求解，获得最终的横波反射系数与最终的横纵波速度比。
图2a、图2b分别示出了根据本发明的一个实施例的目标函数的误差值与部分放大的目标函数的误差值的示意图，从图中可以看出：随着迭代次数的增加，目标函数的误差值逐渐减小，基本上从第15次开始趋于稳定，稳定后的误差值很小，几乎可被忽略。
图3a、图3b分别示出了根据本发明的一个实施例的横波反射系数的误差值与横纵波速度比的误差值的示意图。虽然由于非线性量在求解过程中陷入局部寻优的概率不同，横波反射系数rs和横纵波速度比rk的误差变化规律有所差别，但整体处于很小的数值范围。相比rs而言，rk的非线性程度更高，但其获得了相对较小数值的误差程度，这是因为在本发明的相对误差的计算方式下，横纵波速度比相对反射系数较大的数值会对相对降低误差的统计额度。
图4a、图4b、图4c、图4d分别示出了根据本发明的一个实施例的纵波反射系数、横波反射系数、密度梯度、横纵波速度比与理论值的对比示意图。
进一步对比本发明的实施例所得的四项参数最终结果与理论值的差别，如图4a、图4b、图4c、图4d所示，纵波反射系数rp和密度梯度rd的误差几乎为零，而横波反射系数rs和横纵波速度比rk的误差也十分小，非常准确地描述了理论值的趋势和细节。
图5a、图5b、图5c分别示出了根据叠前线性反演、叠前非线性反演与本发明的一个实施例的纵波反射系数的示意图，图5d示出了三者的误差对比示意图。
图6a、图6b、图6c分别示出了根据叠前线性反演、叠前非线性反演与本发明的一个实施例的横波反射系数的示意图，图6d示出了三者的误差对比示意图。
图7a、图7b、图7c分别示出了根据叠前线性反演、叠前非线性反演与本发明的一个实施例的密度梯度的示意图，图7d示出了三者的误差对比示意图。
图8a、图8b、图8c分别示出了根据叠前线性反演、叠前非线性反演与本发明的一个实施例的横纵波速度比的示意图，图8d示出了三者的误差对比示意图。
叠前线性反演和本发明的叠前混合非线性反演方法均取得了非常准确的rp和rd结果，因为这两种方法在求解过程中均采用线性办法求解纵波反射系数，其结果为完全独立的线性量。另外，对于非线性量rs和横纵波速度比rk的各种方法求解结果均存在不同程度的误差，叠前非线性反演的误差程度最大。一方面这是因为横波参量的参数可解性本身较低，其在非线性方程组式中又是与横纵波速度比密切相关的非线性变量，因而求解难度较大；但更主要的原因还是由于四参量的非线性寻优路径或策略严重降低了单一弹性参量的求解精度。由于对横纵波速度比直接为0.5的初始运算近似，叠前线性反演的rs和横纵波速度比结果在幅值上与理论值有所偏差；而本发明的叠前混合非线性反演方法只非线性求解两个参量，且其初始值与约束条件以相对准确的迭代反演结果为基础，因而最终结果品质最高。
图9示出了根据叠前线性反演、叠前非线性反演与本发明的一个实施例的运算时间的对比示意图。叠前非线性反演方法的运算时间最多，为233.35秒，叠前线性反演耗时最少，为0.134秒。本发明的叠前混合非线性反演方法耗时4.25秒，相比于其他两类方法的精度和时间，其工作效率最高。
综上所述，本方法将传统四参数非线性问题降解为两参数非线性问题，实现高效率、高精度地求取各种弹性参数，提高了运算效率与精度。
根据本发明提出了一种计算机可读存储介质，其上存储有计算机程序，其中，程序被处理器执行时实现以下步骤：对于叠前地震数据进行叠前线性反演，获得横波反射系数rs的初始值与横纵波速度比rk的初始值；基于横波反射系数rs、横纵波速度比rk、密度梯度rd，构造组合弹性参量ks；将叠前纵波反射系数公式表示为关于密度梯度rd、纵波反射系数rp和组合弹性参量ks的线性方程，对叠前纵波反射系数公式进行线性求解，获得密度梯度rd、纵波反射系数rp和组合弹性参量ks；构造关于横波反射系数rs和横纵波速度比rk的目标函数，基于横波反射系数rs的初始值与横纵波速度比rk的初始值对目标函数进行非线性求解，进而获得最终的横波反射系数与最终的横纵波速度比。
在一个示例中，叠前纵波反射系数公式为：
其中，rpp(θ)表示叠前纵波反射系数，θ表示反射角与入射角的平均值。
在一个示例中，组合弹性参量ks为：
在一个示例中，目标函数为：
其中，f(rs,rk)表示目标函数。
在一个示例中，目标函数为：
其中，f(rs,rk)表示目标函数。
在一个示例中，通过标准粒子群算法对目标函数进行非线性求解。
以上已经描述了本发明的各实施例，上述说明是示例性的，并非穷尽性的，并且也不限于所披露的各实施例。在不偏离所说明的各实施例的范围和精神的情况下，对于本技术领域的普通技术人员来说许多修改和变更都是显而易见的。本文中所用术语的选择，旨在最好地解释各实施例的原理、实际应用或对市场中的技术的改进，或者使本技术领域的其它普通技术人员能理解本文披露的各实施例。
- 还没有人留言评论。精彩留言会获得点赞！