近似解析離散化方法的粘彈聲波方程數(shù)值模擬及波場(chǎng)特征分析

汪 勇，段焱文，王 婷，桂志先，高 剛

（油氣資源與勘探技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室（長(zhǎng)江大學(xué)），湖北武漢430100）

近似解析離散化方法的粘彈聲波方程數(shù)值模擬及波場(chǎng)特征分析

汪 勇，段焱文，王 婷，桂志先，高 剛

（油氣資源與勘探技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室（長(zhǎng)江大學(xué)），湖北武漢430100）

地層的粘彈性對(duì)地震波產(chǎn)生的吸收和衰減規(guī)律非常復(fù)雜，所以研究地震波在粘彈性介質(zhì)中的傳播規(guī)律對(duì)地震勘探有著重要的意義。近似解析離散化方法是近年來(lái)出現(xiàn)的一種有限差分?jǐn)?shù)值模擬方法，籍此探討了幾種近似解析離散化方法的原理和過(guò)程。通過(guò)與一維粘彈聲波方程解析解的對(duì)比，認(rèn)為改進(jìn)的近似解析離散化方法（INAD）的精度、穩(wěn)定性和計(jì)算效率更具優(yōu)勢(shì)，并將其首次應(yīng)用于二維粘彈聲波方程數(shù)值模擬。通過(guò)數(shù)值模擬，重點(diǎn)分析了聲波振幅、頻率、頻帶寬度隨品質(zhì)因子及傳播距離的衰減規(guī)律，粘彈介質(zhì)中的地震波不僅因波前擴(kuò)散導(dǎo)致振幅衰減，還因吸收作用導(dǎo)致振幅和頻率同時(shí)衰減，利用其衰減特征可以指導(dǎo)含油氣儲(chǔ)層的預(yù)測(cè)。

粘彈介質(zhì)；近似解析離散方法；數(shù)值模擬；品質(zhì)因子；吸收衰減

Keywords：viscoelastic media，nearly－analytic discrete method，numerical simulation，Qvalue，absorption and attenuation

地震波場(chǎng)數(shù)值模擬是由已知的巖石結(jié)構(gòu)和物理參數(shù)建立地震地質(zhì)模型，依據(jù)不同的算法模擬地震波在地下介質(zhì)中的傳播過(guò)程，計(jì)算在地面或地下各觀測(cè)點(diǎn)的地震記錄的一種方法，是勘探地震學(xué)和地震波理論的重要基礎(chǔ)。目前，常用的數(shù)值模擬方法主要有射線追蹤法［1］和波動(dòng)方程法，其中，波動(dòng)方程法有偽譜法［23］、有限元法［4］、邊界元法［5］、譜元法［6］和有限差分法。自20世紀(jì)60年代ALTERMAN等［7］首先將有限差分方法用于地震波場(chǎng)的數(shù)值模擬以來(lái)，該方法已經(jīng)成為地震波場(chǎng)數(shù)值模擬中應(yīng)用最為廣泛的一種方法，差分格式也由早期的低階差分發(fā)展到高階差分，由常規(guī)網(wǎng)格發(fā)展到交錯(cuò)網(wǎng)格、旋轉(zhuǎn)交錯(cuò)網(wǎng)格、可變網(wǎng)格等［8－16］。近似解析離散化（Nearly－Analytic Discrete，以下簡(jiǎn)稱(chēng)NAD）［17－18］方法是近些年發(fā)展起來(lái)的一種有限差分?jǐn)?shù)值模擬方法，基本思想是在時(shí)間上采用泰勒公式展開(kāi)，在空間上利用截?cái)嗟奶├展綐?gòu)造高階插值函數(shù)來(lái)逼近空間偏導(dǎo)數(shù)，與其它有限差分法的根本區(qū)別是該方法利用空間節(jié)點(diǎn)的位移和梯度來(lái)共同逼近變量的空間偏導(dǎo)數(shù)，從而提高計(jì)算精度，在粗網(wǎng)格下可有效地壓制數(shù)值頻散，提高計(jì)算效率。YANG等［19－22］、盧明［23］、宋國(guó)杰［24－25］和TONG等［26］對(duì)該方法進(jìn)行了研究和改進(jìn)，使得該方法逐步完善。

粘彈性介質(zhì)是介于流相粘滯性介質(zhì)和固相完全彈性介質(zhì)之間的一種介質(zhì)。開(kāi)爾文粘彈性介質(zhì)模型是地震學(xué)中一種重要的模型，側(cè)重于固相的特征，對(duì)它的研究和數(shù)值模擬受到廣大學(xué)者的重視。CARCIONE等［27－28］對(duì)粘彈地層中傳播的聲波進(jìn)行了數(shù)值模擬。ROBERTSSON等［29］利用速度－應(yīng)力方程組，采用交錯(cuò)網(wǎng)格有限差分法對(duì)二維粘彈性介質(zhì)中的彈性波進(jìn)行了數(shù)值模擬研究，并且將該算法推廣到三維介質(zhì)中。杜啟振［30］利用偽譜法對(duì)各向異性粘彈介質(zhì)中的波場(chǎng)進(jìn)行了數(shù)值模擬。奚先等［31］利用波動(dòng)方程交錯(cuò)網(wǎng)格有限差分?jǐn)?shù)值模擬方法，對(duì)地震波在二維粘彈性隨機(jī)介質(zhì)中的傳播進(jìn)行了模擬并研究了其波場(chǎng)特征。宋常瑜等［32］采用交錯(cuò)網(wǎng)格高階有限差分?jǐn)?shù)值模擬方法進(jìn)行了井間地震粘彈性波場(chǎng)數(shù)值模擬。唐啟軍等［33］將Von Karman型的隨機(jī)各向同性背景引入粘彈性單斜各向異性波動(dòng)方程，并應(yīng)用交錯(cuò)網(wǎng)格技術(shù)進(jìn)行模擬。

