任意各向异性介质多方位反射P波走时计算方法与流程

本发明涉及石油地球物理勘探，特别是涉及到一种任意各向异性介质多方位反射p波走时计算方法。

背景技术：

1、在地震勘探领域，反射波走时或时差近似方法被广泛应用于反射波数据处理的多个分支中，其中，双曲时差近似公式被广泛应用于各向同性介质反射数据处理方法之中。然而，在各向异性介质中，时差往往是非双曲线型的，这就需要我们基于物理学近似和各向异性参数研究新的走时近似公式。

2、各向异性介质中，应用最广泛的介质是横向各向同性介质，而随着地震采集、处理技术的不断发展，更多更加复杂、对称性较差的各向异性介质越来越受到关注，如正交各向异性介质、单斜各向异性介质等。在传统的各向异性模型构建过程中，需要首先对介质各向异性的类型做出假设，并对对称方向进行估计，这往往是一个复杂的过程。同时，模型构建的过程往往还基于另一个假设，即介质的一个或几个对称轴或对称面始终与坐标轴平行或垂直。在具有明显应力构造的区域或者多组不同方向裂隙的高裂隙带这类复杂各向异性地层中，上述假设的应用会遇到问题。同时，上述模型构建过程还依赖先验信息的可靠性，如果各向异性介质的对称类型、对称方向等先验信息时错误的，那么模拟的结果也必然时错误的。psenc ik(2017)提出了基于弱各向异性(weak ani sotrop ic，wa)参数的任意各向异性介质二维反射p波走时近似方法，该方法消除了对各向异性类型的假设以及估计对称轴/面的方向的依赖，但还是存在着纵横波走时不解耦的问题，且仅针对二维局部坐标系下的反射p波走时进行了研究。

3、在申请号：cn202110340994.1的中国专利申请中，涉及到一种考虑层间多次波的弹性反射波走时反演方法，包括以下步骤：输入反射地震记录以及建立的线性增加的线性初始速度模型、开展最小二乘真振幅偏移获取弹性波偏移剖面、进行震源侧非线性born正演和检波点侧非线性borm正演获取含多次波的背向扰动波场、建立时移互相关目标函数，并利用模型参数梯度计算公式构建反演参数梯度、将所得反演参数梯度对线性初始速度模型进行迭代更新等步骤。该方法给出了弹性多次反射波走时反演技术框架，并对结构复杂的marmous i模型进行了测试，多次反射波反演结果能够大致刻画marmous i模型轮廓，以此开展的常规fwi结果更加稳健。

4、在申请号：cn201611099983.4的中国专利申请中，涉及到一种基于p波到时与波形振幅的震源定位方法，包括以下步骤：步骤1：通过微震监测系统获取震源参数：步骤2：在未知p波波速v的情况下，求解一次震源坐标(x0,y0,z0)及相应速度v；步骤3：计算波形振幅趋势线的斜率步骤4：以logistic函数拟合波形振幅趋势线的斜率k与每个震源的p波波速v之间的函数关系式；步骤5：计算修正后的各个传感器对应的p波波速步骤6：利用修正后的各个传感器对应的p波波速求解二次震源坐标(x,y,z)。该发明方法具有适用性强、准确性高等特点。

5、在申请号：cn202211185705.6的中国专利申请中，涉及到一种地震面波走时和重力异常联合反演方法，涉及地震数据处理技术领域。该技术方案利用双差层析成像方法来获取高分辨率的研究区地壳三维p波、s波速度结构及地震重定位结果，然后应用体波走时与重力联合反演的方法进一步改善地壳三维p波速度和密度结构，根据反演所得的地壳p波、s波速度、密度结构及地震重定位结果，根据修正后的线性化反演模型建立联合反演目标函数；利用反演目标函数的解对待反演的地震面波进行反演成像。该技术方案利用地震面波数据和高分辨率卫星重力数据对地震面波进行联合反演成像，不仅可以大大提高计算的效率，而且也可以提高反演结果的可靠性，提高了结构成像的横向分辨率。

6、以上现有技术均与本发明有较大区别，未能解决我们想要解决的技术问题，为此我们发明了一种新的任意各向异性介质多方位反射p波走时计算方法。

技术实现思路

