高國(guó)超 夏鵬 吉健維

基金項(xiàng)目：貴州大學(xué)引進(jìn)人才科研項(xiàng)目（貴大人基合字（2021）36號(hào)）；貴州省基礎(chǔ)研究（自然科學(xué)）項(xiàng)目（黔科合基礎(chǔ)-ZK[2024]一般088）；貴州省級(jí)地質(zhì)勘查資金項(xiàng)目（52000021MGQSE7S7K6PRP）

第一作者簡(jiǎn)介：高國(guó)超（1988-），男，博士，實(shí)驗(yàn)師。研究方向?yàn)榈卣鸩〝?shù)值模擬、工程物探。

DOI：10.19981/j.CN23-1581/G3.2024.13.002

摘? 要：數(shù)值仿真作為勘探地震學(xué)中研究地震波的重要手段，為野外實(shí)際數(shù)據(jù)采集和處理提供有力的科學(xué)依據(jù)。作為反射地震成像的重要補(bǔ)充，繞射波可以提高成像分辨率。因此，有效識(shí)別和確認(rèn)數(shù)據(jù)中繞射波成分顯得尤為重要。以正弦形粗糙界面為研究對(duì)象，采用譜元法數(shù)值模擬技術(shù)，分析繞射波傳播規(guī)律，確定繞射波空間分布特征。結(jié)果表明，地震波與正弦界面作用后產(chǎn)生不同階繞射波，其運(yùn)動(dòng)學(xué)特征滿足光柵方程，且高階繞射波具有更加嚴(yán)重的角頻散。此外，仿真結(jié)果對(duì)學(xué)生深入理解繞射波對(duì)地震數(shù)據(jù)處理的影響具有重要的輔助作用。

關(guān)鍵詞：繞射波；數(shù)值仿真；周期界面；譜元法；角頻散

中圖分類號(hào)：P631? ? ? 文獻(xiàn)標(biāo)志碼：A? ? ? ? ? 文章編號(hào)：2095-2945（2024）13-0005-04

Abstract： As an important means of studying seismic waves in exploration seismology， numerical simulation provides a powerful scientific basis for field data acquisition and processing. As an important supplement to reflection seismic imaging， diffracted waves can improve the imaging resolution. Therefore， it is particularly important to effectively identify and confirm the diffracted wave components in the data. Taking the sinusoidal rough interface as the research object， the spectral element method（SEM） is used to analyze the propagation law of diffracted waves and determine the spatial distribution characteristics of diffracted waves. The results show that different order diffraction waves are produced after the interaction between seismic waves and sinusoidal interfaces， and their kinematic characteristics satisfy the grating equation， and the higher order diffraction waves have more serious angular dispersion. In addition， the simulation results play an important auxiliary role in students' in-depth understanding of the influence of diffracted waves on seismic data processing.

Keywords： diffraction wave; numerical simulation; periodic interface; spectral element method （SEM）; angular dispersion

繞射波在地震勘探中越來(lái)越受到重視，得益于其提高反射地震成像分辨率[1]的能力，以及在高精度恢復(fù)構(gòu)造細(xì)節(jié)方面的巨大潛力[2]，尤其對(duì)于地下小尺度構(gòu)造，如斷層、裂縫、通道和粗糙界面。隨著勘探開(kāi)發(fā)的深入，這些微小地質(zhì)特征成為研究重點(diǎn)，因此除了鏡面反射數(shù)據(jù)，繞射波成為地震數(shù)據(jù)中最有價(jià)值的信息。然而，如果繞射波未得到恰當(dāng)處理，則會(huì)造成最終成像結(jié)果模糊，降低成像質(zhì)量。因此，有效識(shí)別并確認(rèn)地震數(shù)據(jù)中的繞射波成分至關(guān)重要。

