起伏地表?xiàng)l件下頻率-空間域聲波介質(zhì)正演模擬

郭振波，李振春

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

起伏地表?xiàng)l件下頻率-空間域聲波介質(zhì)正演模擬

郭振波，李振春

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

頻率-空間域正演模擬是頻率域及Laplace-Fourier域全波形反演的基礎(chǔ)，起伏地表?xiàng)l件下波形反演算法的關(guān)鍵是正演算法中考慮起伏地表的影響。基于帶PML吸收邊界的聲波波動(dòng)方程，在已有最優(yōu)9點(diǎn)有限差分正演算法的基礎(chǔ)上構(gòu)建了起伏地表?xiàng)l件下頻率-空間域正演算法。通過(guò)應(yīng)用變網(wǎng)格技術(shù)，進(jìn)一步提高算法的計(jì)算效率、降低內(nèi)存開(kāi)銷(xiāo)，使得大規(guī)模起伏地表模型的頻率域正反演問(wèn)題成為可能。理論分析及數(shù)值測(cè)試表明：通過(guò)對(duì)近地表區(qū)域進(jìn)行局部網(wǎng)格加密，可有效地壓制由于矩形網(wǎng)格離散引起的角點(diǎn)散射；結(jié)合變網(wǎng)格技術(shù)可較易獲得5倍以上計(jì)算效率的提高及內(nèi)存占用的降低，且隨著模型尺度的增加及地表起伏高程差的減小，倍數(shù)將顯著增加；在細(xì)網(wǎng)格與粗網(wǎng)格交界處產(chǎn)生的虛假反射振幅幅值控制在原始波場(chǎng)的2%以內(nèi)，滿足地震波場(chǎng)正反演的需求。

起伏地表；變網(wǎng)格；頻率域；數(shù)值正演

0 引言

隨著計(jì)算能力的提高及對(duì)勘探目標(biāo)刻畫(huà)精度要求的提高，全波形反演正逐步成為一種高精度成像工具。全波形反演基于原始的雙程波動(dòng)方程，且能夠同時(shí)考慮地震數(shù)據(jù)的各種波形信息，因此，該方法能夠同時(shí)對(duì)模型的各個(gè)波數(shù)成分進(jìn)行成像，成像精度較高。目前主要的波形反演方法集中在時(shí)間域、頻率域、Laplace-Fourier域、Laplace域等，其中頻率域與Laplace-Fourier域算法由于具有計(jì)算效率較高、較易實(shí)現(xiàn)多尺度反演、更容易估算震源子波等優(yōu)勢(shì)，在科學(xué)研究及實(shí)際應(yīng)用中受到更多的關(guān)注[1-3]。高精度的頻率-空間域正演算法是頻率域與Laplace-Fourier域波形反演算法的基礎(chǔ)，正演算法的效率及精度將直接影響反演的結(jié)果。目前主要的頻率域正反演算法大都基于水平地表，且較少考慮自由邊界，在處理實(shí)際資料時(shí)所用算法受到限制。波形反演算法能否直接處理起伏地表?xiàng)l件下采集數(shù)據(jù)的關(guān)鍵是所用正演算法是否考慮起伏地表的影響，因此，考慮起伏地表的頻率域正演算法對(duì)于頻率域與Laplace-Fourier域波形反演具有重要意義。此外，頻率域正演算法與時(shí)間域算法相比數(shù)值實(shí)現(xiàn)及理論分析存在較大差異，同時(shí)還具有更易引入介質(zhì)的黏滯性效應(yīng)、對(duì)于多源正演效率較高、PML(完全匹配層)類(lèi)吸收邊界更易處理等優(yōu)勢(shì)。因此，考慮起伏地表的頻率域正演算法在繼承這些優(yōu)勢(shì)的同時(shí)可對(duì)起伏地表模型進(jìn)行直接數(shù)值模擬，充分發(fā)揮正演算法的優(yōu)勢(shì)。綜上所述，考慮起伏地表的高精度頻率-空間域正演算法的研究不管是對(duì)于正演算法本身,還是隨后的全波形反演方法都具有重要意義。

考慮起伏地表的正演方法較多，但各有優(yōu)缺點(diǎn)。在處理起伏地表的正演方法中，可大致分為基于不規(guī)則網(wǎng)格的正演方法和基于規(guī)則網(wǎng)格的正演方法2類(lèi)。基于不規(guī)則網(wǎng)格的正演方法為當(dāng)前研究熱點(diǎn)，主要包括有限元法(譜元法[4]、間斷Galerkin法[5]等)、有限體積法[6]、基于曲網(wǎng)格[7-9]或非規(guī)則網(wǎng)格[10-11]的有限差分法等。對(duì)于有限元和有限體積類(lèi)方法，由于求解波動(dòng)方程的弱形式解能夠自動(dòng)滿足自由邊界條件，所以此類(lèi)方法計(jì)算精度較高，但計(jì)算量一般較大；對(duì)于有限差分類(lèi)方法，自由邊界條件需要顯示設(shè)置，計(jì)算精度相對(duì)于有限元類(lèi)方法較低，但計(jì)算效率較高?；诓灰?guī)則網(wǎng)格的正演方法嚴(yán)重依賴網(wǎng)格剖分的精度，由于速度模型在反演過(guò)程中不斷更新，若沒(méi)有一種有效的、實(shí)時(shí)的網(wǎng)格剖分方法，該類(lèi)方法很難有效地應(yīng)用到波形反演算法中。相對(duì)而言，基于規(guī)則網(wǎng)格的正演方法較易移植到波形反演算法中。然而，利用規(guī)則矩形網(wǎng)格并不能有效地描述起伏地表，階梯狀離散易產(chǎn)生角點(diǎn)散射。Jih等[12]通過(guò)將起伏地表進(jìn)行折線近似來(lái)壓制角點(diǎn)散射，但該方法較難擴(kuò)展到三維；隨后Robertsson等[13]提出了通過(guò)在近地表區(qū)域加密網(wǎng)格來(lái)壓制角點(diǎn)散射，Hayshi等[14]、Wang等[15]對(duì)該方法進(jìn)行了發(fā)展。雖然起伏地表?xiàng)l件下的正演算法較多，但是大都在時(shí)間域?qū)崿F(xiàn)，頻率域中算法很少。

