何兵紅 吳國忱

（中國山東青島 266580 中國石油大學(xué)（華東）地球科學(xué)與技術(shù)學(xué)院）

引言

油氣儲層與地下介質(zhì)中的孔隙裂縫、地層的非均質(zhì)性和孔隙流體有關(guān)．基于完全彈性假設(shè)的彈性波理論不能滿足對油氣儲層中地震波能量損失、速度頻散以及波形扭曲等波場特征的表達(dá)．為了更精確地描述實(shí)際油氣儲層中地震波的傳播規(guī)律，需要建立與油氣特征相應(yīng)的介質(zhì)模型．介質(zhì)模型與地震波的傳播特征密切相關(guān)（牛濱華，孫春巖，2004）．基于黏彈性介質(zhì)理論的地震波響應(yīng)中包含了彈性響應(yīng)和黏滯性響應(yīng)，能夠合理地表達(dá)油氣儲層所表現(xiàn)出的固相彈性介質(zhì)和流相黏性介質(zhì)特征，從而成為研究油氣儲層中地震波傳播特征的主要介質(zhì)模型．黏彈性介質(zhì)中的黏滯性特征通常用地層品質(zhì)因子Q進(jìn)行衡量．實(shí)驗(yàn)測試表明，地層品質(zhì)因子Q在地震頻帶內(nèi)幾乎與頻率無關(guān)，但其衰減與頻率密切相關(guān)，故頻率域更有利于引入黏彈性特征的描述．時間域可利用彈簧與阻尼組合建立基于流變學(xué)的黏彈性模型，將黏滯性特征引入波動方程．基本流變體模型包含麥克斯韋模型和開爾文－佛格特模型，其中麥克斯韋模型中Q值隨頻率的增加呈線性關(guān)系；而開爾文－佛格特模型中Q值則與頻率成反比關(guān)系，并且在地震頻帶內(nèi)Q值隨頻率變化較為劇烈．這些基本流變體模型無法滿足對黏彈性介質(zhì)中常Q模型特征的精確模擬，其根本原因在于只包含單個松弛時間特征；而廣義流變體模型將多個基本流變體模型進(jìn)行組合，使其更能適合對常Q黏彈性模型的描述（Carcioneetal，1988）．流變體模型中采用應(yīng)力松弛時間與應(yīng)變延遲時間表示介質(zhì)的黏滯性強(qiáng)度，需要建立地層品質(zhì)因子Q與應(yīng)力松弛時間和應(yīng)變延遲時間之間的轉(zhuǎn)化關(guān)系．常規(guī)的τ值法由于在應(yīng)力松弛時間與應(yīng)變延遲時間相等的假設(shè)條件下確定常Q模型參數(shù)，故得到的常Q模型誤差較大，特別是對于低Q值的情況，地震波衰減程度存在很大的誤差．

時間域黏彈性介質(zhì)數(shù)值模擬方面，國內(nèi)外研究均已成熟．Day和Minster（1984）通過Padé有理近似將時間域的褶積關(guān)系轉(zhuǎn)化為微分關(guān)系，進(jìn)而推動了時間域線性黏彈性體的數(shù)值模擬．Emmerich和Korn（1987）基于廣義麥克斯韋模型確定了具有明確物理含義的有理近似式系數(shù)，解決了Padé有理近似只能適用于弱衰減短波長的問題．Carcione等（1988）和Carcione（1993）在廣義標(biāo)準(zhǔn)線性黏彈體中引入記憶變量的概念，并利用偽譜法計(jì)算空間導(dǎo)數(shù)，實(shí)現(xiàn)了二維和三維黏彈性介質(zhì)的數(shù)值模擬．偽譜法能夠減弱頻散特征的影響，但是離散傅里葉變換的計(jì)算降低了計(jì)算效率．Robertsson等（1994）利用時間－空間域有限差分算法，實(shí)現(xiàn)了基于廣義標(biāo)準(zhǔn)線性固體模型的黏彈性介質(zhì)的數(shù)值模擬．?tekl和Pratt（1998）在頻率域利用旋轉(zhuǎn)有限差分算法開展了黏彈性數(shù)值模擬．Kristek和Moczo（2003）基于廣義麥克斯韋模型重新定義了與黏滯性系數(shù)無關(guān)的黏彈體方程，實(shí)現(xiàn)了三維黏彈性介質(zhì)交錯網(wǎng)格有限差分的正演模擬．楊文采（1986）基于諧Q模型，采用一維頻率域遞推實(shí)現(xiàn)了黏彈性介質(zhì)中反射地震道的合成．畢玉英和楊寶?。?995）推導(dǎo)了頻率域波場遞推關(guān)系，并形成了含有AVO特征的復(fù)反射系數(shù)，利用線源代替點(diǎn)源從一維計(jì)算推廣到二維計(jì)算．崔建軍和何繼善（2001）從二維黏彈性波動方程出發(fā)，在f-k域內(nèi)實(shí)現(xiàn)了基于開爾文介質(zhì)的黏彈性波動方程的正演和偏移．杜啟振等（2006）研究了線性黏彈性各向異性介質(zhì)中速度頻散和衰減特征．郭少華（2004）研究了基于開爾文－佛格特模型的黏彈性波動力學(xué)特征．苑春方等（2005）利用疊加原理推導(dǎo)了在任意震源條件下，開爾文－佛格特模型均勻黏彈性介質(zhì)中三階波動方程地震波傳播的平面波解析解．單啟銅（2007）基于開爾文模型實(shí)現(xiàn)了黏彈性波動方程正演模擬及參數(shù)反演．孫成禹等（2010）分析了矩形網(wǎng)格下不同黏彈模型波動方程有限差分解的穩(wěn)定性，導(dǎo)出了開爾文－佛格特黏彈模型和麥克斯韋黏彈模型在任意空間差分精度下穩(wěn)定性條件的表達(dá)式．

