基于黏彈性Chapman-Kelvin模型的裂縫儲層頻變AVAZ反演研究

2024-02-04 07:07張青廖建平劉和秀周林張學娟

地球物理學報 2024年2期

張青,廖建平,*,劉和秀,周林,張學娟

1 湖南科技大學地球科學與空間信息工程學院,湘潭 411201 2 中國石化地球物理重點實驗室,中石化石油物探技術研究院有限公司,南京 211103 3 重慶科技學院非常規油氣開發研究院,重慶 401331

0 引言

隨著油氣勘探逐漸從構造油氣藏轉向非常規油氣藏,地下裂縫在油氣生產中發揮著重要作用(Zeng et al.,2013;Ghasemi et al.,2020;Lei et al.,2022).Sullivan等(2006)表明全球三分之一以上的油氣儲存在碳酸鹽巖儲層中,并且大部分為裂縫儲層.對于致密儲層,裂縫不僅可以連通孤立孔隙,增加儲層有效孔隙度,提高滲透率,更是儲層中主要的滲流通道(Lei and Wang,2016;Shi et al.,2018;Hou et al.,2021).因此,裂縫描述是包括碳氫化合物和地下水的生產和監測二氧化碳的地質儲存等各種應用的一個關鍵問題(Iding and Ringrose,2010;劉敬壽等,2019;Sheng et al.,2021;Yasin et al.,2022).Liu和Martinez(2012)表明,可以通過露頭、巖心或鉆孔測井觀測和地震資料解釋等方法來描述裂縫的特征.通常需要將所有這些方法與不同尺度的測量數據結合起來,對裂縫儲層進行準確的描述.

裂縫儲層最基本的特征是各向異性(Cao et al.,2018;Song et al.,2019,2020;Jiang et al.,2021).地震各向異性是描述裂縫特征最常用的方法,并且為了建立各向異性參數與裂縫系統特征之間的聯系,使用裂縫儲層的等效介質理論來描述裂縫巖石(Pan et al.,2018;Ding et al.,2020;Wang et al.,2022).對裂縫的假設不同,目前常用的等效介質理論方法有兩種.其一是基于夾雜物的Hudson模型(Hudson,1981),將裂縫視為嵌入各向同性固體中的橢球體孔隙;另一種是線性滑動模型(Schoenberg,1980),認為裂縫具有線性滑動邊界條件的薄弱面.Hudson模型和Schoenberg(1980)提出的線性滑動模型是常用的靜態等效介質模型.當儲層中存在一組垂直排列的裂縫時,該儲層可以看作具有水平對稱軸的橫向各向同性介質,基于這種靜態等效介質模型的裂縫預測方法無法獲取裂縫長度、裂縫密度等裂縫參數信息,也不能預測出儲層的孔隙度、含水飽和度等物性信息.同樣,這類模型也不能解釋在實際裂縫儲層巖石中地震波的速度頻散和衰減現象,這將與真實裂縫儲層巖石中的頻變各向異性理論相違背(李博南等,2017;Kumar et al.,2018;潘新朋,2019).流體在裂縫與孔隙之間流動的動態等效介質理論模型,更加全面地描述出實際裂縫儲層中的頻變各向異性特征(Kobayashi and Mavko,2016;Gao and Deng,2022).這類模型通過地震波致流作用,可以反映出真實裂縫儲層巖石中的地震波頻散和衰減現象.Thomsen(1995)通過考慮裂縫和孔隙之間的流體傳遞,認識到頻散對解釋地震各向異性的重要性;Hudson等(1996)建立了一個動態等效介質模型,該模型考慮了整個頻率范圍內地震波傳播特性,但在低頻極限時與Thomsen(1995)等徑孔隙模型的相應結論不一致,所以不能在全頻段內有效;Chapman(2003)基于Chapman等(2002)的顆粒噴射流的孔彈性理論,結合一組定向排列的介觀尺度硬幣型裂縫,考慮二者之間的相互作用,提出了各向異性模型,通過對比各向異性Gassmann方程和Hudson模型在低頻和高頻極限下的結論,發現與其一致,說明Chapman(2003)動態等效介質模型可以反映出實際裂縫儲層巖石中的地震波的頻散和衰減現象.

近年來,一些學者開展了基于動態等效介質模型的裂縫儲層參數反演的研究(劉倩等,2016;潘新朋等,2018;Li et al.,2017;Zhang et al.,2022).對于裂縫儲層的反演,主要有定性反演和定量反演兩種類型.其中,定量反演是指基于地震各向異性來反演裂縫密度和裂縫走向等裂縫信息,包括疊前AVAZ反演等方法(陳懷震等,2014;劉宇巍等,2016;葛子建等,2018;陳珂磷等,2022).Maultzsch等(2003)根據Chapman(2003)模型,反演了裂縫密度和裂縫半徑,通過反演得到的裂縫參數信息與裂縫鉆孔成像的觀測結果一致,說明Chapman(2003)模型的有效性;Ali和Jakobsen(2011)通過T矩陣方法求取了裂縫孔隙介質的等效剛度矩陣系數,并且通過求解Christoffel方程,建立裂縫密度和裂縫走向與地震波速度頻散和衰減之間的關系,同時反演了裂縫密度和裂縫走向信息;基于Chapman(2003)模型,Al-Harrasi等(2011)使用網格搜索法對裂縫走向、裂縫密度和裂縫長度進行反演,反演得到的裂縫尺寸與地層露頭觀測結果一致;Ali和Jakobsen(2014)結合T矩陣法和HTI介質PP波反射系數的Rüger近似方程(Rüger,1998),計算了含有一組垂直排列的介觀尺度裂縫的裂縫孔隙介質反射PP波頻變AVAZ數據,通過T矩陣法所計算的等效剛度系數建立了裂縫參數與頻變AVAZ數據的關系,結合貝葉斯方法和蒙特卡羅隨機采樣方法對裂縫密度、裂縫走向、裂縫開度和裂縫半徑進行了同時反演,并表明只要在反演過程中采用動態等效介質模型,并提供有效的裂縫長度先驗信息,可以改善裂縫儲層的特征;Liu等(2021)將頻變各向異性用于裂縫儲層的檢測和評價,并研究了裂縫儲層隨方位變化的規律.由于可以通過分析縱波反射率隨方位角和偏移距的變化來確定裂縫方向和裂縫密度,引起了人們的廣泛關注(Rüger,1998;Hall and Kendall,2003;劉喜武等,2015;Xie et al.,2020;Ma et al.,2022).然而,由于傳統的地震各向異性對裂縫長度不敏感(Hudson,1981;Liu等,2000),因此,從縱波反射中估計裂縫大小仍然很具有挑戰性.

