印興耀,劉志國,李春鵬,,李愛山,蓋海洋

(1.中國石油大學(華東)地球科學與技術學院,山東青島266580;2.中海油研究總院,北京100028)

裂縫作為地殼中一種普遍的構造現(xiàn)象,廣泛存在于各類巖石中。到目前為止,已在砂巖、泥頁巖和碳酸鹽巖,甚至火成巖等各類巖石的裂縫型儲層中獲得了大量的工業(yè)油氣流。地震波在裂縫型地層中傳播會發(fā)生橫波分裂現(xiàn)象,Ata等[1]指出快橫波極化方向平行于裂縫走向,慢橫波極化方向垂直于裂縫走向,為利用橫波技術進行裂縫勘探奠定了理論基礎。但是實際橫波地震勘探成本非常昂貴,因此需要發(fā)展成本相對比較低廉的裂縫型地層縱波勘探技術。裂縫型地層具有較強的縱波方位各向異性,梁鍇[2]和思薌等[3]分別推導了二維和三維情況下TI介質反射透射方程;Grechka等[4]研究了方位AVO梯度與裂縫參數之間的關系;劉洋等[5]、范國章等[6]討論了裂縫傾角與AVO隨方位變化的關系,傾角越小,AVO方位變化特征也越小;Liu等[7]研究了薄裂縫層的方位AVO特征,當薄層厚度大于1/4地震波長時可以得到準確的裂縫密度和走向;Gray等[8]認為方位AVO的各向異性程度與裂縫開啟、流體填充和裂縫傾角等因素有關。因此,通過方位縱波地震數據提取地層方位AVO梯度識別地層裂縫的方法是切實可行的。

常規(guī)方位AVO梯度提取方法對各個方位地震道集獨立進行,該方法提取的AVO截距項即自激自收反射系數在各個方位可能會不一致,進而會影響方位AVO梯度的精度。雖然方位AVO梯度反演方法可以將截距項和各方位梯度項同時引入反演方程,得到截距項一致的方位AVO梯度,但該方法需要各向異性先驗信息以穩(wěn)定反演結果。為了克服常規(guī)方位AVO梯度提取方法截距項難以一致和方位AVO梯度反演方法所需各向異性測井信息匱乏的問題,我們發(fā)展了穩(wěn)定的方位AVO梯度無約束反演方法。據此研究了截距項一致情況下4類含氣砂巖方位AVO梯度與裂縫參數之間的關系,并且通過模型試算和實際工區(qū)應用驗證了該反演方法的有效性,為裂縫型儲層預測提供了有利的技術支撐。

1 方位AVO梯度與裂縫參數

近垂直定向分布的裂縫型地層具有方位各向異性特征,它導致AVO梯度隨著方位變化而變化。因此,有必要研究方位AVO截距項一致情況下,方位AVO梯度與地層裂縫參數之間的關系,為利用方位AVO梯度預測裂縫密度和裂縫走向提供理論基礎。

1.1 方位AVO梯度與裂縫密度

Rutherford等[9]和Castagna等[10]描述了4類含氣砂巖的AVO特征,這4類含氣砂巖模型的縱波速度、橫波速度和密度如表1所示。假設這4類含氣砂巖是由石英、孔隙和氣體構成的,根據《巖石物理手冊》[11],相關彈性參數見表2。

表1 4類含氣砂巖彈性參數

表2 常見礦物彈性參數

在以上4類含氣砂巖中嵌入近垂直定向分布的裂縫,裂縫密度e分別為0.05,0.10,0.15,裂縫縱橫比為0.01,裂縫走向為90°。根據李春鵬等[12]推導的HTI介質彈性波反射透射方程計算4類含氣砂巖模型的方位反射系數,并且將反演的方位AVO梯度擬合成橢圓,以研究方位AVO梯度與裂縫密度之間的關系。

圖1至圖5顯示了4類含氣砂巖裂縫密度e分別為0.05,0.10,0.15時的方位AVO梯度和橢圓擬合結果,圖中方位AVO梯度近似為一條周期為π的余弦曲線,且裂縫走向與裂縫傾向的方位AVO梯度差異最大;裂縫密度越大,第Ⅰ,Ⅱ-1,Ⅱ-2和Ⅲ類含氣砂巖裂縫走向與裂縫傾向的方位AVO梯度差異也越大,但是第Ⅳ類含氣砂巖的差異卻先增大后減小。

