基于虚数域Rytov近似的裂缝型地层地震散射波场特征模拟方法与流程

本发明涉及一种地下裂缝地震散射波场特征正演模拟方法，具体是关于一种基于虚数域rytov近似的裂缝型地层地震散射波场特征模拟方法，属于勘探地球物理学领域。

背景技术：

随着地球物理勘探技术的不断发展，勘探目标由常规油气藏转向非常规油气藏，对于地震波正演模拟的要求也越来越高。地下均质地层受构造等运动的影响，易发育定向排列的高角度缝，形成裂缝型地层，且对于非常规勘探领域，由于地层中的裂缝具有油气运移与储集的作用，致使裂缝型储层成为勘探重点目标之一。因此，具有定向排列特征的高角度裂缝介质的正演模拟在新时期油气藏开采具有重要意义。积分法是地震波波动方程正演模拟的主要方法，微分法是利用差分方法求导算得波动方程差分形式解，在推导过程中采用高阶近似。尽管可以保证波场传播的所有信息，但是对于裂缝这种小尺度的非均质体，存在稳定性和频散等问题，影响其计算的稳定与精度。

格林(green)函数积分正演模拟方法可以有效的模拟和刻画地震波在这种小尺度非均匀介质中的波长特征，受到广泛关注。aki提出尾波是介质非均匀性引起的背向散射波，martin和flatte建立了适用于地震波在任何介质中传播的相屏法，wu,r.s利用不同方法得到了格林函数的散射波场和弹性波的背景散射场，积分方程正演模拟理论逐渐成熟。秦雪霏等和蕾蕾在面向地震勘探应用的散射波场正演模拟研究做了深入的研究。但如何在保证精度的前提下，提高green函数积分法的效率有待深入研究。

moczo首次提出可变网格的思路，对于局部变化较大的区域，改变相应的空间和时间步长，以便在满足精度的前提下提高效率。tessmer在波动方程有限差分正演模拟中引入了变网格步长的思想，使得刻画的地下介质的地震记录更加精确。为了达到模拟孔洞缝介质的精度要求，孙林洁提出多级变网格思想，从而使得变网格步长从十倍提高致百倍以上。然而，目前可变网格主要应用于微分模拟方法。

技术实现要素：

针对上述问题，本发明的目的是提供一种基于虚数域rytov近似的裂缝型地层散射波场特征模拟方法，该方法可模拟地下裂缝引起的地震散射波场特征，指导地震识别裂缝型储层，并且该方法可以较好地解决模拟效率和稳定性的问题，且对于提高模拟结果信噪比和精度有较高的优势。

为实现上述目的，本发明采取以下技术方案：一种虚数域rytov近似的裂缝型地层散射波场特征模拟方法，包括如下步骤：

步骤1：建立地面观测系统以确定地面震源激发点位置与接收位置，设计合适的震源信号，作为正演模拟系统的初始波场；

步骤2：对实数域震源信号进行二维傅里叶变换，将实数域震源信号转换至虚数域震源信号，并将虚数域震源信号作为地下裂缝介质的入射波场ui＝u0(x,y,z；ω)，其中u0为背景波场，x、y和z为笛卡尔坐标系的坐标轴，ω为角频率；

步骤3：基于rytov近似建立散射场外推方程：

①根据rytov变换公式得：

②基于rytov近似，得到虚数域rytov近似散射场方程us为：

式中，i为虚数单位；s0为均匀背景慢度；kx、ky和kz分别为x、y和z方向的波数，且其中k0为背景波数；δz为沿垂直方向将地下裂缝介质划分为等厚的薄板厚度；δh为薄板间的水平慢度扰动量；

③基于虚数域rytov近似散射场方程，前向散射场和背向散射场分别表征为：

式中，ft^-1为二维傅里叶反变换算子。

步骤4：首先利用虚数域rytov近似散射波传播方程(3)和(4)分别计算上层地下裂缝介质的背向散射场与前向散射场并对上层地下裂缝介质的背向散射场进行存储；然后基于rytov近似，将上层地下裂缝介质的背景波场u0(z1)与前向散射场代入将计算结果作为下一层地下裂缝介质的入射波场，将新的入射波场代入虚数域rytov近似散射波传播方程(3)和(4)，分别计算下一层地下裂缝介质的前向散射场和背向散射场并且对下一层地下裂缝介质的背向散射场继续进行存储；

步骤5：根据步骤3和步骤4计算虚数域粗网格rytov近似散射波场，当模拟的波场传播到地下裂缝介质的加密边界时，利用插值原理将水平方向上的粗网格波场进行插值加密处理，直至满足网格加密倍数，垂直方向上对粗网格进行分割，并基于波场外推理论，对加密网格中地下裂缝介质的波场逐层进行外推计算；