由于地震勘探野外實(shí)際數(shù)據(jù)采集過(guò)程周期長(zhǎng)、成本高、存在安全風(fēng)險(xiǎn)、現(xiàn)場(chǎng)作業(yè)技術(shù)人員要求高，還會(huì)受工區(qū)地表環(huán)境和設(shè)備種類型號(hào)限制，致使地震實(shí)驗(yàn)在實(shí)驗(yàn)室或者野外開(kāi)展實(shí)物模擬極為困難，進(jìn)而無(wú)法滿足實(shí)驗(yàn)教學(xué)及科學(xué)研究的需求。數(shù)值仿真技術(shù)作為一種運(yùn)用計(jì)算機(jī)編程模仿實(shí)際環(huán)境進(jìn)行科學(xué)實(shí)驗(yàn)的技術(shù)，具有實(shí)驗(yàn)經(jīng)濟(jì)、安全、靈活、可視化和可重復(fù)性等優(yōu)點(diǎn)[3]，可以擺脫真實(shí)地震采集實(shí)驗(yàn)的約束。此外，仿真模擬還可以針對(duì)性地反饋野外實(shí)驗(yàn)設(shè)計(jì)的不足，有效合理地改進(jìn)實(shí)驗(yàn)方案，降低后期數(shù)據(jù)處理成本。因此，數(shù)值仿真技術(shù)已成為地震物探領(lǐng)域一種必不可少的手段。

近年來(lái)計(jì)算機(jī)技術(shù)的迅猛發(fā)展，為地震波仿真提供了有利條件。許多學(xué)者已利用數(shù)值仿真手段研究了地震波繞射現(xiàn)象：Rayleigh[4]開(kāi)創(chuàng)性地研究了周期性表面的繞射波現(xiàn)象。自此，該問(wèn)題得到廣泛關(guān)注[5]，但是大部分工作集中于反射系數(shù)的計(jì)算，繞射波的傳播規(guī)律和分布特征并未引起重視[6]。周期界面的繞射波問(wèn)題尤為重要，一個(gè)原因是在平面波入射情況下，界面周期性在除鏡面反射之外的某些方向上產(chǎn)生非常強(qiáng)的散射，影響處理效果；另一個(gè)原因是由傅里葉分析可知，任何復(fù)雜形態(tài)的界面都可以通過(guò)周期性結(jié)構(gòu)研究，為非周期界面波場(chǎng)分析提供研究借鑒。

本文針對(duì)周期性界面的地震繞射波問(wèn)題，以正弦形起伏界面為研究對(duì)象，采用譜元法數(shù)值仿真技術(shù)[7]，通過(guò)不同的震源組合和模型配置，研究地震波的傳播規(guī)律，識(shí)別地震數(shù)據(jù)中的繞射波成分，分析繞射波的空間分布特征，對(duì)揭示繞射波在地震成像中的作用具有重要意義，為進(jìn)一步研究隨機(jī)粗糙界面的繞射波特征提供思路。

1 周期性界面繞射波理論

單頻波與正弦形周期粗糙界面相互作用時(shí)，由于界面連續(xù)周期發(fā)射的小波相互增強(qiáng)，地震波會(huì)在鏡面反射（亦稱為0階繞射）以外的某些方向上產(chǎn)生非常強(qiáng)的繞射波[4]，如圖1所示。這些離散的方向稱為特征方向，與特征方向?qū)?yīng)的繞射波為繞射階。如果界面位于xy平面，其周期性沿x軸（圖1），且入射波位于xz平面，則繞射階的傳播角度（相對(duì)于粗糙界面的平均平面的法線）滿足下列繞射光柵方程[6]

式中：i是入射角，？茲n表示n階繞射波的角度，？姿是地震波波長(zhǎng)，且滿足？姿=？淄/f，？淄為介質(zhì)中的地震波速度，f為地震波的頻率。參數(shù)d代表周期性起伏界面的空間波長(zhǎng)。請(qǐng)注意，繞射階的符號(hào)由相對(duì)于0階繞射波（即鏡面反射波）的方向定義，如果繞射階位于相對(duì)于0階的順時(shí)針?lè)较?，則繞射階為正，否則，繞射階為負(fù)。

2? 數(shù)值仿真技術(shù)