在頻率-空間域正演方法方面，研究重點(diǎn)大都集中在差分精度的提高上。Jo等[16]結(jié)合常規(guī)網(wǎng)格與旋轉(zhuǎn)網(wǎng)格各自的優(yōu)勢(shì)，利用最小平方的思想求取了頻率-空間域最優(yōu)9點(diǎn)格式的差分系數(shù)；隨后tekl等[17]將其推廣到彈性介質(zhì)正演；Shin等[18]發(fā)展了25點(diǎn)格式的頻域正演算法；Hustedt等[19]詳細(xì)分析并對(duì)比了最優(yōu)9點(diǎn)格式與四階差分格式；最近，曹書(shū)紅等[20]推導(dǎo)了聲波介質(zhì)17點(diǎn)格式正演算法；吳國(guó)忱、梁鍇等[21-22]推導(dǎo)并實(shí)現(xiàn)了TTI介質(zhì)、VTI介質(zhì)中頻率正演的相關(guān)算法；韓立等[23]討論了頻率域正演中PML邊界加載方法及頻率域變網(wǎng)格步長(zhǎng)計(jì)算策略。在頻率域起伏地表處理方面，Brossier等[24]利用有限體積法對(duì)起伏地表?xiàng)l件下彈性波進(jìn)行正演模擬，并隨后將其應(yīng)用到波形反演中。目前在頻率域有限差分正演模擬中，關(guān)于起伏地表的研究較少。

筆者沿用經(jīng)典的最優(yōu)9點(diǎn)格式有限差分算法，在處理起伏地表時(shí)將時(shí)域算法中變網(wǎng)格的思想應(yīng)用到頻率域正演中；在近地表區(qū)域?qū)W(wǎng)格進(jìn)行加密處理，而其余部分應(yīng)用常規(guī)粗網(wǎng)格，在保證正演精度的前提下提高了計(jì)算效率、降低了內(nèi)存開(kāi)銷(xiāo)。除自由邊界外，其余邊界均采用PML類(lèi)吸收邊界條件。筆者首先概述最優(yōu)9點(diǎn)格式頻率-空間域正演算法的基本原理，然后給出起伏地表?xiàng)l件下自由邊界條件及變網(wǎng)格技術(shù)的數(shù)值實(shí)現(xiàn)方法，最后通過(guò)數(shù)值實(shí)例證明方法的精度及對(duì)效率提高方面的優(yōu)勢(shì)。

1 頻率-空間域最優(yōu)9點(diǎn)格式有限差分正演模擬

頻率-空間域帶PML邊界條件的二維聲波波動(dòng)方程為

(1)

其中：p為波場(chǎng)；ρ為密度；K為體積模量；f為震源項(xiàng)；s為PML吸收邊界相關(guān)參數(shù)；下角x、z分別表示x軸與z軸。s的表達(dá)式為

(2)

其中：ω為角頻率；d為吸收因子，主要影響PML層內(nèi)波場(chǎng)振幅的衰減；α為頻移因子，主要影響低頻波的吸收；β為尺度因子，主要影響邊界大角度入射波的吸收。當(dāng)s=1時(shí)，沒(méi)有吸收效應(yīng)，波動(dòng)方程退化為原始波動(dòng)方程；當(dāng)β=1、α=0時(shí)，吸收邊界退化為常規(guī)PML邊界條件；當(dāng)β>1、α>0時(shí)，為CFS-PML(復(fù)頻移完全匹配層)邊界條件。根據(jù)不同的情況與要求，可對(duì)相應(yīng)參數(shù)進(jìn)行人為調(diào)整[25-26]。由于PML類(lèi)吸收邊界條件的定義及相關(guān)理論推導(dǎo)是在頻率域給出的，所以相對(duì)于時(shí)域正演算法，在頻率域?qū)崿F(xiàn)更為簡(jiǎn)單。

圖1 常規(guī)最優(yōu)9點(diǎn)差分格式(a)與網(wǎng)格離散方式(b)Fig. 1 Conventional optimal 9-point finite difference scheme (a) and discrete scheme (b)

在常規(guī)坐標(biāo)系(圖1a中x-z坐標(biāo)系)中，有限差分?jǐn)?shù)值頻散在波傳播方向與x(或z)軸夾角為0°時(shí)最大、45°時(shí)最小，而在旋轉(zhuǎn)坐標(biāo)系(圖1a中x′-z′坐標(biāo)系)中，數(shù)值頻散在夾角為0°時(shí)最小、45°時(shí)最大。最優(yōu)9點(diǎn)差分格式充分利用2種坐標(biāo)系下數(shù)值頻散的互補(bǔ)性，通過(guò)將2種坐標(biāo)系下的空間求導(dǎo)項(xiàng)進(jìn)行加權(quán)平均以減弱數(shù)值頻散的各向異性;然后借鑒頻率域有限元方法中將質(zhì)量加速度項(xiàng)也用周?chē)鷰讉€(gè)點(diǎn)加權(quán)平均表示的方法來(lái)提高正演的精度;最后利用最優(yōu)化方法求取涉及到的加權(quán)系數(shù)以構(gòu)建高精度頻率域正演算法。詳細(xì)的推導(dǎo)及分析可見(jiàn)參考文獻(xiàn)[16,19]。

將波場(chǎng)及介質(zhì)參數(shù)進(jìn)行空間離散，并用最優(yōu)9點(diǎn)差分格式代替空間求導(dǎo)項(xiàng)，可得原始波動(dòng)方程的離散形式：

(3)

