董少群 曾聯(lián)波* Xu Chaoshui曹 菡 王圣嬌 呂文雅(①中國石油大學(xué)(北京)地球科學(xué)學(xué)院，北京 102249； ②油氣資源與探測國家重點實驗室，北京102249； 阿德萊德大學(xué)土木、環(huán)境和采礦工程學(xué)院,澳大利亞阿德萊德 5005； 中國科學(xué)院對地觀測與數(shù)字地球科學(xué)中心，北京 100094)

1 引言

裂縫是指巖石受到構(gòu)造變形發(fā)生破裂作用或物理成巖作用形成的沒有明顯位移的面狀不連續(xù)體，常呈網(wǎng)狀分布[1]。裂縫在很多領(lǐng)域均有涉及，其空間展布研究對于油氣勘探開發(fā)、地?zé)衢_發(fā)、采礦、巖土工程、水文地質(zhì)等領(lǐng)域均有重要意義。在油氣田勘探開發(fā)領(lǐng)域常稱其為“裂縫”，在巖土工程、采礦、水文地質(zhì)等領(lǐng)域則習(xí)慣稱其為“裂隙”，兩者中文名稱雖不同，但其所指一般等同[2]。為便于敘述，本文統(tǒng)一將其稱為裂縫。裂縫是巖體中流體移動的通道，如裂縫網(wǎng)絡(luò)是儲層中油氣、地?zé)嵯到y(tǒng)中熱液、水文地質(zhì)工程中地下水等流體運移的通道。

裂縫建模是一種常用的研究裂縫空間展布的手段，可分為確定性建模和隨機建模。確定性建模是根據(jù)已知信息建立確定的裂縫模型，如通過地震資料解釋出規(guī)模較大的裂縫，該方法一般不適用于規(guī)模較小的裂縫，且不能較好地綜合利用多種資料。裂縫隨機建模則是利用裂縫的先驗信息，通過隨機模擬方式生成可選的相同概率裂縫模型。這類方法不僅滿足已知點的裂縫統(tǒng)計學(xué)特征，而且承認(rèn)未知區(qū)域裂縫發(fā)育的隨機性，較好地尊重了裂縫模擬不確定性的客觀事實[3]。近年來裂縫隨機建模方法應(yīng)用較為廣泛，且已在油氣開采、地?zé)衢_發(fā)、采礦等領(lǐng)域均取得較好的應(yīng)用效果。

目前常用的隨機裂縫建模方法大致可分為五類： ①基于空間剖分的裂縫建模； ②離散裂縫網(wǎng)絡(luò)建模； ③基于變差函數(shù)的裂縫建模； ④基于多點地質(zhì)統(tǒng)計學(xué)的裂縫建模； ⑤基于分形特征迭代的裂縫建模。本文結(jié)合實例介紹五類方法的原理、優(yōu)缺點及適用性等，同時討論這些方法在應(yīng)用中存在的問題，并在此基礎(chǔ)上總結(jié)了裂縫隨機建模方法研究的四個發(fā)展趨勢。

2 基于空間剖分的裂縫建模

假設(shè)巖石中單個裂縫為平面多邊形，所有裂縫多邊形延伸后的平面便可將巖石剖分成多個巖塊，即進行空間剖分。裂縫建模時，若空間剖分合理，則這些平面上的某些多邊形即為裂縫?；诳臻g剖分的裂縫建模方法便是這種建模思路的一種實現(xiàn)，常用模型包括正交模型、馬賽克鑲嵌模型和Veneziano模型等[4]。

2.1 正交模型

正交模型是針對巖石裂縫系統(tǒng)進行幾何模擬的一種建模方法，其早期應(yīng)用于水文地質(zhì)等領(lǐng)域[5]。實際應(yīng)用中的正交模型多為三維模型，初期的正交模型為無邊界正交模型，它將巖石裂縫系統(tǒng)中的裂縫簡化為三個方向上的無限正交面(圖1a的灰色平面)。不同方向裂縫面間距分別為S1、S2、S3，可通過隨機模擬獲得。圖1b為澳大利亞塔斯馬尼亞的棋盤石[6]，通過野外露頭測量統(tǒng)計不同方向的裂縫間距，分別獲得S1、S2、S3的分布函數(shù)，隨后便可通過蒙特卡洛隨機模擬得到隨機正交裂縫模型。該方法適用于正交的高角度裂縫系統(tǒng)建模，且裂縫規(guī)模相對于研究區(qū)范圍足夠大。該方法的優(yōu)點是簡單易行，但實際裂縫系統(tǒng)并非都為如此規(guī)整正交裂縫網(wǎng)絡(luò)(圖1d中的延河剖面露頭、圖1f中的加利福尼亞的碳酸鹽巖露頭剖面[7])。針對此類問題，有學(xué)者[4，8，9]將無邊界正交模型進行了改進，提出有邊界正交模型(圖1c、圖1e)。將圖1a中剖分出的小網(wǎng)格進行編號，然后依據(jù)一定的規(guī)則選擇某些小網(wǎng)格或小網(wǎng)格的一部分作為裂縫。

圖1 正交模型及與其相似露頭剖面(a)無邊界正交模型； (b)塔斯馬尼亞較規(guī)整的棋盤石(正交裂縫)[6]； (c)有邊界正交模型；(d)延河剖面露頭； (e)有邊界正交模型； (f)加利福尼亞圣巴巴拉市碳酸鹽巖露頭剖面[10]

2.2 馬賽克鑲嵌模型

二維馬賽克鑲嵌模型，先在平面上生成一系列散點，然后通過相鄰兩點的中垂線將平面分割成彼此鑲嵌的凸多邊形“馬賽克”[11]，圖2a是其俯視圖。當(dāng)給予每個多邊形一定高度后便可形成三維馬賽克鑲嵌模型(圖2a)，多邊形柱狀體相鄰面即為裂縫，這種方法可較好地模擬柱狀節(jié)理(圖2b)。另一種三維馬賽克實現(xiàn)方式如圖2c所示，它直接在空間中生成一系列散點，然后通過相鄰兩點的中垂面將三維空間剖分維彼此鑲嵌的凸多面體“馬賽克”。此時巖石被剖分成斷塊或基質(zhì)巖塊，按一定規(guī)則選定凸多面體的部分面為裂縫，所形成的裂縫網(wǎng)絡(luò)即為所得的裂縫模型。野外露頭及巖心常見此類裂縫，如圖2d所示。在利用這種方法建立裂縫模型時，生成平面和空間中的散點十分關(guān)鍵，因為它直接決定了模型是否可用。

2.3 Veneziano模型

Veneziano模型也稱泊松平面模型(Poisson flat model)，是由Veneziano于1978年提出，其建立過程如圖3所示，可分為三步[4，12]： ①泊松面過程(圖3a)，在空間中生成隨機平面Ai(i=1,2,…,n)，n為預(yù)定義的參數(shù)； ②泊松線過程(圖3b)，通過泊松過程生成平面Ai上的隨機直線，將平面分割成多邊形集Pi(i=1,2,…,m)，m為第i個平面上的多邊形的個數(shù)； ③依概率選取多邊形作為裂縫面(圖3c)。整個模型生成過程由裂縫密度決定[12]，步驟①反映了主應(yīng)力走向及研究區(qū)的整體裂縫密度[13]，由P30、P31或P33(表1)控制，步驟②、步驟③通過裂縫多邊形的大小和位置反映裂縫密度的變化，步驟②由P20、P21或P22(表1)控制。前面提及的正交模型、馬賽克鑲嵌模型均可視為Veneziano模型的一種特例。

圖2 馬賽克鑲嵌模型及與之相似的露頭剖面(a)二維馬賽克鑲嵌擴展模型； (b)濟州島柱狀節(jié)理； (c)三維馬賽克鑲嵌模型； (d)塔河油田一間房組的巖心

圖3 Veneziano模型(a)泊松面過程確定隨機面； (b)泊松線過程確定隨機線； (c)確定裂縫面； (d)三維裂縫模型

表1 裂縫密度、強度、孔隙度參數(shù)

該模型提出后，一些學(xué)者陸續(xù)對該模型做了改進。Dershowitz等[8，9]將該模型的后兩步進行了改進，提出了Dershowitz模型。Einstein、Ivanova等[4]在Veneziano模型(圖3d)的基礎(chǔ)上進行了改進，通過將多邊形進行旋轉(zhuǎn)及篩選等操作將Veneziano模型轉(zhuǎn)換為離散裂縫網(wǎng)絡(luò)模型。Ivanova等也嘗試對步驟②進行改進，使用1.2節(jié)中介紹的Voronoi圖法將平面進行剖分[13]，提高了生成符合條件的裂縫多邊形的效率。Mosser等[14]也對該模型進行了研究，認(rèn)為直接利用泊松線過程建立的裂縫模型過度估計了裂縫的連通性，并且針對此問題提出了STIT方法。STIT方法在泊松線過程的基礎(chǔ)上，對其剖分出的網(wǎng)格進行進一步剖分，Mosser等[14]認(rèn)為該方法可以緩解過度估計裂縫連通性的問題。

3 離散裂縫網(wǎng)絡(luò)建模

