一种时空双变正演模拟方法与流程

本发明公开一种时空双变正演模拟方法，属于地震勘探基础应用领域。

背景技术：

有限差分正演模拟方法是目前处理常规非均匀介质应用较为广泛的一种正演方法，将地质模型网格化得到数学模型，通过数值计算，求解描述地震波传播的微分方程获取地震响应波场，能够有效模拟地震波传播过程中的绕射波、多次波等全波场信息。除了对地质模型进行空间网格化，在正演模拟过程中还需要进行时间域的网格离散，通常来说，时空域网格尺度越小，正演模拟精度越高，计算稳定性越强，但计算量和占用内存会增大，正是这种模拟精度、稳定性与计算效率、占用内存之间的矛盾一直推动着正演模拟技术的发展。

现今，勘探对象逐步向强纵横向变速区域、低速带、复杂构造区域及碳酸盐岩储层的小型孔缝洞等地质体转变，上述模拟精度与计算效率之间的矛盾愈加突出。1989年，Moczo提出了可变网格的思想，即模拟不同区域采用不同的空间网格尺度，该方法被证明是保持模拟精度同时又降低内存需求的可行方法。随后，国内外大量学者就该空间变网格思想做了深入的研究，增强了变网格算法的适应性和实用性(Jastram，1992，1994；Pitarka，1999；Aoi，1999；Wang，2001；李胜军，2007；朱生旺，2007；赵海波，2007；孙成禹，2008；李振春，2008)。

为了保证算法的稳定性，空间小网格尺度要求较小的时间网格步长，但在非变网格区域采用小时间步长会造成时间过采样，增加计算波场延拓次数，同时研究表明，这不仅不会提高精度，还会在一定程度上引入频散误差(Collino，2003)。为了克服这一不足，Falk(1998)提出了局部可变时间步长算法，允许变时间步长为2的幂级数倍。Tessmer(2000)基于二阶声波方程做了进一步的改进，时间步长变化倍数可为任意整数。

专利《高精度空间和时间任意倍数可变网格有限差分正演方法(CN105277980)》公开了一种基于二阶声波运动方程的高精度空间和时间任意倍数可变网格有限差分正演方法，可适应低速层、裂缝介质、生物礁等各类复杂模型。但二阶方程常规网格离散解的频散条件较一阶交错网格的更严格，需要更小的网格尺度，占用内存和计算量较大。

针对这一问题，专利《基于时空双变网格的弹性波正演模拟技术(CN102183790)》公开了一种基于交错网格的双变弹性波正演模拟技术，利用一阶速度-应力方程，对复杂模型采用局部精细网格剖分方案，并对空间精细区域进行精细时间步长延拓，提高了计算效率。该技术未考虑高倍网格变化带来的不稳定，同时，当存在多个相距较远的需要局部精细化的目标区时，覆盖区统一变网格剖分方案对计算效率的改善不明显。

地球物理学进展杂志29卷第3期公开了一种《裂缝型储层时空双变正演模拟研究》，该方法将Lanczos滤波算子引入到时空双变正演模拟算法中，利用优化的正演模拟算法模拟裂缝储层地震响应，解决了常规变网格算法长时间采样下的不稳定问题，对深层微构造有较好的模拟精度及效率，但该方法同样未考虑高倍网格变化下的稳定性及多目标区复杂构造问题。

Applied Geophysics(应用地球物理(英文版))杂志58卷第1期公开了一种《基于时空双变网格的起伏地表变坐标系正演模拟方法》，该方法通过坐标变换法将起伏地表转化为水平地表,同时将物理空间的波动方程转化为计算空间的波动方程,在计算空间采样变网格技术完成数值模拟，较全局细网格算法能显著节约计算内存，对起伏地表构造具有较高的模拟精度和一定的适应性，拓展了时空双变技术的应用范围。

综上所述，现有的变网格技术在网格变化倍数较高或模拟时间较长时，易产生不稳定，此外，当研究的地质构造中存在相距较远的多目标区时，采用同一变网格剖分方案对计算效率和内存的提升不明显，因此，变网格技术在其稳定性和计算效率等方面仍有待改善。

技术实现要素：

