張 博 吳國(guó)忱 李青陽(yáng) 楊凌云 單俊臻

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

引言

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

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

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

1 基礎(chǔ)理論

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

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

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

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

圖2 多級(jí)散射原理示意圖（a）一級(jí)散射；（b）二級(jí)散射；（c）三級(jí)散射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 地震多級(jí)散射波方程推導(dǎo)

2.1 基于波恩近似的散射波場(chǎng)分級(jí)

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

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

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

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

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

為方便表示，多級(jí)散射的波恩近似表達(dá)式則為：

即

則有

式中，上標(biāo)1，2，···，N為散射波場(chǎng)級(jí)數(shù).

2.2 多級(jí)散射波方程推導(dǎo)

在均勻介質(zhì)中，二維標(biāo)量波方程（無(wú)震源項(xiàng)）為

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

背景波場(chǎng)滿(mǎn)足標(biāo)量波方程

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

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

以此類(lèi)推，可以得到，每一級(jí)獨(dú)立的散射波方程

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

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

式中等號(hào)右端整體可視為虛震源，其中二階導(dǎo)數(shù)項(xiàng)?2uN－1/?t2無(wú)法直接計(jì)算獲得，本文以空間導(dǎo)數(shù)來(lái)替代時(shí)間導(dǎo)數(shù)項(xiàng).

2.3 多級(jí)散射波方程的有限差分離散

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

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

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

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

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

3.1 兩點(diǎn)強(qiáng)散射體模型

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

圖3 兩點(diǎn)強(qiáng)散射體模型Fig. 3 A model with two strong scatterers

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

圖4 不同級(jí)數(shù)的散射波波場(chǎng)快照（a）多級(jí)散射波；（b）一級(jí)散射波；（c）二級(jí)散射波；（d）三級(jí)散射波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

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

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

圖6 各級(jí)散射記錄與參考記錄第100道的波形對(duì)比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 兩點(diǎn)弱散射體模型

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

圖7 兩點(diǎn)弱散射體模型（a）及多級(jí)散射記錄與參考記錄第 100 道的波形對(duì)比（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 各級(jí)散射記錄與參考記錄第100道的波形對(duì)比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ù)雜散射體模型

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

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

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

圖10 多級(jí)散射記錄Fig. 10 Multi-order scattering records

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

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

4 結(jié)論

針對(duì)非均勻散射體中地震波傳播正演模擬問(wèn)題，本文推導(dǎo)了多級(jí)和各級(jí)散射波方程.通過(guò)算法分析和模型測(cè)試得到以下認(rèn)識(shí)：

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

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

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

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

審稿專(zhuān)家提出有益的修改意見(jiàn)，作者在此表示感謝.