非均勻黏彈性介質(zhì)中聲波模擬的非規(guī)則網(wǎng)格方法

豆輝，張劍鋒

1 中國(guó)科學(xué)院地質(zhì)與地球物理研究所，中國(guó)科學(xué)院油氣資源研究重點(diǎn)實(shí)驗(yàn)室，北京 1000292 中國(guó)科學(xué)院大學(xué)，北京 100049

?

非均勻黏彈性介質(zhì)中聲波模擬的非規(guī)則網(wǎng)格方法

豆輝1,2，張劍鋒1

1 中國(guó)科學(xué)院地質(zhì)與地球物理研究所，中國(guó)科學(xué)院油氣資源研究重點(diǎn)實(shí)驗(yàn)室，北京 1000292 中國(guó)科學(xué)院大學(xué)，北京 100049

本文研究了滿足衡Q模型的黏彈性介質(zhì)中聲波模擬的時(shí)間域數(shù)值解法.通過采用有理式基函數(shù)描述頻率相關(guān)的體積模量，使模型滿足地震波黏性吸收的衡Q要求.結(jié)合彈性波模擬的非規(guī)則、非結(jié)構(gòu)網(wǎng)格方法——格子法，本文提出了一種非均勻黏彈性介質(zhì)中聲波模擬的非規(guī)則網(wǎng)格方法.非規(guī)則、非結(jié)構(gòu)網(wǎng)格的使用，可以精細(xì)地刻畫地下介質(zhì)中復(fù)雜速度、Q值的變化，及界面的形狀.通過和解析解的比較及復(fù)雜模型算例分析，驗(yàn)證了該方法的精度及有效性.

時(shí)間域；地震波模擬；非規(guī)則網(wǎng)格；黏彈性；衡Q模型；有理函數(shù)

1 引言

地震波數(shù)值模擬對(duì)于定量地分析研究復(fù)雜地下介質(zhì)中的地震波成像非常重要.傳統(tǒng)的地震波模擬基于彈性模型假設(shè)，可以滿足一般的地震勘探精度要求，得到很好的數(shù)值模擬結(jié)果(Gazdag，1981).但近年來勘探技術(shù)的發(fā)展對(duì)數(shù)值模擬的精度提出了更高的要求，需要考慮地下介質(zhì)黏性對(duì)數(shù)值模擬結(jié)果的影響.同時(shí)，計(jì)算機(jī)技術(shù)的飛速發(fā)展為高精度黏彈性數(shù)值計(jì)算提供了可能.前人關(guān)于黏彈性數(shù)值模擬的研究分為頻率域和時(shí)間域兩類.頻率域數(shù)值模擬可以適用任意的衰減準(zhǔn)則描述介質(zhì)的衰減特性，得到廣泛的應(yīng)用(Brossier et al.，2010；廖建平等，2011).但該方法在處理大型數(shù)據(jù)時(shí)可能面臨著內(nèi)存需求過大的問題.而時(shí)間域數(shù)值模擬時(shí)，面臨著黏性應(yīng)力-應(yīng)變之間的褶積運(yùn)算會(huì)增加數(shù)值計(jì)算難度的問題.為改善時(shí)間域數(shù)值計(jì)算的困難，前人們提出了很多描述黏性衰減的策略，如引入流變模型或改變波動(dòng)方程表示等.

目前描述衰減問題常用的流變模型，是通過不同形式的串并聯(lián)線性彈簧體和阻尼器得到的一系列簡(jiǎn)單的線性黏彈性固體模型，如Kelvin-Voigt固體模型，Maxwell固體模型等(Blake et al.，1982).這些流變模型基于相應(yīng)的物理模型，描述介質(zhì)的蠕變和松弛響應(yīng)，及介質(zhì)的應(yīng)力應(yīng)變關(guān)系.但也存在一定的局限性.首先，和實(shí)驗(yàn)所測(cè)材料的應(yīng)力應(yīng)變關(guān)系相比，Maxwell模型和Kelvin-Voigt模型都只能部分(蠕變或松弛過程)和實(shí)驗(yàn)結(jié)果相吻合.其次，該類方法沒有考慮模型的Q(ω)函數(shù)是否符合衡Q要求，不能得到地震波場(chǎng)的真實(shí)振幅，與實(shí)際地球介質(zhì)屬性不相符.

另一種描述衰減的策略是對(duì)彈性波方程做變動(dòng)(吳玉等，2014)，以克服傳統(tǒng)計(jì)算中的內(nèi)存需求過大的問題.Zhu和Harris(2013，2014)、Zhu(2014)基于Kjartansson(1979)的衡Q模型和Chen和Holm(2004)的分?jǐn)?shù)式Laplace算子方法，通過在彈性波波動(dòng)方程上加上黏性近似項(xiàng)，實(shí)現(xiàn)黏彈性介質(zhì)的時(shí)間域數(shù)值模擬.該方法可以描述近乎衡Q的衰減頻散響應(yīng)，但所加的黏性項(xiàng)是近似項(xiàng)，容易造成算法的不穩(wěn)定，導(dǎo)致計(jì)算精度受限.