本发明的目的是提供一种克服上述变网格技术在高倍变网格和长时间采样下稳定性较低，多目标复杂地质构造条件下计算量和内存消耗较大等缺陷的时空双变正演模拟方法。

本发明利用待评价地质构造的深度域速度场表征的物性参数，对深度域速度场进行分块、分级，在建立深度域速度场的背景网格模型、单级变网格分块的变网格模型和多级变网格分块的多级变网格模型的基础上，设定人工震源，通过时空双变正演方法，获得由每个分块不同变网格尺度下的波场特征确定的背景网格地震响应波场，具体包括以下步骤：

1、利用待评价地质构造的地震、地质、测井资料，确定深度域速度场、深度域速度场的背景网格尺度和背景时间步长，建立深度域速度场的背景网格模型；

2、利用待评价地质构造的深度域速度场，确定深度域速度场的各单级变网格分块的单级变网格尺度和变网格时间步长，多级变网格分块中各级变网格尺度和变网格时间步长，建立深度域速度场的单级变网格分块的单级变网格模型和多级变网格分块的多级变网格模型。

2.1利用步骤1获得的深度域速度场表征的待评价地质构造的物性参数及目标区数量，确定待评价地质构造的分块数量、各分块变网格尺度和变网格倍数；

2.2利用步骤2.1得到的各分块变网格尺度和变网格倍数，将所有分块分为单级变网格分块和多级变网格分块，并确定多级变网格分块的变网格级数和各级变网格倍数；

2.3利用步骤2.2确定的各单级变网格分块的变网格倍数、多级变网格分块的变网格级数和各级变网格倍数，和步骤1确定的深度域速度场的背景网格尺度和背景时间步长，确定深度域速度场的各单级变网格分块的单级变网格尺度和变网格时间步长，多级变网格分块中各级变网格尺度和变网格时间步长，建立深度域速度场的各单级变网格分块的单级变网格模型和多级变网格分块的多级变网格模型。

3、设定人工震源，利用二维声波压力-速度波动方程离散差分式，正演模拟其在深度域速度场的背景网格模型、各单级变网格模型、多级变网格模型中的波场传播，确定包含每个分块不同变网格尺度下的波场特征的背景网格地震响应波场。

3.1利用步骤1生成的深度域速度场的背景网格尺度和背景时间步长，以及步骤2得到的深度域速度场的各单级变网格分块的单级变网格尺度和变网格时间步长，多级变网格分块中各级变网格尺度和变网格时间步长，建立背景网格、各单级变网格以及多级变网格中各级变网格的时间二阶，空间偶数阶差分精度的二维声波压力-速度波动方程离散差分式；

3.2设定人工震源，利用步骤3.1得到的背景网格的时间二阶，空间偶数阶差分精度的二维声波压力-速度波动方程离散差分式，更新计算背景网格尺度和时间网格步长上的压力场和速度场值，获得背景网格模型产生的地震响应波场；

3.3将步骤3.2得到的背景网格尺度和时间网格步长上的压力场和速度场值传递给各分块变网格模型边界上的对应点，同时，分别判断波场是否传播到该分块所在区域，若未传播到，则该区域仍按步骤3.2进行更新，若传播到单级变网格分块所在区域，则进入步骤3.4，若传播到多级变网格分块所在区域，则进入步骤3.5；

3.4设定步骤3.3获得的单级变网格模型边界上的压力场和速度场值为变网格正演的初始值和边界值，利用步骤3.1得到的该单级变网格时间二阶，空间偶数阶差分精度的二维声波压力-速度波动方程离散差分式，更新计算该分块区域内变网格尺度和时间网格步长上的压力场和速度场，获得该单级变网格分块模型产生的地震响应波场，进入步骤3.6；

3.5进入步骤3.3确定的波场已传播到的多级变网格分块所在区域，利用步骤3.4相同的原理，获得该多级变网格分块第一级变网格模型产生的地震响应波场，同时，判断波场是否传播到该多级变网格分块的第二级变网格区域，若未传播到，则仍按第一级变网格模型进行波场更新，若传播到，则利用步骤3.4相同的原理，获得该多级变网格分块第二级变网格模型产生的地震响应波场，依次类推，直至获得该多级变网格分块最后一级变网格模型产生的地震响应波场；