其中：pi,j為波場(chǎng)在離散網(wǎng)格點(diǎn)(i,j)上的波場(chǎng)值，如圖1b所示；Ci與Ri為由介質(zhì)參數(shù)、頻率、差分系數(shù)、PML層參數(shù)、網(wǎng)格間距等共同控制的在相應(yīng)點(diǎn)的加權(quán)系數(shù)；fi,j為震源項(xiàng)。與時(shí)域算法不同，在數(shù)值求解的過(guò)程中需要將笛卡爾坐標(biāo)系下介質(zhì)參數(shù)與波場(chǎng)映射到線性坐標(biāo)系下以便將其寫(xiě)成向量的形式，如圖1b中數(shù)字編號(hào)所示。此時(shí)，離散后的波動(dòng)方程可以整理成矩陣的形式：

(4)

其中：A為阻抗矩陣；P為波場(chǎng)向量；F為震源項(xiàng)向量。

通過(guò)空間離散及有限差分近似，將原始偏微分方程的求解問(wèn)題轉(zhuǎn)化成為一稀疏線性方程組的求解問(wèn)題。若模型離散后總點(diǎn)數(shù)為M，則P為大小為M×1的向量，A為大小為M×M的矩陣。式(4)可以通過(guò)直接法與迭代法2種方法求解，不管哪種方法，計(jì)算效率都隨M的增加而顯著增加。利用直接法求解時(shí)需要顯示求取矩陣A的逆。雖然矩陣A是稀疏的，可壓縮存儲(chǔ)，但是矩陣A的逆并不一定是稀疏的，其大小與模型離散后總點(diǎn)數(shù)的平方成正比。在模型較大時(shí)，利用直接解法內(nèi)存存儲(chǔ)將是一個(gè)主要問(wèn)題。綜上所述，頻率域正演算法的計(jì)算效率與模型離散后的總點(diǎn)數(shù)M有直接關(guān)系。在差分算法一定的情況下，提高頻率域正演的計(jì)算效率、較低內(nèi)存開(kāi)銷(xiāo)的最直接途徑就是減少模型離散后的總點(diǎn)數(shù)。筆者借助MUMPS軟件包進(jìn)行稀疏線性方程組的求解[27]。

2 起伏地表?xiàng)l件下自由地表邊界條件及變網(wǎng)格技術(shù)的引入

對(duì)于聲波介質(zhì)，在自由地表需滿足應(yīng)力為0的邊界條件。在水平地表情況下，由于矩形網(wǎng)格與地表一致，較易處理自由邊界；當(dāng)?shù)乇砥鸱鼤r(shí)，由于矩形網(wǎng)格并不能完整地描述起伏地表，階梯狀的離散方式易在正演波場(chǎng)中產(chǎn)生角點(diǎn)散射。當(dāng)網(wǎng)格間距相對(duì)于地震波波長(zhǎng)足夠小時(shí)，角點(diǎn)散射可得到有效壓制。研究表明，當(dāng)一個(gè)波長(zhǎng)內(nèi)有15個(gè)以上采樣時(shí)，角點(diǎn)散射才能得到有效壓制[13]。本節(jié)首先概述起伏地表?xiàng)l件下頻率域正演自由地表邊界條件的加載方式，然后引入變網(wǎng)格技術(shù)，只對(duì)近地表附近網(wǎng)格進(jìn)行加密處理，提高算法的計(jì)算效率并降低內(nèi)存開(kāi)銷(xiāo)。

2.1 起伏地表?xiàng)l件下頻率-空間域自由地表邊界條件數(shù)值實(shí)現(xiàn)

雖然最優(yōu)9點(diǎn)有限差分算法的推導(dǎo)是以一階速度-應(yīng)力方程為基礎(chǔ)，并在數(shù)值離散的過(guò)程中用到了交錯(cuò)網(wǎng)格的相關(guān)算法，但是頻域正演的最終形式是同位網(wǎng)格的形式。相對(duì)于交錯(cuò)網(wǎng)格算法，該方法加入自由邊界較為容易；但相對(duì)于時(shí)域算法，在頻域正演中自由邊界條件的引入是隱式的，算法復(fù)雜性相對(duì)更高。

以圖2中黃色方框所示區(qū)域?yàn)槔M(jìn)行分析。圖2中綠色實(shí)線為地表實(shí)際位置，藍(lán)色階梯狀實(shí)線為利用矩形網(wǎng)格對(duì)地表進(jìn)行離散后的結(jié)果?？紤]粗網(wǎng)格情況，若求取A點(diǎn)的空間導(dǎo)數(shù)需要用到B-I周?chē)?個(gè)點(diǎn)的信息，由于E、F點(diǎn)位于自由表面以上，此處波場(chǎng)值為0。由于頻率域正演最終的形式是求取式(4)所示的線性方程組，波場(chǎng)向量P中均為未知量；因此，首先需要將E、F點(diǎn)從波場(chǎng)向量中去除。同時(shí)結(jié)合式(3)可知，E、F點(diǎn)波場(chǎng)值為0等價(jià)于在差分時(shí)不考慮該點(diǎn)的影響，也就是說(shuō)，在構(gòu)建波場(chǎng)向量P及阻抗矩陣A時(shí)不考慮自由表面以上的點(diǎn)。在實(shí)際處理過(guò)程中，需要重新定義笛卡爾坐標(biāo)系與線性坐標(biāo)系之間的映射關(guān)系，并需做好自由表面以上網(wǎng)格點(diǎn)的識(shí)別。

對(duì)比圖1b與圖2左側(cè)數(shù)字編號(hào)可知，在起伏地表?xiàng)l件下笛卡爾坐標(biāo)系與線性坐標(biāo)系的映射關(guān)系不能利用簡(jiǎn)單的函數(shù)關(guān)系表示，需要提前建立兩者之間完整的映射關(guān)系。自由表面以上網(wǎng)格點(diǎn)的識(shí)別在笛卡爾坐標(biāo)系中實(shí)現(xiàn)。

2.2 變網(wǎng)格技術(shù)及數(shù)值實(shí)現(xiàn)