本文以一維開(kāi)爾文粘彈模型為例，首先討論4種近似解析離散化方法的基本思想，推導(dǎo)了差分方程的建立過(guò)程，介紹了空間導(dǎo)數(shù)的求取和方程穩(wěn)定性，然后進(jìn)行數(shù)值模擬，依據(jù)數(shù)值模擬精度和計(jì)算效率，確定了適合粘彈性聲波方程數(shù)值模擬的方法。并首次將近似解析離散化方法應(yīng)用到二維粘彈聲波方程的數(shù)值模擬中，通過(guò)對(duì)波場(chǎng)快照和地面地震記錄的分析，討論了粘彈介質(zhì)中聲波的衰減規(guī)律。

1 方法原理

地震勘探的非完全彈性介質(zhì)的研究中，經(jīng)常使用一種稱(chēng)為粘彈性體或開(kāi)爾文佛各特（Kelvin－Voigt）體的模型。這種模型假定，介質(zhì)的應(yīng)力包括兩個(gè)部分，一部分是彈性應(yīng)力，另一部分是粘滯應(yīng)力，其中彈性應(yīng)力和應(yīng)變成正比，粘滯應(yīng)力和應(yīng)變的時(shí)間變化率成正比。也就是說(shuō)，應(yīng)力與應(yīng)變不再是簡(jiǎn)單的瞬時(shí)關(guān)系，而與應(yīng)變的歷史有關(guān)，其位移滿(mǎn)足的矢量波動(dòng)方程表示為：

式中：ρ為介質(zhì)密度；U＝（u，v，w）T為介質(zhì)質(zhì)點(diǎn)x，y和z方向的位移；grad表示梯度算子；div表示散度算子2為拉普拉斯算子；λ和μ為介質(zhì)的拉梅系數(shù)；λ′和μ′為描述粘滯性質(zhì)的兩個(gè)參數(shù)，稱(chēng)為粘滯系數(shù)。

1．1 4種NAD方法

近似解析離散化是一種新的數(shù)值模擬方法，下面以一維粘彈介質(zhì)聲波方程為例，討論該類(lèi)方法的基本思想及4種NAD差分方程建立過(guò)程。

由公式（1）可以得到一維粘彈介質(zhì)聲波波動(dòng)方程：

其中，u表示質(zhì)點(diǎn)x方向的位移，η＝λ＋2μ，η′＝λ′＋2μ′。

設(shè)滿(mǎn)足聲波方程的解為U，有如下定義：式中：u表示質(zhì)點(diǎn)位移大??；ux表示位移u在x方向的梯度；V表示U對(duì)時(shí)間的一階偏導(dǎo)數(shù)；P表示U對(duì)時(shí)間的二階偏導(dǎo)數(shù)。U，V和P之間滿(mǎn)足：

式中：m，k分別表示U對(duì)x和t的偏導(dǎo)階數(shù)。

時(shí)間上利用截?cái)嗟奶├展綄?duì)U和V展開(kāi)可以得到：

式中：i表示U和V的空間網(wǎng)格點(diǎn)坐標(biāo)；n和n＋1表示U和V的時(shí)間網(wǎng)格坐標(biāo)。P和P對(duì)時(shí)間的一階偏導(dǎo)數(shù)可以由波動(dòng)方程（2）轉(zhuǎn)化為U和V對(duì)空間的高階偏導(dǎo)，即：

由公式（5）至公式（8）可知，只需要知道U和V及它們的空間偏導(dǎo)數(shù)在第n時(shí)間層的值，即可計(jì)算第n＋1時(shí)間層的值，其中，i表示空間網(wǎng)格坐標(biāo)。在公式（7）和公式（8）中，有位移U和V的空間高階偏導(dǎo)數(shù)，如根據(jù)NAD類(lèi)方法的思想，這些空間導(dǎo)數(shù)是利用待求空間網(wǎng)格點(diǎn)周?chē)?7個(gè)點(diǎn)的位移及梯度共同逼近，而不是僅僅采用位移近似，這里只列舉其中兩個(gè)u對(duì)x的高階導(dǎo)數(shù)（（9）式，（10）式）進(jìn)行說(shuō)明，其余高階偏導(dǎo)數(shù)的近似函數(shù)見(jiàn)文獻(xiàn)［19］。

從公式（9）和公式（10）可以看出，要計(jì)算u對(duì)空間x的高階導(dǎo)數(shù)，不僅需要x方向左右兩個(gè)點(diǎn)的位移，還需要這兩個(gè)點(diǎn)的位移梯度ux，這正是NAD類(lèi)方法與其它方法的不同之處。由于加上了梯度值，使得近似精度更高，更適合大尺度模型下粗網(wǎng)格的數(shù)值模擬［19］。