本文研究目的是基于已提出的黏彈性Chapman-Kelvin的動態等效介質模型,深入研究在黏彈性裂縫儲層中地震波的傳播特性,分析裂縫儲層參數(主要為裂縫密度、裂縫長度、孔隙度和含水飽和度)對地震波傳播特性的影響,以便為復雜裂縫的儲層參數反演問題研究提供有力的理論依據.我們考慮到傳統的等效介質模型未充分考慮非完全彈性介質理論和頻變各向異性理論的雙相或多相流體假設,也無法對地震波在實際裂縫儲層中的頻散和衰減給出準確合理的解釋.基于黏彈性Chapman-Kelvin的動態等效介質模型(廖建平等,2023a;Liao et al.,2023a,b),我們將關于黏彈性、噴射流以及斑塊效應耦合的頻變彈性系數結合到Schoenberg和Protázio(1992)概括的Zoeppritz方程,計算出反射系數,同時考慮到發生地震頻散時,反射系數和頻率產生關系,構建了在角度、方位和時間域內的新型正演方程.最后,基于PP波頻變反射系數對裂縫儲層參數變化有較好的敏感性,開展了兩種反演方法研究:其一是基于貝葉斯理論直接反演方法;其二是基于MCMC的隨機反演方法.通過頻變AVAZ進行儲層參數的反演研究,結果表明,MCMC隨機反演方法在缺失先驗的儲層參數信息時,反演結果的不確定性較強.當存在有效且足夠的先驗信息時,反演結果的可靠性進一步提升.而基于貝葉斯理論直接反演方法,由于在方位角AVO分析中對地震頻散的考慮,則顯示出區別大尺度裂縫和微尺度裂隙的潛力,進而說明基于頻變反射系數進行儲層參數反演的可行性.

1 裂縫儲層的等效介質建模

1.1 基于黏彈性Chapman-Kelvin動態等效介質模型

為了描述實際的裂縫儲層介質在全頻帶范圍內的地震響應規律,Chapman等(2002)提出一種微觀結構的噴射理論模型,該理論模型在低頻時符合Gassmann方程,在高頻時符合與噴射相關的頻散,并且描述了在含有微裂隙和孔隙的巖石中,流體在微裂隙和孔隙之間,以及微裂隙之間的流動.Chapman(2003)認為,頻變各向異性理論的有效性不僅取決于裂縫尺度,還取決于顆粒尺度的流體效應.所以,在基于Chapman等(2002)的噴射理論基礎上,描述了流體流動引起的噴射流產生的各向異性頻散效應,引入了一套垂直排列的介觀尺度裂縫.該模型考慮了介觀和微觀兩種尺度的流體流動,能夠對地震頻段地震波的頻散和衰減給出合理的解釋.為了準確描述模型,需要考慮一系列參數,包括孔隙度、裂縫密度、流體參數、裂縫參數以及時間尺度等.Chapman(2003)模型中存在橢球狀微裂隙、球形孔隙和定向排列的橢球狀介觀裂縫.圖1為Chapman(2003)模型的示意圖,其中微裂隙、孔隙與巖石顆粒均為微觀尺度,而裂縫為介觀尺度.

圖1 Chapman裂縫孔隙介質示意

考慮流體在整個孔隙空間內交換,Chapman(2003)模型的頻變剛度矩陣形式為:

(1)

(2)

式(2)的矩陣形式,如式(3)所示:

(3)

其中,C′11、C′33、C′44、C′23、C′13、C′66為兩種非混相流體飽和裂縫巖石的有關噴射流和斑塊效應的頻變各向異性彈性剛度系數,且黏彈性Chapman-Kelvin模型的等效剛度張量在本文的研究中主要取決于以下參數:頻率、裂縫密度、裂縫長度、孔隙度和含水飽和度.

1.2 等效黏彈性剛度張量

五個獨立的等效黏彈性剛度系數C11、C13、C33、C44、C66是儲層參數的函數,儲層參數是由基于黏彈性Chapman-Kelvin動態等效介質模型給出.我們對基于Chapman-Kelvin模型的等效黏彈性剛度參數的響應特征進行分析,這里的剛度參數分別是頻率f、裂縫密度d、裂縫長度l、孔隙度por和含水飽和度sw的函數.圖2為各剛度參數的實部曲線和虛部曲線,其中設定裂縫密度為0.05和孔隙度為0.1,可以觀察到剛度張量C11、C13、C33以及C66的實部分量隨頻率的增大,先增大,當達到高頻104Hz以后,逐漸趨于平緩,達到最大值(圖2a、c、e、i),而C44的實部分量隨頻率逐漸減小,當達到高頻104Hz以后,逐漸趨于平緩,達到最小值(圖2g);剛度張量C11、C13的虛部分量表現為先隨頻率的增大而增大,但在100Hz以后出現了顯著的下降現象,如圖2b、d所示.觀察到地震頻段為0到100 Hz的頻散和衰減現象,其中較為明顯的是較大尺度裂縫,即l為1 m和較高含水飽和度,即sw為0.6,也反映出基于黏彈性Chapman-Kelvin模型的頻變AVAZ裂縫儲層參數反演,與裂縫尺寸和含水飽和度關系不唯一.

圖2 剛度參數的實部和虛部隨頻率的變化

1.3 基于黏彈性Chapman-Kelvin模型的裂縫儲層介質地震波頻變特征的分析

裂縫儲層中地震波的頻散、衰減和頻變各向異性,以及分析出重要的裂縫儲層參數對地震波在裂縫介質中傳播特性的影響,將為裂縫儲層的儲層參數反演提供有力的理論基礎.在基于黏彈性Chapman-Kelvin模型中,我們研究了裂縫密度、裂縫長度、孔隙度、含水飽和度對地震波頻變特征的影響.試驗中考慮了一塊被水和CO2氣體飽和的巖石,體積模量和剪切模量代表不存在裂縫時的測量值,表1給出了各個參數的測量值.其中,附錄A中給出了橫波速度、橫波分裂、橫波衰減和準縱波速度的計算公式.

