魏欣偉,薛 姣,羅 霞
(1.中國石油化工股份有限公司勝利油田分公司物探研究院,山東東營257000;2.中國地質大學(武漢)地球物理與空間信息學院,湖北武漢430074;3.地球內部多尺度成像湖北省重點實驗室,湖北武漢430074)
裂縫是油氣儲層的儲集空間和滲濾通道,非常規(guī)致密儲層存在孔隙度低、非均質性強等特點,尋找有效的裂縫發(fā)育帶對致密儲層的勘探和開發(fā)尤為重要。由于地層上覆載荷的壓實作用,儲層中廣泛發(fā)育中、高角度裂縫,造成了巖石的方位各向異性。通常把垂直定向排列裂縫儲層等效為HTI(具有水平對稱軸的橫向各向同性)介質,利用地震波在HTI介質中傳播的方位各向異性信息可以進行裂縫預測。
地震縱波在HTI介質中傳播的方向特性表現(xiàn)為:對于固定偏移距,地震屬性F與測線方位角φ的關系為F(φ)=A+Bcos[2(φ-φ0)],其中系數(shù)B反映了裂縫發(fā)育強度,角度φ0包含了裂縫走向信息[1]。RüGER等[2]在HTI介質線性近似反射系數(shù)研究的基礎上,提出利用AVO(反射振幅隨偏移距變化)梯度、動校正(NMO)速度和橫波分裂分析計算裂縫走向和各向異性參數(shù)。對裂縫方位各向異性比較敏感的地震屬性有NMO速度[2]、反射振幅[3]、縱波旅行時[4]、頻率和衰減[5]等。國內的多個碳酸鹽巖、火山巖裂縫儲層綜合分析實例證明了利用P波地震屬性隨方位角變化的橢圓擬合方法進行裂縫檢測的有效性[6-8]。宋維琪等[9]利用傅里葉級數(shù)形式的P波方位各向異性反射系數(shù)擬合進行裂縫各向異性估計,進一步提高了裂縫檢測的準確性。
隨著裂縫等效介質模型各向異性理論研究的發(fā)展,以及寬方位地震采集和處理技術研究的深入,基于AVOA數(shù)據(jù)反演裂縫參數(shù)成為裂縫儲層預測研究的重點之一。AVOA裂縫參數(shù)反演中常用的裂縫等效介質模型是Schoenberg線性滑動模型[10],利用法向和切向裂縫弱度來表征裂縫,具有類似于HTI介質的彈性剛度矩陣形式。SHAW等[11]基于線性滑動裂縫等效介質模型推導了反射系數(shù)與裂縫弱度之間的關系,提出利用AVOA數(shù)據(jù)反演法向和切向裂縫弱度的方法?;诹芽p介質巖石物理模型,利用疊前地震數(shù)據(jù)AVOA進行裂縫參數(shù)反演能夠充分利用疊前方位地震信息,獲得裂縫等效介質模型參數(shù),實現(xiàn)裂縫儲層定量描述[12-14]。XUE等[15]提出了一種新流體因子表達式,以及方位各向異性反射振幅與法向和切向裂縫弱度的矩陣表達式,利用疊前AVOA技術反演了裂縫密度和新的流體因子。PAN等[16]和CHEN等[17]利用貝葉斯理論McMc(馬爾科夫鏈蒙特卡洛)算法反演法向和切向裂縫弱度。
疊前AVOA數(shù)據(jù)反演在實際地震資料中的應用得益于寬方位地震數(shù)據(jù)采集和處理的發(fā)展。OVT技術[18-20]是一種有效的寬方位三維地震數(shù)據(jù)處理技術,OVT域道集具有豐富的炮檢距和方位角信息,基于OVT域地震數(shù)據(jù)可以有效地實現(xiàn)方位各向異性分析與反演研究,提高裂縫儲層預測的精度與可靠性[18]。目前,OVT域地震數(shù)據(jù)主要用于分方位地震屬性橢圓擬合[7,19]或Thomsen各向異性參數(shù)反演[20]。
常規(guī)AVOA裂縫參數(shù)反演方法先反演法向和切向裂縫弱度,然后根據(jù)線性滑動模型與Hudson扁圓幣狀模型[21]的等價關系計算裂縫密度。本文基于含油和含氣條件下的裂縫等效介質模型開展AVOA裂縫密度直接反演方法研究,建立基于OVT域疊前地震數(shù)據(jù)的AVOA裂縫密度反演技術流程,最后基于SK地區(qū)陸上勘探OVT域地震數(shù)據(jù)開展了AVOA裂縫密度反演和裂縫儲層預測,以驗證方法的有效性。
裂縫等效介質模型建立了裂縫參數(shù)與彈性參數(shù)之間的關系。最簡單、最常用的裂縫等效介質模型是Schoenberg線性滑動模型,忽略裂縫的形狀和結構,把裂縫當作無限薄的柔軟平面,或者以線性滑動邊界條件表示的弱度平面。線性滑動模型的基礎是Backus平均[22-23],即裂縫介質的等效柔度矩陣等于各向同性背景介質的柔度矩陣與裂縫柔度矩陣之和。在線性滑動模型中裂縫被當作一個應力連續(xù)的位移間斷面,位移間斷和連續(xù)的應力之間的關系是線性的,因此稱為線性滑動模型。線性滑動模型的剛度矩陣可以用巖石骨架的拉梅常數(shù)λ和μ,法向裂縫弱度ΔN和切向裂縫弱度ΔT來表示[11]:
C=Cb-(λ+2μ)·
(1)
式中:Cb是各向同性背景介質彈性剛度矩陣;g是背景介質橫波與縱波速度比的平方,g=μ(λ+2μ);r=λ/(λ+2μ)。常規(guī)的HTI介質具有3個各向異性參數(shù)ε,γ和δ,其中ε和γ分別用于表征縱波和橫波的方位各向異性。線性滑動模型是一種特殊的HTI介質模型,其中,切向裂縫弱度代表橫波各向異性,與裂縫密度呈正比,法向裂縫弱度代表縱波各向異性,與裂縫密度呈正比,同時也與流體充填情況有關,當裂縫中飽和流體且流體無法流動時,法向裂縫弱度為零。
為了更加直觀地表征裂縫發(fā)育程度,引入Hudson扁圓幣狀裂縫介質模型[21]。Hudson模型假設各向同性彈性介質中嵌入稀疏定向排列的扁圓幣狀(或稱扁橢球狀)裂縫,模擬的是超聲波實驗室條件下飽和巖石的性質[24],裂縫密度e可以表示為:
(2)
式中:N/V表示單位體積內的裂縫個數(shù);a是裂縫半徑(即橢球長軸半徑);φ是扁圓幣狀裂縫孔隙度;α是裂縫扁度(即橢球短軸與長軸之間的比值)。Hudson模型的一階等效彈性剛度矩陣為[23]:
(3)
式中:U11和U33是與裂縫性質有關的無量綱常量。線性滑動模型與Hudson扁圓幣狀裂縫模型對裂縫的描述有很大差別,但是又存在等價性。線性滑動模型的彈性剛度矩陣與扁圓幣狀裂縫介質模型的一階近似具有相同的形式[23]。
(4)
對于流體充填裂縫,流體剪切模量μf=0,切向裂縫弱度簡化為[23]:
(5)
特別地,當裂縫中飽含油、水時,流體體積模量Kf≠0,裂縫的扁度α→0,法向裂縫弱度為[23]:
ΔN=0
(6)
對于含氣裂縫介質,流體體積模量Kf=0,法向裂縫弱度簡化為[23]:
(7)
垂直定向排列裂縫造成了巖石性質的各向異性,利用常規(guī)橢圓擬合方法得到的橢圓扁率是一種定性的裂縫發(fā)育強度表示方法,基于裂縫巖石物理模型進行裂縫密度和裂縫弱度反演則是一種定量裂縫參數(shù)反演方法。切向裂縫弱度與裂縫密度呈正相關關系,法向裂縫弱度與裂縫密度呈正相關,與充填流體的體積模量呈負相關關系。
裂縫等效介質分界面PP波反射系數(shù)可以分為與反射界面兩側各向同性背景介質物性差異有關的反射系數(shù)各向同性部分Riso(θ)以及與定向排列裂縫有關的反射系數(shù)各向異性部分Rani(θ,φ)[14]:
RPP(θ,φ)=Riso(θ)+Rani(θ,φ)
(8a)
(8b)
Rani(θ,φ)=-[g(1-2g)cos2(φ-φ0)sin2θ(1+tan2θ)+g2cos4(φ-φ0)sin2θtan2θ]δΔN+gcos2(φ-φ0)sin2θ[1-sin2(φ-φ0)tan2θ]δΔT
(8c)
式中:θ和φ分別表示入射角和測線方位角;φ0表示裂縫對稱軸方向;vP是縱波速度;Z和G分別表示縱波阻抗和橫波模量;符號“-”和δ分別表示反射界面兩側參數(shù)的平均值和差值,例如δΔN=ΔN2-ΔN1。當入射角θ<30°時,忽略高階項sin2θtan2θ,則有:
(9a)
Rani(θ,φ)=-g(1-2g)cos2(φ-φ0)·sin2θδΔN+gcos2(φ-φ0)sin2θδΔT
(9b)
對于含油裂縫介質,法向裂縫弱度ΔN=0,反射系數(shù)公式簡化為:
(10)
對于含氣裂縫介質,根據(jù)法向裂縫弱度與裂縫密度之間的關系(公式(7)),反射系數(shù)公式簡化為:
(11)
假設已知裂縫對稱軸方向φ0,兩個不同測線方位φ1和φ2的反射系數(shù)之差為:
ΔR(θ)=R(θ,φ2)-R(θ,φ1)=Ce(θ)Δe
(12)
對于含油裂縫,模型參數(shù)正演系數(shù)為:
(13)
對于含氣裂縫,模型參數(shù)正演系數(shù)為:
(14)
對于K+1層裂縫介質,引入地震子波矩陣,將N個入射角不同方位的地震道集差值ΔS=S(φ2)-S(φ1)與裂縫參數(shù)之間的關系表示成矩陣的形式:
(15)
式中:S(θi)=[S(θi,t1)S(θi,t2)S(θi,tK)]T表示第i個入射角對應的地震數(shù)據(jù);W(θi)是子波褶積矩陣;Ce(θi)是由裂縫參數(shù)正演系數(shù)組成的對角矩陣:
Ce(θi)=
(16)
假設已知裂縫走向φ0,對于含油裂縫介質,利用公式(13)計算Ce(θi)矩陣對角線上的正演系數(shù);對于含氣裂縫,利用公式(14)計算Ce(θi)矩陣對角線上的正演系數(shù)。正演方程(公式(15))是一個線性方程,可以簡寫為:
d=Gm
(17)
式中:G是由子波矩陣,模型參數(shù)正演系數(shù)矩陣和一階差分算子組成的線性正演算子;m=[e(t0),e(t1),,e(tK)]T表示裂縫密度向量。
假設地震數(shù)據(jù)噪聲相互獨立且滿足高斯分布,具有相同的方差,數(shù)據(jù)的似然函數(shù)也滿足高斯分布,將數(shù)據(jù)似然函數(shù)與模型先驗分布結合起來得到后驗概率分布,將后驗概率最大化,得到目標函數(shù):
J(m)=(d-Gm)T(d-Gm)+μR(m)
(18)
式中:R(m)是由模型參數(shù)先驗分布得到的正則化項,通過改變系數(shù)μ改變模型約束條件在反演求解過程中占的比重。目標函數(shù)對模型參數(shù)m的導數(shù)為0即可求得目標函數(shù)的最小值,得到反演問題的解。本文采用Tikhonov正則化,并加入非負約束進行迭代反演。
常規(guī)AVAZ(振幅隨方位角變化)分析或裂縫介質各向異性橢圓擬合方法中,通常通過定義疊加量版,對OVT道集數(shù)據(jù)部分偏移距、部分方位角進行疊加,從而進一步提高數(shù)據(jù)信噪比。而在AVOA反演中,需要利用偏移速度從OVT道集中抽取分方位-入射角道集,即從偏移距轉換到(入射角)角度域,加強振幅能量和穩(wěn)定性。
AVOA反演利用的是相同反射點的振幅隨方位角和入射角的變化,而疊前道集中可能存在方位各向異性時差,以及遠偏移距數(shù)據(jù)仍然存在動校正未校平的殘余時差。方位各向異性和動校正誤差引起的方位時差會影響反演結果的準確性,因此,需要進一步對AVOA道集做時差校正。利用滑動時窗波形匹配剩余時差校正技術,將疊后數(shù)據(jù)作為參考道,選定時窗內某一道集的所有道與參考道進行互相關,求取每道最大相關系數(shù)對應的時移量用于時差校正。
AVOA反演中需要已知裂縫走向或裂縫對稱軸方向,首先利用井信息或者地震數(shù)據(jù)的橢圓擬合或AVAZ反演獲取裂縫走向,繼而進行AVOA裂縫參數(shù)反演[15-16,25-26]。
基于OVT道集數(shù)據(jù)的AVOA裂縫參數(shù)反演的流程如圖1所示。
圖1 基于OVT數(shù)據(jù)的AVOA反演流程
建立裂縫等效介質模型,其中各向同性背景介質縱橫波速度、密度和裂縫密度參數(shù)分別為R1,R2和R3,如圖2所示。該模型中存在3段裂縫儲層,裂縫走向為北東向45°,裂縫密度均為0.1。假設裂縫扁度α=10-6,法向裂縫弱度隨充填流體體積模量的增大而減小(圖3),當裂縫中充填油或水時可忽略法向弱度,切向裂縫弱度不隨充填流體變化。
圖2 裂縫介質模型縱橫波速度(a)、密度(b)和裂縫密度(c)
圖3 裂縫弱度隨充填物體積模量的變化(紅圈表示含氣、含油和含水裂縫介質的法向裂縫弱度)
利用本文提出的AVOA裂縫密度反演算法對裂縫儲層模型的合成地震記錄進行試算。假設裂縫儲層中飽和充填油,利用Zoeppritz方程計算反射系數(shù),并將反射系數(shù)與30Hz雷克子波進行褶積得到合成疊前地震記錄,然后分別添加20%,50%和100%的隨機噪聲,得到的合成地震記錄信噪比分別為RS/N=5,RS/N=2和RS/N=1。首先,對合成疊前地震數(shù)據(jù)進行AVAZ裂縫走向反演,圖4顯示,當信噪比大于等于2時,儲層段裂縫走向反演結果與裂縫模型走向相吻合,非裂縫儲層段(裂縫密度為0的位置)裂縫走向反演受噪聲影響較大。利用與裂縫對稱軸夾角分別為0和90°的疊前地震數(shù)據(jù)(圖5和圖6)進行AVOA裂縫密度反演,其結果如圖7所示。由于裂縫密度不可能是負值,因此在反演過程中加入非負約束條件。圖7顯示,當信噪比大于等于2時滿足反演要求,當信噪比為1時反演結果與裂縫模型真實值有較大差異。
圖4 不同信噪比疊前地震數(shù)據(jù)裂縫走向AVAZ反演a 不含噪聲; b 含20%噪聲; c 含50%噪聲; d 含100%噪聲
圖5 不同信噪比疊前地震數(shù)據(jù)(Ⅰ)a 不含噪聲; b 含20%噪聲
圖6 不同信噪比的疊前地震數(shù)據(jù)(Ⅱ)a 含50%噪聲; b 含100%噪聲
圖7 不同信噪比疊前地震數(shù)據(jù)裂縫密度反演結果(紅線)與模型真實值(藍線)對比a 不含噪聲; b 含20%噪聲; c 含50%噪聲; d 含100%噪聲
假設裂縫模型為干裂縫或者含氣裂縫,依然用含油情況下的AVOA裂縫密度反演方法進行反演,則其裂縫密度反演結果小于真實值(圖8)。
圖8 干裂縫介質模型裂縫密度反演結果(紅線)與模型真實值(藍線)對比a 不含噪聲; b 含20%噪聲; c 含50%噪聲; d 含100%噪聲
中國中東部某探區(qū)致密砂巖油氣資源潛力大,巖性致密、基質孔隙結構復雜,具有較強的非均質性,裂縫的發(fā)育程度對致密砂巖儲層的控制作用非常明顯,該探區(qū)中SK工區(qū)多層系含油,裂縫較為發(fā)育,某些井試油日產達百噸以上。實際地震數(shù)據(jù)為SK工區(qū)陸上勘探得到的OVT道集。
SK工區(qū)OVT道集覆蓋次數(shù)為93~312次,不同位置覆蓋次數(shù)分布不均勻。方位角-偏移距變化關系分析圖(圖9)顯示OVT道集偏移距范圍為300~6000m,南北方向炮檢距分布均勻,覆蓋次數(shù)較高,然而橫縱比較小,方位分布不均勻;南北方向最大偏移距達6000m,而東西方向最大偏移距不足2500m,東西方向遠偏移距信息缺失。偏移距大于2500m的部分,方位角主要分布在0~60°和120°~180°之間,如果使用偏移距大于2500m的疊前地震數(shù)據(jù)進行AVOA反演則會產生較大誤差。經過分析,有效的偏移距范圍是300~2500m,其覆蓋次數(shù)和方位角分布較為均衡。
圖9 方位角-偏移距分析
對OVT數(shù)據(jù)進行分方位-入射角道集抽取時,需要考慮各面元覆蓋次數(shù)以及方位角和入射角個數(shù)之間的平衡。理論上講,方位角和入射角個數(shù)越多,反演的精度越高,但隨著方位角和入射角的增多,覆蓋次數(shù)變少,AVOA道集數(shù)據(jù)的信噪比逐漸降低。
雖然目的層T6和T7界面最大入射角約為50°,但由于不同方位炮檢距分布不均,在抽取分方位-入射角道集時應充分考慮最大有效炮檢距(2500m)對應的入射角作為最大入射角(約30°)。經過試驗最終確定劃分為4個方位和3個入射角,生成中心方位角分別為22.5°,67.5°,112.5°和157.5°,中心入射角分別為10°,17°和24°的分方位-入射角道集(圖10a)。方位角劃分范圍分別為0~45°,45°~90°,90°~135°和135°~180°,入射角劃分范圍分別為6°~13°,13°~20°,和20°~27°。選用合適的滑動時窗參數(shù)進行剩余時差校正后(圖10b),方位各向異性和動校正誤差引起的時差被校正,同相軸被拉平。
圖10 時差校正前(a)、后(b)分方位-入射角道集
在疊前OVT道集資料分析的基礎上,抽取分方
位-入射角道集,利用滑動時窗波形匹配時差校正技術消除方位各向異性和動校正剩余時差,然后求取相互垂直的角道集之差:
(19)
式中:S(θi,φj)表示第i個入射角、第j個方位角對應的地震道;對于SK實際地震數(shù)據(jù),S(θi,φ3)-S(θi,φ1)表示入射角為θi時第3個方位(φ3=112.5°)與第1個方位(φ1=22.5°)地震道之間的差異。首先,利用AVAZ反演估計裂縫走向,圖11顯示研究區(qū)目的層(T6上10ms到T6下10ms)裂縫走向主要為近東西向,裂縫走向與斷裂帶的延展方向一致。然后,應用本文提出的AVOA裂縫參數(shù)反演方法估計裂縫密度,達到裂縫儲層預測的目的。測井、鉆井、巖心和試油資料表明SK工區(qū)多層系含油,為裂縫型含油致密砂巖儲層,因此將含油裂縫的模型參數(shù)正演系數(shù)(公式(13))代入模型正演方程(公式(15))中,然后利用目標函數(shù)(公式(18))最小化進行AVOA反演計算。裂縫密度沿T6層切片(圖12)顯示在研究區(qū)西北部yi177井到y(tǒng)i178井之間裂縫密度較大,為裂縫發(fā)育帶。裂縫的分布主要受斷裂帶控制,裂縫的分布與斷裂系統(tǒng)的分布具有一定相關性,裂縫密度較大的位置位于斷裂帶周圍和斷裂帶交會處。
圖11 沿T6層裂縫走向切片(藍線走向表示裂縫走向,藍線長度表示各向異性程度)
圖12 沿T6層裂縫密度切片
將測井資料解釋成果與AVOA裂縫密度反演結果進行對比分析,驗證基于AVOA裂縫密度反演的裂縫儲層預測效果。yi176井T6~T7含油層段為粉細砂巖,綜合解釋為油層,常規(guī)測井曲線上具有“高聲波時差、高中子、低電阻、低密度和自然電位負異?!钡奶卣?。圖13是過yi176井(線號702)裂縫密度反演剖面,其中,粉紫色曲線顯示了正交多極子聲波陣列測井得到的各向異性參數(shù),黑色曲線是井旁地震道的裂縫密度反演結果。井旁地震道的裂縫密度反演結果(黑色曲線)與聲波陣列各向異性(粉紫色曲線)異常相吻合,裂縫密度反演結果顯示在T6到T7層中間及T7層上部裂縫較為發(fā)育,測井解釋為裂縫儲層。裂縫密度反演結果顯示裂縫儲層在橫向上呈似層狀和團塊狀分布。
圖13 過yi176井裂縫密度剖面
本文基于裂縫等效介質模型,分別推導了含油和含氣情況下反射系數(shù)與裂縫密度之間的關系,在此基礎上提出了基于疊前AVOA反演的裂縫密度估計方法和基于OVT域疊前數(shù)據(jù)的裂縫密度反演技術流程。SK工區(qū)裂縫參數(shù)反演結果與測井資料解釋成果相吻合,顯示了OVT域AVOA裂縫參數(shù)反演的有效性及其在SK工區(qū)裂縫儲層預測中的適用性?;贠VT域地震數(shù)據(jù)的AVOA裂縫參數(shù)反演結果能夠定量刻畫裂縫分布特征,為SK工區(qū)致密砂巖含油裂縫儲層預測和開發(fā)提供了重要依據(jù)。