還有一種方法是Liu等(1976)提出的用松弛機(jī)制譜的概念定義本構(gòu)關(guān)系.Liu認(rèn)為黏性Q在地震波頻帶范圍內(nèi)近似為常數(shù)，可由多個(gè)松弛機(jī)制的組合近似表示.Day和Minster(1984)更進(jìn)一步提出，借助Padé近似將時(shí)間域的褶積算子化為一系列微分算子，進(jìn)而由一個(gè)n階有理函數(shù)近似表示黏彈性模量，并解析地求出其中的參數(shù).Emmerich和Korn(1987)在此基礎(chǔ)上，賦予了該黏彈性模量相應(yīng)的物理意義——可由廣義的Maxwell流變模型表示.該廣義模型由一個(gè)彈簧體和n個(gè)經(jīng)典的Maxwell體并聯(lián)組成，n個(gè)經(jīng)典的Maxwell體分別對(duì)應(yīng)有理式中的n個(gè)疊加項(xiàng).該類方法可以近似地滿足衡Q模型要求，但需要尋得盡可能低階的n，以減少計(jì)算時(shí)對(duì)內(nèi)部變量的存儲(chǔ)需求，保證計(jì)算效率.

本文描述衰減問題的方法就是基于上述策略展開的.為了快速地搜索到低階的又能滿足衡Q要求的模型參數(shù)，本文以無窮范數(shù)構(gòu)建目標(biāo)函數(shù)，結(jié)合全局搜索的模擬退火方法尋求最優(yōu)解.最終確定，當(dāng)階數(shù)n=2或3時(shí)，我們就能找到滿足要求的衡Q模型，從而極大地減少了數(shù)值計(jì)算量和所需的存儲(chǔ)空間.另一方面，為進(jìn)一步提高計(jì)算效率和精度，本文結(jié)合了非規(guī)則非結(jié)構(gòu)網(wǎng)格的格子法(Zhang and Liu，1999)作為黏彈性介質(zhì)的數(shù)值模擬算法.相比于常見的有限差分法(Robertsson et al.，1994；Saenger and Bohlen，2004；Operto et al.，2007；王忠仁等，2014)、有限元法(Marfurt，1984)和偽譜法(Carcione，2010)等，格子法可以自然地滿足自由邊界條件，更精細(xì)地刻畫復(fù)雜介質(zhì)的速度變化及Q值變化.為驗(yàn)證本文算法的精度和有效性，文中給出了簡(jiǎn)單模型算例中數(shù)值解和解析解的比較，測(cè)試了算法在氣云構(gòu)造和Marmousi模型中的模擬效果，都取得了很好的結(jié)果.

2 非規(guī)則、非結(jié)構(gòu)網(wǎng)格聲波模擬的格子法

本文數(shù)值解法是基于聲波模擬的格子法(Zhang and Liu，1999，2002)展開工作.格子法采用非規(guī)則、非結(jié)構(gòu)化網(wǎng)格離散速度模型，既可以精細(xì)地刻畫復(fù)雜邊界條件、速度及Q值的變化，也具有與有限差分法相當(dāng)?shù)挠?jì)算效率.

這里簡(jiǎn)單地介紹一下聲波模擬下格子法的基本原理.考慮非均勻介質(zhì)中的無源聲波方程(徐義和張劍鋒，2008)：

(1)

格子法的理論核心是在非規(guī)則網(wǎng)格下，給出式(1)的基于積分平衡的微分方程弱形式.如圖1所示，在節(jié)點(diǎn)k的鄰域，定義虛線輪廓線1-2-3-…-12-1所圍區(qū)域?yàn)棣?，?jiǎn)記該輪廓線為?Ω，對(duì)式(1)兩邊作面積分，并對(duì)等式右邊應(yīng)用Green定理得到

(2)

式中m為圍繞節(jié)點(diǎn)k的三角形單元數(shù)目，Skl對(duì)應(yīng)節(jié)點(diǎn)k周圍第l個(gè)三角形單元中的虛線段，α與β分別是邊界Ω外法線沿x與z軸的方向余弦.利用動(dòng)力學(xué)計(jì)算的集中質(zhì)量原理及三角形線性插值方法可得到式(2)的離散形式：

(3)

式中Mk是節(jié)點(diǎn)k各鄰域三角形單元中的?Ω1/ρv2dΩ值總和的三分之一，(ρ)l表示第l個(gè)單元的密度，Dxp,Dzp是三角形差分算子，可表述為

(4)