基于全波波動方程的黏彈性介質(zhì)地震波不僅包含反射波，還包含折射波、潛水波以及多次波，能夠有效地模擬復(fù)雜構(gòu)造中地震波的傳播特征．但在目前地震勘探中，地震資料的解釋仍主要利用一次反射波的信息．特別是目前應(yīng)用最為廣泛的地震反射波反演方法對波場特征的敏感性強(qiáng)，而多種波場同時存在則容易引起反演結(jié)果的不準(zhǔn)確性．單程波地震波數(shù)值模擬提供了一種適用于簡單介質(zhì)的一次反射波的構(gòu)建方法（何兵紅等，2009，2014；何兵紅，2011）．但對于鹽丘等復(fù)雜介質(zhì)，單程波法對于其地下地質(zhì)構(gòu)造的模擬精度較低．基于地震波散射理論的一階波恩近似波動方程在弱擾動假設(shè)下忽略多級散射，可以得到連續(xù)界面的一次反射波．與單程波波動方程相比，一階波恩近似波動方程對于高角度波場信息模擬精度更高．鄭江峰（2006）基于開爾文－佛格特模型推導(dǎo)了頻率域黏彈性介質(zhì)的波恩近似．本文基于廣義流變體模型推導(dǎo)了時間域黏滯性聲波介質(zhì)的速度－應(yīng)力方程，根據(jù)散射波理論得到了其一階波恩近似方程；提出了改進(jìn)的τ值法擬合常Q模型，由于該模型不再受應(yīng)力松弛時間與應(yīng)變延遲時間相等條件的約束，從而提高了低Q情況下常Q模型的擬合精度．卷積完全匹配層（convolutional perfect matched layer，簡寫為CPML）對大角度入射邊界反射吸收效果要優(yōu)于完全匹配層（perfect matched layer，簡寫為PML）（張魯新等，2010），為消除邊界反射進(jìn)一步得到了含有CPML邊界的時間域黏滯性介質(zhì)的一階波恩近似方程；數(shù)值試驗(yàn)對比了黏彈性介質(zhì)中一階波恩近似方程與全波波動方程以及單程波波動方程的模擬精度，并分析了一階波恩近似方程對速度擾動和地層Q擾動的適應(yīng)性．

1 基于流變體模型改進(jìn)的τ值法擬合常Q特征

實(shí)驗(yàn)證明地震波在介質(zhì)中傳播時的衰減與頻率相關(guān)，在地震頻帶內(nèi)地層品質(zhì)因子Q與頻率無關(guān)．為了在地震頻帶內(nèi)構(gòu)建與頻率無關(guān)的黏彈性模型，以麥克斯韋模型和開爾文－佛格特模型為基礎(chǔ)的廣義流變體模型得以發(fā)展．Moczo和Kristek（2005）、曹丹平（2008）以及Cao和Yin（2014）等證明了廣義麥克斯韋模型與廣義標(biāo)準(zhǔn)線性體模型（基于麥克斯韋模型的標(biāo)準(zhǔn)線性體固體模型（generalized Maxwell standard linear solid model，簡寫為GMSLS）、基于開爾文－佛格特的標(biāo)準(zhǔn)線性體固體模型（generalized Kelvin-Voigt standard linear solid model，簡寫為GKSLS））是等價的．

基于GMSLS模型的Q值可以表示為（Blanchetal，1995）

式中，ω為角頻率，L為標(biāo)準(zhǔn)線性黏彈體的個數(shù)，τεl和τσl分別為應(yīng)變松弛時間和應(yīng)力松弛時間．

為了獲得Q值的參數(shù)，定義松弛時間為

則式（1）重新整理為

強(qiáng)衰減低Q介質(zhì)中不滿足假設(shè)條件τl＋1≈1，對于常Q模型擬合存在很大誤差；常規(guī)τ值法穩(wěn)定性雖得到提高，但其精度不高．本文為了提高Q值擬合精度，定義函數(shù)

則式（3）可表示為