利用矩形網(wǎng)格對(duì)起伏地表進(jìn)行階梯狀空間離散，在正演過(guò)程中會(huì)因?yàn)殡A梯狀離散產(chǎn)生嚴(yán)重的角點(diǎn)散射。這是矩形網(wǎng)格處理起伏地表的共同問(wèn)題，但角點(diǎn)散射會(huì)隨著網(wǎng)格間距的減小而減弱。如圖2所示，綠色曲線為實(shí)際地表高程，藍(lán)色階梯狀曲線為利用粗網(wǎng)格離散后的結(jié)果，紅色階梯狀曲線為利用細(xì)網(wǎng)格離散后的結(jié)果。通過(guò)對(duì)比可以看到，利用細(xì)網(wǎng)格離散后的結(jié)果不僅更接近真實(shí)地表形態(tài)，而且階梯狀幅度更小，將產(chǎn)生更少的虛假繞射波。

由前述可知，頻率域正演的計(jì)算效率、內(nèi)存開(kāi)銷(xiāo)與總共的網(wǎng)格點(diǎn)數(shù)有直接關(guān)系，全局細(xì)網(wǎng)格的正演算法在實(shí)際應(yīng)用中無(wú)法接受。引入變網(wǎng)格技術(shù)，只在近地表區(qū)域加密網(wǎng)格，可在較精確模擬近地表效應(yīng)的同時(shí)提高計(jì)算效率、降低內(nèi)存開(kāi)銷(xiāo)。

變網(wǎng)格離散策略如圖2所示。由于最優(yōu)9點(diǎn)格式有限差分方法計(jì)算時(shí)需要周?chē)粋€(gè)網(wǎng)格間距內(nèi)網(wǎng)格點(diǎn)上的波場(chǎng)值，細(xì)網(wǎng)格與粗網(wǎng)格的交界需要定義在地表最小高程一個(gè)粗網(wǎng)格間距以下的位置。在粗網(wǎng)格與細(xì)網(wǎng)格區(qū)域分別利用各自的網(wǎng)格間距進(jìn)行求解。在過(guò)渡區(qū)域，粗網(wǎng)格中空間求導(dǎo)要用到該網(wǎng)格點(diǎn)上方網(wǎng)格的波場(chǎng)值，該波場(chǎng)值可直接從細(xì)網(wǎng)格中對(duì)應(yīng)位置處取出；細(xì)網(wǎng)格中空間求導(dǎo)要用到該網(wǎng)格下方網(wǎng)格點(diǎn)的波場(chǎng)值，而該波場(chǎng)值大部分在粗網(wǎng)格中并沒(méi)有定義，需利用粗網(wǎng)格點(diǎn)上波場(chǎng)值插值得到，如圖2中空心圓圈所示。

圖2 起伏地表?xiàng)l件下基于變網(wǎng)格的網(wǎng)格離散策略Fig. 2 Strategy of grid discretization based on variable grid when consider irregular topography

細(xì)網(wǎng)格點(diǎn)分別位于左邊界(a)、右邊界(b)、粗網(wǎng)格內(nèi)部(c)、內(nèi)部粗網(wǎng)格上(d)。圖3 4種網(wǎng)格點(diǎn)內(nèi)插示意圖Fig. 3 Four different situations for grid interpolation

由于在粗網(wǎng)格上內(nèi)插得到的點(diǎn)并未在阻抗矩陣中顯式表示，而是隱式地表示為對(duì)應(yīng)粗網(wǎng)格點(diǎn)上的權(quán)系數(shù),所以需要考慮不同的情況分別進(jìn)行處理。對(duì)于如圖2所示的變網(wǎng)格離散策略，網(wǎng)格點(diǎn)內(nèi)插時(shí)共有如圖3所示的4種情況。當(dāng)細(xì)網(wǎng)格點(diǎn)在左右邊界及粗網(wǎng)格內(nèi)部時(shí)，共需要粗網(wǎng)格2個(gè)網(wǎng)格點(diǎn)的值(圖3a-c)；當(dāng)細(xì)網(wǎng)格在x方向與粗網(wǎng)格重合且不在左右邊界上時(shí)，共需要3個(gè)粗網(wǎng)格點(diǎn)上的值(圖3d)。具體實(shí)現(xiàn)方式以圖3d中情況為例進(jìn)行說(shuō)明。對(duì)于圖3d中標(biāo)號(hào)為1的網(wǎng)格點(diǎn)，構(gòu)建阻抗時(shí)需要其本身的值及周?chē)?個(gè)點(diǎn)的信息。由于標(biāo)號(hào)為7的網(wǎng)格點(diǎn)與粗網(wǎng)格點(diǎn)重合，所以只有標(biāo)號(hào)為6和8的網(wǎng)格點(diǎn)需要通過(guò)粗網(wǎng)格上的值插值得到。假設(shè)其插值用下式表示：

(5)

其中：w1、w2、w3、w4為插值權(quán)系數(shù)。將式(5)帶入式(3)可得到該情況下過(guò)渡區(qū)域的有限差分格式：

C1p1+C2p5+C3p9+C4p3+

(C5+R3w2+R4w3)p7+

(6)

筆者利用線性插值對(duì)缺失網(wǎng)格點(diǎn)進(jìn)行波場(chǎng)插值，該方法使得在網(wǎng)格過(guò)渡區(qū)域產(chǎn)生的虛假反射控制在真實(shí)波場(chǎng)的2%以內(nèi)，滿足波形反演的需求。通過(guò)應(yīng)用更適當(dāng)?shù)牟逯捣椒蛇M(jìn)一步提高正演的精度。此外，變網(wǎng)格技術(shù)引入以后，笛卡爾坐標(biāo)系與線性坐標(biāo)系之間的映射關(guān)系將更為復(fù)雜，進(jìn)一步增加了程序的復(fù)雜性。

a. 常規(guī)網(wǎng)格； b. 變網(wǎng)格。圖5 15 Hz頻域波場(chǎng)的實(shí)部Fig. 5 Real part of wavefield in frequency domain at 15 Hz

3 模型試算