圖1 非規(guī)則三角網(wǎng)格局部示意圖Fig.1 Local irregular triangle grid

式中，Δ為三角形ijk的面積，下標(biāo)r代表三角形中的各頂點(diǎn)，(bk)l，(ck)l分別表示圍繞節(jié)點(diǎn)k的第l個(gè)單元對(duì)應(yīng)的幾何參數(shù)，以圖1中的三角形單元ijk為例，可簡(jiǎn)單地表述為bk=(zj-zi)/2,ck=(xj-xi)/2.

最后，利用中心差分格式離散方程式(3)左端中對(duì)時(shí)間的二階偏導(dǎo)數(shù)，我們就可以實(shí)現(xiàn)聲波波場(chǎng)的迭代更新.

3 體積模量的有理函數(shù)展開和衡Q模型

3.1 體積模量的有理函數(shù)展開

黏彈性介質(zhì)在頻率域中有如下應(yīng)力應(yīng)變關(guān)系：

(5)

其中，M(ω)是與頻率有關(guān)的復(fù)體積模量，σ(ω)是應(yīng)力，ε(ω)為應(yīng)變，經(jīng)由Fourier變換，得到如下褶積關(guān)系：

(6)

這意味著在時(shí)間域計(jì)算時(shí)，需要存儲(chǔ)歷史時(shí)刻的應(yīng)變值，會(huì)導(dǎo)致極大的計(jì)算負(fù)擔(dān)，是數(shù)值計(jì)算很不樂見的.所以，我們需要考慮從其他方面對(duì)其做簡(jiǎn)化近似處理.

已知實(shí)驗(yàn)測(cè)得實(shí)際介質(zhì)的應(yīng)力松弛響應(yīng)關(guān)系為：在常應(yīng)變作用下，應(yīng)力呈指數(shù)衰減變化，即實(shí)際黏性介質(zhì)的松弛函數(shù)可用一個(gè)負(fù)指數(shù)衰減的函數(shù)表示.不同黏性介質(zhì)的松弛響應(yīng)對(duì)應(yīng)不同的幅值及衰減系數(shù)表示.我們總可以找到這樣的一組基函數(shù)——由一系列幅值和衰減系數(shù)組合表示任意的衰減函數(shù)，并用其表示介質(zhì)的松弛響應(yīng).我們定義未松弛模量MU及松弛模量MR，分別表示介質(zhì)衰減前、后的體積模量，并基于此構(gòu)建如下負(fù)指數(shù)衰減的松弛函數(shù)：

(7)

至此，可將體積模量表示為

(8)

我們通過Fourier變換，將其變換到頻率域(e-ωjt→FFT→1/(iω+ωj))，有

(9)

黏彈性介質(zhì)中，頻率描述的體積模量可以由n個(gè)松弛頻率的有理式函數(shù)展開表示.

3.2 衡Q模型

實(shí)際地球介質(zhì)中，黏性Q(ω)是隨頻率變化的函數(shù).但是在地震頻帶范圍內(nèi)Q可近似為常數(shù)，或隨頻率微弱變化.地球物理勘探中的信號(hào)頻帶一般在3～100 Hz內(nèi)，我們可以采用在該頻帶范圍內(nèi)滿足衡Q要求的模型描述黏彈性介質(zhì)的高精度勘探問題.在介質(zhì)黏性吸收滿足衡Q要求時(shí)，我們就可以確定體積模量Mn(ω)中有理式函數(shù)的系數(shù).

針對(duì)式(9)中有理函數(shù)表示的體積模量Mn(ω)，品質(zhì)因子Q可表示為

(10)

為尋求滿足衡Q模型要求的最優(yōu)解，我們用Q0表示目標(biāo)Q值，以無窮范數(shù)構(gòu)建目標(biāo)函數(shù)，得到下列表示：

(11)

圖2 階數(shù)(a)n=2，(b)n=3時(shí)，Q-1=0.05的目標(biāo)值(實(shí)線)和Q-1(ω)曲線擬合結(jié)果(虛線)Fig.2 Target values (solid lines) and curve fitting results (dash lines) of the Q-1(ω),with different order:(a) n=2,(b) n=3