譜元法（Spectral Element Method，SEM）作為有限元方法（Finite Element Method，F(xiàn)EM）的變種方法，采用了高階分段多項(xiàng)式逼近波動(dòng)方程的弱公式，既具有偽譜法的準(zhǔn)確性，又具有有限元法的靈活性[7]?；贕auss-Lobatto-Legendre數(shù)值積分和高階Lagrange插值，譜元法公式中的質(zhì)量矩陣呈現(xiàn)對(duì)角性，且時(shí)間遞推格式是全顯式的，因此該方法便于優(yōu)化存儲(chǔ)成本、易于通過(guò)并行方案高效實(shí)現(xiàn)。由于譜元法繼承了有限元法優(yōu)點(diǎn)，擅長(zhǎng)處理復(fù)雜幾何界面：粗糙的地形、傾斜或彎曲的界面，甚至扭曲的網(wǎng)格[8]，因此，為了準(zhǔn)確模擬粗糙界面情況下的地震波傳播，采用譜元法數(shù)值仿真技術(shù)更加合適。

圖1? 平面波與正弦形周期界面相互作用的示意圖

3? 數(shù)值仿真實(shí)驗(yàn)

在數(shù)值仿真過(guò)程中，使用帶寬為15 Hz的突發(fā)源，其主頻為100 Hz（圖2）。由其振幅譜可知，突發(fā)源信號(hào)是一個(gè)窄帶寬信號(hào)，其最小和最大頻率分別約為85 Hz和115 Hz。根據(jù)光柵方程（1），為了減少不同繞射階的重疊覆蓋，震源信號(hào)的帶寬應(yīng)盡可能窄。因此本文采用的窄帶寬突發(fā)源將有助于識(shí)別仿真結(jié)果中的各階繞射波。

圖2? 模擬仿真中采用的突發(fā)源

對(duì)于模型配置（圖3），研究底部具有周期性界面的聲學(xué)介質(zhì)，其聲波速度為1 500 m/s，因此主頻處對(duì)應(yīng)的波長(zhǎng)為？姿0=15 m。該模型的水平和垂直尺寸分別為3 000 m和1 605 m，底部為正弦形周期性粗糙界面，振幅為a=0.5？姿0=7.5 m，空間周期設(shè)置為d=■？姿0≈21.2 m。此外，模型底部周期性界面是壓力自由邊界，其他邊界設(shè)置為完美匹配層吸收邊界（PML）[9]。在使用點(diǎn)源進(jìn)行數(shù)值模擬的情況下，突發(fā)源放置在（1 500、105 m）處。為了更準(zhǔn)確地描述彎曲界面，模型網(wǎng)格剖分采用9個(gè)控制點(diǎn)定義四邊形譜元的幾何形狀。

圖3? 模型配置示意圖

3.1? 平面波入射-全起伏界面情況

為了觀察繞射階在空間上的離散分布特征，首先使用近平面波作為入射波，該平面波由751個(gè)突發(fā)源組合均勻排列在直線上并同時(shí)激發(fā)產(chǎn)生。設(shè)定入射角i為30°（相對(duì)于z軸），根據(jù)繞射光柵方程（1），理論上計(jì)算出帶寬內(nèi)（即85～115 Hz范圍）每個(gè)頻率的繞射階角度，結(jié)果見(jiàn)表1。

表1? 帶寬信號(hào)入射情況下不同階的繞射波角度

從表1中發(fā)現(xiàn)，在平面波入射情況下，除鏡面反射之外，只存在-1階和-2階繞射波。此外，對(duì)于-1階繞射波，由？茲n（fmax）-？茲n（fmin）定義的繞射角范圍約為13°（從-19.38°到-6.6°），而-2階繞射波約為43°（從-90°到-46.87°），因此-2階繞射波的角頻散現(xiàn)象比-1階繞射波更加嚴(yán)重。此外，根據(jù)光柵方程（1），圖4繪制了帶寬信號(hào)入射角 i與繞射角？茲n之間的關(guān)系。當(dāng)入射角僅為30°（圖5中黑色虛線表示）時(shí)，出現(xiàn)了0階、-1階和-2階繞射波，且-2階繞射角度范圍明顯大于-1階繞射波，這些結(jié)果與表1中結(jié)果一致。需要注意的是，由于模型尺寸限制了最大入射角度為45°，所以圖4中并不會(huì)出現(xiàn)-3階繞射波曲線。

