謝俊法 孫成禹 伍敦仕 喬志浩

1) 中國(guó)山東青島266580中國(guó)石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院2) 中國(guó)蘭州730020中國(guó)石油勘探開發(fā)研究院西北分院

?

基于無(wú)分裂復(fù)頻移卷積完全匹配層邊界的黏彈介質(zhì)勒夫波模擬

謝俊法1,2)孫成禹1),伍敦仕1)喬志浩1)

1) 中國(guó)山東青島266580中國(guó)石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院2) 中國(guó)蘭州730020中國(guó)石油勘探開發(fā)研究院西北分院

本文建立了無(wú)分裂復(fù)頻移卷積完全匹配層(CFS-CPML)吸收邊界條件， 利用交錯(cuò)網(wǎng)格下的高精度有限差分格式對(duì)黏彈性介質(zhì)中的勒夫波場(chǎng)進(jìn)行了數(shù)值模擬； 分析了松弛機(jī)制個(gè)數(shù)對(duì)品質(zhì)因子擬合精度的影響， 驗(yàn)證了CFS-CPML邊界條件對(duì)大角度掠射波的吸收效果． 數(shù)值結(jié)果表明： 本文方法所使用的5個(gè)松弛機(jī)制和空間4階差分精度， 即可在保證計(jì)算效率的前提下滿足目前理論研究的需要； 隨著品質(zhì)因子的減小， 頻散特征曲線的相速度逐漸向增高的方向偏離理論頻散特征曲線的相速度， 且各模式的高頻能量也隨之減弱． 本文結(jié)果可為發(fā)展高精度的面波反演方法提供必要的理論依據(jù)．

勒夫波 正演模擬 復(fù)頻移卷積完全匹配層(CFS-CPML) 吸收邊界條件 頻散特征

引言

勒夫波是在低速層分界面附近傳播的一種SH型不均勻平面波， 包含了近地表的構(gòu)造信息， 因此可以利用勒夫波獲取近地表的橫波速度、 品質(zhì)因子等地層信息(Lee， Solomon， 1975， 1978； Xiaetal， 2002， 2013； Winsborrowetal， 2003， 2005； Safanietal， 2005， 2006； Dongetal， 2013)．

地震波在實(shí)際地層的傳播過(guò)程中， 伴隨著能量的吸收與衰減， 因而采用線性黏彈性的假設(shè)能夠更準(zhǔn)確地表征地震波尤其是面波在近地表衰減作用下的傳播特征． 一般利用品質(zhì)因子Q來(lái)表征介質(zhì)的黏彈性， 且其在地震頻帶范圍內(nèi)可近似為常數(shù)(Kanamori， Anderson， 1977)． 時(shí)間域常Q黏彈性介質(zhì)中的地震波模擬， 通常采用廣義標(biāo)準(zhǔn)線性體(Genera-lized Zener body， 簡(jiǎn)寫為GZB)模型(Carcioneetal， 1988； Carcione， 1993)或廣義麥克斯韋爾(generalized Maxwell body define by Emmerich and Korn， 簡(jiǎn)寫為GMB-EK)模型(Emmerich， Korn， 1987； Kristek， Moczo， 2003)實(shí)現(xiàn)常Q近似， Moczo和Kristek(2005)已證明這兩種模型是等價(jià)的． 關(guān)于常Q黏彈性介質(zhì)中體波模擬的研究較多(Liuetal， 1976； 孫成禹等， 2010)， 而針對(duì)面波的研究則相對(duì)較少. Carcione(1992)曾采用偽譜法對(duì)黏彈性半空間情形下的瑞雷波進(jìn)行正演， 并采取廣義標(biāo)準(zhǔn)線性體模型進(jìn)行常Q擬合， 但其松弛函數(shù)遺漏了1/L(L為標(biāo)準(zhǔn)線性體個(gè)數(shù)或松弛機(jī)制個(gè)數(shù))因子， 該疏漏后來(lái)在其專著中得以修正(Carcione， 2001)． Zhang等(2011)同樣采用偽譜法結(jié)合廣義標(biāo)準(zhǔn)線性體模型進(jìn)一步分析了兩層介質(zhì)情況下瑞雷面波的頻散特征變化， 但是由于其采用的模型較為簡(jiǎn)單， 高階模式的變化規(guī)律不明顯． 高靜懷等(2012)基于位移-應(yīng)力旋轉(zhuǎn)交錯(cuò)網(wǎng)格有限差分法研究了黏彈性介質(zhì)中瑞雷波的正演， 并與完全彈性情形進(jìn)行了對(duì)比， 但其采用的松弛函數(shù)同樣也沒有1/L因子． 就勒夫波而言， Boore(1970)采用常規(guī)網(wǎng)格有限差分法實(shí)現(xiàn)了完全彈性橫向非均勻介質(zhì)中的勒夫波數(shù)值模擬； Luo等(2010)基于標(biāo)準(zhǔn)交錯(cuò)網(wǎng)格有限差分法研究了完全彈性介質(zhì)中3種典型地層的勒夫波頻散特征． 而關(guān)于黏彈性常Q介質(zhì)中勒夫波數(shù)值模擬的研究， 特別是Q值對(duì)勒夫波頻散特征的影響程度， 尚不多見．

