多松弛时间格子Boltzmann方法的黏滞吸收边界

姜春涛 周 辉* 夏木明 唐瑾璇 王 颖

(①中国石油大学(北京)油气资源与探测国家重点实验室，北京 102249；②中国石油大学(北京)CNPC物探重点实验室，北京 102249；③中国科学院地质与地球物理研究所，北京 100029；④中国科学院地球科学研究院，北京 100029；⑤中国科学院油气资源研究重点实验室，北京 100029；⑥北京航天控制仪器研究所，北京 100094；⑦北京市光纤传感系统工程技术研究中心，北京 100094)

0 引言

地震波场模拟在地震勘探中占有举足轻重的地位，既是人们了解地震波在地下介质中传播规律的主要手段，也是地震反演、偏移成像和地质解释等[1-3]的基础。传统地震波场模拟一般是基于声波方程或弹性波方程[4]。而实际地下介质是非均匀、多相、非弹性甚至各向异性的，为了能更准确地描述地震波的传播规律，相继发展了黏弹性、多相以及各向异性波动方程模拟方法。时间域有限差分法和波数域伪谱法因为计算效率高、易于实现等优点而备受人们青睐。但这两种方法面对速度低且变化剧烈的复杂构造时，往往需要更小的时空采样步长，以解决数值频散和稳定性问题。除此之外，波动方程是一种宏观方程，因而必须要作介质连续性假设，对包含多相强间断面、多孔、多尺度介质，应用往往受限，模拟结果可信度降低。因此，很有必要发展一种更灵活、更适用于复杂介质的波场模拟技术。

格子Boltzmann方法(LBM)，是一种基于Boltzmann方程，以当代统计物理学为理论依据，求解Navier-Stokes方程的简单流体数值模拟方法[5]。它通过追踪介观尺度上离散粒子的相互作用，试图模拟宏观层面的复杂物理现象。鉴于该方法稳定性好、计算精度高、物理意义明确、不考虑非线性项、算法简单易实现且扩展迁移性强、天然的并行特性以及内部边界条件处理灵活等优势，且对复杂非线性问题和复杂边界适用性强，在计算流体力学、气动声学、热力学、电磁学以及地震学等领域[6-8]都有一定的应用。LBM自诞生以来，主要经历了细胞自动机(CA)、格子气自动机(LGA)、声格固体模型(PLS)和LBM模型等[9]几个发展过程。

在地震学领域，姚姚[10]及李幼铭等[11]使用CA模型进行地震模拟，刘舒考等[12]将LGA模型用于地震数据正演。针对CA模型在地震正演中的计算量大的问题，王真理等[13]提出了相应的并行算法。刘劲松等[14]利用CA模型进行地震模拟、偏移试验，取得了一定效果。贺艳晓[15]实现了PLS和LBM模型在简单非均匀介质中的波场模拟。针对黏弹介质波场模拟，Xia等[8]根据松弛时间与品质因子之间的定量关系，建立了LBM模型和地震波衰减参数的联系。Jiang等[16]用数值实验进一步探索了LBM在地震波场正演时的稳定性条件。凭借对微观力学模型研究成果，夏木明[17]结合格子Boltzmann模型和弹簧网络模型实现了流—固耦合复杂介质波场的正演模拟。

众所周知，在数值模拟过程中，由于计算区域有限，不可避免地会引入人工截断边界，如何施加吸收效果好的边界条件一直是一个热点问题[18]。目前关于LBM的无反射边界条件(NRBC)的研究主要集中在计算流体力学和气动声学领域[19-20]。在地震学领域，刘劲松等[21]提出在一定范围内的网格点上的准粒子数乘以一个衰减系数以实现边界反射的压制；Jiang等[22]引入了基于松弛时间的海绵层吸收边界，通过改变衰减系数及函数类型，在LBM波场数值模拟时取得了较好的边界反射压制效果。

以上关于LBM的吸收边界大多是基于单松弛时间格子Boltzmann模型(SRT-LBM)，依赖Bhatnagar-Gross-Krook(BGK)碰撞算子，且只考虑一个松弛时间。精度更高、物理意义更明确、能实现各物理模式解耦、稳定性更好，但碰撞过程更加复杂的多松弛时间格子Boltzmann模型(MRT-LBM)，用于地震模拟的吸收边界研究还未见到。相对于SRT-LBM，MRT-LBM的碰撞更加灵活，控制方程更复杂，吸收边界的加载就更困难。

