吳國忱 秦海旭

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

引言

全球發(fā)現(xiàn)的油氣藏中有大量的裂縫性油氣藏，包括碳酸鹽裂縫型油氣藏、致密砂巖裂縫型油氣藏及基巖裂縫型油氣藏（毛寧波等，2008；楊曉等，2010）.裂縫通常是裂縫性油氣藏的儲存空間或運(yùn)移通道，因此裂縫介質(zhì)中地震波場的研究越來越受到關(guān)注（謝桂生，1994；金抒辛，何樵登，2005；潘保芝等，2006）.研究發(fā)現(xiàn)，在含有定向分布的裂縫介質(zhì)中地震波的傳播呈現(xiàn)各向異性的特征（張中杰，何樵登，1989；郭建，1993；Schoenberg，Sayers，1995）.因此，利用地震波場來研究裂縫發(fā)育帶位置、確定儲層物性并進(jìn)行儲層預(yù)測具有重要的意義.

目前研究裂縫的方法之一是將裂縫看成散射體，用散射理論研究裂縫介質(zhì)的地震響應(yīng)（刁順等，1994a，b），但此方法的正演效率非常低.本文用等效理論將裂縫介質(zhì)等效為各向異性介質(zhì)，運(yùn)用各向異性介質(zhì)理論研究裂縫的地震響應(yīng).各向異性是指介質(zhì)的彈性參數(shù)隨方向而變化，該方法效率高，能得到較好的效果.裂縫對地震波的影響已研究了多年.例如：Hudson（1981）提出裂縫介質(zhì)中P波和S波速度與衰減的公式；Thomsen（1986）提出利用彈性弱各向異性理論研究定向裂縫與等徑孔隙度對地震速度的綜合影響；Schoenberg和Sayers（1995）使用裂縫法向與切向柔度推導(dǎo)裂縫介質(zhì)中的有效彈性常數(shù)；Bakulin等（2000a，b，c）回顧了幾種不同的裂縫等效理論，并比較了它們對裂縫參數(shù)估計(jì)的可行性.

地震波場數(shù)值模擬是地球物理學(xué)領(lǐng)域的一個研究熱點(diǎn)，既可以提高我們對各種復(fù)雜介質(zhì)中地震波傳播規(guī)律的認(rèn)識，又可以作為檢測工具檢驗(yàn)各種方法技術(shù)的應(yīng)用效果.前人研究結(jié)果表明，地震波在定向裂縫介質(zhì)中的傳播具有橫向各向同性的特征（Hudson，1981；Thomsen，1986；Bakulin et al，2000a，b，c）.Crampin（1985）提出構(gòu)造應(yīng)力產(chǎn)生的具有水平對稱軸的橫向各向同性（horizontal transverse isotropy，簡稱為HTI）介質(zhì)，可用來模擬高陡傾角裂縫介質(zhì)（孔麗云等，2012）.目前常用的裂縫介質(zhì)正演模擬方法主要是交錯網(wǎng)格有限差分方法（董良國等，2000；裴正林，2004，2006；裴正林，王尚旭，2005；李雪靜，劉定進(jìn)，2008）.普通交錯網(wǎng)格能有效地模擬地震波在裂縫介質(zhì)中的傳播，但也存在一定的缺點(diǎn).首先，當(dāng)模擬非均勻性比較強(qiáng)的介質(zhì)或各向異性介質(zhì)時，為了獲得滿意的計(jì)算結(jié)果，需要對介質(zhì)參數(shù)進(jìn)行平均或內(nèi)插；其次，在遇到自由界面時，必須對自由界面進(jìn)行單獨(dú)的顯式處理（嚴(yán)洪勇，劉洋，2012）.旋轉(zhuǎn)交錯網(wǎng)格將相同的物理量定義在相同的網(wǎng)格點(diǎn)上，計(jì)算其沿對角線的差分，然后利用對角線差分組合計(jì)算其沿坐標(biāo)軸的差分.該方法不需要對介質(zhì)參數(shù)進(jìn)行平均或內(nèi)插，在遇到自由邊界時不需要單獨(dú)的顯示處理（Saenger et al，2000），從而解決了普通交錯網(wǎng)格存在的缺點(diǎn)，是目前正演中比較準(zhǔn)確的方法.本文采用旋轉(zhuǎn)交錯網(wǎng)格正演模擬方法對裂縫性儲集層進(jìn)行正演模擬.

1 裂縫介質(zhì)的等效理論