假設(shè)真實(shí)品質(zhì)因子為Qtrue，根據(jù)最小二乘法建立目標(biāo)函數(shù)：

式中，ωs和ωe分別為積分范圍的開始角頻率和截止角頻率．假設(shè)每個標(biāo)準(zhǔn)線性黏彈體τl相等，即τl＝τ．對τ求導(dǎo)并令其為零，即

得到τ的計(jì)算表達(dá)式為

即為改進(jìn)的τ值法．Blanch等（1995）提出的常規(guī)τ值表達(dá)式為

改進(jìn)的τ值計(jì)算式中，H（ω）／F（ω）隨頻率變化的特征能夠?qū)（ω）進(jìn)行與頻率有關(guān)的修正．當(dāng)L＝5時，常規(guī)τ值法與改進(jìn)的τ值法在地震頻帶內(nèi)均能夠?qū)Τ模型進(jìn)行較高精度的擬合．本文采用L＝4對改進(jìn)的τ值法進(jìn)行測試．由測試結(jié)果（圖1）可知，改進(jìn)的τ值法所確立的Q模型更能接近常Q模型，特別是在中高頻部分對常規(guī)的τ值法進(jìn)行了很大的修正．

圖1 改進(jìn)的τ值法與常規(guī)τ值法常Q模型擬合精度對比（a）Q＝30；（b）Q＝50Fig．1 Comparison of constant Qfitting by improvedτmethod with that by conventionalτmethod（a）Q＝30；（b）Q＝50

2 基于廣義流變體模型的位移方程的建立

根據(jù)玻爾茲曼疊加原理，在線性黏彈性介質(zhì)中應(yīng)力－應(yīng)變方程為

式中，σ為應(yīng)力，ε為應(yīng)變，＊為卷積算子，ψ為松弛函數(shù)．黏滯性聲波介質(zhì)中松弛函數(shù)可以分別用松弛模量和未松弛模量來表示：

本文采用未松弛模量表達(dá)式．

根據(jù)定義，時間域模量是松弛函數(shù)對時間的導(dǎo)數(shù)，可得

將式（15）帶入式（12）得到用模量表示的應(yīng)力方程為

為了解決時間域卷積形式求解問題，Day和Minster（1984）基于Padé近似引入記憶變量的概念，該思想被逐步接納并廣泛應(yīng)用．

定義記憶變量為

并對時間t求導(dǎo)可得

則式（16）可表示為

將應(yīng)力、應(yīng)變與位移之間的關(guān)系分別代入式（18）和式（19）中，可得黏滯性位移方程為

式中，體積模量K＝c2ρ，c為地震波傳播速度．

3 基于地震波散射理論的波恩近似方程的建立

攝動理論假設(shè)下，地下介質(zhì)可以劃分為背景介質(zhì)和擾動介質(zhì)，其對應(yīng)的地震波場分解為背景波場和擾動波場．在黏滯性聲波介質(zhì)中，介質(zhì)參數(shù)包含彈性參數(shù)和黏滯性參數(shù)．本文中彈性參數(shù)包含體積模量K和密度ρ．GMSLS模型中黏滯性參數(shù)用τεl和τσl表示．τσl和τεl與Q值的對應(yīng)關(guān)系通過式（2）和式（9）獲得，因此在τσl固定的情況下，所對應(yīng)的有3個背景介質(zhì)參數(shù)和3個擾動介質(zhì)參數(shù)．本文定義

黏彈性介質(zhì)中地震波場除了總壓力波場，還需要計(jì)算黏滯性波場．其對應(yīng)的波場分為2個背景波場和2個擾動波場．根據(jù)地震波散射理論，總波場可以分解為入射波場和散射波場，即

對于變密度介質(zhì)，位移方程式（20）重新表示為

其中可分離出背景波場

以及包含多級散射的散射波場

在弱散射條件下忽略多級散射，用背景場代替總波場，背景介質(zhì)參數(shù)代替總介質(zhì)參數(shù)，可得到基于一階波恩近似的散射波場為

波場對于時間一階導(dǎo)數(shù)的存在，導(dǎo)致常規(guī)網(wǎng)格數(shù)值求解可執(zhí)行性差，通常利用速度、應(yīng)力與位移三者之間的關(guān)系，即

將位移方程式（20）降階為應(yīng)力－速度方程：

該方程將所有關(guān)于時間和空間的二階導(dǎo)數(shù)降低為一階導(dǎo)數(shù)，從而可采用二維交錯網(wǎng)格進(jìn)行數(shù)值求解．但同時在該方程中也需引入變量vx和vz來分別表示x方向和z方向的速度分量，因此黏滯性聲波介質(zhì)中實(shí)際存在4個波場．

根據(jù)地震波散射理論，總波場可以分解為入射波場和散射波場，即

本文定義如下變量：

