姜春濤 周 輝* 夏木明 唐瑾璇 王 穎

(①中國石油大學(北京)油氣資源與探測國家重點實驗室，北京 102249；②中國石油大學(北京)CNPC物探重點實驗室，北京 102249；③中國科學院地質與地球物理研究所，北京 100029；④中國科學院地球科學研究院，北京 100029；⑤中國科學院油氣資源研究重點實驗室，北京 100029；⑥北京航天控制儀器研究所，北京 100094；⑦北京市光纖傳感系統(tǒng)工程技術研究中心，北京 100094)

0 引言

地震波場模擬在地震勘探中占有舉足輕重的地位，既是人們了解地震波在地下介質中傳播規(guī)律的主要手段，也是地震反演、偏移成像和地質解釋等[1-3]的基礎。傳統(tǒng)地震波場模擬一般是基于聲波方程或彈性波方程[4]。而實際地下介質是非均勻、多相、非彈性甚至各向異性的，為了能更準確地描述地震波的傳播規(guī)律，相繼發(fā)展了黏彈性、多相以及各向異性波動方程模擬方法。時間域有限差分法和波數(shù)域偽譜法因為計算效率高、易于實現(xiàn)等優(yōu)點而備受人們青睞。但這兩種方法面對速度低且變化劇烈的復雜構造時，往往需要更小的時空采樣步長，以解決數(shù)值頻散和穩(wěn)定性問題。除此之外，波動方程是一種宏觀方程，因而必須要作介質連續(xù)性假設，對包含多相強間斷面、多孔、多尺度介質，應用往往受限，模擬結果可信度降低。因此，很有必要發(fā)展一種更靈活、更適用于復雜介質的波場模擬技術。

格子Boltzmann方法(LBM)，是一種基于Boltzmann方程，以當代統(tǒng)計物理學為理論依據(jù)，求解Navier-Stokes方程的簡單流體數(shù)值模擬方法[5]。它通過追蹤介觀尺度上離散粒子的相互作用，試圖模擬宏觀層面的復雜物理現(xiàn)象。鑒于該方法穩(wěn)定性好、計算精度高、物理意義明確、不考慮非線性項、算法簡單易實現(xiàn)且擴展遷移性強、天然的并行特性以及內部邊界條件處理靈活等優(yōu)勢，且對復雜非線性問題和復雜邊界適用性強，在計算流體力學、氣動聲學、熱力學、電磁學以及地震學等領域[6-8]都有一定的應用。LBM自誕生以來，主要經(jīng)歷了細胞自動機(CA)、格子氣自動機(LGA)、聲格固體模型(PLS)和LBM模型等[9]幾個發(fā)展過程。

在地震學領域，姚姚[10]及李幼銘等[11]使用CA模型進行地震模擬，劉舒考等[12]將LGA模型用于地震數(shù)據(jù)正演。針對CA模型在地震正演中的計算量大的問題，王真理等[13]提出了相應的并行算法。劉勁松等[14]利用CA模型進行地震模擬、偏移試驗，取得了一定效果。賀艷曉[15]實現(xiàn)了PLS和LBM模型在簡單非均勻介質中的波場模擬。針對黏彈介質波場模擬，Xia等[8]根據(jù)松弛時間與品質因子之間的定量關系，建立了LBM模型和地震波衰減參數(shù)的聯(lián)系。Jiang等[16]用數(shù)值實驗進一步探索了LBM在地震波場正演時的穩(wěn)定性條件。憑借對微觀力學模型研究成果，夏木明[17]結合格子Boltzmann模型和彈簧網(wǎng)絡模型實現(xiàn)了流—固耦合復雜介質波場的正演模擬。

眾所周知，在數(shù)值模擬過程中，由于計算區(qū)域有限，不可避免地會引入人工截斷邊界，如何施加吸收效果好的邊界條件一直是一個熱點問題[18]。目前關于LBM的無反射邊界條件(NRBC)的研究主要集中在計算流體力學和氣動聲學領域[19-20]。在地震學領域，劉勁松等[21]提出在一定范圍內的網(wǎng)格點上的準粒子數(shù)乘以一個衰減系數(shù)以實現(xiàn)邊界反射的壓制；Jiang等[22]引入了基于松弛時間的海綿層吸收邊界，通過改變衰減系數(shù)及函數(shù)類型，在LBM波場數(shù)值模擬時取得了較好的邊界反射壓制效果。

