沈天晶，胡葉正

（西南石油大學 地球科學與技術學院，四川成都 610500）

地震正演是在已知地下巖層空間結構與介質參數(shù)的情況下，通過數(shù)值模擬來模擬地震波的傳播。模擬結果中包括了走時信息與波場信息，反映了地下結構的空間規(guī)律與介質信息。正演數(shù)值模擬不僅是地球物理的基礎，也是其重要組成部分。通過對正演模擬的研究發(fā)現(xiàn)，一階的聲波波動方程的有限差分模擬是最常用的方法[1-2]。盡管許多學者對其求解方法與模型離散進行了優(yōu)化，但也增加了計算量，降低了計算效率?？紤]到上述問題，本文利用一階聲波波動方程進行正演模擬，同時利用交錯網(wǎng)格對其進行求解，利用CPML進行邊界處理。在計算效率的提升方面，本文利用了OpenMP與mpi對程序進行并行化。

1 CPML吸收邊界

利用交錯網(wǎng)格對一階聲波波動方程進行求解，將質點震動速度，波場應力與模型屬性定義在網(wǎng)格上。其中比較非常關鍵的一點就是選取了CPML吸收邊界，原理如下：

二階CPML吸收邊界的吸收因子為：

在時域中，在拉伸坐標中的偏導數(shù)的求取公式如下：

上述公式的卷積運算可以利用遞歸卷積法進行運算，公式如下：

式中，n為時間步長，ai與bi為系數(shù)，由下式給出：

由上述公式，二維二階聲波方程在拉伸坐標下的公式為：

由公式（3）可知：

將公式（7）與公式（3）帶入公式（6），可得：

公式（7）即為二階聲波方程的CPML吸收邊界公式。

一階的CPML吸收邊界的拉伸因子為：

i為計算點到吸收邊界內的距離，L為吸收層厚度，vmax為最大聲波速度，αmax是與地震波主頻有關的系數(shù)。

進一步推導可得出：

上式即為一階聲波方程CPML吸收邊界的公式。

2 模型實驗

通過Marmousi模型實驗來驗證以上方法，利用該模型的一塊區(qū)域，其中包括了斷層與細層，模型較為復雜。首先，利用二階的聲波方程進行正演模擬，震源主頻為30Hz，空間采樣為5m，時間采樣為0.000 5s，CPML吸收邊界厚度為16，差分精度為8，得到結果如圖1（a）所示，因為模型較為復雜，反射波雜亂，分辨率較低，箭頭位置甚至有頻散的出現(xiàn)。然后在所有參數(shù)不變的情況下，利用一階的聲波方程，利用交錯網(wǎng)格對其求解實現(xiàn)正演模擬，圖1（b）為質點震動速度的x方向的分量。圖1（c）為質點震動速度y方向的分量。這兩種結果模擬了野外實際采集數(shù)據(jù)的方式，沿著水平與豎直的方向分開采集。圖1（d）代表了合成的炮記錄，主要反映了縱波的傳播方式?？梢钥吹较啾榷A的常網(wǎng)格方程的結果，這種方法的精度更高，不容易頻散，更能夠模擬野外實際情況的采集。同時，使用了交錯網(wǎng)格，不僅沒有使計算量增加，還讓結果更加精確。

圖1 正演模擬結果

3 結論

利用交錯網(wǎng)格與一階波動方程進行數(shù)值模擬，并且利用Marmousi模型進行了模型實驗。同時對二階波動方程與一階波動方程的模擬效果的比較可以看到，一階波動方程模擬效果更好，精度更高，更加契合實際的采集數(shù)據(jù)。同時，一階的波動方程更穩(wěn)定，不容易頻散。利用交錯網(wǎng)格，在不提高計算量的情況下，使模擬更加精確。