波場(chǎng)數(shù)值模擬需要高效的邊界反射處理技術(shù)． Bérenger(1994)在進(jìn)行電磁波模擬時(shí)提出完全匹配層技術(shù)(perfectly matched layer， 簡(jiǎn)寫為PML). 由于PML具有良好的吸收效果， 目前已成為地震波數(shù)值模擬中最流行的邊界反射處理技術(shù)． 早期的PML多采用波場(chǎng)分裂方式， 稱為分裂式完全匹配層(Split-PML)(Collino， Tsogka， 2011)， 相對(duì)于旁軸近似、 指數(shù)衰減等吸收邊界處理技術(shù)(Clayton， Engquist， 1977； Cerjanetal， 1985)， 其吸收性能有了很大的提高， 但仍然存在某些局限性， 例如， 當(dāng)震源離PML層很近或接收排列過(guò)長(zhǎng)時(shí)， 不能有效地吸收以掠射角入射的波． 為了克服該問(wèn)題， 無(wú)需波場(chǎng)分裂的PML技術(shù)得以發(fā)展， 其中比較典型的為復(fù)頻移卷積完全匹配層技術(shù)(complex frequency shifted convolutional perfectly matched layer， 簡(jiǎn)寫為CFS-CPML)． Komatitsch和Martin(2007)比較了完全彈性介質(zhì)中CFS-CPML與Split-PML在吸收掠射角入射波時(shí)的特性， Martin和Komatitsch(2009)進(jìn)一步比較了黏彈性介質(zhì)中兩種邊界的吸收特性． 國(guó)內(nèi)多位研究人員對(duì)吸收邊界問(wèn)題也進(jìn)行了大量的研究工作(張顯文等， 2009； 張魯新等， 2010； 田坤等， 2013； Wangetal， 2013)， 結(jié)果表明， 采用CFS-CPML能夠有效解決掠射情況下PML吸收效果差的問(wèn)題， 而且實(shí)現(xiàn)過(guò)程中無(wú)需分裂變量， 提高了計(jì)算效率．

鑒于實(shí)際地層并非完全彈性， 尤其是近地表介質(zhì)， 其黏彈性特征較為明顯. 為了模擬黏彈性介質(zhì)中的勒夫波， 本文將首先利用5個(gè)松弛機(jī)制的廣義麥克斯韋爾模型來(lái)刻畫黏彈性特征， 采取無(wú)需分裂的復(fù)頻移卷積完全匹配層技術(shù)處理人工邊界反射， 結(jié)合時(shí)間2階、 空間4階精度的標(biāo)準(zhǔn)交錯(cuò)網(wǎng)格有限差分法， 以實(shí)現(xiàn)黏彈性常Q介質(zhì)中的勒夫波正演模擬． 其次， 分析不同Q值對(duì)勒夫面波頻散特征的影響， 以及速度遞增、 夾高速層和夾低速層的復(fù)雜地層的勒夫波頻散特征， 為利用勒夫波獲取近地表信息提供必要的理論依據(jù)．

1 勒夫波一階速度-應(yīng)力-記憶變量方程推導(dǎo)

在完全彈性條件下， 應(yīng)力與應(yīng)變之間是一種瞬時(shí)對(duì)應(yīng)關(guān)系， 不隨時(shí)間變化； 而在黏彈性條件下， 應(yīng)力與應(yīng)變均為時(shí)間的函數(shù)， 三維各向同性線性黏彈性介質(zhì)的本構(gòu)方程可表示為(Robertssonetal， 1994)：

(1)

式中， “*”表示時(shí)域卷積，ψπ(t)和ψμ(t)為應(yīng)力松弛函數(shù)，σ為應(yīng)力，ε為應(yīng)變，t為時(shí)間， δ為克羅內(nèi)克函數(shù)， 下標(biāo)遵循愛因斯坦求和約定．

通常基于廣義麥克斯韋爾模型(Emmerich， Korn， 1987； Kristek， Moczo， 2003)來(lái)實(shí)現(xiàn)常Q近似， 基于該模型得到的松弛函數(shù)為

(2)

(3)

式中:μu為未松弛模量;ξxyl和ξzyl為對(duì)應(yīng)第l個(gè)松弛機(jī)制的記憶變量；v為速度分量， 單位為m/s;ρ為密度， 單位為kg/m3．

2 勒夫波黏彈方程的CFS-CPML吸收邊界條件

高效的吸收邊界是進(jìn)行波場(chǎng)數(shù)值模擬必不可少的條件， 而采用無(wú)分裂波場(chǎng)的CFS-CPML吸收邊界， 避免了非物理的波場(chǎng)分裂所引起的數(shù)值誤差， 同時(shí)引入了復(fù)頻移技術(shù)， 提高了對(duì)高角度掠射波場(chǎng)的吸收效果． 對(duì)式(3)進(jìn)行傅里葉變換， 得到其頻域形式：

(4)

式中 “^”表示對(duì)應(yīng)的頻域物理量． 經(jīng)典完全匹配層(perfectly matched layer， 簡(jiǎn)寫為PML)技術(shù)實(shí)際上是對(duì)PML層內(nèi)的參數(shù)進(jìn)行復(fù)坐標(biāo)拉伸， 即