運(yùn)用各向異性理論進(jìn)行裂縫介質(zhì)正演模擬，需要將裂縫介質(zhì)等效為各向異性介質(zhì).目前廣泛存在的裂縫等效理論主要有Thomsen等效理論、Hudson等效理論及線性滑動理論等.Thomsen等效理論相當(dāng)于低頻模型，假設(shè)流體流動無限慢，該方法限制條件較多，不利于研究（Thomsen，1986）.Hudson等效理論與Thomsen等效理論正好相反，相當(dāng)于高頻模型，其假設(shè)流體流動無限快，也不利于實(shí)際研究（Hudson，1999）.線性滑動理論（Bakulin et al，2000a，b，c）是將裂縫看作地層中的一個柔性面，該柔性面可以用裂縫的切向柔度和法向柔度進(jìn)行表示，柔度符合線性滑動邊界條件便于求解.對于多組裂縫可以對裂縫柔度進(jìn)行矢量疊加便于實(shí)際求解.因此本文采用線性滑動等效理論將裂縫介質(zhì)等效為各向異性介質(zhì).

含有裂縫巖石的總?cè)岫鹊扔诹芽p柔度與圍巖柔度的和，即

式中，S為巖石總的柔度矩陣.Sf為裂縫的柔度矩陣.Sb為圍巖的柔度矩陣.其中柔度矩陣為彈性矩陣C的倒數(shù)，即S＝C－1.根據(jù)線性滑動理論，裂縫介質(zhì)的彈性矩陣為

式中，λ和μ為背景巖石的拉梅常數(shù)，ΔN和ΔT分別為裂縫介質(zhì)的法向柔度與切向柔度.當(dāng)介質(zhì)中存在幾組裂縫時，先求出每組裂縫介質(zhì)的切向柔度與法向柔度，然后對切向柔度與法向柔度進(jìn)行矢量疊加求出裂縫介質(zhì)總的柔度矩陣.裂縫介質(zhì)等效示意圖如圖1所示.

由于裂縫介質(zhì)的法向柔度與切向柔度需要滿足線性滑動邊界條件，不方便求解.而采用Hudson等效理論與線性滑動理論相結(jié)合，既彌補(bǔ)了Hudson等效理論的局限性，又可以方便計(jì)算裂縫的法向柔度與切向柔度.實(shí)際計(jì)算中采用線性滑動理論與Hudson理論相結(jié)合的等效理論將裂縫介質(zhì)等效為各向異性介質(zhì).此時，裂縫的法向柔度和切向柔度分別為

圖1 裂縫介質(zhì)等效示意圖圖中ΔN和ΔT分別為裂縫介質(zhì)的法向柔度與切向柔度.（a）單縫模型；（b）裂縫組合模型Fig.1 The sketch of fracture equivalence ΔNdenotes the normal weakness，and ΔTdenotes the tangential weakness（a）Sigle fracture model；（b）Combination model of fractures

式中，e為裂縫密度，k′和μ′分別為裂縫充填物的體積模量與切變模量，k為背景介質(zhì)的體積模量，a為裂縫的長度，c為裂縫的寬度，g＝v2S／v2P.彈性介質(zhì)模型的性質(zhì)由剛度矩陣C確定，C確定了應(yīng)力與應(yīng)變之間的關(guān)系，但其確定彈性波動方程系數(shù)的物理意義很不直觀.Thomsen（1986）提出了一整套表征各向異性的參數(shù)，具體表達(dá)如下：

式中VP0，VS0分別為qP波和qS波沿各向異性介質(zhì)對稱軸方向的相速度；ε，δ和γ是表示HTI介質(zhì)各向異性強(qiáng)度的3個無量綱因子.根據(jù)式（1）－（4）可得裂縫參數(shù)與各向異性參數(shù)之間的關(guān)系為

將線性滑動理論與Hudson等效理論相結(jié)合，對4組典型的裂縫介質(zhì)進(jìn)行簡單的等效，其結(jié)果如表1所示.裂縫層埋深3 000m，巖層厚度為100m，背景巖石的縱波速度為6 200m／s，橫波速度為3 500m／s，密度為2 800kg／m3.

表1 裂縫介質(zhì)參數(shù)Table 1 Fracture parameters

2 HTI介質(zhì)一階速度應(yīng)力方程

式（2）—（3）可將裂縫參數(shù)等效為表征裂縫介質(zhì)的彈性矩陣參數(shù)，即將裂縫介質(zhì)等效為各向異性介質(zhì).下面推導(dǎo)HTI介質(zhì)的一階速度應(yīng)力方程.由廣義胡克定律、運(yùn)動微分方程及幾何方程可得HTI介質(zhì)的波動方程.

廣義胡克定律

運(yùn)動微分方程

幾何方程

式中，U＝（ux，uy，uz）為地震波場，F(xiàn)＝（fx，fy，fz）為震源，L為偏微分算子，ρ為密度，

根據(jù)式（6）—（8）可以求得HTI介質(zhì)一階速度應(yīng)力方程的表達(dá)式為

3 裂縫介質(zhì)旋轉(zhuǎn)交錯差分方程

3.1 時間2階差分方程近似

3.2 空間高階差分方程近似

式中，Δx為x方向的采樣間隔，cn為有限差分系數(shù)，具體如表2所示.