表1 模型計算中輸入的參數值

裂縫密度是裂縫儲層中的一個重要參數,它描述了儲層中裂縫發育的程度.圖3顯示了裂縫密度對準縱波速度的影響.首先分別取頻率為1 Hz、10 Hz、100 Hz、1000 Hz,從圖3a可以看出:準縱波速度隨裂縫密度的增大而增大.然后分別取裂縫密度為0.02、0.04、0.06和0.08,研究了在不同的裂縫密度下準縱波速度與頻率的關系,從圖3b可以看出,準縱波速度隨頻率的增大,先增大,后緩慢減小;準縱波產生的地震波速度頻散的頻率范圍與裂縫密度的變化無關,而頻散的幅度隨著裂縫密度的增大而增大.

圖3 準縱波速度隨裂縫密度(a)和頻率(b)的變化

裂縫與孔隙之間流體流動的弛豫時間τf與裂縫長度l成正比關系,因而裂縫長度將會影響頻散出現的頻段(劉喜武等,2015).圖4顯示了裂縫長度對橫波分裂的影響.首先從圖4a可以看出,橫波分裂隨裂縫長度的增大,而非線性增大,其中在低頻即f=1 Hz和10 Hz時,變化的幅度小,但是在高頻即f=100 Hz和1000 Hz時,變化的幅度變大.然后分別取裂縫長度為0.15 m、0.25 m、0.35 m和0.45 m,研究了不同的裂縫長度下橫波分裂與頻率的關系,從圖4b可以看出,橫波分裂產生的速度頻散的頻率范圍不隨裂縫長度變化,而頻散的幅度隨著裂縫長度的增大而增大.

圖4 橫波分裂隨裂縫長度(a)和頻率(b)的變化

圖5顯示了孔隙度對橫波速度的影響.從圖5a中,可以看出,橫波速度隨孔隙度的增大而增大,還可以發現橫波速度隨頻率變化的敏感性較低.分別取孔隙度為0.05、0.1、0.15和0.2,研究了不同的孔隙度下橫波速度與頻率的關系,從圖5b可以看出,橫波產生的速度頻散的頻率范圍不隨孔隙度變化,而頻散的幅度隨著孔隙度的增大而增大.

圖5 橫波速度隨孔隙度(a)和頻率(b)的變化

圖6顯示了含水飽和度對橫波衰減的影響.從圖6a中研究發現,橫波衰減在頻率為1000 Hz時,隨著含水飽和度的增大而減小,在其他頻率下(1 Hz、10 Hz、1000 Hz),橫波衰減隨含水飽和度的增大,先增大后減小.為了研究不同含水飽和度下橫波衰減隨頻率的變化,分別取含水飽和度為0、0.25、0.5、0.75和1,從圖6b可以看出,含水飽和度為零時,橫波衰減隨頻率先增大后減小;在含水飽和度不為零時,橫波產生衰減的頻率范圍也與含水飽和度的變化無關,而衰減的幅度隨含水飽和度的增大而增大.

2 基于黏彈性Chapman-Kelvin模型的AVAZ正問題與分析

可以利用地震屬性檢測裂縫導致的巖石彈性性質變化,如振幅隨方位角和入射角的變化、地震波傳播速度隨方位的變化等(Liu et al.,2012;Zhang,2020).概況地說,目前可以分為兩類,其一是根據AVO屬性隨方位角的變化進行橢圓擬合,定性估計裂縫發育方位和強度;另一種則是根據相應的裂縫等效介質理論,將反射的地震波振幅與偏移距和方位角聯系起來,從而建立出反射系數與裂縫儲層參數的關系,進而反演出裂縫儲層參數(劉喜武等,2015).上述研究均是基于靜態等效介質理論,其剛度系數矩陣與頻率無關.并且在基于靜態等效介質的AVO或AVAZ分析中,忽略了頻散和衰減的影響.而動態等效介質模型可以描述具有頻變特性的裂縫介質,如在1.1小節中的基于黏彈性Chapman-Kelvin動態等效介質模型.速度頻散和衰減的出現將使得上覆地層與裂縫儲層分界面的反射系數具有頻變性,這種頻變性與裂縫儲層參數密切相關.因此,基于動態等效介質理論計算的頻變反射系數與基于靜態等效介質理論計算的不具有頻變的反射系數相比,具有更為接近實際反射系數的值,即包含了更多裂縫儲層信息.

本部分采用黏彈性Chapman-Kelvin動態等效介質模型描述裂縫孔隙介質,結合Schoenberg和Protázio(1992)概括的Zoeppritz方程,給出了頻變反射系數的計算方法,并分析了重要的儲層參數與頻變反射系數的關系.最后,給出了正演方程的建立過程.

2.1 PP波和PS波反射系數的計算與頻變AVAZ分析

R=(Y′-1Y+X′-1X)-1(Y′-1Y-X′-1X)

(4)

從而得到PP波和PS波的反射系數分別為R(1,1)和R(2,1),這個反射系數是入射角、方位角、頻率以及儲層參數矢量m的函數:

m=[l,fd,por,sw],

(5)

其中,l為裂縫長度;fd為裂縫密度;por為孔隙度;sw為含水飽和度.

為了研究裂縫孔隙介質與上覆介質分界面處的PP波和PS波反射系數具有頻變性的同時,也隨著入射角和方位角的變化,我們設計一個兩層介質模型,即上層為頁巖,并且假定為各向同性;下層為帶有垂直裂縫排列的碳酸鹽巖(各向異性,被80%水和20%CO2所飽和,體積模量、剪切模量等其他參數與表1中的參數一致),該模量的具體參數值,如表2所示.

表2 水和氣飽和模型的參數

采用表1和表2中的參數,通過式(4)計算PP波和PS波的反射系數,并分析PP波和PS波的頻變AVAZ特性.