圖4? 帶寬信號(hào)入射情況下入射角與繞射角之間的關(guān)系

為了驗(yàn)證上述理論結(jié)果，對(duì)該模型進(jìn)行了數(shù)值模擬實(shí)驗(yàn)。圖5展示了在t=1.32 s時(shí)刻的波場(chǎng)快照。從快照中，不僅可以觀察到0階（即鏡面反射）、-1階和-2階的繞射波，還可以直觀地看到-2階繞射波的分布范圍遠(yuǎn)大于-1階繞射波，該數(shù)值結(jié)果與表1所示的理論計(jì)算結(jié)果一致。同時(shí)發(fā)現(xiàn)-2階繞射波同相軸的彎曲曲率更大，這是由于在入射波為帶寬信號(hào)時(shí)，高階繞射波的角頻散效應(yīng)比低階繞射波更為嚴(yán)重。

注：其中直線代表751個(gè)突發(fā)源組合的位置。

圖5? t=1.32 s時(shí)刻的波場(chǎng)快照

3.2? 點(diǎn)源入射-半水平半起伏界面情況

為了仿真地震數(shù)據(jù)實(shí)際采集情形，對(duì)點(diǎn)源激發(fā)情況開(kāi)展了數(shù)值測(cè)試。與圖3所示模型不同，設(shè)計(jì)了一個(gè)特殊的模型：底部界面一半是平的，一半是正弦起伏的（圖6），其主要是為了避免來(lái)自負(fù)入射角的繞射波與來(lái)自正入射角的繞射波的干擾，從而更方便追蹤和識(shí)別不同階繞射波場(chǎng)。圖6是在t=1.74 s時(shí)刻的波場(chǎng)快照，一方面，觀察到了0階、+1階、-1階和-2階繞射波（圖6），與理論曲線（圖4）結(jié)果吻合。與平面波入射的情況相比（圖5），繞射波場(chǎng)形態(tài)差異很大，這是因?yàn)槭褂命c(diǎn)源時(shí)，入射角并不是固定不變，而是隨著入射界面位置變化。因此，對(duì)于帶寬信號(hào)的點(diǎn)源，多頻和多入射角共同影響了繞射波場(chǎng)的最終形態(tài)。另一方面，-1階與+1階繞射波相互干擾，且兩者緊隨0階繞射波，但由于傳播方向相反，-1階繞射波與+1階繞射波隨時(shí)間逐漸分離。另外，從該數(shù)值結(jié)果上可以看到由有限長(zhǎng)度粗糙界面的邊緣產(chǎn)生的虛假繞射波（圖6）。

圖6? t=1.74 s時(shí)刻的波場(chǎng)快照

3.3? 點(diǎn)源入射-全起伏界面情況

同樣是點(diǎn)源激發(fā)，對(duì)于底部界面是完全起伏情況，與正入射角入射波場(chǎng)類似，負(fù)入射角入射波場(chǎng)同樣產(chǎn)生了不同階繞射波，但由于點(diǎn)源兩側(cè)入射角具有對(duì)稱性，負(fù)入射角與正入射角產(chǎn)生的繞射波場(chǎng)本質(zhì)上是相同的，只是符號(hào)相反而已。與半水平-半起伏模型的波場(chǎng)快照（圖6）相比，完全起伏模型的負(fù)入射角與正入射角產(chǎn)生的各階繞射波在時(shí)空上相互覆蓋疊加，造成整個(gè)地震波場(chǎng)形態(tài)更加復(fù)雜，無(wú)法直觀辨別不同階繞射波場(chǎng)特征（圖7）。在對(duì)半水平-半起伏模型繞射波場(chǎng)特征分析的基礎(chǔ)上，發(fā)現(xiàn)完全起伏模型除了產(chǎn)生+1階、-1階以及-2階繞射波外，同時(shí)出現(xiàn)了+2階繞射波；由于負(fù)入射角的引入，-1階與+1階的繞射角范圍明顯比半水平-半起伏模型情況下增加。