以上關于LBM的吸收邊界大多是基于單松弛時間格子Boltzmann模型(SRT-LBM)，依賴Bhatnagar-Gross-Krook(BGK)碰撞算子，且只考慮一個松弛時間。精度更高、物理意義更明確、能實現(xiàn)各物理模式解耦、穩(wěn)定性更好，但碰撞過程更加復雜的多松弛時間格子Boltzmann模型(MRT-LBM)，用于地震模擬的吸收邊界研究還未見到。相對于SRT-LBM，MRT-LBM的碰撞更加靈活，控制方程更復雜，吸收邊界的加載就更困難。

為采用MRT-LBM進行地震波場模擬，本文提出一種基于多松弛參數(shù)的黏滯吸收邊界。通過數(shù)值模擬實驗優(yōu)化衰減參數(shù)組合，利用均勻模型和非均質模型檢驗了該黏滯吸收邊界的效果和適用性。

1 多松弛時間格子Boltzmann方法

1.1 離散速度模型

LBM是一種基于介觀領域的數(shù)值求解方法，它描述了每個格子點上的粒子速度分布函數(shù)的演化過程。與傳統(tǒng)宏觀方程的離散方法不同，LBM不僅在物理空間和時間上進行離散，同時也需要在速度空間上進行離散。

LBM主要包含兩個要素：離散速度模型和演化方程。根據(jù)維度不同，最常用的離散速度模型依次為D1Q3、D2Q9和D3Q19，其中D后的數(shù)字表示空間維度，Q后的數(shù)字表示離散速度對的個數(shù)，本文研究是基于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)

其元素是松弛參數(shù)；f(x,t)和feq(x,t)分別表示在t時刻、x處的粒子速度分布函數(shù)向量和相應的平衡態(tài)分布函數(shù)向量

(6)

m和meq分別為速度矩向量和平衡態(tài)速度矩向量，與f和feq的轉換關系為

(7)

m的每個元素對應于宏觀物理量[5]，可表示為

(8)

式中：ρ為介質密度；e為能量；ε為能量的平方；jx和jy是動量的x和y方向分量；qx和qy是能量通量的x和y方向分量；pxx和pxy是壓力張量的對角線和非對角線成分。meq可以通過式(7)求得，即

(9)

式(6)中的平衡態(tài)分布函數(shù)[8]為

(10)

式中：u為介質速度；cs和wi分別為D2Q9模型對應的格子聲速和權系數(shù)

(11)

LBM通過下式建立了微觀量fi與宏觀量ρ和u的關系[8]

(12)

使用MRT-LBM進行波場數(shù)值模擬主要分為以下五個步驟：

①初始化密度、速度、松弛參數(shù)、粒子數(shù)密度等物理量；

②計算平衡態(tài)分布函數(shù)；

③在速度矩空間進行碰撞過程以及外部吸收層區(qū)域的數(shù)值計算；

④在速度空間進行遷徙過程以及應用周期、鏡面、外推等內部邊界條件；

⑤根據(jù)粒子數(shù)密度計算宏觀量密度、速度。

1.3 LBM與常規(guī)波動方程對比

有限差分法(FDM)與LBM在實現(xiàn)算法、數(shù)值離散方式和震源加載這三方面的對比如下。

(1)實現(xiàn)算法：LBM是一種介觀尺度的數(shù)值算法，求解的是Navier-Stokes方程，其關鍵步驟是平衡態(tài)分布函數(shù)的求取、執(zhí)行碰撞過程和遷徙過程、宏觀量的求取。而FDM是一種宏觀尺度的數(shù)值算法，可以求解常規(guī)的波動方程，其實現(xiàn)算法的關鍵步驟是將波動方程的時間或者空間偏導數(shù)用差分代替，迭代求取壓力項或振動速度項。

(2)數(shù)值離散方式：LBM除了離散空間和時間外，還離散速度；FDM只離散空間和時間。

(3)震源加載：LBM震源一般加載在介觀量粒子數(shù)密度上；FDM通常加載在宏觀量壓力或質點速度上。

2 基于多松弛參數(shù)的黏滯吸收邊界

式(5)中的9個松弛參數(shù)對MRT-LBM的數(shù)值計算尤為重要，它們直接影響模擬波場的頻散性、衰減性、各向異性、準確性、穩(wěn)定性等[5]。松弛參數(shù)的大小表明了相應的速度矩向其平衡態(tài)矩線性松弛的速率