離散裂縫網(wǎng)絡(luò)(Discrete fracture networks，DFN)最早由Baecher等[12，15]引入。如圖4a所示，Baecher模型假設(shè)裂縫為一圓盤(圓形或橢圓形)，各圓盤之間可以交叉，是一種基于示性點過程的隨機建模方法[16，17]。其實現(xiàn)過程分為兩步：點過程、示性過程。點過程確定裂縫的中心位置； 示性過程確定點的屬性，如裂縫大小、傾角、傾向、開度、滲透率等屬性[17]。

自Baecher模型被引入以來，便得到了廣泛應(yīng)用，并且出現(xiàn)了不同形式的改進，主要體現(xiàn)在以下三個方面。

(1)裂縫面形狀表示方法的改進。三維裂縫面最初假設(shè)為圓盤，因為裂縫產(chǎn)生時最初近似為圓形或橢圓形，隨著破裂過程的繼續(xù)，受層面等多方面限制，其形狀變得相對不規(guī)則，故后續(xù)有學(xué)者相繼將其改進為多邊形等，并對裂縫的期次(如裂縫多邊形之間的的切割、限制關(guān)系)做了相關(guān)改進(圖4b～圖4e)。圖4b為增強Baecher模型，它將橢圓換成近似橢圓的多邊形，而且裂縫可被其他裂縫截斷； 圖4c為BART(Baecher Algorithm Revised Terminations)模型，其對裂縫終止過程做了進一步改進[9]； 圖4d為泊松矩形模型，是一種簡化的增強Baecher模型，其裂縫通過矩形表示，而不再是原來近似圓形的多邊形[4]； 圖4e中模型為隨機多邊形模型，它將裂縫面定義為邊數(shù)固定的隨機多邊形[17，20]，目前很多裂縫建模商業(yè)軟件中已采用該方案，可預(yù)設(shè)多邊形邊數(shù)，如FracMan、Petrel等。

(2)裂縫中心位置的產(chǎn)生規(guī)則的改進。裂縫位置的確定一般通過穩(wěn)態(tài)泊松過程來實現(xiàn)， 為擴展其適用范圍，其他方式的裂縫位置產(chǎn)生規(guī)則相繼被引入，如基于非穩(wěn)態(tài)泊松過程、馬爾科夫鏈蒙特卡洛(Markov Chain MonteCarlo，MCMC)[21]、分形[18，22]、Cluster過程[17]、Cox過程[17]等(圖4f～圖4j)。圖4f中模型為最近鄰模型，是一種半概率的、基于模式的模型。裂縫易分布于斷層或大規(guī)模裂縫附近，該模型模擬該趨勢，優(yōu)先在斷層或大規(guī)模裂縫附近生成裂縫[4]。圖4g中模型為戰(zhàn)區(qū)模型，也是一種半概率的、基于模式的模型。它通過優(yōu)先在近似平行的裂縫之間生成裂縫來模擬剪切帶的幾何形態(tài)[23]。圖4h～圖4j分別對應(yīng)Levy-Lee分形模型[24]、母女模型(Parent-daughter model)[25，26]、二元分形裂縫網(wǎng)絡(luò)(Binary fractal fracture network，BFFN)模型[18，22]，它們是考慮裂縫位置的分形特征的改進DFN模型。

圖4 離散裂縫網(wǎng)絡(luò)模型[4，9，11，13，17-19](a)Baecher模型； (b)增強Baecher模型； (c)BART模型； (d)泊松矩形模型； (e)隨機多邊形模型； (f)最近鄰模型； (g)戰(zhàn)區(qū)模型；(h)Levy-Lee分形模型； (i)母女模型； (j)二元分形裂縫網(wǎng)絡(luò)模型； (k)非平面帶模型； (l)GEOFRAC模型

(3)裂縫屬性分布的改進。很多學(xué)者對裂縫的走向、傾向、傾角、開度、孔隙度、滲透率、間距等參數(shù)的分布做了大量研究工作。如描述裂縫長度常使用指數(shù)分布[27]、對數(shù)正態(tài)分布[15，28]及伽馬分布[29]等； 描述裂縫走向的有von-Mises分布[17]； 描述裂縫傾向、傾角的有Fisher分布； 描述開度的有指數(shù)分布、對數(shù)正態(tài)分布[30]、與裂縫長度呈線性或非線性關(guān)系[31]等； 描述裂縫間距的有指數(shù)分布[12]、對數(shù)正態(tài)分布等； 描述裂縫跡長的有對數(shù)正態(tài)分布[12]、 指數(shù)分布等。

對DFN模型進行改進時，常兼顧上述多個方面，如圖4h中的非平面帶模型，是一種半隨機模型，其裂縫位置和走向隨設(shè)定曲面變化[4]，對裂縫面的形狀及裂縫中心位置的產(chǎn)生規(guī)則均有所改進。

建模方法的主要原理相同，即示性點過程。除了基于示性點過程的離散裂縫網(wǎng)絡(luò)建模外，Ivanova等[4，13，32]在Veneziano模型的基礎(chǔ)上進行了改進，也可以生成離散裂縫網(wǎng)絡(luò)模型，稱為GEOFRAC模型(圖4l)。限于篇幅，本文僅選取幾組有代表性的DFN模型進行介紹。

3.1 Baecher模型

Baecher模型也稱為隨機圓盤模型(Random-disk model)[12]，它通過點過程和示性過程共同確定裂縫網(wǎng)絡(luò)中的各條裂縫。原始Baecher模型有五個假設(shè)條件[4，12]： ①單個裂縫為平面圓盤； ②裂縫中心位置符合隨機分布，且各裂縫中心位置彼此獨立； ③圓盤半徑符合對數(shù)正態(tài)分布； ④裂縫半徑與傾角不相關(guān)(統(tǒng)計角度上獨立)； ⑤裂縫半徑與裂縫中心位置不相關(guān)(統(tǒng)計角度上獨立)。由于構(gòu)造應(yīng)力作用，實際儲層的裂縫常具有一個或多個優(yōu)勢產(chǎn)狀(傾向、傾角)?；诖?，Baecher模型建立過程如圖5所示，分為四步。

(1)利用極點圖劃分裂縫組系。如圖5a所示，半徑對應(yīng)裂縫法線與鉛垂線夾角(即裂縫傾角)的正弦，方位角對應(yīng)裂縫法線在水平面的投影與正北方的夾角(即傾向)，顏色代表裂縫的數(shù)量，暖色調(diào)表示裂縫數(shù)量多，可以看出存在兩組高角度裂縫，走向分別為北北東和北北西。

(2)點過程確定裂縫位置。以第①組裂縫為例，通過下式中均勻分布[33]確定隨機數(shù)(x,y,z)，即為裂縫中心位置坐標(biāo)

(1)

式中V為研究區(qū)體積。該組裂縫圓盤的中心位置如圖5b所示。

(3)示性過程確定每條裂縫的屬性(如傾向、傾角、走向、大小、開度等)。通過蒙特卡洛隨機模擬產(chǎn)生的隨機數(shù)作為分位數(shù)，從累計分布函數(shù)上獲得分位數(shù)對應(yīng)的裂縫屬性值。裂縫產(chǎn)狀(傾向、傾角)模擬可通過模擬裂縫面法線來實現(xiàn)[17，33-36]，確定裂縫面法線后再將其轉(zhuǎn)換為裂縫面的傾向和傾角。除此之外，產(chǎn)狀也可通過直接模擬裂縫面的傾向、傾角來實現(xiàn)[17，19，33-35]，即首先擬合裂縫傾向、傾角所符合的分布，然后直接對兩個分布進行隨機模擬，獲取相應(yīng)的傾向和傾角數(shù)值。

下面將以通過裂縫面法線確定裂縫產(chǎn)狀的方式進行介紹。利用von-Mises分布(式(2))擬合裂縫傾向，利用均勻分布(式(3))擬合裂縫傾角，利用對數(shù)正態(tài)分布(式(4))擬合裂縫圓盤半徑和裂縫開度。

(2)

(3)

式中a、b分別表示裂縫傾角的最大值和最小值。

(4)

式中：σ體現(xiàn)圓盤半徑或開度的變化程度，σ越大，分布函數(shù)越平緩；μ決定半徑或開度分布的峰值位置。

通過圖5a中第①組裂縫的參數(shù)統(tǒng)計結(jié)果進行概率密度估計可獲得式(2)～式(4)中的參數(shù)。隨后便可得到對應(yīng)的累積概率密度曲線，如圖5c～圖5f所示。

通過蒙特卡洛隨機模擬產(chǎn)生四個隨機數(shù)作為四個分布的分位數(shù)，然后找到對應(yīng)累積概率密度曲線的橫坐標(biāo)即為裂縫屬性，如圖11中虛線所示，傾向、傾角、半徑、開度分別為61°、84°、1.2m、0.0192mm。重復(fù)上述操作3次后，生成的裂縫如圖5g所示，顏色表示開度的大小。繼續(xù)重復(fù)上述過程直至滿足終止條件即可生成離散裂縫網(wǎng)絡(luò)模型。終止條件多與裂縫密度有關(guān)，如表1所示。常用的終止條件有： ①單位體積中裂縫的條數(shù)P30達(dá)到預(yù)設(shè)值； ②單位體積中裂縫的面積P32達(dá)到預(yù)設(shè)值； ③單位體積中裂縫的體積P33達(dá)到預(yù)設(shè)值。這些預(yù)設(shè)值常通過巖心、測井、露頭等測得的裂縫密度給出。完成第一組裂縫模擬后，用同樣的方式可以實現(xiàn)另一組裂縫的建模。