3.6利用Lanczos滤波公式，将步骤3.4获得的单级变网格分块模型产生的地震响应波场，和步骤3.5获得的多级变网格分块各级变网格模型产生的地震响应波场，传递给对应的背景网格，获得同时包含单级变网格分块的单级变网格尺度和多级变网格分块中各级变网格尺度下的波场特征的背景网格地震响应波场，进入步骤3.7；

3.7利用步骤3.6获得的同时包含单级变网格分块的单级变网格尺度和多级变网格分块中各级变网格尺度下的波场特征的背景网格地震响应波场，按照步骤3.2-3.6迭代，完成所有背景时间步长的地震响应波场更新，获得由每个分块不同变网格尺度下的波场特征确定的背景网格地震响应波场。

本发明的有益效果是：本发明采用Lanczos滤波公式完成不同尺度网格间的波场传递，能够保证长时间采样下波场延拓的稳定性、减小变网格界面的虚假反射误差；同时，通过多级交错变网格技术解决了高倍网格变化带来的不稳定，保证了对微小尺度地质目标体的模拟精度；此外，分块变网格思想可同时对包含若干个尺度级数不同的目标区进行模拟，每个目标区采用不同的变网格倍数，最大限度的提高正演模拟的效率，增强变网格技术的适用性。

附图说明

图1本发明技术方案流程图；

图2深度域速度场的背景网格模型；

图3(a)深度域速度场的地表起伏带分块对应的背景网格模型；

图3(b)深度域速度场的地表起伏带分块的单级3倍变网格模型；

图4(a)深度域速度场的深部低速层分块对应的背景网格模型；

图4(b)深度域速度场的深部低速层分块的第一级3倍变网格模型；

图4(c)深度域速度场的深部低速层的第二级3*5倍变网格模型；

图4(d)深度域速度场的深部低速层的第三级3*5*11倍变网格模型；

图5(a)本发明方法得到的不同变网格尺度下的波场特征确定的背景网格地震响应波场；

图5(b)深度域速度场的深部低速层分块内裂缝带的地震响应波场；

图6常规网格正演模拟地震响应波场；

图7三个偏移距下，两种不同正演地震响应波场的单道波形对比；

图8(a)、(b)、(c)分别为采用直接传递法、九点加权法和Lanczos滤波法三种滤波方法正演模拟6-8s的地震响应波场；

图9(a)三种建模方式占用内存对比；

图9(b)三种正演方法计算时间对比。

具体实施方式

以一个待评价的地质构造问题为实例并结合附图，对本发明做进一步描述，由图1可知本发明实施例如下：

1、获取待评价地质构造的地震、地质、测井资料，确定深度域速度场：深度域速度场大小为1.8km*1.8km，包含两个目标区：地表起伏带和深部低速层(厚度为12m)，同时深部低速层局部发育5m厚的微小尺度裂缝带，目标区之外为两个水平地层，其最大速度为4000m/s，最小速度为3000m/s，地震主频为30Hz，利用公式(1)和公式(2)得到深度域速度场的背景网格尺度为6m，时间步长为0.2ms，由此建立深度域速度场的背景网格模型，如图2所示，该深度域速度场的背景网格模型的网格大小为301*301，第一层为地表起伏带，纵向网格点241处为深部低速层，由于深度域速度场的背景网格尺度大于裂缝带厚度，无法刻画裂缝带，背景网格模型中不能识别裂缝带。

Δx:深度域速度场的背景网格尺度，Δt:深度域速度场的背景时间步长，v_min：深度域速度场的最小速度，v_max：深度域速度场的最大速度，feq：地震主频。

2、利用待评价地质构造的深度域速度场，确定深度域速度场的各单级变网格分块的单级变网格尺度和变网格时间步长，多级变网格分块中各级变网格尺度和变网格时间步长，建立深度域速度场的单级变网格分块的变网格模型和多级变网格分块的变网格模型。

2.1利用步骤1获得的深度域速度场表征的待评价地质构造的地震主频30Hz，两个目标区地表起伏带和深部低速层的速度分别为3000m/s、2500m/s，裂缝开度0.36cm，裂缝内含油的速度为1300m/s，确定变网格区为：网格范围(1-301)*(1-81)的地表起伏带和网格范围(1-301)*(220-260)的含裂缝的深部低速层两个分块，利用上述目标区的物性参数、公式(1)和公式(2)，确定深度域速度场的地表起伏带分块和含裂缝的深部低速层分块的网格尺度分别为2m、0.36cm，变网格倍数分别为3倍、165倍。

