基于體視學(xué)技術(shù)的二維數(shù)據(jù)模擬三維裂隙網(wǎng)絡(luò)

倪春中，劉春學(xué)，張世濤，李?lèi)?ài)敏，李雨健，蒼寶海，劉 雷

(1.昆明理工大學(xué) 國(guó)土資源工程學(xué)院，云南 昆明 650093; 2.云南財(cái)經(jīng)大學(xué) 城市與環(huán)境學(xué)院，云南 昆明 650221;3.云南錫業(yè)集團(tuán)有限責(zé)任公司，云南 個(gè)舊 661021; 4.中國(guó)石油 遼河油田分公司 特種油開(kāi)發(fā)公司，遼寧 盤(pán)錦 124010;5.中國(guó)地震局 地震預(yù)測(cè)研究所，北京 100038)

基于體視學(xué)技術(shù)的二維數(shù)據(jù)模擬三維裂隙網(wǎng)絡(luò)

倪春中1，劉春學(xué)2，張世濤1，李?lèi)?ài)敏3，李雨健1，蒼寶海4，劉 雷5

(1.昆明理工大學(xué) 國(guó)土資源工程學(xué)院，云南 昆明 650093; 2.云南財(cái)經(jīng)大學(xué) 城市與環(huán)境學(xué)院，云南 昆明 650221;3.云南錫業(yè)集團(tuán)有限責(zé)任公司，云南 個(gè)舊 661021; 4.中國(guó)石油 遼河油田分公司 特種油開(kāi)發(fā)公司，遼寧 盤(pán)錦 124010;5.中國(guó)地震局 地震預(yù)測(cè)研究所，北京 100038)

CT掃描；體視學(xué)；裂隙網(wǎng)絡(luò)；巖石裂隙；個(gè)舊錫礦

近年來(lái)國(guó)內(nèi)外對(duì)于裂隙空間分布的研究雖然很多，可所用的數(shù)據(jù)大多是通過(guò)鉆孔或者露頭等獲取的一維或二維樣品資料，而真正的三維裂隙網(wǎng)絡(luò)分布樣品資料是很難觀測(cè)的，因此不少學(xué)者嘗試通過(guò)統(tǒng)計(jì)學(xué)和應(yīng)用概率等方法將一維或二維的數(shù)據(jù)推算三維巖石裂隙網(wǎng)絡(luò)的屬性。地下巖石裂隙的勘察方法主要有鉆孔巖心法、超聲成像法、鉆孔電視法以及鉆孔流量測(cè)井法等[1]，小尺度的觀測(cè)方有如核磁共振(NMR)、X-射線掃描(CT)、γ-射線掃描和同步加速器(synchrotron)等。裂隙空間分布具有很多屬性，各屬性之間大部分不是相互獨(dú)立的，存在多種形式的相關(guān)關(guān)系，因此確定裂隙的長(zhǎng)度、密度等重要的屬性后[2-3]，其他變量可以隨之確定。Berkowitz和Adler發(fā)現(xiàn)裂隙網(wǎng)絡(luò)與觀測(cè)截面交嵌的圓盤(pán)直徑以及其他變量均服從Power Law分布、lognormal Law分布和Exponential Law分布[4]。Adler和Thovert運(yùn)用二維圖像中的痕跡長(zhǎng)度及數(shù)量估算三維裂隙網(wǎng)絡(luò)中裂隙密度和概率密度等屬性[5]。Yeomans認(rèn)為裂隙的個(gè)別屬性(如長(zhǎng)度、位移、開(kāi)口寬度等)以及整個(gè)裂隙網(wǎng)絡(luò)分布也服從Gamma Law趨勢(shì)[6]，而Dershowitz和Einstein則認(rèn)為裂隙方向可以用Bivariate Normal分布和Fisher分布進(jìn)行描述[7]。