圖5 Baecher模型建模流程(a)極點圖確定裂縫組系； (b)點過程確定裂縫中心位置； (c)示性過程確定裂縫傾向、傾角；(d)示性過程確定裂縫半徑長度； (e)示性過程確定裂縫開度； (f)第①組裂縫模擬結(jié)果

(4)質(zhì)量控制(QC)。裂縫網(wǎng)絡(luò)模型建完后，很重要的一步就是質(zhì)量控制，而裂縫密度則是常用的一類衡量參數(shù)，如統(tǒng)計井軌跡通過P10(單位長度上裂縫的條數(shù))或P11(單位長度上裂縫的長度)是否與井?dāng)?shù)據(jù)相對吻合，如果不吻合，則需要對模型進行調(diào)整。

3.2 最近鄰模型

按裂縫形成時構(gòu)造應(yīng)力場作用范圍，可將裂縫分為區(qū)域構(gòu)造裂縫和局部構(gòu)造裂縫。局部構(gòu)造裂縫又可分為與斷層相關(guān)的裂縫和與褶皺相關(guān)的裂縫。與斷層相關(guān)的裂縫存在于斷層附近，與之伴生，其發(fā)育的應(yīng)力場與斷層一致，且隨著與斷層距離的增大，裂縫的線密度明顯降低[1，37](圖6[38])。最近鄰模型是一種改進的Baecher模型，它對裂縫位置做了改進，假設(shè)裂縫密度與裂縫和主裂縫的距離有關(guān)，它隨著距離的增大呈指數(shù)減小，從而使裂縫更易于分布于斷層或大規(guī)模裂縫附近[4]。

針對上述現(xiàn)象，高志勇等[37]對庫車前陸盆地裂縫與斷層關(guān)系進行了定量分析。圖7顯示裂縫線密度與距斷面距離成較好的指數(shù)關(guān)系，與最近鄰模型假設(shè)相近，故可用于此類裂縫建模及最后的裂縫模型效果檢驗。

圖6 裂縫分布與距斷層距離相關(guān)的剖面[38]

圖7 庫車前陸盆地裂縫與斷層關(guān)系[38，40]

使用最近鄰法建立與斷層相關(guān)的裂縫模型時，裂縫的性質(zhì)可以參考與斷層相關(guān)的裂縫模式[39]。通過分析斷層的應(yīng)力特征及發(fā)育裂縫的力學(xué)特征可以獲得裂縫的走向、傾向、傾角等相關(guān)信息，也可以輔助檢驗輸入模型參數(shù)合理性及建立的裂縫網(wǎng)絡(luò)模型的有效性。

在建立與斷層相關(guān)的裂縫模型時，也可以使用母女模型[25，26]，同樣能達(dá)到與最近鄰模型類似的效果，母女模型除了考慮斷層附近裂縫密度高的特性外，還考慮了裂縫的分形特征，本文將在3.3節(jié)中進行介紹。

3.3 分形模型

利用了裂縫分形特征的DFN模型均可看作是分形模型，認(rèn)為裂縫長度符合冪律分布[33，41]。除考慮裂縫長度的分形特征外，裂縫位置的分形特征也常被引入到DFN模型中，如Levy-Lee分形模型[24](圖4h)、母女模型[25，26](圖4i)、二元分形裂縫網(wǎng)絡(luò)模型[18，22](圖4j)等。

實際應(yīng)用中，裂縫的分形特征也被學(xué)者廣泛認(rèn)同[42-45]。很多學(xué)者認(rèn)為晶質(zhì)巖中的裂縫系統(tǒng)常呈現(xiàn)分形特征，具體體現(xiàn)在兩個方面：裂縫空間分布符合冪律分布； 裂縫長度符合冪律分布。BFFN模型正是基于上述兩點提出的。模型中裂縫中心位置按照“分形破碎”的方法產(chǎn)生，裂縫長度利用隨機過程結(jié)合冪律定理產(chǎn)生[18，22]。其具體步驟如下。

(1)利用“分形破碎”的方法，根據(jù)D2(式(5))確定裂縫中心位置

(5)

(2)示性過程確定裂縫屬性，如圖8h所示。裂縫長度根據(jù)式(6)確定，傾向和傾角依據(jù)Fisher分布(式(7))獲得

lgN=alglmax-algl

(6)

(7)

式中：N為裂縫長度不小于l的裂縫數(shù)量；lmax為最大裂縫的長度；φ、θ分別為裂縫面的傾向、傾角。

3.4 Veneziano改進模型

Ivanova等在離散裂縫建模方法研究和改進方面未沿用Baecher模型中示性點過程的方法，而是在Veneziano模型[4，12]基礎(chǔ)上進行了改進，將圖3d中生成的裂縫模型進行了平移、旋轉(zhuǎn)，從而反映裂縫位置和走向的局部變化，并將其稱為GEOFRAC模型[13]，如圖4l所示。原始GEOFRAC模型在裂縫生成過程中存在一個問題，泊松線過程產(chǎn)生的多邊形中有大量因其形狀不能較好地反映裂縫形態(tài)而被刪除，即裂縫多邊形生成效率低，如圖3b、圖3c所示。針對這一問題，Ivanova等用Delauney三角剖分代替泊松線過程對原始模型進行了改進，有效地改善了裂縫多邊形生成效率低的問題，他們將改進后的模型稱為新GEOFRAC模型[13]。新GEOFRAC模型的生成過程如圖9所示。

圖8 三維BFFN模型構(gòu)建過程[18](a)～(g)裂縫中心位置產(chǎn)生過程； (h)基于示性點過程生成裂縫模型

圖9 新GEOFRAC模型中裂縫生成過程[13](a)通過泊松過程生成平面； (b)平面網(wǎng)格剖分； (c)多邊形隨機平移； (d)多邊形隨機旋轉(zhuǎn)

4 基于變差函數(shù)的裂縫建模

基于變差函數(shù)的建模方法在沉積相、儲層物性等建模中廣泛應(yīng)用[46]，但在裂縫建模中應(yīng)用相對較少。通過變差函數(shù)可以表征空間中不同位置裂縫的相關(guān)性[16，47-49]，因此不少學(xué)者也嘗試將其應(yīng)用于裂縫建模。大部分工作是在2.1節(jié)中Baecher模型的框架內(nèi)展開的，即利用變差函數(shù)模擬裂縫密度或確定裂縫位置，即代替點過程確定裂縫位置[16]。Dowd等[16]、Xu等[47]利用截斷高斯指示模擬(Truncated Gaussian indicator simulation)和截斷多高斯場截斷模擬(Truncated Pluri Gaussian simulation)[50]對裂縫位置進行模擬，并結(jié)合示性過程進行裂縫建模，實驗結(jié)果表明截斷高斯指示模擬適用于裂縫系統(tǒng)相對簡單的模擬，截斷多高斯場截斷模擬可更好地對復(fù)雜裂縫系統(tǒng)中裂縫位置進行模擬。

變差函數(shù)除了可以用于代替點過程外，也可以用于代替示性過程。Koike等[48，49，51，52]利用通過指示變換和變差函數(shù)來模擬裂縫走向、傾向，其文章中將其方法也稱為GEOFRAC[48，49]，但根據(jù)發(fā)表先后，3.4節(jié)中方法要早，因此本文中GEOFRAC如果沒有特別注釋，則代指2.3節(jié)中的方法。Koike等提出的方法的操作步驟如圖10所示，可分為以下三步。

(1)確定小平面中心，如圖10a所示。該方法不直接模擬裂縫，而是先生成尺度較小的小平面，然后再將小平面連接成最終裂縫面。在確定裂縫中心時，先將儲層進行網(wǎng)格化，然后利用序貫高斯模擬(Sequential Gaussian Simulation，SGS)模擬每個網(wǎng)格的裂縫密度，結(jié)合裂縫密度通過隨機過程在各網(wǎng)格中確定是否有小平面中心，若有，則確定其中心位置。

(2)確定小平面走向和方位角(圖10b～圖10e)。將走向和傾角劃分為16個區(qū)域，并進行指示變換，I=(g1,g2,…,gL)，此時L為16，當(dāng)走向和方位角位于第i個區(qū)域時，gi=1，其他元素為0。此處Koike等定義的走向符合右手定則，范圍為[0,2π]。通過普通克里金(Ordinary Kriging，OK)計算某個小平面中心位置處的I,如I=(0.2,0.9,…,0.05)，則其屬于第2個區(qū)域。隨后再根據(jù)該區(qū)域內(nèi)裂縫的走向和方位的分布函數(shù)隨機生成裂縫走向和方位。在普通克里金插值時，為了減少計算量，該方法引入了主成分分析。