表2 幾種常用有限差分系數(shù)Table 2 Several kinds of finite difference coefficients

3.3 一階速度應(yīng)力差分方程

交錯網(wǎng)格主體網(wǎng)格定義的變量與交錯網(wǎng)格半定義的變量不同，速度、應(yīng)力和彈性參數(shù)的定義規(guī)則如圖2所示.其中變量分別定義在網(wǎng)格1，2，3，4點(diǎn)上.

在普通交錯網(wǎng)格有限差分中，通常將介質(zhì)參數(shù)定義在法向應(yīng)力所在的網(wǎng)格點(diǎn)上.本文中法向應(yīng)力及所有的背景參數(shù)都定義在主體網(wǎng)格1點(diǎn)上，z方向速度定義在2點(diǎn)上，x方向速度定義在3點(diǎn)上，切應(yīng)力定義在4點(diǎn)上.應(yīng)用普通交錯網(wǎng)格需要定義1套主網(wǎng)格點(diǎn)和3套半網(wǎng)格點(diǎn)，在模擬裂縫介質(zhì)傳播時需要內(nèi)插相關(guān)的彈性參數(shù).應(yīng)力速度及背景參數(shù)的定義規(guī)則如表3所示.其中：正應(yīng)力τxx，τzz以及彈性參數(shù)c11，c13，c55，c33，ρ定義在主體網(wǎng)格1點(diǎn)上；切應(yīng)力τxz定義在4點(diǎn)所示的半網(wǎng)格點(diǎn)上，z方向速度vz定義在2點(diǎn)所示的半網(wǎng)格點(diǎn)上；x方向速度vx定義在3點(diǎn)所示的半網(wǎng)格點(diǎn)上；半網(wǎng)格3點(diǎn)的密度需要用密度差值移位或者取平均得到.

圖2 標(biāo)準(zhǔn)交錯網(wǎng)格示意圖Fig.2 The sketch of standard staggered grid

圖3 旋轉(zhuǎn)交錯網(wǎng)格示意圖Fig.3 The sketch of rotated staggering grid

表3 標(biāo)準(zhǔn)交錯網(wǎng)格中彈性波場分量和彈性參數(shù)的空間位置Table 3 The positions of elastic wavefield components and elastic parameters in standard staggered grid

表4 旋轉(zhuǎn)交錯網(wǎng)格中彈性波場分量和彈性參數(shù)的空間位置Table 4 The positions of elastic wavefield components and elastic parameters in rotated staggering grid

當(dāng)介質(zhì)非均質(zhì)非常強(qiáng)時，普通交錯網(wǎng)格有限差分得不到滿意的結(jié)果，此時需要用旋轉(zhuǎn)交錯網(wǎng)格模擬.旋轉(zhuǎn)交錯網(wǎng)格首先將圖2沿對角線旋轉(zhuǎn)，使得所有的應(yīng)力分量位于同一點(diǎn)上，所有的速度位于同一點(diǎn)上.本文的變量定義規(guī)則如圖3所示.每個點(diǎn)代表的變量含義如表4所示.

旋轉(zhuǎn)交錯網(wǎng)格中所有速度分量和應(yīng)力分量均位于相同的網(wǎng)格點(diǎn)上，因此介質(zhì)密度與彈性參數(shù)均位于速度所在的網(wǎng)格點(diǎn)或者應(yīng)力所在的網(wǎng)格點(diǎn)上.旋轉(zhuǎn)交錯網(wǎng)格中只需要定義一套主體網(wǎng)格和一套半網(wǎng)格，在模擬裂縫介質(zhì)地震波場傳播時不需要對彈性參數(shù)進(jìn)行內(nèi)插，從而提高了穩(wěn)定性.旋轉(zhuǎn)交錯網(wǎng)格中速度應(yīng)力參數(shù)及背景參數(shù)如表4所示.

表4給出了應(yīng)力速度及背景參數(shù)的定義規(guī)則.其中正應(yīng)力τxx和τzz、切應(yīng)力τxz以及彈性參數(shù)c11，c13，c55，c33，ρ均定義在主體網(wǎng)格1點(diǎn)上，x方向速度vx與z方向速度vz定義在2點(diǎn)所示的半網(wǎng)格點(diǎn)上.旋轉(zhuǎn)交錯網(wǎng)格將相同的物理量定義在相同的網(wǎng)格點(diǎn)上，計(jì)算其沿對角線的差分，再利用對角線差分組合計(jì)算其沿坐標(biāo)軸的差分.此方法不需要對介質(zhì)參數(shù)進(jìn)行平均或者差值.下面推導(dǎo)二維情況下旋轉(zhuǎn)交錯網(wǎng)格沿坐標(biāo)軸方向的表達(dá)式.

根據(jù)矢量合成有

將式（12）、（13）帶入式（9），得到旋轉(zhuǎn)交錯網(wǎng)格下裂縫介質(zhì)的一階速度應(yīng)力方程與一階應(yīng)力速度方程為