本研究以云南個(gè)舊錫礦區(qū)高松礦田為例，所用的小尺度裂隙數(shù)據(jù)取自樣本巖石的CT掃描圖像，應(yīng)用體視學(xué)技術(shù)的方法將二維圖像中的裂隙延展及其連接等數(shù)據(jù)反演到三維空間，然后在三維空間相關(guān)特征的約束條件下，反演推斷出三維空間中裂隙密度分布，最終實(shí)現(xiàn)應(yīng)用二維裂隙數(shù)據(jù)模擬三維裂隙網(wǎng)絡(luò)。

1 地質(zhì)概況

高松礦田是云錫公司的主要生產(chǎn)礦山之一，位于個(gè)舊錫礦東礦區(qū)北部，處于礦區(qū)的五子山復(fù)背斜北段，

夾持于南北向個(gè)舊斷裂、甲界山斷裂與東西向個(gè)松斷裂、背陰山斷裂之間，地質(zhì)條件復(fù)雜。高松礦區(qū)主要出露地層為二疊系和三疊系，以白云巖、白云質(zhì)灰?guī)r為主，另有古近系與新近系泥巖、第四系殘坡積物零星出露。

褶皺構(gòu)造有近EW向?qū)﹂T(mén)山-阿西寨向斜、近SN向五子山背斜等。斷裂構(gòu)造可以分為EW，NE，NW，SN四組，其中以近EW向和NE向?yàn)橹?，NW向次之。EW向斷裂組自北而南有個(gè)松斷裂、麒麟山斷裂、馬吃水?dāng)嗔?、高阿斷裂、背陰山斷裂，呈近等距離分布，其中個(gè)松和背陰山斷裂規(guī)模較大，麒麟山斷裂次之。NE向斷裂組自西而東有蓮花山斷裂、蘆塘壩斷裂、麒阿西斷裂。NW向斷裂組主要為大箐東斷裂、黑碼石斷裂、駝峰山斷裂、阿西寨斷裂。SN向斷裂組不發(fā)育(圖1)。

2 樣本裂隙提取

2.1 巖石樣本

本研究在個(gè)舊錫礦高松礦田1 360，1 540，1 720和1 920 m四個(gè)不同中段采取巖石樣品35件，從中選取了12件巖石樣品進(jìn)行CT掃描，用于微觀尺度下三維裂隙網(wǎng)絡(luò)的模擬，各樣品采樣平面位置如圖1所示，在中段的位置分布如表1。

圖1 個(gè)舊礦區(qū)高松礦田構(gòu)造地質(zhì)略圖(文獻(xiàn)[8]，修改)Fig.1 Tectonic geological sketch of Gaosong mine in GejiuT2g1.中三疊統(tǒng)個(gè)舊組卡房段白云巖灰?guī)r；T2g2.中三疊統(tǒng)個(gè)舊組馬拉格段白云巖；T2g3.中三疊統(tǒng)個(gè)舊組白泥硐段灰?guī)r

樣品編號(hào)001002003004005006007008009010011012采樣位置/m192019201720172015401540154015401360136013601360

2.2 裂隙跡線提取

CT圖像中提取的信息主要包括：掃描定位的上下、前后、左右以及圖像中特殊點(diǎn)的位置，掃描定位后，CT機(jī)給出設(shè)定的圖像放大倍數(shù)，圖像中以X/Y的數(shù)值給出光標(biāo)指示點(diǎn)的二維坐標(biāo)；試樣的幾何特征，包括試驗(yàn)中試樣直徑、長(zhǎng)度、兩點(diǎn)間的距離和角度[9-10]，進(jìn)行幾何測(cè)量和內(nèi)部變形的研究；圖像的分布特征，從CT數(shù)及其方差正負(fù)變化的組合情況，可以判斷試件整體或某個(gè)部位處于壓密還是裂隙擴(kuò)展等；巖土試樣中某些節(jié)理構(gòu)造，在CT圖像中表現(xiàn)為某些方向的明暗條紋。