(5)

(6)

式中sx和sz為復(fù)拉伸函數(shù)， 具體表示為

(7)

經(jīng)典PML技術(shù)通常采用波場(chǎng)分裂的方式， 對(duì)質(zhì)點(diǎn)的速度分量和應(yīng)力分量沿坐標(biāo)軸方向進(jìn)行分裂(Qinetal， 2009) ． 當(dāng)震源離PML層很近或接收排列過(guò)長(zhǎng)時(shí)， 不能有效地吸收以掠射角入射的波． 為解決該問(wèn)題， CFS-CPML技術(shù)(Komatitsch， Martin， 2007； Martin， Komatitsch， 2009； 張顯文等， 2009； Zhang， Yang， 2010； 張魯新等， 2010； 田坤等， 2013； 王華等， 2013)引入兩個(gè)頻移因子αx，αz和兩個(gè)收縮因子kx，kz， 將式(7)推廣為

(8)

式中，kx≥1，kz≥1，αx≥0，αz≥0． 頻移因子αx和αz的引入有效避免了衰減系數(shù)出現(xiàn)奇異值； 收縮因子kx和kz的引入使得地震波進(jìn)入PML區(qū)域后， 傳播方向沿該區(qū)域的法方向彎曲， 這樣高角度入射波能夠更加深入PML區(qū)域中并逐漸衰減(Zhang， Yang， 2010)． 進(jìn)一步整理式(8)為

(9)

將式(9)代入式(4)， 經(jīng)傅里葉反變換， 得到基于CFS-CPML技術(shù)實(shí)現(xiàn)的一階速度-應(yīng)力-記憶變量時(shí)域方程如下：

(10)

其中

(11)

為4個(gè)內(nèi)部變量. 沿x方向的衰減系數(shù)dx， 收縮因子kx， 頻移因子αx分別為(Zhang， Yang， 2010)

(12)

圖1 不同物理參數(shù)在交錯(cuò)網(wǎng)格 單元中的布局方式 Fig.1 Distribution of different physical parameters within the staggered-grid cells

式中：x表示網(wǎng)格節(jié)點(diǎn)到PML層內(nèi)界面的距離；L為PML層的厚度；d0，k0，α0為常數(shù)， 通常取d0=-3cSlogRc/2L，cS為橫波速度，Rc為反射系數(shù)． 當(dāng)PML層的厚度為10個(gè)網(wǎng)格節(jié)點(diǎn)時(shí)， 可取Rc=0.1%，k0=2，α0=πf0，f0為震源子波的主頻(Collino， Tsogka， 2001)． 沿z方向的處理辦法與x方向類似．

由于CFS-CPML技術(shù)處理邊界不需要進(jìn)行波場(chǎng)分裂， 降低了存儲(chǔ)和求解的復(fù)雜程度， 其計(jì)算效率比使用Split-PML技術(shù)高．

3 標(biāo)準(zhǔn)交錯(cuò)網(wǎng)格離散差分方式

根據(jù)圖1所示的網(wǎng)格交錯(cuò)方式， 對(duì)式(10)進(jìn)行離散， 即可實(shí)現(xiàn)黏彈性SH波的數(shù)值模擬． 限于篇幅， 這里僅給出速度分量和內(nèi)部變量的時(shí)間2階、 空間4階離散形式：

(13)

(14)

其中

(15)

式(10)中的其它變量， 如σxy，σzy，ξxyl，ξzyl， 可類似導(dǎo)出．

利用顯式有限差分進(jìn)行數(shù)值模擬時(shí)， 須考慮計(jì)算過(guò)程的穩(wěn)定性． 以黏彈介質(zhì)中的未松弛(高頻極限)速度代替彈性速度(田坤， 2014)， 得到

(16)

可以將此作為黏彈性介質(zhì)穩(wěn)定性的判斷依據(jù)． 式中Δx和Δz為橫向和縱向的空間網(wǎng)格采樣步長(zhǎng)， Δt為時(shí)間采樣間隔，cu為高頻極限速度，an為相應(yīng)的交錯(cuò)網(wǎng)格差分算子系數(shù)．

4 數(shù)值模擬算例與波場(chǎng)分析

4.1 松弛機(jī)制個(gè)數(shù)的選取

利用式(2)的松弛函數(shù)進(jìn)行常Q近似時(shí)， 松弛機(jī)制個(gè)數(shù)對(duì)擬合精度具有較大的影響， 松弛機(jī)制個(gè)數(shù)越多， 擬合精度越高， 但相應(yīng)的計(jì)算量也越大； 因此， 在保證擬合精度足夠高的前提下， 需要盡量減少松弛機(jī)制的個(gè)數(shù)， 從而減小計(jì)算量． 圖2給出了分別采用3個(gè)和5個(gè)松弛機(jī)制在1—10 Hz范圍內(nèi)的常Q擬合結(jié)果. 可以看出， 采用3個(gè)松弛機(jī)制的擬合曲線與Q=20的曲線相差較大， 而采用5個(gè)松弛機(jī)制的擬合曲線則與Q=20的曲線吻合得較好．

圖2 松弛機(jī)制個(gè)數(shù)與常Q擬合精度的關(guān)系

number and constant Q fitting precision