在公式（7）和公式（8）中，有V對(duì)空間的高階導(dǎo)數(shù)，如果采用公式（11）表示的一階向后差分近似，則稱(chēng)為近似解析離散化方法（NADM）［19］。式中：m表示空間偏導(dǎo)階數(shù)。如果質(zhì)點(diǎn)速度V也采用類(lèi)似（9）式和（10）式的展開(kāi)式，則稱(chēng)為改進(jìn)的近似解析離散化（Improving Nearly－Analytic Discrete，簡(jiǎn)稱(chēng)INAD）方法［22］。如果將公式（5）改為公式（12），則稱(chēng)為精化的近似解析離散化方法（Refined Nearly－Analytic Discrete，簡(jiǎn)稱(chēng)RNAD）［24］，僅僅是將公式（5）中上一時(shí)刻的速度場(chǎng)替換為當(dāng)前速度場(chǎng)，用來(lái)計(jì)算新時(shí)間層的位移。

如果計(jì)算位移場(chǎng)U時(shí)的速度項(xiàng)V由n和n＋1時(shí)刻速度場(chǎng)的加權(quán)組合來(lái)代替，則稱(chēng)為加權(quán)近似解析離散化方法（Weighted Nearly Analytic Discrete，簡(jiǎn)稱(chēng)WNAD）［25］，即采用公式（13）對(duì)位移場(chǎng)進(jìn)行時(shí)間層推進(jìn)。

式中：φ為加權(quán)系數(shù)，取值范圍為0～1。

1．2 算法穩(wěn)定性

穩(wěn)定性條件是有限差分?jǐn)?shù)值模擬中一個(gè)非常重要的問(wèn)題，如果不滿(mǎn)足穩(wěn)定性條件，則前面各層的舍入誤差將會(huì)影響到后面各層的值，如果誤差的影響越來(lái)越大，會(huì)使得有限差分的結(jié)果完全錯(cuò)誤。文獻(xiàn)［19，22，24，25］利用Fourier方法研究了以上4種方法的穩(wěn)定性，并給出了穩(wěn)定性條件。假設(shè)Δt為時(shí)間步長(zhǎng)，Δx為空間網(wǎng)格間距，最大的速度值為vmax，穩(wěn)定性條件為：

穩(wěn)定性條件方面，NADM方法最為嚴(yán)格，空間步長(zhǎng)相同時(shí)，所需的時(shí)間網(wǎng)格最小。INAD方法最為寬松，所需要滿(mǎn)足穩(wěn)定性條件的時(shí)間網(wǎng)格最大，更有利于進(jìn)行長(zhǎng)采樣時(shí)間的數(shù)值模擬。

1．3 算法精度和效率分析

為了比較NADM，RNAD，INAD和WNAD 4種近似解析離散化方法的數(shù)值模擬精度，設(shè)置了一維粘彈介質(zhì)模型，模型長(zhǎng)度16km，聲波速度4 000m／ s。模型品質(zhì)因子，其中，ω表示縱波源中心圓頻率。品質(zhì)因子與吸收系數(shù)成反比，其值越大，表示能量耗散越小。當(dāng)Q趨近于無(wú)窮大時(shí)，介質(zhì)變?yōu)橥耆珡椥越橘|(zhì)，粘彈聲波方程則變?yōu)榱送耆珡椥月暡ǚ匠獭?/p>

公式（2）表示的一維粘彈介質(zhì)聲波方程，其解析解可以表示為：式中：A為振幅；波數(shù)k和吸收系數(shù)α都與頻率ω有關(guān)。

數(shù)值模擬中，取dt＝0．8ms，dx＝20m。在模型的8km處加載公式（15）表示的震源，可得圖1所示的2s時(shí)刻的波場(chǎng)快照和相對(duì)誤差曲線，相對(duì)誤差定義如下：

式中：u（tn，xi）表示精確解析解表示數(shù)值解。采用WNAD方法時(shí)，公式（13）中的加權(quán)系數(shù)φ取0．5［25］。

圖1 一維粘彈聲波方程數(shù)值模擬波場(chǎng)快照（a）和相對(duì)誤差（b）

從圖1可以看出：①I(mǎi)NAD和WNAD兩種方法數(shù)值模擬的結(jié)果與解析解的波剖面基本重合，模擬精度較高，而NADM和RNAD的相對(duì)誤差較大；②數(shù)值模擬結(jié)果顯示，隨著傳播距離的增加，振幅隨之發(fā)生衰減，而波長(zhǎng)隨之增加。采用同樣的方法，分析不同模型參數(shù)下的相對(duì)誤差，結(jié)果如表1所示。

表1 采用不同計(jì)算方法和模型時(shí)位移的平均相對(duì)誤差

從表1可以看出，4種方法的相對(duì)誤差都隨時(shí)間步長(zhǎng)或空間步長(zhǎng)的增加而增大。整體而言，INAD結(jié)果的相對(duì)誤差最小，在空間步長(zhǎng)為50m時(shí)，平均相對(duì)誤差為10．31%，說(shuō)明了該方法的可靠性和實(shí)用性。

