構(gòu)造物理模擬和PIV有限應(yīng)變分析對構(gòu)造裂縫預(yù)測的啟示

沈禮，賈東，尹宏偉，魏東濤，陳竹新，孫闖，崔鍵.南京大學(xué)能源科學(xué)研究院，地球科學(xué)與工程學(xué)院，南京2002 2.中國石油天然氣集團(tuán)公司油藏描述重點實驗室，蘭州70020 .中國石油勘探開發(fā)研究院，北京0008 4.中國石油勘探開發(fā)研究院西北分院，蘭州70020

?

構(gòu)造物理模擬和PIV有限應(yīng)變分析對構(gòu)造裂縫預(yù)測的啟示

沈禮1，2，賈東1*，尹宏偉1，魏東濤2，4，陳竹新3，孫闖1，崔鍵1
1.南京大學(xué)能源科學(xué)研究院，地球科學(xué)與工程學(xué)院，南京210023 2.中國石油天然氣集團(tuán)公司油藏描述重點實驗室，蘭州730020 3.中國石油勘探開發(fā)研究院，北京100083 4.中國石油勘探開發(fā)研究院西北分院，蘭州730020

摘要：油氣儲層中構(gòu)造裂縫發(fā)育與有限應(yīng)變狀態(tài)關(guān)系密切，為了探索有限應(yīng)變分析與構(gòu)造裂縫預(yù)測的新技術(shù)方法，此次研究設(shè)計完成了一組單側(cè)擠壓收斂模型的物理模擬實驗，并引入粒子圖像測速（PIV，Particle Image Velocimetry）技術(shù)對實驗過程進(jìn)行了定量化分析。實驗?zāi)Ｐ驮诖瓜蛏蠟楹承詫拥亩鄬咏Y(jié)構(gòu)，實驗結(jié)果形成了一個肉眼可見的箱狀褶皺。通過PIV技術(shù)可以獲取實驗?zāi)Ｐ妥冃窝莼^程中各階段的位移場數(shù)據(jù)，計算出各階段的增量應(yīng)變，實現(xiàn)從初始狀態(tài)到褶皺形成之后整個變形過程的有限應(yīng)變分析，探討構(gòu)造裂縫成因機(jī)制和分布規(guī)律，進(jìn)行定量化裂縫預(yù)測。擠壓變形過程初期，應(yīng)變分布范圍很廣，有限應(yīng)變較弱（約4%～8 %），在擠壓方向上的線應(yīng)變表現(xiàn)為弱壓應(yīng)變，在垂向上的線應(yīng)變表現(xiàn)為弱張應(yīng)變，這種現(xiàn)象是褶皺和斷層產(chǎn)生前平行層縮短和層增厚的純剪變形結(jié)果，也是區(qū)域型張裂縫和剪裂縫形成的主要機(jī)制。褶皺和斷層即將發(fā)育之時至發(fā)育之后，應(yīng)變局限在斷層發(fā)育的剪切帶及附近區(qū)域，有限應(yīng)變表現(xiàn)為較強(qiáng)（達(dá)20%）的剪切應(yīng)變和剪切張應(yīng)變，是斷層面附近簡單剪切變形作用的結(jié)果，也是局部型剪裂縫和張剪裂縫形成的主要機(jī)制。

關(guān)鍵詞：物理模擬；粒子圖像測速（PIV）；有限應(yīng)變分析；構(gòu)造裂縫預(yù)測

1　引言

天然裂縫是儲層中油氣運移的重要通道，按成因可分為構(gòu)造裂縫和非構(gòu)造裂縫。構(gòu)造裂縫的發(fā)育狀況與構(gòu)造事件相關(guān)聯(lián)，包括與褶皺、斷層等局部構(gòu)造事件相關(guān)的裂縫和在弱構(gòu)造變形區(qū)形成的裂縫（Nelson，1985）；非構(gòu)造裂縫的發(fā)育狀況與構(gòu)造作用無關(guān)，包括成巖裂縫、收縮裂縫以及與表面相關(guān)的裂縫（Lorenz，1991）。在致密儲層中，裂縫控制著油氣藏的分布（周新桂等，2003）。由于致密儲層中非構(gòu)造裂縫的孔徑小、滲透率低以及橫向連通性差等因素，構(gòu)造裂縫成為了儲層中主要的滲流通道，控制著儲層中的油氣成藏過程（王鵬威等，2014；曾聯(lián)波等，2008；Zeng and Li，2009）。構(gòu)造裂縫定量預(yù)測可以優(yōu)化井位設(shè)置，對于油氣藏開采具有指導(dǎo)意義。

構(gòu)造裂縫定量預(yù)測的主要方法有：（1）構(gòu)造主曲率法，該方法假設(shè)與撓曲有關(guān)的裂縫發(fā)生在構(gòu)造面主曲率或傾角變化率極大處（曾錦光等，1982；Murray，1968）。（2）應(yīng)力場數(shù)值模擬與巖石破裂法，利用有限元數(shù)值模擬確定地質(zhì)模型中的應(yīng)力分布，再根據(jù)巖石脆性破裂準(zhǔn)則預(yù)測脆性破裂發(fā)生的位置（陳波和田崇魯，1998；沈國華，2008；宋惠珍等，1999；文世鵬和李德同，1996）。（3）應(yīng)變能法，通過計算地層中單位體積應(yīng)變能來預(yù)測裂縫發(fā)育程度，該方法認(rèn)為同樣厚度的巖石中應(yīng)變能越高，裂縫發(fā)育越多（郭迎春等，2010；Price，1966）。（4）統(tǒng)計學(xué)方法，通過統(tǒng)計學(xué)規(guī)律推測地下不同厚度巖層中裂縫間距的大小，該方法也適用于非構(gòu)造裂縫（Narr，1991； Narr and Lerche，1984）。構(gòu)造裂縫的成因和影響因素很多，每種方法的依據(jù)可能只是裂縫的一種成因機(jī)制或分布規(guī)律，所以研究者們也會綜合利用多種方法來進(jìn)行構(gòu)造裂縫預(yù)測（丁中一等，1998；高帥等，2015；周新桂等，2007），預(yù)測效果較之單一方法有所改善，但仍然存在一些局限，還不能夠全面地解決構(gòu)造裂縫定量預(yù)測問題。