圖3給出了1—10 Hz范圍內(nèi)， 采用5個(gè)松弛機(jī)制進(jìn)行常Q擬合的數(shù)值相速度與使用Futterman(1962)公式計(jì)算所得理論相速度的對(duì)比結(jié)果， 其中相速度的變化范圍為2560—2680 m/s， 其最大相對(duì)誤差為0.15%， 這充分說(shuō)明采用5個(gè)松弛機(jī)制得到的擬合結(jié)果具有較高的精度， 故本文的正演模擬中選取5個(gè)松弛機(jī)制．

圖4 均勻各向同性線性黏彈性介質(zhì)模型示意圖Fig.4 Schematic diagram of the homogeneous isotropic linear viscoelastic medium

4.2 邊界條件的吸收效果對(duì)比

圖5 基于黏彈性介質(zhì)模型采用Split-PML吸收邊界所獲取的0.14 s (上)，0.3 s (中)， 0.54 s (下)時(shí)刻的垂直分量波場(chǎng)快照

圖7給出了采用Split-PML和CFS-CPML兩種吸收邊界獲得的炮記錄對(duì)比， 兩側(cè)各40道為吸收邊界， 考慮到邊界反射較弱， 制圖時(shí)對(duì)圖7a和圖7b均加了4%的增益． 可以看到： 圖7a中炮檢距較大的地震道存在頂部邊界的反射(箭頭所示區(qū)域)， 而炮檢距較小的地震道則沒有， 這是由于大炮檢距處波入射到匹配層的入射角過(guò)大， 匹配層吸收不完全所致； 圖7b中沒有頂部邊界的反射， 說(shuō)明CFS-CPML吸收邊界對(duì)大角度掠射到邊界的地震波具有較好的吸收效果．

圖6 基于黏彈介質(zhì)模型采用CFS-CPML吸收邊界所獲取的0.14 s (上)，0.3 s (中)， 0.54 s (下)時(shí)刻的垂直分量波場(chǎng)快照

圖7 炮記錄的對(duì)比(a) 采用Split-PML技術(shù)， 箭頭所指為頂部邊界的反射； (b) 采用CFS-CPML技術(shù)

將圖4所示的介質(zhì)模型擴(kuò)大成6000 m×6000 m， 即網(wǎng)格點(diǎn)數(shù)為3000×3000， 激發(fā)點(diǎn)設(shè)在(3000 m， 3000 m)點(diǎn)處， 接收線的深度設(shè)為3002 m． 在這種條件下， 當(dāng)計(jì)算時(shí)間tcal≤0.54 s時(shí)， 直達(dá)波未到達(dá)邊界， 屬于無(wú)邊界反射記錄， 本文將該記錄作為評(píng)價(jià)邊界吸收效果的標(biāo)準(zhǔn)． 選取圖7a和圖7b中炮檢距為零的地震數(shù)據(jù)進(jìn)行波形顯示， 如圖8a所示， 可以看出， 采用Split-PML和CFS-CPML吸收邊界的波形與無(wú)邊界反射記錄非常一致． 圖8a下部放大圖中所顯示的0.05—0.1 s內(nèi)振幅為10-3數(shù)量級(jí)， 存在邊界反射的波形與無(wú)邊界反射的波形幾乎一致， 說(shuō)明地震波以小角度入射到邊界時(shí)， 采用Split-PML和CFS-CPML吸收邊界均有較好的吸收效果． 選取圖7a和圖7b中炮檢距為400 m的地震數(shù)據(jù)進(jìn)行波形顯示， 如圖8b所示. 從圖8b下圖的直達(dá)波尾部放大圖可以看出， 采用CFS-CPML吸收邊界的波形更接近無(wú)邊界反射的波形， 說(shuō)明對(duì)于大角度掠射的地震波， 采用CFS-CPML吸收邊界的處理效果較好．

圖8 單道波形的對(duì)比

4.3 波場(chǎng)模擬與分析

4.3.1 雙層介質(zhì)模型的正演模擬與分析

由于勒夫波不存在于均勻半空間中， 其正演模擬需采用多層介質(zhì)模型． 建立最簡(jiǎn)單的兩層介質(zhì)模型， 參數(shù)見表1. 由于密度對(duì)勒夫波的影響較小， 所以模型中兩層介質(zhì)的密度均設(shè)為2000 kg/m3． 整個(gè)模型由801×200個(gè)網(wǎng)格點(diǎn)構(gòu)成， 網(wǎng)格大小為Δx=Δz=0.5 m， 時(shí)間步長(zhǎng)為Δt=0.5 ms， CFS-CPML層由10個(gè)網(wǎng)格點(diǎn)構(gòu)成； 震源函數(shù)采用主頻為f0=40 Hz的雷克子波， 反射系數(shù)Rc=0.1%， 收縮因子k0=2， 頻移因子α0=πf0， 采用5個(gè)松弛機(jī)制在1—100 Hz范圍內(nèi)進(jìn)行常Q擬合． 鑒于Q>150時(shí)， 黏彈性衰減效應(yīng)很小， 因此當(dāng)?shù)貙観≥200時(shí)， 該地層可視為完全彈性介質(zhì)． 保持底層半空間介質(zhì)為完全彈性， 依次使第一層的橫波品質(zhì)因子QS為200， 40， 20和10， 采用時(shí)間2階、 空間4階的有限差分法進(jìn)行正演得到4個(gè)勒夫波記錄， 通過(guò)傾斜疊加算法(Parketal， 1999)計(jì)算其所對(duì)應(yīng)的頻散特征(圖9a--d)， 并將其與理論頻散特征進(jìn)行對(duì)比以考察表層Q值變化對(duì)勒夫波頻散特征的影響．

