基于流變體模型的波恩近似方程地震正演模擬*

何兵紅 吳國(guó)忱

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

引言

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

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

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

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

實(shí)驗(yàn)證明地震波在介質(zhì)中傳播時(shí)的衰減與頻率相關(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，簡(jiǎn)寫為GMSLS）、基于開爾文－佛格特的標(biāo)準(zhǔn)線性體固體模型（generalized Kelvin-Voigt standard linear solid model，簡(jiǎn)寫為GKSLS））是等價(jià)的．

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

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

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

則式（1）重新整理為

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

則式（3）可表示為

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

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

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

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

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

圖1 改進(jìn)的τ值法與常規(guī)τ值法常Q模型擬合精度對(duì)比（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í)間域模量是松弛函數(shù)對(duì)時(shí)間的導(dǎo)數(shù)，可得

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

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

定義記憶變量為

并對(duì)時(shí)間t求導(dǎo)可得

則式（16）可表示為

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

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

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

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

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

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

其中可分離出背景波場(chǎng)

以及包含多級(jí)散射的散射波場(chǎng)

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

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

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

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

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

本文定義如下變量：

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

4 非分裂式CPML邊界條件

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

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

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

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

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

通過遞歸得到迭代解為

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

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

5 數(shù)值測(cè)試

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

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

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

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

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

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

圖6 弱速度擾動(dòng)及弱Q擾動(dòng)模型中由全波波動(dòng)方程、一階波恩近似方程與單程波波動(dòng)方程得到的地震波形對(duì)比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)速度擾動(dòng)模型Fig．7 The model with strong velocity perturbation

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

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

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

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

6 討論與結(jié)論

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

畢玉英，楊寶俊．1995．二維粘彈介質(zhì)中完全波場(chǎng)正演模擬［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］．青島：中國(guó)石油大學(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．粘彈性波動(dòng)方程正演和偏移［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ì)波場(chǎng)數(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模型的粘彈性波動(dòng)力學(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）．

何兵紅，吳國(guó)忱，梁鍇，白曉寅．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］．青島：中國(guó)石油大學(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）．

何兵紅，吳國(guó)忱，許沖．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．粘彈性波動(dòng)方程正演模擬與參數(shù)反演［D］．青島：中國(guó)石油大學(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ì)波動(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）．

吳如山．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］．中國(guó)科學(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)交錯(cuò)網(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．