CN109613608A - 一种任意复杂介质地震波传播矩阵模拟方法 - Google Patents

Abstract

本发明涉及一种任意复杂介质地震波传播矩阵模拟方法，该方法包括以下步骤：1)将水平地表的非均匀地质模型剖分成厚度为Δz的若干水平薄层；2)在每个水平薄层内对广义Lippmann‑Schwinger积分方程进行傅里叶变换，推导反射/透射矩阵；3)利用水平薄层之间的边界连续条件，通过反射/透射矩阵推导广义反射/透射矩阵；4)使用层间迭代关系，使用广义反射/透射矩阵在频率域下进行地震波场的求解；5)使用傅里叶反变换将地震波场转换到时间域，达成所求目标。本发明将传统广义反射/透射矩阵法(推广到界面起伏的分层非均匀地质结构模型，推导的反射/透射矩阵中含有层内非均匀介质的影响，可以在保证精度的同时大大节省计算资源。

Description

一种任意复杂介质地震波传播矩阵模拟方法

技术领域

本发明涉及一种任意复杂介质地震波传播透射与反射广义传播矩阵模拟方法，属于地震勘探技术领域。

背景技术

地球内部介质的复杂性和多样性制约地震勘探的发展。目前，地震波场数值模拟方法是地球物理常用的正演方法之一，该地震波场数值模拟方法是通过使用设定模型模拟地下介质，计算地震波场的正演结果，从而推测地下介质的分布。

沉积盆地表现为半空间下不规则界面的多层地质结构，每层内为随机非均质，地质模型的地震记录合成在近几十年来已经被广泛的研究和发展。目前较为成熟的地震记录合成方法包含有限元法、有限差分法、边界元-体积元法等全数值的模拟方法，对于这种界面起伏的分层非均质结构，需要进行较为精细的模型剖分，计算量与均匀模型相比将大大提升，急需发展一种半解析的高精度快速模拟技术。

发明内容

针对上述问题，本发明的目的是提供一种任意复杂介质地震波传播透射与反射广义传播矩阵模拟方法，该方法可以在保证精度的同时大大节省计算资源。

为实现上述目的，本发明采取以下技术方案：一种任意复杂介质地震波传播矩阵模拟方法，其特征在于，该方法包括以下步骤：1)将水平地表的非均匀地质模型剖分成厚度为Δz的若干水平薄层；2)在每个水平薄层内对广义Lippmann-Schwinger积分方程进行傅里叶变换，推导反射/透射矩阵；3)利用水平薄层之间的边界连续条件，通过反射/透射矩阵推导广义反射/透射矩阵；4)使用层间迭代关系，使用广义反射/透射矩阵在频率域下进行地震波场的求解；5)使用傅里叶反变换将地震波场转换到时间域，达成所求目标。

所述的任意复杂介质地震波传播矩阵模拟方法，优选的，在所述步骤2)中，具体包括如下过程：

①通过格林函数的频率域表达将广义Lippmann-Schwinger积分方程转化到频率-波数域，并使用Born近似，将体扰动的计算转化为边界扰动计算，根据广义Lippmann-Schwinger积分方程在频率-波数域的表达式，分别讨论观测点所在位置，第j水平薄层的位移与应力满足以下关系：

式中，为第j水平薄层的下行波的振幅；为第j+1水平薄层的上行波的振幅；E^(j)为第j水平薄层的厚度参数；E^(j+1)为第j+1水平薄层的厚度参数；和分别为观测点在第j水平薄层以下时上边界的位移与应力扰动参数；与分别为观测点在第j+1水平薄层以上时下边界的位移与应力扰动参数；为观测点在第j水平薄层以下时上边界的体扰动参数；为观测点在第j+1水平薄层以上时下边界的体扰动参数；α^(j)和β^(j)分别为第j水平薄层中的位移及其法相导数的参数；

同时，任意层中的地震波场可以解耦为上行波和下行波，有：

式中，和分别为观测点在第j水平薄层以上时上边界的位移与应力扰动参数；和为分别观测点在第j+1水平薄层以下时下边界的位移与应力扰动参数；为观测点在第j水平薄层以上时上边界的体扰动参数；为观测点在第j+1水平薄层以下时下边界的体扰动参数；为第j水平薄层的上行波的振幅；为第j+1水平薄层的下行波的振幅；

同时，层间的上下行波与水平薄层的反射透射系数矩阵满足以下关系：

式中，为反射/透射矩阵；为第j水平薄层下行波向上反射的反射系数；为第j水平薄层上行波向下反射的反射系数；为第j水平薄层上行波向上透射的透射系数；为第j水平薄层下行波向下透射的透射系数；