1、本发明的目的是提供一种解决了对介质各向异性对称类型、对称方向等先验信息的依赖问题的任意各向异性介质多方位反射p波走时计算方法。

2、本发明的目的可通过如下技术措施来实现：任意各向异性介质多方位反射p波走时计算方法，该任意各向异性介质多方位反射p波走时计算方法包括：

3、步骤1、建立局部坐标系和全局坐标系；

4、步骤2、进行局部坐标系下反射p波走时计算；

5、步骤3、进行局部坐标系下各向异性参数和全局坐标系下各向异性参数之间的转换；

6、步骤4、进行全局坐标系下不同方位反射p波走时的计算。

7、本发明的目的还可通过如下技术措施来实现：

8、在步骤1，假设一笛卡尔坐标系作为全局坐标系，其x3轴是垂直方向的，x1，x2轴是水平方向且相互正交；背景介质为一均匀的地层，地层底部的水平反射面与(x1，x2)平面平行；在这种情况下，p波的入射下行和反射上行射线位于垂直于反射面的平面内，并且相对于反射面法线对称。

9、在步骤1，针对地震波观测系统，定义局部坐标系；地震波从全局坐标系的原点处的震源s激发，在点o处反射，传播至检波器r；将点s和r所在的轴定义为局部坐标系的x′1轴，与x′1轴在同一水平面且与x′1轴垂直的轴为x′2轴，垂直于x′1轴、x′2轴的轴为x′3轴；局部坐标系的原点与震源s一致，并且还表示全局坐标系xi的原点，垂直轴x3在局部和全局两个坐标系中也是相同的；反射射线位于垂直平面(x′1，x′3)内。

10、在步骤2，通过局部坐标系下的a参数εx、εz、ηy以及偏移距x、反射面深度h、参考纵波速度α计算反射波走时。

11、步骤2包括：

12、步骤21，建立反射走时基本计算公式；

13、步骤22，运用时差分析中的表示法改写公式；

14、步骤23，根据弱各向异性近似求射线矢量；

15、步骤24，根据相速度平方近似求射线速度；

16、步骤25，线性化处理简化反射走时近似公式。

17、在步骤21，在局部坐标系中，反射波走时t可通过以下公式进行计算：

18、

19、其中，x为偏移距，即震源s和检波器r之间的距离，h是水平反射面的深度，与地层的对称平面重合；t＝t(x)表示反射p波走时v＝v(n)表示射线速度，射线速度v(n)在入射射线和反射射线中的大小相同。

20、在步骤22，根据时差分析中广泛应用的表示法，公式(1)可以改写为：

21、

22、其中，为层厚度归一化偏移距，α为参考各向同性介质中的p波相速度，t0为零偏移距双程走时；运用公式(2)转化后，公式(1)变为：

23、

24、其中公式(3)中的相位矢量n可能会偏离反射射线所在的垂直平面(x′1，x′3)。

25、在步骤23，引入弱各向异性近似，在弱各向异性介质条件下，射线矢量与相位矢量的方向相同且位于平面(x′1，x′3)内，射线矢量n在局部坐标系中通过表示的分量为：

26、

27、其中n′3表达式中的+表示入射方向，-表示反射方向；

28、根据上行波和下行波之间的区别，可将公式(3)分为两部分：沿参考射线下行方向的走时td和沿上行方向的走时tu；

29、

30、式中，nd和nu分别代表参考射线的上行方向和下行方向；根据t2＝(td+tu)2，可以得到：

31、

32、为求公式(6)中的反射波走时，需要得到射线速度v(n)的表达式。

33、在步骤24，引入相速度平方近似，即通过相同方向相速度平方的一阶近似c2(n)来近似表示射线速度的平方v2(n)；局部坐标系下的p波射线速度公式：

34、

35、公式(7)中，α是预设的参考速度，α的值不影响相速度的值。εx，εz分别代表介质在坐标轴x1、x3方向p波相速度与参考速度α的差异；ηy代表介质在坐标轴x2方向p波速度的变化梯度；χy代表在平面(x1，x3)内p波偏离轴向的程度；ξ15，ξ35代表在平面(x1，x3)内p波方向与快、慢s波方向的偏离程度。