下面從所需存儲(chǔ)空間和計(jì)算時(shí)間兩方面分析以上4種方法的計(jì)算效率。由于RNAD，WNAD和INAD 3種方法在速度場(chǎng)V的時(shí)間層推進(jìn)上沒(méi)有區(qū)別，只是在計(jì)算位移場(chǎng)U 時(shí)采用的速度場(chǎng)不同，RNAD用的是n＋1時(shí)刻的速度場(chǎng)，INAD用的是n時(shí)刻的速度場(chǎng)，WNAD用的是加權(quán)的速度場(chǎng)，所以這3種方法存儲(chǔ)數(shù)組的數(shù)量沒(méi)有差別。而NADM方法由于不需要按公式（6）進(jìn)行速度場(chǎng)時(shí)間層推進(jìn)，而是用位移場(chǎng)的向后差分進(jìn)行計(jì)算，所以存儲(chǔ)數(shù)組個(gè)數(shù)減少了2個(gè)。從計(jì)算時(shí)間來(lái)看，由于INAD具有最大的庫(kù)郎數(shù)，在取同樣空間步長(zhǎng)時(shí)，它可采用較大的時(shí)間步長(zhǎng)，所以在同樣的采樣時(shí)間內(nèi)，其計(jì)算次數(shù)和時(shí)間最少。例如計(jì)算上述模型1時(shí)，NADM，RNAD，WNAD和INAD可取的時(shí)間步長(zhǎng)分別為0．8，1．2，1．7，2．2ms，它們的時(shí)間推進(jìn)層數(shù)分別約為2 500，1 667，1 176和909次，INAD計(jì)算時(shí)間約是NADM，RNAD和WNAD的36%，54%和77%。

綜合來(lái)看，對(duì)粘彈聲波方程進(jìn)行數(shù)值模擬時(shí)，INAD方法具有較高的精度和計(jì)算效率，所以本文利用該方法進(jìn)行數(shù)值模擬，并進(jìn)一步分析一維粘彈介質(zhì)中聲波的傳播規(guī)律。模型長(zhǎng)度設(shè)置為4km，聲波速度4 000m／s，空間步長(zhǎng)dx＝10m，時(shí)間步長(zhǎng)dt＝0．5ms。數(shù)值模擬時(shí)，改變品質(zhì)因子和信號(hào)頻率：①Q(mào)→∞，f＝50Hz；②Q＝100，f＝50Hz；③Q＝100，f＝20Hz；④Q＝10，f＝50Hz。在x＝2km處加載雷克子波震源信號(hào)，得到2 000，2 500，3 000，3 500m處質(zhì)點(diǎn)的振動(dòng)曲線和振幅譜如圖2所示，振動(dòng)曲線的主頻及振幅見(jiàn)表2。

圖2 一維粘彈聲波方程數(shù)值模擬振動(dòng)曲線（a，b，c，d）及其對(duì)應(yīng)的頻譜（e，f，g，h）

一維模型中傳播的聲波為平面波，不存在波前擴(kuò)散，振幅和頻率的衰減都是由介質(zhì)的粘彈性導(dǎo)致的，分析圖2和表2，結(jié)論如下。

1）當(dāng)Q→∞時(shí)，隨傳播距離和時(shí)間的增加，振幅和頻率沒(méi)有衰減，粘彈聲波方程等價(jià)于完全彈性聲波方程，同時(shí)也證明了INAD方法模擬結(jié)果的可靠性。

2）比較②和④說(shuō)明，當(dāng)震源主頻不變時(shí)，Q的大小決定了振幅和頻率的衰減速度，當(dāng)Q＝10時(shí)，振幅發(fā)生最大的衰減，衰減至0．004；頻率也降低得最快，降至8Hz。

3）比較②和③可以看出，當(dāng)Q一定時(shí)，衰減系數(shù)近似地與頻率的平方成正比，高頻信號(hào)衰減得快，低頻衰減得慢。②中的頻率降低了26Hz，振幅降低了90%，而③中的頻率降低了8Hz，振幅降低了74%。所以在地震勘探中，地震波隨著傳播距離的增大，高頻成分很快被衰減，深部地震信號(hào)只保留較低的頻率成分。

4）同公式（14）表示的含義相同，振幅隨傳播距離呈指數(shù)衰減，開(kāi)始衰減速度快，隨傳播距離的增大，衰減速率降低。頻率的衰減也是如此。如模型④中的最初500m，頻率降低了36Hz，振幅降低了98%，而最后500m，頻率只降低了2Hz，振幅降低50%。

表2 不同Q及f時(shí)的信號(hào)主頻及振幅

2 二維模型試算

二維粘彈介質(zhì)聲波波動(dòng)方程為：

與一維情況相同，二維粘彈介質(zhì)聲波方程INAD方法的差分公式依然為公式（5）和公式（6），只是公式中P和有所變化，如下式所示：

其中，

且品質(zhì)因子Q與λ，μ，λ′和μ′有關(guān)，其大小由1．3節(jié)中公式給出。

2．1 均勻粘彈聲波方程模擬

設(shè)置模型長(zhǎng)度和深度均為2 000m，vP＝4 000m／s。為了比較，計(jì)算Q＝20，50，100和趨近于∞（即完全彈性聲波模擬）4種情況。時(shí)間步長(zhǎng)Δt＝0．5ms，空間網(wǎng)格Δx＝Δz＝10m，模型中間（x＝1 000m，z＝1 000m）位置激發(fā)主頻50Hz雷克子波，獲得的200ms波場(chǎng)快照如圖3所示，圖4為地表接收到的地震記錄，記錄長(zhǎng)度500ms。

圖3 200ms時(shí)刻波場(chǎng)快照

圖4 z＝0處地震記錄