將式（30）中的各個變量帶入應(yīng)力－速度方程式（28）中，可將總波場分解為背景波場和散射波場．同時在散射波場計(jì)算中用背景波場代替總波場，背景介質(zhì)參數(shù)代替總介質(zhì)參數(shù)，可得到基于一階波恩近似的散射波場滿足的應(yīng)力－速度方程：

4 非分裂式CPML邊界條件

地震波數(shù)值模擬對實(shí)際地球介質(zhì)進(jìn)行了有限區(qū)域假設(shè)，在區(qū)域邊界通常會產(chǎn)生人為的邊界反射．Berenger（1994）提出的PML吸收邊界條件從電磁場模擬引入到地震波模擬中，因其能夠有效地消除邊界反射的影響而成為廣泛應(yīng)用的一種邊界條件．研究表明，PML邊界條件對于地震波的吸收具有方向性，垂直入射的波能夠在較短的時間內(nèi)衰減一個數(shù)量級，而對于大角度入射的波則吸收效果欠佳．目前，地震波的研究多數(shù)是基于地面地震反射波法．當(dāng)震源位于地層表面時，遠(yuǎn)偏移距波場以大入射角到達(dá)PML吸收邊界區(qū)域，邊界反射部分被吸收．Kuzuoglu和Mittra（1996）提出了復(fù)頻移PML邊界條件，其在時間域表現(xiàn)為卷積的形式，因此又稱為CPML邊界條件．

PML邊界的實(shí)施是通過將常規(guī)坐標(biāo)x用含有阻尼項(xiàng)dx的復(fù)坐標(biāo)來代替，即

在CPML邊界條件中，引入實(shí)數(shù)項(xiàng)κx和αx：

式中，α≥0，κ≥1．該式中對復(fù)坐標(biāo)求導(dǎo)可得

通過傅里葉反變換得到時間域形式為

通過遞歸得到迭代解為

將式（35）和（38）帶入式（31），得到基于CPML邊界條件的黏滯性介質(zhì)一階波恩近似方程：

采用CPML邊界條件避免了將原始方程進(jìn)行分解的問題，只需要增加輔助變量φ即可大大提高計(jì)算效率．

5 數(shù)值測試

為研究黏滯性聲波波動方程吸收衰減特征和基于一階波恩近似的黏滯性聲波介質(zhì)散射波波場特征，本文建立了含有不整合界面的理論斷層速度場．該模型縱測線長度為3 640m，深度為2 380m，采樣間隔為10m．炮點(diǎn)位于地表測線方向2 000m處，檢波器以10m間隔均勻分布于地表．弱速度擾動模型中最低速度為2 600m／s，最高速度為3 900m／s（圖2a）．弱Q擾動模型中Q值在李氏經(jīng)驗(yàn)公式（李慶忠，1993）計(jì)算所得結(jié)果的基礎(chǔ)上加入深度加權(quán)系數(shù)以增強(qiáng)淺層地層吸收特征（圖2b），淺層低Q值特征表示地表疏松介質(zhì)的強(qiáng)衰減特征．在檢波點(diǎn)號范圍為180—220、深度約為750m的斷層構(gòu)造中，其所對應(yīng)的地層低Q特征表示了含氣特征．

圖2 弱速度擾動（a）和弱Q擾動（b）模型Fig．2 The model with weak velocity perturbation（a）and weak Qperturbation（b）

根據(jù)散射理論和量子力學(xué)理論，在弱散射條件下采用波恩近似用背景波場代替總波場得到散射波場．一階波恩近似方程是對一級散射波場的描述，某個散射點(diǎn)處的散射波場可當(dāng)做二次信號的激發(fā)源（吳如山，1993）．由于地球介質(zhì)中的散射問題滿足輻射邊界條件，即當(dāng)?shù)卣鸩▊鞑サ綗o限遠(yuǎn)處時波場能量為零，因此可以利用格林函數(shù)求解方法得到其精確解的解析表達(dá)式．基于格林函數(shù)的散射波場求解既可以采用雙程波思路，也可以采用波場分解理論將散射波場分為前向散射和背向散射，再利用單程波思路計(jì)算散射波場．本文采用時間域有限差分法利用雙程波思路求解波恩近似方程，其計(jì)算簡單且無須進(jìn)行傅里葉變換．?dāng)?shù)值實(shí)驗(yàn)中全波波動方程包含了背景波場和散射波場（包含多級散射）（圖3a，b），在連續(xù)界面上基于波恩近似的一級散射波場表現(xiàn)為一次反射波（圖3c，d）．由全波波動方程與一階波恩近似方程所得到的波場差值中包含了透射波場和多級散射波場（圖3e，f）．在反射波偏移成像中廣泛使用的單程波法波場延拓能夠?qū)崿F(xiàn)主要反射波的數(shù)值模擬，是目前工業(yè)生產(chǎn)中主要使用的偏移成像方法之一．本文基于分步傅里葉算法得到了基于單程波波動方程的黏滯性介質(zhì)波場（圖4a，b）．?dāng)?shù)值結(jié)果表明，單程波法在低角度入射條件下能夠獲得較高的精度．但是對于高角度入射的情況，單程波波動方程與全波波動方程和一階波恩近似方程相比，其精度較低．