3.1 均勻介質(zhì)模型

圖4 均勻速度模型Fig. 4 Homogeneous velocity model

如圖4所示，設(shè)計(jì)速度為2 500 m/s的均勻介質(zhì)模型，網(wǎng)格采樣間隔為10 m，縱橫向采樣點(diǎn)數(shù)均為150個(gè)。地表為一傾斜界面，高程由左到右逐漸減小，最大高差為149 m。震源位于模型中心位置，震源子波為15 Hz的雷克子波，檢波器位于地表以下30 m。正演共用80個(gè)頻率，頻率采樣率為0.5 Hz，除自由地表外其余3個(gè)邊界均應(yīng)用PML邊界條件，PML層厚度為15個(gè)網(wǎng)格點(diǎn)。

除在原始網(wǎng)格上進(jìn)行頻域正演外，還利用基于變網(wǎng)格技術(shù)的正演算法進(jìn)行正演模擬，如圖4所示。加密網(wǎng)格的厚度為500 m，細(xì)網(wǎng)格間距為2 m，相應(yīng)的加密倍數(shù)為5倍。正演所用其他參數(shù)與原始網(wǎng)格正演一致。

圖5a、b分別為應(yīng)用常規(guī)網(wǎng)格與變網(wǎng)格方法得到的15 Hz波場(chǎng)的實(shí)部。通過(guò)對(duì)比可知，2種網(wǎng)格正演波場(chǎng)基本一致，只是在自由表面處有明顯差異。由于頻率域波場(chǎng)對(duì)應(yīng)震源子波為正弦波的時(shí)域波場(chǎng)，在頻率域波場(chǎng)上較難進(jìn)行波場(chǎng)的對(duì)比，故在隨后的分析中將頻率域波場(chǎng)變換到時(shí)間域進(jìn)行分析。

0.316 s時(shí)刻的波場(chǎng)快照如圖6所示。對(duì)比圖6a與b可知，2種方法所得波場(chǎng)基本相同，細(xì)網(wǎng)格與粗網(wǎng)格邊界處產(chǎn)生的虛假反射在圖6b中完全看不到。將2種波場(chǎng)做差后的結(jié)果如圖6c所示。從振幅幅值對(duì)比上可知，虛假反射的振幅小于原始波場(chǎng)的2%，遠(yuǎn)小于原始波場(chǎng)的振幅。從波形上分析可知，圖6c中波場(chǎng)快照的差剖面主要由兩部分構(gòu)成，其中：標(biāo)號(hào)①所示部分，即能量較強(qiáng)的部分為細(xì)網(wǎng)格與粗網(wǎng)格交界處的虛假反射；標(biāo)號(hào)②所示部分，即能量較弱但分布范圍較廣的波場(chǎng)殘差是由自由邊界引起的，雖然直達(dá)波還未到達(dá)自由邊界，但是由于傅里葉變換的周期性，圖中波場(chǎng)殘差為較大時(shí)刻的波場(chǎng)殘差。通過(guò)兩者的定量對(duì)比分析可知，本文所述方法與其他變網(wǎng)格方法一樣，雖然在細(xì)網(wǎng)格與粗網(wǎng)格交界處產(chǎn)生虛假反射，但是相對(duì)于原始波場(chǎng)能量較弱，在波場(chǎng)的正反演中可以接受。

a. 常規(guī)網(wǎng)格正演結(jié)果；b. 變網(wǎng)格正演結(jié)果；c. 兩者之差。圖6 0.316 s時(shí)刻波場(chǎng)快照對(duì)比Fig. 6 Comparison of wavefield snaps at 0.316 s

0.48 s時(shí)刻波場(chǎng)快照：a. 常規(guī)網(wǎng)格；b. 變網(wǎng)格。0.62 s時(shí)刻波場(chǎng)快照：c. 常規(guī)網(wǎng)格；d. 變網(wǎng)格。圖7 不同時(shí)刻波場(chǎng)快照對(duì)比Fig. 7 Comparison of wavefield snaps at different time

圖7為2個(gè)較大時(shí)刻波場(chǎng)快照對(duì)比圖。由圖7a與b可知，引入變網(wǎng)格技術(shù)并沒(méi)有改變PML邊界吸收的效果，也沒(méi)有因?yàn)榫W(wǎng)格變化產(chǎn)生不穩(wěn)定(白色虛框所示區(qū)域)。同時(shí)，從圖7a與b中可以看到PML吸收邊界良好的吸收效果。對(duì)于起伏情況下自由邊界的處理效果，在圖7c與d中可以更明顯地看到。在圖7c中由于粗網(wǎng)格對(duì)起伏地表的階梯狀離散更為明顯，產(chǎn)生嚴(yán)重的角點(diǎn)散射(箭頭所指區(qū)域)。通過(guò)在近地表區(qū)域加密網(wǎng)格，有效地壓制了角點(diǎn)散射，如圖7d所示。在圖8所示的炮記錄上,也可以明顯地看到變網(wǎng)格技術(shù)的引入對(duì)角點(diǎn)散射的壓制。

a. 基于常規(guī)網(wǎng)格正演方法；b. 基于變網(wǎng)格正演方法。圖8 炮集對(duì)比Fig. 8 Comparison of shot gathers

3.2 復(fù)雜介質(zhì)模型

通過(guò)對(duì)復(fù)雜介質(zhì)模型進(jìn)行數(shù)值測(cè)試驗(yàn)證本文所述方法對(duì)模型復(fù)雜性的適應(yīng)性，同時(shí)定量分析變網(wǎng)格方法的引入對(duì)效率提高等方面的優(yōu)勢(shì)。