從圖3和圖4所示的波場(chǎng)快照和地震記錄可以看出：①模擬波場(chǎng)非常清晰，沒(méi)有數(shù)值頻散；②同一時(shí)刻，粘彈聲波模擬和完全彈性（Q→∞）聲波模擬的波前位置相同，證實(shí)了方法的可靠性；③隨著Q的減小，地震波波長(zhǎng)增加，波數(shù)降低；頻率降低，同相軸變波的振幅歸一化到0～1，得到炮檢距 振幅曲線（圖5）。從圖5可以看出：①當(dāng)采用完全彈性聲波模擬時(shí)（Q→∞），與一維情況不同，振幅隨地震波的傳播發(fā)生波前擴(kuò)散，振幅與傳播距離呈反比關(guān)系；②當(dāng)采用粘彈聲波模擬時(shí)，除存在波前擴(kuò)散外，還存在粘彈介質(zhì)對(duì)地震波的吸收衰減，使得衰減加??；③衰減速度由Q決定，Q越小，衰減越快。寬，子波延續(xù)時(shí)間增加；④粘彈性介質(zhì)中，地震波隨傳播距離的增加，振幅呈指數(shù)衰減，地震記錄中遠(yuǎn)道的振幅要明顯小于近道。

為了分析地震波衰減規(guī)律，將圖4中的4個(gè)反射

圖5 炮檢距 振幅曲線

對(duì)圖4所示的地震記錄，提取x＝1 000m位置處質(zhì)點(diǎn)的單道記錄并進(jìn)行頻譜分析，得到不同Q情況下的振動(dòng)曲線的振幅譜（圖6）。圖6顯示了不同Q時(shí)頻率衰減大小，完全彈性聲波模擬時(shí)記錄信號(hào)的主頻依然為50Hz，頻帶范圍為32～70Hz。隨著Q的減小，粘彈聲波信號(hào)的高頻成分減少得越多，Q＝20時(shí)的頻率降低至14Hz，頻帶范圍減小為8～19Hz，這與實(shí)際地震記錄中的頻譜特征一致。

圖6 振幅譜曲線（x＝1 000m，z＝0）

2．2 多層水平層狀模型模擬

利用INAD方法對(duì)4層層狀模型分別進(jìn)行粘彈和完全彈性聲波方程模擬試算，模型長(zhǎng)度和深度均為2 000m，各層巖石物理參數(shù)如表3所示。

表3 4層水平層狀介質(zhì)模型參數(shù)

為了比較粘彈性和完全彈性介質(zhì)中地震波場(chǎng)的差異，第2層和第3層的Q分別取50，100。時(shí)間步長(zhǎng)Δt＝0．5ms，空間網(wǎng)格Δx＝Δz＝10m，模型地表中間（x＝1 000m，z＝0）位置激發(fā)主頻為40Hz的雷克子波，獲得不同時(shí)刻的波場(chǎng)快照（圖7，圖8，圖9）和地表接收的單炮記錄（圖10），記錄長(zhǎng)度600ms。在數(shù)值模擬中，采用PML吸收邊界處理，算法見(jiàn)文獻(xiàn)［33］和文獻(xiàn)［34］，這里不再贅述。

圖7至圖10中，P1，P2和P3代表來(lái)自3個(gè)反射界面的反射P波。圖7中，因?yàn)榈?層設(shè)置為完全彈性介質(zhì)，所以粘彈聲波和彈性聲波模擬快照中的P1相同。而透射進(jìn)入第2層、第3層地層的P波則出現(xiàn)了差異，正是由于粘彈介質(zhì)對(duì)地震波的吸收衰減，導(dǎo)致了粘彈聲波模擬結(jié)果中的P2和P3反射波振幅明顯衰減，同相軸變寬，子波延續(xù)時(shí)間增加，波長(zhǎng)增大。

圖7 t＝300ms波場(chǎng)快照

圖8 t＝400ms波場(chǎng)快照

圖9 t＝500ms波場(chǎng)快照

圖10 地震記錄

圖11 x＝1 500m處單道地震記錄

圖11是從上述兩次模擬地面記錄中抽取的單道地震記錄（x＝1 500m），由于直達(dá)波能量太強(qiáng)，所以顯示時(shí)間從300ms開(kāi)始。波前擴(kuò)散導(dǎo)致彈性聲波記錄的振幅隨傳播時(shí)間和距離逐步減弱，而粘彈聲波記錄的振幅衰減速度增大，這與波場(chǎng)快照中的規(guī)律一致。圖11顯示了第1層反射波（350ms左右）振幅，可以看出，彈性聲波的振幅要略小于粘彈聲波的振幅，這是由于第2層變?yōu)檎硰椊橘|(zhì)后，導(dǎo)致P波的垂直反射系數(shù)增大。為了進(jìn)一步分析兩次模擬記錄在頻譜上的差別，利用廣義S變換，對(duì)這兩道地震記錄進(jìn)行時(shí)頻分析，結(jié)果如圖12和圖13所示。