圖3 弱速度擾動及弱Q擾動模型中由全波波動方程（a，b）和一階波恩近似方程（c，d）得到的黏滯性介質(zhì)波場以及由全波波動方程與一階波恩近似方程得到的波場差值（e，f）（a），（c），（e）中t＝420ms；（b），（d），（f）中t＝558msFig．3 The viscoacoustic wave fields obtained by different wave equations in the model with weak velocity perturbation and weak Qperturbation Figs．（a）and（b）are the wave fields obtained by full wave equation；Figs．（c）and（d）are the wave fields obtained by first-order Born approximation equation；Figs．（e）and（f）are the differences between the wave fields in（a），（b）and those in（c），（d）．In Figs．（a），（c），（e）t＝420ms and in Figs．（b），（d），（f）t＝558ms

圖4 黏滯性介質(zhì)單程波波動方程得到的波場．（a）t＝420ms；（b）t＝558msFig．4 The wave fields obtained by one-way wave equation in viscoacoustic media（a）t＝420ms；（b）t＝558ms

圖5 弱速度擾動及弱Q擾動模型中由全波波動方程（a）、一階波恩近似方程（b）和單程波波動方程（c）得到的地震記錄Fig．5 The seismic records obtained by full wave equation（a），first-order Born approximation equation（b）and one-way wave equation（c）in the model with weak velocity perturbation and weak Qperturbation

圖5為弱速度擾動及弱Q擾動模型中得到的單炮地震記錄，其中全波波動方程與一階波恩近似方程所得地震記錄在相同幅值范圍內(nèi)．由于單程波波動方程所得地震記錄振幅能量范圍與前兩者相差較大，我們對其進(jìn)行整體加權(quán)，使炮點(diǎn)位置處的檢波器接收到的第一個反射界面的反射波振幅與對應(yīng)的全波波動方程所得該反射界面的反射波振幅一致．由于利用一階波恩近似方程得到的波場不包含透射波，因此地面地震接收記錄中不包含直達(dá)波信息．由于直達(dá)波能量較強(qiáng)，當(dāng)其傳播到邊界處時的反射能量較強(qiáng)．特別是對于弱散射介質(zhì)，直達(dá)波在邊界處的反射能量甚至能與一次散射波場能量相當(dāng)（圖5a），從而增加了反射波地震勘探的復(fù)雜性．

由于本文單程波數(shù)值模擬中未考慮反射系數(shù)與偏移距及入射角的關(guān)系，故單程波波場在近炮點(diǎn)位置能夠取得比較理想的效果．但在遠(yuǎn)偏移距處，采用單程波波動方程得到的反射波能量強(qiáng)度與全波波動方程及一階波恩近似方程的明顯不一致．因此如何準(zhǔn)確地計(jì)算不同入射角下復(fù)雜介質(zhì)界面的反射系數(shù)是單程波法準(zhǔn)確描述反射波場要解決的關(guān)鍵問題．

圖6為近炮點(diǎn)處的檢波點(diǎn)所接收到的一道地震記錄．全波波動方程中包含直達(dá)波信息，當(dāng)?shù)谝粋€反射界面較淺時，直達(dá)波與反射波發(fā)生干涉．在常規(guī)處理中通常用切除直達(dá)波的方法得到反射波信息，但這種切除方法也損失了重要的反射波信息．特別是在淺層疏松介質(zhì)中直達(dá)波與反射波干涉尤為嚴(yán)重時．當(dāng)切除直達(dá)波，有效的反射波信息也基本被切除．而利用基于一階波恩近似的波動方程則能夠得到只含有一次反射波的波場，從而避免了直達(dá)波與反射波分離的問題．從能量分布上看，基于一階波恩近似方程得到的反射波場與全波波動方程得到的波場吻合較好．單程波法能夠直接得到反射波信息，但在能量匹配上與全波波動方程不一致，特別是對于深層陡傾角的界面，其反射波能量誤差較大．

圖6 弱速度擾動及弱Q擾動模型中由全波波動方程、一階波恩近似方程與單程波波動方程得到的地震波形對比Fig．6 Comparison of the seismic waveforms obtained by full wave equation（blue lines）with that by first-order Born approximation wave equation（red lines）and one-way wave equation（black lines）in the model with weak velocity perturbation and weak Qperturbation

圖7 強(qiáng)速度擾動模型Fig．7 The model with strong velocity perturbation

