侯 凱， 李 黎， 羅 澤， 楊小江， 趙紅娟， 趙偉超

(中海石油(中國)有限公司 深圳分公司研究院，深圳 南山 518000)

0 引言

裂縫性油氣藏作為一種非常復雜的隱蔽油氣藏，廣泛發(fā)育在礁灰?guī)r、泥頁巖以及火成巖中。據(jù)統(tǒng)計，世界上已探明的總儲量的一半為裂縫性油氣藏[1]。裂縫性油氣藏研究最突出的問題是裂縫的尺度非常小，遠遠低于其他地質(zhì)體的地震響應，地震采集的數(shù)據(jù)中雖然包含裂縫型儲層的響應信息，但是裂縫能量的響應很弱，識別具有很大的困難。在20世紀80年代以前，裂縫性儲層的研究依據(jù)野外露頭、巖芯和測井資料等，從構(gòu)造地質(zhì)學角度定性研究裂縫的類型、成因機制和分布特征。后來根據(jù)應力場數(shù)值模擬和巖石破裂準則，進行半定量研究。隨著計算機技術(shù)的快速發(fā)展，裂縫研究逐漸定量化[2]。Wu Huang等[3]利用De Wolf近似和相位屏算子計算三維背景下的背向散射場，符力耘等[4]建立了均勻背景介質(zhì)的Green函數(shù)基本解體積分方程，尹軍杰等[11]對于大量的地質(zhì)模型進行散射波數(shù)值模擬，系統(tǒng)分析了散射波與反射波的差異，蕾蕾[5]根據(jù)散射積分算法，對于縫洞模型進行了可變網(wǎng)格的積分模擬，分析了散射波長的傳播特征。裂縫儲層的正演模擬雖然取得了很大地進步，但是仍然無法精確刻畫地下厘米和毫米級別的裂縫響應特征。

本次研究從散射理論出發(fā)，采用多級可變網(wǎng)格的正演模擬方法對于裂縫進行了大量的模擬，實現(xiàn)了厘米甚至毫米級別裂縫的地震正演模擬，并統(tǒng)計和定量分析了裂縫群的發(fā)育密度、速度、傾角和分布范圍等不同時，地震波場的傳播特征，準確地找出裂縫的各項參數(shù)與地震接收記錄的波場對應關(guān)系，大大提高了裂縫模擬的計算效率，為裂縫型油氣藏的準確預測提供了有力根據(jù)。

1 地震散射理論

散射理論最早是在量子力學中提出的，定義為由傳播介質(zhì)的不均勻性引起的光線向四周射去的現(xiàn)象。假設介質(zhì)是完全均勻的，地震波從震源出發(fā)將按照球面擴散的形式以不變的速度傳播。若介質(zhì)中存在與均勻背景不一樣的非均質(zhì)體時，地震波的傳播發(fā)生改變，此時的非均質(zhì)體可以看作二次震源，我們稱這種現(xiàn)象為散射現(xiàn)象[6-8]。根據(jù)散射理論，我們將介質(zhì)分為背景介質(zhì)和擾動介質(zhì)，在背景介質(zhì)中傳播的波場，稱之為背景波場。如果介質(zhì)的性質(zhì)發(fā)生變化，那么這部分介質(zhì)稱為擾動介質(zhì)。由介質(zhì)的擾動部分產(chǎn)生的波場稱為散射波場。

地震波在均勻介質(zhì)中聲波的傳播符合以下規(guī)律：

(1)

其中：C(r)為均勻介質(zhì)中的傳播速度；u0(r,t)為產(chǎn)生的波場，將式(1)轉(zhuǎn)換到頻率域得到式(2)。

(2)

在索莫菲輻射邊界條件的假設下：

r→∞,r=|r|

(3)

ε(r)為相對于C0(r)的擾動波場，并滿足C0(r)橫向不變，C0(r)可以表示為C0(z)。

(4)

根據(jù)擾動理論，地震波場u可以分解為均勻介質(zhì)中的波場u0與擾動介質(zhì)產(chǎn)生的波場的疊加，即：

u(x,y,z;ω)=u0(x,y,z;ω)+U(x,y,z;ω)

(5)

均勻介質(zhì)中的波場可以表示成式(6)。

(6)

將式(5)、式(6)帶入到式(2)中得式(7)。

u(x,y,z;ω)

(7)

用格林函數(shù)法求解求解方程(7)，最終得到波場的積分表達行式如下：

u(x,y,z;ω)=u0(x,y,z;ω)+us(x,y,z;ω)=

u0(x,y,z;ω)+

(8)

方程(8)就是Lippmann-Swinger方程[9]，此方程是一個非線性方程，在實際中需要使用born公式簡化。

2 多級可變網(wǎng)格