本文以物理模擬實驗和粒子圖像測速（PIV）技術(shù)為基礎(chǔ)，從構(gòu)造裂縫發(fā)育和有限應(yīng)變狀態(tài)的關(guān)系出發(fā)，探索構(gòu)造裂縫定量預(yù)測的新技術(shù)方法。物理模擬實驗可以提供褶皺、斷層等構(gòu)造發(fā)育形成過程的真實可見物理模型（Koyi，1997）。運用PIV技術(shù)可以實現(xiàn)對模型中構(gòu)造發(fā)育演化過程的定量監(jiān)測（沈禮等，2012；Adam et al.，2005），分析構(gòu)造變形發(fā)育之前和形成之后的有限應(yīng)變狀態(tài)，再結(jié)合構(gòu)造裂縫發(fā)育與有限應(yīng)變狀態(tài)的關(guān)系進(jìn)行定量化構(gòu)造裂縫預(yù)測。

2　構(gòu)造物理模擬實驗

2.1實驗?zāi)Ｐ?/p>

本文旨在探索有限應(yīng)變分析和構(gòu)造裂縫預(yù)測的新技術(shù)方法，因此在進(jìn)行構(gòu)造物理模擬時設(shè)計了較為簡單的實驗?zāi)Ｐ汀嶒災(zāi)Ｐ偷某跏汲叽鐬殚L600 mm（x方向）、寬400 mm（y方向）、厚45 mm （z方向）。為了考察粘性層上下變形的差異，模型垂向上設(shè)計為含粘性層的多層結(jié)構(gòu)。模型底部設(shè)置了5 mm厚的微玻璃珠以減少基底摩擦力，使構(gòu)造變形更容易沿基底滑脫面向前傳遞。底部微玻璃珠之上從底面往頂面依次設(shè)置了17 mm厚的石英砂、5 mm厚的硅膠（靠近擠壓端一段為微玻璃珠）和18 mm厚的石英砂（圖1）。

圖1　實驗?zāi)Ｐ推拭媸疽鈭DFig.1 Schematic diagram of the experimental model

實驗?zāi)Ｐ椭胁吭O(shè)置的硅膠是一種粘度系數(shù)為1×104Pa·s的粘性材料，變形規(guī)律符合牛頓流體變形規(guī)則，是模擬自然界巖鹽層、膏巖層等軟弱滑脫層的理想材料（Bonini，2007；Costa and Vendeville，2002）。脆性材料石英砂和微玻璃珠的變形遵循莫爾-庫侖破裂準(zhǔn)則，具有較小的內(nèi)聚強(qiáng)度，石英砂內(nèi)摩擦角約30°，微玻璃珠內(nèi)摩擦角約25°，都是模擬地殼脆性變形的理想材料（Cotton and Koyi，2000；Massoli et al.，2006）。為了提高PIV計算識別精度，實驗中使用了顆粒感更加明顯的雜色石英砂，由白色石英砂和彩色石英砂均勻混合而成。實驗中還采用彩色石英砂設(shè)置了標(biāo)志層以便于觀察模型的變形情況，硅膠層之下用紅色石英砂作標(biāo)志層，硅膠層之上用綠色石英砂作標(biāo)志層。不同顏色石英砂的力學(xué)性質(zhì)是一致的。實驗材料的具體物理性質(zhì)參數(shù)見表1。

實驗動力來源為x方向上的單側(cè)擠壓，擠壓總縮短量為18 mm，縮短速率為0.0015 mm/s。實驗過程中，側(cè)面相機(jī)每隔5分鐘拍攝一幅剖面彩色圖像，PIV系統(tǒng)的高精度CCD鏡頭每隔兩分鐘采集一幅剖面灰度圖像。

2.2模型相似性

表1中列出了實驗?zāi)Ｐ秃妥匀唤缭椭械母魑锢砹繀?shù)以及通過這些參數(shù)計算出來的相似比例參數(shù)。實驗?zāi)Ｐ团c自然界原型的幾何學(xué)相似性由實驗?zāi)Ｐ团c自然界原型的長度相似比例給出：

表1　材料物理性質(zhì)及模型相似性參數(shù)Tab.1 Parameters of the materials and the scaling of the model

l*=2.5×10-6（即實驗中1 cm代表自然界4 km）

動力學(xué)相似性由應(yīng)力相似比例給出，根據(jù)Hubbert(1937)最早提出的相似理論，應(yīng)表示為：

其中，ρ*表示密度相似比例，g*表示重力加速度相似比例。實驗中所用材料與自然界巖石的密度相似比例ρ*≈0.54，實驗在正常重力加速度場中進(jìn)行，所以重力加速度相似比例g*=1。因此，模型與自然界原型的應(yīng)力相似比例σ*= 1.4×10-6。設(shè)自然界滑脫層粘度系數(shù)為1×1019Pa·s（Bonini，2007；Cotton and Koyi，2000；Couzens-Schultz et al.，2003；Leturmy et al.，2000），那么模型與自然界原型的粘度系數(shù)相似比例η*=1×10-15，因此，模型與原型之間的應(yīng)變速率相似比例：

時間相似比例是應(yīng)變速率相似比例的倒數(shù)，所以時間相似比例t*=7.1×10-10，即實驗中的1小時代表了自然界的0.16 Ma，本文實驗擠壓過程耗時200分鐘，代表了自然界中的0.53 Ma。運動學(xué)相似性由實驗?zāi)Ｐ团c自然界原型的縮短速率相似比例給出：

因此，實驗中推擠速率0.0015 mm/s就代表了自然界縮短速率為14 mm/a。

2.3實驗結(jié)果及討論

實驗?zāi)Ｐ妥冃窝莼^程如圖2所示，擠壓縮短量s從0 mm增加到5.4 mm的過程中，模型物質(zhì)沿基底滑脫面向前傳遞，但剖面上并沒有明顯的變形（圖2a-c）。當(dāng)擠壓縮短量達(dá)到8.1 mm時，模型物質(zhì)發(fā)生了微弱的形變（圖2d）。當(dāng)擠壓縮短量達(dá)到10.8mm時，靠近擠壓端的位置形成了一個箱狀褶皺，形態(tài)由一對共軛的剪切帶所控制（圖2e）。擠壓縮短量從10.8 mm增加到18.0 mm的過程中，箱狀褶皺沿共軛剪切帶繼續(xù)發(fā)育、抬升，剪切帶寬度明顯減?。▓D2f）。整個實驗過程中，粘性層上下未見明顯變形。

