地震多級散射波正演模擬方法＊

張 博 吳國忱 李青陽 楊凌云 單俊臻

（中國山東青島 266580 中國石油大學（華東）地球科學與技術(shù)學院）

引言

廣義而言，一切由地球三維非均勻體擾動所引起地震波的變化均可稱為地震波散射；狹義上，地震波散射現(xiàn)象通常是指摒除廣義概念中可以用幾何光學（射線）理論表征的由大尺度非均勻體擾動引起的地震波旅行時和振幅變化的部分，由剩余地球三維非均勻體擾動所引起的地震波場畸變現(xiàn)象（吳如山等，1993；李燦蘋等，2005）.而實際地下地質(zhì)構(gòu)造十分復(fù)雜，各種尺度的非均質(zhì)體共同存在，地面檢波器接收到的往往是多波互相干涉的復(fù)雜地震波場信息，地震散射波作為地震波的一種，它是由入射波與地下的非均勻體相互作用而產(chǎn)生的波，因此地面接收到的散射波中攜帶了大量源于散射體的構(gòu)造和巖性信息，具有對地下縫、洞等小尺度非均質(zhì)體進行精細刻畫的潛力.而對裂縫-孔洞型等非均勻體散射體進行散射波場正演模擬，分析此類介質(zhì)的散射波場特征是利用散射波分析研究地下散射體的基礎(chǔ).

對散射波場的正演模擬是識別散射波的基礎(chǔ).Wu和Huang （1992）采用相位屏算子，模擬了二維垂直變背景情況下的散射場；Wu等（1995）引入De Wolf近似和相屏算子，計算了三維常背景情況下的背向散射場；符力耘等（1998）利用擾動理論，建立了用均勻介質(zhì)格林函數(shù)作為基本解的體積分方程，給出了配置法求解體積分方程的數(shù)值方法；田麗花（2007）采用相位移法在f-k域?qū)Φ卣鹕⑸洳▓鲞M行了波動方程正演模擬；劉鐵華（2012）設(shè)計了一種基于微擾論的f-k域積分法，在散射場的二次震源和空間能量衰減處理兩方面進行了改進，提高了計算精度和效率；Eaton和David （1999）采用波恩（Born）近似和射線理論近似方法，計算了格林函數(shù)的三維彈性波弱散射場；奚先（2018，2020）采用卷積神經(jīng)網(wǎng)絡(luò)（convolution neural network，縮寫為CNN）對散射波進行識別，大致識別出測試模型中散射點的準確位置及復(fù)雜偏移剖面中的各散射體的位置.

上述研究主要是對一次散射波進行模擬，本文則推導多級散射波和各級散射波微分方程，并基于此采用高階有限差分對能夠產(chǎn)生多級散射波的散射體模型進行正演模擬，分析多級散射波的地震響應(yīng)特征，以推斷地下散射體的分布情況和性質(zhì).

1 基礎(chǔ)理論

根據(jù)介質(zhì)分解理論，將地下介質(zhì)分解為背景介質(zhì)和擾動介質(zhì)（圖1）.震源激發(fā)產(chǎn)生的入射波在背景介質(zhì)中的波場稱為背景波場，背景介質(zhì)所對應(yīng)的速度稱為背景速度c0.入射波與擾動介質(zhì)相互作用所產(chǎn)生的波場稱為散射波場u（sr，ω），散射波場可以看作是觀測波場與背景波場之差（符力耘等，1998；雷蕾等，2011），即u（sr，ω）＝u（r，ω）－u（0r，ω）.

圖1 地震散射波產(chǎn)生原理示意圖Fig. 1 Schematic diagram of generating seismic scattered wave

波恩近似理論.波恩近似利用波場的振幅，當介質(zhì)產(chǎn)生的散射波場遠小于背景波場時，即u（sr，ω）?u（r，ω），用背景波場代替總波場將方程線性化，波恩近似適用于弱散射當圍巖和散射體的速度差與圍巖速度的比值ΔV＜15%的小擾動情況，或傳播距離短的情況（尹軍杰等，2005）.