所用模型如圖9所示，橫向與縱向范圍均為2 400 m，速度為2 500～3 800 m/s。模型中既有大尺度的背斜、斷層等構(gòu)造，也有小尺度的繞射點(diǎn)，可測(cè)試算法對(duì)復(fù)雜模型的適應(yīng)性。起伏地表左側(cè)設(shè)計(jì)為連續(xù)變化的類(lèi)正弦狀起伏，右側(cè)為不連續(xù)變化的折線狀起伏，可測(cè)試算法對(duì)不同起伏情況的適應(yīng)性。繪圖時(shí)將起伏地表以上利用速度為1 500 m/s的均勻介質(zhì)填充以獲取較好的圖形顯示效果。為對(duì)比算法的計(jì)算效率與精度，常規(guī)網(wǎng)格選取間距為4 m的細(xì)網(wǎng)格，橫縱向均為600個(gè)采樣點(diǎn)。震源位于橫向1 400 m、地表以下40 m處，正演所用子波為主頻、為15 Hz的雷克子波，共用到110個(gè)頻率，頻率采樣間隔為0.416 67 Hz。PML吸收層厚度為20個(gè)采樣點(diǎn)。在基于變網(wǎng)格的正演方法中，粗網(wǎng)格處網(wǎng)格間距為12 m，表層加密后的網(wǎng)格間距為4 m，加密網(wǎng)格厚度為576 m，其他正演參數(shù)與細(xì)網(wǎng)格下正演參數(shù)一致。

圖9 復(fù)雜介質(zhì)模型Fig. 9 Complex velocity model

0.418 s時(shí)刻波場(chǎng)快照：a. 細(xì)網(wǎng)格；b. 變網(wǎng)格；c. 兩者之差。0.837 s時(shí)刻波場(chǎng)快照：d. 細(xì)網(wǎng)格；e. 變網(wǎng)格；f. 兩者之差。圖10 不同時(shí)刻波場(chǎng)快照對(duì)比Fig. 10 Comparison of wavefield snaps at different time

a. 常規(guī)方法在細(xì)網(wǎng)格下正演記錄；b. 本文所述方法所得正演記錄。圖11 炮集對(duì)比Fig. 11 Comparison of shot gather

利用2種方法所得不同時(shí)刻波場(chǎng)快照如圖10所示。圖10a與d分別為0.418 s與0.837 s時(shí)刻的波場(chǎng)快照，在波場(chǎng)快照上可以較清晰地看到起伏地表對(duì)波場(chǎng)的影響。標(biāo)號(hào)②處為直達(dá)波與虛反射的疊加；由于震源距離地表較近，兩者在旅行時(shí)上差異較少以致在波形上很難分開(kāi)。標(biāo)號(hào)①處為在右側(cè)高角度起伏地表產(chǎn)生的反射波；該類(lèi)型波使得隨后的波場(chǎng)傳播過(guò)程更為復(fù)雜，這也在一定程度上說(shuō)明了考慮起伏地表的重要性。標(biāo)號(hào)③處為地下第一個(gè)界面反射波在自由地表產(chǎn)生的表層多次波。通過(guò)波場(chǎng)分析可知，由于起伏地表的影響使得地下波場(chǎng)更為復(fù)雜，這種波場(chǎng)的復(fù)雜性在接收記錄上也可以清晰地看到，如圖11所示。由于地表的起伏主要反射軸已不完全是雙曲形態(tài)，給波場(chǎng)的分析帶來(lái)很大困難，這也在一定程度上驗(yàn)證了在反問(wèn)題中考慮起伏地表的重要性。

分別對(duì)比分析圖10a與b及圖10c與d中波場(chǎng)快照可知：基于變網(wǎng)格的正演方法與基于細(xì)網(wǎng)格的常規(guī)方法具有相同的精度，并沒(méi)有因?yàn)榫W(wǎng)格的變化而引起波場(chǎng)不穩(wěn)定等問(wèn)題；同時(shí)細(xì)網(wǎng)格與粗網(wǎng)格交界處的虛假反射由于能量相對(duì)很弱，在波場(chǎng)快照中基本看不到其存在。圖10c與f分別為2個(gè)時(shí)刻波場(chǎng)快照在2種方法下的波場(chǎng)差異。對(duì)比分析圖10a-c可以看到：由于地震波還沒(méi)有傳播到地下反射層，波場(chǎng)殘差主要是細(xì)網(wǎng)格與粗網(wǎng)格交界處產(chǎn)生的虛假反射(圖中白色箭頭所示)，虛假反射波振幅是原始波場(chǎng)的0.5%，對(duì)波場(chǎng)影響較小。與圖10c不同，圖10f由于是較大時(shí)刻的波場(chǎng)殘差，其中除了細(xì)網(wǎng)格與粗網(wǎng)格交界處的虛假反射外，更多的是由于576 m以下粗網(wǎng)格區(qū)域與細(xì)網(wǎng)格對(duì)介質(zhì)參數(shù)的離散采樣不同引起的，這是有限差分方法共同的問(wèn)題。然而通過(guò)能量對(duì)比可以看到，即使波場(chǎng)殘差較為復(fù)雜，但是整體能量相對(duì)原始波場(chǎng)較弱，在實(shí)際波場(chǎng)的正演與反演過(guò)程中是可以接受的。圖11為2種方法的觀測(cè)記錄。兩者波場(chǎng)基本一致，在圖11b上基本看不到由網(wǎng)格變化引起的虛假反射。

本模型試算所用程序利用Fortran語(yǔ)言編寫(xiě)，程序在Linux操作系統(tǒng)下利用ifort編譯，編譯時(shí)未加優(yōu)化選項(xiàng)，運(yùn)算CPU主頻為2.4 GHz的Intel雙核處理器(型號(hào)為E6600)，運(yùn)算時(shí)只用單線程?；诩?xì)網(wǎng)格的常規(guī)方法共用時(shí)15 min 31 s，內(nèi)存占用在500 MB左右；而基于變網(wǎng)格的方法用時(shí)2 min 36 s，內(nèi)存占用在90 MB左右?？梢钥吹?，通過(guò)引入變網(wǎng)格技術(shù)，計(jì)算效率提高了6倍以上，內(nèi)存占用降低到原先的1/5以內(nèi)。當(dāng)所用模型橫縱向范圍進(jìn)一步變大時(shí)或起伏地表高程差進(jìn)一步縮小時(shí)，在內(nèi)存占用及計(jì)算效率等方面的優(yōu)勢(shì)將更明顯。