(13)

松弛參數(shù)s1、s7、s8與介質的運動剪切黏度ν和運動體積黏度μ有關[8]

(14)

該式表明，s1越小，μ越大，衰減越嚴重；同樣地，s7越小，ν越大，衰減越嚴重。

鑒于此，為了實現(xiàn)MRT-LBM地震波場模擬時截斷反射的壓制，可以考慮在模擬區(qū)域外加一定厚度的吸收層，通過同時改變吸收層松弛參數(shù)s1以及s7和s8的取值，以改變黏滯性，從而實現(xiàn)邊界反射波的“黏滯”吸收。

模擬區(qū)域設計如圖2所示，分為內部區(qū)域和上、下、左、右四個梯形吸收區(qū)域。內部區(qū)域的松弛參數(shù)保持不變，吸收區(qū)域的松弛參數(shù)會隨坐標x或者z而變化。以左邊吸收區(qū)域(0≤x≤L-1)為例，松弛參數(shù)為

圖2 模擬區(qū)域示意圖

(15)

為簡單起見，本文假設

(16)

很明顯，參數(shù)L、n、a的選擇對吸收效果影響很大。不失一般性，L越大，對邊界反射壓制更好，但是會增加計算負擔；當L一定時，如何選取n和a的參數(shù)組合則顯得尤為重要。

3 衰減參數(shù)組合優(yōu)化

衰減系數(shù)a可以看作是吸收層厚度L和冪函數(shù)指數(shù)n的函數(shù)。對于MRT-LBM的黏滯吸收邊界，由于沒有經(jīng)驗公式可參考，因而首先需要進行參數(shù)優(yōu)化。策略是通過選取不同的L、n、a參數(shù)組合在均勻模型中進行波場正演，記錄直達波和截斷邊界反射的能量，進而計算信噪比[22]

SNR=10 lg(Ed/Er)

(17)

式中：Ed是直達波能量；Er是截斷邊界反射能量。參數(shù)組合的取值如下

(18)

均勻模型的參數(shù)設置如下：模擬網(wǎng)格數(shù)為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)為源點的離散坐標；區(qū)域外部加上厚度為L的吸收層；內部區(qū)域的松弛參數(shù)設置為

(19)

為了保持模擬的穩(wěn)定，各個松弛參數(shù)始終需滿足

0

(20)

MRT-LBM在式(18)的L、n、a參數(shù)組合下進行若干次波場正演并通過四個檢波點記錄直達波和邊界反射的能量。圖3展示了不同參數(shù)組合時上述四個檢波點信噪比的算術平均，很容易發(fā)現(xiàn)衰減系數(shù)的選取對邊界的吸收效果影響很明顯。實際參數(shù)優(yōu)化的過程如下：

(1)先選擇L。通過對比圖3的6個子圖，人為定義信噪比大于45dB且L在可接受的范圍內，最終選擇L=30。

(2)其次選擇n和a。通過對比圖3e的六條曲線的最大值，發(fā)現(xiàn)指數(shù)為3.0時信噪比最大，因此選擇n=3.0；然后找到該曲線最大信噪比對應的衰減系數(shù)，可得a=0.001。

圖3 不同參數(shù)組合下信噪比曲線

理論上，MRT-LBM的黏滯吸收邊界采用30層厚度的計算量會大于常規(guī)PML厚度10～20層的計算量，然而此吸收邊界形式簡單，兼容性強，同樣的層數(shù)不會帶來過多的計算量。

進一步地，將該均勻模型擴大到801×801的網(wǎng)格數(shù)，時間記錄1001個采樣點，比較LBM與求解一階壓力—速度形式的二維波動方程的FDM在計算復雜度和計算效率方面的差異。計算復雜度可以分為時間復雜度和空間復雜度，一般地，這可分別由程序的實際運行時間(Time)與占用的最大物理內存百分比(Mem)表示。所有的模擬均在3.0GHz的英特爾處理器、內存為15.980GB的電腦上實施，采用C語言編程，結果如表1所示。從表1可以看出，LBM的計算效率比FDM低，其所需計算時間大約是FDM的兩倍；LBM所需的計算內存大約是FDM的2～3倍。這主要是因為LBM在速度模型上進行了9個方向的離散，開辟的動態(tài)數(shù)組以及實際需要計算的物理量會相應地增加。

