吳國忱 秦海旭

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

引言

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

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

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

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

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

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

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

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

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

圖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)系，但其確定彈性波動(dòng)方程系數(shù)的物理意義很不直觀.Thomsen（1986）提出了一整套表征各向異性的參數(shù)，具體表達(dá)如下：

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

將線性滑動(dòng)理論與Hudson等效理論相結(jié)合，對(duì)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)動(dòng)微分方程及幾何方程可得HTI介質(zhì)的波動(dòng)方程.

廣義胡克定律

運(yùn)動(dòng)微分方程

幾何方程

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

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

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

3.1 時(shí)間2階差分方程近似

3.2 空間高階差分方程近似

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

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

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

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

在普通交錯(cuò)網(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)用普通交錯(cuò)網(wǎng)格需要定義1套主網(wǎng)格點(diǎn)和3套半網(wǎng)格點(diǎn)，在模擬裂縫介質(zhì)傳播時(shí)需要內(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)交錯(cuò)網(wǎng)格示意圖Fig.2 The sketch of standard staggered grid

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

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

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

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

旋轉(zhuǎn)交錯(cuò)網(wǎng)格中所有速度分量和應(yīng)力分量均位于相同的網(wǎng)格點(diǎn)上，因此介質(zhì)密度與彈性參數(shù)均位于速度所在的網(wǎng)格點(diǎn)或者應(yīng)力所在的網(wǎng)格點(diǎn)上.旋轉(zhuǎn)交錯(cuò)網(wǎng)格中只需要定義一套主體網(wǎng)格和一套半網(wǎng)格，在模擬裂縫介質(zhì)地震波場傳播時(shí)不需要對(duì)彈性參數(shù)進(jìn)行內(nèi)插，從而提高了穩(wěn)定性.旋轉(zhuǎn)交錯(cuò)網(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)交錯(cuò)網(wǎng)格將相同的物理量定義在相同的網(wǎng)格點(diǎn)上，計(jì)算其沿對(duì)角線的差分，再利用對(duì)角線差分組合計(jì)算其沿坐標(biāo)軸的差分.此方法不需要對(duì)介質(zhì)參數(shù)進(jìn)行平均或者差值.下面推導(dǎo)二維情況下旋轉(zhuǎn)交錯(cuò)網(wǎng)格沿坐標(biāo)軸方向的表達(dá)式.

根據(jù)矢量合成有

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

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

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

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

3.4.1 穩(wěn)定性條件

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

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

3.4.2 邊界條件

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

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

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

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

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

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

3.4.3 網(wǎng)格頻散

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

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

2）求n時(shí)刻的擴(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)定性條件.裂縫引起的反射波能量相對(duì)于背景介質(zhì)較小，因而正演模擬中裂縫介質(zhì)引起的反射系數(shù)變化經(jīng)常被淹沒在背景介質(zhì)中.為了研究裂縫介質(zhì)的波場，本文設(shè)計(jì)一個(gè)均勻模型裂縫位于均勻模型中間，如圖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，時(shí)間采樣率為1ms，時(shí)間長度為5s.正演模擬結(jié)果如圖5所示.

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

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

圖5 不同時(shí)刻地震波場示意圖.（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的雷克子波，時(shí)間采樣率為0.1ms，時(shí)間長度為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)裂縫延伸為縱波半個(gè)波長時(shí)，頂?shù)追瓷渫耆珠_；當(dāng)裂縫延伸為縱波1／4波長時(shí)，可以找出頂?shù)追瓷浣缑妫藭r(shí)頂?shù)追瓷洳ㄉ形赐耆珠_.由圖6b可知，當(dāng)裂縫延伸為橫波半個(gè)波長時(shí)，頂?shù)追瓷渫耆珠_；當(dāng)裂縫延伸為橫波1／4波長時(shí)，可以找出頂?shù)追瓷浣缑?，此時(shí)頂?shù)追瓷洳ㄉ形赐耆珠_.

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

4.1 裂縫密度對(duì)地震波場的影響

采用圖4所示的模型示意圖.裂縫長度為2m，寬度為1mm，裂縫充填物為油；裂縫密度分別為0.5，1，2，4和5條／m；采用雷克子波進(jìn)行正演模擬，子波的主頻為25Hz，時(shí)間采樣率為1ms，時(shí)間長度為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 裂縫縱橫比對(duì)地震波場的影響

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

4.3 裂縫充填物對(duì)地震波場的影響

采用圖4所示的模型示意圖.裂縫密度為1條／m，長度為2m，寬度為1mm，裂縫充填物分別為油、水和泥巖；采用雷克子波進(jìn)行正演模擬，子波的主頻為25Hz，時(shí)間采樣率為1ms，時(shí)間長度為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)交錯(cuò)網(wǎng)格一階速度應(yīng)力方程及所對(duì)應(yīng)的差分方程，實(shí)現(xiàn)了裂縫介質(zhì)的正演模擬，分析了裂縫介質(zhì)中的地震波場響應(yīng)特征，討論了不同裂縫密度、裂縫縱橫比和裂縫充填物的地震響應(yīng)特征.結(jié)果表明：裂縫的存在相當(dāng)于人為增加反射界面；裂縫密度越大，裂縫縱橫比越小，裂縫充填物與背景介質(zhì)彈性性質(zhì)差別越大，引起的反射波能量變化越大.本文模擬結(jié)果為利用地震數(shù)據(jù)進(jìn)行裂縫介質(zhì)參數(shù)反演與儲(chǔ)層識(shí)別及預(yù)測提供了依據(jù).需要說明的是：本文結(jié)果僅僅針對(duì)一組高陡傾角裂縫介質(zhì)，而實(shí)際地層中裂縫的存在形式是多種多樣的，需進(jìn)一步研究實(shí)際裂縫介質(zhì)中的地震波場特征；本文僅僅對(duì)裂縫充填物、裂縫密度、裂縫縱橫比等裂縫物性參數(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.一階彈性波方程交錯(cuò)網(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ì)彈性波交錯(cuò)網(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.裂縫儲(chǔ)層地震方位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.裂縫和孔洞型儲(chǔ)層孔隙模型的理論進(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ì)中彈性波方程交錯(cuò)網(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ì)中彈性波波場交錯(cuò)網(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ì)彈性波傳播交錯(cuò)網(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.波動(dòng)方程有限差分法中的頻散與假頻分析［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)交錯(cuò)網(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.裂縫型儲(chǔ)層地震檢測方法綜述［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)動(dòng)學(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.