三维地表下地震波走分法

开篇：润墨网以专业的文秘视角，为您筛选了一篇三维地表下地震波走分法范文，如需获取更多写作素材，在线客服老师一对一协助。欢迎您的阅读与分享！

1引言

地震数据处理技术效果不好是起伏地表地区地震勘探经常会面临的一个挑战，其主要原因是现有的地震数据处理技术很多都基于水平地表假设，所以研究专门针对起伏地表的三维地震数据处理技术是很有必要的［1］．地震波走时计算技术是很多地震数据处理技术［2－7］中必不可少的正演模拟工具，因此对起伏地表条件下的三维地震波走时计算技术进行研究，可以为研发能很好适应起伏地表条件的三维地震数据处理技术打下基础，进而有助于提高起伏地表地区的三维地震勘探的效果．计算地震波走时的方法主要有两大类：射线追踪和数值求解Eikonal方程．其中，射线追踪方法主要包括：两点射线追踪法［8］、波前构建法［9－10］、最短路径射线追踪法［11－12］、插值法［13－14］、辛几何算法［15－16］等；数值求解Eikonal方程的方法主要是指有限差分法［17－22］．关于三维起伏地表条件下的地震波走时计算，国内外学者主要采用射线追踪类方法（包括：两点射线追踪法［23］、插值法［24］和最短路径射线追踪法［25］等）直接或间接地研究过该问题，并取得了一些研究成果．实际上，由于各种差分格式和各种波前扩展方式的提出和发展，数值求解Eikonal方程的方法已具有计算效率高、精度易调控（可以通过差分格式阶数的改变来调控算法的精度）、适用于复杂介质模型、稳定性好、易于编程实现等优点［18－22］．综合利用有限差分法所具备的优点，笔者曾基于对该方法的改进求解过二维起伏地表条件下的地震波走时计算问题［26－27］．与二维空间中的算法相比，有限差分法在处理三维空间中的起伏地表问题时，必然会遇到一些新的问题，但算法的基本原理决定了其具有处理三维空间问题的能力．因此，本文提出了三维空间中算法的具体内容：①三维空间中起伏地表模型的不等距网格建立方法；②综合运用不等距差分格式、迎风差分格式、Huygens原理和Fermat原理推导建立三维空间中适应于不等距网格的局部走时计算公式；③基于对常规窄带技术的改进来获得三维起伏地表条件下算法的整体实现步骤；④对方法进行精度分析，并通过引入网格加密技术在不过多牺牲计算效率的情况下实现计算精度的大幅度提高；⑤通过一个计算实例检验算法在处理三维起伏地表、地下复杂地质构造时的有效性和适应性．

2三维起伏地表模型的网格剖分

为了准确地描述三维起伏地表的形态，并为后续数值计算打下基础，采用不等距网格剖分三维起伏地表模型．如图1所示，不等距网格在进行三维起伏地表模型剖分时，在地表附近的区域采用纵向上网格间距不等的网格，在除地表附近以外的空间采用规则的立方体网格．不等距网格剖分三维起伏地表模型具有如下特点：①它能准确定位地表上采样点的高程位置，而不用作“利用就近的规则网格节点代替地表采样点”的近似处理；②与地表附近的四面体网格和其他网格形态相比，不等距网格相对简洁，这也会让地表附近的局部走时算法变得相对简洁3局部走时计算公式的推导为了让有限差分法能够适应如图1所示的不等距网格，首先需要建立该网格中的局部走时计算公式，当然该局部走时计算公式首先必须满足地震波的传播规律，同时还需要在面对不等距网格时具有无条件稳定性．综合分析有限差分法中各种差分格式后发现迎风差分格式具有无条件稳定性［18］，但是该差分格式不能直接应用于不等距网格中．为了解决该问题，本文通过在迎风差分格式中融入不等距差分格式、Huygens原理和Fermat原理来建立三维起伏地表条件下的局部走时计算公式．其中迎风差分格式的思想用于控制算法的稳定性，不等距差分格式用于适应不等距网格，而Huygens原理和Fermat原理则用于保证算法在地震波传播规律的控制下进行．下面首先简要地回顾一下常规迎风差分格式的基本原理．为了便于阐述，首先引入三维空间中的Eikonal方程：其中，max（A，B，C）为求取A，B，C中的最大值的函数，D－lnti，j，k、D＋lnti，j，k（l＝x、y、z）分别为空间中一点（i，j，k）处l方向上的函数T的n阶向后、向前差分算子（为了表述简单，此处和后续的ti，j，k均用t表示，max（D－l1ti，j，k，－D＋l1ti，j，k，0）简写为max（l）），如图3所示，根据地表附近的被算点在不等距网格中所处位置的不同，可将处于不规则网格中的节点类型分为：①位于地表表面上的网格节点（如图3中的点H）；②与地表表面上的网格节点直接相邻的网格节点（该节点到地表的距离小于一个网格间距，即如图3中的M点）．下面将分别推导这两类网式中：θ＝h2／d2．公式（8）即为计算M点的走时值时的局部走时计算公式，该公式是通过在迎风差分格式中融入不等距差分格式，并综合应用Fermat原理而获得的．②当H点为当前被算点时．设点H1、H2、H3、H4的走时值分别为tH1、tH2、tH3、tH4．此时H点处的局部走时公式的建立问题将变得更加复杂，因为在该点处无法建立迎风差分格式，同时H点在X、Y两个方向上的网格节点均为地表上的点（如图3所示的点H1、H2、H3、H4），并且绝大多数情况下与H点不在同一高程上．为了处理这些问题，首先，同样在M点处将不等距差分格式融入迎风差分格式中得公式（6），不过此时的未知量变为了tH，而解此以tH为未知量的一元二次方程有：公式（11）即为计算H点的走时值时的局部走时计算公式，该公式是通过在迎风差分格式中融入不等距差分格式，并综合应用Fermat原理和Huygens原理而获得的．综合上述推导可知，三维起伏地表条件下走时计算的局部计算公式，因为被算点所处的位置不同而不同．当计算规则网格节点的走时值时，采用公式（5）；当计算地表附近的不规则网格中的节点的走时值时则采用公式（8）和（11），其中公式（8）用于计算地表表面上的网格节点的走时值，而公式（11）则用于计算与地表表面上节点直接相邻的网格节点的走时值．

