多松弛時間格子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計算量大的問題。