圖7給出了不同入射角情況下,PP波、PS波反射系數隨方位角和頻率的變化.當頻率一定時,從圖7a、b可以看出,PP波和PS波的反射系數幅值均隨方位角的變化,而且這種變化隨著入射角的增大變得更為顯著;當固定方位角為50°,從圖7c、d我們可以發現PP波反射系數幅值隨頻率的增大,先減小后增大;且PP波反射系數幅值在地震頻段(0～100 Hz)內隨頻率的變化最為顯著,尤其在大角度(40°)入射情況下;PS波反射系數幅值隨頻率的增大,先增大后減小;PS波反射系數幅值隨入射角的增大而增大.

圖7 不同入射角情況下PP波和PS波反射系數隨頻率和方位角的變化

2.2 儲層參數變化對PP波頻變反射系數的影響

研究裂縫儲層參數(裂縫密度fd、裂縫長度l、孔隙度por和含水飽和度sw)的變化對PP波頻變反射系數的影響.分析PP波頻變AVO特性對這些儲層參數的敏感性,以便于我們探究利用頻變AVO方法反演這些參數的可行性.其中,選擇表2所示的兩層模型,并利用式(4)計算反射系數.

圖8顯示了裂縫密度變化對PP波頻變反射系數的影響.分別取裂縫密度為0.02、0.05、0.08和0.1,其他模型參數保持不變.不同裂縫密度情況下,反射系數幅值均隨著入射角的增大而增大,在入射角大于20°的范圍內,這一現象較為顯著;反射系數幅值隨著裂縫密度的增大而增大.反映出PP波頻變AVO特性對裂縫密度的變化較為敏感.

圖8 不同裂縫密度的PP波反射系數幅值隨入射角和頻率的變化

圖9顯示了裂縫長度變化對PP波頻變反射系數的影響.分別取裂縫長度為0.2 m、0.5 m、1 m和2 m,其他模型參數保持不變.我們可以發現頻變反射系數的變化主要在入射角30°到45°之間,隨著裂縫長度的增大,反射系數幅值隨頻率的變化幅度減小,而不同頻率的反射系數幅值均有所增大(如圖9所示);裂縫長度的增大導致不同入射角的反射系數幅值均有所增大,而反射系數幅值隨入射角的變化幅度減小(如圖9所示).反映出PP波頻變AVO特性對裂縫長度的變化具有一定的敏感性,尤其在大角度即30°到45°、低頻即10 Hz到40 Hz情況下,但其敏感性不如裂縫密度.

圖9 不同裂縫長度的PP波反射系數幅值隨入射角和頻率的變化

圖10顯示了PP波頻變反射系數隨孔隙度的變化.分別取孔隙度為0.05、0.1、0.15和0.2,其他模型參數保持不變.從圖10中可以看出,孔隙度的變化引起PP波頻變AVO特性的顯著變化,隨著孔隙度的增大,不同頻率、不同入射角對應的反射系數幅值均明顯減小.從而說明,PP波頻變AVO特性對孔隙度的變化非常敏感.因而,可以利用頻變AVO反演技術較為準確地反演出孔隙度.

圖10 不同孔隙度的PP波反射系數幅值隨入射角和頻率的變化

圖11顯示了PP波頻變反射系數隨含水飽和度的變化.分別取含水飽和度為0.2、0.4和0.6,其他模型參數保持不變.從圖11可以看出,只有在低頻時,反射系數幅值隨著入射角的增大而增大,并且隨著含水飽和度的增大,反射系數的幅值在不斷地變小.所以,PP波頻變AVO特性對含水飽和度的變化具有一定的敏感性,但要遠遠弱于裂縫密度、裂縫長度和孔隙度.

圖11 (a) sw=0.2; (b) sw=0.4; (c) sw=0.6; (d) 不同含水飽和度的PP波反射系數幅值隨入射角和頻率的變化

2.3 正演方程的建立

傳統的合成地震記錄是通過震源與時間序列的反射率的褶積來生成,但是當地震頻散發生時,反射系數與頻率產生關系,使得我們在時間域上難以應用褶積.于是,我們考慮頻率域中的地震子波(雷克子波)和傅里葉反變換方法來合成疊前地震記錄,模型參數m和疊前地震數據的關系為:

d(θ,φ,t)=F-1[W(f)·R(θ,φ,f;m)]+n,

(6)

式(6)可以簡化為:

d=F-1G(m)+n,

(7)

其中,dT=(d1,d2,…,dM)T為疊前觀測地震數據;mT=(m1,m2,…,ml)T為模型參數;F-1表示傅里葉反變換;G為正演算子;W(f)為頻率域中的地震子波;R為反射系數,如式(4)所示;f為頻率;θ為入射角;φ為方位角;n為觀測地震記錄中的噪聲.

3 裂縫孔隙介質中儲層參數的反演

裂縫儲層參數反演是裂縫儲層評價的重要內容,對于裂縫油氣藏的勘探和開發有著重要的意義.用于描述裂縫儲層介質的傳統裂縫等效介質理論是基于靜態等效介質.由于沒有考慮裂縫儲層介質彈性響應的頻變性,在含有少量較大尺寸裂縫的介質與含有較多小尺寸微裂隙的介質中,它們的彈性常數和各向異性參數可能相等,所以靜態等效介質理論無法區別出地震各向異性是由微裂隙引起,還是由較大尺寸的裂縫引起(Maultzsch et al.,2003;Al-Harrasi et al.,2011).與微裂隙相比,一般來說,大尺度裂縫在裂縫儲層中主導著流體流動.因此,裂縫尺寸是裂縫儲層的一個重要參數,也是儲層評價標準的一個重要參考.又根據2.2小節可以得知,動態等效介質模型所描述的地震響應的頻變特征對裂縫密度、裂縫長度、孔隙度、含水飽和度等裂縫儲層參數具有很好的敏感性.

在已知裂縫走向的條件下,利用貝葉斯理論,通過疊前頻變AVAZ反演得到裂縫儲層參數,即裂縫密度、裂縫長度、孔隙度和含水飽和度,以及顯示出區別大尺度裂縫和微尺度裂隙的潛力.

3.1 貝葉斯理論

式(7)所示的反問題在地球物理領域中總是不適定的非線性問題,這種不適定性是由噪聲干擾或不適當的正演算子等因素引起的(黃廣譚,2019).因此,求解相應的反問題需要引入盡可能多的先驗約束信息來緩解這種不適定性,這通常被稱為正則化,本文我們采用貝葉斯理論來減輕這種不適定性問題.