表1 兩層介質(zhì)模型參數(shù)

當(dāng)?shù)谝粚咏橘|(zhì)的橫波品質(zhì)因子QS值為200時(shí)， 兩層介質(zhì)均為完全彈性， 其勒夫波地震記錄及相應(yīng)的頻散特征如圖9a所示． 從勒夫波記錄中可以看出CFS-CPML邊界的吸收效果較好． 從圖9a給出的頻散特征(彩色區(qū)域)可以看到： 該記錄中包含基階、 第一高階、 第二高階和第三高階模式的勒夫波； 對(duì)每個(gè)模式而言， 隨著頻率的增大， 相速度逐漸減小并趨于表層橫波速度300 m/s； 隨著頻率的減小， 相速度趨于底層半空間橫波速度600 m/s． 圖9a中的白色曲線為根據(jù)矩陣傳遞算法(Haskell， 1953)正演得到的理論勒夫波頻散曲線， 很容易看到， 各模式的頻散曲線與其吻合得較好， 說(shuō)明正演算法能夠準(zhǔn)確地模擬該地層條件下的勒夫波． 此外， 理論上已證明， 在兩層完全彈性介質(zhì)的情況下， 勒夫波高階模式存在截止頻率， 即當(dāng)頻率低于截止頻率時(shí)不存在高階模式(Aki， Richards， 2002)， 這在圖9a中也得以體現(xiàn)， 即當(dāng)頻率低于20 Hz時(shí)只存在基階模式， 不存在高階模式．

圖9b給出了當(dāng)?shù)谝粚咏橘|(zhì)QS值為40時(shí)的勒夫波地震記錄及其相應(yīng)頻散特征． 可以看出， 4個(gè)模式的高頻成分損耗十分明顯， 從約80 Hz降低到近似70 Hz， 這是由于吸收衰減作用所致． 此外， 由于黏彈性所導(dǎo)致的固有頻散， 各模式頻散曲線的最低相速度均高于完全彈性情況下的理論相速度． 當(dāng)?shù)谝粚拥腝S值進(jìn)一步減小到20(圖9c)甚至10(圖9d)時(shí)， 從勒夫波記錄中可以看出高頻成分較完全彈性情況下的地震記錄損失較嚴(yán)重， 特別是第601道之后的遠(yuǎn)偏移距部分更為嚴(yán)重． 從圖9c頻散圖中也可以看到， 此時(shí)第三高階模式高頻能量已近乎完全衰減， 第二高階模式的能量也僅集中在40—60 Hz之間的狹窄頻帶內(nèi)， 基階和第一高階模式的最高頻率降低至60 Hz附近． 若表層介質(zhì)QS值減小到10(圖9d)， 第三高階模式高頻能量則完全衰減， 第二高階模式的最高頻率降低至50 Hz左右， 基階和第一高階模式的最高頻率也降低， 并且相速度在高頻時(shí)趨于340 m/s， 與完全彈性時(shí)的相速度300 m/s的相對(duì)偏差可達(dá)13.3%．

圖9 雙層介質(zhì)模型中的勒夫波記錄(左)及其頻散特征圖(右)頻散圖中的白色曲線為根據(jù)矩陣傳遞算法(Haskell， 1953)正演所得的理論勒夫波頻散曲線(a) 上層QS=200； (b) 上層QS=40； (c) 上層QS=20; (d) 上層QS=10

綜上可見， 勒夫波頻散曲線中各模式的最低相速度隨第一層介質(zhì)QS值的減小而增大． 然而， 實(shí)際地層并不是完全彈性的， 近地表介質(zhì)尤其如此， 因此， 在利用面波獲取近地表速度時(shí)， 第一層介質(zhì)的速度受頻散曲線最低相速度的影響較大， 導(dǎo)致利用面波所獲取的第一層速度與真實(shí)速度存在一定的誤差， 且QS值越小， 該誤差越大．

4.3.2 4層介質(zhì)模型的正演模擬與分析

近地表介質(zhì)復(fù)雜多變， 垂向上的速度并非總是遞增， 且近地表壓實(shí)程度較低， 速度和Q值通常較?。?為便于對(duì)比， 研究多層黏彈介質(zhì)的勒夫波特征時(shí)， 保持底層半空間介質(zhì)為完全彈性， 對(duì)3種黏彈介質(zhì)進(jìn)行組合， 構(gòu)成4層介質(zhì)模型， 具體參數(shù)如表2所示． 橫波速度為400， 600， 800 m/s和1000 m/s時(shí)， 橫波品質(zhì)因子QS分別為10， 20， 50， 200． 除模型參數(shù)設(shè)置不同外， 4層介質(zhì)模型的數(shù)值模擬參數(shù)與表1所示的雙層介質(zhì)模型一致． 從圖9a可知彈性介質(zhì)的勒夫波頻散曲線與理論曲線吻合， 因此僅分析黏彈介質(zhì)的勒夫波頻散曲線特征， 彈性介質(zhì)的頻散曲線特征可根據(jù)理論曲線作出判斷． 下面分速度遞增模型、 夾高速層模型和夾低速層模型等3種情況逐一討論．