實驗?zāi)Ｐ蛿D壓縮短過程中形成的箱狀褶皺亦可以描述為斷層端點褶皺，其演化過程從初始階段的平行層縮短開始，接著產(chǎn)生沿基底滑脫面滑動的滑脫褶皺，最后形成斷層端點褶皺（Bernard et al.，2007；Storti et al.，1997）。這與本文實驗中觀察到的現(xiàn)象是一致的。如果縮短量繼續(xù)增加，斷面所在的剪切帶將逐漸變窄，模型最后將演化為斷層斷坡褶皺（Storti et al.，1997）。對物理模擬實驗中這類現(xiàn)象的進(jìn)一步分析可以讓構(gòu)造地質(zhì)學(xué)家更深入地理解斷層相關(guān)褶皺運動學(xué)演化過程，幫助復(fù)雜地震數(shù)據(jù)的解釋和構(gòu)造演化的平衡恢復(fù)（（Daeron et al.，2007；Simoes et al.，2007））。本文接下來將運用PIV技術(shù)分析實驗?zāi)Ｐ椭旭薨櫋鄬樱羟袔В┊a(chǎn)生前、后的有限應(yīng)變狀態(tài)，探討油氣儲層中構(gòu)造裂縫的形成機(jī)制和分布規(guī)律。

3　粒子圖像測速（PIV）技術(shù)

3.1 PIV工作的基本原理

粒子圖像測速，英文為ParticleImageVelocimetry，通常簡稱為PIV。該測速技術(shù)最初是在流體力學(xué)實驗的研究領(lǐng)域發(fā)展起來的（Adrian，1991），后被引入到巖土工程實驗中計算顆粒材料的位移和變形情況（White et al.，2003；Wolf et al.，2003），最近十年才被應(yīng)用到構(gòu)造物理模擬實驗中定量化監(jiān)測實驗?zāi)Ｐ偷臉?gòu)造演化過程（Adam et al.，2005）。與物理模擬實驗傳統(tǒng)分析方法相比，PIV技術(shù)既節(jié)省時間又能提高位移變形監(jiān)測的時間、空間分辨率和精準(zhǔn)度（至少一個數(shù)量級以上），所以近年來被廣泛地應(yīng)用到物理模擬實驗中分析諸如應(yīng)變集中過程（Boutelier et al.，2014；Dooley et al.，2015；Leever et al.，2011；Santimano et al.，2015）、地表過程（Cruz et al.，2008；Reiter et al.，2011；Wu et al.，2015）、地震活動（Knap?pett et al.，2006；Rosenau et al.，2009）以及斷層的粘滑運動（Yamadaet al.，2006）等地質(zhì)過程。值得一提的是，與PIV技術(shù)工作原理類似的數(shù)字圖像相關(guān)（DIC,Digital Image Correlation）技術(shù)與光學(xué)流動技術(shù)（Optical Flow Technique）也經(jīng)常被應(yīng)用到物理模擬實驗中監(jiān)測模型的運動和變形過程（Bernard et al.，2007；Bonini et al.，2015；Grave?leau et al.，2008；Hall et al.，2010）。

PIV技術(shù)是一種可以實現(xiàn)非線性位移變形可視化的光學(xué)、非接觸式圖像相關(guān)技術(shù)。其工作的基本原理如圖3，通過PIV系統(tǒng)高分辨率鏡頭獲取實驗?zāi)Ｐ妥冃芜^程中的一系列灰度圖像（圖3A），每幅圖像都可以用對應(yīng)的灰度坐標(biāo)函數(shù)表示，不同圖像函數(shù)可以通過互相關(guān)函數(shù)聯(lián)系起來。將時間間隔為dt的兩幅圖像劃分為多個問詢域（圖3B），通過基于快速傅里葉變化的互相關(guān)運算可以計算出局部問詢域內(nèi)顆粒在dt時間內(nèi)的位移，而互相關(guān)函數(shù)的峰值代表了局部問詢域內(nèi)顆粒的平均位移（圖3C）。這樣每個問詢域都獲得了一個位移矢量，所有問詢域的位移矢量就構(gòu)成了整幅圖像在dt時間內(nèi)的位移矢量場。對于位移和變形比較大的地方，通過逐步縮小問詢域的迭代算法和變形窗口算法可以提高位移矢量的空間分辨率（Fincham and Spedding，1997；Scarano and Riethmuller，2000；Wieneke，2001）。以增量位移矢量為基礎(chǔ)，便可以進(jìn)一步獲取模型的增量應(yīng)變，如線應(yīng)變、剪切應(yīng)變和面積應(yīng)變等。將每個階段的增量應(yīng)變按順序求和，也就獲得了模型變形演化過程中的有限應(yīng)變狀態(tài)。

圖3　PIV技術(shù)工作基本原理Fig.3 Fundamental of PIV technique (modified after Adam et al.,2005)

3.2 PIV計算結(jié)果的精度

