1.本发明属于隧道地球物理探测技术领域，具体涉及一种地质预报探测断层二维数值模拟方法。
背景技术：2.我国每年新掘进的隧道进尺可以贯穿整个地球，总长超过12000km。隧道开挖前进行超前地质预报意义重大，对隧道超前地质预报进行研究十分必要。由于施工过程中，隧道周围围岩形状是起伏不定的，而传统的隧道超前地质预报探测有限差分数值模拟时，将隧道形状当成规则形状进行处理，如圆柱、矩形等，与实际条件是不符合的，直接影响数值模拟的结果，目前针对隧道起伏条件下的数值模拟方法有限，因此亟需开展关于隧道起伏地震超前地质预报探测数值模拟技术方法，进一步提高数值模拟精度。
技术实现要素：3.为了克服现有技术的上述缺点，本发明的目的是提供一种地质预报探测断层二维数值模拟方法，能够有效的提高隧道超前地质预报探测数值模拟的准确性。
4.本发明解决其技术问题所采用的技术方案是：
5.一种地质预报探测断层二维数值模拟方法，包括以下步骤：
6.步骤一：测量隧道掌子面反方向的隧道的起伏四周的起伏情况，构建起伏隧道模型；
7.步骤二：建立均匀介质速度的三维模型，速度参数为v
p1
,v
s1
,ρ1，其中v
p1
,v
s1
,ρ1分别为三维模型的纵波速度、横波速度、密度；
8.步骤三：将隧道高程参数信息插入速度模型，建立起伏隧道的速度模型，隧道模型速度参数为v
p2
,v
s2
,ρ2，其中v
p2
,v
s2
,ρ2分别为隧道空间的纵波速度、横波速度、密度；
9.步骤四：在x方向处设置一个走向90度的断层，则断层部分的速度参数为v
p3
,v
s3
,ρ3，其中v
p3
,v
s3
,ρ3分别为断层部分的纵波速度、横波速度、密度；
10.步骤五：沿隧道建立二维坐标系，x方向为指向掌子面，y方向为垂直于隧道侧壁，以模型边界上隧道面的中心为原点建立二维坐标系，建立隧道开挖超前地质预报观测系统；
11.步骤六：给定网格间距和时间步长等参数，需要满足采样定律和稳定性条件；
12.步骤七：通过真实高程计算真实地形点和隧道表面的网格点的距离；
13.步骤八：在隧道处采用不等距有限差分格式，其余部分采用常规有限差分格式；
14.步骤九：在起伏隧道的侧壁布置检波器，从掌子面往隧道方向进行依次均布n个三分量检波器和m个震源，其中三分量检波器号设为rn，震源号为sm，构建基底地震超前地质预报观测系统；
15.步骤十：采用雷克子波、震源频率为z hz，采样间隔为l，采样时间为t；
16.步骤十一：采用常规的pml吸收边界进行处理；
17.步骤十二：迭代速度和应力，输出波场快照和地震记录。
18.优选的，在所述步骤一中：
19.(1)根据隧道情况，从掌子面往隧道方向测量90m范围内的整个隧道高程参数；
20.(2)将高程参数转换成相对高程，构建整个隧道的基本形态；
21.(3)通过隧道的基本形态，构建起伏隧道模型。
22.优选的，在所述步骤二中：所述三维模型中，x方向、y方向长度分别为300m、100m；
23.在步骤三中：v
p2
＝340m/s,v
s2
＝0m/s,ρ2＝1.29kg/m2。
24.优选的，在所述步骤六中：
25.(1)利用波动方程对网格进行网格化，大小为h；
26.(2)设置数值模拟参数：频率为f，单位为hz，采样间隔为l，单位为ms；
27.(3)为避免数值频散，要满足采样定理和稳定性条件，具体为网格大小和时间步长：
[0028]ⅰ.网格间距大小满足：
[0029]
其中，v
min
/(2fn)代表模拟的最小波长，v
min
为数值模拟参数中最小的速度参数，n为所取的差分阶数相关因子；
[0030]ⅱ.时间步长大小满足：
[0031]
其中，m是与所取差分阶数相关的因子，v
max
为为数值模拟参数中最小的速度参数。
[0032]
优选的，在所述步骤八中：
[0033]
(1)对模型进行网格化后，在起伏空气隧道侧壁外围采用常规等距离有限差分方法处理；
[0034]
(2)在起伏隧道附近进行不等距差分，主要分为应力和速度的不等距差分，具体如下：
[0035]ⅰ.通过u'＝r
β
u将起伏隧道坐标旋转为水平条件，其中u＝(v
x
,vy,σ
xx
,σ
xy
,σ
yy
)
t
,u'＝(v'
x
,v'y,σ'
xx
,σ'
xy
,σ'
yy
)
t
，v
x
,vy分别代表的速度分量，σ
xx
,σ
xy
,σ
yy
为应力分量，v'
x
,v'y为转换后的速度分量，σ'
xx
,σ'
xy
,σ'
yy
为转换后的应力分量，r
α
为坐标旋转矩阵，β为旋转角度，具体为：
[0036][0037]ⅱ.计算应力时，不等距差分共有8种情况：
[0038]
①
起伏条件下σ
xx
、σ
xy
、σ
yy
距离隧道起伏位置四周均有四个网格点时，则不存在不等距差分，则σ
xx
、σ
xy
、σ
yy
表达式如下：
[0039][0040]
其中，δl、δh分别是采样间隔和网格大小，i,j分别代表x方向和y方向的网格点数，λ、μ为拉梅参数；
[0041]
②
起伏条件下σ
xx
、σ
yy
距离隧道起伏位置四周均有四个网格点，但σ
xy
距离隧道起伏位置四周只有三个网格点时，计算σ
xx
、σ
yy
采用常规差分格式计算；而计算σ
xy
时存在不等距差分计算，则σ
xx
、σ
xy
、σ
yy
表达式如下：
[0042][0043]
其中，y0(i)为隧道起伏的坐标参数，下述y1(i)为隧道坐标向下取整；
[0044]
③
起伏条件下σ
xx
、σ
yy
距离隧道起伏位置四周均有三个网格点，σ
xy
距离隧道起伏位置四周只有两个网格点时，则差分格式均为不等距差分格式，则σ
xx
、σ
xy
、σ
yy
表达式如下：
[0045][0046]
④
起伏条件下σ
xx
、σ
yy
距离隧道起伏位置四周均有三个网格点，该分量上σ
xy
距离隧道起伏位置四周只有一个网格点时，差分格式改用不等距差分格式，则只需要计算σ
xx
、σ
yy
，表达式如下：
[0047][0048]
⑤
起伏条件下σ
xx
、σ
yy
距离隧道起伏位置四周均有三个网格点，该分量上σ
xy
距离隧道起伏位置四周只有四个网格点时，计算σ
xx
、σ
yy
采用不等距差分格式计算；而计算σ
xy
采用常规差分方法，则σ
xx
、σ
xy
、σ
yy
表达式如下：
[0049][0050]
⑥
起伏条件下σ
xx
、σ
yy
距离隧道起伏位置四周均有两个网格点，该分量上σ
xy
距离隧道起伏位置四周只有三个网格点时，计算σ
xx
、σ
xy
、σ
yy
采用不等距差分格式计算，则σ
xx
、σ
xy
、σ
yy
表达式如下：
[0051][0052]
⑦
起伏条件下σ
xx
、σ
yy
距离隧道起伏位置四周均有两个网格点，该分量上σ
xy
距离隧道起伏位置四周只有两个网格点时，计算σ
xx
、σ
xy
、σ
yy
采用不等距差分格式计算，则σ
xx
、σ
xy
、σ
yy
表达式如下：
[0053][0054]
⑧
起伏条件下σ
xx
、σ
yy
距离隧道起伏位置四周均有两个网格点，该分量上σ
xy
距离隧道起伏位置四周只有一个网格点时，只需要计算σ
xx
、σ
yy
采用不等距差分格式计算，则σ
xx
、σ
yy
表达式如下：
[0055][0056]ⅲ.当计算速度时，在水平或者垂直方向上共存在5种不同表达，具体为：
[0057]
①
在计算v
x
和vy时，水平和垂直方向上均可正常计算，采用常规差分格式计算，具体如下：
[0058][0059]
②
在计算vy时采用常规差分格式；而在计算v
x
时在垂直方向上存在不等距，改用不等距差分格式，具体如下：
[0060][0061]
③
在计算vy时采用常规差分格式；而计算v
x
时在水平和垂直方向上均存在不等距，改用不等距差分格式，具体如下：
[0062][0063]
④
在计算vy时采用常规差分格式；而v
x
在地表外时不需计算，具体如下：
[0064][0065]
⑤
计算vy时在水平方向上存在不等距，改用不等距差分格式；而v
x
在地表外时不需计算，具体如下：
[0066][0067]ⅳ.隧道的另一侧壁起伏情况按照上述步骤
ⅰ‑ⅲ
处理；
[0068]
ⅴ
.隧道掌子面处的起伏情况按照上述步骤
ⅰ‑ⅲ
处理。
[0069]
优选的，在所述步骤九中：相邻的所述三分量检波器之间的间隔为三米，所述掌子面设有用于设置所述三分量检波器的第一孔位，所述第一孔位的孔深为两米，所述三分量检波器布置在所述第一孔位的侧壁的中心位置。
[0070]
优选的，在所述步骤九中：所述三分量检波器的数量不大于31个。
[0071]
优选的，在所述步骤九中：相邻的所述震源之间的间隔为五米，所述掌子面设有用于设置所述震源的第二孔位，所述第二孔位的孔深为三米，所述震源布置在所述第二孔位的侧壁的中心位置。
[0072]
优选的，在所述步骤九中：所述震源的数量不大于17个。
[0073]
优选的，在所述步骤十二中：
[0074]
(1)本文采用各向同性介质的一阶速度应力二维弹性波波动方程进行：
[0075][0076][0077][0078][0079][0080]
上式中的σ
xx
、σ
yy
分别为沿x、y方向的正应力，λ、μ为拉梅参数，ρ为密度，v
x
,vy分别为质点沿着x、y方向的震动速度，t为时间；
[0081]
(2)交错网格下用高阶差分求解，具体差分方式如步骤八所示；
[0082]
(3)不断的进行迭代速度和应力，输出波场快照和地震记录文件。
[0083]
与现有技术相比，本发明的有益效果包括有：
[0084]
本方案基于不等距差分的隧道起伏地震进行超前地质预报探测，通过实际测量的隧道参数情况进行建模，对不同的网格点位置进行速度赋值，根据隧道的起伏边界不同，利用不等距网格差分法对模型进行数值计算，采用pml吸收边界，不断迭代速度和应力进行数值模拟，输出地震波场和地震记录，该方法能够有效的提高隧道超前地质预报探测数值模拟的准确性。
附图说明
[0085]
为了更清楚地说明本发明实施例技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
[0086]
图1为本发明方案应用的场景。
[0087]
图2是网格剖分的结构；
[0088]
图3为不等距差分应力表达示意图；
[0089]
图4为不等距差分速度表达示意图。
[0090]
其中：
[0091]
1-起伏隧道，2-断层构造，3-差分算子，4-起伏隧道的差分算子，5-不等距差分中心点，6-等距差分中心点，7-隧道起伏表面，8-x方向网格间距，9-z方向网格间距，10-高程z0。
具体实施方式
[0092]
为了能够更清楚地理解本发明的上述目的、特征和优点，下面结合附图和具体实施方式对本发明进行详细描述。需要说明的是，在不冲突的情况下，本技术的实施方式及实施方式中的特征可以相互组合。在下面的描述中阐述了很多具体细节以便于充分理解本发明，所描述的实施方式仅仅是本发明一部分实施方式，而不是全部的实施方式。基于本发明中的实施方式，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施方式，都属于本发明保护的范围。
[0093]
除非另有定义，本文所使用的所有的技术和科学术语与属于本发明的技术领域的技术人员通常理解的含义相同。本文中在本发明的说明书中所使用的术语只是为了描述具体的实施方式的目的，不是旨在于限制本发明。
[0094]
实施例：
[0095]
如图1-4所示，本实施例中提供一种地质预报探测断层二维数值模拟方法，包括以下步骤：
[0096]
步骤一：测量隧道掌子面反方向的隧道的起伏四周的起伏情况，构建起伏隧道模型；
[0097]
步骤二：建立均匀介质速度的三维模型，速度参数为v
p1
,v
s1
,ρ1，其中v
p1
,v
s1
,ρ1分别为三维模型的纵波速度、横波速度、密度；
[0098]
步骤三：将隧道高程参数信息插入速度模型，建立起伏隧道的速度模型，隧道模型速度参数为v
p2
,v
s2
,ρ2，其中v
p2
,v
s2
,ρ2分别为隧道空间的纵波速度、横波速度、密度；
[0099]
步骤四：在x方向处设置一个走向90度的断层，则断层部分的速度参数为v
p3
,v
s3
,ρ3，其中v
p3
,v
s3
,ρ3分别为断层部分的纵波速度、横波速度、密度；
[0100]
步骤五：沿隧道建立二维坐标系，x方向为指向掌子面，y方向为垂直于隧道侧壁，以模型边界上隧道面的中心为原点建立二维坐标系，建立隧道开挖超前地质预报观测系统；
[0101]
步骤六：给定网格间距和时间步长等参数，需要满足采样定律和稳定性条件；
[0102]
步骤七：通过真实高程计算真实地形点和隧道表面的网格点的距离；
[0103]
步骤八：在隧道处采用不等距有限差分格式，其余部分采用常规有限差分格式；
[0104]
步骤九：在起伏隧道的侧壁布置检波器，从掌子面往隧道方向进行依次均布n个三分量检波器和m个震源，其中三分量检波器号设为rn，震源号为sm，构建基底地震超前地质预报观测系统；
[0105]
步骤十：采用雷克子波、震源频率为z hz，采样间隔为l，采样时间为t；
[0106]
步骤十一：采用常规的pml吸收边界进行处理；
[0107]
步骤十二：迭代速度和应力，输出波场快照和地震记录。
[0108]
具体的，在步骤一中：
[0109]
(1)根据隧道情况，从掌子面往隧道方向测量90m范围内的整个隧道高程参数；
[0110]
(2)将高程参数转换成相对高程，构建整个隧道的基本形态；
[0111]
(3)通过隧道的基本形态，构建起伏隧道模型。
[0112]
具体的，在步骤二中：三维模型中，x方向、y方向长度分别为300m、100m；
[0113]
在步骤三中：v
p2
＝340m/s,v
s2
＝0m/s,ρ2＝1.29kg/m2。
[0114]
具体的，在步骤六中：
[0115]
(1)利用波动方程对网格进行网格化，大小为h；
[0116]
(2)设置数值模拟参数：频率为f，单位为hz，采样间隔为l，单位为ms；
[0117]
(3)为避免数值频散，要满足采样定理和稳定性条件，具体为网格大小和时间步长：
[0118]ⅰ.网格间距大小满足：
[0119]
其中，v
min
/(2fn)代表模拟的最小波长，v
min
为数值模拟参数中最小的速度参数，n为所取的差分阶数相关因子；
[0120]ⅱ.时间步长大小满足：
[0121]
其中，m是与所取差分阶数相关的因子，v
max
为为数值模拟参数中最小的速度参数。
[0122]
具体的，在步骤八中：
[0123]
(1)对模型进行网格化后，在起伏空气隧道侧壁外围采用常规等距离有限差分方法处理；
[0124]
(2)在起伏隧道附近进行不等距差分，主要分为应力和速度的不等距差分，具体如下：
[0125]ⅰ.通过u'＝r
β
u将起伏隧道坐标旋转为水平条件，其中u＝(v
x
,vy,σ
xx
,σ
xy
,σ
yy
)
t
,u'＝(v'
x
,v'y,σ'
xx
,σ'
xy
,σ'
yy
)
t
，v
x
,vy分别代表的速度分量，σ
xx
,σ
xy
,σ
yy
为应力分量，v'
x
,v'y为转换后的速度分量，σ'
xx
,σ'
xy
,σ'
yy
为转换后的应力分量，r
α
为坐标旋转矩阵，β为旋转角度，具体为：
[0126][0127]ⅱ.计算应力时，不等距差分共有8种情况，参见图3：
[0128]
①
起伏条件下σ
xx
、σ
xy
、σ
yy
距离隧道起伏位置四周均有四个网格点时，则不存在不等距差分，则σ
xx
、σ
xy
、σ
yy
表达式如下：
[0129][0130]
其中，δl、δh分别是采样间隔和网格大小，i,j分别代表x方向和y方向的网格点数，λ、μ为拉梅参数；
[0131]
②
起伏条件下σ
xx
、σ
yy
距离隧道起伏位置四周均有四个网格点，但σ
xy
距离隧道起伏位置四周只有三个网格点时，计算σ
xx
、σ
yy
采用常规差分格式计算；而计算σ
xy
时存在不等距差分计算，则σ
xx
、σ
xy
、σ
yy
表达式如下：
[0132][0133]
其中，y0(i)为隧道起伏的坐标参数，下述y1(i)为隧道坐标向下取整；
[0134]
③
起伏条件下σ
xx
、σ
yy
距离隧道起伏位置四周均有三个网格点，σ
xy
距离隧道起伏位置四周只有两个网格点时，则差分格式均为不等距差分格式，则σ
xx
、σ
xy
、σ
yy
表达式如下：
[0135][0136]
④
起伏条件下σ
xx
、σ
yy
距离隧道起伏位置四周均有三个网格点，该分量上σ
xy
距离隧道起伏位置四周只有一个网格点时，差分格式改用不等距差分格式，则只需要计算σ
xx
、σ
yy
，表达式如下：
[0137][0138]
⑤
起伏条件下σ
xx
、σ
yy
距离隧道起伏位置四周均有三个网格点，该分量上σ
xy
距离隧道起伏位置四周只有四个网格点时，计算σ
xx
、σ
yy
采用不等距差分格式计算；而计算σ
xy
采用常规差分方法，则σ
xx
、σ
xy
、σ
yy
表达式如下：
[0139][0140]
⑥
起伏条件下σ
xx
、σ
yy
距离隧道起伏位置四周均有两个网格点，该分量上σ
xy
距离隧道起伏位置四周只有三个网格点时，计算σ
xx
、σ
xy
、σ
yy
采用不等距差分格式计算，则σ
xx
、σ
xy
、σ
yy
表达式如下：
[0141][0142]
⑦
起伏条件下σ
xx
、σ
yy
距离隧道起伏位置四周均有两个网格点，该分量上σ
xy
距离隧道起伏位置四周只有两个网格点时，计算σ
xx
、σ
xy
、σ
yy
采用不等距差分格式计算，则σ
xx
、σ
xy
、σ
yy
表达式如下：
[0143][0144]
⑧
起伏条件下σ
xx
、σ
yy
距离隧道起伏位置四周均有两个网格点，该分量上σ
xy
距离隧道起伏位置四周只有一个网格点时，只需要计算σ
xx
、σ
yy
采用不等距差分格式计算，则σ
xx
、σ
yy
表达式如下：
[0145][0146]ⅲ.当计算速度时，在水平或者垂直方向上共存在5种不同表达，参见图4，具体为：
[0147]
①
在计算v
x
和vy时，水平和垂直方向上均可正常计算，采用常规差分格式计算，具体如下：
[0148][0149]
②
在计算vy时采用常规差分格式；而在计算v
x
时在垂直方向上存在不等距，改用不等距差分格式，具体如下：
[0150][0151]
③
在计算vy时采用常规差分格式；而计算v
x
时在水平和垂直方向上均存在不等距，改用不等距差分格式，具体如下：
[0152][0153]
④
在计算vy时采用常规差分格式；而v
x
在地表外时不需计算，具体如下：
[0154][0155]
⑤
计算vy时在水平方向上存在不等距，改用不等距差分格式；而v
x
在地表外时不需计算，具体如下：
[0156][0157]ⅳ.隧道的另一侧壁起伏情况按照上述步骤
ⅰ‑ⅲ
处理；
[0158]
ⅴ
.隧道掌子面处的起伏情况按照上述步骤
ⅰ‑ⅲ
处理。
[0159]
具体的，在步骤九中：相邻的三分量检波器之间的间隔为三米，掌子面设有用于设置三分量检波器的第一孔位，第一孔位的孔深为两米，三分量检波器布置在第一孔位的侧壁的中心位置。
[0160]
具体的，在步骤九中：三分量检波器的数量不大于31个。
[0161]
具体的，在步骤九中：相邻的震源之间的间隔为五米，掌子面设有用于设置震源的第二孔位，第二孔位的孔深为三米，震源布置在第二孔位的侧壁的中心位置。
[0162]
具体的，在步骤九中：震源的数量不大于17个。
[0163]
具体的，在步骤十二中：
[0164]
(1)本文采用各向同性介质的一阶速度应力二维弹性波波动方程进行：
[0165][0166][0167][0168][0169][0170]
上式中的σ
xx
、σ
yy
分别为沿x、y方向的正应力，λ、μ为拉梅参数，ρ为密度，v
x
,vy分别为质点沿着x、y方向的震动速度，t为时间；
[0171]
(2)交错网格下用高阶差分求解，具体差分方式如步骤八所示；
[0172]
(3)不断的进行迭代速度和应力，输出波场快照和地震记录文件。
[0173]
图1为本方案应用的场景，包括起伏隧道1和断层构造2；
[0174]
图2为网格剖分的结构，具体的为不等距差分网格，其中包括差分算子3(v
x
,vz)，起伏隧道的差分算子4(v
x
,vz)，不等距差分中心点5，等距差分中心点6，隧道起伏表面7，x方向网格间距8和z方向网格间距9，z方向网格的高程z0标号为10。
[0175]
隧道开挖过程中对围岩的破坏是不均一的，地震超前地质预报探测时检波器布置在参差不齐的隧道壁上，当震源在隧道侧壁进行激发时，复杂的隧道条件对地震波场的影响不容忽视。与在起伏地表条件不同，地面三维地震是半空间探测，隧道超前地质预报探测
为全空间条件下进行的，隧道空间侧壁起伏。本发明提出的模拟方法，通过实际测量的隧道参数情况进行建模，对不同的网格点位置进行速度赋值，根据隧道的起伏边界不同，利用不等距网格差分法对模型进行数值计算，采用pml吸收边界，不断迭代速度和应力进行数值模拟，输出地震波场和地震记录，该方法能够有效的提高隧道超前地质预报探测数值模拟的准确性。
[0176]
以上所述，仅是本发明的较佳实施例而已，并非对本发明作任何形式上的限制，故凡是未脱离本发明技术方案内容，依据本发明的技术实质对以上实施例所作的任何修改、等同变化与修饰,均仍属于本发明技术方案的范围内。
技术特征：1.一种地质预报探测断层二维数值模拟方法，其特征在于，包括以下步骤：步骤一：测量隧道掌子面反方向的隧道的起伏四周的起伏情况，构建起伏隧道模型；步骤二：建立均匀介质速度的三维模型，速度参数为v
p1
,v
s1
,ρ1，其中v
p1
,v
s1
,ρ1分别为三维模型的纵波速度、横波速度、密度；步骤三：将隧道高程参数信息插入速度模型，建立起伏隧道的速度模型，隧道模型速度参数为v
p2
,v
s2
,ρ2，其中v
p2
,v
s2
,ρ2分别为隧道空间的纵波速度、横波速度、密度；步骤四：在x方向处设置一个走向90度的断层，则断层部分的速度参数为v
p3
,v
s3
,ρ3，其中v
p3
,v
s3
,ρ3分别为断层部分的纵波速度、横波速度、密度；步骤五：沿隧道建立二维坐标系，x方向为指向掌子面，y方向为垂直于隧道侧壁，以模型边界上隧道面的中心为原点建立二维坐标系，建立隧道开挖超前地质预报观测系统；步骤六：给定网格间距和时间步长等参数，需要满足采样定律和稳定性条件；步骤七：通过真实高程计算真实地形点和隧道表面的网格点的距离；步骤八：在隧道处采用不等距有限差分格式，其余部分采用常规有限差分格式；步骤九：在起伏隧道的侧壁布置检波器，从掌子面往隧道方向进行依次均布n个三分量检波器和m个震源，其中三分量检波器号设为r
n
，震源号为s
m
，构建基底地震超前地质预报观测系统；步骤十：采用雷克子波、震源频率为z hz，采样间隔为l，采样时间为t；步骤十一：采用常规的pml吸收边界进行处理；步骤十二：迭代速度和应力，输出波场快照和地震记录。2.根据权利要求1所述的一种地质预报探测断层二维数值模拟方法，其特征在于，在所述步骤一中：(1)根据隧道情况，从掌子面往隧道方向测量90m范围内的整个隧道高程参数；(2)将高程参数转换成相对高程，构建整个隧道的基本形态；(3)通过隧道的基本形态，构建起伏隧道模型。3.根据权利要求1所述的一种地质预报探测断层二维数值模拟方法，其特征在于，在所述步骤二中：所述三维模型中，x方向、y方向长度分别为300m、100m；在步骤三中：v
p2
＝340m/s,v
s2
＝0m/s,ρ2＝1.29kg/m2。4.根据权利要求1所述的一种地质预报探测断层二维数值模拟方法，其特征在于，在所述步骤六中：(1)利用波动方程对网格进行网格化，大小为h；(2)设置数值模拟参数：频率为f，单位为hz，采样间隔为l，单位为ms；(3)为避免数值频散，要满足采样定理和稳定性条件，具体为网格大小和时间步长：ⅰ.网格间距大小满足：其中，v
min
/(2fn)代表模拟的最小波长，v
min
为数值模拟参数中最小的速度参数，n为所取的差分阶数相关因子；ⅱ.时间步长大小满足：
其中，m是与所取差分阶数相关的因子，v
max
为为数值模拟参数中最小的速度参数。5.根据权利要求1所述的一种地质预报探测断层二维数值模拟方法，其特征在于，在所述步骤八中：(1)对模型进行网格化后，在起伏空气隧道侧壁外围采用常规等距离有限差分方法处理；(2)在起伏隧道附近进行不等距差分，主要分为应力和速度的不等距差分，具体如下：ⅰ.通过u'＝r
β
u将起伏隧道坐标旋转为水平条件，其中u＝(v
x
,v
y
,σ
xx
,σ
xy
,σ
yy
)
t
,u'＝(v'
x
,v'
y
,σ'
xx
,σ'
xy
,σ'
yy
)
t
，v
x
,v
y
分别代表的速度分量，σ
xx
,σ
xy
,σ
yy
为应力分量，v'
x
,v'
y
为转换后的速度分量，σ'
xx
,σ'
xy
,σ'
yy
为转换后的应力分量，r
α
为坐标旋转矩阵，β为旋转角度，具体为：ⅱ.计算应力时，不等距差分共有8种情况：
①
起伏条件下σ
xx
、σ
xy
、σ
yy
距离隧道起伏位置四周均有四个网格点时，则不存在不等距差分，则σ
xx
、σ
xy
、σ
yy
表达式如下：其中，δl、δh分别是采样间隔和网格大小，i,j分别代表x方向和y方向的网格点数，λ、μ为拉梅参数；
②
起伏条件下σ
xx
、σ
yy
距离隧道起伏位置四周均有四个网格点，但σ
xy
距离隧道起伏位置四周只有三个网格点时，计算σ
xx
、σ
yy
采用常规差分格式计算；而计算σ
xy
时存在不等距差分计算，则σ
xx
、σ
xy
、σ
yy
表达式如下：其中，y0(i)为隧道起伏的坐标参数，下述y1(i)为隧道坐标向下取整；
③
起伏条件下σ
xx
、σ
yy
距离隧道起伏位置四周均有三个网格点，σ
xy
距离隧道起伏位置四周只有两个网格点时，则差分格式均为不等距差分格式，则σ
xx
、σ
xy
、σ
yy
表达式如下：
④
起伏条件下σ
xx
、σ
yy
距离隧道起伏位置四周均有三个网格点，该分量上σ
xy
距离隧道起伏位置四周只有一个网格点时，差分格式改用不等距差分格式，则只需要计算σ
xx
、σ
yy
，表达式如下：
⑤
起伏条件下σ
xx
、σ
yy
距离隧道起伏位置四周均有三个网格点，该分量上σ
xy
距离隧道起伏位置四周只有四个网格点时，计算σ
xx
、σ
yy
采用不等距差分格式计算；而计算σ
xy
采用常规差分方法，则σ
xx
、σ
xy
、σ
yy
表达式如下：
⑥
起伏条件下σ
xx
、σ
yy
距离隧道起伏位置四周均有两个网格点，该分量上σ
xy
距离隧道起伏位置四周只有三个网格点时，计算σ
xx
、σ
xy
、σ
yy
采用不等距差分格式计算，则σ
xx
、σ
xy
、σ
yy
表达式如下：
⑦
起伏条件下σ
xx
、σ
yy
距离隧道起伏位置四周均有两个网格点，该分量上σ
xy
距离隧道起伏位置四周只有两个网格点时，计算σ
xx
、σ
xy
、σ
yy
采用不等距差分格式计算，则σ
xx
、σ
xy
、σ
yy
表达式如下：
⑧
起伏条件下σ
xx
、σ
yy
距离隧道起伏位置四周均有两个网格点，该分量上σ
xy
距离隧道起伏位置四周只有一个网格点时，只需要计算σ
xx
、σ
yy
采用不等距差分格式计算，则σ
xx
、σ
yy
表达式如下：ⅲ.当计算速度时，在水平或者垂直方向上共存在5种不同表达，具体为：
①
在计算v
x
和v
y
时，水平和垂直方向上均可正常计算，采用常规差分格式计算，具体如下：
②
在计算v
y
时采用常规差分格式；而在计算v
x
时在垂直方向上存在不等距，改用不等距差分格式，具体如下：
③
在计算v
y
时采用常规差分格式；而计算v
x
时在水平和垂直方向上均存在不等距，改用不等距差分格式，具体如下：
④
在计算v
y
时采用常规差分格式；而v
x
在地表外时不需计算，具体如下：
⑤
计算v
y
时在水平方向上存在不等距，改用不等距差分格式；而v
x
在地表外时不需计算，具体如下：ⅳ.隧道的另一侧壁起伏情况按照上述步骤
ⅰ‑ⅲ
处理；
ⅴ
.隧道掌子面处的起伏情况按照上述步骤
ⅰ‑ⅲ
处理。6.根据权利要求1所述的一种地质预报探测断层二维数值模拟方法，其特征在于，在所述步骤九中：相邻的所述三分量检波器之间的间隔为三米，所述掌子面设有用于设置所述三分量检波器的第一孔位，所述第一孔位的孔深为两米，所述三分量检波器布置在所述第一孔位的侧壁的中心位置。7.根据权利要求6所述的一种地质预报探测断层二维数值模拟方法，其特征在于，在所述步骤九中：所述三分量检波器的数量不大于31个。8.根据权利要求1所述的一种地质预报探测断层二维数值模拟方法，其特征在于，在所述步骤九中：相邻的所述震源之间的间隔为五米，所述掌子面设有用于设置所述震源的第二孔位，所述第二孔位的孔深为三米，所述震源布置在所述第二孔位的侧壁的中心位置。9.根据权利要求1所述的一种地质预报探测断层二维数值模拟方法，其特征在于，在所述步骤九中：所述震源的数量不大于17个。10.根据权利要求1所述的一种地质预报探测断层二维数值模拟方法，其特征在于，在所述步骤十二中：(1)本文采用各向同性介质的一阶速度应力二维弹性波波动方程进行：(1)本文采用各向同性介质的一阶速度应力二维弹性波波动方程进行：(1)本文采用各向同性介质的一阶速度应力二维弹性波波动方程进行：
上式中的σ
xx
、σ
yy
分别为沿x、y方向的正应力，λ、μ为拉梅参数，ρ为密度，v
x
,v
y
分别为质点沿着x、y方向的震动速度，t为时间；(2)交错网格下用高阶差分求解，具体差分方式如步骤八所示；(3)不断的进行迭代速度和应力，输出波场快照和地震记录文件。
技术总结本发明属于隧道地球物理探测技术领域，提供一种地质预报探测断层二维数值模拟方法，利用不等距网格差分法对模型进行数值计算，采用PML吸收边界，不断迭代速度和应力进行数值模拟，输出地震波场和地震记录，该方法能够有效的提高隧道超前地质预报探测数值模拟的准确性。性。性。
技术研发人员：郑方坤 刘学 孙英
受保护的技术使用者：中交第四航务工程局有限公司
技术研发日：2022.06.17
技术公布日：2022/11/1