(3)將小平面連接形成裂縫面(圖10c～圖10f)。首先判斷兩個小平面是否連接，它可以通過兩個小平面的走向差、傾向差和小平面中心距離判斷是否連接，也可以通過投影后的走向差和中心距離來判斷兩個小平面是否連接。隨后將孤立的小平面轉(zhuǎn)換成圓盤，作為孤立裂縫，將相互連接的小平面連接成裂縫面，如圖10g，可以直接將小平面中心連接成不規(guī)則三角形格網(wǎng)(Triangulated Irregular Network，TIN)，也可以根據(jù)小平面的走向、傾角的均值計算裂縫面，并用小平面的包絡(luò)作為裂縫面的邊界。最終生成的裂縫模型，如圖10h所示。

圖10 基于變差函數(shù)的裂縫建模(據(jù)文獻(xiàn)[48， 49]，有改動)(a)確定小平面中心； (b)～(e)確定小平面走向和方位角； (f)連接小平面形成裂縫面； (g)、(h)離散裂縫網(wǎng)絡(luò)模型

5 基于多點地質(zhì)統(tǒng)計學(xué)的裂縫建模

多點地質(zhì)統(tǒng)計學(xué)較以變差函數(shù)為核心的兩點地質(zhì)統(tǒng)計學(xué)而言，易于表征復(fù)雜的空間結(jié)構(gòu)、再現(xiàn)復(fù)雜目標(biāo)的幾何形態(tài)[46]。多點地質(zhì)統(tǒng)計學(xué)包括迭代的和非迭代的兩大類方法。迭代方法有基于模擬退火的方法[53]、基于吉布斯取樣的后處理迭代方法等，但這類方法普遍存在迭代收斂問題[46]。文獻(xiàn)中提及的多點地質(zhì)統(tǒng)計學(xué)方法多指非迭代的方法，該類方法最早由Guardiano等提出，其通過訓(xùn)練圖像掃描獲取待估點的概率分布，進而模擬沉積相[54]，但在掃描圖像時運算量及內(nèi)存消耗都較大，不同學(xué)者對其嘗試進行改進。目前常用的改進算法有SNESIM(Single Normal Equation Simulation)算法[55， 56]、SIMPAT算法[57-59]、FILTERSIM 算法[60]、DisPAT算法[61,62]等。SNESIM使用一個平衡方程進行概率估計，一次性將訓(xùn)練圖像的條件概率存儲于“搜索樹”中，極大地減少了機時。SIMPAT利用訓(xùn)練圖像獲取地下結(jié)構(gòu)的模式，采用相似性方法對地下儲層進行圖像恢復(fù)和再現(xiàn)。FILTERSIM為基于模式濾波的多點地質(zhì)統(tǒng)計模擬方法。DisPAT算法構(gòu)建數(shù)據(jù)樣式的距離矩陣，應(yīng)用多維尺度分析進行降維，再用K-means聚類建立樣式聚類，后續(xù)模擬與SIMPAT相同[62]。這些方法從不同角度對原始方法改進，以提高效率和拓展應(yīng)用范圍。

多點地質(zhì)統(tǒng)計學(xué)早期主要用于沉積相模擬，近年來正逐漸被應(yīng)用于裂縫建模，對于裂縫建模，野外露頭剖面、地震解釋的斷層等都可以作為裂縫幾何形態(tài)的先驗認(rèn)識，即訓(xùn)練圖像，同時測井解釋或巖心觀察的裂縫可以為硬數(shù)據(jù)。Dowd等[16]嘗試使用SNESIM算法對尤卡山二維裂縫剖面進行模擬，并與其他方法進行對比，認(rèn)為多點地質(zhì)統(tǒng)計學(xué)模擬可以體現(xiàn)裂縫的空間相關(guān)性，并且能模擬復(fù)雜的裂縫網(wǎng)絡(luò)。謝青[63]嘗試使用多線程并行運算提高SNESIM算法效率，并對頁巖儲層天然裂縫進行了三維空間建模。Mohammadmoradi[64]、張烈輝等[65]、Jia等[66]針對FILTERSIM算法在裂縫建模中的問題分別進行了改進，并進行了二維裂縫建模。Liu等[67]使用SIMPAT算法對二維裂縫剖面進行建模并數(shù)模，結(jié)果表明所建立模型可以較好地描述裂縫系統(tǒng)。利用多點地質(zhì)統(tǒng)計學(xué)進行裂縫建模，普遍存在穩(wěn)態(tài)的問題，針對這一問題，Chugunova等[68]提出一種非穩(wěn)態(tài)多點地質(zhì)統(tǒng)計學(xué)的模擬方法，該方法除了使用裂縫圖像作為訓(xùn)練圖像外，還加入了裂縫密度圖作為輔訓(xùn)練圖像進行裂縫建模。隨后，Chugunova等[69]利用多點地質(zhì)統(tǒng)計學(xué)建立裂縫網(wǎng)絡(luò)模型，并結(jié)合動態(tài)數(shù)值模擬對模型進行了驗證。其建模具體流程如圖11所示，具體分為三步：①訓(xùn)練圖像構(gòu)建(圖11a～圖11c),通過遙感圖像獲取巖性及裂縫，并選擇某一相同巖性的區(qū)域作為訓(xùn)練圖像; ②確定硬數(shù)據(jù)； 利用遙感數(shù)據(jù)解釋明顯的斷層作為硬數(shù)據(jù)； ③結(jié)合訓(xùn)練圖像和硬數(shù)據(jù)對待估點進行預(yù)測； ④通過生產(chǎn)數(shù)據(jù)驗證模型的有效性。