圖8是地震波檢波器（其位置如圖7中直線所示）記錄的共炮點(diǎn)地震數(shù)據(jù)。基于繞射光柵方程（1）和模型介質(zhì)聲波速度，可以理論獲得不同偏移距下各階繞射波初至?xí)r間，結(jié)合圖8中模擬的共炮點(diǎn)數(shù)據(jù)發(fā)現(xiàn)-1階與+1階繞射波在時(shí)間上呈現(xiàn)出完全重疊的特點(diǎn)，-2階與+2階繞射波亦是如此，但是由于傳播方向相反，它們只在某些時(shí)刻重疊（如圖8中的-2與+2階繞射波），該結(jié)論與波場(chǎng)快照展示的結(jié)果（圖7）一致，如果傳播時(shí)間足夠長(zhǎng)，重疊的繞射波場(chǎng)在空間上將會(huì)相互獨(dú)立。需要指出的是，為避免邊界條件對(duì)數(shù)值結(jié)果的影響，采集長(zhǎng)度往往小于模型尺寸，致使模型范圍內(nèi)入射角在粗糙界面產(chǎn)生的繞射波并非都能夠被圖7中的檢波器記錄到。

注：其中檢波器記錄位置由圖7中直線所示。

圖8? 單炮地震數(shù)據(jù)

4? 結(jié)束語(yǔ)

本文通過(guò)不同的震源組合和模型配置，使用譜元法數(shù)值仿真技術(shù)，研究了正弦形周期性起伏界面對(duì)地震波傳播的影響，并借助相關(guān)理論計(jì)算和關(guān)系曲線，有效識(shí)別了波場(chǎng)快照和地震數(shù)據(jù)中不同階的繞射波，揭示了地震繞射波的空間分布形態(tài)。雖然實(shí)際激發(fā)的地震信號(hào)并不具有模擬中的窄帶寬特征，但研究結(jié)果為在地震成像中利用繞射波奠定了基礎(chǔ)，對(duì)進(jìn)一步探索更加復(fù)雜起伏界面-隨機(jī)起伏界面提供了思路。另一方面，用數(shù)值模擬仿真，可以直觀顯示繞射波的傳播規(guī)律和分布特點(diǎn)，對(duì)輔助學(xué)生深入理解繞射波對(duì)地震數(shù)據(jù)處理的影響具有重要意義。

參考文獻(xiàn)：

[1] 盛同杰，趙驚濤，彭蘇萍.地震繞射波弱信號(hào)U-net網(wǎng)絡(luò)提取方法[J].地球物理學(xué)報(bào)，2023，66（3）：1192-1204.

[2] 欒錫武，楊佳佳.地震繞射波波場(chǎng)分離與成像方法綜述[J].石油物探，2022，61（5）：761-770.

[3] 向世明，王傳雷，張學(xué)強(qiáng).計(jì)算機(jī)虛擬技術(shù)在物探仿真實(shí)驗(yàn)教學(xué)中的應(yīng)用及教學(xué)實(shí)踐[J].中國(guó)地質(zhì)教育，2011，20（4）：48-50.

[4] RAYLEIGH， J W S. The theory of sound[M]. Volume 2. New York： Dover Publications， 1945.

[5] 孫成禹，杜世通.平面聲波在粗糙界面上的反射特征研究[J].地球物理學(xué)報(bào)，2006，49（3）：903-907.

[6] GAO G， CRISTINI P， FAVRETTO C N， et al. On the reflection of time-domain acoustic spherical waves by a sinusoidal diffraction grating[J]. Acta Acustica， 2021（5）：6.

[7] 孟雪莉，劉少林，楊頂輝，等.基于優(yōu)化數(shù)值積分的譜元法模擬地震波傳播[J].石油地球物理勘探，2022，57（3）：602-612.

[8] 劉玉柱，劉偉剛，吳世林，等.基于譜元法的全波形反演及其在海底地震數(shù)據(jù)中的應(yīng)用[J].地球物理學(xué)報(bào)，2020，63（8）：3063-3077.

[9] 廉西猛，單聯(lián)瑜，隋志強(qiáng).地震正演數(shù)值模擬完全匹配層吸收邊界條件研究綜述[J].地球物理學(xué)進(jìn)展，2015，30（4）：1725-1733.