本次試驗(yàn)按照巖石的實(shí)際產(chǎn)狀對(duì)巖樣進(jìn)行CT掃描，得到了巖樣的三維空間網(wǎng)格點(diǎn)上的CT數(shù)，用Vo-xelCalc NDT Software進(jìn)行了初步數(shù)據(jù)處理，得到[xyzCTN]格式的元數(shù)據(jù)。然后用Matlab編寫(xiě)了數(shù)據(jù)處理程序，對(duì)CT掃描得到的元數(shù)據(jù)進(jìn)行了進(jìn)一步處理。通過(guò)調(diào)整各個(gè)巖樣CT數(shù)的閾值，使CT數(shù)在空間的變化更加明晰，便于觀察裂隙的展布(圖2)。

從50個(gè)水平掃描剖面中提取了相關(guān)的裂隙數(shù)據(jù)。裂隙數(shù)據(jù)除了實(shí)際觀測(cè)的跡線外，另外一部分資料利用Matlab編寫(xiě)的輔助程序從CT掃描圖像。巖樣水平剖面的大小多數(shù)為60 mm×80 mm。從中提取了284條裂隙跡線，其中最長(zhǎng)為64.1 mm，最短為2.2 mm，平均值為16.69 mm，方差為12.33 mm，集中在10 mm左右。

3 基于體視學(xué)方法由二維跡線模擬裂隙長(zhǎng)度及密度

體視學(xué)利用形態(tài)學(xué)和數(shù)學(xué)的分析，借助計(jì)算機(jī)及數(shù)據(jù)處理系統(tǒng)和顯微鏡及顯做成像系統(tǒng)，將二維平面

經(jīng)過(guò)成像及計(jì)算機(jī)分析處理得到三維形態(tài)，準(zhǔn)確地對(duì)物體進(jìn)行定量及形態(tài)結(jié)構(gòu)分析。本文應(yīng)用體視學(xué)技術(shù)，研究圓盤(pán)裂隙三維空間與二維剖面之間的相互聯(lián)系，解析裂隙長(zhǎng)度在不同維數(shù)之間的聯(lián)系和轉(zhuǎn)化，以處理二維剖面中的樣本裂隙跡線數(shù)據(jù)。

定義直徑為φ的圓盤(pán)，其中心點(diǎn)與觀測(cè)面之間的距離為z，其單位法線矢量與平面法線之間的夾角為α(圖3)。

假設(shè)φ為常數(shù)，直徑為φ的圓盤(pán)與觀測(cè)面交切的比例為S(z,φ)。從圖3可以看出，只有當(dāng)|z|<φ/2時(shí)才會(huì)有交切。對(duì)于給定的z，S(z,φ)可以表示為

(1)

式中：z為裂隙圓盤(pán)中心到觀測(cè)面S的距離，mm；φ為圓盤(pán)直徑，mm。

下述公式應(yīng)用冪次定律可以計(jì)算裂隙直徑分布[11-13]：

(2)

式中：α為分形維數(shù)，對(duì)于裂隙跡線長(zhǎng)度來(lái)說(shuō)，1≤α≤3；a為待定參數(shù)，由式(3)確定。

如果裂隙直徑長(zhǎng)度假定在φm≤φ≤φM的范圍內(nèi)，那么常數(shù)α可以表示為：

(3)

(4)

式中：φM,φm分別為裂隙最大直徑和最小直徑，mm；α為分形維數(shù)，根據(jù)裂隙分布試驗(yàn)確定。

根據(jù)直徑的冪次定律，假定裂隙直徑最大值φM是無(wú)限的，跡線c的概率可表示為：

圖2 巖樣CT圖像裂紋提取Fig.2 Crack extraction from rock CT image

圖3 利用體視學(xué)技術(shù)研究圓盤(pán)裂隙三維空間與二維剖面之間相互關(guān)系原理Fig.3 Study on the relationship between 3D space and 2D profiles of disc fissures used Stereologya.直徑為φ的裂紋圓盤(pán)中心點(diǎn)與觀測(cè)面S的距離為z；b.裂紋圓盤(pán)與觀測(cè)面S的交切示意圖

(5)

(6)

(7)

由于跡線長(zhǎng)度總小于或等于裂隙直徑，為了便于計(jì)算，設(shè)定c=φ，則有：

(8)

因?yàn)閏1≤φ≤φM,c0=φM時(shí)h(φ)=h1，可得：

