一种面向隧道探测的快速地震波走时计算方法

本发明涉及地震波走时计算，具体涉及一种面向隧道探测的快速地震波走时计算方法。

背景技术：

1、反射地震波成像方法是主要的隧道超前探测方法，地震波走时是其核心参数，常被广泛的应用于偏移与层析成像等技术中。现有的走时计算方法基本都是基于规则模型提出的，而现实中处于施工期的隧道常常拥有不规则的界面，此时要求地震波走时算法不仅要对复杂地表有良好的适应性而且在处理不规则边界时仍然有较高的计算效率，因此，研究面向隧道探测的快速地震波走时计算方法具有重要意义。

2、《勘探地球物进展》2007年第30卷第5期公开了孙章庆等“波前快速推进法起伏地表地震波走时计算”介绍了波前快速推进法(fmm)是一种准确而稳定的地震波走时计算方法，并且以fmm为基础，提出了一种复杂地表条件的计算地震波走时的计算方法，计算取得了良好的效果，其算法的计算复杂度为o(nlog2n),n表示计算网格点总数。

3、中国发明专利申请201110397445.4中公开了《vti介质中地震波的走时计算方法》，该算法通过输入地下地质速度模型以及观测系统参数，然后计算相应中间结果，主要采用多项式求解的方式，获取地震波走时，通过不同模型的验证计算结果取得了良好的效果。

4、中国发明专利申请201810077621.8中公开了《一种混合二维地震波走时计算方法》，该算法通过读入速度模型与相关计算参数，从炮点沿不同方向追踪一定距离的射线信息，利用波前构建法计算射线范围内地震波走时，然后利用快速推进法计算剩余网格节点的走时，通过数值模拟证明了该方法良好的计算效果。

5、通过以上的例子可以看出，现有的地震波走时计算方法针对常规水平模型具有良好的计算效果，但是这些方法在面对复杂起伏地表问题时为了获取较高的计算精度，常需要较长的计算时间，或者为了提高计算效率而相应降低了算法的计算精度。

技术实现思路

1、本发明要解决的技术问题是提供一种面向隧道探测的快速地震波走时计算方法，过采用迎风差分方案离散程函方程梯度项，求解程函方程从而获取地震波走时，基于快速推进法的窄带技术，采用一种更快速且更接近实际地震波前传播方式，即在地震波前扩展时每次在窄带中选取一组点作为扩展点，在改善原有走时计算方法的精度的同时，大大提高了计算效率。

2、为解决上述技术问题，本发明采用如下技术方案：

3、设计一种面向隧道探测的快速地震波走时计算方法，包括以下步骤：

4、s1、读入地质模型以及相关参数信息，所述参数包含速度模型的网格点数，网格间距，震源位置和起伏地表信息。

5、s2、对网格点进行属性划分：将所有网格点分为接收点、窄带点和远离点；接收点指所有网格点中被认为已经走时计算的点；窄带点是指所有窄带内的点；远离点是指除接收点和窄带点以外的所有网格节点，计算过程中用tt表示该计算点的走时值。

6、s3、震源初始化、窄带初始化和最小走时点选取；

7、s31、以震源点为起始点将震源点属性设置为idtt＝2，tt＝0，选取震源点的邻点，即上下左右四个点，判断其属性，idtt不为2，则求取其解析值，公式如下：

8、ti±1,j＝si±1,jh,ti,j±1＝si,j±1h

9、其中，h表示网格间距，s表示该点慢度；t表示该点的地震波走时值，i表示该点的横坐标，j表示该点的纵坐标；

10、s32、将震源临近网格点属性设置为idtt＝1，并将其纳入窄带构成初始窄带，同时将各邻近点的走时tt存入数组gt中；

11、s33、设置gt数组内最小走时值为tm。

12、s4、群组g选取：群组g表示在数组gt中走时属于(tm，tm+dt)范围内的点；群组g的所有点都作为窄带中的更新计算点，进行下一步更新计算。

13、所述步骤s4中群组g的选择标准为：

14、在二维情况下，设h＝δx＝δz,定义

15、

16、

17、其中，h表示网格间距，x、z分别表示横向与纵向的单位步长，w表示计算波前；表示波前w中的一点，表示慢度函数，表示走时函数；

18、则群组g的选择标准如下：

19、

20、其中，dt表示时间，h表示网格间距，tw,min表示计算波前w中的最小走时点，sw,min表示计算波前w中最小走时点对应的慢度值。这样确定时间的数值，是为了保证，群组g中的元素个数大于1。

21、s5、gmm采用迎风差分公式离散程函方程的走时梯度项，来求取地震波走时。

22、所述步骤s5中，求取地震波走时的过程为：

23、网格节点走时是通过迎风差分法求解程函方程获得的，在二维情况下，地震波传播波前满足程函方程：

24、

25、其中t(x,z)为走时函数，s(x,z)为模型介质慢度函数；

26、gmm采用迎风差分公式离散程函方程的走时梯度项，来求取地震波走时；

27、

28、其中分别为x和z方向前向和后向差分算子，以一阶差分格式为例，有

29、

30、

31、s6、窄带扩展演化；

32、通过逆序遍历数组gt和顺序遍历数组gt，根据数组中的每个点(i,j)的走时，对其所有临近点(x,z)进行属性判断，更新该邻近点(x,z)的走时，将该临近点(x,z)纳入或移出gt数组，直至gt数组为空，输出地震波走时计算结果。

33、所述步骤s6中窄带扩展演化的具体过程为：

34、s61、设置

35、其中，dt表示时间，h表示网格间距，sw,min表示计算波前w中最小走时点对应的慢度值；

36、s62、设置tm＝tm+dt；

37、s63、以逆序遍历数组gt，对于数组中的每个点(i,j)如果其走时满足tt(i,j)＜＝tm，则对其所有临近点(x,z)进行属性判断，若idtt<＝1，则更新该邻近点(x,z)的走时；

38、s64、以顺序遍历数组gt，对于数组gt中的每个点(i,j)如果其走时满足tt(i,j)＜＝tm，则

39、(ⅰ)对其所有临近点(x,z)进行属性判断，若idtt<＝1，则更新该邻近点(x,z)的走时；

40、(ⅱ)若其临近点(x,z)属性为idtt＝0，则将其属性改为idtt＝1，并将其临近点(x,z)的走时纳入gt数组中；

41、(ⅲ)将该点(i,j)的属性设置为idtt＝2，即将该点(i,j)移出数组gt；

42、s7、判断数组gt是否为空，若为空则输出地震波走时计算结果，若不为空则返回s6继续计算。

43、本发明的有益效果在于：本发明提供的一种面向隧道探测的快速地震波走时计算方法，通过利用迎风有限差分方案求解程函方程来获取地震波走时，并且基于窄带技术，在地震波前扩展时选择一组点作为扩展点同时进行扩展计算，在保证了计算精度的同时，大大提高了地震波走时算法的计算效率。

44、隧道探测的环境通常比较复杂，对方法的稳定性要求较高，本发明具有无条件稳定性，有很强的适应能力。另外隧道探测对探测结果的时效性要求比较高，本方法具有很高的计算效率，可以提高隧道地震勘探处理的计算效率。