常規(guī)的有限差分正演模擬方法考慮到模擬的精度，又要考慮計算的效率，空間采樣間隔通常是5 m～10 m。Moczo[16]首次提出了可變網(wǎng)格的思想，根據(jù)非均勻體的尺度，對空間步長調(diào)整，可以在一定程度上，保證模擬效率的情況下提高精度。然而對裂縫儲層來說，仍然不能準確刻畫，為了克服可變網(wǎng)格倍數(shù)的限制，筆者首次提出并采用了多級變網(wǎng)格的方法對于裂縫進行模擬，成功的將變網(wǎng)格區(qū)域繼續(xù)放大。

如圖1所示，常規(guī)網(wǎng)格無法對于小尺度模型準確刻畫，經(jīng)過可變網(wǎng)格加密后，可以提高大尺度溶洞的刻畫程度，但仍然無法準確地模擬裂縫。在這種情況下，在可變網(wǎng)格的基礎上，對于微裂縫繼續(xù)空間加密，經(jīng)過兩次的加密可以達到了模型放大100倍的效果，準確地刻畫微裂縫，這種局部加密的方法提高了模擬的精度，又保證了正演的效率。

圖1 多級變網(wǎng)格的模型剖分示意圖Fig.1 The schematic diagram of multistage variable mesh model

圖2 過渡區(qū)域空間導數(shù)網(wǎng)格計算點Fig.2 Spatial derivative grid computing point of transition region

圖3 含裂縫的斷層模型Fig.3 A fractured fault model

1)多級可變網(wǎng)格空間過渡。為了保證變網(wǎng)格區(qū)域的波場穩(wěn)定過渡，需要重新定義過渡區(qū)域網(wǎng)格差分格式，讓加密區(qū)域的波場得到耦合。如圖2所示，變網(wǎng)格區(qū)域參與導數(shù)計算點是不對稱的，為了解決此問題，小網(wǎng)格區(qū)域取點不是等距的。對差分格式和差分系數(shù)做如下修正：

u(x-mNdx,z)-2u(x,z) ]

(9)

u(x,z-mNdz)-2u(x,z) ]

(10)

式(9)和式(10)分別給出了加密區(qū)域x方向和z方向的網(wǎng)格空間二階導數(shù)M階精度的計算公式。dx表示加密網(wǎng)格步長，m表示空間導數(shù)的點到參與計算點的距離。

圖4 含裂縫斷層模型的接收記錄Fig.4 The receiving record of the fractured model(a)單程波的接收記錄及裂縫響應；(b)雙程波的接收記錄及裂縫響應

2)多級變網(wǎng)格時間過渡方法。在使用多級可變網(wǎng)格時，模型的時間采樣步長同樣也發(fā)生了變化。波動方程的空間采用間隔與時間采樣間隔必須滿足穩(wěn)定性條件，才能保證算法是穩(wěn)定性的。當空間采樣間距減小時，時間采樣的間距必須相應的減小。多級變網(wǎng)格的時間采樣步長與空間采樣步長倍數(shù)的變化保持一致。設變網(wǎng)格的倍數(shù)是M時，t時刻到t+1時刻之間的波場值用線性插值算法計算M-1時刻的波場值，在小網(wǎng)格步長內(nèi)計算t時刻到t+MΔt的波場(Δt是小網(wǎng)格的時間步長)，如果大網(wǎng)格的時間步長為ΔT，則Δt=1/M(ΔT)。最后將t+MΔt時刻的波場值回歸到大網(wǎng)格的網(wǎng)格上。若在一階變網(wǎng)格區(qū)域內(nèi)需要進行二級變網(wǎng)格時，需要對于二級變網(wǎng)格的時間步長和空間步長做與同樣方法的改變即可。經(jīng)過測試可知，改變時間采樣步長不會引起明顯的人為反射，因此不做特殊的處理。

2.1 積分法與有限差分微分法的波場特征對比

圖3為一個斷層模型，在斷嵌位置處發(fā)育垂直裂縫群。模型設置的參數(shù)如下：裂縫的寬度為0.5 m，斷開深度為50 m，模型大小設置為600×500，間隔為5 m。激發(fā)主頻30 Hz的雷克子波，時間間隔是0.5 ms。震源位于地面的中心位置?？勺兙W(wǎng)格積分法使用高斯衰減邊界條件，有限差分微分法使用PML邊界條件。

圖4為兩種裂縫模擬方法的接受記錄。從接受記錄情況看，波動方程微分法可以接收到豐富的波場信息[10]，同時也存在裂縫響應弱，容易被其他信息掩蓋的缺點。波動方程積分法得到的地震記錄不含直達波、多次反射波等信息，如圖4所示，相對于雙程波的模擬來說，單程波模擬的裂縫響應更清晰，且噪音成分少，可以直觀地找到裂縫非均質(zhì)體的響應特征。綜合分析，對于含裂縫的儲層使用單程波進行裂縫響應的研究，更能觀測到裂縫的響應，滿足定量分析的要求，同時也能滿足微裂縫模擬對計算機的要求[11-13]。