從時(shí)頻分析結(jié)果中提取3層反射波主頻和振幅，結(jié)果如表4所示。表4中數(shù)據(jù)表明：①當(dāng)采用彈性聲波模擬時(shí)，子波信號(hào)的主頻和頻帶范圍基本不變，實(shí)際上這并不能準(zhǔn)確地描述地震波在地下的傳播規(guī)律。在傳播過(guò)程中發(fā)生的振幅衰減是由波前擴(kuò)散導(dǎo)致的，在2 000m厚度模型中，振幅衰減了約50%。②采用粘彈聲波模擬能更加精確地描述地震波在粘彈介質(zhì)中的傳播特征，子波主頻和頻帶寬度均隨傳播深度的增加而減小，衰減速度與介質(zhì)的品質(zhì)因子和子波頻率有關(guān)。振幅除了波前擴(kuò)散外，衰減會(huì)變得更劇烈，本模型中，振幅衰減了約98%。③雖然子波頻率降低和頻帶變窄會(huì)導(dǎo)致地震資料縱向分辨率降低，但這些頻率衰減特征正好可以指導(dǎo)含油氣儲(chǔ)層的預(yù)測(cè)工作（因?yàn)閹r石含流體后，會(huì)使其粘彈性質(zhì)更加明顯）。

圖12 x＝1 500m處單道地震記錄時(shí)頻分析結(jié)果（完全彈性聲波模擬）

圖13 x＝1 500m處單道地震記錄時(shí)頻分析結(jié)果（粘彈聲波模擬）

表4 彈性和粘彈性聲波模擬的3層反射波主頻及振幅

3 結(jié)束語(yǔ)

粘彈介質(zhì)地震波場(chǎng)的數(shù)值模擬研究對(duì)地震勘探具有重要的意義，有助于人們更清楚地認(rèn)識(shí)地震波在巖石和油氣儲(chǔ)集區(qū)的傳播規(guī)律，具有較強(qiáng)的理論和現(xiàn)實(shí)意義。本文將一種新的數(shù)值模擬方法應(yīng)用到粘彈介質(zhì)聲波方程的數(shù)值模擬中，研究結(jié)果表明：①改進(jìn)的近似解析離散化（INAD）有限差分算法在粗網(wǎng)格下具有較高的精度，更有利于大尺度模型的數(shù)值模擬，為粘彈介質(zhì)的地震波場(chǎng)數(shù)值模擬提供了新的思路和方法；②采用完全彈性聲波模擬時(shí)，振幅只隨波前擴(kuò)散發(fā)生衰減，而子波信號(hào)的頻率和頻帶范圍基本不變；③采用粘彈聲波模擬時(shí)，振幅在波前擴(kuò)散的基礎(chǔ)上，衰減幅度增加；同時(shí)，子波主頻和頻帶寬度隨傳播深度的增加而減小，衰減速度與介質(zhì)的品質(zhì)因子和子波頻率有關(guān)。

由于近似解析離散化方法在空間上可求得位移的最高5階空間導(dǎo)數(shù)，且粘彈介質(zhì)波動(dòng)方程中存在速度項(xiàng)，使得本文中采用的INAD方法的差分精度為時(shí)間二階和空間四階，在下一步研究中，需要進(jìn)一步研究時(shí)間差分精度更高的近似解析離散化方法，提高粘彈介質(zhì)的模擬精度。其次，可將該類(lèi)方法在其它粘彈介質(zhì)模型的全波場(chǎng)數(shù)值模擬中進(jìn)一步推廣，擴(kuò)大該類(lèi)方法的適用范圍。

［1］ CERVENY V．Seismic rays and ray intensities in inhomogeneous anisotropic media［J］．Geophysical Journal International，1972，29（1）：1－33

［2］ KOSLOFF D D，BAYSAL E．Forward modeling by a Fourier method［J］．Geophysics，1982，47（10）：1402－1412

［3］ RESHEF M，KOSLOFF D D，EDWARDS M，et al．Three dimensional elastic modeling by the Fourier method［J］．Geophysics，1988，53（9）：1184－1193

［4］ DRAKE L A．Rayleigh waves at a continental boundary by the finite element method［J］．Bulletin of the Seismological Society of America，1972，62（5）：1259－1268

［5］ BOUCHON M，COUTANT O．Calculation of synthetic seismograms in a laterally varying medium by the boundary element－discrete wavenumber method［J］．Bulletin of the Seismological Society of America，1994，84（6）：1869－1881

［6］ KOMATISSCH D，BARNES C，TROMP J．Simulation of anisotropic wave propagation based upon a spectral element method［J］．Geophysics，2000，65（4）：1251－1260

［7］ ALTERMAN Z，KARAL F C．Propagation of seismic wave in layered media by finite difference methods［J］．Bulletin of the Seismological Society of America，1968，58（1）：367－398

［8］ KELLY K R，WARD R W，TREITEL S，et al．Synthetic seismograms：a finite－difference approach［J］．Geophysics，1976，41（1）：2－27

［9］ VIRIEUX J．P－SV wave propagation in heterogeneous media：velocity－stress finite－difference method［J］．Geophysics，1986，51（4）：889－901

［10］ GRAVES R W．Simulation seismic wave propagation in 3Delastic media using staggered－grid finite differences［J］．Bulletin of the Seismological Society of America，1996，86（4）：1091－1106

［11］ 董良國(guó)，馬在田，曹景忠，等．一階彈性波方程交錯(cuò)網(wǎng)格高階差分解法［J］．地球物理學(xué)報(bào)，2000，43（3）：411－419 DONG L G，MA Z T，CAO J Z，et al．The staggeredgrid high－order difference method of first－order elastic equation［J］．Chinese Journal of Geophysics，2000，43（3）：411－419