一般低階(n≤8)時(shí)有2n+1個(gè)模型參數(shù)，還比較少，我們可以采用蒙特卡洛、模擬退火、遺傳算法等全局搜索方法.本文采用的是優(yōu)化的模擬退火方法(Goffeetal.，1994)，可實(shí)現(xiàn)全局的快速搜索，求得最優(yōu)解.然后將最優(yōu)解對(duì)應(yīng)的參數(shù)結(jié)果列表記錄下來，便于以后模擬計(jì)算時(shí)查表使用.最終得到的結(jié)果和目標(biāo)值的擬合效果，可以參見圖2所示，圖中實(shí)線表示目標(biāo)值Q0=20.0，虛線為Q-1(ω)擬合結(jié)果.圖2中，分別對(duì)應(yīng)階數(shù)n=2，n=3時(shí)Q(ω)擬合得到的曲線結(jié)果和目標(biāo)值的比較.可以看出，在頻帶3～100 Hz內(nèi)，n=2時(shí)即可得到不錯(cuò)的擬合效果；n=3時(shí)，Q(ω)曲線可以更加貼合目標(biāo)值，擬合效果更優(yōu).本文數(shù)值模擬的簡(jiǎn)單算例，采用的就是階數(shù)n=3時(shí)的參數(shù)結(jié)果.

假設(shè)δM/MU?1，則式(10)可近似為

(12)

4 非均勻黏彈性介質(zhì)中聲波模擬的非規(guī)則網(wǎng)格方法

將式(9)表示的體積模量Mn(ω)代入式(5)中，并引入中間變量φj(ω)，得到

(13)

式中，中間變量φj(ω)滿足關(guān)系式

(14)

對(duì)式(13)和式(14)做Fourier反變換，得到時(shí)間域黏彈性介質(zhì)的波動(dòng)方程：

(15)

由此，我們得到如下非均勻黏彈性介質(zhì)中的聲波方程組：

(16)

(17)

(18)

式中的系數(shù)aj,ωj,及δM可由本文3.2節(jié)中的計(jì)算列表查詢得到.當(dāng)δM=0時(shí)，介質(zhì)不含黏性衰減，則式(16)—(18)可合并成式(1)所示的聲波方程.

為便于模擬復(fù)雜介質(zhì)下波傳播情況，本文采用格子法中的三角網(wǎng)格單元，即非規(guī)則、非結(jié)構(gòu)網(wǎng)格離散速度模型，并延用格子法原理到黏聲波方程組(式(16)—(18))，最終得到非規(guī)則網(wǎng)格的黏聲波數(shù)值模擬解法.

仍以圖1為例，延用聲波格子法的積分近似推導(dǎo)，在虛線包圍的k節(jié)點(diǎn)鄰域?qū)κ?17)做面積分，并利用Green定理，得到

(20)

式中，ck,bk和式(4)中三角單元的幾何參數(shù)有相同的意義，Δ為三角單元的面積.最終，對(duì)式(16)、(20)和(18)使用時(shí)間域的中心差分近似，可得

(21)

(22)

(23)

式中上標(biāo)q代表時(shí)間，下標(biāo)l表示節(jié)點(diǎn)k周圍第l個(gè)三角形單元.式(21)-(23)構(gòu)成了非規(guī)則網(wǎng)格黏聲波方法的應(yīng)變和壓強(qiáng)更新基礎(chǔ).由q時(shí)刻的壓強(qiáng)p代入式(22)可求得q+0.5時(shí)刻各節(jié)點(diǎn)處的形變?chǔ)?，再結(jié)合式(23)求得q+0.5時(shí)刻各節(jié)點(diǎn)處的中間變量φj，最后由式(21)求得q+1時(shí)刻的壓強(qiáng)波場(chǎng)p,完成對(duì)壓強(qiáng)p在時(shí)間上的更新.

相比于聲波方程，該黏聲波方程需要額外計(jì)算n個(gè)一階時(shí)間導(dǎo)方程，儲(chǔ)存n+1個(gè)中間變量ε(x,z,t),φj(x,z,t)(j=1,…,n).為保證計(jì)算效率，我們需要使用盡可能低階的階數(shù)n，同時(shí)又要滿足地震波黏性吸收的衡Q要求.本文3.2節(jié)已經(jīng)確定在n=2或3時(shí)找到很好的滿足衡Q要求的參數(shù).另一方面，格子法中的非規(guī)則網(wǎng)格的使用，允許在數(shù)值計(jì)算時(shí)，黏彈性介質(zhì)區(qū)域參與黏聲波方程的計(jì)算，彈性介質(zhì)仍按聲波方程模擬計(jì)算，也可大大縮減計(jì)算存儲(chǔ)需求，提高計(jì)算效率.

數(shù)值模擬時(shí)四周的邊界條件，延用徐義和張劍鋒(2008)的非規(guī)則PML吸收方法.在吸收帶內(nèi)，本文采用隨邊界位置變化的局部坐標(biāo)系，見圖1中的(w,v)坐標(biāo)系，結(jié)合引入的中間變量φj，構(gòu)建基于局部坐標(biāo)系的黏聲波PML分裂方程，這里不再詳述.該方法可以很好地吸收或透射不同方向的入射波，避免邊界反射問題.

5 數(shù)值算例