圖11 基于多點地質(zhì)統(tǒng)計裂縫建模(據(jù)文獻(xiàn)[69]，有改動)(a)遙感圖像； (b)提取的不同巖性區(qū)域的裂縫； (c)訓(xùn)練圖像砂巖4對應(yīng)的裂縫； (d)硬數(shù)據(jù)(以(c)中

6 基于分形特征迭代的裂縫建模

基于分形特征[70，71]迭代的裂縫建模的常用方法有迭代函數(shù)系統(tǒng)(Iterated Function System， IFS)[72-74]、L-System、謝爾賓斯基三角形等，本章重點介紹前兩種。

IFS將待生成的圖像看作由許多與整體相似的(自相似)或經(jīng)過一定變換與整體相似的(自仿射)小塊拼貼而成。而裂縫網(wǎng)絡(luò)常表現(xiàn)出自相似性或自仿射性，所以IFS也被用于裂縫網(wǎng)絡(luò)模擬。Acuna等[72，75]利用IFS對二維和三維裂縫系統(tǒng)進行建模，如圖12a、圖12b所示，該方法首先對初始點群進行一系列迭代數(shù)值變換，隨后每次迭代中，應(yīng)用該系統(tǒng)中的函數(shù)對點群進行轉(zhuǎn)換、映射、旋轉(zhuǎn)、收縮、扭曲。當(dāng)圖像中的點群符合分形目標(biāo)時便終止迭代，此時即可獲得最終裂縫分布圖像[46，75]。周德華等[74]將IFS應(yīng)用于油氣藏天然裂縫網(wǎng)絡(luò)建模，方法實施包括三部分： ①裂縫網(wǎng)絡(luò)的迭代生成，根據(jù)生成模型的分形維數(shù)Df與實際地層裂縫參數(shù)空間Dr的分形維數(shù)的差異決定迭代發(fā)生器的選擇和迭代的終止； ②裂縫網(wǎng)絡(luò)中流體流動的有限元模擬； ③通過與試井資料的匹配關(guān)系，利用模擬退火算法優(yōu)化IFS參數(shù)，進而優(yōu)化天然裂縫網(wǎng)絡(luò)的幾何結(jié)構(gòu)。

L-System是描述植物生長的數(shù)學(xué)模型，其基本思想可解釋為理想化的樹木生長過程：從一條樹枝(種子)開始，發(fā)出新的芽枝，而發(fā)過芽枝的枝干又全部發(fā)新芽枝，直至最后長出葉子。該方法與裂縫生長過程較為相近，故可嘗試應(yīng)用于裂縫建模。根據(jù)L-System 的基本思想，先建立二維樹的數(shù)學(xué)模型：從樹干開始，然后沿著樹干逐漸擴展到連接的樹枝，再以遞歸的方式進行同樣的過程，該過程持續(xù)到最終分枝。Zhou等[76，77]結(jié)合微地震數(shù)據(jù)、生產(chǎn)數(shù)據(jù)對非常規(guī)儲層中的裂縫系統(tǒng)利用L-system進行模擬，該方法將微地震數(shù)據(jù)解釋為與之匹配的分形網(wǎng)絡(luò)，進而解釋復(fù)雜裂縫網(wǎng)絡(luò)的分叉及多尺度特征。注水誘導(dǎo)裂縫的形成過程近似“從樹干逐漸擴展到連接的樹枝”，因此也有學(xué)者嘗試將其應(yīng)用于注水誘導(dǎo)縫的模擬。Fan等[78]基于分形特征對分形誘導(dǎo)裂縫網(wǎng)絡(luò)進行模擬，并結(jié)合數(shù)值模擬進行了驗證。

圖12 基于分形特征迭代的裂縫模型(據(jù)文獻(xiàn)[72， 75， 76， 78]，有改動)(a)二維IFS模型； (b)三維IFS模型； (c)二維L-system模型； (d)三維分形誘導(dǎo)裂縫網(wǎng)絡(luò)示意圖

7 結(jié)論和討論

儲層裂縫分布的復(fù)雜性使得儲層裂縫建模成為一大難點，目前隨機裂縫建模方法是較常應(yīng)用的一類建模方法。本文系統(tǒng)地調(diào)研了現(xiàn)有的裂縫隨機建模方法，將其分為五大類，他們的優(yōu)缺點、適用性、發(fā)展趨勢概括如下。

(1)基于空間剖分的裂縫建模

基于空間剖分的裂縫建模原理簡單，便于操作，適用于高角度正交裂縫、柱狀節(jié)理、巖石破裂等，但不適用于復(fù)雜裂縫建模。無邊界正交模型和有邊界正交模型適合模式較為簡單的正交裂縫系統(tǒng)，如區(qū)域構(gòu)造裂縫。馬賽克鑲嵌模型簡單易行，可較好地應(yīng)用于玄武巖柱狀節(jié)理建模及巖石破碎模擬，但此類方法初始隨機點的確定準(zhǔn)則是一個難點。Veneziano模型包含大部分的空間剖分的裂縫模型，但其生成符合裂縫形態(tài)的多邊形的效率較低。基于空間剖分的裂縫建模近年來發(fā)展緩慢，但仍有應(yīng)用潛力。

(2)離散裂縫網(wǎng)絡(luò)建模

離散裂縫網(wǎng)絡(luò)模型中的Baecher模型是目前最常用的裂縫模型之一，簡單易行，便于后續(xù)數(shù)模。它可綜合多種信息將已有裂縫模式在研究區(qū)進行裂縫空間位置拓展，所建立的模型可包含較全的裂縫屬性(傾向、傾角、長度、開度等)，這也是其他模型所不具備的。該方法易于與地質(zhì)模式相結(jié)合，進而形成更加符合地質(zhì)認(rèn)識的裂縫模型。例如在斷層附近的裂縫系統(tǒng)建模時可以通過改進點過程，使其符合遠(yuǎn)離斷層裂縫密度低的規(guī)律；晶質(zhì)巖中的裂縫系統(tǒng)常呈現(xiàn)分形特征，將分形特性融入到點過程和示性過程，即可較好地體現(xiàn)此類裂縫系統(tǒng)的特性；裂縫密度約束點過程可以使生成裂縫網(wǎng)絡(luò)更加逼近此類裂縫系統(tǒng)；此外也可將研究區(qū)應(yīng)力狀態(tài)導(dǎo)致的裂縫模式融入。實際裂縫系統(tǒng)中裂縫展布常受層控，因為巖性變化對裂縫終止有較強影響，如沉積巖儲層中裂縫，DFN模型可以通過分層位、不同組系、不同規(guī)模裂縫的疊加來實現(xiàn)多尺度復(fù)雜裂縫系統(tǒng)的建模。

然而DFN模型在實際應(yīng)用中也存在三個主要問題： ①需要預(yù)設(shè)裂縫的形狀、大小、傾向、傾角等的分布及其分布中的參數(shù)。這類參數(shù)可以通過露頭分析、巖心觀察、測井解釋等獲取，但常存在采樣不均一、不充分、數(shù)據(jù)需要校正等一系列問題； ②裂縫系統(tǒng)中裂縫在空間上常具有相關(guān)性，DFN模型未考慮這一點，這是該模型的不足，而這一點正是基于變差函數(shù)和多點地質(zhì)統(tǒng)計學(xué)的裂縫建模的一大優(yōu)點。③不同組系裂縫數(shù)量比例如何分配，這一問題也是前三類方法都面臨的問題。

近年來，DFN建模研究趨勢主要有以下幾個方面：①裂縫中心位置確定方法的研究，如點過程的改進(非穩(wěn)態(tài)泊松過程、MCMC、Cluster過程、Cox過程等)、裂縫密度約束下的裂縫中心位置產(chǎn)生方法研究等；②示性過程中裂縫屬性的研究，如裂縫形狀的選擇、模型中裂縫參數(shù)(傾向、傾角、大小、開度等)的分布函數(shù)研究；③DFN模型中參數(shù)獲取問題的研究。獲取裂縫參數(shù)的手段越來越多，如通過遙感衛(wèi)星數(shù)據(jù)獲取大尺度裂縫及斷裂帶信息[79，80]；基于三維激光掃描(如Light Detection and Ranging， LiDAR)技術(shù)對露頭進行掃描獲取三維數(shù)據(jù)[81-83]，并通過手動或自動的方式提取裂縫參數(shù)[80，84]；基于掃描線法的露頭裂縫信息獲取[85]；利用無人駕駛飛機(Unmanned Aerial Vehicle，UAV)獲取多尺度裂縫模式、長度等信息[86]；利用成像測井獲取單井裂縫信息[87]，如位置、傾角、開度等；通過巖心觀察統(tǒng)計單井裂縫信息；利用CT掃描獲取小尺度裂縫信息[88，89]，如裂縫幾何形態(tài)等。這些裂縫規(guī)模數(shù)千米到幾納米，不可避免的是裂縫多尺度耦合的問題，如何合理地將其轉(zhuǎn)換成研究問題模型中的參數(shù)，即參數(shù)校正成為一項重點工作[90]；④地質(zhì)信息的融合，針對具體地質(zhì)模式對DFN的改進，如與褶皺伴生的裂縫系統(tǒng)建模[91，92]、與斷層伴生的裂縫系統(tǒng)建?？紤]斷距的影響因素[93，94]等；⑤構(gòu)建裂縫密度約束體的研究。建模時，裂縫參數(shù)的統(tǒng)計分布、露頭統(tǒng)計數(shù)據(jù)等可以反映裂縫的統(tǒng)計信息，但不能有效地指出裂縫更為發(fā)育的區(qū)域，而裂縫密度約束體正是起到這一作用，對裂縫中心位置的生成過程或裂縫方位確定過程進行約束。目前構(gòu)建約束體有：通過疊前三維地震數(shù)據(jù)預(yù)測裂縫密度及方位[95，96]、弱度[97，98]、扁率、通過疊后地震數(shù)據(jù)邊緣檢測獲取裂縫帶的信息[99-101]、通過應(yīng)力場數(shù)值模擬構(gòu)建裂縫密度約束體[94，102]，井震結(jié)合地質(zhì)建模等實現(xiàn)[93，103]；⑥裂縫模型有效性檢驗，即質(zhì)量控制(QC)。如對生成的裂縫模型進行抽樣，統(tǒng)計一維密度、裂縫長度和二維密度，并與實測裂縫進行對比[104]，將裂縫連通性分析[20，105，106]結(jié)果與示蹤劑數(shù)據(jù)進行對比，并與實際生產(chǎn)數(shù)據(jù)進行對比[93]。根據(jù)對比結(jié)果對模型參數(shù)等進行調(diào)整，校正模型。重復(fù)上述過程直至滿足質(zhì)量控制條件，否則模型不可用。

除上述基于示性點過程的離散裂縫網(wǎng)絡(luò)模型外，利用其他數(shù)學(xué)方法生成離散裂縫網(wǎng)絡(luò)的模型也逐漸得到重視，如Veneziano改進模型。

(3)基于變差函數(shù)的裂縫建模

基于變差函數(shù)的裂縫建模有兩類方法：一類是用沉積相建模流程直接對井上測井解釋的裂縫用地質(zhì)統(tǒng)計學(xué)的方法進行，但由于裂縫非均質(zhì)性遠(yuǎn)超過沉積相的變化，因此不能較好地模擬裂縫形態(tài)，同時對裂縫屬性也不能較好地估計； 另一類研究多在DFN模型的框架下開展，因此它繼承了DFN模型的大部分特點。它用指示變換和變差函數(shù)來模擬裂縫屬性，其功能近似于Baecher模型中示性過程，但其考慮了裂縫的空間相關(guān)性，是對DFN模型的改進，而兩點地質(zhì)統(tǒng)計學(xué)對復(fù)雜裂縫系統(tǒng)中裂縫的空間相關(guān)性表征不足是其面臨的一大問題。從發(fā)表成果看，近年來該方法在實用性方面有了較大改進。其研究趨勢除了DFN模型中提到的六點外，復(fù)雜裂縫系統(tǒng)空間相關(guān)性的描述也是亟待解決的問題。

(4)基于多點地質(zhì)統(tǒng)計學(xué)的裂縫建模

基于多點地質(zhì)統(tǒng)計學(xué)的裂縫建模，易于結(jié)合先驗裂縫模式信息，且忠實于硬數(shù)據(jù)，并在后續(xù)數(shù)值模擬方面具有較強優(yōu)勢，適用于靜態(tài)裂縫建模。但訓(xùn)練圖像獲取及算法效率提高是其發(fā)展過程中必須要解決的兩個問題。目前訓(xùn)練圖像多為二維圖像，常通過衛(wèi)星、露頭等獲得，如何獲取適合研究區(qū)裂縫建模的三維訓(xùn)練圖像是近年來逐漸被關(guān)注的一個研究方向，如通過三維地震解釋斷層、通過應(yīng)力場數(shù)值模擬三維斷層、通過三維CT獲取裂縫。獲取的訓(xùn)練圖像普遍存在尺度耦合的問題，這也是該類方法需要解決的問題。同時訓(xùn)練圖像穩(wěn)定性問題是應(yīng)用多點地質(zhì)統(tǒng)計學(xué)進行模擬時要考慮的一個問題。算法效率提高也是多點地質(zhì)統(tǒng)計學(xué)近年來的研究熱點，數(shù)據(jù)結(jié)構(gòu)、降維、并行計算等手段相繼被引入，雖有所緩解，但仍待改進。