②将式(1)和式(2)联立，求得反射/透射矩阵的计算公式为：

③求解反射/透射矩阵，获得反射/透射矩阵的表达式：

式(5)-(8)即为式(4)的解，其中，I为单位矩阵。

所述的任意复杂介质地震波传播矩阵模拟方法，优选的，在所述步骤3)中，具体包括如下过程：

对比反射/透射矩阵和广义反射/透射矩阵与层间上下行波的关系，获得两个矩阵间的转化关系，广义反射/透射矩阵的计算公式如下：

对于震源之上的水平薄层，有：

对于震源之下的水平薄层，有：

式中，为自由表面上行波向下反射的广义反射系数；为第j水平薄层上行波向上透射的广义透射系数，为第j水平薄层上行波向下反射的广义反射系数；为第j水平薄层下行波向下透射的广义透射系数；为第j水平薄层下行波向上反射的广义反射系数。

所述的任意复杂介质地震波传播矩阵模拟方法，优选的，在所述步骤4)中，具体包括如下过程：

震源层的总上行波的振幅表达式如下：

类似地，震源层的总下行波的表达式为：

式中，s为震源所在层的编号；为第s-1层上行波向下反射的广义反射系数；为第s层上行波向下反射的广义反射系数；和分别为震源所在层的总上行波和总下行波的振幅；s_↑(z)和s_↓(z)为震源项的上行波和下行波的振幅；ξ^(s)为震源所在层下界面的深度；ξ^(s-1)为震源所在层上界面的深度；

则任意水平地表模型中的剖分的任意水平薄层的上行波及下行波的振幅按以下方式通过广义反射/透射矩阵进行迭代：

式中，为第s-1水平薄层上行波向上透射的广义透射系数；为第s水平薄层下行波向下透射的广义透射系数；

特别地，对于自由表面的位移的频率域表达式为：

式中u⁽⁰⁾(k_x)为自由表面的位移在频率域的表达形式；k_x为水平波数，为离散后的自由表面位移参数；i为虚数单位；x为水平位置；m为水平剖分节点编号；k_(x,m)为离散后的水平波数。

所述的任意复杂介质地震波传播矩阵模拟方法，优选的，在所述步骤5)中，达成的所求目标为：

式中，u⁽⁰⁾(x)为时间域下的地表位移；为x方向上的傅里叶反变换。

本发明由于采取以上技术方案，其具有以下优点：本发明将传统广义反射/透射矩阵法(只适用于界面起伏的分层均匀地质结构模型)推广到界面起伏的分层非均匀地质结构模型，推导的反射/透射矩阵中含有层内非均匀介质的影响。本发明可以在保证精度的同时大大节省计算资源。

附图说明

图1是本发明水平地表的非均匀地质模型剖分方式图；

图2是一个起伏地下界面的均匀地质模型图；

图3是一个起伏地下界面的非均匀地址模型图；

图4是本发明与边界元方法计算的均匀地质模型的地震合成记录比较图；

图5(a)-(c)是本发明与边界元方法计算的非均匀地质模型的地震合成记录比较图。

具体实施方式

下面结合附图和实施例对本发明进行详细的描述。然而应当理解，附图的提供仅为了更好地理解本发明，它们不应该理解成对本发明的限制。

本发明提出了一种任意复杂介质地震波传播矩阵模拟方法，包括以下步骤：

1)如图1所示，将水平地表的非均匀地质模型剖分成厚度为Δz的若干水平薄层；

2)在每个水平薄层内对广义Lippmann-Schwinger积分方程进行傅里叶变换，推导反射/透射矩阵，具体过程如下：

①通过格林函数的频率域表达将广义Lippmann-Schwinger积分方程转化到频率-波数域，并使用Born近似，将体扰动的计算转化为边界扰动计算，根据广义Lippmann-Schwinger积分方程在频率-波数域的表达式，分别讨论观测点所在位置，第j水平薄层的位移与应力满足以下关系：

式中，为第j水平薄层的下行波的振幅；为第j+1水平薄层的上行波的振幅；E^(j)为第j水平薄层的厚度参数；E^(j+1)为第j+1水平薄层的厚度参数；和分别为观测点在第j水平薄层以下时上边界的位移与应力扰动参数；与分别为观测点在第j+1水平薄层以上时下边界的位移与应力扰动参数；为观测点在第j水平薄层以下时上边界的体扰动参数；为观测点在第j+1水平薄层以上时下边界的体扰动参数；α^(j)和β^(j)分别为第j水平薄层中的位移及其法相导数的参数。