此次從50個(gè)水平面中提取的裂隙跡線長(zhǎng)度的頻率直方圖見(jiàn)圖4?？梢?jiàn)巖樣剖面中的裂隙跡線長(zhǎng)度分布遵從對(duì)數(shù)正態(tài)分布，與一般研究中的發(fā)現(xiàn)一致。在考慮取樣偏差后，為了研究的簡(jiǎn)便，可以認(rèn)為跡線長(zhǎng)度分布服從冪律分布，即將較短的裂隙跡線長(zhǎng)度省去，與一般的裂隙觀察一致，從而可以用前述的理論方法進(jìn)行分析。

在運(yùn)用公式(9)計(jì)算裂隙直徑時(shí)，依據(jù)實(shí)際觀測(cè)數(shù)據(jù)，最小直徑φm設(shè)定為1 mm；最大直徑φM在1 mm到200 mm之間變化，本文取最大值200 mm；指數(shù)α一般選擇兩個(gè)值，分別是1.5和2.5，根據(jù)與現(xiàn)場(chǎng)裂隙分布相比較，取α=2.5，可以發(fā)現(xiàn)該值與實(shí)際情況相吻合。依據(jù)公式(2)—公式(6)，可分別求出α=1.500 5，β=0.422 6，b=1.5。

圖4 巖樣裂隙跡線長(zhǎng)度分布頻率Fig.4 Frequency distribution of trace length of rock fissures

4 結(jié)論

表2 裂紋直徑和密度計(jì)算結(jié)果

致謝：在野外工作期間，得到云南錫業(yè)公司松礦朱文捷高工的大力幫助，在此表示感謝！

[1] 張連明，周峰，張兆坤.地下巖層裂隙的幾種勘察研究方法[J].水文地質(zhì)工程地質(zhì)，2003，4：111-114. Zhang Lianming,Zhou Feng,Zhang Zhaokun.Several methods for investigating crack in underground rock stratum[J].Hydrogeology and Engineering Geology，2003,4:111-114.

[2] 白斌，鄒才能，朱如凱,等.四川盆地九龍山構(gòu)造須二段致密砂巖儲(chǔ)層裂縫特征、形成時(shí)期與主控因素[J].石油與天然氣地質(zhì)，2012，33(4)：526-535. Bai Bin，Zou Caineng，Zhu Rukai,et al.Characteristics,timing and controlling factors of structural fractures in tight sandstones of the 2nd member of Xujiahe Formation in Jiulong Mountain structure,Sichuan Basin[J].Oil & Gas Geology，2012，33(4)：526-535.

[3] 張仲培,劉士林,楊子玉,等.塔里木盆地麥蓋提斜坡構(gòu)造演化及油氣地質(zhì)意義[J].石油與天然氣地質(zhì),2011,32(6):909-919. Zhang Zhongpei,Liu Shilin,Yang Ziyu,et al.Tectonic evolution and its petroleum geological significances of the Maigaiti Slope,Tarim Basin[J].Oil & Gas Geology，2011,32(5):909-919.

[4] Berkowitz B,Adler P M.Stereological analysis of fracture network structure in geological formations[J].Geophysics Research,1998，B103:15339-15360.

[5] Adler P M，Thovert.Fractures and fracture network[M].USA:Klu-wer Academic Publishers,1999.

[6] Acuna J,Yeomans Y.Application of fractal geometry to the study of networks of fractures and their pressure transient[J].Water Resources

Research,1995,31(3):527-540.

[7] Dershowitz W S，Einstein H H.Characterizing rock joint geometry with joint system models[J].Rock Mechanics Rock Engineering,1988,21:21-51.

[8] 莊永秋，王任重，鄭樹(shù)培，等.云南個(gè)舊錫銅多金屬礦床[M].北京：地震出版社，1996. Zhuang Yongqiu,Wang Renzhong,Zheng Shupei,et al.Geology of Gejiu tin-copper polymetallic depositl(in China)[M].Beijing:Seismological Press,1996.

