基于馬爾科夫鏈的裂縫期次分析
——以重慶天府地區上二疊統長興組露頭為例

袁 帥， 胡 明,3*， 李 瑞， 王 豪， 胡 峰

(1.西南石油大學地球科學與技術學院，成都 610500；2.中國地質大學構造與油氣資源教育部重點實驗室，武漢 430074；3.天然氣地質四川省重點實驗室，成都 610500；4.東方地球物理公司西南物探分公司，成都 610213)

裂縫研究對于了解巖體強度、巖體變形性、巖體穩定性、流體流動、油氣儲集與開采、斷裂機制和構造歷史至關重要[1]。目前有多種研究構造裂縫發育期次的方法,其中基于巖心資料的裂縫觀察、裂縫充填物同位素分析、包裹體測溫、巖石聲發射等方法都是研究裂縫期次的有效方法。然而，由于裂縫具有多期成因、多期改造、多期充填的特點，每一種技術方法都具有一定的局限性或多解性[2]，且實驗方法涉及采樣、制樣等繁瑣的實驗步驟，對樣品、操作以及設備的要求極高，費用相當昂貴且容易出現誤差。外國學者以概率學為理論基礎提出了一種基于馬爾科夫鏈的時間順序研究方法，這種方法曾用于沉積學中地層旋回模式、復理石韻律結構等研究[3-6]，Snyder將其公式修改之后首次運用到裂縫期次的研究之中[7-8]。為此，在前人的研究基礎之上，以重慶天府地區上二疊統長興組大量出露的裂縫為研究對象，用馬爾科夫鏈分析該地區的裂縫發育期次,再利用裂縫充填物碳氧同位素測試結果對結論加以驗證，意在為日后油氣勘探或基礎地質研究提供一種經濟有效的裂縫期次研究方法。

1 研究區概況

1.1 地質概況

研究區(圖1)位于重慶天府地區，研究對象為觀音峽背斜核部或靠近核部的長興組地層。觀音峽背斜是川東高陡褶皺帶西南方向的分支。整體呈北東-南西向延伸，出北碚后轉向南延，經中梁山后轉向東南伸達九龍坡，被長江切穿軸部，過江以后逼近綦江而傾伏。背斜為兩翼不對稱的斜歪褶曲，軸面傾斜(先東南傾斜，后西北傾斜)，軸線扭轉彎曲(反“S”形)，樞紐鞍狀起伏[9]。

1.2 裂縫發育特征

裂縫的基本參數包括裂縫的寬度(張開度)、長度、間距、密度、產狀、充填情況以及溶蝕改造等[10]。通過大量的野外觀察和測量，發現本地區天然裂縫極為發育，主要以構造成因的張性裂縫為主[圖2(a)]，裂縫面較為粗糙且不平整，有少量的X形剪切縫[圖2(b)]以及雁列式節理組[圖2(c)]；裂縫面較寬(即張開度較大)，可高達3～4 cm，且充填情況較好，充填深度可達20～25 cm[圖2(d)]；單條裂縫長度最長可達20～30 cm，而最短不過數厘米；裂縫發育較為密集，平均為5條/m；由于靠近背斜轉折端，裂縫所在巖層產狀均較為平緩，傾角為20°～30°，且裂縫多為高角度裂縫或垂直縫；研究區裂縫發育多而雜亂，交切關系復雜，沒有明顯的交切關系能說明裂縫發育期次，但總體而言有四個優勢方位，分別是北東向、近東西向、南東東向以及南南東向(圖3)。

圖1 研究區地理位置圖

圖2 研究區裂縫特征

圖3 研究區裂縫露頭及主要方向

2 馬爾科夫鏈綜合分析

野外觀察發現，裂縫展布極為復雜，無法簡單地根據若干條裂縫之間的分期配套得出整個區域的裂縫發育期次，因此需要一種方法來判斷是否某一方向的裂縫頻繁地與其他方向裂縫顯示出時間的先后關系，從而評估研究區裂縫的發育期次。馬爾科夫鏈就是這樣一種研究時間序列中事件之間關系的統計方法。