表2 4層介質(zhì)模型參數(shù)

圖10a--c分別給出了速度遞增模型、 夾高速層模型和夾低速層模型的勒夫波記錄及其相應(yīng)的頻散曲線， 其中， 彩色區(qū)域是從勒夫波記錄中采用傾斜疊加算法所提取的頻散曲線， 白色線條是理論頻散曲線． 可以看出： 速度遞增模型和夾高速層模型的頻散曲線都只存在基階模式和第一高階模式， 基階模式的能量最強(qiáng)， 且隨著頻率的增加， 每個(gè)模式的頻散曲線與理論頻散曲線的吻合程度逐漸變差； 與速度遞增模型的頻散曲線(圖10a)相比， 夾高速層模型頻散曲線(圖10b)的基階模式彎曲程度更小， 第一高階模式的能量更弱； 夾低速層模型頻散曲線(圖10c)尚存在第二高階模式， 但其能量較弱， 總體頻率范圍和各模式的能量較速度遞增模型增加， 且基階模式的頻率范圍更小， 與理論頻散曲線的吻合程度更高， 此外， 第一高階模式的頻率范圍更大．

由于黏彈性介質(zhì)對(duì)高頻能量的吸收較其對(duì)低頻能量的吸收強(qiáng)， 因此頻散曲線中各模式高頻部分與理論曲線的吻合程度相對(duì)較差． 勒夫波是由一次波和多次反射的SH波相長(zhǎng)干涉產(chǎn)生的， 基于幾種介質(zhì)不同組合形成的模型得到的勒夫波頻散曲線， 其各階模式存在差異， 即速度遞增模型、 夾高速層模型和夾低速層模型的基階模式和第一高階在能量分布、 頻率范圍和彎曲程度上都存在差異． 反之， 勒夫波包含了地下構(gòu)造信息， 其頻散曲線中各模式的能量分布、 頻率范圍和彎曲程度是地下構(gòu)造信息的一種反映， 采用適當(dāng)?shù)姆椒軌驈睦辗虿ㄖ蝎@取地下的構(gòu)造信息， 目前常用的做法是拾取頻散特征曲線， 再結(jié)合最小二乘或者遺傳算法等優(yōu)化方法來(lái)反演得到地下橫波速度等底層信息(Xiaetal， 2002， 2013； Safanietal， 2005， 2006)．

圖10 4層介質(zhì)中勒夫波記錄(左)及其頻散特征圖(右)頻散圖中的白色曲線為根據(jù)矩陣傳遞算法(Haskell， 1953)正演所得的理論勒夫波頻散曲線(a) 速度遞增模型； (b) 夾高速層模型； (c) 夾低速層模型

5 討論與結(jié)論

本文以無(wú)分裂復(fù)頻移卷積完全匹配層為吸收邊界, 利用高精度交錯(cuò)網(wǎng)格有限差分， 推導(dǎo)了黏彈性介質(zhì)波動(dòng)方程的一階速度-應(yīng)力-記憶變量方程， 并給出了數(shù)值解法， 實(shí)現(xiàn)了黏彈性勒夫波的數(shù)值模擬. 通過(guò)模型試算， 得到結(jié)論如下：

1) 相對(duì)于傳統(tǒng)的分列式完全匹配層邊界條件， CFS-CPML邊界條件可以更加有效地吸收掠射到邊界的地震波， 且無(wú)需進(jìn)行波場(chǎng)分裂， 從而提高了計(jì)算效率．

2) 本文試驗(yàn)表明， 使用5個(gè)松弛機(jī)制和4階差分精度， 能夠在保證計(jì)算效率的前提下滿足目前理論研究的需要．

3) 介質(zhì)的黏彈性使頻散特征曲線的相速度高于理論頻散特征曲線的相速度， 各模式的高頻成分能量減弱， 且Q值越小， 頻散特征曲線與理論頻散特征曲線的偏差越大， 各模式高頻成分的能量越弱．

實(shí)際的近地表地層具有較強(qiáng)的黏彈性， 本文發(fā)展了基于無(wú)分裂復(fù)頻移卷積完全匹配層邊界的黏彈性勒夫波數(shù)值模擬方法， 模擬的頻散特征更加接近真實(shí)地層情況， 可為面波反演提供理論依據(jù)．

高靜懷， 何洋洋， 馬逸塵. 2012. 黏彈性與彈性介質(zhì)中Rayleigh面波特性對(duì)比研究[J]. 地球物理學(xué)報(bào)， 55(1): 207--218.

Gao J H， He Y Y， Ma Y C. 2012. Comparison of the Rayleigh wave in elastic and viscoelastic media[J].ChineseJournalofGeophysics， 55(1): 207--218 (in Chinese).

孫成禹， 肖云飛， 印興耀， 彭洪超. 2010. 黏彈介質(zhì)波動(dòng)方程有限差分解的穩(wěn)定性研究[J]. 地震學(xué)報(bào)， 32(2): 147--156.