圖5 模型示意圖Fig.5 A fractured model diagram(a)常規(guī)網(wǎng)格;(b)一級變網(wǎng)格;(c)二級變網(wǎng)格

圖6 地面記錄的接收結(jié)果Fig.6 Seismic record(a)常規(guī)網(wǎng)格的模擬結(jié)果；(b)多極變網(wǎng)格的模擬區(qū)域

圖7 接收記錄的波形分析Fig.7 Wave figure from cmp

圖8 背斜模型Fig.8 The anticline model(a)無裂縫模型;(b)含裂縫模型

2.2 多級可變網(wǎng)格積分法與常規(guī)網(wǎng)格積分法對比

采用同一模型，分別使用常規(guī)網(wǎng)格積分和多級可變網(wǎng)格的方法進行模擬(圖5)。震源子波是雷克子波，主頻為40 Hz，采樣間距是1 ms。為了測試多級可變網(wǎng)格，我們將模型中的第二個矩形溶洞進行二級變網(wǎng)格的加密，第一次選取模型中的虛形框部分。橫向加密10倍得到一級變網(wǎng)格的模型，變網(wǎng)格區(qū)域內(nèi)的空間采樣間隔變?yōu)? m，時間采樣間距為0.1 ms。然后對虛形框部分繼續(xù)加密10倍，得到二級變網(wǎng)格模型，這時候的網(wǎng)格間距變?yōu)?.1 m，時間采樣間距為0.01 ms。

對比圖6(a)、圖6(b)可看出，多級變網(wǎng)格的模擬區(qū)域的記錄與常規(guī)網(wǎng)格的記錄一致，多級變網(wǎng)格沒有帶來明顯的邊界反射。抽取如圖7所示位置出的波形圖做對比分析，可以清晰地看出,常規(guī)網(wǎng)格的模擬方法與多級變網(wǎng)格的模擬方法得到波形基本一致?？梢钥吹蕉嗉壙勺兙W(wǎng)格可以實現(xiàn)高倍數(shù)的變網(wǎng)格模擬，而且穩(wěn)定性和精度很高。

2.3 網(wǎng)格積分法與網(wǎng)格微分法計算效率對比

圖8給出的背斜模型，模型的速度在圖中已經(jīng)標出。在圖(b)中發(fā)育一組垂直的裂縫群，裂縫的速度為1 500 m/s。模型的網(wǎng)格大小為120×100，空間采樣間距為10 m，震源使用雷克子波，主頻為30 Hz，時間采樣步長是0.5 ms，采樣點是1 500 個。對這兩個模型分別使用多級變網(wǎng)格積分法和多級變網(wǎng)格微分法正演模擬，兩種方法給定的模型參數(shù)以及震源參數(shù)是完全相同的。

圖9 地面接收記錄Fig.9 Seismic record(a)單程波無裂縫模擬結(jié)果;(b)單程波有裂縫模擬結(jié)果;(c)裂縫響應;(d)雙程波無裂縫模擬結(jié)果;(e)雙程波有裂縫模擬結(jié)果;(f)裂縫響應

由圖9可以看出，雙程波模擬結(jié)果中有直達波的成分，而單程波模擬僅含有一次反射波。兩種方法都需要進行高倍數(shù)的變網(wǎng)格來模擬裂縫。多級可變網(wǎng)格采用的10倍的加密，采用微分法正演模擬的計算時間為6 h。采用格林函數(shù)積分法的單程波模擬，進行正演模擬的計算時間為20 min。因此可以得出，對于含有裂縫非均質(zhì)體的模型，格林函數(shù)積分法要比多級變網(wǎng)格的雙程波模擬效率高20倍左右。

3 微裂縫地震響應定量分析

3.1 微裂縫響應特征

通過圖10中給定2個模型分析裂縫的接收記錄響應特征，給定裂縫的寬度為5 mm，傾角為60°，裂縫間的橫向間距為0.1 m,圖10(a)中發(fā)育3組裂縫，裂縫的切深為10 m，范圍是35 m。圖10(b)中不發(fā)育裂縫。

通過對比圖10可以看出，裂縫在地層中的能量很弱，裂縫群能量會對界面反射有干涉作用，在地震記錄中以尾波的形式出現(xiàn)。

3.2 震源主頻度裂縫能量影響