PIV計算位移矢量數(shù)據(jù)的空間分辨率由PIV系統(tǒng)數(shù)字鏡頭的光學(xué)分辨率和互相關(guān)算法的精度決定（Adam et al.，2005）。普通的四百萬像素數(shù)碼相機(jī)照片就能夠滿足PIV分析的需求，但是其采集的圖像通常為彩色圖像，而在進(jìn)行PIV互相關(guān)運算前必須轉(zhuǎn)換為灰度圖像，這個轉(zhuǎn)換和壓縮過程將會損失部分圖像信息，降低PIV分析的精度。PIV系統(tǒng)的鏡頭可以直接采集灰度圖像，比較完整的保留圖像的灰度梯度。根據(jù)前人的研究，PIV系統(tǒng)計算的位移矢量標(biāo)準(zhǔn)差小于0.1個像素（White et al.，2003）。本文采用的PIV分析系統(tǒng)為北京立方天地科技發(fā)展有限責(zé)任公司所開發(fā)，其技術(shù)報告中對人工旋轉(zhuǎn)場粒子圖像的計算結(jié)果進(jìn)行了誤差分析，用實例證明了其計算結(jié)果的標(biāo)準(zhǔn)差小于0.1像素。本文PIV系統(tǒng)的CCD（Charge-coupled Device）鏡頭光學(xué)分辨率為四百萬像素（2048×2048），用于PIV計算分析的實驗?zāi)Ｐ烷L度為42 cm，所以本文計算結(jié)果位移矢量長度的空間分辨率為0.02 mm。應(yīng)變測量的相對誤差與問詢域的尺寸和計算區(qū)域重疊率相關(guān)，本文在進(jìn)行互相關(guān)計算時采用的問詢域大小為16× 16像素，步長為12×12像素，重疊率為75%，計算結(jié)果相對誤差小于1%。時間分辨率受限于數(shù)字鏡頭采集圖像的頻率（幀轉(zhuǎn)移的速度），普通的數(shù)碼相機(jī)采集圖像時需要對圖像進(jìn)行壓縮，采集頻率大約為0.1 Hz，而本文PIV系統(tǒng)的CCD數(shù)字鏡頭采集圖像時不需要對圖像進(jìn)行壓縮，采集頻率最大可達(dá)8 Hz，并且還不會丟失圖像信息。

圖4　褶皺發(fā)育前的增量位移和增量應(yīng)變狀態(tài)（底圖對應(yīng)圖2　c）Fig.4 Incremental displacement field and incremental strain state before the fold initiated (Background corresponding to Fig.2c)

3.3褶皺發(fā)育前后的增量應(yīng)變分析

增量位移量和總的位移量數(shù)據(jù)可以提供實驗?zāi)Ｐ驼麄€構(gòu)造變形過程的應(yīng)變演化歷史。圖4a展示了實驗?zāi)Ｐ婉薨櫚l(fā)育前（實驗時間t從60分鐘到80分鐘，dt=20分鐘；縮短量s從5.4 mm到7.2 mm，ds=1.8 mm）的增量位移矢量場，底圖為實驗進(jìn)行到60分鐘，縮短量為5.4 mm時的圖像（圖2c）。圖4b是根據(jù)位移矢量長度所成的增量位移彩色云圖，右側(cè)的色標(biāo)給出了每種顏色代表的具體位移量。在靠近擠壓端的位置，增量位移矢量場顯示出一個斜邊并不十分明顯的紅色直角三角區(qū)域，其位移方向與實驗?zāi)Ｐ蛿D壓縮短方向一致，且位移量大小約為1.8 mm，與模型擠壓端的縮短量相當(dāng)。在這個斜邊的右側(cè)，位移矢量方向幾乎均為右上方向，僅在底面附近和粘性層右端表現(xiàn)為純向右，這在水平位移量（圖4c）和垂向位移量（圖4b）的分布范圍上也可以表現(xiàn)出來。緊挨著紅色直角三角區(qū)，呈現(xiàn)出一個倒立的黃色等腰三角區(qū)，位移量約1.3 mm，再往右是一個近似平行四邊形的綠色區(qū)域，位移量約0.8～1.0 mm。繼續(xù)向右來到模型前端含粘性層的地方，可以發(fā)現(xiàn)位移量在粘性層之上比之下覆蓋范圍更大，但位移量都小于0.6 mm。位移場在模型剖面上表現(xiàn)出了從擠壓端向前逐漸減小的梯度。

增量應(yīng)變可以根據(jù)增量位移場計算出來，本文主要分析了與構(gòu)造裂縫預(yù)測相關(guān)的水平線應(yīng)變（exx）、垂向線應(yīng)變（ezz）、面積應(yīng)變（exx+ezz）和剪切應(yīng)變分量（exz）。褶皺發(fā)育之前，水平線應(yīng)變、垂向線應(yīng)變和面積應(yīng)變都比較分散，分布范圍都很廣，沒有表現(xiàn)出應(yīng)變集中的位置，也都較弱，大小約為3%～4 %（圖4e-g）。水平線應(yīng)變?yōu)槿鯄簯?yīng)變，垂向線應(yīng)變?yōu)槿鯊垜?yīng)變，面積應(yīng)變既有張應(yīng)變也有壓應(yīng)變，規(guī)律性不強(qiáng)。剪切應(yīng)變則主要集中在底面和粘性層面之上，強(qiáng)度較大，均大于7.5 %（圖4h）。

褶皺發(fā)育之后（實驗時間t從180分鐘到200分鐘，dt=20分鐘；縮短量s從16.2 mm到18.0 mm，ds=1.8 mm），靠近擠壓端的直角三角區(qū)域增量位移大小和方向都仍與擠壓端一致，斜邊邊界更加清晰，對應(yīng)反沖斷層（BT）所在的位置，右側(cè)倒立等腰三角區(qū)域位移方向為右上方向，位移大小在1.5 mm左右，前沖斷層（FT）所在的剪切帶位移明顯減小約為0.8～1.0 mm，前沖斷層下盤位移全部為零（圖5a-d）。褶皺發(fā)育之后增量位移梯度非常明顯，增量應(yīng)變也幾乎全部集中在前沖斷層和反沖斷層所在的剪切帶位置，前沖斷層所在剪切帶應(yīng)變強(qiáng)度明顯強(qiáng)于反沖斷層所在剪切帶（圖5e-h）。

增量應(yīng)變分析顯示，緊挨著模型擠壓端的直角三角區(qū)域在褶皺發(fā)育前后都與模型擠壓端的位移保持一致，也幾乎沒有內(nèi)部應(yīng)變，它的作用可以等同于不變形的后擋板。褶皺發(fā)育前后最明顯的區(qū)別在于褶皺發(fā)育前應(yīng)變比較分散，分布范圍大，強(qiáng)度較弱；而褶皺發(fā)育后應(yīng)變主要集中在反沖斷層和前沖斷層附近的剪切帶，強(qiáng)度較大。從應(yīng)變產(chǎn)生的機(jī)制上，褶皺發(fā)育前位移梯度主要靠平行層縮短和層增厚吸收，應(yīng)變屬于純剪切變形的結(jié)果；而褶皺發(fā)育后位移梯度明顯反應(yīng)出兩條剪切帶的活動，應(yīng)變屬于簡單剪切變形的結(jié)果。另外，模型底面和粘性層面與其上覆層之間的剪切應(yīng)變也應(yīng)該屬于簡單剪切變形的結(jié)果。