地下介質(zhì)可以劃分為背景介質(zhì)和擾動介質(zhì)，在研究散射特征中，擾動介質(zhì)相較于周圍介質(zhì)其體積較小，地震波在穿越背景介質(zhì)遇到小擾動體時，散射體會作為一次源激發(fā)散射，相應(yīng)稱其為散射波.散射波響應(yīng)過程如圖2a所示，當散射體僅為一個時，稱之為一級散射，其中S1表示一次散射波.當?shù)叵陆橘|(zhì)散射體（擾動體）增多后，不同位置的散射體會作為二次源激發(fā)其它一次散射源產(chǎn)生的散射波，如圖2b所示，散射體V1激發(fā)的一次散射波在散射體V2處發(fā)生二次散射（Alkhalifah，2015，2016），稱為二級散射；而在散射體V2產(chǎn)生的二次散射傳播到散射體V1（也可以是散射體V3，為了方便描述未展示）處再次產(chǎn)生三次散射，稱之為三級散射.由此可知，當?shù)叵麓嬖诖罅可⑸潴w時，檢波器最終接收到的散射波信息數(shù)量龐大，非常復(fù)雜，是地震波在各個散射體之間錯綜復(fù)雜的傳播和被激發(fā)所導致的.因此，為了研究這些散射波的特征，將其進行分級處理，散射波對應(yīng)的級數(shù)即為入射波經(jīng)過幾次散射體后所產(chǎn)生的.

圖2 多級散射原理示意圖（a）一級散射；（b）二級散射；（c）三級散射Fig. 2 Schematic diagram of multi-order scattering principle（a）The first-order scattering；（b）The second-order scattering；（c）The third-order scattering

2 地震多級散射波方程推導

2.1 基于波恩近似的散射波場分級

李普曼-施溫格（Lippmann-Swinger）方程（Wu，Huang，1992；Wuet al，1995）的總波場計算公式為：

式中： ω為角頻率，k0＝ω/C0為背景介質(zhì)中波數(shù)，r＝（x，z）為 散射點位置，r1＝（x1，z1）為積分空間的任一點源，G為背景介質(zhì)中的格林函數(shù)， ε（r1）為背景速度上疊加的一個擾動量.當u（sr，ω）≤u（r，ω）時 ，上式右邊被積函數(shù)中的u用u0代替，則波恩一級近似表達式為（Born，Wolf，1999）

若用u1代 替右邊被積函數(shù)中u，則可得到波恩二級近似表達式

將式（2）帶入式（3）可得：

式中：等式右端第二項積分是關(guān)于散射體的一級散射，第三項是關(guān)于散射體的二級散射，依次進行替換迭代，每項與前一項的迭代關(guān)系式為

為方便表示，多級散射的波恩近似表達式則為：

即

則有

式中，上標1，2，···，N為散射波場級數(shù).

2.2 多級散射波方程推導

在均勻介質(zhì)中，二維標量波方程（無震源項）為

式中：C（r）為 波在均勻介質(zhì)中的傳播速度，為拉普拉斯算子，t為時間.為方便推導，將式（9）轉(zhuǎn)換至頻率域，表示為

背景波場滿足標量波方程

將式（8）和式（11）帶入式（10），可得到頻率域多級散射波方程

隨后，根據(jù)介質(zhì)擾動原理及波恩近似理論，對多級散射波方程進行分級處理.

以此類推，可以得到，每一級獨立的散射波方程

式中，N（N≥1）表 示散射波的級數(shù)，當N＝1時，等式右端為背景場.

將式（16）轉(zhuǎn)換到時間域，得到時間域各級散射波方程：

式中等號右端整體可視為虛震源，其中二階導數(shù)項?2uN－1/?t2無法直接計算獲得，本文以空間導數(shù)來替代時間導數(shù)項.

2.3 多級散射波方程的有限差分離散

采用規(guī)則網(wǎng)格有限差分法進行離散，時間二階導數(shù)的離散格式為

式中，u為壓力場， Δt為時間步長.空間二階導數(shù)的空間2N階離散格式為

式中, ωn表示二階導數(shù)的 2N階精度差分系數(shù)（趙茂強，2010；劉慶敏，2007）.規(guī)則網(wǎng)格有限差分方程的穩(wěn)定性條件為（梁展源，2016；李青陽等，2018）