當一個隨機過程在給定現在狀態及所有過去狀態情況下，其未來狀態的條件概率分布僅依賴于當前狀態，稱這個隨機過程具有馬爾可夫性質，是一個馬爾科夫過程[11]。地質過程的馬爾科夫鏈模擬作為地質分析的一種手段，其有效性已經被中外地質界所應用。構造裂縫的形成就是一個隨機過程，且下一條裂縫的狀態概率分布僅與當前裂縫有關，與上一條裂縫無關，因此理論上該隨機過程具有馬爾科夫性質，可使用馬爾科夫鏈進行研究。

2.1 基本方法

馬爾科夫鏈是一個基于概率統計的方法，在運用到裂縫研究中時，先通過大量的野外觀察和測量，將各個方向裂縫之間的時間先后關系統計到矩陣(實際轉移頻數或觀測值)，再通過公式計算出一個隨機的時間先后關系(隨機轉移頻數或估計值)，最后計算兩者的標準化殘差。意在通過標準殘差值來對比是否某一方向的裂縫會比預期更頻繁地與另一方向的裂縫產生時間關系，并最終得出一個總體的時間序列。具體方法如下。

(1)數據準備：馬爾科夫鏈分析是一種統計學方法，對樣本數據的準確度有著較高的要求，如何減少樣本數據采集時的誤差是必須要考慮的問題，現采用圓形測線法來收集樣本數據[12]:在研究區地層露頭上用粉筆畫直徑為1 m的圓，然后記錄圓內裂縫的產狀，交切關系等特征(圖4)。由于研究區地層較為平緩(巖層傾角均小于20°)，則可不考慮巖層對裂縫產狀測量的影響。研究共使用測線50組，測得裂縫797條，記錄了346次交切關系。

地層中展布的各個時期形成的裂縫稱之為狀態，某一時刻裂縫從一種狀態轉移到另一種狀態，稱之為狀態轉移。狀態的選擇是應用馬爾科夫鏈分析的關鍵。狀態的選擇要考慮到實際的裂縫交切關系，若狀態區域劃分太大，則容易忽略掉同一個區域內裂縫的轉移。現以玫瑰花圖中顯示的“極小值”為邊界，分別是25°、65°、95°以及145°，共劃分了4個區域(圖5)，結合野外觀察，設立了四種狀態用于馬爾科夫鏈分析，分別是以45°為主的北東向裂縫、以165°為主的南東向裂縫，以80°為主的近東西向裂縫和以120°為主的南東東向裂縫，其中北東向裂縫和近東西向裂縫雖發育數量較少，但與其他裂縫之間存在時間關系，故為不可忽略的狀態。

(2)轉移頻數：根據裂縫的交錯和限制關系，將某一方位的裂縫晚于另一方位裂縫的次數作為轉移頻數，并記錄在表格中(表1)，其縱列表示較晚的裂縫，橫列表示較早的裂縫(如表1中B行A列(記作OB,A)數據表示B區域的裂縫晚于A區域的裂縫的次數為1次)，并將總的轉移頻數記錄在表格的最右下角。

圖4 使用圓形測線的野外數據統計

圖5 節理玫瑰花圖及分區

表1 轉移頻數矩陣

在計算出系統的狀態轉移頻數矩陣的基礎上，檢驗該系統是否具有馬爾科夫性質可選用統計量x2進行假設檢驗。Snyder根據Cochran的理論進行x2檢驗[13]，然而這一檢驗方法并沒有結合馬爾科夫鏈的特征，現在前人研究的基礎上對檢驗方法進行改善[14]，采用公式：

(1)

式(1)中：fij為轉移頻數矩陣的值；pj為轉移頻數矩陣第j列之和除以總轉移頻數的值，稱為邊際概率；pij為fij除以總轉移頻數的值，即轉移概率(表2)。

(3)隨機轉移值計算：根據觀測的實際轉移數據即轉移頻數，可以計算出狀態的隨機轉移數據。通常情況下，其公式為

EA,B=SASB/N

(2)

式(2)中：EA,B為A行B列的隨機轉移值；SA為A區域觀測到的總時間關系數；SB為B區域觀測到的總時間關系數，N為所有測線內時間關系總數。