4 結(jié)論

筆者建立并數(shù)值實(shí)現(xiàn)了頻率-空間域起伏地表?xiàng)l件下的數(shù)值正演算法，通過(guò)將時(shí)域正演中的變網(wǎng)格技術(shù)應(yīng)用到頻域正演，有效地提高了計(jì)算效率、減少了內(nèi)存開(kāi)銷(xiāo)。數(shù)值試驗(yàn)證明由網(wǎng)格間距變化引起的虛假反射波振幅幅值控制在原始波場(chǎng)的2%以內(nèi)，滿足波場(chǎng)正反演的要求；借助變網(wǎng)格技術(shù)可在保證精度的情況下較易得到5倍以上的加速比，且加速比隨著模型尺度的增加或起伏地表高程差的減小而增加。此外對(duì)內(nèi)存占用的大幅減少，使得對(duì)大尺度模型的頻域正演成為可能。

通過(guò)加入起伏地表自由邊界條件，進(jìn)一步擴(kuò)展了頻域正演方法的應(yīng)用范圍，可充分發(fā)揮頻率正演較易引入介質(zhì)黏滯性、對(duì)多源正演效率較高、PML類(lèi)吸收邊界精度更高等優(yōu)勢(shì)。為了提高計(jì)算效率，時(shí)域算法在空間變網(wǎng)格的同時(shí)往往還需要進(jìn)行局部時(shí)間變步長(zhǎng)，而頻域算法則不需要，相對(duì)而言減少了誤差來(lái)源，正演精度更高。此外，本文所述正演方法為起伏地表下的頻率域與Laplace-Fourier域波形反演奠定了堅(jiān)實(shí)的基礎(chǔ)。

本文在細(xì)網(wǎng)格與粗網(wǎng)格的過(guò)渡區(qū)域，通過(guò)線性插值對(duì)缺失網(wǎng)格點(diǎn)上的波場(chǎng)進(jìn)行插值。雖然得到了可接受的插值精度，但并不是最優(yōu)的方式；推導(dǎo)并應(yīng)用更精確的插值方法以取得更高精度的波場(chǎng)將是隨后的研究重點(diǎn)。

[1] Virieux J, Operto S. An Overview of Full-Waveform Inversion in Exploration Geophysics[J]. Geophysics, 2009,74(6):1-26.

[2] Pratt R G. Seismic Waveform Inversion in the Frequency Domain:Part 1: Theory and Verification in a Physical Scale Model[J]. Geophysics, 1999, 64(3): 888-901.

[3] Shin C, Cha Y H. Waveform Inversion in the Laplace-Fourier Domain[J]. Geophys J Int, 2009, 177(3): 1067-1079.

[4] Komatitsch D, Tromp J. Introduction to the Spectral-Element Method for 3D Seismic Wave Propagation[J]. Geophys J Int, 1999, 139(3):806-822.

[5] K?ser M, Dumbser M. An Arbitrary High-Order Discontinuous Galerkin Method for Elastic Waves on Unstructured Meshes: I： The Two-Dimensional Isotropic Case with External Source Terms[J]. Geophys J Int, 2006, 166(2): 855-877.

[6] Dumbser M, K?ser M, Puente J D L. Arbitrary High-Order Finite Volume Schemes for Seismic Wave Propagation on Unstructured Meshes in 2D and 3D[J]. Geophys J Int, 2007, 171(2)：665-694.

[7] Tessmer E, Kosloff D, Behle A. Elastic Wave Propagation Simulation in the Presence of Surface Topography[J]. Geophys J Int, 1992,108(2):621-632.

[8] Wang X C, Liu X W. 3-D Acoustic Wave Equation Forward Modeling with Topography[J]. Applied Geophysics, 2007, 4(1):8-15.

[9] 丘磊. 正交曲線坐標(biāo)系下的地震波數(shù)值模擬技術(shù)研究[D]. 杭州:浙江大學(xué)，2011. Qiu Lei. Study on Seismic Wave Simulations in Orthogonally Curvilinear Coordinate System[D]. Hangzhou: Zhejiang University, 2011.

[10] 張劍鋒. 彈性波數(shù)值模擬的非規(guī)則網(wǎng)格差分法[J]. 地球物理學(xué)報(bào)，1998,41(增刊1):357-366. Zhang Jianfeng. Non-Orthogonal Grid Finite-Diffe-rence Method for Numerical Simulation of Elastic Wave Propagation[J]. Acta Geophysica Sinica, 1998, 41(Sup.1): 357-366.

[11] K?ser M, Igel H. Numerical Simulation of 2D Wave Propagation on Unstructured Grids Using Explicit Differential Operators[J]. Geophysical Prospecting, 2001, 49(5) :607-619.

[12] Jih R, McLaughlin K L, Der Z A. Free-Boundary Conditions of Arbitrary Polygonal Topography in a Two-Dimensional Explicit Elastic Finite-Difference Scheme[J]. Geophysics, 1988, 53(8):1045-1055.

[13] Robertsson J O A, Holliger K. Modeling of Seismic Wave Propagation near the Earth’s Surface[J]. Physics of the Earth and Planetary Interiors, 1997, 104(1/2/3): 193-211.

[14] Hayashi K, Bruns D R, Toks?z M N. Discontinuous-Grid Finite-Difference Seismic Modeling Including Surface Topography[J]. Bulletin of the Seismological Society of America, 2001, 91(6): 1750-1764.

[15] Wang Y, Xu J, Schuster G. Viscoelastic Wave Simulation in Basin by a Variable-Grid Finite-Difference Method[J]. Bulletin of the Seismological Society of America, 2001, 91(6): 1741-1749.

[16] Jo C, Shin C, Suh J H. An Optimal 9-Point, Finite-Difference,Frequency-Space 2-D Scalar Wave Extra-polator[J]. Geophysics, 1996, 61(2):529-537.