为采用MRT-LBM进行地震波场模拟，本文提出一种基于多松弛参数的黏滞吸收边界。通过数值模拟实验优化衰减参数组合，利用均匀模型和非均质模型检验了该黏滞吸收边界的效果和适用性。

1 多松弛时间格子Boltzmann方法

1.1 离散速度模型

LBM是一种基于介观领域的数值求解方法，它描述了每个格子点上的粒子速度分布函数的演化过程。与传统宏观方程的离散方法不同，LBM不仅在物理空间和时间上进行离散，同时也需要在速度空间上进行离散。

LBM主要包含两个要素：离散速度模型和演化方程。根据维度不同，最常用的离散速度模型依次为D1Q3、D2Q9和D3Q19，其中D后的数字表示空间维度，Q后的数字表示离散速度对的个数，本文研究是基于D2Q9模型。

图1为D2Q9离散速度模型的基本结构，表示在中心点处的粒子每一时刻有9种不同的迁徙方向。图1的离散速度对为

(1)

1.2 演化方程

MRT-LBM的演化方程[5]为

f(x+ciΔt,t+Δt)-f(x,t)=

-M-1S[m(x,t)-meq(x,t)]

(2)

该式左边表示粒子在速度空间中的迁徙过程，右边表示粒子在速度矩空间中的碰撞过程。式中：Δt为时间采样间隔；M是转换矩阵

M=

(3)

其逆为

(4)

S是一个对角矩阵

S=diag[s0,s1,s2,s3,s4,s5,s6,s7,s8]

(5)

其元素是松弛参数；f(x,t)和feq(x,t)分别表示在t时刻、x处的粒子速度分布函数向量和相应的平衡态分布函数向量

(6)

m和meq分别为速度矩向量和平衡态速度矩向量，与f和feq的转换关系为

(7)

m的每个元素对应于宏观物理量[5]，可表示为

(8)

式中：ρ为介质密度；e为能量；ε为能量的平方；jx和jy是动量的x和y方向分量；qx和qy是能量通量的x和y方向分量；pxx和pxy是压力张量的对角线和非对角线成分。meq可以通过式(7)求得，即

(9)

式(6)中的平衡态分布函数[8]为

(10)

式中：u为介质速度；cs和wi分别为D2Q9模型对应的格子声速和权系数

(11)

LBM通过下式建立了微观量fi与宏观量ρ和u的关系[8]

(12)

使用MRT-LBM进行波场数值模拟主要分为以下五个步骤：

①初始化密度、速度、松弛参数、粒子数密度等物理量；

②计算平衡态分布函数；

③在速度矩空间进行碰撞过程以及外部吸收层区域的数值计算；

④在速度空间进行迁徙过程以及应用周期、镜面、外推等内部边界条件；

⑤根据粒子数密度计算宏观量密度、速度。

1.3 LBM与常规波动方程对比

有限差分法(FDM)与LBM在实现算法、数值离散方式和震源加载这三方面的对比如下。

(1)实现算法：LBM是一种介观尺度的数值算法，求解的是Navier-Stokes方程，其关键步骤是平衡态分布函数的求取、执行碰撞过程和迁徙过程、宏观量的求取。而FDM是一种宏观尺度的数值算法，可以求解常规的波动方程，其实现算法的关键步骤是将波动方程的时间或者空间偏导数用差分代替，迭代求取压力项或振动速度项。

(2)数值离散方式：LBM除了离散空间和时间外，还离散速度；FDM只离散空间和时间。

(3)震源加载：LBM震源一般加载在介观量粒子数密度上；FDM通常加载在宏观量压力或质点速度上。

2 基于多松弛参数的黏滞吸收边界

式(5)中的9个松弛参数对MRT-LBM的数值计算尤为重要，它们直接影响模拟波场的频散性、衰减性、各向异性、准确性、稳定性等[5]。松弛参数的大小表明了相应的速度矩向其平衡态矩线性松弛的速率

(13)

松弛参数s1、s7、s8与介质的运动剪切黏度ν和运动体积黏度μ有关[8]

(14)

该式表明，s1越小，μ越大，衰减越严重；同样地，s7越小，ν越大，衰减越严重。

鉴于此，为了实现MRT-LBM地震波场模拟时截断反射的压制，可以考虑在模拟区域外加一定厚度的吸收层，通过同时改变吸收层松弛参数s1以及s7和s8的取值，以改变黏滞性，从而实现边界反射波的“黏滞”吸收。