Sun C Y， Xiao Y F， Yin X Y， Peng H C. 2010. Study on the stability of finite difference solution of visco-elastic wave equations[J].ActaSeismologicaSinica， 32(2): 147--156 (in Chinese).

田坤. 2014. 黏性介質(zhì)正演及逆時(shí)偏移成像方法研究[D]. 青島: 中國(guó)石油大學(xué)(華東): 35--39.

Tian K. 2014.StudyonMethodofForwardModelingandReverseTimeMigrationinViscousMedia[D]. Qingdao: China University of Petroleum: 35--39 (in Chinese).

田坤， 黃建平， 李振春， 李娜， 孔雪， 李慶洋. 2013. 實(shí)現(xiàn)復(fù)頻移完全匹配層吸收邊界條件的遞推卷積方法[J]. 吉林大學(xué)學(xué)報(bào)： 地球科學(xué)版， 43(3): 1022--1032.

Tian K， Huang J P， Li Z C， Li N， Kong X， Li Q Y. 2013. Recursive convolution method for implementing complex frequency-shifted PML absorbing boundary condition[J].JournalofJilinUniversity：EarthScienceEdition， 43(3): 1022--1032 (in Chinese).

張魯新， 符力耘， 裴正林. 2010. 不分裂卷積完全匹配層與旋轉(zhuǎn)交錯(cuò)網(wǎng)格有限差分在孔隙彈性介質(zhì)模擬中的應(yīng)用[J]. 地球物理學(xué)報(bào)， 53(10): 2470--2483.

Zhang L X， Fu L Y， Pei Z L. 2010. Finite difference modeling of Biot’s poroelastic equations with unsplit convolutional PML and rotated staggered grid[J].ChineseJournalofGeophysics， 53(10): 2470--2483 (in Chinese).

張顯文， 韓立國(guó)， 黃玲， 李翔， 劉前坤. 2009. 基于遞歸積分的復(fù)頻移PML彈性波方程交錯(cuò)網(wǎng)格高階差分法[J]. 地球物理學(xué)報(bào)， 52(7): 1800--1807.

Zhang X W， Han L G， Huang L， Li X， Liu Q K. 2009. A staggered-grid high-order difference method of complex frequency-shifted PML based on recursive integration for elastic wave equation[J].ChineseJournalofGeophysics， 52(7): 1800--1807 (in Chinese).

Aki K， Richards P G. 2002.QuantitativeSeismology[M]. 2nd ed. Sausalito: University Science Books: 251--254.

Bérenger J P. 1994. A perfectly matched layer for the absorption of electromagnetic waves[J].JComputPhys， 114(2): 185--200.

Boore D M. 1970. Love waves in nonuniform wave guides: Finite difference calculations[J].JGeophysRes， 75(8): 1512--1527.

Carcione J M， Kosloff D， Kosloff R. 1988. Viscoacoustic wave propagation simulation in the earth[J].Geophysics， 53(6): 769--777.

Carcione J M. 1992. Modeling anelastic singular surface waves in the earth[J].Geophysics， 57(6): 781--792.

Carcione J M. 1993. Seismic modeling in viscoelastic media[J].Geophysics， 58(1): 110--120.

Carcione J M. 2001.WaveFieldsinRealMedia：TheoryandNumericalSimulationofWavePropagationinAnisotropic，AnelasticandPorousMedia[M]. Oxford: Pergamon Press: 45--78.

Cerjan C， Kosloff D， Kosloff R， Reshef M. 1985. A nonreflecting boundary condition for discrete acoustic and elastic wave equations[J].Geophysics， 50(4): 705--708.

Clayton R， Engquist B. 1977. Absorbing boundary conditions for acoustic and elastic wave equations[J].BullSeismolSocAm， 67(6): 1529--1540.

Collino F， Tsogka C. 2001. Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media[J].Geophysics， 66(1): 294--307.

Dong H F， Ke G P， Duffaut K. 2013. Seabed shear-wave velocity estimation from dispersion curves of Love waves[J].ActaAcustunitedAc， 99(6): 890--897.

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

Futterman W I. 1962. Dispersive body waves[J].JGeophysRes， 67(13): 5279--5291.

Haskell N A. 1953. The dispersion of surface waves on multilayered media[J].BullSeismolSocAm， 43(1): 17--34.

Kanamori H， Anderson D L. 1977. Importance of physical dispersion in surface wave and free oscillation problems: Review[J].RevGeophysSpacePhys， 15(1): 105--112.

Komatitsch D， Martin R. 2007. An unsplit convolutional perfectly matched layer improved at grazing incidence for the seismic wave equation[J].Geophysics， 72(5): SM155--SM167.

Kristek J， Moczo P. 2003. Seismic-wave propagation in viscoelastic media with material discontinuities: A 3D fourth-order staggered-grid finite-difference modeling[J].BullSeismolSocAm， 93(5): 2273--2280.

Lee W B， Solomon S C. 1975. Inversion schemes for surface wave attenuation andQin the crust and the mantle[J].GeophysJRastrSoc， 43(1): 47--71.

Lee W B， Solomon S C. 1978. Simultaneous inversion of surface wave phase velocity and attenuation: Love waves in western North America[J].JGeophysRes， 83(B7): 3389--3400.