為驗(yàn)證本文提出的非規(guī)則網(wǎng)格黏聲波數(shù)值模擬解法的正確性和精度，首先將其應(yīng)用于均勻速度模型(對(duì)此模型可容易得到黏聲波傳播的解析解)，并將其計(jì)算結(jié)果與解析解對(duì)比分析.然后，給出一個(gè)強(qiáng)Q衰減擾動(dòng)情況下的數(shù)值模擬解算例，以驗(yàn)證該計(jì)算方法的穩(wěn)定性.最后應(yīng)用到Marmousi模型考察該方法對(duì)復(fù)雜橫向速度、Q值變化情況的適應(yīng)能力.

5.1 均勻介質(zhì)模型

為了評(píng)估本文解法的精度，我們分別從該算法對(duì)炮檢距和Q值的響應(yīng)這兩個(gè)方面討論、比較數(shù)值模擬得到的數(shù)值解與解析解.用于對(duì)比的解析解，是根據(jù)對(duì)應(yīng)原理(Carcione et al.，1988)，在頻率域中用其中的黏彈性參數(shù)代替彈性解析解中的彈性參數(shù)，再反變換到時(shí)間域得到.數(shù)值解方面，我們?cè)O(shè)計(jì)了一個(gè)橫向4 km，縱向2 km的各向同性均勻黏彈性介質(zhì)模型，密度、速度分別為ρ=2.0 g·cm-3,v=2000 m·s-1，炮點(diǎn)位于模型的中心位置.文中我們采用主頻為30Hz的Ricker子波作為震源，時(shí)間采樣步長(zhǎng)為0.33 ms，利用非規(guī)則網(wǎng)格離散該速度模型共得到4527531個(gè)網(wǎng)格單元參與數(shù)值計(jì)算.

圖3 不同炮檢距時(shí)，數(shù)值解(虛線)和解析解(實(shí)線)的比較.檢波器距離震源距離依次為(a) 600，(b) 1600，(c) 3200 m.二者結(jié)果吻合較好Fig.3 Comparison of numerical solutions (dash line) with the analytical solutions (solid line).The receivers are at (a) 600,(b) 1600,(c) 3200 m from the source.The two results are matched with each other well

5.1.1 不同炮檢距時(shí)

本文首先研究Q=20時(shí)，數(shù)值解對(duì)炮檢距的敏感度——能否在遠(yuǎn)炮檢距處保證相應(yīng)的計(jì)算精度.我們分別在距震源600 m，1600 m，3200 m處，放置一個(gè)檢波器記錄波場(chǎng)值，并與相同位置的解析解對(duì)比，最終結(jié)果如圖3所示(圖中實(shí)線表示解析解，虛線表示數(shù)值解).從圖3中可以看出，該算法的數(shù)值解與解析解都能較好的吻合.為進(jìn)一步定量地評(píng)估本文算法的精度，本文利用數(shù)值解與解析解的均方根誤差來衡量，并定義均方根誤差E：

(24)

5.1.2 不同Q值時(shí)

在這部分，我們討論數(shù)值模擬結(jié)果對(duì)Q值的敏感度，尤其在低Q時(shí)，模擬結(jié)果的精度問題.仍采用上面的均勻速度模型和震源，在點(diǎn)(2.6 km，1.0 km)處放置一檢波器，記錄四個(gè)不同Q值時(shí)波傳播情況.圖4為不同Q值時(shí)，非規(guī)則網(wǎng)格數(shù)值模擬結(jié)果在T=200 ms時(shí)刻的波場(chǎng)快照?qǐng)D：從左到右，從上到下，Q值分別為1000，100，30，10.從圖4中可以看到，Q值的減小伴隨著波場(chǎng)幅值更強(qiáng)的衰減及相位延遲.為了更詳細(xì)地了解模擬結(jié)果對(duì)Q的敏感度，我們?nèi)〔煌琎值在同一檢波器處的波場(chǎng)記錄，并與相應(yīng)的解析解對(duì)比分析，如圖5所示，圖中實(shí)線表示解析解，虛線為數(shù)值解.我們?nèi)岳檬?24)得到相應(yīng)的均方根誤差分別為3.8×10-3，3.9×10-3，5.4×10-3，1.3×10-2，說明Q值很小時(shí)，本文解法仍可得到穩(wěn)定的高精度的計(jì)算結(jié)果.

圖4 均勻衰減介質(zhì)中，T=200 ms時(shí)刻的波場(chǎng)快照?qǐng)D，圖中的四個(gè)部分分別對(duì)應(yīng)不同的Q值：(a)1000，(b)100，(c)30，(d)10.震源位于模型的中間位置.隨著Q變小，地震波波場(chǎng)有更強(qiáng)的幅值衰減和相位延遲Fig.4 Four snapshot parts corresponding to four Q values:(a) 1000,(b) 100,(c) 30,(d) 10.A point source located at the center of the model.Snapshots are recorded at 200 ms in homogeneous attenuating media.The amplitude of wavefields decreases rapidly with the Q values decline