2.2利用步骤2.1得到的深度域速度场的地表起伏带分块和含裂缝的深部低速层分块的变网格尺度和变网格倍数，将地表起伏带定义为单级变网格分块，含裂缝的深部低速层定义为多级变网格分块，分三级，第一级3倍变网格，第二级3*5倍变网格，第三级3*5*11倍变网格；

2.3利用步骤2.2确定的深度域速度场的地表起伏带分块的单级变网格倍数、深度域速度场的含裂缝的深部低速层分块的多级变网格级数和各级变网格倍数，和步骤1确定的深度域速度场的背景网格尺度和背景时间步长，利用公式(2)确定深度域速度场的地表起伏带分块的单级变网格时间步长为0.067ms，深度域速度场的含裂缝的深部低速层分块中第一级变网格尺度和变网格时间步长分别为2m、0.067ms，第二级变网格尺度和变网格时间步长分别为0.4m、0.013ms，第三级变网格尺度和变网格时间步长分别为0.36cm、0.0012ms，分别建立深度域速度场的地表起伏带分块的单级变网格模型，如图3(b)所示，该模型中的地表起伏带构造更平滑，消除了背景网格模型中出现的起伏界面阶梯毛刺，提高了构造表征精度；含裂缝的深部低速层分块的多级变网格模型，如图4所示，图4(b)为第一级3倍变网格模型，该模型中的深部低速层比图4(a)背景网格表征下的更精细，但由于裂缝开度很小，无法表征深部低速层内的裂缝带，图4(c)为第二级3*5倍变网格模型，该模型中的裂缝带刻画比较模糊，图4(d)第三级3*5*11倍变网格模型，该模型中可以清晰分辨裂缝的倾角和缝长等产状，达到了对裂缝的表征精度。

3、设定人工震源，利用二维声波压力-速度波动方程离散差分式，正演模拟其在深度域速度场的背景网格模型、地表起伏带分块的单级变网格、含裂缝的深部低速层分块的多级变网格模型中的波场传播过程，确定包含每个分块不同变网格尺度下的波场特征的背景网格地震响应波场。

3.1利用步骤1生成的深度域速度场的背景网格尺度和背景时间步长，以及步骤2得到的深度域速度场的地表起伏带分块的单级变网格尺度和变网格时间步长，含裂缝的深部低速层分块的多级变网格中各级变网格尺度和变网格时间步长，建立背景网格、地表起伏带分块的单级变网格以及含裂缝的深部低速层分块的多级变网格中各级变网格的时间二阶，空间偶数阶差分精度的二维声波压力-速度波动方程离散差分式：

τ:压力场，ν_x:x方向速度场,ν_z:z方向速度场，υ_d:不同网格模型速度，S_d:不同网格模型变网格尺度，n_d:不同网格模型时间点坐标，i_d:不同网格模型x方向网格坐标，j_d:不同网格模型z方向网格坐标，D_x,D_z：x,z方向的四阶差分算子：

f(x,z):压力场或x、z方向速度场，Δx，Δz：不同网格模型网格尺度，c_m：四阶Taylor中心差分系数。

3.2设定人工震源，利用步骤3.1得到的背景网格的时间二阶，空间偶数阶差分精度的二维声波压力-速度波动方程离散差分式，更新计算背景网格尺度和时间网格步长上的压力场和速度场值，获得背景网格模型产生的地震响应波场；

3.3将步骤3.2得到的背景网格尺度和时间网格步长上的压力场和速度场值分别传递给地表起伏带分块和含裂缝的深部低速层分块变网格模型边界上的对应点，同时，分别判断波场是否传播到各个分块所在区域，若未传播到，则该区域仍按步骤3.2进行更新，若传播到单级变网格分块所在区域，则进入步骤3.4，若传播到多级变网格分块所在区域，则进入步骤3.5；