4　有限應(yīng)變分析

實驗?zāi)Ｐ妥冃稳^程的有限應(yīng)變是所有間斷過程增量應(yīng)變的總和。本文實驗過程中，PIV系統(tǒng)每隔2分鐘采集一次圖像，便可以計算一次增量位移和增量應(yīng)變。由于實驗?zāi)Ｐ蛷鸟薨櫚l(fā)育之時至實驗結(jié)束的過程中，變形方式幾乎沒有發(fā)生變化，所以本文僅選取了模型初始狀態(tài)到褶皺發(fā)育之后（實驗進(jìn)行時間為120分鐘時）的過程進(jìn)行了有限應(yīng)變分析。圖6-8展示了這120分鐘內(nèi)實驗?zāi)Ｐ推拭孀冃慰偽灰屏亢陀邢迲?yīng)變的累積過程，即60個連續(xù)間斷過程（dt=2分鐘）增量位移和增量應(yīng)變的累加過程。

實驗?zāi)Ｐ推拭孀冃挝灰屏坷鄯e過程顯示，實驗從初始狀態(tài)到擠壓縮短量為8.1mm的過程中，累積位移量在剖面上大部分地方均逐漸增加（圖6a-d），而在擠壓縮短量從8.1 mm增加至10.8 mm的過程中，累積位移量只在前沖斷層左側(cè)區(qū)域增加，這說明在前沖斷層即將發(fā)育之時至發(fā)育之后，其下盤的位移量就不再增加。在靠近擠壓端的直角三角區(qū)域內(nèi)累積位移量始終與擠壓縮短量保持大致相等（圖6e），說明該區(qū)域內(nèi)的物質(zhì)在整個過程中都與擠壓端的運動基本一致。

有限剪切應(yīng)變的累積過程顯示，擠壓縮短量達(dá)到8.1 mm之前，有限剪切應(yīng)變的分布比較廣泛，應(yīng)變強(qiáng)度隨著擠壓縮短過程逐漸增大（圖7a-d）。但除了在粘性層所在位置較為集中，強(qiáng)度較大（在擠壓縮短量為8.1 mm時約20%）外，其余地方都較弱（在擠壓縮短量為8.1 mm時約4%）。在擠壓縮短量從8.1 mm增加至10.8 mm的過程中，有限剪切應(yīng)變在前沖斷層所在位置增加最為明顯，強(qiáng)度達(dá)到20%，前沖斷層左側(cè)區(qū)域有限剪切應(yīng)變只有少量增加，而其下盤有限剪切應(yīng)變完全沒有增加（圖7e）。有限面積應(yīng)變累積過程顯示，在擠壓縮短量達(dá)到8.1 mm之前，剖面上有限面積應(yīng)變在比較寬廣的區(qū)域內(nèi)逐漸增大，且有向模型前端（右側(cè)）擴(kuò)展的趨勢，最遠(yuǎn)可以到達(dá)離擠壓端約30 cm處，但有限面積應(yīng)變的強(qiáng)度較弱，在擠壓縮短量為8.1 mm時約4%～8 %（圖8a-d）。另外，分散的有限面積應(yīng)變還出現(xiàn)了一定程度的分層現(xiàn)象。擠壓縮短量從8.1 mm增加到10.8 mm的過程中，有限面積應(yīng)變在前沖斷層左側(cè)有少量增加，在反沖斷層和前沖斷層所在位置及鄰近區(qū)域增加比較明顯。在縮短量為10.8 mm時，有限面積應(yīng)變強(qiáng)度在反沖斷層及附近位置約為12 %，在前沖斷層及附近位置達(dá)±20 %（圖8e）。

結(jié)合本文前一部分對褶皺發(fā)育前后的增量應(yīng)變分析，我們可以認(rèn)識到，從初始狀態(tài)到褶皺發(fā)育形成之后的整個過程中，分布廣泛的有限應(yīng)變主要來源于褶皺發(fā)育之前，強(qiáng)度較弱；分布較為集中的有限應(yīng)變主要來源于褶皺和斷層即將發(fā)育之時至發(fā)育之后，主要集中在斷層及鄰近區(qū)域，強(qiáng)度較大。賈東等（2007）在采用磁組構(gòu)方法研究斷層相關(guān)褶皺有限應(yīng)變狀態(tài)的工作中也得到了相同的認(rèn)識。說明褶皺和斷層的發(fā)育過程中地層經(jīng)歷了從分散應(yīng)變到應(yīng)變集中這兩種變形方式（Adam et al.，2005；Bernard et al.，2007）。分散應(yīng)變主要反映了純剪切變形造成的平行層縮短和層增厚過程，應(yīng)變集中主要反映的是斷層及附近位置的簡單剪切變形過程。另外，粘性層所在位置集中的有限剪切應(yīng)變反映了粘性層之上的變形比粘性層之下的變形傳遞得更遠(yuǎn)，是粘性層上下簡單剪切變形的結(jié)果。

圖6　總位移量的累積過程Fig.6 Cumulative process of the total displacement

5　構(gòu)造裂縫預(yù)測與展望

天然裂縫按成因機(jī)制可分為構(gòu)造裂縫和非構(gòu)造裂縫兩大類（Lorenz，1991；Nelson，1985）。按裂縫性質(zhì)可以分為與應(yīng)變能相關(guān)的剪裂縫、與巖層褶皺相關(guān)的張裂縫以及與斷層相關(guān)的剪裂縫和張剪裂縫（McQuillan，1973；Murray，1968；Price，1966）。秦啟榮和蘇培東（2006）以構(gòu)造裂縫預(yù)測為目的，依據(jù)構(gòu)造裂縫的形成方式、空間分布以及與構(gòu)造運動的時間關(guān)系將其進(jìn)一步劃分為區(qū)域型、局部型和復(fù)合型構(gòu)造裂縫3大類。