式中，C（x，z）為 速度， Δt為時間步長， Δx和Δz為空間步長， ωn表示差分系數(shù).最后，加入完美匹配層（Perfectly matched layer，縮寫為PML）邊界條件（何燕，2008；吳國忱等，2014，2020），對多級散射波方程采用有限差分法進行正演模擬.

3 數(shù)值模擬結(jié)果及分析

3.1 兩點強散射體模型

為驗證多級和各級散射波方程的正演模擬精度，設(shè)計一簡單兩點散射模型（圖3），模型大小為2 000 m×2 000 m，背景速度為2 500 m/s，在深度1 000 m和1 400 m處分別設(shè)置一個點散射體，散射體的速度為2 000 m/s，縱、橫向網(wǎng)格間距均為10 m，震源采用主頻為25 Hz的雷克子波，炮點位于地表1 000 m處，采樣時間點數(shù)為2 000，采樣間隔為0.8 ms，全排列接收.

圖3 兩點強散射體模型Fig. 3 A model with two strong scatterers

兩點強散射模型對應(yīng)的三個時刻的多級散射波波場和每一級獨立的散射波波場如圖4所示.圖4a中二級散射波和三級散射波的能量弱于一級散射波，為方便展示，對t＝960 ms和t＝1 200 ms的波場快照添加增益，因此整體未加色標，圖4b?d可直觀地看出各級散射波方程所模擬的波場快照僅包含自身級數(shù)的散射波，不包含其它級的散射干擾.各級散射波的傳播均符合惠更斯-菲涅爾（Huygens-Fresnel）原理，且隨著散射波的級數(shù)增加，其能量衰減.

圖4 不同級數(shù)的散射波波場快照（a）多級散射波；（b）一級散射波；（c）二級散射波；（d）三級散射波Fig. 4 Snapshots of scattered wavefields at different orders（a）The multi-order scattering；（b）The first-order scattering；（c）The second-order scatterings；（d）The third-order scattering

為討論各級散射波方程穩(wěn)定性，抽取各級散射波第100道記錄與參考地震記錄進行對比，結(jié)果如圖5所示.圖中參考記錄為觀測記錄和背景速度正演模擬記錄之差，多級散射記錄為本文方程正演所得記錄，可以看出：兩條記錄曲線在相位和振幅上完全吻合，驗證了多級散射波方程的準確性.各級散射記錄與參考記錄第100道的波形對比如圖6a所示，其中，t為0.832—0.960 s時，參考記錄中的一級散射信息與一級散射記錄的相位信息一致，振幅信息存在較小的誤差，該誤差主要是速度模型中擾動體與圍巖的速度差未滿足波恩弱散射近似（圖6b）；而t為1.240—1.360 s時的誤差與一級散射信息的誤差來源相同（圖6c）.對擾動體與圍巖的速度差需滿足波恩弱散射近似的速度模型做出如下分析.

圖5 多級散射記錄與參考記錄第100道的波形對比Fig. 5 Comparison of the 100th record of the multi-order scattering records with that of the reference records

圖6 各級散射記錄與參考記錄第100道的波形對比Fig. 6 Comparison of the 100th record of each scattering records and the reference records（a）t＝0—1.6 s；（b）t＝0.832—0.960 s；（c）t＝1.240—1.360 s

3.2 兩點弱散射體模型

設(shè)計兩點弱散射模型，設(shè)散射體的速度為2 400 m/s，其余參數(shù)均與兩點強散射模型相同.兩點弱散射模型如圖7a所示，其對應(yīng)的多級散射記錄與參考記錄第100道的波形對比如圖7b所示，由該圖可見：兩條記錄曲線完全吻合，驗證了多級散射波方程的準確性.兩點弱散射模型對應(yīng)的各級散射記錄與參考記錄第100道的波形對比如圖8所示，當t為0.848—0.944 s時，參考記錄中的一級散射信息與一級散射記錄的相位和振幅信息基本一致，驗證一級散射波的準確性（圖8b）；t為1.296—1.318 s時，參考記錄中的二級散射信息與二級散射記錄的相位和振幅信息基本一致，驗證二級散射波的準確性（圖8c）.由于散射的波的能量逐級遞減，三級散射波的能量十分微弱，在單道記錄上振幅信息不明顯，此處不作討論.