式（14）中只需要用到對角線方向的差分格式，將對角線的差分格式進(jìn)行線性組合即可得到最終的結(jié)果.將式（11）帶入式（14）得

式（15）即為變量u對?x，?z的差分，將變量u換成速度分量和應(yīng)力分量便可得到式（14）的差分格式，進(jìn)行正演模擬即可得到裂縫介質(zhì)的波場響應(yīng)特征.

3.4 正演模擬中需要注意的問題

3.4.1 穩(wěn)定性條件

對于裂縫介質(zhì)一階速度應(yīng)力方程與一階應(yīng)力速度差分方程而言，得到其穩(wěn)定性的解析式比較困難.通過矩陣特征分析（Igel et al，1995）可得到其穩(wěn)定性的一種表達(dá)形式.通過矩陣特征分析法得到裂縫介質(zhì)一階速度應(yīng)力彈性波方程時間（2階）、空間（10階）的穩(wěn)定性條件為

式中，Δt為時間采樣間隔，Cmax為彈性矩陣最大值，ρ為密度，Δx和Δz為空間采樣間隔，ci為差分系數(shù)（表2）.若選取參數(shù)不正確，程序不穩(wěn)定，會使得正演模擬結(jié)果錯誤，消除的方法是根據(jù)模型參數(shù)通過式（16）選取合適的時間空間采樣間隔.須說明，時間采樣點(diǎn)不是越小越好，時間采樣點(diǎn)越小則需要的時間點(diǎn)數(shù)越多，累計(jì)誤差越大，也會引起不穩(wěn)定且計(jì)算時間太長，不利于實(shí)際研究.

3.4.2 邊界條件

采用有限差分進(jìn)行裂縫介質(zhì)正演模擬的過程中，由于模擬區(qū)域不是無限大，這種人為限定模型計(jì)算區(qū)域的方法會引起邊界問題.這種邊界引起的反射波常常比物理邊界更強(qiáng)，在地震記錄中會對有效信息產(chǎn)生干擾，因此針對邊界問題必須采用一定的邊界條件消除這種邊界反射（何燕，2008）.消除邊界反射的方法有很多，其中具有代表性的有：Clayton和Engquist（1977）基于旁軸近似提出了吸收邊界條件，其在特定的入射角和頻率范圍內(nèi)具有較好的吸收效果；Reynolds（1978）提出了透明邊界條件，其特點(diǎn)是零角度入射時反射系數(shù)為零；Marfurt（1984）提出了海綿邊界條件；Berenger（1994）針對電磁波傳播情況提出了最佳匹配層（perfectly matcher layer，簡稱PML）吸收邊界條件，其在海綿吸收邊界的基礎(chǔ)上作了進(jìn)一步改進(jìn)，并從理論上證明了該方法可以完全吸收來自各個方向的各種頻率波，且不產(chǎn)生任何反射，是目前應(yīng)用最廣泛的邊界條件.本文采用PML邊界條件消除邊界反射.該方法的原理是在研究區(qū)域的邊界上增加一個吸收層，使邊界上傳入吸收層的波隨傳播距離的增加成指數(shù)衰減，而不產(chǎn)生任何邊界反射，可以得到裂縫介質(zhì)一階速度應(yīng)力帶有PML邊界的差分格式.旋轉(zhuǎn)交錯網(wǎng)格的PML邊界條件與普通交錯網(wǎng)格的PML邊界條件一致.下面以普通交錯網(wǎng)格方程推導(dǎo)PML邊界條件，然后再應(yīng)用于旋轉(zhuǎn)交錯網(wǎng)格.以公式（9）中其中一個方程（省略震源）為例，

根據(jù)復(fù)數(shù)坐標(biāo)變換

對式（17）兩邊進(jìn)行傅里葉變換得

根據(jù)波場分解原理，式（19）可分解為

將式（18）帶入式（20），并進(jìn)行傅里葉反變換得

式（9）中其余方程的推導(dǎo)過程如上，不再一一重復(fù).

3.4.3 網(wǎng)格頻散

采用有限差分法存在固有缺陷，即在利用有限差分算子逼近微分算子時會產(chǎn)生誤差項(xiàng)，它使得具有不同頻率的地震波表現(xiàn)為不同的相速度，從而使地震波發(fā)生頻散，影響數(shù)值模擬的精度和偏移成像效果.數(shù)值頻散實(shí)質(zhì)上是一種因離散化求解波動方程而產(chǎn)生的偽波動，它是指組成地震波的不同頻率分量以不同的速度傳播，隨著傳播時間的推移地震波擴(kuò)展為更長波列的現(xiàn)象（孫成禹等，2009）.時間網(wǎng)格頻散表現(xiàn)為相速度超前，空間網(wǎng)格頻散表現(xiàn)為相速度滯后.本文采用通量傳輸校正方法（Book et al，1975；Boris，Book，1973）消除頻散，因?yàn)槭褂每臻g高階有限差分方法會壓制頻散且計(jì)算量非常大.通量傳輸校正方法主要分為以下幾步：