(5)基于分形特征迭代的裂縫建模

基于分形特征迭代的裂縫建模簡單易行，既可用于靜態(tài)裂縫模擬，也可用于動態(tài)模擬(如巖石破裂的模擬、注水誘導(dǎo)縫模擬等)。該方法僅適用于分形特征明顯的裂縫系統(tǒng)模擬。如何高效地生成符合約束條件(如微地震數(shù)據(jù)、裂縫密度數(shù)據(jù))的分形模型是目前該方法的一個研究重點，如通過約束條件和優(yōu)化算法改進迭代規(guī)則生成符合地質(zhì)認(rèn)識的裂縫模型。

綜上所述，近年來裂縫隨機建模方法的研究趨勢可歸納為四點： ①原有數(shù)學(xué)算法的改進及新算法的引入，如空間剖分方法中Delaunay三角剖分算法的改進、離散裂縫網(wǎng)絡(luò)模型中隨機過程穩(wěn)定性的考慮、多點地質(zhì)統(tǒng)計學(xué)中掃描效率的提高、分形特征迭代模型中迭代算法的改進等； ②將地質(zhì)認(rèn)識融合到數(shù)學(xué)模型中，建立針對某一特定地質(zhì)體的裂縫系統(tǒng)，如與斷層伴生的裂縫系統(tǒng)建模、具有明顯分形特征的晶質(zhì)巖裂縫等； ③確定模型參數(shù)的相關(guān)研究。首先是資料獲取技術(shù)的改進，如疊前三維地震、成像測井、CT掃描、三維激光掃描、無人駕駛飛機等手段的引入。其次是所獲取的不同尺度裂縫信息的耦合，將所獲取的裂縫信息轉(zhuǎn)換成模型中的參數(shù)； ④裂縫模型有效性驗證，即質(zhì)量控制。

感謝毛哲、王小垚、曹東升等在本文成文過程中提供了幫助！

參考文獻(xiàn)

[1] 曾聯(lián)波．低滲透砂巖儲層裂縫的形成與分布．北京：科學(xué)出版社，2008．

Zeng Lianbo．Formation and Distribution of Fracture in Low-permeability Reservoir．Science Press，Beijing，2008.

[2] 邱殿明．?dāng)嗔?、斷層、?jié)理、劈理、裂隙、裂縫之間的關(guān)系小結(jié)．吉林大學(xué)學(xué)報(地球科學(xué)版)，2013，43(5)：1392.

[3] 潘欣．構(gòu)造裂縫識別與建模研究——以鄂爾多斯盆地中西部定邊—華池地區(qū)延長組為例[學(xué)位論文]．陜西西安：西北大學(xué)，2010．

Pan Xin．Identifying and Modeling Research of Structural Fracture — A Case Study of the Yanchang Formation in Dingbian-Huachi Region of the Central-West Ordos Basin [D]．Northwest University，Xi’an，Shaanxi，2010.

[4] Ivanova V M．Three-dimensional Stochastic Modeling of Rock Fracture Systems [D]．Massachusetts Institute of Technology， 1995.

[5] Childs E C．The anisotropic hydraulic conductivity of soil．European Journal of Soil Science，1957，8(1)：42-47.

[6] Robinson N I，Sharp Jr J M，Kreisel I．Contaminant transport in sets of parallel finite fractures with fracture skins．Journal of Contaminant Hydrology，1998，31(1-2)：83-109.

[7] Dunne W M，North C P．Orthogonal fracture systems at the limits of thrusting：an example from southwestern Wales．Journal of Structural Geology，1990，12(2)：207-215.

[8] Dershowitz W S．Rock Joint Systems [D]．Massachusetts Institute of Technology，1984.

[9] Dershowitz W S，Einstein H H．Characterizing rock joint geometry with joint system models．Rock Mechanics and Rock Engineering，1988，21(1)：21-51.

[10] Bai T，Maerten L，Gross M R et al．Orthogonal cross joints：do they imply a regional stress rotation?Journal of Structural Geology，2002，24(1)：77-88.

[11] Alghalandis Y F．Stochastic Modelling of Fractures in Rock Masses [D]．University of Adelaide，2014.

[12] Baecher G B，Lanney N A，Einstein H H．Statistical description of rock properties and sampling．Proceedings of the 18th US Symposium on Rock Mechanics,Golden,CO,1977.

[13] Ivanova V M，Sousa R，Murrihy B et al．Mathematical algorithm development and parametric studies with the GEOFRAC three-dimensional stochastic model of natural rock fracture systems．Computers & Geosciences，2014，67：100-109.

[14] Mosser L J，Matth?i S K．Tessellations stable under iteration．International Discrete Fracture Network Engineering Conference，2014.

[15] Baecher G B．Statistical analysis of rock mass fracturing．Journal of the International Association for Mathematical Geology，1983，15(2)：329-348.

[16] Dowd P A，Xu C，Mardia K V et al．A comparison of methods for the stochastic simulation of rock fractures．Mathematical Geology，2007，39(7)：697-714.

[17] Xu C，Dowd P．A new computer code for discrete fracture network modelling．Computers & Geosciences，2010，36(3)：292-301.

[18] Nakaya S，Nakamura K．Percolation conditions in fractured hard rocks：A numerical approach using the three-dimensional binary fractal fracture network (3D-BFFN) model．Journal of Geophysical Research，2007，112(B12)：1-15.

[19] Alghalandis F Y．ADFNE：Open source software for discrete fracture network engineering，two and three dimensional applications．Computers & Geosciences，2017，102：1-11.

[20] Dong S，Zeng L，Dowd P et al．A fast method for fracture intersection detection in discrete fracture networks．Computers and Geotechnics，2018，98：205-216.

[21] Mardia K V，Nyirongo V B，Walder A N et al．Markov chain Monte Carlo implementation of rock fracture modelling．Mathematical Geology，2007，39(4)：355-381.

[22] Nakaya S，Yoshida T，Shioiri N．Percolation conditions in binary fractal fracture networks：Applications to rock fractures and active and seismogenic faults．Journal of Geophysical Research：Solid Earth，2003，108(B7)：1-13.

[23] Black J H．Hydrogeology of fractured rocks-a question of uncertainty about geometry．Applied Hydrogeology，1994，2(3)：56-70.

[24] Geier J，Lee K，Dershowitz W S．Field validation of conceptual models for fracture geometry．Transactions of American Geophysical Union，1988，69(44)：1177.

[25] Chiles J P．Fractal and geostatistical methods for modeling of a fracture network．Mathematical Geology，1988，20(6)：631-654.

[26] Billaux D，F(xiàn)uller P．An algorithm for mesh simplification applied to fracture hydrology．Mathematical Geology，1989，21(2)：221-232.

[27] Kulatilake P H S W，Um J，Wang M et al．Stochastic fracture geometry modeling in 3-D including validations for a part of Arrowhead East Tunnel，California，USA．Engineering Geology，2003，70(1)：131-155.

[28] Zadhesh J，Jalali S-M E，Ramezanzadeh A．Estimation of joint trace length probability distribution function in igneous， sedimentary，and metamorphic rocks．Arabian Journal of Geosciences，2014，7(6)：2353-2361.

[29] Noroozi M，Kakaie R，Jalali S E．3D geometrical-stochastical modeling of rock mass joint networks：case study of the right bank of rudbar Lorestan dam plant．Journal of Geology and Mining Research，2015，7(1)：1-10.

[30] Watanabe N，Ishibashi T，Hirano N et al．Precise 3D numerical modeling of fracture flow coupled with X-ray computed tomography for reservoir core samples．SPE Journal，2011，16(3)：683-691.

[31] Gong J，Rossen W R．Modeling flow in naturally fractured reservoirs：Effect of fracture aperture distribution on dominant sub-network for flow．Petroleum Science，2017，14(1)：138-154.

[32] Ivanova V M．Geologic and Stochastic Modeling of Fracture Systems in Rocks [D]．Massachusetts Institute of Technology，1998.

[33] 鄭松青，姚志良．離散裂縫網(wǎng)絡(luò)隨機建模方法．石油天然氣學(xué)報，2009，31(4)：106-110．

Zheng Qingsong，Yao Zhiliang．Discrete fracture network stochastic modeling．Journal of Oil and Gas Technology，2009，31(4)：106-110.

[34] 于青春，大西有三．巖體三維不連續(xù)裂隙網(wǎng)絡(luò)及其逆建模方法．地球科學(xué)，2003，28(5)：522-526．

Yu Qingchun，Ohnishi Yuzo．Three-dimensional discrete fracture network model and its inverse method．Earth Science—Journal of China University of Geosciences，2003，28(5)：522-526.

[35] 宋曉晨，徐衛(wèi)亞．裂隙巖體滲流模擬的三維離散裂隙網(wǎng)絡(luò)數(shù)值模型(Ⅰ)：裂隙網(wǎng)絡(luò)的隨機生成．巖石力學(xué)與工程學(xué)報，2004，23(12)：2015-2020．