36、根据公式(7)可得，沿参考射线下行方向的射线速度和由上行方向的走时分别为：

37、

38、其中，pd和pu分别为：

39、

40、

41、将公式(8)代入公式(6)中，即可得到反射走时的近似公式：

42、

43、在步骤25，对公式(11)进行进一步线性化处理，在忽略公式(11)中的a参数χy，ξ15和ξ35的影响之后，公式(11)可以改写为：

44、

45、其中，

46、

47、在步骤3，基于bond变换和任意各向异性理论的不同坐标系下各向异性参数直接转换，包括坐标系转换和a参数转换。

48、在步骤3，坐标系转换矩阵为：

49、

50、其中，为局部坐标系的走向角，θ为局部坐标系的倾向角。

51、在步骤3，将全局坐标系中的a参数使用上标′来与局部坐标系中的a参数进行区分：

52、

53、

54、

55、公式(17)中，εx，εy，εz分别代表介质在坐标轴x1、x2、x3方向p波相速度与参考速度α的差异；ηx，ηy，ηz代表介质在坐标轴x1、x2、x3方向p波速度的变化梯度；χx，χy，χz代表在平面(x2，x3)，(x1，x3)，(x1，x2)内p波偏离轴向的程度；ξ15，ξ35代表在平面(x1，x3)内p波方向与快、慢s波方向的偏离程度；ξ24，ξ34代表在平面(x2，x3)内p波方向与快、慢s波方向的偏离程度；ξ16，ξ26代表在平面(x1，x2)内p波方向与快、慢s波方向的偏离程度。rij、dij为坐标转换因子由局部坐标系的走向角与局部坐标系的倾向角θ的三角函数组成。rij公式见权利要求12，dij公式如下：

56、d11＝r22r33+r23r32，d12＝r21r33+r23r31，d13＝r22r31+r21r32，

57、d21＝r12r33+r13r32，d22＝r11r33+r13r31，d23＝r12r31+r11r32，

58、d31＝r12r23+r13r22，d32＝r11r23+r13r21，d33＝r12r21+r11r22。

59、在步骤4，通过全坐标系下的a参数εx、εz、ηy以及偏移距x、反射面深度h、参考纵波速度α和测线方位角计算反射波走时。

60、步骤4包括：

61、步骤41，确立转换矩阵；

62、步骤42，将局部坐标系下a参数转换为全局坐标系下a参数；

63、步骤43，建立多方位任意各向异性介质反射p波走时计算公式。

64、在步骤41，假设在平面上沿不同方向分布的地震测线，测线的走向角为以测线走向为x轴方向可以构建局部坐标系；当此坐标系向全局坐标系转换时，其转换矩阵r为：

65、

66、在步骤42，局部坐标系中的a参数ε′x，η′y和ε′z可以通过以下公式转换为全局坐标系：

67、

68、

69、ε′z＝εz    (21)。

70、在步骤43，将公式(19)-(21)代入公式(13)中，可得全局坐标系下分方位反射p波走时计算公式：

71、

72、

73、公式(19)-(22)中，εx，εy，εz分别代表介质在坐标轴x1、x2、x3方向p波相速度与参考速度α的差异；ηx，ηy，ηz代表介质在坐标轴x1、x2、x3方向p波速度的变化梯度；χz代表在平面(x1,x2)内p波偏离轴向的程度；ξ16，ξ26代表在平面(x1,x2)内p波方向与快、慢s波方向的偏离程度。

74、本发明中的任意各向异性介质多方位反射p波走时计算方法，利用任意各向异性理论及a参数计算任意各向异性介质多方位反射p波走时。经过弱各向异性近似、相速度平方近似等近似方法，简化了局部坐标系下反射p波的计算公式。利用bond变换推导a参数转换公式，计算全局坐标系下分方位反射p波走时。

75、本发明的优点是解决了对介质各向异性对称类型、对称方向等先验信息的依赖问题；而且在适用于最复杂的任意各向异性介质的同时，也能够向下兼容各向异性较为简单的介质。并且在保障精度的同时，简化后的公式能够大大提高计算效率，并为反演方法打下基础。地球物理数据处理人员可依据本方法模拟复杂各向异性介质反射p波走时，或依据本发明提供的正演方法开发反演方法。