方位AVO梯度可以擬合成橢圓,表3顯示了4類含氣砂巖裂縫密度e分別為0.05,0.10,0.15時的橢圓率,橢圓率=(橢圓長軸-橢圓短軸)/橢圓長軸。表3中裂縫密度越大,第Ⅰ,Ⅱ-1,Ⅱ-2和Ⅲ類含氣砂巖橢圓率也越大,但是第Ⅳ類含氣砂巖的橢圓率先增大后減小。大體來說,裂縫密度越大相應的方位AVO梯度橢圓率也越大,據此可以用于實際工區(qū)裂縫密度預測研究。

圖1 第Ⅰ類含氣砂巖裂縫密度e分別為0.05,0.10,0.15時的方位AVO梯度(a)和方位AVO梯度的橢圓擬合結果(b)

圖2 第Ⅱ -1類含氣砂巖裂縫密度e分別為0.05,0.10,0.15時的方位AVO梯度(a)和方位AVO梯度的橢圓擬合結果(b)

圖3 第Ⅱ -2類含氣砂巖裂縫密度e分別為0.05,0.10,0.15時的方位AVO梯度(a)和方位AVO梯度的橢圓擬合結果(b)

砂巖模型方位AVO梯度橢圓率e=0.05e=0.10e=0.15Ⅰ0.02400.07500.1193Ⅱ-10.03380.08520.1297Ⅱ-20.04050.09430.1404Ⅲ0.05080.10990.1594Ⅳ0.39460.63810.1211

圖4 第Ⅲ類含氣砂巖裂縫密度e分別為0.05,0.10,0.15時的方位AVO梯度(a)和方位AVO梯度的橢圓擬合結果(b)

圖5 第Ⅳ類含氣砂巖裂縫密度e分別為0.05,0.10,0.15時的方位AVO梯度(a)和方位AVO梯度的橢圓擬合結果(b)

1.2 方位AVO梯度與裂縫走向

上述第Ⅳ類含氣砂巖的第1層為均勻各向同性介質,第2層為裂縫介質,裂縫密度是0.05,裂縫縱橫比為0.005,裂縫走向為90°,假設第2層裂縫介質之下有一個均勻各向同性介質且彈性參數與第1層相同。通過模型試算研究方位AVO梯度橢圓半軸方向與裂縫走向之間的關系。圖6a顯示了第Ⅳ類含氣砂巖頂、底界面的方位AVO梯度,圖中頂、底界面的方位AVO梯度呈反相關關系,頂界面裂縫走向處的AVO梯度最大,底界面裂縫走向處的AVO梯度最小。圖6b顯示了第Ⅳ類含氣砂巖頂、底界面方位AVO梯度橢圓擬合結果,圖中頂界面橢圓長軸指示裂縫走向,底界面橢圓短軸指示裂縫走向。因此,方位AVO梯度橢圓半軸可以指示裂縫走向,但是卻存在90°不確定性。

圖6 第Ⅳ類含氣砂巖頂、底界面方位AVO梯度(a)和方位AVO梯度橢圓擬合結果(b)

2 穩(wěn)定方位AVO梯度反演方法

常規(guī)方位AVO梯度提取方法是在各個方位獨立進行的,該方法提取的AVO截距項即自激自收反射系數在各個方位可能會不一致,進一步會影響方位AVO梯度的精度。為了克服各個方位AVO截距不一致的問題和實際工區(qū)各向異性測井信息匱乏的問題,我們發(fā)展了穩(wěn)定的方位AVO梯度無約束反演方法。

2.1 假設方位地震數據是方位反射系數

入射角小于30°時,方位反射系數公式為

(1)

式中:R表示方位反射系數;A表示AVO截距;B表示方位AVO梯度;φ表示方位角;θ表示入射角。地震數據S的正演方程為

(2)

式中:W表示子波矩陣;R表示方位反射系數序列;G表示反射系數正演算子矩陣;m=[AB1B2…BN]T表示待反演參數,其中,B1,B2,…,BN表示各個方位的AVO梯度序列,N表示方位角數目。

假設方位地震數據是方位反射系數,則方位地震數據的正演方程為

(3)

(2)式的最小二乘反演方程為

(4)

(3)式的最小二乘反演方程為

(5)

(5)式相對于(4)式消除了地震子波的影響,增強了反演方程的穩(wěn)定性。

2.2 最少反射層假設降低反褶積維度

