公开了一种基于横向约束的分步叠前弹性参数反演方法及系统，该方法及系统包括：基于经过同相轴匹配一致的PP波和PS波角道集，获得PP波和PS波角度域反射系数；利用Stewart反射系数公式进行多波联合反演，并且利用卡尔曼滤波算法对反演方程进行横向约束，获得纵波速度变化率和横波速度变化率；基于所述纵波速度变化率和横波速度变化率，代入Aki-Richards三参数反射系数公式，并同时利用卡尔曼滤波算法对反演方程进行横向约束，得到密度变化率；对获得的反演结果进行低频补偿获得最终的叠前弹性参数的绝对值。本发明通过基于横向约束的分步AVO反演方法，对地下储层识别及流体预测有重要意义。
1.一种基于横向约束的分步叠前弹性参数反演方法，其特征在于，包括：
1)基于经过同相轴匹配一致的PP波和PS波角道集，获得PP波和PS波角度域反射系数；
2)利用Stewart反射系数公式进行多波联合反演，并且利用卡尔曼滤波算法对反演方程进行横向约束，获得纵波速度变化率和横波速度变化率；
3)基于所述纵波速度变化率和横波速度变化率，代入Aki-Richards三参数反射系数公式，并同时利用卡尔曼滤波算法对反演方程进行横向约束，得到密度变化率；
4)对步骤2)和步骤3)中获得的反演结果进行低频补偿获得最终的叠前弹性参数的绝对值。
2.根据权利要求1所述的基于横向约束的分步叠前弹性参数反演方法，其中，在步骤2)中利用以下Stewart纵波反射系数公式和Stewart转换波反射系数公式获得纵波速度变化率 和横波速度变化率
其中，RPP、RPS分别为纵波、转换波反射系数；VP＝(VP1+VP2)/2，ΔVP＝VP2-VP1，VP1、VP2分别为界面上、下层的纵波速度；VS＝(VS1+VS2)/2，ΔVS＝VS2-VS1，VS1、VS2分别为界面上、下层的SV波速度；ρ＝(ρ1+ρ2)/2，Δρ＝ρ2-ρ1，ρ1、ρ2分别为界面上、下层的密度；θ＝(θ1+θ2)/2，θ1、θ2分别为P波入射角和透射角； 分别为转换波反射角和透射角。
3.根据权利要求1所述的基于横向约束的分步叠前弹性参数反演方法，其中，在步骤3)中通过如下公式得到密度变化率
其中，R'PP为新推导得到的反射系数，VP＝(VP1+VP2)/2，ΔVP＝VP2-VP1，VP1、VP2分别为界面上、下层的纵波速度；VS＝(VS1+VS2)/2，ΔVS＝VS2-VS1，VS1、VS2分别为界面上、下层的SV波速度，θ＝(θ1+θ2)/2，θ1、θ2分别为P波入射角和透射角。
4.根据权利要求1中所述的基于横向约束的分步叠前弹性参数反演方法，其中，在步骤
2-1)根据第k道反演的纵、横波速度变化率预测第k+1道的反演结果；
2-2)基于卡尔曼增益利用测量方程对步骤2-1)中预测的第k+1道反演结果进行修正进而获得最优化反演结果。
5.根据权利要求1中所述的基于横向约束的分步叠前弹性参数反演方法，其中，在步骤
3-1)根据第k道反演的密度变化率预测第k+1道的反演结果；
3-2)基于卡尔曼增益利用测量方程对步骤3-1)中预测的第k+1道反演结果进行修正进而获得最优化反演结果。
6.根据权利要求4或5中所述的基于横向约束的分步叠前弹性参数反演方法，其中，通过如下公式预测第k+1道的反演结果：
mk+1,k＝Ak,kmk,k+rk (4)其中，mk+1,k为由第k道预测的第k+1道的反演参数，Ak为状态变换矩阵，mk,k为由第k道反演的纵、横波速度变化率和密度变化率的最优结果，rk为mk,k和mk+1,k之间的状态转换误差。
7.根据权利要求4或5中所述的基于横向约束的分步叠前弹性参数反演方法，其中，所述测量方程为：
其中，dk、mk分别为AVO反演方程的地震记录与弹性参数，lk为地震数据反演结果的误差值，G为系数矩阵。
8.根据权利要求4或5中所述的基于横向约束的分步叠前弹性参数反演方法，其中，所述卡尔曼增益表达式为：
Zk+1＝Yk+1,kGk+1T(Gk+1Yk+1,kGk+1T+B)-1 (6)其中，Zk+1是第k+1道的卡尔曼增益，Yk+1,k为mk+1,k的预测误差协方差矩阵，B是lk的协方差矩阵，Gk+1为系数矩阵。
9.根据权利要求4或5中所述的基于横向约束的分步叠前弹性参数反演方法，其中，通过如下公式对获得最优化反演结果：
mk+1,k+1＝mk+1,k+Zk+1(dk+1-Gk+1mk+1,k) (7)其中，mk+1,k+1为第k+1道的最终优化反演结果，Z为第k+1道卡尔曼增益，Gk+1为系数矩阵，mk+1,k为由第k道预测的第k+1道的反演参数，dk+1是第k+1道地震数据。
10.一种基于横向约束的分步叠前弹性参数反演系统，其上存储有计算机程序，其中，所述程序被处理器执行时实现以下步骤：
步骤1：基于经过同相轴匹配一致的PP波和PS波角道集，获得PP波和PS波角度域反射系数；
步骤2：利用Stewart反射系数公式进行多波联合反演，并且利用卡尔曼滤波算法对反演方程进行横向约束，获得纵波速度变化率和横波速度变化率；
步骤3：基于所述纵波速度变化率和横波速度变化率，代入Aki-Richards三参数反射系数公式，并同时利用卡尔曼滤波算法对反演方程进行横向约束，得到密度变化率；
步骤4：对步骤2和步骤3中获得的反演结果进行低频补偿获得最终的叠前弹性参数的绝对值。
一种基于横向约束的分步叠前弹性参数反演方法及系统
技术领域
[0001] 本发明涉及地震资料AVO叠前弹性参数反演与流体识别领域，更具体地，涉及一种基于横向约束的分步叠前弹性参数反演方法及系统。
背景技术
[0002] 叠前地震反演是储层预测和流体识别的主要技术之一，纵、横波速度和密度等弹性参数可以研究储层岩性和孔隙流体性质，其中由反演得到的密度参数是检测储层含油气
性的重要参数。常规的纵波叠前地震反演方法仅使用了反射纵波地震数据，由纵波叠前地
震数据提取纵、横波速度和密度参数存在多解性、反演过程不适定等问题，因此人们试图寻找一种稳定的高精度叠前地震反演方法。
[0003] AVO反演作为一种叠前反演方法已经成为复杂油气藏储层预测的关键技术，弹性参数反演结果的准确性是AVO叠前反演技术的关键难题。部分实际地震资料信噪比、分辨率不高，PP波反演结果不稳定、反演过程不适定等问题亟需解决。因此，有必要开发一种基于横向约束的分步叠前弹性参数反演方法及系统。
[0004] 公开于本发明背景技术部分的信息仅仅旨在加深对本发明的一般背景技术的理解，而不应当被视为承认或以任何形式暗示该信息构成已为本领域技术人员所公知的现有
技术。
发明内容
[0005] 本发明提出了一种基于横向约束的分步叠前弹性参数反演方法及系统，其能够通过对密度反演结果的准确性问题，提出了基于卡尔曼滤波约束下的分步多波AVO反演方法
研究，解决部分实际地震资料信噪比、分辨率不高，PP波反演结果不稳定、反演过程不适定等问题。
[0006] 根据本发明的一方面，提出了一种基于横向约束的分步叠前弹性参数反演方法。所述方法可以包括：
[0007] 1)基于经过同相轴匹配一致的PP波和PS波角道集，获得PP波和PS波角度域反射系数；
[0008] 2)利用Stewart反射系数公式进行多波联合反演，并且利用卡尔曼滤波算法对反演方程进行横向约束，获得纵波速度变化率和横波速度变化率；
[0009] 3)基于所述纵波速度变化率和横波速度变化率，代入Aki-Richards三参数反射系数公式，并同时利用卡尔曼滤波算法对反演方程进行横向约束，得到密度变化率；
[0010] 4)对步骤2)和步骤3)中获得的反演结果进行低频补偿获得最终的叠前弹性参数的绝对值。
[0011] 优选地，在步骤2)中利用以下Stewart纵波反射系数公式和Stewart转换波反射系数公式获得纵波速度变化率 和横波速度变化率
[0012]
[0013]
[0014] 其中，RPP、RPS分别为纵波、转换波反射系数；VP＝(VP1+VP2)/2，ΔVP＝VP2-VP1，VP1、VP2分别为界面上、下层的纵波速度；VS＝(VS1+VS2)/2，ΔVS＝VS2-VS1，VS1、VS2分别为界面上、下层的SV波速度；ρ＝(ρ1+ρ2)/2，Δρ＝ρ2-ρ1，ρ1、ρ2分别为界面上、下层的密度；θ＝(θ1+θ2)/2，θ1、θ2分别为P波入射角和透射角； 分别为转换波反射角和透射角。
[0015] 优选地，在步骤3)中通过如下公式得到密度变化率
[0016]
[0017] 其中，R'PP为新推导得到的反射系数，VP＝(VP1+VP2)/2，ΔVP＝VP2-VP1，VP1、VP2分别为界面上、下层的纵波速度；VS＝(VS1+VS2)/2，ΔVS＝VS2-VS1，VS1、VS2分别为界面上、下层的SV波速度，θ＝(θ1+θ2)/2，θ1、θ2分别为P波入射角和透射角。
[0018] 优选地，在步骤2)中，通过如下步骤进行横向约束：
[0019] 2-1)根据第k道反演的纵、横波速度变化率预测第k+1道的反演结果；
[0020] 2-2)基于卡尔曼增益利用测量方程对步骤3-1)中预测的第k+1道反演结果进行修正进而获得最优化反演结果。
[0021] 优选地，在步骤3)中，通过如下步骤进行横向约束：
[0022] 3-1)根据第k道反演的密度变化率预测第k+1道的反演结果；
[0023] 3-2)基于卡尔曼增益利用测量方程对步骤3-1)中预测的第k+1道反演结果进行修正进而获得最优化反演结果。
[0024] 优选地，通过如下公式预测第k+1道的反演结果：
[0025] mk+1,k＝Ak,kmk,k+rk (4)
[0026] 其中，mk+1,k为由第k道预测的第k+1道的反演参数，Ak为状态变换矩阵，mk,k为由第k道反演的纵、横波速度变化率和密度变化率的最优结果，rk为mk,k和mk+1,k之间的状态转换误差。
[0027] 优选地，所述测量方程为：
[0028] dk＝Gmk+lk (5)
[0029] 其中，dk、mk分别为AVO反演方程的地震记录与弹性参数，lk为地震数据反演结果的误差值，G为系数矩阵。
[0030] 优选地，所述卡尔曼增益表达式为：
[0031] Zk+1＝Yk+1,kGk+1T(Gk+1Yk+1,kGk+1T+B)-1 (6)
[0032] 其中，Zk+1是第k+1道的卡尔曼增益，Yk+1,k为mk+1,k的预测误差协方差矩阵，B是lk的协方差矩阵，Gk+1为系数矩阵。
[0033] 优选地，通过如下公式对获得最优化反演结果：
[0034] mk+1,k+1＝mk+1,k+Zk+1(dk+1-Gk+1mk+1,k) (7)
[0035] 其中，mk+1,k+1为第k+1道的最终优化反演结果，Z为第k+1道卡尔曼增益，Gk+1为系数矩阵，mk+1,k为由第k道预测的第k+1道的反演参数，dk+1是第k+1道地震数据。
[0036] 根据本发明的另一方面，提出了一种基于横向约束的分步叠前弹性参数反演系统，其上存储有计算机程序，其中，所述程序被处理器执行时实现以下步骤：
[0037] 步骤1：基于经过同相轴匹配一致的PP波和PS波角道集，获得PP波和PS波角度域反射系数；
[0038] 步骤2：利用Stewart反射系数公式进行多波联合反演，并且利用卡尔曼滤波算法对反演方程进行横向约束，获得纵波速度变化率和横波速度变化率；
[0039] 步骤3：基于所述纵波速度变化率和横波速度变化率，代入Aki-Richards三参数反射系数公式，并同时利用卡尔曼滤波算法对反演方程进行横向约束，得到密度变化率；
[0040] 步骤4：对步骤2和步骤3中获得的反演结果进行低频补偿获得最终的叠前弹性参数的绝对值。
[0041] 本发明具有其它的特性和优点，这些特性和优点从并入本文中的附图和随后的具体实施方式中将是显而易见的，或者将在并入本文中的附图和随后的具体实施方式中进行
详细陈述，这些附图和具体实施方式共同用于解释本发明的特定原理。
附图说明
[0042] 通过结合附图对本发明示例性实施例进行更详细的描述，本发明的上述以及其它目的、特征和优势将变得更加明显，其中，在本发明示例性实施例中，相同的参考标号通常代表相同部件。
[0043] 图1示出了根据本发明的一种基于横向约束的分步叠前弹性参数反演方法的步骤的流程图；
[0044] 图2(a)示出了深度道集CDP 1000点的纵波角道集；
[0045] 图2(b)示出了深度道集CDP 1000点的转换横波角道集；
[0046] 图3(a)示出了AVO反演后经过低频补偿的M测线目标层纵波速度剖面；
[0047] 图3(b)示出了AVO反演后经过低频补偿的M测线目标层横波速度剖面；
[0048] 图3(c)示出了AVO反演后经过低频补偿的M测线目标层密度剖面。
具体实施方式
[0049] 下面将参照附图更详细地描述本发明。虽然附图中显示了本发明的优选实施例，然而应该理解，可以以各种形式实现本发明而不应被这里阐述的实施例所限制。相反，提供这些实施例是为了使本发明更加透彻和完整，并且能够将本发明的范围完整地传达给本领
域的技术人员。
[0050] 图1示出了根据本发明的一种基于横向约束的分步叠前弹性参数反演方法的步骤的流程图。
[0051] 在该实施例中，根据本发明的一种基于横向约束的分步叠前弹性参数反演方法的可以包括：
[0052] 步骤101，基于经过同相轴匹配一致的PP波和PS波角道集，获得PP波和PS波角度域反射系数；
[0053] 步骤102，利用Stewart反射系数公式进行多波联合反演，并且利用卡尔曼滤波算法对反演方程进行横向约束，获得纵波速度变化率和横波速度变化率；
[0054] 在一个示例中，利用以下Stewart纵波反射系数公式和Stewart转换波反射系数公式获得纵波速度变化率 和横波速度变化率
[0055]
[0056]
[0057] 其中，RPP、RPS分别为纵波、转换波反射系数；VP＝(VP1+VP2)/2，ΔVP＝VP2-VP1，VP1、VP2分别为界面上、下层的纵波速度；VS＝(VS1+VS2)/2，ΔVS＝VS2-VS1，VS1、VS2分别为界面上、下层的SV波速度；ρ＝(ρ1+ρ2)/2，Δρ＝ρ2-ρ1，ρ1、ρ2分别为界面上、下层的密度；θ＝(θ1+θ2)/2，θ1、θ2分别为P波入射角和透射角； 分别为转换波反射角和透射角。
[0058] 在一个示例中，通过如下步骤进行横向约束：
[0059] 2-1)根据第k道反演的纵、横波速度变化率预测第k+1道的反演结果；
[0060] 2-2)基于卡尔曼增益利用测量方程对步骤2-1)中预测的第k+1道反演结果进行修正进而获得最优化反演结果。
[0061] 具体地，在步骤2-1)中通过如下公式预测第k+1道的反演结果：
[0062] mk+1,k＝Ak,kmk,k+rk (4)
[0063] 其中，mk+1,k为由第k道预测的第k+1道的反演参数，Ak为状态变换矩阵，mk,k为由第k道反演的纵、横波速度变化率最优结果，rk为mk,k和mk+1,k之间的状态转换误差。
[0064] 步骤2-2)中所述测量方程为：
[0065] dk＝Gmk+lk (5)
[0066] 其中，dk、mk分别为AVO反演方程的地震记录与弹性参数，lk为地震数据反演结果的误差值，G为公式(1)、(2)、(3)式中方程的系数矩阵。
[0067] 在步骤2-2)所述卡尔曼增益表达式为：
[0068] Zk+1＝Yk+1,kGk+1T(Gk+1Yk+1,kGk+1T+B)-1 (6)
[0069] 其中，Zk+1是第k+1道的卡尔曼增益，Yk+1,k为mk+1,k的预测误差协方差矩阵，B是lk的协方差矩阵，Gk+1为(1)，(2)，(3)式中方程的系数矩阵。
[0070] 在步骤2-2)中通过如下公式对获得最优化反演结果：
[0071] mk+1,k+1＝mk+1,k+Zk+1(dk+1-Gk+1mk+1,k) (7)
[0072] 其中，mk+1,k+1为第k+1道的最终优化反演结果，Zk+1为第k+1道卡尔曼增益，Gk+1为系数矩阵，mk+1,k为由第k道预测的第k+1道的反演参数，dk+1是第k+1道地震数据。
[0073] 步骤103，基于所述纵波速度变化率和横波速度变化率，代入
[0074] Aki-Richards三参数反射系数公式，并同时利用卡尔曼滤波算法对反演方程进行横向约束，得到密度变化率；
[0075] 在一个示例中，通过如下公式得到密度变化率
[0076]
[0077] 其中，R'PP为推导得到的新反射系数方程，VP＝(VP1+VP2)/2，ΔVP＝VP2-VP1，VP1、VP2分别为界面上、下层的纵波速度；VS＝(VS1+VS2)/2，ΔVS＝VS2-VS1，VS1、VS2分别为界面上、下层的SV波速度，θ＝(θ1+θ2)/2，θ1、θ2分别为P波入射角和透射角。
[0078] 在一个示例中，通过如下步骤进行横向约束：
[0079] 3-1)根据第k道反演的密度变化率预测第k+1道的反演结果；
[0080] 3-2)基于卡尔曼增益利用测量方程对步骤3-1)中预测的第k+1道反演结果进行修正进而获得最优化反演结果。
[0081] 具体地，在步骤3-1)中通过如下公式预测第k+1道的反演结果：
[0082] mk+1,k＝Ak,kmk,k+rk (4)
[0083] 其中，mk+1,k为由第k道预测的第k+1道的反演参数，Ak为状态变换矩阵，mk,k为由第k道反演的密度变化率最优结果，rk为mk,k和mk+1,k之间的状态转换误差。
[0084] 在步骤3-2)中所述测量方程为：
[0085] dk＝Gmk+lk (5)
[0086] 其中，dk、mk分别为AVO反演方程的地震记录与弹性参数，lk为地震数据反演结果的误差值，G为公式(1)、(2)、(3)式中方程的系数矩阵。
[0087] 在步骤3-2)所述卡尔曼增益表达式为：
[0088] Zk+1＝Yk+1,kGk+1T(Gk+1Yk+1,kGk+1T+B)-1 (6)
[0089] 其中，Zk+1是第k+1道的卡尔曼增益，Yk+1,k为mk+1,k的预测误差协方差矩阵，B是lk的协方差矩阵，Gk+1为公式(1)、(2)、(3)式中方程的系数矩阵。
[0090] 在步骤3-2)中通过如下公式对获得最优化反演结果：
[0091] mk+1,k+1＝mk+1,k+Zk+1(dk+1-Gk+1mk+1,k) (7)
[0092] 其中，mk+1,k+1为第k+1道的最终优化反演结果，Z为第k+1道卡尔曼增益，Gk+1为系数矩阵，mk+1,k为由第k道预测的第k+1道的反演参数，dk+1是第k+1道地震数据。
[0093] 步骤104，对步骤102和步骤103中获得的反演结果进行低频补偿获得最终的叠前弹性参数的绝对值。
[0094] 在一个示例中，低频补偿是通过对反演获得的弹性参数结果的反射率计算即可获取参数的绝对值，这个计算的过程将低频成分包含进来。
[0095] 应用示例
[0096] 为便于理解本发明实施例的方案及其效果，以下给出一个具体应用示例。本领域技术人员应理解，该示例仅为了便于理解本发明，其任何具体细节并非意在以任何方式限
制本发明。
[0097] 本发明选取了中国LJ地区A油田M测线的纵波和转换横波地震数据进行测试。该区在地质上位于凹陷南缓斜坡，受燕山运动和喜马拉雅运动的影响，形成了向北倾没的大型
鼻状构造，油气藏类型主要为受构造控制的断块油气藏，兼有岩性、泥岩裂缝等特殊油气
藏。如图2(a)所示，为原始地震记录的深度道集CDP 1000点的纵波角道集，图2(b)为深度道集CDP 1000点的转换波角道集，采用的地震资料PP波主频接近35Hz，PS波主频接近15Hz。在时间域内将纵波地震数据与转换横波地震数据匹配一致得到纵波角道集和转换横波角道
集。
[0098] 经过本发明中的基于横向约束的基于横向约束的分步叠前弹性参数反演方法反演后获得如图3(a)中所示的经过低频补偿的M测线目标层段纵波剖面，如图3(b)中所示的
经过低频补偿的M测线目标层段横波剖面和如图3(c)中所示的经过低频补偿的M测线目标
层段密度剖面。根据图中可以清晰的看出反演的目标层，经与生产井信息进行对比发现，纵波速度，即图3(a)中黑色线段标注的生产井段内的油层和密度反演结果，即图3(c)中黑色
线段标注的生产井段内的油层接近真实值范围，实现了较为准确的预测。
[0099] 本发明还提出了一种基于横向约束的分步叠前弹性参数反演系统，其上存储有计算机程序，其中，所述程序被处理器执行时实现以下步骤：
[0100] 步骤1：基于经过同相轴匹配一致的PP波和PS波角道集，获得PP波和PS波角度域反射系数；
[0101] 步骤2：利用Stewart反射系数公式进行多波联合反演，并且利用卡尔曼滤波算法对反演方程进行横向约束，获得纵波速度变化率和横波速度变化率；
[0102] 步骤3：基于所述纵波速度变化率和横波速度变化率，代入Aki-Richards三参数反射系数公式，并同时利用卡尔曼滤波算法对反演方程进行横向约束，得到密度变化率；
[0103] 步骤4：对步骤2和步骤3中获得的反演结果进行低频补偿获得最终的叠前弹性参数的绝对值。
[0104] 综上所述，本发明通过针对部分实际地震资料信噪比、分辨率不高，PP波反演结果不稳定、反演过程不适定等问题，本发明提出了一种基于横向约束的分步叠前弹性参数反演方法及系统，能够较好的反演地下储层中密度等弹性参数，其基于卡尔曼滤波算法则增
强横向连续性，理论严谨可靠且操作流程简单实用，为AVO反演识别提供了帮助。
[0105] 本领域技术人员应理解，上面对本发明的实施例的描述的目的仅为了示例性地说明本发明的实施例的有益效果，并不意在将本发明的实施例限制于所给出的任何示例。
[0106] 以上已经描述了本发明的各实施例，上述说明是示例性的，并非穷尽性的，并且也不限于所披露的各实施例。在不偏离所说明的各实施例的范围和精神的情况下，对于本技术领域的普通技术人员来说许多修改和变更都是显而易见的。