1）輸運(yùn)計(jì)算，計(jì)算待修正值

2）求n時刻的擴(kuò)散通量

式中η1為擴(kuò)散系數(shù).

3）擴(kuò)散計(jì)算

4）反擴(kuò)散通量計(jì)算

式中η2為反擴(kuò)散系數(shù).

5）限制反擴(kuò)散通量條件

6）反擴(kuò)散計(jì)算

4 裂縫介質(zhì)正演模擬分析

采用式（14）、（15）的差分方程，編寫程序進(jìn)行裂縫介質(zhì)正演模擬，可以得到裂縫介質(zhì)的正演模擬記錄.其中模擬過程中需要注意邊界條件、頻散壓制及穩(wěn)定性條件.裂縫引起的反射波能量相對于背景介質(zhì)較小，因而正演模擬中裂縫介質(zhì)引起的反射系數(shù)變化經(jīng)常被淹沒在背景介質(zhì)中.為了研究裂縫介質(zhì)的波場，本文設(shè)計(jì)一個均勻模型裂縫位于均勻模型中間，如圖4所示.模型尺寸為4 000m×4 000 m，縱向、橫向網(wǎng)格間距均為10m，縱波速度為3 000m／s，橫波速度為1 700m／s，密度為2 000 kg／m3.黑色區(qū)域?yàn)榱芽p所在的區(qū)域，其中：裂縫長度為5m，寬度為1mm，密度為2條／m，裂縫充填物為油；裂縫橫向延伸為4 000m、縱向延伸為100m.采用式（14）、（15）的差分方程進(jìn)行正演模擬.正演模擬采用雷克子波，子波的主頻為25Hz，時間采樣率為1ms，時間長度為5s.正演模擬結(jié)果如圖5所示.

圖4 單層模型示意圖Fig.4 The sketch of single layer model

圖5 中由于裂縫引起的各向異性較小，qSV波和qSH波在地震波場中無法區(qū)分，統(tǒng)一用qS波表示.從圖5a中可以看出，t＝0.75s時波還沒有傳播到裂縫處，地震波場中含有兩種波，即qP波和qS波.其中qP和qS分別為準(zhǔn)P波和準(zhǔn)S波.從圖5b可以看出，t＝1.0s時qP波已經(jīng)穿過裂縫層，地震波場中含有qS波、qPRqP波（qP波入射qP反射波）、qPRqS波（qP波入射qS反射波）、qPTqS波（qP波入射qS透射波）及qPTqP波（qP波入射qP透射波）.地震波場比較復(fù)雜，可以看出裂縫介質(zhì)頂?shù)锥紩a(chǎn)生反射波、透射波.從圖5c中可以看出，t＝1.6s時qS波已經(jīng)穿過裂縫層，地震記錄中含有qPRqP波（qP波入射qP反射波）、qPRqS波（qP波入射qS反射波）、qSRqP波（qS波入射qP反射波）、qSRqS波（qS波入射qS反射波）、qPTqS波（qP波入射qS透射波）、qPTqP波（qP波入射qP透射波）、qSTqS波（qS波入射qS透射波）及qSTqP波（qS波入射qP透射波）.由此可見，t＝1.6s時地震波場非常復(fù)雜，裂縫介質(zhì)頂?shù)锥紩a(chǎn)生反射波、透射波.地震波穿過裂縫介質(zhì)時地震波場非常復(fù)雜，便于我們利用地震波場的信息反演裂縫的參數(shù).

圖5 不同時刻地震波場示意圖.（a）t＝0.75s；（b）t＝1.0s；（c）t＝1.6sFig.5 The sketch of seismic wavefield in different moments.（a）t＝0.75s；（b）t＝1.0s；（c）t＝1.6s

同樣采用圖4模型示意圖，模型尺寸改為400m×400m，模型采樣間隔為10m×10m.采用主頻為75Hz的雷克子波，時間采樣率為0.1ms，時間長度為0.5s，炮點(diǎn)位于模型中間.模型縱波速度為3 000m／s，橫波速度為1 700m／s.黑色區(qū)域?yàn)榱芽p所在的區(qū)域，其中：裂縫長度為5m，寬度為1mm，密度為2條／m，裂縫充填物為油；裂縫橫向延伸為400m，縱向延伸為1—24m.從地震記錄中提取第50道的PP波和SS波記錄，如圖6所示.

圖6 不同裂縫延伸長度的單道地震記錄.（a）PP波；（b）SS波Fig.6 Seismic records of single trace in different crack growth.（a）PP wave；（b）SS wave