圖5 震源位于點(diǎn)(2.0,1.0)km，檢波器位于(2.6,1.0)km處，不同的Q值：(a)1000，(b)100，(c)30，(d)10時(shí)，對(duì)應(yīng)的數(shù)值解(虛線)和解析解(實(shí)線)的比較.它們的結(jié)果都能很好的吻合Fig.5 Comparison between simulated solutions (dash lines) and analytical solutions (solid lines) corresponding to four Q values:(a) 1000,(b) 100,(c) 30,and (d) 10 at point (2.6,1.0) km from a source (2.0,1.0) km .The results show an excellent agreement

以上兩個(gè)數(shù)值算例結(jié)果都驗(yàn)證了本文數(shù)值解法在遠(yuǎn)炮檢距、低Q處均具有穩(wěn)定的高精度解.

圖6 (a)370 ms,(b)450 ms時(shí)刻的波場(chǎng)快照.在圖中白色圓圈代表的氣云周圍有反射產(chǎn)生Fig.6 Snapshots of pressures at (a)370 ms and(b)450 ms propagation time.There are reflections at interface of the gas-pocket which is represented by the white circle in the figures

5.2 氣云模型

具有強(qiáng)衰減(即低Q)特征的氣云構(gòu)造常給地震成像帶來很大的困擾，應(yīng)用本文數(shù)值算法到氣云模型，可以驗(yàn)證算法的穩(wěn)定性.模型設(shè)計(jì)如下：模型大小為2000 m×1000 m，地下介質(zhì)為速度v=2000 m·s-1的黏彈性介質(zhì).在點(diǎn)(600 m，600 m)處有一半徑為50 m的圓形氣云(圖6中的白色圓圈)，氣云內(nèi)Q=10，氣云外為Q=∞(即彈性介質(zhì)).數(shù)值模擬時(shí)，取主頻為30 Hz的Ricker子波作為爆炸源置于點(diǎn)(1000 m，10 m)處(圖6中的星號(hào))，采樣時(shí)間步長(zhǎng)為0.333 ms.圖6a和6b分別對(duì)應(yīng)傳播時(shí)間為370 ms和450 ms時(shí)的地震波場(chǎng)快照.從圖中可以看出，在氣云界面有弱的反射波產(chǎn)生，透過氣云的波場(chǎng)存在幅值的衰減和相位的略微滯后現(xiàn)象.說明本文方法在強(qiáng)Q變化模型中可得到穩(wěn)定的計(jì)算結(jié)果.另一方面，鑒于非規(guī)則、非結(jié)構(gòu)網(wǎng)格的優(yōu)勢(shì)，計(jì)算時(shí)氣云內(nèi)使用黏聲波方程，氣云外的區(qū)域按彈性波方程計(jì)算，可以有效節(jié)省黏聲波計(jì)算時(shí)的存儲(chǔ)空間，提高計(jì)算效率，是其他規(guī)則網(wǎng)格方法所不具備的.

5.3 Marmousi模型

圖7 (a)Marmousi速度模型;(b)根據(jù)速度計(jì)算得到的Q值模型Fig.7 (a)Marmousi velocity model；(b)Q model derived from the velocity model

圖8 傳播時(shí)刻T=1.2 s時(shí)的波場(chǎng)快照(a)黏聲波結(jié)果;(b)聲波結(jié)果(左半部分)與黏聲波結(jié)果(右半部分)在同一時(shí)刻的對(duì)比.黏聲波波場(chǎng)的幅值有明顯的衰減.Fig.8 The snapshots of pressure at T=1.2 s(a) The viscoacoustic wave;(b) The comparisons of the acoustic wave (the left part)and the viscoacoustic wave (the right part)at the same time.The amplitude of viscoacoustic wave field is significant reduced by the attenuation.

圖9 Marmousi模型模擬結(jié)果的聲波頻譜與黏聲波頻譜對(duì)比圖.圖中實(shí)線代表聲波頻譜，虛線代表黏聲波頻譜.黏性介質(zhì)下地震波能量向低頻方向移動(dòng)，高頻成分被強(qiáng)烈吸收Fig.9 The comparison of acoustic spectrum and viscoacoustic spectrum derived from the Marmousi Model.The solid line represents the acoustic spectrum,and the dash line represents the viscoacoustic spectrum.In the attenuation medium,the seismic energy moves toward to the low frequencies,and strongly absorbs the high frequency components

6 結(jié)論

