技术领域
[0001] 本
发明属于石油勘探领域，具体涉及一种基于非规则扇形网格剖分的起伏隧道空间正演模拟系统及方法
背景技术
[0002] 石油是保障国家安全的重要战略资源，随着我国石油供需对外依存度逐年快速升高，亟需加强
地震勘探理论研究以提高我国油气储量及产量。我国主要油气探区地下构造复杂且储层岩性差异较大，为
地震勘探带来了巨大挑战，主要困难表现如下：
地震波传播规律复杂，地震有效信息微弱，高
精度地震成像困难。针对我国探区地质构造的复杂性，需要系统地开展针对复杂介质的地震波传播规律及高精度成像理论研究。
[0003] 随着油气勘探开发的深入，所面临地质构造越来复杂，为地震勘探带来了巨大挑战。针对我国探区地质构造的复杂性，需要系统开展介震波传播规律及高精度成像理论与方法的研究，从而更好地为油气藏的识别、描述挖潜等提供更加重要的地球物理依据。
[0004] 在隧道、
测井、钻孔等特殊地质环境中，传统的直
角坐标基于矩形网格的弹性波方程不能准确地模拟地震波场。
发明内容
[0005] 针对
现有技术中存在的上述技术问题，本发明提出了一种基于非规则扇形网格剖分的起伏隧道空间正演模拟系统，以及基于该正演模拟系统的模拟方法，该正演模拟系统及方法设计合理，克服了现有技术的不足，具有良好的效果。
[0006] 本发明所采用的技术解决方案是：
[0007] 一种基于非规则扇形网格剖分的起伏隧道空间正演模拟系统，包括顺序连接的输入模
块、网格生成模块、坐标变换模块、波场延拓模块、吸收边界模块、坐标反变换模块和输出模块；
[0008] 输入模块，被配置为用于输入起伏隧道空间高程函数、
震源及
检波器的大地坐标，隧道空间模型的纵横
波速度和
密度模型；
[0009] 网格生成模块，被配置为用于起伏隧道空间模型进行非规则扇形网格剖分；
[0010] 坐标变换模块，被配置为用于将起伏隧道空间模型映射到曲扇
坐标系下的矩形网格剖分的模型；
[0011] 波场延拓模块，被配置为用于推导曲扇坐标系弹性波一阶速度-应
力方程，利用曲扇坐标系下的全交错网格机制进
行波场更新；
[0012] 吸收边界模块，被配置为用于曲扇坐标系下三方向吸收边界条件吸收人工边界反射；
[0013] 坐标反变换模块，被配置为用于将波场快照变换到笛卡尔坐标系下；
[0014] 输出模块，输出起伏隧道空间模型的炮记录和波场快照。
[0015] 一种基于非规则扇形网格剖分的起伏隧道空间正演模拟方法，采用上述的基于非规则扇形网格剖分的起伏隧道空间正演模拟系统，包括以下步骤：
[0016] 步骤1：通过输入模块，输入起伏隧道空间高程函数、震源及检波器的大地坐标、隧道空间模型的纵横波速度和密度模型；
[0017] 步骤2：通过网格生成模块，对起伏隧道空间模型进行非规则扇形网格剖分；
[0018] 步骤3：通过坐标变换模块，将起伏隧道空间模型映射到曲扇坐标系下的矩形网格剖分的模型，通过如下所示的坐标变换公式进行变换：
[0019]
[0020] 其中，θ表示极坐标系下的角度变量，r表示极坐标系下的半径变量，和γ分别表示曲扇坐标系下的角度变量和半径变量； 和 分别定义 处扇形坐标系和曲扇坐标系下的最大半径；
[0021] 步骤4：通过波场延拓模块，推导曲扇坐标系弹性波一阶速度-
应力方程，曲扇坐标系下的方程可表示为：
[0022]
[0023] 其中，t为时间；ρ为密度，λ和μ为
拉梅常数；vθ和vr分别为角度分量和半径分量的速度；τrr，τθθ，τrθ为扇形坐标系下辅助应力场；γ0和r0分别表示扇形坐标系和曲扇坐标系下的最大半径；偏导数 和 可由下式求得：
[0024]
[0025] 步骤5：利用曲扇坐标系下的全交错网格机制进行波场更新，在二阶上使用二阶精度O(Δt2)，空间上使用高阶精度O(Δθ2M,Δγ2M)的有限差分格式来求解方程(2)，其中M表示空间精度；用f代表任意波场值，那么f的时间偏导数和两分量空间偏导数可以写为：
[0026]
[0027] 其中， 表示 方向的空间导数，Dγ表示γ方向的空间导数，Δt为时间步长，和Δγ分别为 方向和γ方向的网格步长，am表示有限差分算子，上标n表示时间坐标，下标i，j表示空间坐标；
[0028] 为了提高曲扇坐标系下的正演模拟的
稳定性，使用曲扇坐标系下的全交错网格机制进行数值离散；
[0029] 波场更新的递推格式为：
[0030]
[0031] 其中，
[0032]
[0033] 步骤6：通过吸收边界模块，曲扇坐标系下三方向吸收边界条件吸收人工边界反射，曲扇坐标系下包含三个边界：外边界，内边界和角边界；
[0034] 对于外边界，采用曲扇坐标系下的半径方向的复频移完美匹配层进行人工边界反射的压制；对于内边界，在隧道空间的起伏地表处：
[0035]
[0036] 将方程(7)带入方程(2)可得起伏自由地表边界条件：
[0037]
[0038] 对于角边界存在两种情况，第一种情况为封闭模型(θ∈[0,2π])，则边界条件满足：
[0039]
[0040] 第二种情况为非封闭模型(θ∈[0,θmax],θmax＜2π)，采用曲扇坐标系下的角度方向的复频移完美匹配层进行人工边界反射的压制；
[0041] 步骤7：通过坐标反变换模块，将波场快照变换到笛卡尔坐标系下；
[0042] 步骤8：通过输出模块，输出起伏隧道空间模型的炮记录和波场快照。
[0043] 优选的，在步骤4中，笛卡尔坐标系下二维弹性速度应力
波动方程为：
[0044]
[0045] 其中，vx和vz分别为
水平分量和垂直分量的速度，τxx，τzz，τxz为应力场，t为时间；ρ为密度，λ和μ为拉梅常数；在一些特殊的地质环境如隧道空间，传统笛卡尔坐标系下的正演模拟方法无法对地震波进行准确模拟；为了克服这一难题，首先将方程(10)变换到扇形坐标系下，在扇形坐标系中的弹性波波动方程为：
[0046]
[0047] 其中，vθ和vr分别为角度分量和半径分量的速度；τrr，τθθ，τrθ为扇形坐标系下辅助应力场；应用链式法则，曲扇坐标系下的方程(11)可表示为：
[0048]
[0049] 其中，偏导数 和 可由下式求得：
[0050]
[0051] 本发明的有益技术效果是：
[0052] 相比于传统正演模拟方法，本发明能够更准确地模拟某些特殊地质环境下的地震波场传播，如隧道、测井和钻孔空间。为了处理包含不规则表面形貌的隧道空间，本发明开发了一种过渡坐标系来改进扇形坐标下有限差分正演模拟方法的精度。在该方法中，首先将包含不规则表面形貌的隧道空间划分为非规则扇形网格，然后在曲扇坐标下波场计算，其中起伏地表被映射到一个平面上；采用曲扇坐标系下全交错网格数值离散格式求解曲扇坐标系中的波动方程。此外，曲扇坐标系下复频移完全匹配层、柱状自由表面和圆周边界条件用来吸收曲扇坐标周围不需要的人工反射。
[0053] 本发明的基于非规则扇形网格剖分的起伏隧道空间正演模拟方法为存在剧烈起伏地表的隧道空间提供准确的正演模拟波场，有助于分析地震波在该探区的传播规律，更好地为起伏隧道空间下油气藏的识别、描述挖潜等提供更加重要的地球物理依据。
附图说明
[0054] 图1为本发明的基于非规则扇形网格剖分的起伏隧道空间正演模拟方法的
流程图；
[0055] 图2为曲扇坐标系下全交错网格示意图；
[0056] 图3示出曲扇坐标系下三方向吸收边界条件；
[0057] 图4示出不规则隧道空间模型；其中(a)笛卡尔坐标系；(b)曲坐标系；(c)扇形坐标系；(d)曲扇坐标系；
[0058] 图5为采用本发明方法得到的曲扇坐标系下的波场快照；其中(a)，(c)，(e)示出角度方向；(b),(d),(f)示出半径方向；(a)，(b)为500ms；(c)，(d)为1000ms；(e),(f)为1500ms；
[0059] 图6为采用本发明方法得到的笛卡尔坐标系下的波场快照；其中(a)，(c)，(e)示出角度方向；(b),(d),(f)示出半径方向；(a)，(b)为500ms；(c)，(d)为1000ms；(e),(f)为1500ms；
[0060] 图7为采用本发明方法得到的炮记录；(a)为角度方向；(b)为半径方向；
[0061] 图8为采用常规方法得到的曲扇坐标系下的波场快照；其中(a)，(c)，(e)示出角度方向；(b),(d),(f)示出半径方向；(a)，(b)为500ms；(c)，(d)为1000ms；(e),(f)为1500ms；
[0062] 图9为采用常规方法得到的曲扇坐标系下的波场快照；其中(a)，(c)，(e)示出角度方向；(b),(d),(f)示出半径方向；(a)，(b)为500ms；(c)，(d)为1000ms；(e),(f)为1500ms；
[0063] 图10为本发明中基于非规则扇形网格剖分的起伏隧道空间正演模拟系统的结构示意图。
具体实施方式
[0064] 在隧道、测井、钻孔等特殊地质环境中，传统的直角坐标基于矩形网格的弹性波方程不能准确地模拟地震波场。基于上述技术问题，本发明提出一种通过过渡坐标系来改进扇形坐标下有限差分正演模拟的系统及方法，以能够适应复杂隧道等特殊地质环境。
[0065] 本发明基于非规则扇形网格剖分的起伏隧道空间正演模拟方法，首先将包含不规则表面形貌的隧道空间划分为非规则扇形网格，然后波场计算是在曲扇坐标下计算的，其中起伏地表被映射到一个平面上。采用曲扇坐标系下全交错网格数值离散格式求解曲扇坐标系中的波动方程。此外，曲扇坐标系下复频移完全匹配层、柱状自由表面和圆周边界条件用来吸收曲扇坐标周围不需要的人工反射。
[0066] 本发明为存在剧烈起伏地表的隧道空间提供了准确的正演模拟波场，并分析地震波在该探区的传播规律，更好地为起伏隧道空间下油气藏的识别、描述挖潜等提供更加重要的地球物理依据。
[0067] 下面结合附图及具体
实施例对本发明进行更为具体的说明。
[0068] 实施例1
[0069] 一种基于非规则扇形网格剖分的起伏隧道空间正演模拟系统，其结构如图10所示，包括输入模块、网格生成模块、坐标变换模块、波场延拓模块、吸收边界模块、坐标反变换模块和输出模块；以上模块顺序连接；
[0070] 输入模块，被配置为用于输入起伏隧道空间高程函数、震源及检波器的大地坐标、隧道空间模型的纵横波速度和密度模型；
[0071] 网格生成模块，被配置为用于起伏隧道空间模型进行非规则扇形网格剖分；
[0072] 坐标变换模块，被配置为用于将起伏隧道空间模型映射到曲扇坐标系下的矩形网格剖分的模型；
[0073] 波场延拓模块，被配置为用于推导曲扇坐标系弹性波一阶速度-应力方程，利用曲扇坐标系下的全交错网格机制进行波场更新；
[0074] 吸收边界模块，被配置为用于曲扇坐标系下三方向吸收边界条件吸收人工边界反射；
[0075] 坐标反变换模块，被配置为用于将波场快照变换到笛卡尔坐标系下；
[0076] 输出模块，输出起伏隧道空间模型的炮记录和波场快照。
[0077] 实施例2：
[0078] 在上述实施例的
基础上，本发明还提到一种基于非规则扇形网格剖分的起伏隧道空间正演模拟方法，其流程如图1所示，具体包括如下步骤：
[0079] 步骤1：输入起伏隧道空间高程函数、震源及检波器的大地坐标、隧道空间模型的纵横波速度和密度模型。
[0080] 步骤2：起伏隧道空间模型进行非规则扇形网格剖分。
[0081] 步骤3：将起伏隧道空间模型映射到曲扇坐标系下的矩形网格剖分的模型，通过如下所示的坐标变换公式进行变换：
[0082]
[0083] 其中，θ表示极坐标系下的角度变量，r表示极坐标系下的半径变量，和γ分别表示曲扇坐标系下的角度变量和半径变量。 和 分别定义扇形坐标系和曲扇坐标系下的最大半径。
[0084] 步骤4：推导曲扇坐标系弹性波一阶速度-应力方程，笛卡尔坐标系下二维弹性速度应力波动方程为
[0085]
[0086] 其中，vx和vz分别为水平分量和垂直分量的速度，τxx，τzz，τxz为应力场，t为时间；ρ为密度，λ和μ为拉梅常数。在一些特殊的地质环境如隧道空间，传统笛卡尔坐标系下的正演模拟方法无法对地震波进行准确模拟。为了克服这一难题，首先将方程(2)变换到扇形坐标系下，在扇形坐标系中的弹性波波动方程为：
[0087]
[0088] 其中，vθ和vr分别为角度分量和半径分量的速度；τrr，τθθ，τrθ为扇形坐标系下辅助应力场。应用链式法则，曲扇坐标系下的方程(3)可表示为：
[0089]
[0090] 其中，偏导数 和 可由下式求得：
[0091]
[0092] 步骤5：利用曲扇坐标系下的全交错网格机制进行波场更新，本发明在二阶上使用二阶精度O(Δt2)，空间上使用高阶精度O(Δθ2M,Δγ2M)的有限差分格式来求解方程(4)，其中M表示空间精度。用f代表任意波场值，那么f的时间偏导数和两分量空间偏导数可以写为：
[0093]
[0094] 其中Δt为时间步长， 和Δγ分别为 方向和γ方向的网格步长，am表示有限差分算子，上标n表示时间坐标，下标i，j表示空间坐标。
[0095] 为了提高曲扇坐标系下的正演模拟的稳定性，本发明使用曲扇坐标系下的全交错网格机制进行数值离散(图2)。
[0096] 波场更新的递推格式为：
[0097]
[0098] 其中，
[0099]
[0100] 步骤6：曲扇坐标系下三方向吸收边界条件吸收人工边界反射，曲扇坐标系下包含三个边界：外边界，内边界和角边界(图3)。
[0101] 对于外边界，本发明采用曲扇坐标系下的半径方向的复频移完美匹配层进行人工边界反射的压制；对于内边界，在隧道空间的起伏地表处：
[0102]
[0103] 将方程(9)带入方程(4)可得起伏自由地表边界条件:
[0104]
[0105] 对于角边界存在两种情况。第一种情况为封闭模型(θ∈[0,2π])，则边界条件满足：
[0106]
[0107] 第二种情况为非封闭模型(θ∈[0,θmax],θmax＜2π)，采用曲扇坐标系下的角度方向的复频移完美匹配层进行人工边界反射的压制。
[0108] 步骤7：将波场快照变换到笛卡尔坐标系下。
[0109] 步骤8：输出起伏隧道空间模型的炮记录和波场快照。
[0110] 本发明基于非规则扇形网格剖分的起伏隧道空间正演模拟方法，为存在剧烈起伏地表的隧道空间提供准确的正演模拟波场，有助于分析地震波在该探区的传播规律，更好地为起伏隧道空间下油气藏的识别、描述挖潜等提供更加重要的地球物理依据。
[0111] 应用实验
[0112] 本发明基于非规则扇形网格剖分的起伏隧道空间正演模拟方法，应用于不规则隧道空间模型数据，取得了理想的计算效果。
[0113] 输入起伏隧道空间高程函数、震源及检波器的大地坐标、隧道空间模型的纵横波速度和密度，不规则隧道空间模型的纵波速度如图4(a)所示，纵横波速度满足泊松体，密度模型通过ρ＝0.31vp0.25得到；起伏隧道空间模型进行非规则扇形网格剖分，网格剖分如图4(a)中实线所示；将起伏隧道空间模型映射到曲扇坐标系下的矩形网格剖分的模型，图4(b)和4(c)分别为曲坐标系和扇形坐标系下的模型，图4(d)为最终曲扇坐标系下的模型；推导曲扇坐标系弹性波一阶速度-应力方程；利用曲扇坐标系下的全交错网格机制(图1)进行波场更新；曲扇坐标系下三方向吸收边界条件(图2)吸收人工边界反射；得到的曲扇坐标系下的波场快照如图5所示，将波场快照变换到笛卡尔坐标系下，如图6所示，输出起伏隧道空间模型的炮记录(如图7所示)和波场快照。
[0114] 从图5、6、7可以看出，采用本发明提出的基于非规则扇形网格剖分的起伏隧道空间正演模拟方法能够准确模拟地震波在起伏隧道空间中的波场传播特征。剧烈起伏地表的影响得到了消除，直达波、反射波、转换波及面波都能很好地模拟出来。作为对比，本发明还采用常规曲网格正演模拟方法对该起伏隧道空间模型进行模拟。得到的曲网格坐标系和笛卡尔坐标系下的波场快照如图8和图9所示。相比于本发明的方法，传统方法产生了异常噪音，如图中箭头所示。因此该结果证明了本发明能够得到更加准确、实用的正演模拟结果。
[0115] 上述方式中未述及的部分采取或借鉴已有技术即可实现。
[0116] 当然，上述说明并非是对本发明的限制，本发明也并不仅限于上述举例，
本技术领域的技术人员在本发明的实质范围内所做出的变化、改型、添加或替换，也应属于本发明的保护范围。