本文通過物理模擬實驗?zāi)M了褶皺和斷層構(gòu)造發(fā)育的整個過程，并運用PIV技術(shù)實現(xiàn)了整個過程的有限應(yīng)變分析，分析結(jié)果可以用來探討構(gòu)造裂縫的成因機(jī)制和分布規(guī)律。在擠壓變形初始階段褶皺和斷層尚未開始發(fā)育之前，地層中廣泛地分布著純剪切變形作用引起的弱應(yīng)變，其分布范圍在擠壓方向上可達(dá)120 km（模型中30 cm），是構(gòu)造事件初期分布廣泛的區(qū)域型張裂縫和剪裂縫形成的主要機(jī)制。在褶皺和斷層即將發(fā)育之時至發(fā)育之后，有限應(yīng)變主要集中在斷層所在位置及鄰近區(qū)域，應(yīng)變強(qiáng)度較大，是斷層面附近簡單剪切變形作用的結(jié)果，也是局域型剪裂縫和張剪裂縫發(fā)育的主要機(jī)制。另外，褶皺和斷層發(fā)育之前產(chǎn)生的分布式有限面積應(yīng)變出現(xiàn)的分層現(xiàn)象，很有可能反映了區(qū)域型水平張裂縫的成層發(fā)育狀況，但也有可能與PIV系統(tǒng)對本文實驗?zāi)Ｐ椭胁噬⑸皹?biāo)志層的識別有關(guān)。

圖8　有限面積應(yīng)變的累積過程Fig.8 Cumulative process of the finite area strain

裂縫發(fā)育程度與有限應(yīng)變強(qiáng)度有關(guān)，構(gòu)造事件初期地層中廣泛分布的弱應(yīng)變約為4%～8%，這些應(yīng)變?nèi)羧糠从碁閺埩芽p，那么儲層中的油氣儲集空間將大大增加，尤其對于致密儲層，這樣的變化甚至可以改變儲層的性質(zhì)，其很有可能演化為比較好的儲層。分布在斷層面附近的強(qiáng)應(yīng)變達(dá)20 %，這些應(yīng)變反映的剪裂縫和張剪裂縫在斷層活動時期可以作為油氣運移的通道，而在斷層不活動的時期則可能起到封堵作用。

以物理模擬實驗和PIV技術(shù)為基礎(chǔ)，定量分析褶皺和斷層構(gòu)造發(fā)育過程中有限應(yīng)變演化歷史，通過有限應(yīng)變狀態(tài)與構(gòu)造裂縫發(fā)育的關(guān)系，探討裂縫成因機(jī)制和分布規(guī)律，是實現(xiàn)構(gòu)造裂縫定量預(yù)測的一種切實可行的辦法。本文定量地分析了構(gòu)造發(fā)育前后的有限應(yīng)變狀態(tài)，探討了構(gòu)造裂縫的成因機(jī)制和分布規(guī)律，但還沒有具體地給出有限應(yīng)變狀態(tài)與構(gòu)造裂縫發(fā)育程度之間的定量關(guān)系，所以還有很多方面值得改進(jìn)。比如，將本文研究結(jié)果與實際研究區(qū)域的地面裂縫、鉆井巖芯裂縫發(fā)育情況的統(tǒng)計結(jié)果結(jié)合起來可以為研究有限應(yīng)變狀態(tài)與構(gòu)造裂縫發(fā)育程度之間的定量關(guān)系提供線索；有限元數(shù)值模擬（Maerten and Maerten，2006）和離散元數(shù)值模擬（尹宏偉等，未發(fā)表成果）進(jìn)行實驗?zāi)Ｐ蛻?yīng)力應(yīng)變分析的方法也可以與本文物理模擬實驗有限應(yīng)變分析的方法結(jié)合起來進(jìn)行構(gòu)造裂縫的定量預(yù)測；利用CT掃描和數(shù)字體相關(guān)（DVC,Digital Volume Correlation）技術(shù)進(jìn)行物理模擬實驗的三維有限應(yīng)變分析（Adam et al.，2013）可以實現(xiàn)三維構(gòu)造裂縫定量預(yù)測。

參考文獻(xiàn)(References)：

陳波，田崇魯.1998.儲層構(gòu)造裂縫數(shù)值模擬技術(shù)的應(yīng)用實例[J].石油學(xué)報，19(4): 50-54.

丁中一，錢祥麟，霍紅，等.1998.構(gòu)造裂縫定量預(yù)測的一種新方法——二元法[J].石油與天然氣地質(zhì)，19(1): 1-7,14.

高帥，曾聯(lián)波，馬世忠，等.2015.致密砂巖儲層不同方向構(gòu)造裂縫定量預(yù)測[J].天然氣地球科學(xué)，26(3): 427-434.

郭迎春，薛峰，汪必峰.2010.儲層應(yīng)力應(yīng)變與裂縫體積密度定量關(guān)系研究[J].油氣地質(zhì)與采收率，17(4): 87-90.

賈東，陳竹新，羅良，等.2007.斷層相關(guān)褶皺的磁組構(gòu)與有限應(yīng)變:川西岷江沖斷構(gòu)造的實例分析[J].自然科學(xué)進(jìn)展，17(2): 188-195.

秦啟榮，蘇培東.2006.構(gòu)造裂縫類型劃分與預(yù)測[J].天然氣工業(yè)，26 (10): 33-36.

沈國華.2008.有限元數(shù)值模擬方法在構(gòu)造裂縫預(yù)測中的應(yīng)用[J].油氣地質(zhì)與采收率，15(4): 24-29.

沈禮，賈東，尹宏偉，等.2012.基于粒子成像測速（PIV）技術(shù)的褶皺沖斷構(gòu)造物理模擬[J].地質(zhì)論評，58(3): 471-480.

宋惠珍，曾海容，孫君秀，等.1999.儲層構(gòu)造裂縫預(yù)測方法及其應(yīng)用[J].地震地質(zhì)，21(3): 205-213.

王鵬威，陳筱，龐雄奇，等.2014.構(gòu)造裂縫對致密砂巖氣成藏過程的控制作用[J].天然氣地球科學(xué)，25(2): 185-191.

文世鵬，李德同.1996.儲層構(gòu)造裂縫數(shù)值模擬技術(shù)[J].石油大學(xué)學(xué)報（自然科學(xué)版），20(5): 17-24.

曾錦光，羅元華，陳太源.1982.應(yīng)用構(gòu)造面主曲率研究油氣藏裂縫問題[J].力學(xué)學(xué)報，(2): 202-206.