表1 FDM與LBM計算復雜度及計算效率對比

4 吸收效果驗證

4.1 均勻介質模型

將上面優(yōu)化后的衰減參數(shù)組合的黏滯吸收邊界應用到基于MRT-LBM的均勻介質的波場模擬，模型參數(shù)設置同上。

模擬的波場切片及地震記錄如圖4a和圖4b所示，未見明顯的截斷邊界反射。第4個檢波點記錄的振動曲線如圖4c所示，并將擴大的模型的相應記錄作為參考解(邊界反射出現(xiàn)在考察的時間之外)進行對比。由圖4c可見，施加吸收邊界的振幅曲線與參考解只有微小的差異，說明邊界反射吸收得比較好。圖5為不同吸收層厚度的地震波場總能量(某一時刻波場所有點的振幅的平方和)隨時間的動態(tài)變化，模擬總時長為1s。由圖5可見，選取優(yōu)化后的衰減參數(shù)組合，波場能量在通過吸收層時迅速衰減(大約5個數(shù)量級)，且當剩余反射重新到達邊界時會繼續(xù)得到大幅度的壓制。

圖4 均勻模型的波場切片及地震記錄和單道波形

4.2 非均質介質模型

進一步通過測試非均質模型的波場模擬檢驗本文的黏滯吸收邊界優(yōu)化參數(shù)的有效性。若無特殊說明，黏滯吸收邊界采用的衰減參數(shù)組合均為L=30、n=3、a=0.001。

首先考察兩個簡單的非均質模型：層狀模型(圖6a)和孔洞模型(圖6b)的波場模擬結果。層狀模型的速度為1500、3000m/s，以20 Hz的雷克子波作為震源放在地表中心(五星所示)，其他參數(shù)設置同均質模型。圖7為t=325ms的波場切片以及x=50m和z=150m的地震記錄?？梢园l(fā)現(xiàn)，相對于直達波和反射波，截斷邊界反射波并不明顯，證實了本文的黏滯吸收邊界對MRT-LBM波場模擬時產生的截斷邊界反射波有較好的衰減效果。

圖6 簡單非均質模型

圖7 層狀模型的波場切片及地震記錄

孔洞模型的網(wǎng)格數(shù)為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之間的下邊界反射很嚴重，已經(jīng)干擾了有效層間反射波的識別，而z=0～200m的上邊界反射被吸收得比較徹底。

圖8 孔洞模型t=362.5ms時的波場切片

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

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

最后，用沿深度方向拉伸的BP模型(圖11)驗證本文的黏滯吸收邊界對復雜模型的適用性。模型網(wǎng)格數(shù)為483×398，最大速度為4500m/s，最小速度為1500m/s，25Hz主頻的雷克子波作為震源加在地表中心向下60個網(wǎng)格處(紅五星位置)，其他參數(shù)設置同均勻模型。

圖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地震波動正演模擬的特性，本文提出了一種基于多松弛參數(shù)的黏滯吸收邊界，該邊界算法簡單，易實現(xiàn)且可擴展性強?？紤]到衰減參數(shù)組合對邊界反射吸收效果影響顯著，因而通過數(shù)值實驗對其進行了優(yōu)化。

用均均勻模型、層狀模型、孔洞模型和復雜的BP模型進行了測試，結果表明本文的黏滯吸收邊界條件較好地壓制了邊界反射，且適用復雜模型地震波場模擬。

需要注意的是，本文衰減參數(shù)組合只嘗試了幾種冪函數(shù)，衰減系數(shù)以10的整數(shù)倍為間隔，因而得到的未必是最優(yōu)的參數(shù)組合。還可以通過嘗試新的衰減函數(shù)，比如Sigmoid函數(shù)、正弦函數(shù)、反正切函數(shù)等[23]，以及二次加密衰減系數(shù)的采樣間隔等方式進一步優(yōu)化。另外，為了簡單起見，本文實現(xiàn)邊界衰減的三個松弛參數(shù)采用的是同一套衰減體系，為了使邊界吸收效果更好，各個松弛參數(shù)可以嘗試不同的衰減過程，但是要注意在松弛參數(shù)變化的過程中要滿足MRT-LBM的穩(wěn)定性條件。最后，可以采用并行方式處理MRT-LBM計算量大的問題。