[9] 蘇建政，黃志文，龍秋蓮，等.基于ANSYS軟件的降低破裂壓力機(jī)理模擬[J].石油與天然氣地質(zhì)，2012，33(4)：640-645. Su Jianzheng,Huang Zhiwen,Long Qiulian,et al.ANSYS-based si-mulation of fracturing pressure reducing mechanism[J].Oil & Gas Geology,2012，33(4)：640-645.

[10] 蘇娜，段永剛，于春生.微CT掃描重建低滲氣藏微觀孔隙結(jié)構(gòu)——以新場(chǎng)氣田上沙溪廟組儲(chǔ)層為例[J].石油與天然氣地質(zhì),2011，32(5)：792-796. Su Na,Duan Yonggang,Yu Chunsheng.Reconstruction of microsco-pic pore structure in low permeability gas reservoir by micro-CT scanning:an example from the Upper Shaximiao Formation in Xinchang gas field[J].Oil & Gas Geology，2011，32(5)：792-796.

[11] Barton C A，Zoback M D.Self-similar distribution and properties of macroscopic fractures at depth in crystalline rocks in the Cajon Pass scientific drill hole[J].International Joural of Rock Mechanics and Ming Science，1993，30(1)：5181-5200.

[13] 倪春中，劉春學(xué)，張世濤.從巖石露頭裂隙跡線估算裂隙三維空間方向[J].石油與天然氣地質(zhì)，2013，34(1)：102-106. Ni Chunzhong,Liu Chunxue,Zhang Shitao.Estimation of three-dimensional distribution of fissures according to fissure traces on outcrops[J].Oil & Gas Geology,2013,34(1)：102-106.

(編輯 董 立)

Stereology-basedsimulationof3Dfracturenetworkwith2Ddata

Ni Chunzhong1，Liu Chunxue2，Zhang Shitao1，Li Aimin3，Li Yujian1，Cang Baohai4，Liu Lei5

(1.EngineeringInstituteofLandandResources,KunmingUniversityofScienceandTechnology,Kunming,Yunnan650093,China；2.SchoolofUrbanandEnvironment,YunnanUniversityofFinanceandEconomics,Kunming,Yunnan650221,China；3.YunnanTinGroupCo.Ltd.,GejiuYunnan661021,China；4.SpecialOilDevelopmentCompany,PetroChinaLiaoheOilfieldCompany,Panjin,Liaoning124010,China; 5.InstituteofEarthquakePrediction,ChinaSeismologicalBureau,Beijing100038,China)

Various attribute information of fissures such as length,width and density can be extracted from CT scanning images of oriented rock samples that are systematically collected through field work.These information can play a important role in defining spatial distribution of fissure networks within rocks.Both trace distribution and fissure network distribution follow the power-law distribution,thus statistical analysis can be performed on the field-collected trace data.The extension length and interlink of fissures on 2D images can be inversed to 3D space by applying stereological techniques,so as to deduce 3D fissure length and density distribution and realize conversion of rock fissure attributes from 2D to 3D.This approach was applied in Gaosong ore field.Based on analysis of the traces on the CT scanning images,we obtained the maximum trace length of 64.1 mm,the minimum of 2.2 mm,a mean of 16.69 mm and a variance of 12.33 mm.The fissures were divided into 9 groups,and the probability distribution of diameters of each group was obtained.Simulation was performed on the extracted trace data when the fractal dimension(α)is 1.5 and 2.5 respectively.The mean of fissure diameter was 27.60 mm,and fracture network density was 0.30/cm3.This approach can be popularized to other mines to obtain attributes of 3D fissure network when only 2D data are available.

CT section image,stereology technology,fissure network,rock fissure,Gejiu tin deposit

2013-07-02;

：2013-12-01。

倪春中(1979—)，男，博士、講師，構(gòu)造地質(zhì)學(xué)。

國(guó)家自然科學(xué)基金項(xiàng)目(40902058);云南省應(yīng)用基礎(chǔ)研究計(jì)劃項(xiàng)目(2013FZ026);云南省教育廳科學(xué)研究基金項(xiàng)目(2011Y363)。

0253-9985(2014)01-0148-05

10.11743/ogg20140119

TE19

：A