步骤6：当对加密网格中地下裂缝介质的波场外推计算至地下裂缝介质的加密边界时，基于转换算法，将水平方向上加密网格的波场变换至粗网格中，垂直方向上则继续基于波场外推理论对加密边界下粗网格中地下裂缝介质的波场逐层进行外推计算；

步骤7：当波场传播直到地下裂缝介质的底部时，地下裂缝介质的波场外推计算完毕；此时，从地下裂缝介质的底部开始，利用地下裂缝介质底部计算出的背向散射场作为上一层地下裂缝介质的入射波场，利用式(4)计算向后传播的背向散射场

步骤8：将向后传播的背向散射场与步骤4存储下来的n-1层地下裂缝介质的背向散射场进行叠加处理，将叠加结果作为上一层地下裂缝介质的入射波场，利用式(3)计算向后传播的前向散射场

步骤9：重复步骤4～8，波场向后传播，当到达地下裂缝介质的接收位置时，地下裂缝介质的波场的反向外推计算完毕；最后利用二维傅里叶反变换将虚数域地震散射记录转换至实数域地震散射记录，对时间-空间域的波场特征进行表征，将反变换后的波场结果进行存储，最终获得接收地下裂缝介质地震散射记录。

所述的裂缝型地层散射波场特征模拟方法，优选地，在所述步骤5中，波场由粗网格插值得到加密网格的算法为：

波场由加密网格转换到粗网格的算法为：

式中，k为水平方向网格加密的倍数；m为粗网格点次序；j为加密网格点次序；n为粗网格点数；un^m为粗网格中第m个点的波场值；ukn^km为加密网格波场；ukn^km+j为km+j点在加密网格中的波场；ukn^k(m+1)为k(m+1)点在加密网格中的波场；为l点在加密网格中的波场。

本发明由于采取以上技术方案，其具有以下优点：1、本发明基于波场弱变化的rytov近似，虚数域波场计算更加稳定，避免born近似计算波场时的累计误差。2、本发明在保证精度和稳定性的前提下，以构建的虚数域rytov近似散射场方程为基础，利用多级可变网格格林函数积分法对裂缝地层的地震散射响应特征进行模拟，其精度可达到毫米尺度。3、本发明对于相同尺度的裂缝模型，多级可变网格格林函数积分法要比多级可变网格有限差分法的模拟效率高数十倍，而且对计算机内存消耗要低很多。4、本发明不受直达波和多次波等影响，裂缝地震散射响应可以清晰的分辨出来。

附图说明

为了更清楚的说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单的介绍，下面描述中的附图仅仅是本发明的一些实施例。

图1是含裂缝地层的网格剖分示意图；

图2是本发明的流程示意图；

图3是用于正演模拟的含裂缝地层速度模型示意图；

图4(a)-(f)是不同裂缝倾角情况下的地面接收记录图；

图5(a)-(d)是裂缝在不同位置处时的地面接收记录图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明的基本思路是首先通过虚数域rytov近似散射场方程，得到背景波场与散射波场的解析式，在此基础上构建格林函数积分算子；然后利用薄板近似和小角度近似条件推导三维地震波向前传播的总波场和背向散射场相屏传播算子；最后采用变网格对地下裂缝介质进行剖分(如图1所示)，确定加密边界，对地下裂缝介质中含有裂缝的区域采用加密网格格林函数积分进行正演模拟，对不含裂缝的区域采用粗网格格林函数积分进行正演模拟。

如图2所示，本发明提供的虚数域rytov近似的裂缝型地层散射波场特征模拟方法，其具体步骤如下：

步骤1：建立地面观测系统以确定地面震源激发点位置与接收位置，设计合适的震源信号，作为正演模拟系统的初始波场。

步骤2：对实数域震源信号进行二维傅里叶变换(即对实数域震源信号的频率与波数进行表征)，将实数域震源信号转换至虚数域震源信号，并将虚数域震源信号作为地下裂缝介质的入射波场ui＝u0(x,y,z；ω)，其中u0为背景波场，x、y和z为笛卡尔坐标系的坐标轴，ω为角频率；根据波场的叠加原理，地下裂缝介质的总波场u可看作入射波场ui与散射波场us的叠加：u＝ui+us。

