CN100545679C - 起伏地表的地震叠后正演方法 - Google Patents

Abstract

起伏地表的地震叠后正演方法，对起伏地表的地质模型进行矩形封装，起伏地表上方采用任意速度值或密度值进行充填，再采用斯奈尔定理计算反射系数模型；消除起伏地表的反射影响并变换到频率域作为初始波场，利用单程波动方程的傅立叶有限差分延拓算子沿着深度轴从下往上进行波场延拓，且保持输入基准面(起伏地表)和输出基准面之间的波场为零；将延拓后的频率域的波场变换到时间域进行成像，通过常规的显示软件将成像结果处理为地震剖面图像，对于起伏地表的地质模型具有很好的地震波场模拟效果。

Description

起伏地表的地震叠后正演方法

技术领域

本发明涉及反射波地震波场数值模拟过程中的地震叠后正演技术范畴，是一种起伏地表地震波场数值模拟的叠后正演方法。

背景技术

地震波场的数值模拟就是在假定已知地下介质的结构模型和相应物理参数的情况下，模拟地震波在地下各种介质中的传播规律，并计算在地面或地下各观测点所应观测到的数值地震记录的一种地震模拟方法。它一方面可以为地震数据采集、处理、解释提供理论依据，评估方法的科学性、可行性和先进性；另一方面可以用来检验各种解释成果的可信性，以及各种反演算法的正确性和可靠性。从地震波传播的角度，可以将常规的地震叠后正演技术分为射线追踪方法和波动方程方法。射线追踪方法可以灵活地处理起伏地表，但是它以地震波的运动学为基础，将地震波波动理论简化为射线理论，主要考虑地震波传播的运动学特征，缺少地震波的动力学信息，计算结果在一定程度上是近似的，应用有一定限制条件，对复杂的地质构造成像精度低，有时会出现盲区。波动方程方法通过数值方法直接求解波动方程，包含了地震波传播的动力学特征，波场信息丰富，模拟结果较为准确，但是它只能处理水平地表的情况，波场延拓的基准面必须为一个水平界面。

发明内容

本发明是要提供一种基于单程波动方程的适用于起伏地表的地震叠后正演方法，它不但能够适应起伏的地表条件，而且采用单程波动方程进行计算，具有很快的计算速度和很高的复杂构造的成像精度。

本发明的基于单程波动方程的起伏地表的地震叠后正演技术，具体步骤包括：

(1)采用矩形框对目标地质模型进行整体封装，矩形框的顶面即为输出的基准面，起伏地表上方采用任意速度值或密度值进行充填，并采用规则网格对速度模型和密度模型进行网格化。

(2)采用斯奈尔定理计算反射系数，并利用雷克子波进行褶积处理，得到反射系数模型。

(3)消除起伏地表的反射影响，并变换到频率域作为初始延拓波场。

(4)利用单程波动方程的傅立叶有限差分延拓算子沿着深度轴从下往上进行波场延拓，且保持输入基准面(起伏地表)和输出基准面之间的波场为零。

(5)将延拓后的频率域的波场变换到时间域进行成像。

(6)通过常规的显示软件将成像结果处理为地震剖面图像。

本发明的基于单程波动方程的起伏地表的地震叠后正演方法，在输入基准面和输出基准面之间保持一个数值为零的波场层来解决波动方程地震叠后正演中的起伏地表问题，吸取了射线追踪方法灵活处理起伏地表的优点。

本发明的基于单程波动方程的起伏地表的地震叠后正演方法，利用傅立叶有限差分延拓算子在频率域进行波场延拓计算，具有较强的处理横向变速的能力，对于复杂构造具有很高的成像精度，改善了射线追踪方法对于复杂构造成像精度低的缺点。

本发明的具体实现原理如下：

选定的矩形框的顶面为波场延拓的基准面，其为一个水平界面，且位于起伏地表的最高点之上；矩形框的右边界、左边界和底边界与输入的地质模型完全重合。对封装后的地质模型进行规则的网格剖分，于是起伏的地表形态被离散化后成为规则网格的一部分，起伏地表之上的网格采用任意的速度值和密度值进行充填。利用斯奈尔定理计算反射系数，