貝葉斯反演基于概率和統計的思想,即利用先驗信息和似然函數獲得后驗信息,從而把求取反演目標函數的極小值問題,轉化為獲得模型參數的后驗概率分布的最大化問題.貝葉斯條件概率的公式為:

(8)

P(m|d)∝L(d|m)P(m).

(9)

3.2 裂縫孔隙介質的儲層參數反演

基于統計理論的貝葉斯反演,通過引入模型參數的先驗信息作為反演的正則化約束條件,能夠有效地降低反演的非唯一性.模型參數的先驗分布可以服從不同的分布函數(如均勻分布、高斯分布和柯西分布等),進而可以獲得不同的反演結果(黃廣譚,2019).本文選用多元高斯分布來構建先驗模型,則有:

(10)

其中,Nd為模型參數的個數;m=[l,fd,por,sw]T～P(μm,Cm),l為裂縫長度,fd為裂縫密度,por為孔隙度,μm為模型參數的期望,即μm=E{m}=[μl,μfd,μpor,μsw]T;Cm為模型參數m的協方差矩陣,如式(11)所示:

(11)

其中,Cov表示協方差.

似然函數,即失配函數,是指由模型向數據的映射,描述了觀測數據與合成數據之間的誤差,表3給出了5種常用的失配函數(Huang et al.,2020).

表3 5種常用的失配函數

L1范數和L2范數度量作為地震反演中最常用的兩種方法,常被用作失配函數.但是這兩種方法都有其局限性.L1范數方法在零點處存在奇異性問題,L2范數方法因收斂性差而受到限制(黃廣譚,2019).而Huber范數、混合范數和對數絕對范數能有效地克服常規范數的缺陷,其中對數絕對范數最為突出,對數絕對范數失配函數既能解決L1范數零點處的奇異性問題,又能盡可能保持與L1范數的收斂速度.

所以,我們在計算觀測數據d(疊前地震數據)與正演模擬響應的殘差時,引入了對數絕對范數,則殘差表示為:

(12)

我們從殘差中計算出似然函數,則似然函數表示為:

P(d|m)=aexp(-bΔE),

(13)

其中,a為歸一化系數;b為常數.

式(13)是由Ulrych等(2001)分析得出的;Ulrych等(2001)和Mavko和Mukerji(1998)從理論上解釋b的確定如何與觀測的標準誤差有關.然而Kirkpatrick等(1983)解釋式(13)作為一個經驗關系.

最終的似然函數表示為:

(14)

略去無關的歸一化常量,可以得到m的后驗分布如下:

