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

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

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

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

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

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

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

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

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

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

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

1 方法原理

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

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

1．1 4種NAD方法

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

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

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

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

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

時間上利用截斷的泰勒公式對U和V展開可以得到：

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

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

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

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

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

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

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

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

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

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

為了比較NADM，RNAD，INAD和WNAD 4種近似解析離散化方法的數(shù)值模擬精度，設(shè)置了一維粘彈介質(zhì)模型，模型長度16km，聲波速度4 000m／ s。模型品質(zhì)因子，其中，ω表示縱波源中心圓頻率。品質(zhì)因子與吸收系數(shù)成反比，其值越大，表示能量耗散越小。當(dāng)Q趨近于無窮大時，介質(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時刻的波場快照和相對誤差曲線，相對誤差定義如下：

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

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

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

表1 采用不同計算方法和模型時位移的平均相對誤差

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

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

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

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

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

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

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

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

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

表2 不同Q及f時的信號主頻及振幅

2 二維模型試算

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

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

其中，

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

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

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

圖3 200ms時刻波場快照

圖4 z＝0處地震記錄

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

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

圖5 炮檢距 振幅曲線

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

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

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

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

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

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

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

圖7 t＝300ms波場快照

圖8 t＝400ms波場快照

圖9 t＝500ms波場快照

圖10 地震記錄

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

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

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

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

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

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

3 結(jié)束語

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

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

［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］ 董良國，馬在田，曹景忠，等．一階彈性波方程交錯網(wǎng)格高階差分解法［J］．地球物理學(xué)報，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］ 董良國，馬在田，曹景忠．一階彈性波方程交錯網(wǎng)格高階差分解法穩(wěn)定性研究［J］．地球物理學(xué)報，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］ 董良國．彈性波數(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］ 王書強(qiáng)，楊頂輝，楊寬德．彈性波方程的緊致差分方法［J］．清華大學(xué)學(xué)報（自然科學(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］ 楊頂輝，滕吉文，張中杰．三分量地震波場的近似解析離散模擬技術(shù)［J］．地球物理學(xué)報，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)的近似解析離散化方法及彈性波波場模擬［D］．北京：清華大學(xué)，2004

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

［24］ 宋國杰．三維彈性波方程的改進(jìn)近似解析離散化方法及波場模擬［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］ 宋國杰．基于WNAD方法的非一致網(wǎng)格算法及其彈性波場模擬［J］．地球物理學(xué)報，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ì)偽譜法波場模擬［J］．物理學(xué)報，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ì)中的波場特征分析［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］ 宋常瑜，裴正林．井間地震粘彈性波場特征的數(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］ 唐啟軍，韓立國，王恩利，等．基于隨機(jī)各向異性背景的粘彈性單斜介質(zhì)二維三分量正演模擬［J］．西北地震學(xué)報，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］ 單啟銅，樂友喜．PML邊界條件下二維粘彈性介質(zhì)波場模擬［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

（編輯：顧石慶）

新書介紹

《裂縫各向異性地震等效介質(zhì)理論》（精裝），王赟，劉媛媛，張美根著，2017年3月1日由科學(xué)出版社出版。本書系統(tǒng)介紹了地震各向異性領(lǐng)域近30年來發(fā)展相對成熟的4種裂縫誘導(dǎo)地震各向異性的等效介質(zhì)理論Hudson擾動理論、Thomsen等徑孔隙理論、Schoenberg線性滑移理論和Chapman多尺度裂縫等效介質(zhì)理論；并通過典型裂縫介質(zhì)建模與地震波場的有限元數(shù)值模擬、分析，對比了4種理論的異同及裂縫介質(zhì)的地震波場特征，重點討論了地震波在多組多尺度裂縫介質(zhì)中的傳播。該書可作為固體地球物理地震專業(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ù)值模擬及波場特征分析［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)主要從事地震波場特征研究和儲層預(yù)測研究。

國家重點基礎(chǔ)研究發(fā)展計劃（973計劃）（2013CB228605）和中國石油科技創(chuà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）．