$R(x,y,z_i,t)=\frac{\mathrm{Den}(x,y,z_{i+1},t)\×\mathrm{Vel}(x,y,z_{i+1},t)-\mathrm{Den}(x,y,z_i,t)\×\mathrm{Vel}(x,y,z_i,t)}{\mathrm{Den}(x,y,z_{i+1},t)\×\mathrm{Vel}(x,y,z_{i+1},t)+\mathrm{Den}(x,y,z_i,t)\×\mathrm{Vel}(x,y,z_i,t)}$

并采用雷克子波(主频根据目标工区的地层特征选定，一般介于10Hz～90Hz之间)进行褶积处理。

消除起伏地表的反射影响，并采用傅立叶变换将它变换到频率域R(x，y，z，ω)，作为初始延拓的波场。给定频率ω，提取该频率处的初始波场W(x，y，z，ω)，用上行波方程沿着深度轴采用傅立叶有限差分波场延拓算子从下往上进行波场延拓，延拓后的波场为W′(x，y，0，ω)，该过程可以表示为，

W′(x，y，0，ω)＝FFD_z[W(x，y，z，ω)]

其中x，y分别是空间直角坐标系中的两个水平方向的坐标，z是深度方向的坐标，ω是圆频率，FFD_z是傅立叶有限差分波场延拓算子，Den(x，y，z，t)是网格化的密度值，Vel(x，y，z，t)是网格化的速度值。延拓过程中下一深度z_i的输入波场为上一深度z_i+1的延拓结果，同时在输出基准面和输入基准面之间设置一个零波场层参与正常的波场传播和延拓过程。

对所有的频率进行循环计算，得到输出基准面上的频率域波场U(x，y，0，ω)。对U(x，y，0，ω)进行反傅立叶变换变换到时间域成像，就得到最后的模拟波场I(x，y，t)，然后通过常规的显示软件将成像结果处理为地层剖面图像。

上述延拓算子是根据单程波动方程近似逼近得到的，采用频率域的稳定的傅立叶有限差分算子实现(Ristow，1994)。该波场延拓算子可近似为，

$\sqrt{\frac{{\mathrm{\ω}}^2}{v^2}+\frac{{\∂}^2}{\∂x^2}}=A_1+A_2+A_3$

式中，

$A_1=\sqrt{\frac{{\mathrm{\ω}}^2}{c^2}+\frac{{\∂}^2}{\∂x^2}};$ $A_2=(\frac{\mathrm{\ω}}{v}-\frac{\mathrm{\ω}}{c});$ $A_3=\frac{\mathrm{\ω}}{v}[1-\frac{c}{v}]\frac{a\frac{v^2}{{\mathrm{\ω}}^2}\frac{{\∂}^2}{\∂x^2}}{1+b\frac{v^2}{{\mathrm{\ω}}^2}\frac{{\∂}^2}{\∂x^2}}$

其中v是实际速度场，c是参考速度场，ω是角频率，参数a，b的取值影响到波场成像的最大倾角：

本发明采用基于单程波动方程的适应于起伏地表的地震波场叠后数值模拟方法，具有如下特点，主要表现为：

(1)可以灵活地处理地震叠后正演中的起伏地表条件。

(2)将起伏地表、单程波动方程和频率域的傅立叶有限差分波场延拓算子有机结合。

(3)对于起伏地表条件下的复杂构造具有精确的波场成像效果。

附图说明

图1是是一个水平地层的二维起伏地表模型，地表的最大起伏高差为380米。

图2是用本技术对图1中的速度模型进行地震叠后正演模拟后的结果，可以看出水平地层的反射同相轴不再平直，而是与地表形态相反；此外，模拟仅包含一次反射波，信噪比高。

图3是一个复杂地质构造的二维起伏地表模型，近地表处包含两个速度为2000米/秒的低速异常体。

图4是用本技术对图3中的速度模型进行地震叠后正演模拟后的结果，可以看出近地表的两个低速异常体的反射时间基本一致，不受起伏地形的影响；此外，由于起伏地形和低速异常体的影响，位于模型中部的水平地层的反射同相轴发生了严重的扭曲。

具体实施方式

基于单程波动方程的起伏地表的地震叠后正演方法，具体为以下步骤：