Song Xiaochen，Xu Weiya．Numerical model of three-dimensional discrete fracture network for seepage in fractured rocks (Ⅰ)：generation of fracture network．Chinese Journal of Rock Mechanics and Engineering，2004，23(12)：2015-2020.

[36] Wettstein S J，Wittel F K，Araújo N A M et al．From invasion percolation to flow in rock fracture networks．Physica A：Statistical Mechanics and its Applications，2012，391(1)：264-277.

[37] 高志勇，朱如凱，馮佳睿等．中國前陸盆地構(gòu)造—沉積填充響應(yīng)與深層儲層特征．北京：地質(zhì)出版社，2016．

Gao Zhiyong，Zhu Rukai，F(xiàn)eng Jiarui et al．Tectonic-sedimentary Filling Response and Deep Reservoir Characteristics of Foreland Basins in China．Geological Publishing House，Beijing，2016.

[38] Surrette M J，Allen D M．Quantifying heterogeneity in variably fractured sedimentary rock using a hydrostructural domain．Geological Society of America Bulletin，2008，120(1-2)：225-237.

[39] Lyu W，Zeng L，Liao Z et al．Fault damage zone cha-racterization in tight-oil sandstones of the Upper Triassic Yanchang Formation in the southwest Ordos Basin， China： Integrating cores， image logs， and conventional logs．Interpretation，2017，5(4)：SP27-SP39.

[40] Lyu W，Zeng L，Liu Z et al．Fracture responses of conventional logs in tight-oil sandstones：A case study of the Upper Triassic Yanchang Formation in southwest Ordos Basin，China．AAPG Bulletin，2016，100(9)：1399-1417.

[41] Tóth T M．Determination of geometric parameters of fracture networks using 1D data．Journal of Structural Geology，2010，32(7)：878-885.

[42] Davy P，Sornette A，Sornette D．Some consequences of a proposed fractal nature of continental faulting．Nature，1990，348(6296)：56-58.

[43] Bour O，Davy P．Connectivity of random fault networks following a power law fault length distribution．Water Resources Research，1997，33(7)：1567-1583.

[44] Bour O，Davy P．Clustering and size distributions of fault patterns：Theory and measurements．Geophysical Research Letters，1999，26(3)：2001-2004.

[45] Magnusdottir L，Horne R N．Inversion of time-lapse electric potential data to estimate fracture connectivity in geothermal reservoirs．Mathematical Geosciences，2015，47(1)：85-104.

[46] 吳勝和．儲層表征與建模．北京：石油工業(yè)出版社，2010．

Wu Shenghe．Reservoir Description and Modeling．Petroleum Industry Press，Beijing，2010.

[47] Xu C，Dowd P A，Mardia K V et al．A flexible true plurigaussian code for spatial facies simulations．Computers & Geosciences，2006，32(10)：1629-1645.

[48] Koike K，Kubo T，Liu C et al．3D geostatistical modeling of fracture system in a granitic massif to characterize hydraulic properties and fracture distribution．Tectonophysics，2015，660：1-16.

[49] Koike K，Liu C，Sanga T．Incorporation of fracture directions into 3D geostatistical methods for a rock fracture system．Environmental Earth Sciences，2012，66(5)：1403-1414.

[50] Sebacher B，Hanea R，Stordal A S．An adaptive pluri-Gaussian simulation model for geological uncertainty quantification．Journal of Petroleum Science and Engineering，2017，158：494-508.

[51] Koike K，Sakamoto H，Ohmi M．Detection and hydrologic modeling of aquifers in unconsolidated alluvial plains through combination of borehole data sets：A case study of the Arao area，Southwest Japan．Engineering Geology，2001，62(4)：301-317.

[52] Masoud A，Koike K．Tectonic architecture through Landsat-7 ETM+/SRTM DEM-derived lineaments and relationship to the hydrogeologic setting in Siwa region，NW Egypt．Journal of African Earth Sciences，2006，45(4-5)：467-477.

[53] Deutsch C V，Journel A G．Gslb：Geostatistical Software Library and User’s Guide．Oxford University Press，1992.

[54] Guardiano F B，Srivastava R M．Multivariate Geostatistics：Beyond Bivariate Moment．Springer，1993.

[55] Strebelle S．Conditional simulation of complex geolo-gical structures using multiple-point statistics．Mathematical Geology，2002，34(1)：1-21.

[56] 張文彪，段太忠，劉彥鋒等．綜合沉積正演與多點地質(zhì)統(tǒng)計模擬碳酸鹽巖臺地——以巴西Jupiter油田為例．石油學(xué)報，2017，38(8)：925-934．

Zhang Wenbiao，Duan Taizhong，Liu Yanfeng et al．Integrated sedimentary forward modeling and multipoint geostatistics in carbonate platform simulation：A case study of Jupiter oilfield in Brazil．Acta Petrolei Sinica，2017，38(8)：925-934.

[57] Arpat G B．Sequential Simulation with Patterns [D]．Stanford University，2005.

[58] Arpat G B，Caers J．Conditional simulation with patterns．Mathematical Geology，2007，39(2)：177-203.

[59] Yin Y，Qin Z，Zuo R．Simulation of reservoir geometry using simpat method．Acta Geologica Sinica，2013，87(Z1)：607-609.

[60] Zhang T，Switzer P，Journel A．Filter-based classification of training image patterns for spatial simulation．Mathematical Geology，2006，38(1)：63-80.

[61] Honarkhah M，Caers J．Stochastic simulation of pat-terns using distance-based pattern modeling．Mathematical Geosciences，2010，42(5)：487-517.

[62] 喻思羽，李少華，何幼斌等．基于樣式降維聚類的多點地質(zhì)統(tǒng)計建模算法．石油學(xué)報，2016，37(11)：1403-1409．

Yu Siyu，Li Shaohua，He Youbin et al．Multiple-point geostatistics algorithm based on pattern scale-down cluster．Acta Petrolei Sinica，2016，37(11)：1403-1409.

[63] 謝青．基于GPU的頁巖儲層裂縫建模及壓力模擬[學(xué)位論文]．安徽合肥：中國科學(xué)技術(shù)大學(xué)，2014．

Xie Qing．GPU based Fracture Modeling and Pressure Distribution Simulation of Shale Reservoirs [D]．University of Science and Technology of China，Hefei，Anhui，2014.

[64] Mohammadmoradi P．Facies and fracture network modeling by a novel image processing based method．Geomaterials，2013，3(4)：156-164.

[65] 張烈輝，賈鳴，張芮菡等．裂縫性油藏離散裂縫網(wǎng)絡(luò)模型與數(shù)值模擬．西南石油大學(xué)學(xué)報(自然科學(xué)版)，2017，39(3)：121-127．

Zhang Liehui，Jia Ming，Zhang Ruihan et al．Discrete fracture network modeling and numerical simulation of fractured reservoirs．Journal of Southwest Petroleum University(Science & Technology Edition)，2017，39(3)：121-127.

[66] Jia M，Zhang L，Guo J．Combining a connected-component labeling algorithm with FILTERSIM to simulate continuous discrete fracture networks．Environmental Earth Sciences，2017，76(8)：327.

[67] Liu X，Zhang C，Liu Q et al．Multiple-point statistical prediction on fracture networks at Yucca Mountain．Environmental Geology，2009，57(6)：1361-1370.

[68] Chugunova T L，Hu L Y．Multiple-point simulations constrained by continuous auxiliary data．Mathematical Geosciences， 2008，40(2)：133-146.

[69] Chugunova T，Corpel V，Gomez J P．Explicit fracture network modelling：From multiple point statistics to dynamic simulation．Mathematical Geosciences，2017，49(4)：541-553.

[70] Qiao J．Julia sets and complex singularities of free energies．Memoirs of the American Mathematical Society，2015，5：1-3.

[71] Qiao J．Julia sets and complex singularities in diamond-like hierarchical Potts models．Science China Mathematics，2005，48(3)：388-412.

[72] Acuna J A，Yortsos Y C．Application of fractal geome-try to the study of networks of fractures and their pressure transient．Water Resources Research，1995，31(3)：527-540.

[73] 袁杰，李霞麗，白水．迭代函數(shù)系統(tǒng)及其應(yīng)用．中央民族大學(xué)學(xué)報(自然科學(xué)版)，2004，13(4)：317-321．

Yuan Jie，Li Xiali，Bai Shui．Iterated function system and its application．Journal of The Central University For Nationalities (Natural Science Edition)，2004，13(4)：317-321.

[74] 周德華，焦方正，葛家理．一種油氣藏裂縫網(wǎng)絡(luò)模擬優(yōu)化生成技術(shù)．?dāng)鄩K油氣田，2002，9(4)：30-32．

Zhou Dehua，Jiao Fangzheng，Ge Jiali．A new approach to optimistically model fracture networks of hydrocarbon reservoirs．Fault-block Oil & Gas Field，2002，9(4)：30-32.

[75] Acuna J A，Yortsos Y C．Numerical construction and flow simulation in networks of fractures using fractal geometry．Society of Petroleum Engineers，1991:1-7.