本文發(fā)展了基于非規(guī)則化網(wǎng)格的時(shí)間域黏聲波數(shù)值模擬算法.文中采用有理式基函數(shù)描述頻率相關(guān)的體積模量，通過無窮范數(shù)構(gòu)建目標(biāo)函數(shù)，并結(jié)合模擬退火方法，實(shí)現(xiàn)了滿足衡Q要求的低階基函數(shù)參數(shù)的快速搜索，降低了計(jì)算中的內(nèi)存需求.同時(shí)，由引入n個(gè)關(guān)于中間變量的一階偏微分方程，替代常規(guī)時(shí)間域黏聲波數(shù)值模擬中的褶積運(yùn)算，極大地緩和了計(jì)算存儲(chǔ)問題，提高了計(jì)算效率.低階有理式表示的體積模量和格子法的結(jié)合，使得本文的數(shù)值模擬算法，既可精細(xì)地刻畫復(fù)雜邊界條件，也能兼顧計(jì)算精度和效率.簡(jiǎn)單模型中和理論解的分析比較證明了方法的精度及穩(wěn)定性.二維Marmousi模型算例表明，本文算法對(duì)復(fù)雜介質(zhì)模型具有很好適用性.本文方法也可延用到黏彈性波模擬及成像，進(jìn)一步提高地震勘探精度，有很好的應(yīng)用前景.

Blake R J,Bond L J,Downie A L.1982.Advances in numerical studies of elastic wave propagation and scattering.∥Reviewof Progress in Quantitative Nondestructive Evaluation.New York:Plenum Press,157-166.

Brossier R,Operto S,Virieux J,et al.2010.Frequency-domain numerical modelling of visco-acousticwaveswith finite-difference and finite-element discontinuous Galerkinmethods.∥Dissanayake DW.Acoustic Waves.SCIYO.

Carcione J M,Kosloff D,Kosloff R.1988.Wave propagation simulation in a linear viscoelastic medium.Geophysical Journal International,95(3):597-611.

Carcione J M.2010.A generalization of the Fourier pseudo spectral method.Geophysics,75(6):A53-A56.

Chen W,Holm S.2004.Fractional Laplacian time-space models for linear and nonlinear lossymedia exhibiting arbitrary frequency power-law dependency.The Journal of the Acoustical Society of America,115(4):1424-1430.

Day S M,Minster J B.1984.Numerical simulation of attenuated wavefields using a Padéapproximant method.Geophysical Journal International,78(1):105-118.

Emmerich H,Korn M.1987.Incorporation of attenuation into time-domain computations of seismic wave fields.Geophysics,52(9):1252-1264.

Gazdag J.1981.Modeling of the acoustic wave equation with transform methods.Geophysics,46(6):854-859.

Goffe W L,FerrierG D,Rogers J.1994.Global optimization of statistical functions with simulated annealing.Journal of Econometrics,60(1-2):65-99.

Kjartansson E.1979.Constant Q-wave propagation and attenuation.Journal of Geophysical Research:Solid Earth,84(B9):4737-4748.

Li Q Z.1993.High-resolution Seismic Data Processing (in Chinese).Beijing:Petroleum Industry Press,38.

Liao J P,Liu H X,Wang H Z,et al.2011.Two dimensional visco-acoustic wave modeling research in frequency-space domain.Progress in Geophysics (in Chinese),26(6):2075-2081,doi:10.3969/j.issn.1004-2903.2011.06.023.

Liu H P,Anderson D L,Kanamori H.1976.Velocity dispersion due to anelasticity:Implications for seismology and mantle composition.Geophysical Journal International,47(1):41-58.

Marfurt K J.1984.Accuracy of finite-difference and finite-element modeling of the scalar and elastic wave equations.Geophysics,49(5):533-549.

Operto S,Virieux J,Amestoy P,et al.2007.3D finite-difference frequency-domain modeling of visco-acoustic wave propagation using a massively parallel direct solver:A feasibility study.Geophysics,72(5):SM195-SM211.

Robertsson J O,Blanch J O,Symes W W.1994.Viscoelastic finite-difference modeling.Geophysics,59(9):1444-1456.

Saenger E H,Bohlen T.2004.Finite-difference modeling of viscoelastic and anisotropic wave propagation using the rotated staggered grid.Geophysics,69(2):583-591.

Wang R Z,Liu R,Zhao B X.2014.Numerical simulation of wave equation with segmented smooth curve boundaries.Chinese Journal of Geophysics (in Chinese),57(4):1199-1208,doi:10.6038/cjg20140417.

Xu Y,Zhang J F.2008.An irregular-grid perfectly matched layer absorbing boundary for seismic wave modeling.Chinese Journal of Geophysics (in Chinese),51(5):1520-1526,doi:10.3321/j.issn:0001-5733.2008.05.026.

Zhang J F,Liu T L.1999.P-SV-wave propagation in heterogeneous media:Grid method.Geophysical Journal International,136(2):431-438.

Zhang J F,Liu T L.2002.Elastic wave modelling in 3-D heterogeneous media:3D grid method.Geophysical Journal International,150(3):780-799.