同时，任意层中的地震波场可以解耦为上行波和下行波，有：

式中，和分别为观测点在第j水平薄层以上时上边界的位移与应力扰动参数；和为分别观测点在第j+1水平薄层以下时下边界的位移与应力扰动参数；为观测点在第j水平薄层以上时上边界的体扰动参数；为观测点在第j+1水平薄层以下时下边界的体扰动参数；为第j水平薄层的上行波的振幅；为第j+1水平薄层的下行波的振幅。

同时，层间的上下行波与水平薄层的反射透射系数矩阵满足以下关系：

式中，为反射/透射矩阵；为第j水平薄层下行波向上反射的反射系数；为第j水平薄层上行波向下反射的反射系数；为第j水平薄层上行波向上透射的透射系数；为第j水平薄层下行波向下透射的透射系数。

②将式(1)和式(2)联立，可以求得反射/透射矩阵的计算公式为：

③求解反射/透射矩阵，可以获得反射/透射矩阵的表达式：

式(5)-(8)即为式(4)的解，其中，I为单位矩阵。

3)利用水平薄层之间的边界连续条件，通过反射/透射矩阵推导广义反射/透射矩阵：

对比反射/透射矩阵和广义反射/透射矩阵与层间上下行波的关系，获得两个矩阵间的转化关系，广义反射/透射矩阵的计算公式如下：

对于震源之上的水平薄层，有：

对于震源之下的水平薄层，有：

式(9)和(10)即为通过反射/矩阵推导广义反射/透射矩阵的计算公式。式中，为自由表面上行波向下反射的广义反射系数；为第j水平薄层上行波向上透射的广义透射系数，为第j水平薄层上行波向下反射的广义反射系数；为第j水平薄层下行波向下透射的广义透射系数；为第j水平薄层下行波向上反射的广义反射系数。

4)使用层间迭代关系，使用广义反射/透射矩阵在频率域下进行地震波场的求解；

震源层的总上行波的振幅表达式如下：

类似地，震源层的总下行波的表达式为：

式中，s为震源所在层的编号；为第s-1层上行波向下反射的广义反射系数；为第s层上行波向下反射的广义反射系数；和分别为震源所在层的总上行波和总下行波的振幅；s_↑(z)和s_↓(z)为震源项的上行波和下行波的振幅；ξ^(s)为震源所在层下界面(第s层)的深度；ξ^(s-1)为震源所在层上界面(第s-1层)的深度。

则任意水平地表模型中的剖分的任意水平薄层的上行波及下行波的振幅可以按以下方式通过广义反射/透射矩阵进行迭代：

式中，为第s-1水平薄层上行波向上透射的广义透射系数；为第s水平薄层下行波向下透射的广义透射系数。

特别地，对于自由表面的位移的频率域表达式为：

式中u⁽⁰⁾(k_x)为自由表面的位移在频率域的表达形式；k_x为水平波数，为离散后的自由表面位移参数；i为虚数单位；x为水平位置；m为水平剖分节点编号；k_(x,m)为离散后的水平波数。

5)使用傅里叶反变换将地震波场转换到时间域，达成所求目标：

式中，u⁽⁰⁾(x)为时间域下的地表位移；为x方向上的傅里叶反变换。

下面以具体实例来说明本发明的技术效果：

我们将本发明方法分别应用于一个起伏地下界面的均匀地质模型(如图2所示)及非均匀地质模型(如图3所示)，分为4层，速度分别为2700m/s，3200m/s，3700m/s与4500m/s，震源分布在模型正中600m深处，第二层与第三层的分界面起伏。

图4比较了本发明方法计算的地震合成记录结果与边界元方法计算的结果，两种方法有非常好的一致性，证明了方法的可靠性。

我们将上述模型中的第二层加入5％、15％和25％的速度扰动，其他参数不变，然后将本发明方法计算的地震合成记录结果与边界元方法计算的结果进行比较，比较图分别为图5(a)-(c)，两个结果仍有较好的一致性。

上述各实施例仅用于说明本发明，其中各部件的结构、连接方式和制作工艺等都是可以有所变化的，凡是在本发明技术方案的基础上进行的等同变换和改进，均不应排除在本发明的保护范围之外。

Claims (5)

1.一种任意复杂介质地震波传播矩阵模拟方法，其特征在于，该方法包括以下步骤：

1)将水平地表的非均匀地质模型剖分成厚度为Δz的若干水平薄层；

2)在每个水平薄层内对广义Lippmann-Schwinger积分方程进行傅里叶变换，推导反射/透射矩阵；