[76] Zhou Z，Su Y，Wang W et al．Application of the fractal geometry theory on fracture network simulation．Journal of Petroleum Exploration & Production Technology，2017，7(2)：487-496.

[77] Zhou Z，Su Y，Wang W et al．Integration of microseismic and well production data for fracture network calibration with an L-system and rate transient analysis．Journal of Unconventional Oil & Gas Resources，2016，15：113-121.

[78] Fan D，Ettehadtavakkol A．Semi-analytical modeling of shale gas flow through fractal induced fracture networks with microseismic data．Fuel，2017，193：444-459.

[79] Dasho O A，Ariyibi E A，Akinluyi F O et al． Application of satellite remote sensing to groundwater potential modeling in Ejigbo area，Southwestern Nigeria．Modeling Earth Systems and Environment，2017，3(2)：615-633.

[80] Cao T，Xiao A，Wu L et al．Automatic fracture detection based on Terrestrial Laser Scanning data：A new method and case study．Computers & Geosciences，2017，106：209-216.

[81] 曾慶魯，張榮虎，盧文忠等．基于三維激光掃描技術(shù)的裂縫發(fā)育規(guī)律和控制因素研究——以塔里木盆地庫車前陸區(qū)索罕村露頭剖面為例．天然氣地球科學(xué)，2017，28(3)：397-409．

Zeng Qinglu，Zhang Ronghu，Lu Wenzhong et al．Fracture development characteristics and controlling factors based on 3D laser scanning technology：An outcrop case study of Suohan village，Kuqa foreland area，Tarim Basin．Natural Gas Geoscience，2017，28(3)：397-409.

[82] Vazaios I，Vlachopoulos N，Diederichs M S．Integra-tion of lidar-based structural input and discrete fracture network generation for underground applications．Geotechnical and Geological Engineering，2017，35(5)：2227-2251.

[83] Casini G，Hunt D W，Monsen E et al．Fracture characterization and modeling from virtual outcrops．AAPG Bulletin，2016，100(1)：41-61.

[84] Wang X，Zou L，Shen X et al．A region-growing approach for automatic outcrop fracture extraction from a three-dimensional point cloud．Computers & Geosciences，2017，99：100-106.

[85] Santos R F V C，Miranda T S，Barbosa J A et al．Chara-cterization of natural fracture systems：Analysis of uncertainty effects in linear scanline results．AAPG Bulletin，2015，99(12)：2203-2219.

[86] Bisdom K，Nick H M，Bertotti G．An integrated workflow for stress and flow modelling using outcrop-derived discrete fracture networks．Computers & Geosciences，2017，103：21-35.

[87] 肖小玲，靳秀菊，張翔等．基于常規(guī)測井與電成像測井多信息融合的裂縫識別．石油地球物理勘探，2015，50(3)：542-547．

Xiao Xiaoling，Jin Xiuju，Zhang Xiang et al．Fracture identification based on information fusion of conventional logging and electrical imaging logging．OGP，2015，50(3)：542-547.

[88] Lai J，Wang G，F(xiàn)an Z et al．Three-dimensional quantitative fracture analysis of tight gas sandstones using industrial computed tomography．Scientific Reports，2017，7(1)：1825.

[89] Karimpouli S，Tahmasebi P，Ramandi H L et al．Stochastic modeling of coal fracture network by direct use of micro-computed tomography images．International Journal of Coal Geology，2017，179：153-163.

[90] Huang L，Tang H，Tan Q et al．A novel method for correcting scanline-observational bias of discontinuity orientation．Scientific Reports，2016，6(1)：22942.

[91] Sanai L，Chenini I，Ben M A et al．Fracture related-fold patterns analysis and hydrogeological implications：Insight from fault-propagation fold in Northwestern of Tunisia．Journal of African Earth Sciences，2015，101：375-382.

[92] 祖克威，曾聯(lián)波，鞏磊．?dāng)鄬酉嚓P(guān)褶皺概念模型中的裂縫域．地質(zhì)科學(xué)，2013，48(4)：1140-1147．

Zu Kewei,Zeng Lianbo,Gong Lei.Fractures domains in conceptual models of fault-related folds．Chinese Journal of Geology，2013，48(4)：1140-1147.

[93] 郎曉玲，郭召杰．基于DFN離散裂縫網(wǎng)絡(luò)模型的裂縫性儲層建模方法．北京大學(xué)學(xué)報(自然科學(xué)版)，2013，？(6)：964-972．

Lang Xiaoling，Guo Zhaojie．Fractured reservoir modeling method based on discrete fracture network model．Acta Scientiarum Naturalium Universitatis Pekinensis，2013，？(6)：964-972.

[94] 張雨晴，王暉，范廷恩等．花崗巖潛山儲層裂縫建模表征方法——以渤?；◢弾r潛山A油田為例．吉林大學(xué)學(xué)報(地)，2016，46(5)：1312-1320．

Zhang Yuqing，Wang Hui，F(xiàn)an Ting’en et al．Granite buried hill reservoir characterization and modeling：Taking offshore an oilfield in Bohai Bay as an example．Journal of Jilin University (Earth Science Edition)，2016,46(5)：1312-1320.

[95] 熊曉軍，簡世凱，李翔等．基于標(biāo)準(zhǔn)差統(tǒng)計的窄方位角疊前裂縫預(yù)測方法．石油地球物理勘探，2017，52(1)：114-120．

Xiong Xiaojun，Jian Shikai，Li Xiang et al．Fracture prediction on prestack narrow azimuth data with standard deviation statistic.OGP,2017,52(1):114-120.

[96] 李慧瓊，張盟勃，蒲仁海等．黃257井區(qū)疊前縱波方位各向異性裂縫分布預(yù)測．石油地球物理勘探，2017，52(2)：350-359．

Li Huiqiong，Zhang Mengbo，Pu Renhai et al．Late Triassic fracture detection with seismic azimuth anisotropics in Huang257 survey，OrdosBasin．OGP，2017，52(2)：350-359.

[97] Xue J，Gu H，Cai C．Model-based amplitude versus offset and azimuth inversion for estimating fracture parameters and fluid content．Geophysics，2017，82(2)：M1-M17.

[98] 薛姣，顧漢明，蔡成國．準(zhǔn)靜態(tài)孔縫介質(zhì)廣義裂隙弱度研究．石油地球物理勘探，2015，50(6)：1146-1153．

Xue Jiao，Gu Hanming，Cai Chengguo．General fracture weaknesses for quasi-static porous fractured media．OGP，2015，50(6)：1146-1153.

[99] 袁曉宇，李映濤，葉寧等．基于頻譜分解的螞蟻追蹤裂縫檢測技術(shù)在玉北地區(qū)的應(yīng)用研究．石油地球物理勘探，2015，50(4)：665-671．

Yuan Xiaoyu，Li Yingtao，Ye Ning et al．The application of ant-tracking fracture detection based on spectrum decomposition in Yubei area．OGP，2015，50(4)：665-671.

[100] Wilson T H，Smith V，Brown A．Developing a model discrete fracture network，drilling，and enhanced oil recovery strategy in an unconventional naturally fractured reservoir using integrated field，image log，and three-dimensional seismic data．AAPG Bulletin，2015，99(4)：735-762.

[101] 汪忠德，薛詩桂，李紅敬等．油氣地球物理解釋技術(shù)研究新進展．地球物理學(xué)進展，2016，31(3)：1187-1201．

Wang Zhongde，Xue Shigui，Li Hongjing et al．New progress in the technology of oil and gas geophysical interpretation．Progress in Geophysics，2016，31(3)：1187-1201.

[102] 劉建偉，張云銀，曾聯(lián)波等．非常規(guī)油藏地應(yīng)力和應(yīng)力甜點地球物理預(yù)測——渤南地區(qū)沙三下亞段頁巖油藏勘探實例．石油地球物理勘探，2016，51(4)：792-800．

Liu Jianwei，Zhang Yunyin，Zeng Lianbo et al．Geophysical prediction of stress and stress desserts in unconventional reservoirs：An example in Bonan area．OGP，2016，51(4)：792-800.

[103] 薛艷梅，夏東領(lǐng)，蘇宗富等．多信息融合分級裂縫建模．西南石油大學(xué)學(xué)報(自然科學(xué)版)，2014，36(2)：57-63．

Xue Yanmei，Xia Dongling，Su Zongfu et al．Fracture modeling at different scales based on convergent multi-source information．Journal of Southwest Petroleum University (Science & Technology Edition)，2014，36(2)：57-63.

[104] 張成，夏露，李英杰等．柴達(dá)木盆地東部石炭系頁巖三維裂縫網(wǎng)絡(luò)模型研究．地學(xué)前緣，2016，23(5)：184-192．

Zhang Cheng，Xia Lu，Li Yingjie et al．Study on three-dimensional fracture network models of Carboniferous shale in Eastern Qaidam Basin．Earth Science Frontiers，2016，23(5)：184-192.

[105] Xu C，Dowd P A，Mardia K V et al．A connectivity index for discrete fracture networks．Mathematical Geology，2007，38(5)：611-634.

[106] Alghalandis Y F，Dowd P A，Xu C．Connectivity field：a measure for characterising fracture networks．Mathematical Geosciences，2015，47(1)：63-83．