Liu H P， Anderson D L， Kanamori H. 1976. Velocity dispersion due to anelasticity: Implications for seismology and mantle composition[J].GeophysJInt， 47(1): 41--58.

Luo Y H， Xia J H， Xu Y X， Zeng C， Liu J P. 2010. Finite-difference modeling and dispersion analysis of high-frequency Love waves for near-surface applications[J].PureApplGeophys， 167(12): 1525--1536.

Martin R， Komatitsch D. 2009. An unsplit convolutional perfectly matched layer technique improved at grazing incidence for the viscoelastic wave equation[J].GeophysJInt， 179(1): 333--344.

Moczo P， Kristek J. 2005. On the rheological models used for time-domain methods of seismic wave propagation[J].GeophysResLett， 32(1): L01306.

Park C B， Miller R D， Xia J H. 1999. Multichannel analysis of surface waves (MASW)[J].Geophysics， 64(3): 800--808.

Qin Z， Lu M H， Zheng X D， Yao Y， Zhang C， Song J Y. 2009. The implementation of an improved NPML absorbing boundary condition in elastic wave modeling[J].ApplGeophys， 6(2): 113--121.

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

Safani J， O′Neill A， Matsuoka T， Sanada Y. 2005. Applications of Love wave dispersion for improved shear-wave velocity imaging[J].JEnvironEngGeophys， 10(2): 135--150.

Safani J， O′Neill A， Matsuoka T. 2006. Full SH-wavefield modeling and multiple-mode Love wave inversion[J].ExplGeophys， 37(4): 307--321.

Wang H， Tao G， Shang X F， Fang X D， Burns D R. 2013. Stability of finite difference numerical simulations of acoustic logging-while-drilling with different perfectly matched layer schemes[J].ApplGeophys， 10(4): 384--396.

Winsborrow G， Huws D G， Muyzert E. 2003. Acquisition and inversion of Love wave data to measure the lateral variability of geo-acoustic properties of marine sediments[J].JApplGeophys， 54(1): 71--84.

Winsborrow G， Huws D G， Muyzert E. 2005. The estimation of shear-wave statics usinginsitumeasurements in marine near-surface sediments[J].GeophysProsp， 53(4): 557--577.

Xia J H， Miller R D， Park C B， Tian G. 2002. DeterminingQof near-surface materials from Rayleigh waves[J].JApplGeophys， 51(2/3/4): 121--129.

Xia J H， Yin X F， Xu Y X. 2013. Feasibility of determiningQof near-surface materials from Love waves[J].JApplGeophys， 95(8): 47--52.

Zhang K， Luo Y H， Xia J H， Chen C. 2011. Pseudospectral modeling and dispersion analysis of Rayleigh waves in viscoelastic media[J].SoilDynEarthqEng， 31(10): 1332--1337.

Zhang W， Yang S. 2010. Unsplit complex frequency-shifted PML implementation using auxiliary differential equations for seismic wave modeling[J].Geophysics， 75(4): T141--T154.

Love wave modeling in viscoelastic media with unsplit CFS-CPML conditions

Xie Junfa1,2)Sun Chengyu1),Wu Dunshi1)Qiao Zhihao1)

1)SchoolofGeosciences，ChinaUniversityofPetroleum，ShandongQingdao266580，China2)ResearchInstituteofPetroleumExploration&Development-Northwest,PetroChina,Lanzhou730020,China

Numerical simulation is an important way to study the characteristics of seismic wavefield. In this paper， we established the unsplit complex frequency shifted convolutional perfectly matched layer (CFS-CPML) absorbing boundary condition， and staggered grid and high precision finite difference were used to simulate Love wavefield in viscoelastic medium. And then we analyzed the influence of the number of relaxation mechanisms on theQvalue fitting and verified the absorbing effect of CFS-CPML boundary conditions on the waves with large incident angles. The numerical results show that the requirements of the present theoretical wavefield study can be met with higher calculation efficiency if five relaxation mechanisms and finite difference with spatial four-order accuracy are selected for this method. At the same time， the phase velocity of dispersion curve gradually deviates from that of the theoretical curve with the decrease of quality factor， and the high frequency energy of each pattern also weakens. The study provided certain theoretical basis for high precision surface wave inversion method.

Love wave; forward modeling; complex frequency shifted convolutional perfectly matched layer (CFS-CPML); absorbing boundary condition; dispersion characteristic

謝俊法， 孫成禹， 伍敦仕， 喬志浩. 2016. 基于無(wú)分裂復(fù)頻移卷積完全匹配層邊界的黏彈介質(zhì)勒夫波模擬. 地震學(xué)報(bào), 38(2): 244--258. doi:10.11939/jass.2016.02.009.

Xie J F， Sun C Y， Wu D S， Qiao Z H. 2016. Love wave modeling in viscoelastic media with unsplit CFS-CPML conditions.ActaSeismologicaSinica, 38(2): 244--258. doi:10.11939/jass.2016.02.009.

國(guó)家自然科學(xué)基金(41374123)資助.

2015-08-13收到初稿， 2015-12-07決定采用修改稿.

e-mail: suncy@upc.edu.cn

10.11939/jass.2016.02.009

P315.3+1

A