(1)采用矩形框对目标地质模型进行整体封装，矩形框的顶面即为输出的基准面，起伏地表上方采用任意速度值或密度值进行充填，并采用规则网格对速度模型和密度模型进行网格化。如果只有速度模型，则采用密度与速度的经验公式计算密度模型。

(2)采用斯奈尔定理计算反射系数，并利用雷克子波进行褶积处理，得到反射系数模型。雷克子波的主频根据工区的地层特点进行选择(一般介于10Hz～90Hz之间)。

(3)消除起伏地表的反射影响，并变换到频率域作为初始延拓波场。

(4)利用单程波动方程的傅立叶有限差分延拓算子沿着深度轴从下往上进行波场延拓，且保持输入基准面(起伏地表)和输出基准面之间的波场为零。

具体是：进行频率循环内的波场延拓计算，并根据工区的地层特点，参考初始波和终止波的时间，确定波场延拓的频率范围f1和f2，f1一般介于1Hz～5Hz之间，f2一般介于50Hz～100Hz之间；在单个频率循环内，沿深度轴从下往上进行波场延拓计算，确定傅立叶有限差分波场延拓算子 $\sqrt{\frac{{\mathrm{\ω}}^2}{v^2}+\frac{{\∂}^2}{\∂x^2}}=A_1+A_2+A_3,$ 其中：

$A_1=\sqrt{\frac{{\mathrm{\ω}}^2}{c^2}+\frac{{\∂}^2}{\∂x^2}};$ $A_2=(\frac{\mathrm{\ω}}{v}-\frac{\mathrm{\ω}}{c});$ $A_3=\frac{\mathrm{\ω}}{v}[1-\frac{c}{v}]\frac{a\frac{v^2}{{\mathrm{\ω}}^2}\frac{{\∂}^2}{\∂x^2}}{1+b\frac{v^2}{{\mathrm{\ω}}^2}\frac{{\∂}^2}{\∂x^2}}$

其中v是实际速度场，c是参考速度场，ω是角频率，并根据地层的倾斜程度选择参数a，b的取值，参数a，b的取值和偏移成像的最大倾角的对应关系为：

延拓过程中下一深度z_i的输入波场为上一深度z_i+1的延拓结果，同时在输出基准面和输入基准面之间设置一个零波场层参与正常的波场传播和延拓过程。循环计算完f1和f2之间的波场值，对其它的波场值充零处理。

(5)将延拓后的频率域的波场变换到时间域进行成像。

(6)通过常规的显示软件将成像结果处理为地震剖面图像。

图1和图2是一个水平地层的二维起伏地表模型的测试例子。

(1)对矩形封装的速度模型进行网格化，起伏地表上方的速度给定为500米/秒(任意选定的值)，密度模型为常数1；根据道距确定网格X方向的间距为20米，根据深度确定网格Z方向的间距为5米。

(2)采用斯奈尔定理计算反射系数，并利用主频为50Hz的雷克子波进行褶积处理，得到反射系数模型。

(3)消除起伏地表的反射影响，并变换到频率域作为初始延拓波场。

(4)利用单程波动方程的傅立叶有限差分延拓算子沿着深度轴从下往上进行波场延拓，且保持输入基准面(起伏地表)和输出基准面之间的波场为零。具体的说，根据工区的地层特点和参考初始波和终止波的时间，确定循环计算的频率范围f1和f2，一般将它设置为地震数据的有效频带，这里f1＝1，f2＝100。由于地层为水平地层，参数a，b，取值倾角为45度的情况，为a＝0.5，b＝0.25。对f1和f2之间的每个频率的初始波场沿深度轴从下往上进行波场延拓计算，其波场延拓算子为前面叙述的稳定的频率域的傅立叶有限差分算子。延拓过程中下一深度z_i的输入波场为上一深度z_i+1的延拓结果，同时在输出基准面和输入基准面之间设置一个零波场层参与正常的波场传播和延拓过程。循环计算完f1和f2之间的波场值，对其它的波场值充零处理。

(5)将延拓后的频率域的波场变换到时间域进行成像。

(6)通过常规的显示软件将成像结果处理为地震剖面图像。图2是地震叠后模拟的地震剖面。