然而裂縫之間的角度關系影響著裂縫之間的交切次數，如在一個固定區域內，相互垂直的裂縫會比斜交的裂縫出現更多的交切次數(圖6)，因此需要將隨機值乘以一個因子sin(|θA-θB|)以減小角度對交切關系的影響。修改后的公式為

表2 轉移概率矩陣

圖6 角度關系對交切次數的影響[7]

(3)

式(3)中：SA′為A區域總的裂縫數；SB′為B區域總的裂縫數；θA為A區域中點的方位角度數；θB為B區域中點的方位角度數；N′為測線中觀測到的總的裂縫數。

再將所得隨機數據乘以N/SE(SE為隨機數據之和)，使隨機數據之和等于觀測數據之和，以便于下一步的殘差對比，計算結果見表3。

表3 隨機數據矩陣

(4)標準殘差計算：標準化殘差指的是觀察值與隨機值之差的標準化，突出的是表示觀測值偏離隨機值的程度(表4)。計算時利用公式：

(4)

式(4)中：RA,B為標準殘差值；OA,B為A行B列的轉移頻數；EA,B為A行B列的隨機轉移值。

在本研究中，其值表示某一方向裂縫比隨機值更頻繁或者更少晚于另一方向裂縫的程度。根據標準殘差表，比較表格中斜對角正負相反的單元值來確定相對的時間。如果一個單元格值為正，而斜對角的值為負，那么這一對的列模式就發生在行模式之前(如D行C列值為正值而其對角線C行D列為負值，則D行所代表的以165°為主的裂縫晚于C列中以120°為主的裂縫)。在進行這些比較時，若對角單元格的標準殘差值均為負值，則這組值沒有任何意義。

2.2 結果分析

根據標準殘差表中所得數據(表4)，A區域和B區域(A行B列和B行A列)對應的值分別是-1.83和-1.28，均為負值，為無意義數據；A區域與C區域對應的值分別是0.43和-2.13，表示A區域裂縫形成時間晚于C區域裂縫；A區域與D區域對應的值分別是-1.59和-0.82，無意義；B區域與C區域對應的值分別是-0.85和3.19，表示C區裂縫形成時間晚于B區裂縫；B區域與D區域對應的值分別是-2.46和5.08，表示D區域裂縫形成時間晚于B區域裂縫；C區域與D區域對應的值分別是-3.66和3.56，表示D區域裂縫形成時間晚于C區域裂縫。

表4 標準殘差表

根據上述結論得出裂縫形成的時間順序(圖7)：以80°為主的近東西向裂縫為最早生成的，其次是以120°為主的南東東向裂縫，以165°為主的南南東向裂縫和以45°為主的北東向裂縫沒有明顯的先后關系，推測為同期生成，結合野外裂縫觀察，這兩組裂縫可能是一對X形共軛剪切縫或是在形成時期扭動而成。

圖7 裂縫形成的時間順序

3 充填物碳氧同位素特征

裂縫充填物碳、氧同位素分析也是劃分裂縫形成期次常用的方法。當地層水進入裂縫時，或多或少有結晶礦物析出并沉淀在裂縫壁上，結晶方解石時的古地溫、古水介質條件和有機碳的影響決定了碳、氧同位素的豐度，不同時期充填物的碳、氧同位素值必然不同[15-19]。由于研究區裂縫方解石充填程度較好，因此用該方法分析研究區裂縫發育期次，從而對上述結果進行驗證。

裂縫的形成過程屬于“真空擴容”過程，飽和地層水在第一時間進入裂隙中，必然在縫壁富集，受到多期礦化水影響，可能存在多期膠結物。因此，在研究裂縫形成期次時，采集的裂縫充填物必須是靠近縫壁的，且只能對第一期充填的次生礦物進行分析[20]。此次實驗嚴格選取各期次裂縫對應的充填物，采樣時避開被風化的縫面充填物，采深層靠近縫壁的新鮮充填物進行同位素分析，結果如圖8所示。第一期裂縫充填物的δ18O和δ13C集中在Ⅰ區，δ18O平均值為-8.58‰，δ13C平均值為3.32‰；第二期裂縫充填物的δ18O和δ13C集中在Ⅱ區，δ18O平均值為-7.16‰，δ13C平均值為3.79‰；第三期裂縫充填物的δ18O和δ13C集中在Ⅲ區，δ18O平均值為-5.18‰，δ13C平均值為3.70‰。研究區裂縫充填物δ18O和δ13C的集中分布在3個區域內，因此，判斷裂縫有3個形成期次。

