豆輝，張劍鋒

1 中國科學院地質與地球物理研究所，中國科學院油氣資源研究重點實驗室，北京 1000292 中國科學院大學，北京 100049

?

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

豆輝1,2，張劍鋒1

1 中國科學院地質與地球物理研究所，中國科學院油氣資源研究重點實驗室，北京 1000292 中國科學院大學，北京 100049

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

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

1 引言

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

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

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

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

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

2 非規(guī)則、非結構網格聲波模擬的格子法

本文數值解法是基于聲波模擬的格子法(Zhang and Liu，1999，2002)展開工作.格子法采用非規(guī)則、非結構化網格離散速度模型，既可以精細地刻畫復雜邊界條件、速度及Q值的變化，也具有與有限差分法相當的計算效率.

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

(1)

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

(2)

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

(3)

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

(4)

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

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

最后，利用中心差分格式離散方程式(3)左端中對時間的二階偏導數，我們就可以實現聲波波場的迭代更新.

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

3.1 體積模量的有理函數展開

黏彈性介質在頻率域中有如下應力應變關系：

(5)

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

(6)

這意味著在時間域計算時，需要存儲歷史時刻的應變值，會導致極大的計算負擔，是數值計算很不樂見的.所以，我們需要考慮從其他方面對其做簡化近似處理.

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

(7)

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

(8)

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

(9)

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

3.2 衡Q模型

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

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

(10)

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

(11)

圖2 階數(a)n=2，(b)n=3時，Q-1=0.05的目標值(實線)和Q-1(ω)曲線擬合結果(虛線)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)時有2n+1個模型參數，還比較少，我們可以采用蒙特卡洛、模擬退火、遺傳算法等全局搜索方法.本文采用的是優(yōu)化的模擬退火方法(Goffeetal.，1994)，可實現全局的快速搜索，求得最優(yōu)解.然后將最優(yōu)解對應的參數結果列表記錄下來，便于以后模擬計算時查表使用.最終得到的結果和目標值的擬合效果，可以參見圖2所示，圖中實線表示目標值Q0=20.0，虛線為Q-1(ω)擬合結果.圖2中，分別對應階數n=2，n=3時Q(ω)擬合得到的曲線結果和目標值的比較.可以看出，在頻帶3～100 Hz內，n=2時即可得到不錯的擬合效果；n=3時，Q(ω)曲線可以更加貼合目標值，擬合效果更優(yōu).本文數值模擬的簡單算例，采用的就是階數n=3時的參數結果.

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

(12)

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

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

(13)

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

(14)

對式(13)和式(14)做Fourier反變換，得到時間域黏彈性介質的波動方程：

(15)

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

(16)

(17)

(18)

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

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

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

(20)

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

(21)

(22)

(23)

式中上標q代表時間，下標l表示節(jié)點k周圍第l個三角形單元.式(21)-(23)構成了非規(guī)則網格黏聲波方法的應變和壓強更新基礎.由q時刻的壓強p代入式(22)可求得q+0.5時刻各節(jié)點處的形變ε，再結合式(23)求得q+0.5時刻各節(jié)點處的中間變量φj，最后由式(21)求得q+1時刻的壓強波場p,完成對壓強p在時間上的更新.

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

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

5 數值算例

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

5.1 均勻介質模型

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

圖3 不同炮檢距時，數值解(虛線)和解析解(實線)的比較.檢波器距離震源距離依次為(a) 600，(b) 1600，(c) 3200 m.二者結果吻合較好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 不同炮檢距時

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

(24)

5.1.2 不同Q值時

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

圖4 均勻衰減介質中，T=200 ms時刻的波場快照圖，圖中的四個部分分別對應不同的Q值：(a)1000，(b)100，(c)30，(d)10.震源位于模型的中間位置.隨著Q變小，地震波波場有更強的幅值衰減和相位延遲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 震源位于點(2.0,1.0)km，檢波器位于(2.6,1.0)km處，不同的Q值：(a)1000，(b)100，(c)30，(d)10時，對應的數值解(虛線)和解析解(實線)的比較.它們的結果都能很好的吻合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

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

圖6 (a)370 ms,(b)450 ms時刻的波場快照.在圖中白色圓圈代表的氣云周圍有反射產生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 氣云模型

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

5.3 Marmousi模型

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

圖8 傳播時刻T=1.2 s時的波場快照(a)黏聲波結果;(b)聲波結果(左半部分)與黏聲波結果(右半部分)在同一時刻的對比.黏聲波波場的幅值有明顯的衰減.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模型模擬結果的聲波頻譜與黏聲波頻譜對比圖.圖中實線代表聲波頻譜，虛線代表黏聲波頻譜.黏性介質下地震波能量向低頻方向移動，高頻成分被強烈吸收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 結論

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

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.

附中文參考文獻

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

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

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

吳玉,符力耘,陳高祥.2014.基于分數階拉普拉斯算子解耦的粘聲介質地震正演模擬與逆時偏移.∥2014年中國地球科學聯合學術年會-專題19:地震波傳播與成像論文集.北京:中國地球物理學會.

徐義,張劍鋒.2008.地震波數值模擬的非規(guī)則網格PML吸收邊界.地球物理學報,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.非均勻黏彈性介質中聲波模擬的非規(guī)則網格方法.地球物理學報，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.

國家自然科學基金面上項目(41574135)和國家自然科學基金重點項目(41330316)聯合資助.

豆輝，女，1989年生，中國科學院地質與地球物理研究所博士研究生，主要從事時間域數值模擬研究.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.