4算法的整体实现步骤

以上分别建立了三维起伏地表模型和推导了局部走时计算公式，除此之外还需要提出算法的整个实现步骤．不同于水平地表问题，这里算法的整体实现步骤需要在保证满足地震波传播规律（保证算法的合理性，即算法的实现策略不破坏地震传播过程的因果关系）的同时，还能处理起伏地表条件造成的不规则计算边界的问题．综合分析该问题，并充分结合笔者在处理二维空间中的类似问题时的研究经验［26－27］可知，作为有限差分法中的一种稳定、高效且灵活的波前扩展方式，窄带技术完全具有处理三维起伏地表问题的潜力［18－21］．因此下面将以窄带技术的基本思想为基础，以设定新的网格节点类型的方式来建立三维起伏地表条件下的算法的整体实现步骤．如图4a所示，常规窄带技术在实现时把整个计算区域的所有网格节点分为三种类型，包括：①完成点，表示走时计算已经完成的节点，原文为“Accepted”［18］在此意译为“完成点”，即图4a中黑色填充的圆点；②窄带点，表示近似波前上的节点，即图4a中灰色填充的圆点；③远离点，表示走时计算还未进行的节点，即图4a中白色填充的圆点．窄带技术的实现分为初始化和循环计算两个部分．初始化时，仅有震源点被设定为完成点；与震源点直接邻近的几个网格节点被设定为窄带点，它们构成初始窄带（即初始前）；除此之外的那些网格节点均被设定为远离点循环计算包括如下步骤：步骤1：通过比较找出窄带内走时值最小的带点作为子震源点，并将其属性重新设定为完成点步骤2：对子震源直接邻近的6个网格节点类型进行判断，如果其类型为远离点则采用公式（5计算它的走时值，并将其类型重新设定为窄带点；果其类型为窄带点则再次计算它的走时值，并与原来的走时值比较，较小者被保留；如果其类型为成点则不对其作任何处理．步骤3：判断窄带点是否还存在，是则跳回步1重复循环计算，否则计算结束．

基于以上阐述可知，常规窄带技术通过近似述波前的窄带的不断推进，来模拟波前在介质中传播过程．其中窄带技术实现过程中各个网格节的类型的不断交替变化是窄带技术的灵魂，因为是这样的交替变化才实现了窄带的不断推进．很然这样的网格类型的不断交替变化在起伏地表条下也是可以实现的，当然必须对其进行一些必要改进．改进的目的是保证窄带的不断推进只在地以下进行（因为地震波不能穿越地空界面进入空中），同时还能根据节点的位置对应地选取计算公（因为常规的计算公式在地表附近会失效）．为了达到上述目的，在此采用的策略是在窄技术中设定两种新类型的网格节点，即如图4b的表点（图4b中起伏地表表面的点）和地表上点（地以上的网格节点，图4b没有画出，但计算时由于用规则的数组，所以在计算时存在地表以上的点）．当然有别于其它原有的三种类型的网格节点，这里的地表点的属性在计算过程中可能分别为完成点、窄带点或远离点（如图4b所示地表上网格节点的填充颜色有三种可能情况，即黑色、灰色、白色）．与常规窄带技术相比，起伏地表条件下的窄带技术只需要对上述步骤2进行改进：步骤2：对子震源直接邻近的6个网格节点的类型进行判断，如果其类型为远离点则采用公式（5）、（8）或（11）计算它的走时值（公式的选择是通过判断该点是否为地表点或为与地表点直接相邻的节点来实现的），并将其类型重新设定为窄带点；如果其类型为窄带点则再次计算它的走时值，并与其原来的走时值比较，较小者被保留；如果其类型为完成点或地表上点则不对其作任何处理．通过改进步骤2而获得的起伏地表条件下的窄带技术是否能适应不规则计算边界并满足波前传播的规律呢？当然能，因为：①起伏地表条件下的窄带技术保留步骤1不变，即找出波前上的最小走时点作为子震源点的方法是与Huygens原理和Fermat原理相符的；②修改后的步骤2中加入“如果其类型为完成点或地表上点则不对其作任何处理”的操作是为了保持窄带不能扩展进入地表以上区域，这与地震波不能穿越地空界面进入空气中的物理事实是相符的．

