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

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

0 引言

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

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

1 地震散射理論

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

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

(1)

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

(2)

在索莫菲輻射邊界條件的假設(shè)下：

r→∞,r=|r|

(3)

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

(4)

根據(jù)擾動(dòng)理論，地震波場u可以分解為均勻介質(zhì)中的波場u0與擾動(dòng)介質(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)，最終得到波場的積分表達(dá)行式如下：

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

u0(x,y,z;ω)+

(8)

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

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

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

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

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

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

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

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

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

(9)

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

(10)

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

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

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

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

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

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

圖5 模型示意圖Fig.5 A fractured model diagram(a)常規(guī)網(wǎng)格;(b)一級(jí)變網(wǎng)格;(c)二級(jí)變網(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 多級(jí)可變網(wǎng)格積分法與常規(guī)網(wǎng)格積分法對(duì)比

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

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

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

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

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

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

3 微裂縫地震響應(yīng)定量分析

3.1 微裂縫響應(yīng)特征

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

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

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

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

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

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

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

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

圖11 傾斜模型及對(duì)應(yīng)炮記錄Fig.11 Tilt model and corresponding seismic record

20304050607080

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

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

3.4 實(shí)際模型測試

經(jīng)過理論模型的驗(yàn)證，該算法可以快速模型微裂縫的地震特征。圖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)地震接收到的信號(hào)。

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

4 結(jié)論

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

圖12 抽取中心道的波形圖及對(duì)應(yīng)的振幅譜Fig.12 Wave figure from cmp and the corresponding amplitude spectrum

圖13 傾斜模型及對(duì)應(yīng)炮記錄Fig.13 Tilt model and corresponding seismic record

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

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

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

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

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