3)利用水平薄层之间的边界连续条件，通过反射/透射矩阵推导广义反射/透射矩阵；

4)使用层间迭代关系，使用广义反射/透射矩阵在频率域下进行地震波场的求解；

5)使用傅里叶反变换将地震波场转换到时间域，达成所求目标。

2.根据权利要求1所述的任意复杂介质地震波传播矩阵模拟方法，其特征在于，在所述步骤2)中，具体包括如下过程：

①通过格林函数的频率域表达将广义Lippmann-Schwinger积分方程转化到频率-波数域，并使用Born近似，将体扰动的计算转化为边界扰动计算，根据广义Lippmann-Schwinger积分方程在频率-波数域的表达式，分别讨论观测点所在位置，第j水平薄层的位移与应力满足以下关系：

式中，为第j水平薄层的下行波的振幅；为第j+1水平薄层的上行波的振幅；E^(j)为第j水平薄层的厚度参数；E^(j+1)为第j+1水平薄层的厚度参数；和分别为观测点在第j水平薄层以下时上边界的位移与应力扰动参数；与分别为观测点在第j+1水平薄层以上时下边界的位移与应力扰动参数；为观测点在第j水平薄层以下时上边界的体扰动参数；为观测点在第j+1水平薄层以上时下边界的体扰动参数；α^(j)和β^(j)分别为第j水平薄层中的位移及其法相导数的参数；

同时，任意层中的地震波场可以解耦为上行波和下行波，有：

式中，和分别为观测点在第j水平薄层以上时上边界的位移与应力扰动参数；和为分别观测点在第j+1水平薄层以下时下边界的位移与应力扰动参数；为观测点在第j水平薄层以上时上边界的体扰动参数；为观测点在第j+1水平薄层以下时下边界的体扰动参数；为第j水平薄层的上行波的振幅；为第j+1水平薄层的下行波的振幅；

同时，层间的上下行波与水平薄层的反射透射系数矩阵满足以下关系：

式中，为反射/透射矩阵；为第j水平薄层下行波向上反射的反射系数；为第j水平薄层上行波向下反射的反射系数；为第j水平薄层上行波向上透射的透射系数；为第j水平薄层下行波向下透射的透射系数；

②将式(1)和式(2)联立，求得反射/透射矩阵的计算公式为：

③求解反射/透射矩阵，获得反射/透射矩阵的表达式：

式(5)-(8)即为式(4)的解，其中，I为单位矩阵。

3.根据权利要求2所述的任意复杂介质地震波传播矩阵模拟方法，其特征在于，在所述步骤3)中，具体包括如下过程：

对比反射/透射矩阵和广义反射/透射矩阵与层间上下行波的关系，获得两个矩阵间的转化关系，广义反射/透射矩阵的计算公式如下：

对于震源之上的水平薄层，有：

对于震源之下的水平薄层，有：

式中，为自由表面上行波向下反射的广义反射系数；为第j水平薄层上行波向上透射的广义透射系数，为第j水平薄层上行波向下反射的广义反射系数；为第j水平薄层下行波向下透射的广义透射系数；为第j水平薄层下行波向上反射的广义反射系数。

4.根据权利要求3所述的任意复杂介质地震波传播矩阵模拟方法，其特征在于，在所述步骤4)中，具体包括如下过程：

震源层的总上行波的振幅表达式如下：

类似地，震源层的总下行波的表达式为：

式中，s为震源所在层的编号；为第s-1层上行波向下反射的广义反射系数；为第s层上行波向下反射的广义反射系数；和分别为震源所在层的总上行波和总下行波的振幅；s_↑(z)和s_↓(z)为震源项的上行波和下行波的振幅；ξ^(s)为震源所在层下界面的深度；ξ^(s-1)为震源所在层上界面的深度；

则任意水平地表模型中的剖分的任意水平薄层的上行波及下行波的振幅按以下方式通过广义反射/透射矩阵进行迭代：

式中，为第s-1水平薄层上行波向上透射的广义透射系数；为第s水平薄层下行波向下透射的广义透射系数；

特别地，对于自由表面的位移的频率域表达式为：

式中u⁽⁰⁾(k_x)为自由表面的位移在频率域的表达形式；k_x为水平波数，为离散后的自由表面位移参数；i为虚数单位；x为水平位置；m为水平剖分节点编号；k_(x，m)为离散后的水平波数。

5.根据权利要求4所述的任意复杂介质地震波传播矩阵模拟方法，其特征在于，在所述步骤5)中，达成的所求目标为：

式中，u⁽⁰⁾(x)为时间域下的地表位移；为x方向上的傅里叶反变换。