圖11模型范圍是3 000 m*3 000 m，空間采樣間隔為5 m，時間采樣率是1 ms，傾斜裂縫發(fā)育的位置是1 000 m，裂縫橫線的斷開大小只有5 mm，裂縫速度為2 000 m/s，裂縫間的距離是0.25 m，角度為60°，按照表1改變震源的主頻，定量研究了裂縫的響應特征。

通過抽取地震接受記錄，對地震波形和能量變化定量分析(圖12)。由圖12可以看出，裂縫對于不同主頻的震源響應不同，主頻越低地震波形越分散。裂縫能量并不是隨著主頻的升高而升高，而是呈現(xiàn)拋物線的變化特點，當裂縫角度為60°時，主頻為50 Hz的裂縫能量最強，然后隨著主頻的升高而能量降低。

3.3 裂縫發(fā)育密度與能量響應的定量分析

裂縫發(fā)育密度越大，裂縫響應的能量越強，裂縫發(fā)育密度與能量間的響應存在一定的規(guī)律，通過設定模型參數(shù)來定量分析。模型的大小為2 500 m*2 500 m，空間采樣間隔5 m，時間采樣間隔為1 ms，給定裂縫的尺度為0.005 m，裂縫速度為2 000 m/s，縱向切深為10 m，裂縫的傾角為60°。震源主頻為50 Hz。按照表2改變裂縫的發(fā)育密度，獲得裂縫的正演地震記錄。

圖10 微裂縫地震響應Fig.10 The seismic record of micro-fractuew (a)模型1及對應的炮記錄；(b)模型2及對應的炮記錄；(c)中心道炮記錄波形圖

圖11 傾斜模型及對應炮記錄Fig.11 Tilt model and corresponding seismic record

20304050607080

表2 裂縫間距(單位：m)

圖14為正演地震記錄，分析得到裂縫密度越大，裂縫的響應也越強。圖15為裂縫能量與裂縫密度的統(tǒng)計關(guān)系，可以看出裂縫發(fā)育密度與裂縫能量響應之間并不是單純的線性關(guān)系，而是呈現(xiàn)一定的對數(shù)關(guān)系，通過分析裂縫的能量，反過來可以得到裂縫的發(fā)育密度，對于裂縫地識別和定量描述有著十分重要的意義。

3.4 實際模型測試

經(jīng)過理論模型的驗證，該算法可以快速模型微裂縫的地震特征。圖16為某工區(qū)低幅度構(gòu)造的速度模型，模型大小為2 000 m*2 500 m，在模型深度1 850 m～1 870 m處發(fā)育有6組裂縫群，裂縫斷開大小為0.005 m，裂縫速度2 000 m/s，裂縫間的橫向距離是0.1 m，裂縫斷開深度為10 m。震源為30 Hz的雷克子波，分別在模型的左邊、中間和右邊放炮，圖16(b)地震接收到的信號。

圖17對比了有裂縫和無裂縫發(fā)育的地震接收信差別，可以看出，由于裂縫的干涉作用地層界面的反射能量明顯地增強了。同時也驗證了該算法對于毫米級別的裂縫模擬有著很好的適用性。

4 結(jié)論

地震波正演模擬在地震資料采集、處理和解釋過程中具有十分重要的作用，地震波正演的熱點逐漸集中于對復雜微小構(gòu)造的適用性和模擬效率的高效性[14-15]。對于含有裂縫的非常規(guī)油氣藏，由于裂縫的發(fā)育尺度只有厘米和毫米級別，要實現(xiàn)裂縫的準確模擬，常規(guī)的地震正演模擬方法存在模擬效率和精度低、對計算機要求高等不足。

圖12 抽取中心道的波形圖及對應的振幅譜Fig.12 Wave figure from cmp and the corresponding amplitude spectrum

圖13 傾斜模型及對應炮記錄Fig.13 Tilt model and corresponding seismic record

圖15 裂縫密度不同的振幅譜能量統(tǒng)計Fig.15 The amplitude spectrum energy statistics of different fracture density

圖14 抽取中心道的波形圖及對應的振幅譜Fig.14 Wave figure from cmp and the corresponding amplitude spectrum(a)地震波形；(b)頻譜分析的能量譜

圖16 裂縫模型及對應炮記錄Fig.16 Fracture model and corresponding seismic record(a)裂縫模型;(b)地震記錄

圖17 抽取炮記錄中心道的波形圖Fig.17 Wave figure from cmp

為了解決裂縫模擬的問題，筆者基于格林函數(shù)積分法的波動方程，首次提出并采用了多級可變網(wǎng)格的方法，提高模擬效率和準確度，實現(xiàn)了裂縫的準確刻畫。定量化深入研究和總結(jié)了裂縫參數(shù)與裂縫響應特征的對應關(guān)系，分析了裂縫響應特征的變化，對于裂縫的識別和解釋具有重要的意義。