基于快速多极间接边界元法的高频地震波散射模拟方法
【技术领域】
[0001] 本发明涉及地震波散射模拟领域，具体是一种基于快速多极间接边界元法的高频 地震波散射模拟方法。
【背景技术】
[0002] 近年来，随着计算机硬件水平飞速发展，大规模科学计算对于现代科技研发和大 型工程设计愈发重要。其中不规则地形对地震波的散射模拟是很多学科领域的研究热点， 可广泛应用于复杂地形地震动科学确定、地球物理勘探、岩土层参数反演等。理论分析方法 整体上可分为解析法和数值法，其中数值法包括域离散型的有限元法、有限差分法和边界 型的边界单元法等。同域离散型方法相比，边界元法（BEM)具有高精度、降低问题求解维 度、便于处理高梯度应力变化的优点，同时特别便于处理无限域波动问题（无需引入无反 射边界条件），另外还克服了普通有限元法高频弥散的难题。但该方法在实际应用中的主要 缺陷在于待求解方程系数矩阵为非对称稠密满阵，这样对于目前越发流行的大规模科学运 算，其计算效率将大大降低，在三维问题上该弱点将更为明显。
[0003] 针对传统BEM的稠密矩阵特性，近些年来，很多学者开展了快速BEM的研究工作， 包括快速多极子展开（FMM)、小波方法、ACA技术、PFFT方法等。其中快速多极子方法发展 较早，目前理论构架已臻于成熟。其基本原理是将多个源点和场点之间的相互作用运算由 逐粒计算，通过多极展开及局部展开（聚合、转移），转换为粒组与粒组之间的计算，能够极 大地降低大规模多点相互作用运算所需的计算量和存储量。目前该方法已被广泛应用到热 学、电磁学、固体力学、流体动力学和声学等领域的大型复杂问题分析。需指出的是，对于固 体中弹性波散射问题，由于涉及到压缩波和剪切波间的波型转换（在不连续面附近可能还 衍生面波），波动特征将更为复杂，因而FMM-BEM方法在该领域的发展较为缓慢。已有的快 速多极直接边界元法（FMM-DBEM)由于其数学公式推导复杂，实现难度大以及收敛效果不 佳等原因一直处于待完善状态。
【发明内容】
[0004] 本发明针对不规则地形对高频地震波散射问题，提出了一种基于快速多极间接边 界元法的高频地震波散射模拟方法。本方法采用快速多极间接边界元方法（FMM-IBEM)4B 够快速求解其他方法（解析法，有限差分法，有限元法以及常规边界元法）无法完成的高频 地震波散射问题，能在一台普通个人电脑上完成上百万自由度问题的三维散射问题，为高 频或大尺度地震波散射模拟提供了一种有效方法，提高计算速度的同时大大降低电脑资源 占用，为实际工程的快速模拟带来显著的经济效益。
[0005] 为实现上述目的，本发明提供如下技术方案：
[0006] -种基于快速多极间接边界元法的高频地震波散射模拟方法，具体步骤为：
[0007] (1)建立边界元模型
[0008] 考虑平面波入射，基于单层位势理论，在边界表面上施加三个正交方向的虚拟荷 载以构造散射波场；在建立边界元模型中，首先将边界表面离散为三角形或四边形，得到离 散单元；采用斜面圆盘均布荷载近似覆盖在离散单元上；
[0009] ⑵FMM-IBEM实施过程
[0010] (2. 1)边界离散并生成树结构
[0011] 将边界离散后，采用第0层正方体包含整个计算域，将第0层正方体等分为四个第 一层正方体，再将每个第一层正方体等分为四个第二层正方体，这样逐级划分下去；当最小 正方体的边长等于给定的数值时，停止划分；然后删去不含任何单元点的正方体，即可生成 树结构；定义第L层的正方体为父细胞，L彡0 ;而L+1层的正方体为其子结点，1彡1 ;包 含边界单元点的最小正方体称之为叶子结点；第0层正方体为根结点；其他正方体为树结 占.
[0012] (2. 2)近场计算
[0013] 对于叶子结点，其邻居以及自身立方体内场点和源点间相互作用通过间接边界元 法IBEM计算；这里的场点指单元内部的积分点，源点指虚拟荷载和流量源作用圆盘面上的 离散点；
[0014] (2. 3)下行遍历计算多极展开系数
[0015] 从叶子结点开始，将离散单元上虚拟均布荷载的作用传递到叶子结点，形成最初 的多极展开系数，然后通过多极展开中心点间的转移公式Μ2Μ，逐层向上递归操作，直到求 出第二层树结点的多极展开系数。
[0016] (2. 4)上行遍历逐层计算局部展开系数
[0017] 对第二层的每一个树节点，利用多极展开向局部展开转移公式M2L，将相互作用 列表中的所有节点的多极展开系数传递到该树节点的局部展开系数；对于父结点中所有源 点的影响，通过局部展开中心点间的转移公式L2L得到，最后累加到该树节点的局部展开 系数；如此直到求出叶子结点的局部展开系数；
[0018] (2. 5)边界积分和迭代求解
[0019] 当树结构的多极展开系数和局部展开系数计算完毕后，利用树结构对不同单元的 边界积分求和；对某个单元的场点X，根据距离远近，源点类型分为三类：①对其中的近场 源点，即邻居及本身，使用常规方法直接计算；②对该场点所在叶子"相互作用列表"中的 源点，采用源点所在叶子结点的多极展开系数计算；③对于更远的源点，即上层相互作用列 表，它们包含的边界单元积分已通过多极展开系数的转移、叠加和传递，集成到了 X所在的 叶子结点的局部展开系数；这样逐个单元计算下去，使用树结构来代替系数矩阵与迭代向 量相乘所需的存储和计算，即可得到各个单元上虚拟均布荷载的密度；然后采用GMRES算 法迭代求解直到满足精度要求；
[0020] (2. 6)总波场计算
[0021] 得到各个单元上虚拟均布荷载的密度之后，采用多极展开计算各个单元对所关心 区域的位移格林函数，即得到散射波场；最后将散射波场和自由场叠加，即可得到总体波 场。
[0022] 作为本发明进一步的方案：步骤（2. 5)中GMRES算法迭代求解采用预处理技术，即 拟采用近场格林函数解构造一窄带稀疏矩阵，对方程组进行预处理，以加速GMRES方程求 解收敛。
[0023] 作为本发明进一步的方案：该方法替换了常规IBffl的位移格林函数和应力格林 函数的展开计算方式；常规IBEM的位移格林函数和应力格林函数：
[0026] 上式中i, j = 1，2, 3,对应于X，y, z方向；η;为场点X处边界单元单位法向量与i 轴正方向夹角余弦；
[0027] 其中位势函数：
[0029] 其中q为给定常数，后文取为k，h (横波和纵波波数），复数i2= -1，r = I x-y I，X 和y分别为场点和源点；
[0030] 本方法使用树结构作为主要存储和运算对象，对核函数进行展开和传递，借助 GMRES迭代算法，在每一次迭代中以树结构取代系数矩阵与迭代量相乘；通过迭代精度控 制得出结果；三维问题采用八叉树结构，对任意几何形状的弹性体边界面，都会生成与之对 应的、具有特定层次和分布的树结构；
[0031] 根据平面波展开理论，势函数（3)可多极展开为：
?"为[-1，1]上高斯积分的横坐标及权系数，（ = ]^/|(|，r=|M,满足且 |5| < · O)为第2类第1阶球汉克尔函数，勒让德函数；
[0037] 从而，得到了新的位移格林函数和应力格林函数展开式：
[0042] 其中，C；1表示或V为局部展开系数；上述公式⑶和公式（9)采用树 结构存储各局部展开系数及多极展开系数。
[0043] 与现有技术相比，本发明的有益效果是：
[0044] 本发明针对大型实际复杂非水平层状沉积盆地对地震波的散射模拟存在的不足， 例如塔里木盆地，四川盆地等大型盆地存在的低频不稳定性。本发明主要解决了场地较为 坚硬（剪切波速约为2000m/s~4000m/s)的大小大约在一百公里范围内的实际盆地的地 震动模拟、场地较软的中小型盆地模拟和高频地震波散射模拟问题，实现了在普通个人电 脑上快速模拟高频或大型不规则地形地震动问题，使该问题的模拟实现普遍化和经济化。
[0045] 本发明结合快速多极算法（FMM)和间接边界元法（IBEM)，给出了基于对角形式的 平面波展开的快速多极间接边界元法，该方法数学公式推导简洁易懂，实现形式简便直观