模拟区域设计如图2所示，分为内部区域和上、下、左、右四个梯形吸收区域。内部区域的松弛参数保持不变，吸收区域的松弛参数会随坐标x或者z而变化。以左边吸收区域(0≤x≤L-1)为例，松弛参数为

图2 模拟区域示意图

(15)

为简单起见，本文假设

(16)

很明显，参数L、n、a的选择对吸收效果影响很大。不失一般性，L越大，对边界反射压制更好，但是会增加计算负担；当L一定时，如何选取n和a的参数组合则显得尤为重要。

3 衰减参数组合优化

衰减系数a可以看作是吸收层厚度L和幂函数指数n的函数。对于MRT-LBM的黏滞吸收边界，由于没有经验公式可参考，因而首先需要进行参数优化。策略是通过选取不同的L、n、a参数组合在均匀模型中进行波场正演，记录直达波和截断边界反射的能量，进而计算信噪比[22]

SNR=10 lg(Ed/Er)

(17)

式中：Ed是直达波能量；Er是截断边界反射能量。参数组合的取值如下

(18)

均匀模型的参数设置如下：模拟网格数为401×401，空间步长为1m×1m；记录时长为0.4s，时间采样间隔为0.5ms；纵波速度为1500m/s，密度为1000kg/m3；35Hz主频的Ricker子波放在模型中心作为震源；设置四个检波点记录信号，坐标分别为(xs-160，zs-130)、(xs+100，zs+100)、(xs-120，zs+130)和(xs+30，zs-150)，其中(xs，zs)为源点的离散坐标；区域外部加上厚度为L的吸收层；内部区域的松弛参数设置为

(19)

为了保持模拟的稳定，各个松弛参数始终需满足

0

(20)

MRT-LBM在式(18)的L、n、a参数组合下进行若干次波场正演并通过四个检波点记录直达波和边界反射的能量。图3展示了不同参数组合时上述四个检波点信噪比的算术平均，很容易发现衰减系数的选取对边界的吸收效果影响很明显。实际参数优化的过程如下：

(1)先选择L。通过对比图3的6个子图，人为定义信噪比大于45dB且L在可接受的范围内，最终选择L=30。

(2)其次选择n和a。通过对比图3e的六条曲线的最大值，发现指数为3.0时信噪比最大，因此选择n=3.0；然后找到该曲线最大信噪比对应的衰减系数，可得a=0.001。

图3 不同参数组合下信噪比曲线

理论上，MRT-LBM的黏滞吸收边界采用30层厚度的计算量会大于常规PML厚度10～20层的计算量，然而此吸收边界形式简单，兼容性强，同样的层数不会带来过多的计算量。

进一步地，将该均匀模型扩大到801×801的网格数，时间记录1001个采样点，比较LBM与求解一阶压力—速度形式的二维波动方程的FDM在计算复杂度和计算效率方面的差异。计算复杂度可以分为时间复杂度和空间复杂度，一般地，这可分别由程序的实际运行时间(Time)与占用的最大物理内存百分比(Mem)表示。所有的模拟均在3.0GHz的英特尔处理器、内存为15.980GB的电脑上实施，采用C语言编程，结果如表1所示。从表1可以看出，LBM的计算效率比FDM低，其所需计算时间大约是FDM的两倍；LBM所需的计算内存大约是FDM的2～3倍。这主要是因为LBM在速度模型上进行了9个方向的离散，开辟的动态数组以及实际需要计算的物理量会相应地增加。

表1 FDM与LBM计算复杂度及计算效率对比

4 吸收效果验证

4.1 均匀介质模型

将上面优化后的衰减参数组合的黏滞吸收边界应用到基于MRT-LBM的均匀介质的波场模拟，模型参数设置同上。

模拟的波场切片及地震记录如图4a和图4b所示，未见明显的截断边界反射。第4个检波点记录的振动曲线如图4c所示，并将扩大的模型的相应记录作为参考解(边界反射出现在考察的时间之外)进行对比。由图4c可见，施加吸收边界的振幅曲线与参考解只有微小的差异，说明边界反射吸收得比较好。图5为不同吸收层厚度的地震波场总能量(某一时刻波场所有点的振幅的平方和)随时间的动态变化，模拟总时长为1s。由图5可见，选取优化后的衰减参数组合，波场能量在通过吸收层时迅速衰减(大约5个数量级)，且当剩余反射重新到达边界时会继续得到大幅度的压制。