P(m|d)∝exp[-b·[‖d-G(m)‖1

(15)

因此,目標函數J(m)可以簡化為:

(16)

在得到后驗概率分布之后,為了量化與裂縫相關的反演模型參數中的不確定性,可以通過蒙特卡羅隨機采樣方法得到反演結果.蒙特卡羅隨機采樣是通過大量重復試驗得到某一事件的出現概率.在地震勘探領域中,一般認為地層參數具有馬爾科夫性質,即在某一時刻t的模型參數的概率分布只與前一時刻(t-1)的模型參數有關.馬爾科夫鏈是不同狀態之間的轉移,概率矩陣P決定轉移的過程.其中,不同的轉移矩陣的構建方法將產生不同的MCMC方法(Ali and Jakobsen,2014).因此,一般將蒙特卡羅隨機采樣方法和馬爾科夫鏈結合起來進行采樣.

為了有效地避免陷入局部極小的情況,在貝葉斯理論下,采用Metropolis準則對后驗概率密度分布進行采樣.Metropolis準則對任何樣本都是以一定的概率接受或拒絕(Rubino et al.,2013).假設當前所在點為i,隨機選擇的下一個點為j,由i向j傳遞所遵循的Metropolis采樣的接收準則為:

(17)

4 數值試驗

4.1 基于貝葉斯理論的反演試驗

本文使用表4中的參數進行合成研究,并關注于裂縫長度、裂縫密度、孔隙度和含水飽和度的估計以及在方位角AVO分析中,能夠通過對地震頻散的考慮,區分出大尺度裂縫與微尺度裂隙.本文的合成地震記錄考慮一個具有兩個界面的例子,即上層為砂巖,并且假定為各向同性;中間層為帶有垂直排列裂縫的碳酸鹽巖被60%的水和40%的CO2飽和,體積模量、剪切模量等其他參數與表1中的參數一致;下層為砂巖,并且假定為各向同性介質,具體的參數設定如表4所示.本文考慮了兩種裂縫的情況,第一種是微尺度裂隙,即裂隙長度為1 mm,第二種是大尺度裂縫,即裂縫長度為1 m.利用公式(7)合成了角度域和方位角域的兩種裂縫情況的地震道集記錄,并在道集上加入了信噪比為4的隨機噪聲,信噪比SNR的定義為(Chen and Fomel,2015):

表4 水和氣體飽和模型的參數 (合成示例)

(18)

合成的角度域和方位域的地震道集,如圖12所示.其中標注框部分為考慮微尺度裂隙和考慮大尺度裂縫所合成的地震道集中主要的不同部分.

按照3.2小節中的反演步驟,本文首先假設裂縫長度l、裂縫密度fd、孔隙度por和含水飽和度sw滿足高斯分布,均值分別取5、0.04、0.12、0.5,考慮到裂縫長度、裂縫密度、孔隙度以及含水飽和度在不同地區隨巖性變化較大,因此選取較大的方差,即分別取2.5、0.25、0.4、0.35,并將裂縫長度、裂縫密度、孔隙度和含水飽和度控制在一定范圍內,即l∈[0,10],fd∈[0,0.1],por∈[0,0.2],sw∈[0,1],先驗概率分布函數如圖13所示,并假定所有其他參數都已知.在實際應用中,這些先驗信息可以從測井分析或地震解釋中得到(Aster et al.,2005).

通過使用式(12)求得觀測數據和正演模擬數據的殘差,再通過式(14)得到似然函數.然后根據式(15)計算出后驗概率分布P(fd,l,por,sw|d).對于給定問題的概率分布P(x,y,z),每個變量的邊際概率分布為:

(19)

(20)

(21)

因此,可以從邊際概率的角度出發,直接從后驗分布P(fd,l,por,sw|d)中,計算出的裂縫密度fd、裂縫長度l、孔隙度por、含水飽和度sw的邊際概率分布的表達式,分別為:

(22)

(23)

(24)

(25)

我們分別應用于大尺度裂縫和微尺度裂隙情況,即微尺度裂隙,取裂隙密度為0.05,裂隙長度為1 mm,并且在該種情況下孔隙度取0.1,含水飽和度取0.6;大尺度裂縫,取裂縫密度為0.05,裂縫長度為1 m,并且在該種情況下孔隙度取0.1,含水飽和度取0.6.圖14、圖15、圖16、圖17分別為后驗概率P(fd,l,por,sw|d)的裂縫密度、裂縫長度、含水飽和度和孔隙度的邊際概率分布.從圖14的裂縫密度邊際概率分布來看,裂縫密度在大尺度裂縫下反演結果的最大概率分布集中在0.05處,與真實的裂縫密度很吻合,而在微尺度裂隙下的反演結果的最大概率分布集中在0.04,很顯然與真實的裂縫密度相差很大.從圖15的裂縫長度邊際概率分布來看,微尺度裂隙下裂縫長度的最大概率主要集中在10-4m和10-2m之間,不能準確地估計出真實微尺度裂隙的長度,而大尺度裂縫的最大概率在100m處,與真實大尺度裂縫的長度相符合.從圖16的含水飽和度邊際概率分布來看,大尺度裂縫下含水飽和度的最大概率主要集中在0.6處,與真實含水飽和度值相吻合,而微尺度裂隙下含水飽和度的最大概率集中在0.82處,很顯然不符合真實值.最后,從圖17的孔隙度邊際概率分布來看,大尺度裂縫下孔隙度的最大概率分布集中在0.1和0.12區間,而微尺度裂隙下孔隙度的最大概率集中在0.08處,相對大尺度裂縫下的結果來說,精確度不夠.

圖14 微尺度裂隙(a)、大尺度裂縫(b)下裂縫密度的后驗概率分布

圖15 微尺度裂隙(a)、大尺度裂縫(b)下裂縫長度的后驗概率分布

圖16 微尺度裂隙(a)、大尺度裂縫(b)下含水飽和度的后驗概率分布

圖17 微尺度裂隙(a)、大尺度裂縫(b)下孔隙度的后驗概率分布

通過求解裂縫密度、裂縫長度、含水飽和度和孔隙度等儲層參數的后驗概率分布,反演出各個參數最大概率的分布情況,可以很好地區分開大尺度裂縫和微尺度裂隙兩種情況.也可以發現,我們可以精確地反演出大尺度裂縫的長度,而不能夠準確地恢復微尺度裂隙的長度,這是因為地震頻散受到了裂縫長度的控制,當裂縫長度足夠小時,即微尺度裂隙,地震頻散可以忽略不計,所以對于地震頻散可忽略不計的微尺度裂隙情況,不能準確地反演出裂縫密度、含水飽和度和孔隙度.因此,在方位角AVO分析中對地震頻散的考慮,顯示出將大尺度裂縫和微尺度裂隙區分開的潛力.

4.2 基于頻變反射系數的MCMC反演試驗

基于公式(7)的正演方程,若不考慮地震子波的作用,直接利用反射系數,并增加約束,對儲層參數進行反演.我們將儲層參數集M設置為:

M=[fd,l,por,sw],

(26)

即我們關注的儲層參數主要為裂縫密度fd、裂縫長度l、孔隙度por和含水飽和度sw.因此,隨頻率變化的P波反射系數可以表示成:

R(θ,φ,f;M),

(27)

其中,θ和φ分別為入射角和方位角,f為頻率.

對于K個入射角(k=1,2,…,K),Z個方位角(z=1,2,…,Z),F個頻率(f=1,2,…,F),反射系數R可以寫作為:

Rkzf(M)=R(θk,φz,ff;M),

(28)

將反射系數進一步寫成矩陣形式:

R=GM,

(29)

其中,G為基于Chapman-Kelvin模型和Schoenberg和Protázio(1992)概括的Zoeppritz方程所計算的頻變反射系數的核函數.

基于表4中的模型參數,且此時考慮大尺度裂縫的情況,即我們假定真實裂縫密度為0.05,真實裂縫長度為1 m,真實孔隙度為0.1,以及真實含水飽和度為0.8.然后,我們嘗試使用計算得到的合成地震AVAZ資料來估計真實值,即反射系數的大小作為不同入射角的地震頻率(≤100 Hz)和方位角的函數.使用1.1小節中描述的黏彈性Chapman-Kelvin模型和附錄B中的阻抗矩陣(Schoenberg and Protázio,1992)來計算反射系數的大小.假設入射角為10°、20°、30°、40°,方位角為20°、40°、60°、80°,頻率為10 Hz、20 Hz、30 Hz、40 Hz、50 Hz、60 Hz、70 Hz、80 Hz、90 Hz、100 Hz.當缺少先驗信息時,式(16)中的目標函數J(m)變為:

(30)

圖18 不同的入射角下添加信噪比為5的隨機噪聲的PP波反射系數

圖19顯示了在沒有關于任何儲層參數的先驗信息的情況下,使用頻變的地震AVAZ數據,以儲層參數的后驗概率密度函數形式得到的MCMC采樣的反演結果.從圖19可以看出反演結果的不確定性很強,無法得到儲層參數的準確解.

圖19 無先驗信息的儲層參數的反演結果

我們接下來通過假設一個關于裂縫密度的先驗信息來進行反演.本文我們將裂縫密度、裂縫長度、孔隙度和含水飽和度的目標函數J(m),以及關于裂縫密度的先驗信息定義為:

(31)

其中,第一個平方項是在4個不同的入射角、4個不同的方位角和10個頻率(≤100 Hz)上反射系數的觀測值和計算值之間的不匹配;第二個平方項表示裂縫密度的先驗信息,fd0為裂縫密度先驗信息的均值;Δfd為裂縫密度先驗分布的標準差.

圖20顯示了利用頻變的地震AVAZ數據和裂縫密度的先驗信息對儲層參數的后驗概率密度函數形式的MCMC反演結果.我們可以清楚地看到,包含裂縫密度的先驗信息可以在一定程度上提高裂縫長度、孔隙度和含水飽和度的確定性,同時裂縫密度的反演結果的不確定性也大大減小,但含水飽和度的確定性仍然較差.

同樣,假定已知裂縫長度和孔隙度的先驗信息.可以通過觀測和鉆孔數據的插值得到裂縫長度的先驗信息.Odling等(1999)討論了一種從露頭工區中求取裂縫長度的方法,在這種方法中,裂縫長度總體繪制為歸一化的累積頻率分布,即T(l),其中l為裂縫長度.當使用對數軸時,指數為c,裂縫長度的分布如下:

T(l)∝l-c.

(32)

Maultzsch等(2003)通過使用頻率相關的剪切波各向異性和P波衰減的反演也可以獲得關于裂縫長度的先驗信息.

我們將裂縫密度、裂縫長度、孔隙度和含水飽和度的目標函數J(m),以及關于裂縫長度和孔隙度的先驗信息定義為:

(33)

其中,第一個平方項與式(31)中的定義一致;第二個平方項分別表示裂縫長度和孔隙度的先驗信息,l0、por0分別為裂縫長度和孔隙度先驗信息的均值;Δl、Δpor分別為裂縫長度和孔隙度先驗分布的標準差.

圖21顯示了利用頻變地震AVAZ數據和裂縫長度和孔隙度的先驗信息,對儲層參數的后驗概率密度函數形式的MCMC反演結果.從圖21a的裂縫密度反演結果來看,相對于沒有考慮先驗信息的圖19a的裂縫密度反演的準確度有了進一步的提高.但是我們對比圖20d、d,發現關于裂縫長度和孔隙度的先驗信息對含水飽和度的估計沒有太大的改善作用.

圖21 具有裂縫長度和孔隙度先驗信息的儲層參數的反演結果

最后,我們針對具有裂縫長度和孔隙度先驗信息的含水飽和度反演結果確定性仍然不是很理想,所以我們考慮了足夠多的先驗信息,即假定裂縫長度、裂縫密度和孔隙度的先驗信息是已知的.本文將裂縫密度、裂縫長度、孔隙度和含水飽和度的目標函數J(m),以及關于裂縫長度、裂縫密度和孔隙度的先驗信息定義為:

(34)

其中,第一個平方項與式(31)中的定義一致;第二個平方項分別表示裂縫長度、裂縫密度和孔隙度的先驗信息,l0、fd0、por0分別為裂縫長度、裂縫密度和孔隙度先驗信息的均值;Δl、Δfd、Δpor分別為裂縫長度、裂縫密度和孔隙度先驗分布的標準差.

圖22顯示了利用頻變地震AVAZ數據和裂縫長度、裂縫密度和孔隙度的先驗信息對儲層參數的后驗概率密度函數形式的MCMC反演結果.從圖22d含水飽和度的反演結果來看,相比之前,反演結果的確定性有了一定的加強.并且通過對比圖21和圖22中含水飽和度的反演結果,我們也不難發現關于裂縫密度、裂縫長度、孔隙度三者的先驗信息最終對含水飽和度的估計變得更加關鍵.

圖22 具有裂縫長度、裂縫密度和孔隙度先驗信息的儲層參數的反演結果

5 討論

裂縫等效介質模型能夠更精確地模擬裂縫儲層介質的性質,本文根據所提出的基于黏彈性Chapman-Kelvin動態等效介質模型,其能夠適應復雜的儲層情況,并初步討論了裂縫儲層的儲層參數反演問題,但是離實際裂縫儲層中工業應用還有很遠的距離.綜合國內外對裂縫儲層相關問題的研究現狀,對其進一步的研究應包括以下幾個方面:

(1)完善裂縫儲層動態等效介質模型理論.在Chapman(2003)模型的基礎上,考慮了實際地下介質的黏彈性和地震頻散,引入了線性流變Kelvin黏彈性理論模型,構建出黏彈性Chapman-Kelvin動態等效介質模型.但是能否將線性流變理論模型與基于Boltzmann理論的黏彈性模型統一結合到Chapman(2003)模型中,以便適應更加復雜的黏彈性儲層情況還有待后續研究.

(2)裂縫儲層參數反演由于各向異性、地震資料的處理、參數的選擇以及反演策略的原因具有很強的多解性,尤其在考慮地震頻散的情況下,表現得更為顯著.本文基于貝葉斯理論進行裂縫儲層參數的反演,在缺失先驗的儲層參數信息時,反演結果的不確定性較強.當存在有效且足夠的先驗信息時,反演結果的可靠性進一步提升.但是,關于儲層參數的先驗信息很難獲取.

(3)本文所提出的頻變AVAZ反演技術,為裂縫儲層參數的反演和區別大尺度裂縫和微尺度裂隙提供了新的思路.在實際應用中,需要寬方位地震采集數據和高保真去噪技術來獲得高質量的疊前數據,采用多尺度地震分頻技術,按照本文的反演流程進行大尺度裂縫和微尺度裂隙的區分,以及裂縫儲層參數的反演,并結合實際的測井資料,驗證所提出的反演方法的實際應用效果.

6 結論

裂縫儲層介質的等效介質模型建立了地震數據與儲層特征的聯系.考慮到黏彈性介質比彈性介質更能代表地球的性質,以及傳統的裂縫等效介質模型無法對實際裂縫儲層中地震波的頻散和衰減給出合理的解釋.因此,本文基于已經提出的黏彈性Chapman-Kelvin動態等效介質模型,并在此基礎上就裂縫密度、裂縫長度、孔隙度以及含水飽和度對地震波頻變特征的影響做了分析.通過理論分析和數值試驗,我們可得到以下結論和認識:

(1)利用黏彈性Chapman-Kelvin動態等效介質模型與Schoenberg和Protazio概括的Zoeppritz方程計算了反射系數,分析了反射PP波和PS波的頻變AVAZ特性,即當頻率一定時,PP波和PS波的反射系數幅值均隨方位角的變化,而且這種變化隨著入射角的增大變得更為顯著;當固定方位角為50°,我們可以發現PP波反射系數幅值隨頻率的增大,先減小后增大;且PP波反射系數幅值在地震頻段(0～100 Hz)內隨頻率的變化最為顯著,尤其在大角度入射情況下;PS波反射系數幅值隨頻率的增大,先增大后減小,而隨入射角的增大而增大.并就重要的儲層參數對PP波頻變反射系數的影響進行了分析,發現孔隙度的敏感性最強,其他依次為裂縫密度、裂縫長度和含水飽和度.

(2)基于PP波頻變反射系數對裂縫儲層參數變化有較好的敏感性,開展了兩種反演方法研究,其一是基于貝葉斯理論直接反演方法,其中以對數絕對范數作為似然函數和高斯分布即L2范數度量作為先驗約束;其二是基于MCMC的隨機反演方法.通過頻變AVAZ進行儲層參數的反演研究,結果表明,MCMC隨機反演方法在缺失先驗信息時,反演解的不確定性較強;包含裂縫密度的先驗信息在一定程度上提高裂縫長度、孔隙度和含水飽和度的確定性,但含水飽和度的確定性仍然很差;假定裂縫長度、裂縫密度和孔隙度三者先驗信息是已知時,含水飽和度反演的確定性有了進一步的提高.而基于貝葉斯理論直接反演方法,由于考慮到地震頻散的影響,通過后驗概率的邊際概率分布情況,區別出了大尺度裂縫和微尺度裂隙,使得利用頻變方位AVO來區分微尺度裂隙和大尺度裂縫成為可能.

致謝感謝審稿專家提出的寶貴意見.

附錄A 基于黏彈性Chapman-Kelvin模型的橫波速度、橫波分裂、橫波衰減和準縱波速度的計算

基于黏彈性Chapman-Kelvin地震波傳播的模型,如式(A1)所示:

(A1)

式(A1)寫成矩陣形式為:

(A2)

其中,子塊矩陣:

(A3)

mfw=C44+Cmu,

(A4)

其中:

Cmu=2Udry-real(C44),

(A5)

其中,Udry為干燥巖石的剪切模量;real指復數的實部部分.

則橫波速度vs的表達式為:

(A6)

其中,r為巖石密度,橫波衰減qis的表達式為:

(A7)

其中,imag指復數的虛部部分,橫波分裂sws的表達式為:

(A8)

其中,vps為純橫波的速度,表達式為:

(A9)

其中,α為地震波傳播和對稱軸之間的角度(0°對應于沿裂縫的法線傳播);且vqs為準橫波速度,其表達式為:

(A10)

其中,H為常數,表達式為:

H=[(C11-C44)sin2α-(C33-C44)cos2α]2

+(C13+C44)2sin2(2α).

(A11)

最后,準縱波速度vqp的計算表達式為:

(A12)

附錄B 基于Schoenberg和Protázio(1992)阻抗矩陣的計算

Schoenberg和Protázio(1992)利用Zoeppritz方程系數矩陣的子矩陣給出了平面波反射和透射問題的顯示解;Chapman和Liu(2003)在Schoenberg和Protázio(1992)的研究基礎上,推導了與入射角和方位角相關的反射率;Jin等(2017)假設各向異性介質有一個平行于X1-X2平面的水平對稱平面,設X3=0是兩個半空間之間的水平界面.設定上部彈性介質中入射波的入射角為θ,方位角為φ,傳播方向可以表示為n=(sinθcosφ,sinθsinφ,cosθ).當地震波在垂直X1-X3平面內傳播時,位移可以表示成:

(B1)

其中,e1、e2和e3表示極化,s1為X1方向上的水平慢度,s2為X2方向上的水平慢度,s3為垂直慢度;水平慢度s1和s2分別為式(B2)、(B3)所示:

s1=ξsinθcosφ,

(B2)

s2=ξsinθsinφ,

(B3)

Tsvankin(2012)給出了平面波的解析表達:

(B4)

其中,μ為位移;Cijkl為彈性剛度張量.將式(B1)代入式(B4),可以得出:

ρμi=Cijklslsjμk,

(B5)

其中,Cijklslsj定義為Christoffel矩陣(Mavko et al.,2020),該矩陣建立了彈性模量與水平慢度和垂直慢度之間的聯系.

將式(B5)進行等價改寫成:

(ρδik-Cijklslsj)μk=0,

(B6)

其中,若i等于k,則δik為1;若i不等于k,則δik為0.

通過求解式(B7)中的方程,可以得到垂直慢度s3與彈性模量Cijkl之間的關系:

|ρδik-Cijklslsj|=0,

(B7)

其中,特征值分別是準P波、第一個到達的準S波和最后一個到達的準S波的垂直慢度s3p、s3s和s3t的平方解,相關的特征向量為相應的極化ep、es和et.

上部介質中的地震波場包括入射和反射的P波和S波,而下部介質只有透射波(Jin et al.,2018).上部介質的位移場可以寫成:

(B8)

其中,ip,rp,is,rs,it和rt分別是入射準縱波、反射準縱波、入射快準橫波、反射快準橫波、入射慢準橫波和反射慢橫波的幅值;ep1,ep2,ep3;es1,es2,es3;et1,et2,et3分別為準P波、第一個到達準S波和最后一個到達準S波的極化;s3p,s3s,s3t分別為準P波、第一個到達準S波和最后一個到達準S波的垂直慢度;ω為角頻率.

而在下層介質中只有透射的P波和S波,其位移場表示為:

(B9)

其中,tp,ts和tt分別為透射的準縱波、透射的快準橫波和透射的慢準橫波的幅值;所有其他帶撇號的參數具有與式(B8)中不帶撇號的參數相同的含義.

可以看出反射系數計算的關鍵是求出垂直慢度相應的極化值.由Christoffel方程可以從給定的水平慢度條件下,計算出垂直慢度和相應的極化.Zoeppritz方程給出了在兩個彈性介質之間的界面上,入射的縱波和橫波的反射和透射平面波振幅的精確計算.Schoenberg和Protázio(1992)通過引入兩個阻抗矩陣來求解Zoeppritz方程,得到了水平和垂直慢度值,以及相應的極化值,兩個阻抗矩陣的表達式為:

(B10)

(B11)

其中,Cij表示入射介質實數域的剛度矩陣元素.