［12］ 董良國(guó)，馬在田，曹景忠．一階彈性波方程交錯(cuò)網(wǎng)格高階差分解法穩(wěn)定性研究［J］．地球物理學(xué)報(bào)，2000，43（6）：856－864 DONG L G，MA Z T，CAO J Z．The stability study of the staggered－grid high－order difference method of first－order elastic equation，［J］．Chinese Journal of Geophysics，2000，43（6）：856－864

［13］ 董良國(guó)．彈性波數(shù)值模擬中的吸收邊界條件［J］．石油地球物理勘探，1999，34（1）：46－56

DONG L G．Absorptive boundary condition in elasticwave numerical modeling［J］．Oil Geophysical Prospecting，1999，34（1）：46－56

［14］ 王書(shū)強(qiáng)，楊頂輝，楊寬德．彈性波方程的緊致差分方法［J］．清華大學(xué)學(xué)報(bào)（自然科學(xué)版），2002，42（8）：1128－1131

WANG S Q，YANG D H，YANG K D．Compact finite difference scheme for elastic equations［J］．Journal of Tsinghua University（Science and Technology），2002，42（8）：1128－1131

［15］ SAENGER E H，GOLD N，SHAPIRO S A．Modeling the propagation of elastic waves using a modified finitedifference grid［J］．Wave Motion，2000，31（1）：77－92

［16］ SAENGER E H，BOHLEN T．Finite－difference modeling of viscoelastic and anisotropic wave propagation using the rotated staggered grid［J］．Geophysics，2004，69（2）：583－591

［17］ KONDOH Y，HOSAKA Y，ISHII K．Kernel optimum nearly－analytical discrimination algorithm applied to parabolic and hyperbolic equations［J］．Computers ＆Mathematics with Applications，1994，27（3）：59－90

［18］ 楊頂輝，滕吉文，張中杰．三分量地震波場(chǎng)的近似解析離散模擬技術(shù)［J］．地球物理學(xué)報(bào)，1996，39（增刊1）：283－291

YANG D H，TENG G W，ZHANG Z J．Nearly－analytical discretization modeling technique of 3－component seismic wave－fields［J］．Chinese Journal of Geophysics，1996，39（S1）：283－291

［19］ YANG D H，TENG J W，ZHANG Z J，et al．A nearlyanalytic discrete method for acoustic and elastic wave equation［J］．Bulletin of the Seismological Society of A－merica，2003，93（2）：882－890

［20］ YANG D H，LU M，WU R S，et al．An optimal nearlyanalytic discrete method for 2Dacoustic and elastic wave equations［J］．Bulletin of the Seismological Society of America，2004，94（5）：1982－1991

［21］ YANG D H，SONG G J，LU M．Optimally accurate nearly－analytic discrete scheme for wave－filed simulation in 3Danisotropic media［J］．Bulletin of the Seismological Society of America，2007，97（5）：1557－1569

［22］ YANG D H，SONG G J，CHEN S，et al．An improved nearly analytical discrete method：an efficient tool to simulate the seismic response of 2－D porous structures［J］．Journal of Geophysics ＆Engineering，2007，4（1）：40－52

［23］ 盧明．改進(jìn)的近似解析離散化方法及彈性波波場(chǎng)模擬［D］．北京：清華大學(xué)，2004

LU M．The optimum nearly analytic discrete methodand elastic wave numerical simulations［D］．Beijing：Tsinghua University，2004

［24］ 宋國(guó)杰．三維彈性波方程的改進(jìn)近似解析離散化方法及波場(chǎng)模擬［D］．北京：清華大學(xué)，2008

SONG G J．The improved nearly analytic discrete method for 3Delastic equations and its numerical simulation［D］．Beijing：Tsinghua University，2008

［25］ 宋國(guó)杰．基于WNAD方法的非一致網(wǎng)格算法及其彈性波場(chǎng)模擬［J］．地球物理學(xué)報(bào)，2010，53（8）：1985－1992

SONG G J．Non－uniform grid algorithm based on the WNAD method and elastic wave－field simulations［J］．Chinese Journal of Geophysics，2010，53（8）：1985－1992

［26］ TONG P，YANG D H，HUA B L．High accuracy wave simulation－Revised derivation，numerical analysis and testing of a nearly analytic integration discrete method for solving acoustic equation［J］．International Journal of Solids ＆Structures，2011，48（1）：56－70

［27］ CARCIONE J M，KOSLOFF D，KOSLOFF R．Wave propagation simulation in a linear viscoacoustic medium［J］．Geophysical Journal International，1988，93（2）：393－401

［28］ CARCIONE J M，KOSLOFF D，KOSLOFF R．Viscoacoustic wave propagation simulation in the Earth［J］．Geophysics，1988，53（6）：769－777

［29］ ROBERTSSON J O A，BLANCH J O．Viscoelastic finite－difference modeling［J］．Geophysics，1994，59（9）：1444－1456

［30］ 杜啟振．各向異性黏彈性介質(zhì)偽譜法波場(chǎng)模擬［J］．物理學(xué)報(bào)，2004，53（12）：4428－4434

DU Q Z．Wavefield forward modeling with the pseudo－spectral method in viscoelastic and azimuthally anisotropic media［J］．Acta Physica Sinica，2004，53（12）：4428－4434