图4 均匀模型的波场切片及地震记录和单道波形

4.2 非均质介质模型

进一步通过测试非均质模型的波场模拟检验本文的黏滞吸收边界优化参数的有效性。若无特殊说明，黏滞吸收边界采用的衰减参数组合均为L=30、n=3、a=0.001。

首先考察两个简单的非均质模型：层状模型(图6a)和孔洞模型(图6b)的波场模拟结果。层状模型的速度为1500、3000m/s，以20 Hz的雷克子波作为震源放在地表中心(五星所示)，其他参数设置同均质模型。图7为t=325ms的波场切片以及x=50m和z=150m的地震记录。可以发现，相对于直达波和反射波，截断边界反射波并不明显，证实了本文的黏滞吸收边界对MRT-LBM波场模拟时产生的截断边界反射波有较好的衰减效果。

图6 简单非均质模型

图7 层状模型的波场切片及地震记录

孔洞模型的网格数为801×801，中间异常体的尺寸为267m×267m，速度为1500m/s，围岩的速度为3000m/s，模拟时长为0.5s，只在模型的上表面加载了黏滞吸收边界，其他设置同均匀模型。图8为t=362.5ms的波场切片，由图可知，左、右、下三个边界的反射(黑色箭头处)很明显，而加载吸收边界的上部边界的截断边界反射显著衰减。为了更清楚地展示上边界的吸收效果，抽取了x=400m，t=225.0ms和362.5ms的波形曲线。当t=225.0ms时(图9)直达波并未到达模型边界，因而未产生截断边界反射，波形曲线关于z=400m对称；当t=362.5ms时(图10)，在z=600～800m之间的下边界反射很严重，已经干扰了有效层间反射波的识别，而z=0～200m的上边界反射被吸收得比较彻底。

图8 孔洞模型t=362.5ms时的波场切片

图9 孔洞模型x=400m、t=225.0ms的波形曲线

图10 孔洞模型x=400m、t=362.5ms的波形曲线

最后，用沿深度方向拉伸的BP模型(图11)验证本文的黏滞吸收边界对复杂模型的适用性。模型网格数为483×398，最大速度为4500m/s，最小速度为1500m/s，25Hz主频的雷克子波作为震源加在地表中心向下60个网格处(红五星位置)，其他参数设置同均匀模型。

图11 沿深度拉伸的BP速度模型

图12为MRT-LBM模拟的BP模型t=125ms的波场切片及x=150m和z=50m的地震记录。其中图12a未作吸收边界处理，可以看出，直达波产生的反射(青色箭头处)极其严重，明显覆盖了有效信号；对比而言，四周加载黏滞吸收边界后的波场(图12b)中的截断边界反射得以明显压制，有效信号更加清楚；图12c为图12a和图12b的差，可以清晰地看到，除了直达波导致的边界反射外，还有很多反射波引起的边界噪声。图12证实了本文的吸收边界条件适用于复杂模型的MRT-LBM地震波场模拟。

图12 BP模型计算的波场切片及地震记录

5 结论与讨论

考虑到MRT-LBM地震波动正演模拟的特性，本文提出了一种基于多松弛参数的黏滞吸收边界，该边界算法简单，易实现且可扩展性强。考虑到衰减参数组合对边界反射吸收效果影响显著，因而通过数值实验对其进行了优化。

用均均匀模型、层状模型、孔洞模型和复杂的BP模型进行了测试，结果表明本文的黏滞吸收边界条件较好地压制了边界反射，且适用复杂模型地震波场模拟。

需要注意的是，本文衰减参数组合只尝试了几种幂函数，衰减系数以10的整数倍为间隔，因而得到的未必是最优的参数组合。还可以通过尝试新的衰减函数，比如Sigmoid函数、正弦函数、反正切函数等[23]，以及二次加密衰减系数的采样间隔等方式进一步优化。另外，为了简单起见，本文实现边界衰减的三个松弛参数采用的是同一套衰减体系，为了使边界吸收效果更好，各个松弛参数可以尝试不同的衰减过程，但是要注意在松弛参数变化的过程中要满足MRT-LBM的稳定性条件。最后，可以采用并行方式处理MRT-LBM计算量大的问题。