黏滯性介質(zhì)一階波恩近似方程得到的波場與全波波動方程得到的波場匹配的前提條件是弱散射近似，即波恩近似成立的條件為kRδσ／σ0≤1．其中，δσ為散射波場，σ0為背景波場，k為波場波數(shù)，R為非均勻體的尺度（如球半徑）（吳如山，1993）．在以上討論中我們所使用的速度模型在層速度及厚度上均滿足波恩近似成立的條件．為了進(jìn)一步探索在不滿足波恩近似成立條件下的地震波場誤差問題，本文將速度進(jìn)行修改得到強(qiáng)速度擾動模型．修改后的速度模型最高速度為3 900m／s，最低速度降低為1 800m／s，相鄰兩層速度差增大，特別是對于不整合覆蓋層，其速度明顯降低（圖7）．在地震記錄剖面上，黏滯性介質(zhì)全波波動方程與一階波恩近似方程所得到的波場匹配較好，其波形形態(tài)與振幅能量保持一致（圖8）．分別從圖8a，b中抽取檢波點(diǎn)號為200的一道地震記錄，得到如圖9所示的波形對比圖．通過對比可知，一階波恩近似方程與全波波動方程所得到的波場在振幅方面能夠達(dá)到一致．但在旅行時方面，一階波恩近似方程所得到的反射波場發(fā)生漂移，特別是在t＝538ms處的反射波旅行時向上移動約10ms，這種誤差將導(dǎo)致地震層位解釋不準(zhǔn)確，嚴(yán)重影響儲層預(yù)測精度．

圖8 強(qiáng)速度擾動及弱Q擾動模型中由全波波動方程（a）和一階波恩近似方程（b）得到的地震記錄Fig．8 The seismic records obtained by full wave equation（a）and first-order Born approximation equation（b）in the model with strong velocity perturbation and weak Qperturbation

圖9 強(qiáng)速度擾動及弱Q擾動模型中由全波波動方程與一階波恩近似方程得到的地震波形對比Fig．9 Comparison of the seismic waveforms obtained by full wave equation（blue lines）with that by first-order Born approximation equation（red lines）in the model with strong velocity perturbation and weak Qperturbation

在常密度假設(shè)條件下，黏滯性參數(shù)Q是另一個影響波恩近似精度的參數(shù)．在以上討論中我們固定了黏滯性參數(shù)Q，通過改變速度參數(shù)探索該黏滯性介質(zhì)一階波恩近似方程所得地震記錄中的反射波精度．下面進(jìn)一步固定速度參數(shù)，研究黏滯性參數(shù)Q對一階波恩近似方程所得反射波精度的影響．對于全波波動方程，為了將直達(dá)波與反射波區(qū)分開，我們把第一個反射界面深度下移一定深度，Q值最低值為20、最高值為500，兩個平層之間Q值差為40；第二個平層（Q＝60）與相鄰不整合層的最大Q值相差280（圖10a）．圖11為在均勻速度場的情況下，由黏滯性擾動所產(chǎn)生的波場特征．盡管Q值非均質(zhì)性比較強(qiáng)，全波波動方程與一階波恩近似方程所得波場在旅行時上仍然能夠保持一致，其振幅也基本一致．

圖10 強(qiáng)Q擾動模型（a）和修改后的弱速度擾動模型（b）Fig．10 The models with strong Qperturbation（a）and the modified weak velocity perturbation（b）

圖11 均勻速度場下強(qiáng)Q擾動模型中由全波波動方程與一階波恩近似方程得到的地震波形對比Fig．11 Comparison of the seismic waveforms obtained by full wave equation（blue lines）with that by first-order Born approximation equation（red lines）in homogenous velocity model with strong Qperturbation

圖12 弱速度擾動及強(qiáng)Q擾動模型中由全波波動方程與一階波恩近似方程得到的地震波形對比Fig．12 Comparison of the seismic waveforms obtained by full wave equation（blue lines）with that by first-order Born approximation equation（red lines）in the model with weak velocity perturbation and strong Qperturbation

為了消除全波波動方程中直達(dá)波對反射波的影響，以便與一階波恩近似方程所得的反射波進(jìn)行對比，我們將弱速度擾動模型（圖2a）第一個反射界面下移相同的深度得到修改后的弱速度擾動場（圖10b）．通過對比可知，弱速度擾動和強(qiáng)Q擾動的情況下，全波波動方程與一階波恩近似方程所得波場在旅行時方面能夠匹配較好（圖12），但在振幅能量上略有差異．因此在黏滯性介質(zhì)中，一階波恩近似方程的波場精度主要取決于速度，Q值擾動所引起的振幅誤差幾乎可以忽略不計(jì)．

6 討論與結(jié)論