根據(jù)正演采用的參數(shù)，縱波的波長為40m，橫波的波長為23m.由圖6a可知，當(dāng)裂縫延伸為縱波半個波長時，頂?shù)追瓷渫耆珠_；當(dāng)裂縫延伸為縱波1／4波長時，可以找出頂?shù)追瓷浣缑妫藭r頂?shù)追瓷洳ㄉ形赐耆珠_.由圖6b可知，當(dāng)裂縫延伸為橫波半個波長時，頂?shù)追瓷渫耆珠_；當(dāng)裂縫延伸為橫波1／4波長時，可以找出頂?shù)追瓷浣缑?，此時頂?shù)追瓷洳ㄉ形赐耆珠_.

為了更清楚地了解裂縫介質(zhì)對地震波場的影響，從不同的裂縫參數(shù)出發(fā)研究不同裂縫參數(shù)對地震波場的影響.下面分別研究裂縫密度、裂縫縱橫比、裂縫充填物等裂縫參數(shù)對地震波場的影響.

4.1 裂縫密度對地震波場的影響

采用圖4所示的模型示意圖.裂縫長度為2m，寬度為1mm，裂縫充填物為油；裂縫密度分別為0.5，1，2，4和5條／m；采用雷克子波進(jìn)行正演模擬，子波的主頻為25Hz，時間采樣率為1ms，時間長度為5s.正演模擬結(jié)果如圖7所示.可以看出，裂縫密度越大，引起的反射波能量越強(qiáng).

圖7 不同裂縫密度（e）正演模擬的地震記錄（a）e＝0.5條／m；（b）e＝1條／m；（c）e＝2條／m；（d）e＝4條／m；（e）e＝5條／mFig.7 Seismic records in the media with different fracture densities（a）e＝0.5m－1；（b）e＝1m－1；（c）e＝2m－1；（d）e＝4m－1；（e）e＝5m－1

4.2 裂縫縱橫比對地震波場的影響

采用圖4所示的模型示意圖.裂縫密度為1條／m，裂縫充填物為油，裂縫縱橫比分別為1 000，2 000，5 000，7 000和10 000；采用雷克子波進(jìn)行正演模擬，子波的主頻為25Hz，時間采樣率為1ms，時間長度為5s.正演模擬結(jié)果如圖8所示.可以看出，裂縫縱橫比越小，引起的反射波能量越強(qiáng).

4.3 裂縫充填物對地震波場的影響

采用圖4所示的模型示意圖.裂縫密度為1條／m，長度為2m，寬度為1mm，裂縫充填物分別為油、水和泥巖；采用雷克子波進(jìn)行正演模擬，子波的主頻為25Hz，時間采樣率為1ms，時間長度為5s.正演模擬結(jié)果如圖9所示.可以看出，裂縫充填物速度與背景速度差別越大，引起的反射波能量越強(qiáng).

圖8 不同裂縫縱橫比（a／c）正演模擬的地震記錄（a）a／c＝1000；（b）a／c＝2000；（c）a／c＝5000；（d）a／c＝7000；（e）a／c＝10000Fig.8 Seismic records in the media with different fracture aspect ratios（a）a／c＝1000；（b）a／c＝2000；（c）a／c＝5000；（d）a／c＝7000；（e）a／c＝10000

圖9 不同裂縫充填物正演模擬的地震記錄（a）油；（b）水；（c）泥巖Fig.9 Seismic records in the media with different fracture fillings（a）Oil；（b）Water；（c）Mudstone

5 討論與結(jié)論

本文提出一種高陡傾角裂縫介質(zhì)的正演模擬方法.利用裂縫等效理論將高陡傾角裂縫介質(zhì)等效為橫向各向同性（HTI）介質(zhì)，推導(dǎo)出HTI介質(zhì)旋轉(zhuǎn)交錯網(wǎng)格一階速度應(yīng)力方程及所對應(yīng)的差分方程，實(shí)現(xiàn)了裂縫介質(zhì)的正演模擬，分析了裂縫介質(zhì)中的地震波場響應(yīng)特征，討論了不同裂縫密度、裂縫縱橫比和裂縫充填物的地震響應(yīng)特征.結(jié)果表明：裂縫的存在相當(dāng)于人為增加反射界面；裂縫密度越大，裂縫縱橫比越小，裂縫充填物與背景介質(zhì)彈性性質(zhì)差別越大，引起的反射波能量變化越大.本文模擬結(jié)果為利用地震數(shù)據(jù)進(jìn)行裂縫介質(zhì)參數(shù)反演與儲層識別及預(yù)測提供了依據(jù).需要說明的是：本文結(jié)果僅僅針對一組高陡傾角裂縫介質(zhì)，而實(shí)際地層中裂縫的存在形式是多種多樣的，需進(jìn)一步研究實(shí)際裂縫介質(zhì)中的地震波場特征；本文僅僅對裂縫充填物、裂縫密度、裂縫縱橫比等裂縫物性參數(shù)進(jìn)行了正演模擬，需要進(jìn)一步研究裂縫介質(zhì)其它物性參數(shù)的地震響應(yīng)特征；本文使用彈性波進(jìn)行正演模擬，而實(shí)際地震記錄通常是縱波記錄，因此需要進(jìn)一步研究各向異性介質(zhì)中縱橫波分離的問題.