步骤3：散射波场us可等效为前向散射场与背向散射场的线性叠加：以地下裂缝介质的背景速度作为传播速度，基于散射波传播方程可对前向散射场和背向散射场进行定量表征。在波场外推计算时，传统方法多采用born近似，但由于该近似是基于弱散射近似建立定量表征，因此计算过程中必然产生累计误差。针对上述问题，本发明引入rytov近似建立散射场外推方程。地下裂缝会引起相位变化，rytov近似基于散射场的弱变化近似，为了提高计算的稳定性，通过以下方法建立散射场外推方程：

①根据rytov变换公式得：

②基于rytov近似，可以得到虚数域rytov近似散射场方程us为：

式中，i为虚数单位；s0为均匀背景慢度；kx、ky和kz分别为x、y和z方向的波数，且其中k0为背景波数；δz为沿垂直方向将地下裂缝介质划分为等厚的薄板厚度；δh为薄板间的水平慢度扰动量；

③基于虚数域rytov近似散射场方程，前向散射场和背向散射场可分别表征为：

式中，ft^-1为二维傅里叶反变换算子。

步骤4：首先利用虚数域rytov近似散射波传播方程(3)和(4)分别计算上层地下裂缝介质的背向散射场与前向散射场并对上层地下裂缝介质的背向散射场进行存储；然后基于rytov近似，将上层地下裂缝介质的背景波场u0(z1)与前向散射场代入将计算结果作为下一层地下裂缝介质的入射波场，将新的入射波场代入虚数域rytov近似散射波传播方程(3)和(4)，分别计算下一层地下裂缝介质的前向散射场和背向散射场并且对下一层地下裂缝介质的背向散射场继续进行存储。

步骤5：根据步骤3和步骤4计算虚数域粗网格rytov近似散射波场，当模拟的波场传播到地下裂缝介质的加密边界时，利用插值原理将水平方向上的粗网格波场进行插值加密处理，直至满足网格加密倍数，垂直方向上对粗网格进行分割，并基于波场外推理论，对加密网格中地下裂缝介质的波场逐层进行外推计算：

其中，波场由粗网格插值得到加密网格的算法为：

波场由加密网格转换到粗网格的算法为：

式中，k为水平方向网格加密的倍数；m为粗网格点次序；j为加密网格点次序；n为粗网格点数；un^m为粗网格中第m个点的波场值；ukn^km为加密网格波场；ukn^km+j为km+j点在加密网格中的波场；ukn^k(m+1)为k(m+1)点在加密网格中的波场；为l点在加密网格中的波场。

步骤6：当对加密网格中地下裂缝介质的波场外推计算至地下裂缝介质的加密边界时，基于转换算法，将水平方向上加密网格的波场变换至粗网格中，垂直方向上则继续基于波场外推理论对加密边界下粗网格中地下裂缝介质的波场逐层进行外推计算。

步骤7：当波场传播直到地下裂缝介质的底部时，地下裂缝介质的波场外推计算完毕；此时，从地下裂缝介质的底部开始，利用地下裂缝介质底部计算出的背向散射场作为上一层地下裂缝介质的入射波场，利用式(4)计算向后传播的背向散射场

步骤8：将向后传播的背向散射场与步骤4存储下来的n-1层地下裂缝介质的背向散射场进行叠加处理，将叠加结果作为上一层地下裂缝介质的入射波场，利用式(3)计算向后传播的前向散射场

步骤9：重复步骤4～8，波场向后传播，当到达地下裂缝介质的接收位置时，地下裂缝介质的波场的反向外推计算完毕。最后利用二维傅里叶反变换将虚数域地震散射记录转换至实数域地震散射记录，对时间-空间域的波场特征进行表征，将反变换后的波场结果进行存储，最终获得接收地下裂缝介质地震散射记录。

图3展示了用于正演模拟的含裂缝地层速度模型，本发明利用该速度模型分别计算了裂缝角度为从0°-180°变换的合成地震记录(角度为水平线顺时针与裂缝的夹角)，其它参数保持不变。部分结果如图4(a)-(f)所示，图4(a)-(f)依次为裂缝角度为0°、30°、60°、90°、120°和150°时的合成地震记录。当裂缝的走向水平时，由于裂隙存在产生的地震同相轴最为连续，裂缝垂直时，连续性最差且能量相对集中。

图5(a)-(d)展示了裂缝在地层界面不同位置处的地面接收记录，图中裂缝共有三组，每组七条，从图中可以看出，裂缝在模型中地层界面不同位置时，其地面接收记录有所差异。当裂缝在地层反射界面上方时，会增加地层界面的地震响应，当裂缝在地层界面中间和下面，会减弱地层界面的地震响应。

上述各实施例仅用于对本发明的目的、技术方案和有益效果进行示例性描述，并不局限于上述具体实施方式，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。