3.4以步骤3.3获得的地表起伏带分块的单级变网格模型边界上的压力场和速度场值为变网格正演的初始值和边界值，利用步骤3.1得到的该单级变网格时间二阶，空间偶数阶差分精度的二维声波压力-速度波动方程离散差分式，更新计算该分块区域内变网格尺度和时间网格步长上的压力场和速度场，获得地表起伏带分块的单级变网格模型产生的地震响应波场，进入步骤3.6；

3.5进入步骤3.3确定的波场已传播到的含裂缝的深部低速层分块的多级变网格模型所在区域，利用步骤3.4相同的原理，获得该多级变网格分块第一级变网格模型产生的地震响应波场，同时，判断波场是否传播到该多级变网格分块的第二级变网格区域，若未传播到，则仍按第一级变网格模型进行波场更新，若传播到，则利用步骤3.4相同的原理，获得该多级变网格分块第二级变网格模型产生的地震响应波场，依次类推，获得该多级变网格分块第三级变网格模型产生的地震响应波场；

3.6利用Lanczos滤波公式(4)，将步骤3.4获得的地表起伏带分块的单级变网格模型产生的地震响应波场，和步骤3.5获得的含裂缝的深部低速层分块的多级变网格模型中各级变网格模型产生的地震响应波场，传递给对应的背景网格，获得同时包含地表起伏带分块的单级变网格尺度和含裂缝的深部低速层分块的多级变网格尺度下的波场特征的背景网格地震响应波场，进入步骤3.7；

k：变网格倍数，A由确定，F(i,j)：背景网格点上的压力场和速度场值，f(i,j)：变网格点上的压力场和速度场值，ω_mn:Lanczos滤波算子。

3.7利用步骤3.6获得的同时包含地表起伏带分块的单级变网格尺度和含裂缝的深部低速层分块的多级变网格尺度下的波场特征的背景网格地震响应波场，按照步骤3.2-3.6迭代，完成所有背景时间步长的地震响应波场更新，获得由每个分块不同变网格尺度下的波场特征确定的背景网格地震响应波场，如图5(a)所示。

比较例：图6为采用6m网格尺度对待评价地质构造模型进行常规网格正演模拟，获得的待评价地质构造的地震响应波场。通过图5(a)和图6对比，利用本发明方法正演模拟精度更高，可以消除背景网格尺度下的起伏界面阶梯毛刺产生的边界绕射波噪音，同时能够实现对小尺度裂缝的刻画，图5(b)为裂缝带产生的地震响应波场。

图7为三个不同偏移距下，两个地震响应波场的单道波形，两个地震响应波场分别采用本发明方法(目标区采用15倍变网格、0.4m网格尺度和3倍变网格、2m网格尺度，非目标区采用6m网格尺度)和全区0.4m网格尺度的常规网格正演模拟方法获得，从波形对比可以发现：本发明方法得到正演结果和全局0.4m网格尺度正演结果吻合较好，验证了本发明方法的模拟精度。

图8为分别利用现有的直接传递法(a)、九点加权法(b)和本发明采用的Lanczos滤波法(c)，正演模拟在6s-8s大时间采样时的地震响应波场。通过对比，本发明采用的Lanczos滤波方法在大时间采样情况下更稳定。

图9(a)对比了全区165倍变网格、地表起伏带和含裂缝的深部低速层两个目标区不分块看作一个目标区的165倍变网格以及本发明提出的多级、分块变网格，三种网格模型的内存占用量，通过对比，本发明方法可以降低99.97％的占用内存，效果显著。

图9(b)对比了全区3倍变网格、地表起伏带和含裂缝的深部低速层两个目标区不分块看作一个目标区的3倍变网格以及本发明提出的多级、分块变网格，三种网格模型正演模拟的计算耗时，通过对比，本发明方法可以节约98.37％的计算时间，提高计算效率。

上述实施例说明本发明技术可以更好的适应复杂小尺度裂缝、孔洞型储集体、深层碳酸盐岩储层、多目标区、起伏地表低降速带等非均质较强的复杂地质构造问题。引入的Lanczos滤波方法较好的保证了模拟的稳定性，以及在变网格边界处的模拟精度，而多级分块思想相较于常规变网格方法又进一步提高了对微尺度储层以及复杂地质体的计算效率，增强了变网格技术的适应性和实用性。