刁順，楊慧珠，許云.1994a.裂縫及介質(zhì)非均勻散射［J］.石油物探，33（4）：49-55.

Diao S，Yang H Z，Xu Y.1994a.Scattering due to inhomogeneity of media［J］.Geophysical Prospecting for Petroleum，33（4）：49-55（in Chinese）.

刁順，楊慧珠，許云.1994b.TI介質(zhì)裂紋散射研究［J］.石油地球物理勘探，29（5）：545-551.

Diao S，Yang H Z，Xu Y.1994b.Research on crack scattering in transversally isotropic medium［J］.Oil Geophysical Prospecting，29（5）：545-551（in Chinese）.

董良國，馬在田，曹景忠，王華忠，耿建華，雷冰，許世勇.2000.一階彈性波方程交錯網(wǎng)格高階差分解法［J］.地球物理學(xué)報(bào)，43（3）：411-419.

Dong L G，Ma Z T，Cao J Z，Wang H Z，Geng J H，Lei B，Xu S Y.2000.A staggered-grid high-order difference method of one-order elastic wave equation［J］.Chinese Journal of Geophysics，43（3）：411-419（in Chinese）.

郭建.1993.裂隙介質(zhì)中的各向異性研究［J］.石油地球物理勘探，28（3）：348-353.

Guo J.1993.A research into anisotropy of fractured medium［J］.Oil Geophysical Prospecting，28（3）：348-353（in Chinese）.

何燕.2008.正交各向異性彈性波高階有限差分正演模擬研究［D］.東營：中國石油大學(xué)（華東）地球資源與信息學(xué)院：55-58.

He Y.2008.High-Order Finite-Difference Forward Modeling of Elastic-Wave in Orthorhombic Anisotropic Media［D］.Dongying：College of Geo-resources and Information，China University of Petroleum：55-58（in Chinese）.

金抒辛，何樵登.2005.泥巖裂隙的正演模擬方法研究［J］.石油物探，44（2）：119-127.

Jin S X，He Q D.2005.Forward modeling method research on mudstone fractures medium［J］.Geophysical prospecting for Petroleum，44（2）：119-127（in Chinese）.

孔麗云，王一博，楊慧珠.2012.裂縫誘導(dǎo)HTI雙孔隙介質(zhì)中的裂縫參數(shù)分析［J］.地球物理學(xué)報(bào)，55（1）：189-196.

Kong L Y，Wang Y B，Yang H Z.2012.Fracture parameters analyses in fracture-induced HTI double-porosity medium［J］.Chinese Journal of Geophysics，55（1）：189-196（in Chinese）.

李雪靜，劉定進(jìn).2008.各向同性介質(zhì)彈性波交錯網(wǎng)格有限差分正演模擬［J］.小型油氣藏，13（3）：15-20.

Li X J，Liu D J.2008.Isotropic medium elastic wave cross grid finite difference forward simulation［J］.Small Hydrocarbon Reservoirs，13（3）：15-20（in Chinese）.

毛寧波，謝濤，楊凱，金明霞.2008.裂縫儲層地震方位AVO正演模擬研究及應(yīng)用［J］.石油天然氣學(xué)報(bào)，30（5）：59-63.

Mao N B，Xie T，Yang K，Jin M X.2008.Azimuthal AVO forward model for fracture reservoirs and its application［J］.Journal of Oil and Gas Technology，30（5）：59-63（in Chinese）.

潘保芝，張麗華，單剛義，楊雪.2006.裂縫和孔洞型儲層孔隙模型的理論進(jìn)展［J］.地球物理學(xué)進(jìn)展，21（4）：1232-1237.

Pan B Z，Zhang L H，Shan G Y，Yang X.2006.Progress in porosity model for fractured and vuggy reservoirs［J］.Progress in Geophysics，21（4）：1232-1237（in Chinese）.

裴正林.2004.三維各向異性介質(zhì)中彈性波方程交錯網(wǎng)格高階有限差分法數(shù)值模擬［J］.石油大學(xué)學(xué)報(bào)：自然科學(xué)版，28（5）：23-29.

Pei Z L.2004.Three-dimensional numerical simulation of elastic wave propagation in 3-D anisotropic media with staggered-grid high-order difference method［J］.Journal of China University of Petroleum：Edition of Natural Science，28（5）：23-29（in Chinese）.

裴正林，王尚旭.2005.任意傾斜各向異性介質(zhì)中彈性波波場交錯網(wǎng)格高階有限差分法模擬［J］.地震學(xué)報(bào)，27（4）：441-451.