圖7 兩點弱散射體模型（a）及多級散射記錄與參考記錄第 100 道的波形對比（b）Fig. 7 A model with two weak scatterers （a）and comparison of the 100th record of the multiorderscattering records and with that of the reference records （b）

圖8 各級散射記錄與參考記錄第100道的波形對比Fig. 8 Comparison of the 100th record of each scattering records and the reference records（a）t＝0—1.6 s；（b）t＝0.848—0.944 s；（c）t＝1.296—1.318 s

3.3 復(fù)雜散射體模型

為了驗證本文方法對復(fù)雜模型的適用性和穩(wěn)定性，設(shè)計一個同時包含不同傾斜角度裂縫和不同尺度孔洞的復(fù)雜散射體模型，如圖9所示.模型大小2 000 m×2 000 m，背景速度為2 500 m/s，三條裂縫從左到右的傾角依次為90°，0°和45°，速度均為2 300 m/s，裂縫下面四個不同尺度的溶洞速度為2 400 m/s，縱、橫向網(wǎng)格間距均為10 m，震源采用主頻為25 Hz的雷克子波，炮點位于地表1 000 m處，采樣時間點數(shù)為2 000，采樣間隔為0.8 ms，全排列接收.

圖9 復(fù)雜散射體模型Fig. 9 A model with complex scatterers

復(fù)雜散射體模型對應(yīng)的多級散射波地震記錄如圖10所示，由該圖可見：多級散射記錄無直達波，散射波呈現(xiàn)雙曲形態(tài)，且多個散射體同時存在且距離相近時，各級散射波會產(chǎn)生干涉；同一尺度不同傾角的情況下，橫向裂縫的散射能量強于縱向裂縫；當孔洞的尺度小于1/4波長時，炮記錄表現(xiàn)為單條雙曲線，隨著孔洞散射體尺度的增大，其頂?shù)讓?yīng)的兩條雙曲線能量區(qū)分愈加明顯，直至完全分離.

圖10 多級散射記錄Fig. 10 Multi-order scattering records

分別抽取復(fù)雜散射體模型對應(yīng)的多級散射波地震記錄與參考地震記錄的第100道記錄進行對比，如圖11所示，圖中兩條曲線在振幅和相位上完全吻合，表明多級散射波方程對復(fù)雜模型具有適用性，并驗證了多級散射波方程的準確性.

圖11 多級散射記錄與參考記錄的第100道波形對比Fig. 11 Comparison of the 100th record of the multi-order scattering records with that of and the reference records

4 結(jié)論

針對非均勻散射體中地震波傳播正演模擬問題，本文推導了多級和各級散射波方程.通過算法分析和模型測試得到以下認識：

1）本文推導了多級和各級散射波方程，通過有限差分的正演數(shù)值模擬，抽取其中心道與參考記錄中心道進行對比，結(jié)果顯示，多級散射單道記錄與參考單道記錄完全吻合，驗證了多級散射波方程的準確性，各級散射記錄在波恩弱散射近似范圍內(nèi)，其單道記錄與參考單道記錄的振幅及相位基本吻合，驗證了在波恩弱散射近似下各級散射方程的準確性.

2）各級散射波的傳播符合惠更斯-菲涅爾原理，且隨著散射波級數(shù)的增加，能量隨之衰減.

3）在散射能量方面，尺度相同的情況下，橫向裂縫比縱向裂縫能量更強，散射體與圍巖的速度差異越大，其能量也越大.

本文基于標量波方程的散射波方程推導思路理論上能夠推廣至各向異性介質(zhì)，但具體推導內(nèi)容仍需進一步分析；本文采用的背景模型基本為半空間模型，而對于復(fù)雜背景介質(zhì)下的孔洞、縫洞的散射波識別、提取仍需進一步研究.

審稿專家提出有益的修改意見，作者在此表示感謝.