本发明涉及地震波数值模拟，具体涉及到高铁震源的三维声波方程频率域高精度数值模拟方法。
背景技术：
1、传统地震勘探通常利用炸药、可控震源或气枪等震源激发地震波对地下空间进行成像。然而，此类震源成本昂贵，而且具有环境破坏性。现如今，随着交通需求量的骤增，高铁已经成为了人们快速出行的一个主要选择。行进中的高铁，由于其接近匀速的运行速度、确定的长度和载荷使其能够成为一个可以重复的人工震源，并且纵横交错的高铁铁路系统拓宽了这一人工震源的分布范围，从而使得它具有探测地球内部结构和物性参数的巨大潜力。与传统震源相比，高铁震源具有可重复性强、零成本、不破坏环境以及信号频谱信息丰富等优势。因此，利用移动高铁产生的振动信号对地下结构进行成像并估计相关的地下特征具有重要价值。
2、精确的高铁波场模拟技术是利用高铁波场开展近地表构造成像和参数反演的基础。现有的高铁波场模拟技术主要集中在两个方面：第一方面，基于全空间、半空间以及水平层状模型的解析和半解析类方法和第二方面，时间域数值模拟方法，比如有限差分方法和谱元法。而频率域的数值模拟方案主要是针对固定震源，用以合成传统震源激发的波场信号。综上，针对高铁波场的频率域数值模拟技术目前尚且处于空白状态。
3、利用解析或半解析类方法模拟高铁波场可以有效分析高铁波场的多域特征，但是此类方法只能处理均匀或者水平层状等几类简单的高对称性模型。针对复杂地下空间模型，难以找到相应的解析解。
4、针对上述问题，本发明提供了一种实现高精度频率域反演的高铁震源的三维声波方程频率域高精度数值模拟方法。
技术实现思路
1、本发明提供了一种实现高精度频率域反演的高铁震源的三维声波方程频率域高精度数值模拟方法。
2、本发明的目的是提供高铁震源的三维声波方程频率域高精度数值模拟方法，方案包括以下步骤：
3、步骤一、明确数值求解的定解问题，包括控制方程和边界条件，采用的控制方程为三维频率域声波方程，边界条件设定为边界场值为零，方程的震源项为高铁震源；
4、步骤二、采用33点有限差分格式对三维声波方程的空间导数项和质量加速项进行离散；
5、步骤三、不考虑震源项，在均匀无限介质中，求取33点差分格式的频散关系，并通过优化确定加权系数；
6、步骤四、为避免边界反射，在目标模型周边施加pml吸收边界条件；
7、步骤五、对高铁桥墩模型进行简化建模，构建高铁震源函数，并设计高铁震源的加载方式；
8、步骤六、获取高铁频率域单频波场。
9、进一步的，步骤一中的控制方程和边界条件的具体计算步骤如下：
10、步骤一、三维频率域声波方程及边界条件表示为：
11、
12、其中，x＝(x,y,z)表示空间位置，p为压力波场，ω为角频率，v是传播速度，为拉普拉斯算子，s为震源项；
13、步骤二、对于传统的固定震源，s通常可以表达为：
14、s(x,y,z,ω)＝s(ω)δ(x-xs)δ(y-ys)δ(z-zs),
15、其中，s(ω)为频率域的震源函数，(xs,ys,zs)为震源空间位置；
16、步骤三、当高铁沿着x方向匀速移动，则第m桥墩的震源函数可以表达为：
17、
18、其中，(xm,ys,zs)表示第m桥墩的空间位置，v表示高铁移动速度，表示虚数单位，q是调制因子表示为：
19、
20、从q的表达式中可以看出，高铁震源函数受车厢参数a和b以及移动速度v的调制，因子p表示为：
21、
22、其中，n为车厢个数。
23、进一步的，步骤二中的33点有限差分格式集体计算步骤如下：
24、步骤一、对于三维频率域声波方程的离散，我们提出一种新的33点有限差分格式，对于二阶空间导数项33点离散格式表示为：
25、
26、其中，表示平均导数27点算子，表示经典7点算子。
27、步骤二、平均导数27点算子表示为：
28、
29、其中，
30、
31、其中，pm,l,n≈p(mδx,lδy,nδz)，δx，δy以及δz分别表示在x,y以及z方向的网格间隔，aij(i,j＝1,2,3)为加权系数，表示为：
32、
33、步骤二、对于质量加速项，33点加权算子表示为：
34、
35、其中，
36、
37、其中，bi(i＝0,1,2,3,4)表示加权系数，因此，可获得33点有限差分离散方程为：
38、
39、步骤三、当w＝1以及b4＝0时，33点方法即退化为27点平均导数方法。
40、进一步的，步骤三中的频散关系具体计算步骤如下：
41、步骤一、评估33点差分格式的计算精度并确定其加权系数，进行经典的频散分析，考虑三维情况下的经典平面波解，具体表达式为：
42、
43、其中，p0是一个常数，k是波数，θ是传播角，φ是方位角；
44、步骤二、将步骤一的离散格式插入33点有限差分离散方程，可获得归一化相速度，具体表达式为：
45、
46、其中，vph是相速度，表示每个波长的网格点数(ppw)，不失一般性，δx＝max{δx,δy,δz}，
47、
48、
49、步骤三、为了压制频散误差，考虑如下l2范数的误差函数：
50、
51、其中，θ，φ以及的范围分别设置为以及[0,0.36]，用matlab中加约束的非线性优化程序fmincon进行优化。
52、进一步的，步骤四中的边界条件具体计算步骤如下：
53、步骤一、为了压制为限制模型大小而引入的边界所产生的人工反射，在目标模型的外围施加pml吸收边界条件，考虑pml吸收边界条件的频率域三维声波方程为：
54、
55、其中，
56、
57、其中，lx，ly以及lz分别是x，y以及z方向上pml的厚度，x，y以及z方向上的起点是内部模型的外边缘，cx＝cy＝cz＝180。
58、进一步的，步骤五的具体计算步骤如下：
59、步骤一、考虑模型的所有网格点，线性方程组表示为：
60、ap＝-s,
61、其中，a是大小为(nxnynz)×(nxnynz)的稀疏阻抗矩阵，p是大小为nxnynz×1的波场向量，s是大小为nxnynz×1的震源向量，nx，ny以及nz分别是x，y以及z方向的网格点数；
62、步骤二、在有限差分网格模板中，高铁震源可以表示为：
63、
64、其中，p是充零排序算子,确保s是nxnynz×1大小的向量；t表示矩阵转置，m为桥墩数。
65、本发明具有以下优势：本发明实现了频率域高铁波场的高精度数值模拟，填补了高铁波场频率域正演方法的空白，为频率域高铁波场的全波形反演提供了正演模块，这将为高铁数据成像提供新思路，有力推动高铁沿线及周边地区近地表成像的实际应用。
1.高铁震源的三维声波方程频率域高精度数值模拟方法，其特征在于：所述方法包括以下步骤：
2.如权利要求1所述的高铁震源的三维声波方程频率域高精度数值模拟方法，其特征在于：所述步骤一中的控制方程和边界条件的具体计算步骤如下：
3.如权利要求1所述的高铁震源的三维声波方程频率域高精度数值模拟方法，其特征在于：所述步骤二中的33点有限差分格式集体计算步骤如下：
4.如权利要求1所述的高铁震源的三维声波方程频率域高精度数值模拟方法，其特征在于：所述步骤三中的频散关系具体计算步骤如下：
5.如权利要求1所述的高铁震源的三维声波方程频率域高精度数值模拟方法，其特征在于：所述步骤四中的边界条件具体计算步骤如下：
6.如权利要求1所述的高铁震源的三维声波方程频率域高精度数值模拟方法，其特征在于：所述步骤五的具体计算步骤如下：