[18] Shin C, Sohn H. A Frequency-Space 2-D Scalar Wave Extrapolator Using Extended 25-Point Finite-Difference Operator[J]. Geophysics, 1998, 63(1): 289-296.

[19] Hustedt B, Operto S, Virieux J. Mixed-Grid and Staggered-Grid Finite-Difference Methods for Frequency-Domain Acoustic Wave Modeling[J]. Geophys J Int, 2004, 157(3): 1269-1296.

[20] 曹書(shū)紅，陳景波. 聲波方程頻率域高精度正演的17點(diǎn)格式及數(shù)值實(shí)現(xiàn)[J]. 地球物理學(xué)報(bào)， 2012，55(10)：3440-3449. Cao Shuhong, Chen Jingbo. A 17-point Scheme and Its Numerical Implementation for High-Accuracy Modeling of Frequency-Domain Acoustic Equation[J]. Chinese Journal Geophysics, 2012,55(10): 3440-3449.

[21] 吳國(guó)忱，梁鍇. VTI介質(zhì)qP波方程高精度有限差分算子[J]. 地球物理學(xué)進(jìn)展，2007, 22(3)：896-904. Wu Guochen, Liang Kai. High Precision Finite Difference Operators for qP Wave Equation in VTI Media[J]. Progress in Geophysics, 2007, 22(3): 896-904.

[22] 梁鍇，吳國(guó)忱，印興耀. TTI介質(zhì)qP波方程頻率-空間域加權(quán)平均有限差分算子[J]. 石油地球物理勘探，2007,42(5)：516-525. Liang Kai, Wu Guochen, Yin Xingyao. Weighted Mean Finite-Difference Operator of qP Wave Equation in Frequency-Space Domain for TTI Medium[J]. Oil Geophysical Prospecting, 2007, 42(5):516-525.

[23] 韓立，韓立國(guó)，李翔，等. 二階聲波方程頻域PML邊界條件及頻域變網(wǎng)格步長(zhǎng)并行計(jì)算[J]. 吉林大學(xué)學(xué)報(bào)：地球科學(xué)版，2011，41(4)：1226-1232. Han Li, Han Liguo, Li Xiang, et al. PML Boundary Conditions for Second-Order Acoustic Wave Equations and Variable Grid Parallel Computation in Frequency-Domain Modeling[J]. Journal of Jilin University: Earth Science Edition, 2011, 41(4):1226-1232.

[24] Brossier R, Virieux J, Operto S. Parsimonious Finite-Volume Frequency-Domain Method for 2D P-SV Wave Modeling[J]. Geophy J Int, 2009, 175(2): 541-559.

[25] Pan G, Abubakar A, Habashy T M. An Effective Perfectly Matched Layer Design for Acoustic Fourth-Order Frequency-Domain Finite-Difference Scheme[J]. Geophys J Int, 2012, 188(1): 211-222.

[26] 田坤，黃建平，李振春, 等. 實(shí)現(xiàn)復(fù)頻移完全匹配層吸收邊界條件的遞推卷積方法[J]. 吉林大學(xué)學(xué)報(bào)：地球科學(xué)版，2013，43(3)：1022-1032. Tian Kun, Huang Jianping, Li Zhenchun, et al. Recursive Convolution Method for Implementing Complex Frequency-Shifted PML Absorbing Boundary Condition[J]. Journal of Jilin University: Earth Science Edition, 2013, 43(3):1022-1032.

[27] Amestoy P R, Duff I S, Koster J, et al. A Fully Asynchronous Multifrontal Solver Using Distributed Dynamic Scheduling[J]. SIAM Journal of Matrix Analysis and Applications, 2001, 23(1):15-41.

Acoustic Wave Modeling in Frequency-Spatial Domain with Surface Topography

Guo Zhenbo, Li Zhenchun

SchoolofGeosciences,ChinaUniversityofPetroleum(EastChina),Qingdao266580,Shandong,China

Seismic wave modeling in frequency-space domain is the foundation of full waveform inversion in frequency and Laplace-Fourier domain. An accurate seismic wave modeling with surface topography is the kernel of full waveform inversion which can direct process data recorded on relief surface. Based on acoustic wave equation with PML absorption boundary condition, the free surface condition on relief surface is included into the optimized 9-points finite difference algorithm in frequency space domain. Then the variable grid method which refines the grid near the surface and remains the rest unchangeable is introduced to increase the computational efficiency and decrease the memory cost. The theoretical analysis and numerical test show that: 1) the diffraction caused by grid discretization of the relief surface can be suppressed by refining the grid; 2) through the variable grid method, computational efficiency can be improved for several times and the memory cost can also be decreased for several times easily; 3) the energy of artificial reflection caused by the boundary of fine and coarse grid is less than 2% of the energy of the original wavefield which is acceptable for the seismic wavefield modeling and inversion.

surface topography; variable grid; frequency domain; numerical modeling

10.13278/j.cnki.jjuese.201402305.

2013-07-04

國(guó)家自然科學(xué)基金項(xiàng)目(40974073)；國(guó)家“863”計(jì)劃項(xiàng)目(2011AA060301)；中央高?；究蒲袠I(yè)務(wù)費(fèi)專(zhuān)項(xiàng)資金(13CX06014A)

郭振波(1986-)，男，博士研究生，主要從事地震波場(chǎng)正演模擬與波形反演等方面的研究，E-mail:stone_gzb@163.com。

10.13278/j.cnki.jjuese.201402305

P631.4

A

郭振波，李振春. 起伏地表?xiàng)l件下頻率-空間域聲波介質(zhì)正演模擬.吉林大學(xué)學(xué)報(bào):地球科學(xué)版,2014,44(2):683-693.

Guo Zhenbo, Li Zhenchun. Acoustic Wave Modeling in Frequency-Spatial Domain with Surface Topography.Journal of Jilin University:Earth Science Edition,2014,44(2):683-693.doi:10.13278/j.cnki.jjuese.201402305.