通過地震反褶積將上述得到的方位AVO梯度地震數據體wB1,wB2,…,wBN轉換為方位AVO梯度數據體B1,B2,…,BN,然后通過橢圓擬合得到方位AVO梯度橢圓率和橢圓半軸方向。但是,常規(guī)反褶積方法得到的反射系數不是稀疏的,不能真正反映地下地層界面的真實情況。我們采取Perez等[13]提出的最少反射層地震反演方法,搜索準確的地下界面位置,以降低方位AVO梯度地震記錄反褶積方程維數和條件數,得到穩(wěn)定的反演結果。該方法假設地下地層是稀疏分布的,通過降低小子波矩陣w維數和反褶積系數矩陣wTw條件數,以獲得穩(wěn)定的反演結果。最少反射層稀疏脈沖反演的關鍵性技術是模擬退火法尋找反射層位置和線性優(yōu)化方法解目標函數[14-16]。假設單道采樣點數為L,采樣率為dt。模擬退火法實現(xiàn)步驟如下。

1) 給定反射層個數M,反射層位置τj范圍[0,L·dt],子波主頻f范圍[fmin,fmax]和初始相位φ范圍[φmin,φmax]；給定初始溫度T0,回火升溫次數Nb,每次降溫迭代次數Nt,收斂閾值ε,子波尺度因子s,子波時間長度tw。轉到步驟2)。

2) 令i=1,b=0,T=T0,隨機給定τj,f和φ初始值。

(6)

式中:i0是地震道下標矢量;floor(·)是四舍五入函數;frand(·)是均勻分布隨機函數;fsort(·)是從小到大排序函數。轉到步驟3)。

4) 接受方式。

5) 降溫方式。

T=T0e-C(i-1),其中,C是降溫因子,C越小降溫越慢。轉到步驟6)。

6) 模型位置擾動。

令i=i+1。

如果i>Nt,則i=1,b=b+1；

此法是在顯微鏡下直接進行測定，方便快捷并且儀器損耗較小，但在一定的容積中微生物的個體數目包括死活細胞均被計算在內，還有微小雜物也被計算在內，這樣得出結果往往偏高，因此適用于對形態(tài)個體較大的菌體計數。

如果b>Nb,轉到步驟8)；

如果b≡0,則

(7)

否則

(8)

式中:fsign(·)是符號函數;funidrnd(M,1)表示從1到M隨機找一個整數。轉到步驟7)。

7) 子波擾動。

τj(τj<0‖τj>L·dt)=frand·L·dt

τj=fsort(τj)

f(ffmax)=fmin+frand·

(fmax-fmin)

fsign(frand-0.5)·(φmax-φmin)

φ(φ<φmin‖φ>φmax)=φmin+frand·

(φmax-φmin)

(9)

其中,λ是擾動因子,λ越大擾動時間越長,轉到步驟3)。

3 模型試算

為了驗證穩(wěn)定方位AVO梯度反演方法的有效性,設計4層三維模型如圖7所示,第1層、第2層和第4層是各向同性介質,其中第1層介質的縱波速度、橫波速度、密度和厚度分別為2500m/s,1500m/s,2000kg/m3和20m,第2層介質的縱波速度、橫波速度、密度和厚度分別是3000m/s,1800m/s,2300kg/m3和80m,第4層介質的縱波速度、橫波速度和密度分別是3200m/s,1850m/s,2350kg/m3。第3層介質是80m厚的近垂直定向分布的裂縫介質,該裂縫介質的背景彈性參數和物性參數取自于某油田A井的測井信息,其中A井S地層的縱波速度、橫波速度、密度、石英含量、泥質含量和干酪根含量的均值分別是3476m/s,1935m/s,2400kg/m3,11.28760,17.58680和5.24508。假設該模型inline和xline方向各有100條測線,圖8給出了第3層裂縫介質的裂縫密度平面分布、裂縫走向平面分布和裂縫走向玫瑰圖。

根據李春鵬等[12]推導的HTI介質彈性波反射透射方程可以計算該模型的方位縱波反射系數。模型正演采用相移30°后的50Hz雷克子波,入射角分別是15°,25°和35°,方位角分別是0,30°,60°,90°,120°,150°,180°。圖9顯示了inline75,xline50處,入射角為35°,信噪比為10，采樣間隔為1ms的方位縱波地震記錄,圖中第15采樣點左右處的同相軸表示模型第1層和第2層的各向同性界面的方位地震記錄,各向同性界面的地震記錄與方位角無關;圖中第70和第150采樣點左右處的同相軸分別表示第3層裂縫介質頂面和底面的方位地震記錄,裂縫介質頂面和底面的地震記錄隨方位角變化而變化。

圖7 4層介質模型示意圖解

圖8 4層介質模型的第3層裂縫介質參數a 裂縫密度平面分布; b 裂縫走向平面分布(黑色小線段長度指示裂縫密度大小,小線段方向指示裂縫走向); c 裂縫走向玫瑰圖

圖9 4層介質模型inline75,xline50處,入射角為35°,信噪比為10的方位地震記錄