5计算精度分析及实例

为了分析上述提出的算法的计算精度，在此采用如图5a所示的模型．模型的地表为一山峰，整个计算区域的大小为8．0km×8．0km×6．0km，地下介质的地震波速度为1．0km／s，计算时分别采用50．0m、20．0m的网格间距，震源点置于（4．0km，4．0km，6．0km）点处，并且为了使得地表上各个点的走时值可以通过解析公式计算而获得，特设计地表上的各点可以通过直线直接与震源点相连接．图5c1、图5c2分别给出了采用50．0m、20．0m为网格间距时，地表上各个网格节点的走时值的数值计算结果相对于解析值的相对误差．图5d1、图5d2分别给出了采用50．0m、20．0m为网格间距时，剖面Y＝0．0km内的走时值的数值计算结果相对于解析值的相对误差．这些相对误差的分布均体现出一个共同的特点：震源附近区域的相对误差最大，远离震源的区域相对误差小，并且离震源越远的区域误差越小．基于上述误差分布规律，为了在不过多牺牲计算效率的情况下大幅度提高算法的计算精度，在此引入一种震源附近的网格加密技术［28］．如图5b所示，该网格加密技术对震源附近一定区域内的网格单元依次进行2n×2n，2n－1×2n－1，…21×21的细分，这样的细分使得在实际计算时网格在震源附近区域采用很小的网格间距，而随着离震源距离的增加网格间距逐渐增加，直到网格间距变为常规网格间距．本文选择n＝5对常规网格间距为20．0m的情况进行网格加密处理（关于网格加密的基本原理和实现策略请读者详见参考文献［28］），图5c3、图5d3分别给出了地表上和Y＝0．0km剖面上各网格节点的相对误差，与图5c2、图5d2对比分析发现：算法的整体计算精度有大幅度提高，尤其是在震源附近区域．此外，网格加密后算法的整体计算量仅增加了3．71％．以上得出的精度分析结果表明：本文提出的算法能够很好地解决三维起伏地表条件下的地震波走时计算问题，并且算法能够获得很好的计算精度．与此同时网格加密技术能在不过多牺牲计算效率的前提下大幅度提高计算精度．为了检验本文算法能否稳定、灵活及有效地处理三维起伏地表问题，下面将给出一个三维起伏地表模型算例，该模型由Marmousi模型改造而获得，模型具有地表起伏剧烈、地下（包括近地表）地质构造复杂等特点．如图6a所示，模型的大小为9．0km×9．0km×6．0km，计算时将震源置于（0．0km，0．0km，5．0km）处，采用20．0m的网格间距．图6（b—d）分别给出了X＝0．0km、X－Y＝0．0km、Y＝0．0km剖面的等时线分布．分析图6（b—d）的计算结果可以得出：本文算法的计算结果能与地震波在起伏地表条件下的传播规律相符，同时该计算实例也充分表明本文提出的算法在处理三维起伏地表问题时，能够保证很好的稳定性、灵活性和有效性，并且算法能同时处理三维剧烈地表起伏、地下（包括近地表）复杂地质构造等问题．

6结论

三维起伏地表条件下的地震波走时计算技术是研究三维起伏地表地区很多地震数据处理技术的基础性工具．为了获得适应于三维起伏地表且精度高的走时算法，本文提出了一种不等距迎风差分法．算法通过在常规迎风差分格式中综合引入不等距差分格式、Huygens原理和Fermat原理来获取地表附近的局部走时计算公式，同时通过改进常规窄带技术来获取算法的整体实现步骤．针对算法的精度分析及算例分析表明本文提出的算法具有如下特点：①由于采用不等距网格剖分起伏地表模型，并采用有限差分法作为局部走时算法，所以方法简洁且易于编程实现；②方法的实现步骤能够满足地震波传播的基本规律，即Huygens原理、Fermat原理以及地震波在地空边界处的传播特点；③算法能够保证很高的计算精度，同时网格加密技术能在不过多牺牲计算效率的情况下大幅度提高计算精度；④算法能稳定、灵活及有效地处理剧烈起伏的三维地表、地下（包括近地表）复杂地质构造等问题．