圖8 裂縫充填方解石碳氧同位素分布圖

為進一步分析充填物的形成時期，根據前人的研究方法計算了其鹽度指數Z、形成溫度以及埋深(表5)。根據Weber(1964)公式計算其Z用于區分成巖環境[21-22]，計算公式(Peedee Belemnite,即PDB標準)為

Z=2.048(δ13C+50)+0.498(δ18O+50)

(5)

Z在120以上的碳酸鹽巖應歸入海水成因環境，Z在120以下的應歸入淡水成因環境。碳氧同位素中的氧同位素可以指示充填物形成時的古溫度，現采用Epstein等[23]提出的氧同位素測溫方程來計算充填物充填時的溫度，其方程式為

T=31.9-5.55(δ18O-δ13Ow)+0.7(δ18O-δ13Ow)2

(6)

式(6)中：T為方解石礦物形成時的溫度， ℃；δ18O為測定礦物的氧同位素值，‰；δ13Ow為形成方解石礦物時水介質的氧同位素值，‰。

計算時取海水介質δ13Ow為-2‰，可以得出其充填溫度，再根據研究區古地表溫度(15°)和古地溫梯度(3.0 ℃/hm)[24]，可以折算出裂縫形成時的埋藏深度。根據換算，第一期裂縫形成的平均埋深為3 864 m，形成溫度為130 ℃;第二期裂縫形成的平均埋深為3 084 m，形成溫度為107 ℃;第三期裂縫形成的平均埋深為2 147 m，形成溫度為79 ℃。

4 構造演化史及埋藏史分析

結合研究區構造演化史，天府地區上二疊統長興組地層主要經歷了燕山期和喜山期兩期構造應力場的作用[25]。燕山時期，受太平洋板塊向亞洲板塊俯沖的影響，華南板塊內部發生了大規模的由南東向北西的推覆作用，產生南東向北西的強烈擠壓，從雪峰至江南古陸地區經湘西北至川東，沖斷作用的強度逐漸降低，應力傳遞到天府地區，其區域應力場為北西向，形成了南東東向和近東西向構造裂縫。喜山時期，受太平洋板塊俯沖及印度板塊向歐亞板塊碰撞擠入雙重影響，且有扭動活動，產生了北西-南東向構造應力，在此應力的作用下，形成了南南東向和北東向裂縫。

再根據研究區所在地區埋藏史(圖9)，第一期裂縫形成時平均埋深為3 864 m，為燕山早期產物；第二期裂縫形成時平均埋深為3 084 m，為燕山晚期產物；第三期裂縫形成的平均埋深為2 147 m，為喜山早期產物。

圖9 渝西地區埋藏史圖[26]

5 結論

(1)根據馬爾科夫鏈分析結果，以80°為主的近東西向裂縫為最早生成的，其次是以120°為主的南東東向裂縫，以165°為主的南南東向裂縫和以45°為主的北東向裂縫沒有明顯的先后關系，這兩組裂縫可能是一對X形共軛剪切縫或是在形成時期扭動而成。

(2)充填物碳氧同位素實驗顯示，第一期裂縫充填物δ18O平均值為-8.58‰，δ13C平均值為3.32‰，形成時平均埋深為3 864 m，形成溫度為130 ℃；第二期裂縫充填物δ18O平均值為-7.16‰，δ13C平均值為3.79‰，形成時平均埋深為3 084 m，形成溫度為107 ℃；第三期裂縫充填物δ18O平均值為-5.18‰，δ13C平均值為3.70‰，形成時平均埋深為2 147 m，形成溫度為79 ℃。結合研究區構造演化史以及埋藏史，第一期裂縫形成于燕山早期；第二期裂縫形成于燕山晚期；第三期裂縫形成于喜山早期。

(3)充填物碳氧同位素的實驗結果成功驗證了馬爾科夫鏈的分析，證實馬爾科夫鏈為一種可行的裂縫期次研究方法，可用于日后油氣勘探或基礎地質研究之中。