穩(wěn)定方位AVO梯度反演可以分為兩步。

第一步:假設方位地震數據是方位反射系數,根據(5)式的反演方程,4層介質模型的反演矩陣GTG條件數約為527,但若采用(4)式的反演方程,4層介質模型的反演矩陣(WG)T(WG)條件數約為7.3×1017,可見(5)式的反演方程穩(wěn)定性相對較高。4層介質模型的inline75,xline50處反演結果見圖10a和圖10b,圖中反演的AVO截距地震記錄和實際自激自收反射地震記錄吻合得較好,并且各向同性界面的AVO梯度地震記錄基本與方位無關,但第3層裂縫地層頂、底的AVO梯度地震記錄隨方位變化,說明反演結果能夠指示地層裂縫發(fā)育情況。

第二步:方位AVO梯度地震數據體wB1,wB2,…,wBN轉換為方位AVO梯度數據體B1,B2,…,BN,但是反褶積算子wTw的條件數約為1.6×1018,一般采用正則化算法改善反褶積算子,但是這種常規(guī)反褶積方法得到的反射系數不是稀疏的,不能真正地反映下地層界面的真實情況。我們采用最少反射層假設反演C和D數據體,假設4層介質模型有6個地層分界面,inline75,xline50處反演結果見圖10c,圖中各向同性界面的AVO梯度基本與方位無關,但裂縫地層頂、底的AVO梯度隨方位變化,說明反演結果能夠指示地層裂縫發(fā)育情況。其余地層分界面處幾乎為0,說明最少反射層假設方法不僅能搜索真實地層界面,而且在非裂縫地層分界面處得到的值幾乎為0。

對4層介質模型基于穩(wěn)定方位AVO梯度反演方法得到的結果進行方位AVO梯度橢圓分析(圖11)。比較圖11a和圖8a可以發(fā)現(xiàn),方位AVO梯度橢圓率基本可以反映裂縫地層裂縫密度的平面展布;比較圖11b和圖8b可以發(fā)現(xiàn),方位AVO梯度橢圓半軸方向與裂縫地層裂縫走向基本一致或垂直;圖11c中方位AVO梯度橢圓半軸主方向與裂縫走向主方向一致。由此可見,穩(wěn)定方位AVO梯度反演方法得到的方位AVO梯度橢圓率可以指示裂縫密度,方位AVO梯度橢圓半軸方向可以指示裂縫走向。

圖10 4層介質模型inline75,xline50處的穩(wěn)定方位AVO梯度反演結果a 第一步反演AVO截距地震記錄(黑)與真實自激自收地震記錄(紅)對比; b 第一步反演的方位AVO梯度地震記錄; c 第二步反演的方位AVO梯度

圖11 4層介質模型穩(wěn)定方位AVO梯度反演結果的方位AVO梯度橢圓分析a 橢圓率平面分布; b 橢圓半軸方向平面分布(黑色小線段長度指示裂縫密度大小,小線段方向指示裂縫走向); c 橢圓半軸方向玫瑰圖

4 實際裂縫地層預測應用

某油田發(fā)育大段泥頁巖裂縫地層,工區(qū)內地震、測井資料比較完備,特別是裂縫比較發(fā)育的S地層疊后道集上存在著明顯的方位變化,適合于利用地震方位道集數據進行裂縫型儲層預測研究。對該油田三維寬方位觀測系統(tǒng)的所有可能的炮檢-方位角分析可知,所觀測地震數據隨炮檢距和方位角的分布比較均勻,只是在大炮檢距的小方位角和大方位角數據有點缺失。通過地震資料處理,得到分別為0～30°,30°～60°,60°～90°,90°～120°,120°～150°,150°～180°的6個方位疊加資料。利用穩(wěn)定方位AVO梯度反演方法預測S地層的裂縫密度、走向和發(fā)育帶,圖12給出了預測分析結果。圖12中A井附近存在方位AVO梯度異常,穩(wěn)定方位AVO梯度反演結果的方位AVO梯度方向與測井解釋結果一致,但存在90°不確定性。

圖12 某油田S地層穩(wěn)定方位AVO梯度反演裂縫預測結果與測井解釋結果a 各向異性梯度; b 各向異性梯度方向玫瑰圖; c A井測井解釋結果

5 結束語

針對常規(guī)方位AVO梯度提取方法各個方位AVO截距不一致的問題和實際工區(qū)各向異性測井信息匱乏的問題,我們研究了方位AVO截距項一致情況下方位AVO梯度與裂縫參數之間的關系,提出了穩(wěn)定方位AVO梯度無約束反演方法,得到以下3點認識。

