周期性粗糙界面地震繞射波數值仿真與分析

高國超 夏鵬 吉健維

基金項目：貴州大學引進人才科研項目（貴大人基合字（2021）36號）；貴州省基礎研究（自然科學）項目（黔科合基礎-ZK[2024]一般088）；貴州省級地質勘查資金項目（52000021MGQSE7S7K6PRP）

第一作者簡介：高國超（1988-），男，博士，實驗師。研究方向為地震波數值模擬、工程物探。

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

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

關鍵詞：繞射波；數值仿真；周期界面；譜元法；角頻散

中圖分類號：P631? ? ? 文獻標志碼：A? ? ? ? ? 文章編號：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

繞射波在地震勘探中越來越受到重視，得益于其提高反射地震成像分辨率[1]的能力，以及在高精度恢復構造細節方面的巨大潛力[2]，尤其對于地下小尺度構造，如斷層、裂縫、通道和粗糙界面。隨著勘探開發的深入，這些微小地質特征成為研究重點，因此除了鏡面反射數據，繞射波成為地震數據中最有價值的信息。然而，如果繞射波未得到恰當處理，則會造成最終成像結果模糊，降低成像質量。因此，有效識別并確認地震數據中的繞射波成分至關重要。

由于地震勘探野外實際數據采集過程周期長、成本高、存在安全風險、現場作業技術人員要求高，還會受工區地表環境和設備種類型號限制，致使地震實驗在實驗室或者野外開展實物模擬極為困難，進而無法滿足實驗教學及科學研究的需求。數值仿真技術作為一種運用計算機編程模仿實際環境進行科學實驗的技術，具有實驗經濟、安全、靈活、可視化和可重復性等優點[3]，可以擺脫真實地震采集實驗的約束。此外，仿真模擬還可以針對性地反饋野外實驗設計的不足，有效合理地改進實驗方案，降低后期數據處理成本。因此，數值仿真技術已成為地震物探領域一種必不可少的手段。

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

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

1 周期性界面繞射波理論

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

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

2? 數值仿真技術

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

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

3? 數值仿真實驗

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

圖2? 模擬仿真中采用的突發源

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

圖3? 模型配置示意圖

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

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

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

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

圖4? 帶寬信號入射情況下入射角與繞射角之間的關系

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

注：其中直線代表751個突發源組合的位置。

圖5? t=1.32 s時刻的波場快照

3.2? 點源入射-半水平半起伏界面情況

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

圖6? t=1.74 s時刻的波場快照

3.3? 點源入射-全起伏界面情況

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

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

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

圖8? 單炮地震數據

4? 結束語

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

參考文獻：

[1] 盛同杰，趙驚濤，彭蘇萍.地震繞射波弱信號U-net網絡提取方法[J].地球物理學報，2023，66（3）：1192-1204.

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

[3] 向世明，王傳雷，張學強.計算機虛擬技術在物探仿真實驗教學中的應用及教學實踐[J].中國地質教育，2011，20（4）：48-50.

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

[5] 孫成禹，杜世通.平面聲波在粗糙界面上的反射特征研究[J].地球物理學報，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] 孟雪莉，劉少林，楊頂輝，等.基于優化數值積分的譜元法模擬地震波傳播[J].石油地球物理勘探，2022，57（3）：602-612.

[8] 劉玉柱，劉偉剛，吳世林，等.基于譜元法的全波形反演及其在海底地震數據中的應用[J].地球物理學報，2020，63（8）：3063-3077.

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