曾聯(lián)波，高春宇，漆家福，等.2008.鄂爾多斯盆地隴東地區(qū)特低滲透砂巖儲層裂縫分布規(guī)律及其滲流作用[J].中國科學(xué)（D輯），38（增刊I）：41-47.

周新桂，操成杰，袁嘉音.2003.儲層構(gòu)造裂縫定量預(yù)測與油氣滲流規(guī)律研究現(xiàn)狀和進(jìn)展[J].地球科學(xué)進(jìn)展，18(3): 398-404.

周新桂，張林炎，范昆.2007.含油氣盆地低滲透儲層構(gòu)造裂縫定量預(yù)測方法和實例[J].天然氣地球科學(xué)，18(3): 328-333.

Adam J,Klinkmüller M,Schreurs G,et al.2013.Quantitative 3D strain analysis in analogue experiments simulating tectonic deformation: Integration of X-ray computed tomography and digital volume correlation techniques [J].Journal of Structural Geology,55: 127-149.

Adam J,Urai J L,Wieneke B,et al.2005.Shear localisation and strain distribution during tectonic faulting--new insights from granular-flow experiments and high-resolution optical image correlation techniques [J].Journal of Structural Geology,27: 283-301.

Adrian R J.1991.Particle-imaging techniques for experimental fluid mechanics [J].Annual Review of Fluid Mechanics,23: 261-304.

Bernard S,Avouac J-P,Dominguez S,et al.2007.Kinematics of fault-related folding derived from a sandbox experiment [J].Journal of Geophysical Research,112,B03S12,doi: 10.1029/2005JB004149.

Bonini M.2007.Deformation patterns and structural vergence in brittle-ductile thrust wedges: An additional analogue modeling perspective [J].Journal of Structural Geology,29: 141-158.

Bonini L,Basili R,Toscani G,et al.2015.The role of pre-existing discontinuities in the development of extensional faults: An analog modeling perspective [J].Journal of Structural Geology,74: 145-158.

Boutelier D,Oncken O and Cruden A R.2014.Trench-parallel shortening in the forearc caused by subduction along a seaward-concave plate boundary: Insights fromanalogue modelling experiments [J].Tectonophysics,611(1): 192-203.

Costa E and Vendeville B C.2002.Experimental insights on the geometry and kinematics of fold-and-thrust belts above weak,viscous evaporitic decollement [J].Journal of Structural Geology,24: 1929-1739.

Cotton J T and Koyi H A.2000.Modeling of thrust fronts above ductile and frictional detachments: Application to structures in the Salt Range and Potwar Plateau,Pakistan [J].Geological Society of America Bulletin,112(3): 351-363.

Couzens-Schultz B A,Vendeville B C and Wiltschko D V.2003.Duplex style and triangle zone formation: insights from physical modeling [J].Journal of Structural Geology,25: 1623-1644.

Cruz L,Teyssier C,Perg L,et al.2008.Deformation,exhumation,and topography of experimental doubly-vergent orogenic wedges subjected to asymmetric erosion [J].Journal of Structural Geology,30: 98-115.

Daeron M,Avouac J-P and Charreau J.2007.Modeling the shortening history of a fault tip fold using structural and geomorphic records of deformation [J].Journal of Geophysical Research,112,B03S13,doi: 10.1029/2006JB004460.

Dooley T P,Jackson M P A,Jackson C A-L,et al.2015.Enigmatic structures within salt walls of the Santos Basin-Part 2: Mechanical explanation from physical modelling [J].Journal of Structural Geology,75: 163-187.

Fincham A M and Spedding G R.1997.Low cost high resolution DPIV for measurement of turbulent fluids [J].Experiments in Fluids,23: 449-462.

Graveleau F and Dominguez S.2008.Analogue modelling of the interactions between tectonics,erosion and sedimentation in foreland thrust belts [J].Comptes Rendus Geoscience,340(5): 324-333.

Hall A S,Wood D W,Ibraim E,et al.2010.Localised deformation patterning in 2D granular materials revealed by digital image correlation [J].Granular Matter,12: 1-14.

Hubbert M K.1937.Theory of scale models as applied to the study of geologic structure [J].Geological Society of America Bulletin,48(10): 1459-1519.

Knappett J A,Haigh S K and Madabhushi S P G.2006.Mechanisms of failure for shallow foundations under earthquake loading [J].Soil Dynamics and Earthquake Engineering,26: 91-102.

Koyi H.1995.Mode of internal deformation in sand wedges [J].Journal of Structural Geology,17: 293-300.

Koyi H.1997.Analogue modeling: from a qualitative to a quantitative technique--a historical outline [J].Journal of Petroleum Geology,20(2): 223-238.

Leever K A,Gabrielsen R H,Sokoutis D,et al.2011.The effect of convergence angle on the kinematic evolution of strain partitioning in transpressional brittle wedges: Insight from analog modeling and highresolution digital image analysis [J].Tectonics,30(2): 182-187.

Leturmy P,Mugnier J L and Vinour P.2000.Piggyback basin development above a thin-skinned thrust belt with two detachment levels as a function of interactions between tectonic and superficial mass transfer: the case of the Subandean Zone (Bolivia) [J].Tectonophysics,320: 45-67.

Lorenz J C,Lawrence W T and Warpinski N R.1991.Regional fractures I: a mechanism for the formation of regional fractures at depth in flat-lying reservoirs [J].AAPG Bulletin,75(11): 1714-1737.

Maerten L and Maerten F.2006.Chronologic modeling of faulted and fractured reservoirs using geomechanically based restoration: Technique and industry applications [J].AAPG Bulletin,90(8): 1201-1226.

Massoli D,Koyi H A and Barchi M R.2006.Structural evolution of a fold and thrust belt generated by multiple decollement: analogue models and natural examples from the Northern Apennines (Italy) [J].Journal of Structural Geology,28: 185-199.

McQuillan H.1973.Small-scale fracture density in Asmari formation of southwest Iran and its relation to bed thickness and structural setting [J].AAPG Bulletin,57(4): 2367-2385.

Murray Jr G H.1968.Quantitative Fracture Study--Sanish Pool,Mckenzie County,North Dakota [J].AAPG Bulletin,52(1): 57-65.