1) 4類含氣砂巖的巖石物理模型正演研究結果表明,方位AVO梯度橢圓率能夠指示裂縫發(fā)育區(qū);橢圓半軸方向可以指示裂縫走向,但是存在90°不確定性。

2) 穩(wěn)定方位AVO梯度無約束反演方法能夠提高方位AVO梯度反演方程穩(wěn)定性,并能得到準確的地下地層方位AVO梯度信息。

3) 實際工區(qū)應用結果表明,利用穩(wěn)定方位AVO梯度無約束反演方法得到的方位AVO梯度可以較為準確地預測地層裂縫發(fā)育帶,為裂縫型儲層預測提供了有利的技術支撐。

參 考 文 獻

[1] Ata E B,Michelena R J.Mapping distribution of fractures in a reservoir with P-S converted waves[J].The Leading Edge,1995,14(6):664-674

[2] 梁鍇.TI介質地震波傳播特征與正演方法研究[D].山東:中國石油大學(華東),2009

Liang K.The study on propagation feature and forward modeling of seismic wave in TI media[D].Shandong:China University of Petroleum,1996

[3] 司薌,吳國忱.裂隙等效TTI介質qP波反射特征研究[J].地球物理學進展,2012,27(5):2091-2099

Si X,Wu G C.Study on qP-wave reflectivity in fracture equivalent TTI media[J].Progress in Geophysics (in Chinese),2012,27(5):2091-2099

[4] Grechka V,Tsvankin I.NMO-velocity surfaces and Dix-type formulas in anisotropic heterogeneous media[J].Geophysics,2002,67(3):939-951

[5] 劉洋,董敏煜.各向異性介質中的方位AVO[J].石油地球物理勘探,1999,34(3):261-269

Liu Y,Dong M Y.Azimuthal AVO in anisotropic medium[J].Oil Geophysical Prospecting,1999,34(3):261-269

[6] 范國章,牟永光,金之鈞.裂縫介質中地震波方位AVO特征分析[J].石油學報,2002,23(4):42-46

Fan G Z,Mou Y G,Jin Z J.Amplitude versus offset characteristic of Azimuth of seismic wave in fractured reservoir[J].Acta Petrolei Sinica,2002,23(4):42-46

[7] Liu Y L,Ziolkowski A,Liu E R,et al.Characterization of azimuthal anisotropy in the presence of thin layers using P-waves[J].Expanded Abstracts of 71stAnnual Internat SEG Mtg,2001,169-172

[8] Gray D,Roberts G,Conoco H K.Recent advances in determination of fracture strike and crack density from P-wave seismic data[J].The Leading Edge,2002,21(3):280-285

[9] Rutherford S R,Williams R H.Amplitude-versus-offset variations in gas sands[J].Geophysics,1989,54(6):680-688

[10] Castagna J P,Swan H W,Foster D J.Framework for AVO gradient and intercept interpretation[J].Geophysics,1998,63(3):948-956

[11] Mavko G,Mukerji T,Dnorikin J.The rock physics handbook:tools for seismic analysis in porous media[M].Cambridge:Cambridge University Press,1998:260-263

[12] 李春鵬,印興耀,張峰.HTI介質飽和流體特性和裂縫密度對方位反射系數的影響[J].石油物探,2013,52(1):1-10

Li C P,Yin X Y,Zhang F.Influence of HTI medium saturated fluid properties and fracture density on azimuthal reflectivity[J].Geophysical Prospecting for Petroleum,2013,52(1):1-10

[13] Perez O D,Velis R D.Sparse-spike AVO/AVA attributes from prestack data[J].Expanded Abstracts of 81stAnnual Internat SEG Mtg,2011,340-344

[14] 趙改善.用模擬退火法作AVA反演[J].石油物探,1992,31(4):1-10

Zhao G S.AVA inversion by use of simulated annealing[J].Geophysical Prospecting for Petroleum,1992,31(4):1-10

[15] 王振國,陳小宏,王克斌,等.多參數優(yōu)化的模擬退火法波阻抗反演[J].石油物探,2007,46(2):120-125

Wang Z G,Chen X H,Wang K B,et al.Simulated annealing of multi-parameter optimization method of wave impedance inversion[J].Geophysical Prospecting for Petroleum,2007,46(2):120-125

[16] 李大衛(wèi),尹成,謝兵.模擬退火獨立分量分析方法及應用[J].石油物探,2007,46(1):24-27

Li D W,Yin C,Xie B.Independent component analysis method and simulated annealing[J].Geophysical Prospecting for Petroleum,2007,46(1):24-27