Zhu T Y,Harris J M.2013.A constant-Q time-domain wave equation using the fractional Laplacian.∥SEG Technical Program Expanded Abstracts 2013.Society of Exploration Geophysicists,3417-3422.

Zhu T Y.2014.Time-reverse modelling of acoustic wave propagation in attenuating media.Geophysical Journal International,197(1):483-494.

Zhu T Y,Harris J M.2014.Modeling acoustic wave propagation in heterogeneous attenuating media using decoupled fractional Laplacians.Geophysics,79(3):T105-T116.

附中文參考文獻(xiàn)

李慶忠.1993.走向精確勘探的道路.北京：石油工業(yè)出版社，38.

廖建平,劉和秀,王華忠等.2011.二維頻率空間域黏聲波正演模擬研究.地球物理學(xué)進(jìn)展,26(6):2075-2081,doi:10.3969/j.issn.1004-2903.2011.06.023.

王忠仁,劉瑞,趙博雄.2014.分段光滑曲線邊界波動(dòng)方程數(shù)值模擬研究.地球物理學(xué)報(bào),57(4):1199-1208,doi:10.6038/cjg20140417.

吳玉,符力耘,陳高祥.2014.基于分?jǐn)?shù)階拉普拉斯算子解耦的粘聲介質(zhì)地震正演模擬與逆時(shí)偏移.∥2014年中國(guó)地球科學(xué)聯(lián)合學(xué)術(shù)年會(huì)-專題19:地震波傳播與成像論文集.北京:中國(guó)地球物理學(xué)會(huì).

徐義,張劍鋒.2008.地震波數(shù)值模擬的非規(guī)則網(wǎng)格PML吸收邊界.地球物理學(xué)報(bào),51(5):1520-1526,doi:10.3321/j.issn:0001-5733.2008.05.026.

(本文編輯 胡素芳)

An irregular grid method for acoustic modeling in inhomogeneous viscoelastic medium

DOU Hui1,2，ZHANG Jian-Feng1

1 Key Laboratory of Petroleum Resources Research，Institute of Geology and Geophysics， Chinese Academy of Sciences，Beijing 100029，China2 University of Chinese Academy of Sciences，Beijing 100049，China

Absorption and dispersion in earth media is an issue we should consider for wavefield simulation.In this paper,we present a time-domain numerical modeling algorithm for acoustic wave in inhomogeneous viscoelastic medium.On one hand,to describe the attenuation effect we use the constant-Q model,which can be frequency-independent.On the other hand for the wavefield modeling engine,we utilize the grid method,which uses irregular,unstructured grid to simulate acoustic wave propagation in non-uniform viscoelatsic media.

Time-domain；Seismic modeling；Irregular grid；Viscoelasticity；Constant-Q model；Rational function

豆輝，張劍鋒.2016.非均勻黏彈性介質(zhì)中聲波模擬的非規(guī)則網(wǎng)格方法.地球物理學(xué)報(bào)，59(11):4212-4222,

10.6038/cjg20161123.

Dou H，Zhang J F.2016.An irregular grid method for acoustic modeling in inhomogeneous viscoelastic medium.Chinese J.Geophys.(in Chinese),59(11):4212-4222,doi:10.6038/cjg20161123.

國(guó)家自然科學(xué)基金面上項(xiàng)目(41574135)和國(guó)家自然科學(xué)基金重點(diǎn)項(xiàng)目(41330316)聯(lián)合資助.

豆輝，女，1989年生，中國(guó)科學(xué)院地質(zhì)與地球物理研究所博士研究生，主要從事時(shí)間域數(shù)值模擬研究.E-mail:douhui@mail.iggcas.ac.cn

10.6038/cjg20161123

P631

2016-03-17，2016-09-19收修定稿

Specifically to resolve the convolution problem between strain and stress in the time-domain,we approximate the bulk modulus in the frequency domain with a set of rational functions,which meets the constant-Q requirement.To confirm the model parameters for constant-Q,we apply the L-∞ Norm to reconstruct the objective function,and choose the simulated annealing (SA) method to finish the searching.As a result,we find the optimum solution with low order,like n=2 or 3,which simultaneously meets the constant-Q requirement very well,and reduces the memory demand obviously.In addition,the irregular,unstructured grid,which is used to discretize the velocity model and corresponding Q model,can give detailed descriptions and differentiations when complex interfaces are present,without sacrificing the accuracy and efficiency.

As for the numerical examples,we first validate the scheme proposed using simple model and compare the results with analytic solutions.Result with good consistency is obtained regarding the attenuation effect.Moreover the Marmousi model validation proves that the proposed irregular,unstructured gird based simulation method together with the constant-Q model is effective and efficient to describe the acoustic wave propagation in inhomogeneous viscoelastic medium.