針對常規(guī)τ值法計(jì)算常Q模型參數(shù)精度低的問題，本文推導(dǎo)了改進(jìn)的τ值法．該τ值法能夠?qū)ΤＲ?guī)τ值法不同頻率成分的Q值擬合精度進(jìn)行修正，得到更加精確的常Q模型．在此基礎(chǔ)上，推導(dǎo)了包含CMPL邊界條件的基于廣義流變體模型的波恩近似位移方程．由于黏滯性項(xiàng)引入了時間一階導(dǎo)數(shù)項(xiàng)，常規(guī)網(wǎng)格方法將波場參數(shù)和介質(zhì)參數(shù)定義在同一網(wǎng)格點(diǎn)上而不利于黏滯性介質(zhì)位移方程求解．本文根據(jù)應(yīng)力、應(yīng)變及位移之間的關(guān)系，進(jìn)一步推導(dǎo)了基于交錯網(wǎng)格的包含CPML邊界條件的變密度一階波恩近似應(yīng)力－速度方程．在弱散射條件下，基于波恩近似的黏滯性聲波介質(zhì)方程對于散射波場的模擬能夠與全波波動方程得到的波場匹配．對于全波波動方程，在弱散射條件下，直達(dá)波邊界反射能量的強(qiáng)度甚至能與散射波波場相匹配；對于淺層反射界面，由于直達(dá)波與反射波發(fā)生干涉，常規(guī)的直達(dá)波切除方法不僅損失了反射波信息，而且嚴(yán)重影響了地震解釋精度．黏滯性介質(zhì)波恩近似方程由于不含透射波，不存在直達(dá)波的影響，從而有利于反射波勘探理論研究．相對于單程波法，黏滯性介質(zhì)一階波恩近似方程得到的波場無須計(jì)算反射系數(shù)與角度的關(guān)系，能夠自動進(jìn)行能量分配．當(dāng)速度場存在強(qiáng)擾動時，黏滯性介質(zhì)一階波恩近似方程與全波波動方程得到的波場在旅行時上存在較大誤差，需要采用高階多級散射波波動方程；當(dāng)速度場為弱擾動、模型Q值為強(qiáng)擾動時，盡管黏滯性所引起的速度頻散會帶來相速度變化，但黏滯性介質(zhì)一階波恩近似方程與全波波動方程所得波場基本吻合，說明其在弱速度擾動下能夠適應(yīng)地層Q值強(qiáng)擾動．

畢玉英，楊寶俊．1995．二維粘彈介質(zhì)中完全波場正演模擬［J］．石油地球物理勘探，30（3）：351-362．

Bi Y Y，Yang B J．1995．Forward modeling of full wave field in 2-D viscoelastic medium［J］．OilGeophysicalProspecting，30（3）：351-362（in Chinese）．

曹丹平．2008．多尺度地震資料正反演方法研究［D］．青島：中國石油大學(xué)（華東）地球科學(xué)與技術(shù)學(xué)院：16-86．

Cao D P．2008．MethodResearchofMultiscaleSeismicDataModelingandInversion［D］．Qingdao：School of Geosciences，China University of Petroleum （East China）：16-86（in Chinese）．

崔建軍，何繼善．2001．粘彈性波動方程正演和偏移［J］．中南工業(yè)大學(xué)學(xué)報(bào)：自然科學(xué)版，32（5）：441-444．

Cui J J，He J S．2001．Viscoelastic wave equation forward modeling and migration［J］．JournalofCentralSouthUniversity：ScienceandTechnology，32（5）：441-444（in Chinese）．

杜啟振，王延光，付水華．2006．方位各向異性粘彈性介質(zhì)波場數(shù)值模擬［J］．地球物理學(xué)進(jìn)展，21（2）：502-504．

Du Q Z，Wang Y G，F(xiàn)u S H．2006．Wavefield forward modeling with the pseudo-spectral method in viscoelastic and azimuthally anisotropic media［J］．ProgressinGeophysics，21（2）：502-504（in Chinese）．

郭少華．2004．基于 Kelvin-Voigt模型的粘彈性波動力學(xué)的本征化理論［J］．應(yīng)用數(shù)學(xué)和力學(xué)，25（7）：723-728．

Guo S H．2004．Eigen theory of viscoelastic dynamics based on the Kelvin-Voigt model［J］．AppliedMathematicsand Mechanics，25（7）：723-728（in Chinese）．

何兵紅，吳國忱，梁鍇，白曉寅．2009．粘彈性介質(zhì)單程波法非零偏移距地震數(shù)值模擬與偏移［J］．地球物理學(xué)進(jìn)展，24（4）：1299-1312．

He B H，Wu G C，Liang K，Bai X Y．2009．Numerical simulation and migration of non-zero offset seismic wave in viscoelastic medium by one-way wave equation［J］．ProgressinGeophysics，24（4）：1299-1312（in Chinese）．

何兵紅．2011．基于衰減介質(zhì)的地震波數(shù)值模擬及吸收屬性提取方法研究［D］．青島：中國石油大學(xué)（華東）地球科學(xué)與技術(shù)學(xué)院：36-50．

He B H．2011．TheStudyofNumericalSimulationofSeismicWaveinAttenuationMediumandExtractionofAbsorptionProperties［D］．Qingdao：School of Geosciences，China University of Petroleum （East China）：36-50（in Chinese）．

何兵紅，吳國忱，許沖．2014．基于FFD算子的衰減介質(zhì)地震波模擬［J］．石油地球物理勘探，49（1）：153-160．