Narr W.1991.Fracture density in the deep subsurface: techniques with application to Point Arguello Oil Field [J].AAPG Bulletin,75(8): 1300-1323.

Narr W and Lerche L.1984.A method for estimating subsurface fracture density in core [J].AAPG Bulletin,68(5): 637-648.

Nelson R A.1985.Geologic Analysis of Naturally Fractured Reservoirs [M].Houston,Texas: Gulf Publishing Company: 1-320.

Price N J.1966.Fault and Joint Development in Brittle and Semi-Brittle Rock [M].Oxford: Pergamon Press: 1-176.

Reiter K,Kukowski N and Ratschbacher L.2011.The interaction of two indenters in analogue experiments and implications for curved fold-and-thrust belts [J].Earth and Planetary Science Letters,302: 132-146.

Rosenau M,Lohrmann J and Oncken O.2009.Shocks in a box: An analogue model of subductionearthquake cycles withapplicationto seismotectonic forearc evolution [J].Journal of Geophysical Research Atmospheres,114(B1): 10361.

Santimano T,Rosenau M and Oncken O.2015.Intrinsic versus extrinsic variability of analogue sand-box experiments-Insights from statistical analysis of repeated accretionary sand wedge experiments [J].Journal of Structural Geology,75: 80-100.

Scarano F and Riethmuller M L.2000.Advances in iterative multigrid PIV image processing [J].Experiments in Fluids,29: 51-60.

Simoes M,Avouac J P,Chen Y G,et al.2007.Kinematic analysis of the Pakuashan fault tip fold,west central Taiwan: Shortening rate and age of folding inception [J].Journal of Geophysical Research,112,B03S14,doi: 10.1029/2005JB004198.

Storti F,Salvini F and McClay K.1997.Fault-related folding in sandbox analogue models of thrust wedges [J].Journal of Structural Geology,19: 583-602.

White D J,Take W A and Bolton M D.2003.Soil deformation measurement using particle image velocimetry (PIV) and photogrammetry [J].Geotechnique,53(7): 619-631.

Wieneke B.2001.PIV adaptive multi-pass correlation with deformed interrogation windows [J].PIV Challenge,2001: 1-6.

Wolf H,Konig D and Triantafyllidis T.2003.Experimental investigation of shear band patterns in granular material [J].Journal of Structural Geology,25: 1229-1240.

Wu J E,Mcclay K and Frankowicz E.2015.Niger Delta gravity-driven deformation above the relict Chain and Charcot oceanic fracture zones,Gulf of Guinea: Insights from analogue models [J].Marine & Petroleum Geology,65: 43-62.

Yamada Y,Baba K and Matsuoka T.2006.Analogue and numerical modeling of accretionary prisms with a decollement in sediments [G].In: Buiter S J H,Schreurs G.eds.Analogue and Numerical Modelling of Crustal-ScaleProcesses.Geological Society,London,Special Publications,253: 169-183.

Zeng L B and Li XY.2009.Fractures in sandstone reservoirs with ultra-low permeability: A case study of the Upper Triassic Yanchang Formation in the Ordos Basin,China [J].AAPG Bulletin,93(4): 461-477.

Structural Analogue Modelingand PIVFinite Strain Analysis: Implicationsto Tectonic Fracture Prediction

SHEN Li1,2,JIA Dong1*,YIN Hongwei1,WEI Dongtao2,4,CHEN Zhuxin3,SUN Chuang1,CUI Jian1
1.Institute of Energy Sciences,School of Earth Sciences and Engineering,Nanjing University,Nanjing 210023,China 2.Key Laboratory of Reservoir Characterization,China National Petroleum Corporation (CNPC)，Lanzhou 730020,China 3.Research Institute of Petroleum Exploration and Development,Petrochina，Beijing 100083,China 4.Research Institute of Petroleum Exploration and Development-Northwest (NWGI),Petro Lanzhou 730020,China

Abstract:The development of structural fractures and finite strain state are closely related in the oil and gas reservoir.In order to explore a new method of finite strain analysis and structural fracture prediction,this study designed and performed a set of sandbox experiments.The particle image velocimetry (PIV) technique was applied to quantitatively analyze the experimental process.The experimental model was a unilateral indentation model with a viscous layer in the vertical direction,and the results represent abook=172,ebook=175macroscopic box fold.Using the PIV technique,we can obtain the displacement field data during the deformation process in each stage of the experiments,calculate the incremental strain in each step,and analyze the finite strain state during the whole deformation process from the initial stage to the time after the fold was formed.Furthermore,the genesis mechanism and distribution can be discussed to quantiatively predict the fractures.At the beginning of the deformation,the finite strain was weak (about 4-8%) and widely distributed.Weak compressive linear strain was represented in the direction of indentation and weak tensile linear strain in the vertical direction.These phenomena were interpreted as the results of the pure shear deformation,such as the thickening and parallel shortening of the layers before the folding and faulting initiated.It is also the main mechanism of the formation of the regional tension fissures and shear cracks.The finite strain was localized in the fault zone and the adjacent area when the fold and fault formed.Strong shear and tensile shear strain (up to 20%) was revealed.These were the results of the simple shear deformation along the fault plane.It is also the main mechanism of the formation of the local shear fractures and tensile shear cracks.

Key words:analogue modeling; particle image velocimetry (PIV); finite strain analysis; tectonic fracture prediction

Corresponding author:JIA Dong,Professor; E-mail: djia@nju.edu.cn

*通訊作者：賈東，男，1960年生，教授，博士生導(dǎo)師，主要從事含油氣盆地構(gòu)造和活動構(gòu)造研究；E-mail: djia@nju.edu.cn

作者簡介：沈禮，男，1986年生，博士研究生，構(gòu)造地質(zhì)學(xué)專業(yè)；E-mail: nju.shenli2010@gmail.com

基金項目：江蘇省科技項目（BE2013115）；國家科技重大專項（2011ZX05009-001；2011ZX05003-002）聯(lián)合資助

收稿日期：2015-11-30；修回日期：2016-03-07

DOI:10.16108/j.issn1006-7493.2015238

中圖分類號：P618.13

文獻(xiàn)標(biāo)識碼：A

文章編號：1006-7493（2016）01-0171-12