图3和图4是一个水平地层的二维起伏地表模型的测试例子。

(1)对矩形封装的速度模型进行网格化，起伏地表上方的速度给定为300米/秒(任意选定的值)，密度模型为常数1；根据道距确定网格X方向的间距为20米，根据深度确定网格Z方向的间距为5米。

(2)采用斯奈尔定理计算反射系数，并利用主频为50Hz的雷克子波进行褶积处理，得到反射系数模型。

(3)消除起伏地表的反射影响，并变换到频率域作为初始延拓波场。

(4)利用单程波动方程的傅立叶有限差分延拓算子沿着深度轴从下往上进行波场延拓，且保持输入基准面(起伏地表)和输出基准面之间的波场为零。具体的说，根据工区的地层特点和参考初始波和终止波的时间，确定循环计算的频率范围f1和f2，一般将它设置为地震数据的有效频带，这里f1＝1，f2＝100。由于地层为水平地层，参数a，b，取值倾角为45度的情况，为a＝0.5，b＝0.25。对f1和f2之间的每个频率的初始波场沿深度轴从下往上进行波场延拓计算，其波场延拓算子为前面叙述的稳定的频率域的傅立叶有限差分算子。延拓过程中下一深度z_i的输入波场为上一深度z_i+1的延拓结果，同时在输出基准面和输入基准面之间设置一个零波场层参与正常的波场传播和延拓过程。循环计算完f1和f2之间的波场值，对其它的波场值充零处理。

(5)将延拓后的频率域的波场变换到时间域进行成像。

(6)通过常规的显示软件将成像结果处理为地震剖面图像。图4是地震叠后模拟的地震剖面。

Claims (4)

1、起伏地表的地震叠后正演方法，其特征在于采用以下步骤：

(1)采用矩形框对目标地质模型进行整体封装，矩形框的顶面即为输出的基准面，矩形顶面位于起伏地表的最高点之上，起伏地表上方采用任意速度值或密度值进行充填，并采用规则网格对速度模型和密度模型进行网格化，如果只有速度模型，则采用密度与速度的经验公式计算密度模型；

(2)采用斯奈尔定理计算反射系数，并利用雷克子波进行褶积处理，得到反射系数模型，雷克子波的主频根据工区的地层特点进行选择，一般介于10Hz～90Hz之间；

(3)将起伏地表的反射系数赋值为0，消除起伏地表的反射影响，并对反射系数数组进行正傅立叶变换，变换到频率域作为初始延拓波场；

(4)利用单程波动方程的傅立叶有限差分延拓算子沿着深度轴从下往上进行波场延拓，当延拓深度z_i介于输入基准面与输出基准面之间时，在深度z_i+1与z_i之间延拓的波场增量值赋值为零，即保持输入基准面和输出基准面之间的波场延拓增量为零，根据工区的地层特点，参考初始波和终止波的时间，确定波场延拓的频率范围f1和f2，f1一般介于1Hz～5Hz之间，f2一般介于50Hz～100Hz之间，对其它频率的波场值充零处理，

(5)将延拓后的频率域的波场进行反傅立叶变换，拾取时间域波场的实部，得到最后的成像结果；

(6)通过常规的显示软件将成像结果处理为地震剖面图像。

2、根据权利要求1所述的起伏地表的地震叠后正演方法，其特征在于：对目标模型进行矩形封装，矩形顶面位于起伏地表的最高点之上；并通过在输入基准面和输出基准面之间设置一个零波场层来参与正常的波场传播和延拓，即当延拓深度z_i介于输入基准面与输出基准面之间时，在深度z_i+1与z_i之间延拓的波场增量值赋值为零。

3、根据权利要求1所述的起伏地表的地震叠后正演方法，其特征在于：采用斯奈尔定理计算反射系数，并利用雷克子波进行褶积处理，得到反射系数模型，雷克子波的主频根据目标工区的地层特征选定，一般介于10Hz～90Hz之间。

4、根据权利要求1或2所述的起伏地表的地震叠后正演方法，其特征在于：采用频率域的傅立叶有限差分方法算子沿深度轴从下往上进行波场延拓计算。

Images