He B H，Wu G C，Xu C．2014．Seismic wave simulation in attenuation medium based on FFD operator［J］．OilGeophysicalProspecting，49（1）：153-160（in Chinese）．

李慶忠．1993．走向精確勘探的道路：高分辨率地震勘探系統(tǒng)工程剖析［M］．北京：石油工業(yè)出版社：31-44．

Li Q Z．1993．TheWaytoObtainaBetterResolutioninSeismicProspecting：AsystematicalAnalysisofHighResolutionSeismicExploration［M］．Beijing：Petroleum Industry Press：31-44（in Chinese）．

牛濱華，孫春巖．2004．地震波理論研究進(jìn)展：介質(zhì)模型與地震波傳播［J］．地球物理學(xué)進(jìn)展，19（2）：255-263．

Niu B H，Sun C Y．2004．Developing theory of propagation of seismic waves：Medium model and propagation of seismic waves［J］．ProgressinGeophysics，19（2）：255-263（in Chinese）．

單啟銅．2007．粘彈性波動方程正演模擬與參數(shù)反演［D］．青島：中國石油大學(xué)（華東）地球科學(xué)與技術(shù)學(xué)院：13-83．

Shan Q T．2007．SimulationandInversionofWaveEquationinViscoelasticMedia［D］．Qingdao：School of Geosciences，China University of Petroleum：13-83（in Chinese）．

孫成禹，肖云飛，印興耀，彭洪超．2010．黏彈介質(zhì)波動方程有限差分解的穩(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）．

吳如山．1993．地震波的散射與衰減［M］．北京：地震出版社：49-80．

Wu R S．1993．ScatterandAttenuationofSeismicWave［M］．Beijing：Seismological Press：49-80（in Chinese）．

楊文采．1986．粘彈性介質(zhì)中反射地震道的合成［J］．石油地球物理勘探，21（6）：615-623．

Yang W C．1986．The synthesis of reflection seismic traces in visco-elastic medium［J］．OilGeophysicalProspecting，21（6）：615-623（in Chinese）．

苑春方，彭蘇萍，張中杰，劉振寬．2005．Kelvin-Voigt均勻黏彈性介質(zhì)中傳播的地震波［J］．中國科學(xué)：D輯，35（10）：53-58．

Yuan C F，Peng S P，Zhang Z J，Liu Z K．2006．Seismic wave propagating in Kelvin-Voigt homogeneous visco-elastic media［J］．ScienceinChina：SeriesD，49（2）：147-153．

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

Zhang L X，F(xiàn)u 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）．

鄭江峰．2006．粘彈性介質(zhì)中提高地震資料分辨率方法研究［D］．大慶：大慶石油學(xué)院地球科學(xué)學(xué)院：30-32．

Zheng J F．2006．MethodologyandResearchonImprovingSeismicResolutioninViscoelasticMedia［D］．Daqing：Geoscience College，Daqing Petroleum Institute：30-32（in Chinese）．

Berenger J P．1994．A perfectly matched layer for the absorption of electromagnetic waves［J］．JComputPhys，114（2）：185-200．

Blanch J O，Robertsson J O A，Symes W W．1995．Modeling of a constantQ：Methodology and algorithm for an efficient and optimally inexpensive viscoelastic technique［J］．Geophysics，60（1）：176-184．

Cao D P，Yin X Y．2014．Equivalence relations of generalized rheological models for viscoelastic seismic-wave modeling［J］．BullSeismolSocAm，104（1）：260-268．

Carcione J M，Kosloff D，Kosloff R．1988．Wave propagation simulation in a linear viscoacoustic medium［J］．GeophysJ Int，93（2）：393-401．

Carcione J M．1993．Seismic modeling in viscoelastic media［J］．Geophysics，58（1）：110-120．

Day S M，Minster J B．1984．Numerical simulation of attenuated wavefields using a Padéapproximant method［J］．GeophysJInt，78（1）：105-118．

Emmerich H，Korn M．1987．Incorporation of attenuation into time-domain computations of seismic wave fields［J］．Geophysics，52（9）：1252-1264．

Kristek J，Moczo P．2003．Seismic-wave propagation in viscoelastic media with material discontinuities：A 3Dfourthorder staggered-grid finite-difference modeling［J］．BullSeismolSocAm，93（5）：2273-2280．

Kuzuoglu M，Mittra R．1996．Mesh truncation by perfectly matched anisotropic absorbers in the finite-element method［J］．MicrowOptTechnolLett，12（3）：136-140．

Moczo P，Kristek J．2005．On the rheological models used for time-domain methods of seismic wave propagation［J］．GeophysResLett，32（1）：L1306．

Robertsson J O A，Blanch J O，Symes W W．1994．Viscoelastic finite-difference modeling［J］．Geophysics，59（9）：1444-1456．

?tekl I，Pratt R G．1998．Accurate viscoelastic modeling by frequency-domain finite differences using rotated operators［J］．Geophysics，63（5）：1779-1794．