［31］ 奚先，姚姚．二維粘彈性隨機(jī)介質(zhì)中的波場(chǎng)特征分析［J］．地球物理學(xué)進(jìn)展，2004，19（3）：608－615

XI X，YAO Y．The analysis of the wave field characteristics in 2－D viscoelastic random medium［J］．Progress in Geophysics，2004，19（3）：608－615

［32］ 宋常瑜，裴正林．井間地震粘彈性波場(chǎng)特征的數(shù)值模擬研究［J］．石油物探，2006，45（5）：508－513

SONG C Y，PEI Z L．Numerical simulation of viscoelastic wevefield exploration［J］．Geophysical Prospecting for Petroleum，2006，45（5）：508－513

［33］ 唐啟軍，韓立國(guó)，王恩利，等．基于隨機(jī)各向異性背景的粘彈性單斜介質(zhì)二維三分量正演模擬［J］．西北地震學(xué)報(bào)，2009，30（1）：35－39

TANG Q J，HAN L G，WANG E L．2－D three－component seismic modeling for viscoelastic monoclinic media based on background of random isotropic media［J］．Northwestern Seismological Journal，2009，30（1）：35－39

［34］ 單啟銅，樂(lè)友喜．PML邊界條件下二維粘彈性介質(zhì)波場(chǎng)模擬［J］．石油物探，2007，46（2）：126－130

SHAN Q T，YUE Y X．Wavefield simulation of 2－D viscoelastic medium in perfectly matched layer boundary［J］．Geophysical Prospecting for Petroleum，2007，46（2）：126－130

（編輯：顧石慶）

新書(shū)介紹

《裂縫各向異性地震等效介質(zhì)理論》（精裝），王赟，劉媛媛，張美根著，2017年3月1日由科學(xué)出版社出版。本書(shū)系統(tǒng)介紹了地震各向異性領(lǐng)域近30年來(lái)發(fā)展相對(duì)成熟的4種裂縫誘導(dǎo)地震各向異性的等效介質(zhì)理論Hudson擾動(dòng)理論、Thomsen等徑孔隙理論、Schoenberg線性滑移理論和Chapman多尺度裂縫等效介質(zhì)理論；并通過(guò)典型裂縫介質(zhì)建模與地震波場(chǎng)的有限元數(shù)值模擬、分析，對(duì)比了4種理論的異同及裂縫介質(zhì)的地震波場(chǎng)特征，重點(diǎn)討論了地震波在多組多尺度裂縫介質(zhì)中的傳播。該書(shū)可作為固體地球物理地震專(zhuān)業(yè)本科生、研究生的教材，也可供工程技術(shù)人員參考。

Numerical simulation and the wave field characteristics analysis of viscoelastic acoustic wave equation based on the nearly－analytic discrete method

WANG Yong，DUAN Yanwen，WANG Ting，GUI Zhixian，GAO Gang
（Key Laboratory of Exploration Technologies for Oil and Gas Resources，Ministry of Education（Yangtze University），Wuhan 430100，China）

Viscoelasticity of the media has a very complicated relationship with the absorption and attenuation of seismic wave energy．Therefore，the research on the propagation law of seismic waves in the stratigraphic equivalent media carries great significance to the seismic exploration．Nearly－analytic discrete method is a finite－difference numerical simulation method developed in recent years．In this paper，the author firstly explored the principles and process of four different nearly－analytic discrete methods．In comparison with the one－dimensional analytical solutions，the improved nearly－analytic discrete method（INAD）was proved to excel in precision，stability and computational efficiency．Then INAD was applied to the numerical simulation of 2Dviscoelastic acoustic wave equation for the first time．With the simulated viscoelastic wave field，the author focused on analyzing the propagation law of the amplitude，frequency and bandwidth of the acoustic wave with the change of Qvalues and propagation distance The amplitude of the seismic wave in the viscoelastic medium reduces not only due to the wavefront diffusion，but also the absorption which attenuates the frequency at the same time．These attenuation characteristics can guide the prediction of the oil and gas reservoir．The amplitude of the seismic wave in the viscoelastic medium reduces not only due to the wavefront diffusion，but also the absorption which attenuates the frequency at the same time．These attenuation characteristics can guide the prediction of the oil and gas reservoir．

P631

A

1000－1441（2017）03－0362－11

10．3969／j．issn．1000－1441．2017．03．006

汪勇，段焱文，王婷，等．近似解析離散化方法的粘彈聲波方程數(shù)值模擬及波場(chǎng)特征分析［J］．石油物探，2017，56（3）：362－372

WANG Yong，DUAN Yanwen，WANG Ting，et al．Numerical simulation and the wave field characteristics analysis of viscoelastic acoustic wave equation based on the nearly－analytic discrete method［J］．Geophysical Prospecting for Petroleum，2017，56（3）：362－372

2016－03－24；改回日期：2016－10－18。

汪勇（1979—），男，副教授，博士，現(xiàn)主要從事地震波場(chǎng)特征研究和儲(chǔ)層預(yù)測(cè)研究。

國(guó)家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃（973計(jì)劃）（2013CB228605）和中國(guó)石油科技創(chuàng)新基金項(xiàng)目（2015D－5006－0301）聯(lián)合資助。

This research is financially supported by the National Key Basic Research Project of China（973Program）（Grant No．2013CB228605）and the Science and Technology Innovation Foundation of CNPC（Grant No．2015D－5006－0301）．