Pei Z L，Wang S X.2005.A staggered-grid high-order finite-difference modeling for elastic wave field in arbitrary tilt anisotropic media［J］.Acta Seismologica Sinica，27（4）：441-451（in Chinese）.

裴正林.2006.雙相各向異性介質(zhì)彈性波傳播交錯網(wǎng)格高階有限差分法模擬［J］.石油地球物理勘探，41（2）：137-143.

Pei Z L.2006.Staggered grid high-order finite-difference modeling of elastic wave transmission in biphase anisotropic medium［J］.Oil Geophysical Prospecting，41（2）：137-143（in Chinese）.

孫成禹，宮同舉，張玉亮，張文穎.2009.波動方程有限差分法中的頻散與假頻分析［J］.石油地球物理勘探，44（1）：43-48.

Sun C Y，Gong T J，Zhang Y L，Zhang W Y.2009.Analysis on dispersion and alias in finite-difference solution of wave equation［J］.Oil Geophysical Prospecting，44（1）：43-48（in Chinese）.

謝桂生.1994.裂隙裂縫介質(zhì)的彈性波模擬［J］.石油地球物理勘探，29（5）：537-544.

Xie G S.1994.Simulation of elastic waves in fractured medium［J］.Oil Geophysical Prospecting，29（5）：537-544（in Chinese）.

嚴(yán)洪勇，劉洋.2012.黏彈TTI介質(zhì)中旋轉(zhuǎn)交錯網(wǎng)格高階有限差分?jǐn)?shù)值模擬［J］.地球物理學(xué)報(bào)，55（4）：1354-1365.

Yan H Y，Liu Y.2012.Rotated staggered grid high-order finite-difference numerical modeling for wave propagation in viscoelastic TTI media［J］.Chinese Journal of Geophysics，55（4）：1354-1365（in Chinese）.

楊曉，王真理，喻岳鈺.2010.裂縫型儲層地震檢測方法綜述［J］.地球物理學(xué)進(jìn)展，25（5）：1785-1794.

Yang X，Wang Z L，Yu Y Y.2010.The overview of seismic techniques in prediction of fracture reservoir［J］.Progress in Geophysics，25（5）：1785-1794（in Chinese）.

張中杰，何樵登.1989.含裂隙介質(zhì)中地震波運(yùn)動學(xué)問題的正演模擬［J］.石油地球物理勘探，24（3）：290-299.

Zhang Z J，He Q D.1989.Forward modeling of kinematic problems of seismic wave in fractured medium［J］.Oil Geophysical Prospecting，24（3）：290-299（in Chinese）.

Bakulin A，Grechka V，Tsvankin L.2000a.Estimation of fracture parameters from reflection seismic data.Part 1：HTI model due to a single fracture set［J］.Geophysics，65（6）：1788-1802.

Bakulin A，Grechka V，Tsvankin L.2000b.Estimation of fracture parameters from reflection seismic data.Part 2：Fractured models with orthorhombic symmetry［J］.Geophysics，65（6）：1803-1817.

Bakulin A，Grechka V，Tsvankin L.2000c.Estimation of fracture parameters from reflection seismic data.Part 3：Fractured models with monoclinic symmetry［J］.Geophysics，65（6）：1818-1830.

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

Boris J P，Book D L.1973.Flux-corrected transport.Ⅰ.SHASTA，a fluid transport algorithm that works［J］.J Computat Phys，11（1）：38-69.

Book D L，Boris J P，Hain K.1975.Flux-corrected transportⅡ：Generalizations of the method［J］.J Computat Phys，18（3）：248-283.

Clayton R，Engquist B.1977.Absorbing boundary conditions for acoustic and elastic wave equations［J］.Bull Seismol Soc Am，67（6）：1529-1540.

Crampin S.1985.Evaluation of anisotropy by shear-wave splitting［J］.Geophysics，50（1）：142-152.

Hudson J A.1981.Wave speeds and attenuation of elastic waves in material containing cracks［J］.Geophysics，64（1）：133-150.

Igel H，Mora P，Riollet B.1995.Anisotropic wave propagation through finite-difference grids［J］.Geophysics，60（4）：1203-1216.

Marfurt K J.1984.Accuracy of finite-difference and finite-element modeling of the scalar and elastic wave equations［J］.Geophysics，49（5）：533-549.

Reynolds A C.1978.Boundary conditions for the numerical solution of wave propagation problems［J］.Geophysics，43（6）：1099-1110.

Saenger E H，Gold N，Shapiro S A.2000.Modeling the propagation of elastic waves using a modified finite-difference grid［J］.Wave Motion，31（1）：77-92.

Schoenberg M，Sayers C M.1995.Seismic anisotropy of fractured rock［J］.Geophysics，60（1）：204-221.

Thomsen L.1986.Weak elastic anisotropy［J］.Geophysics，51（10）：1954-1966.