基于體視學技術的二維數據模擬三維裂隙網絡

倪春中，劉春學，張世濤，李愛敏，李雨健，蒼寶海，劉 雷

(1.昆明理工大學 國土資源工程學院，云南 昆明 650093; 2.云南財經大學 城市與環境學院，云南 昆明 650221;3.云南錫業集團有限責任公司，云南 個舊 661021; 4.中國石油 遼河油田分公司 特種油開發公司，遼寧 盤錦 124010;5.中國地震局 地震預測研究所，北京 100038)

基于體視學技術的二維數據模擬三維裂隙網絡

倪春中1，劉春學2，張世濤1，李愛敏3，李雨健1，蒼寶海4，劉 雷5

(1.昆明理工大學 國土資源工程學院，云南 昆明 650093; 2.云南財經大學 城市與環境學院，云南 昆明 650221;3.云南錫業集團有限責任公司，云南 個舊 661021; 4.中國石油 遼河油田分公司 特種油開發公司，遼寧 盤錦 124010;5.中國地震局 地震預測研究所，北京 100038)

CT掃描；體視學；裂隙網絡；巖石裂隙；個舊錫礦

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

本研究以云南個舊錫礦區高松礦田為例，所用的小尺度裂隙數據取自樣本巖石的CT掃描圖像，應用體視學技術的方法將二維圖像中的裂隙延展及其連接等數據反演到三維空間，然后在三維空間相關特征的約束條件下，反演推斷出三維空間中裂隙密度分布，最終實現應用二維裂隙數據模擬三維裂隙網絡。

1 地質概況

高松礦田是云錫公司的主要生產礦山之一，位于個舊錫礦東礦區北部，處于礦區的五子山復背斜北段，

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

褶皺構造有近EW向對門山-阿西寨向斜、近SN向五子山背斜等。斷裂構造可以分為EW，NE，NW，SN四組，其中以近EW向和NE向為主，NW向次之。EW向斷裂組自北而南有個松斷裂、麒麟山斷裂、馬吃水斷裂、高阿斷裂、背陰山斷裂，呈近等距離分布，其中個松和背陰山斷裂規模較大，麒麟山斷裂次之。NE向斷裂組自西而東有蓮花山斷裂、蘆塘壩斷裂、麒阿西斷裂。NW向斷裂組主要為大箐東斷裂、黑碼石斷裂、駝峰山斷裂、阿西寨斷裂。SN向斷裂組不發育(圖1)。

2 樣本裂隙提取

2.1 巖石樣本

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

圖1 個舊礦區高松礦田構造地質略圖(文獻[8]，修改)Fig.1 Tectonic geological sketch of Gaosong mine in GejiuT2g1.中三疊統個舊組卡房段白云巖灰巖；T2g2.中三疊統個舊組馬拉格段白云巖；T2g3.中三疊統個舊組白泥硐段灰巖

樣品編號001002003004005006007008009010011012采樣位置/m192019201720172015401540154015401360136013601360

2.2 裂隙跡線提取

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

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

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

3 基于體視學方法由二維跡線模擬裂隙長度及密度

體視學利用形態學和數學的分析，借助計算機及數據處理系統和顯微鏡及顯做成像系統，將二維平面

經過成像及計算機分析處理得到三維形態，準確地對物體進行定量及形態結構分析。本文應用體視學技術，研究圓盤裂隙三維空間與二維剖面之間的相互聯系，解析裂隙長度在不同維數之間的聯系和轉化，以處理二維剖面中的樣本裂隙跡線數據。

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

假設φ為常數，直徑為φ的圓盤與觀測面交切的比例為S(z,φ)。從圖3可以看出，只有當|z|<φ/2時才會有交切。對于給定的z，S(z,φ)可以表示為

(1)

式中：z為裂隙圓盤中心到觀測面S的距離，mm；φ為圓盤直徑，mm。

下述公式應用冪次定律可以計算裂隙直徑分布[11-13]：

(2)

式中：α為分形維數，對于裂隙跡線長度來說，1≤α≤3；a為待定參數，由式(3)確定。

如果裂隙直徑長度假定在φm≤φ≤φM的范圍內，那么常數α可以表示為：

(3)

(4)

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

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

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

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

(5)

(6)

(7)

由于跡線長度總小于或等于裂隙直徑，為了便于計算，設定c=φ，則有：

(8)

因為c1≤φ≤φM,c0=φM時h(φ)=h1，可得：

此次從50個水平面中提取的裂隙跡線長度的頻率直方圖見圖4。可見巖樣剖面中的裂隙跡線長度分布遵從對數正態分布，與一般研究中的發現一致。在考慮取樣偏差后，為了研究的簡便，可以認為跡線長度分布服從冪律分布，即將較短的裂隙跡線長度省去，與一般的裂隙觀察一致，從而可以用前述的理論方法進行分析。

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

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

4 結論

表2 裂紋直徑和密度計算結果

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

[1] 張連明，周峰，張兆坤.地下巖層裂隙的幾種勘察研究方法[J].水文地質工程地質，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] 白斌，鄒才能，朱如凱,等.四川盆地九龍山構造須二段致密砂巖儲層裂縫特征、形成時期與主控因素[J].石油與天然氣地質，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] 張仲培,劉士林,楊子玉,等.塔里木盆地麥蓋提斜坡構造演化及油氣地質意義[J].石油與天然氣地質,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] 莊永秋，王任重，鄭樹培，等.云南個舊錫銅多金屬礦床[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].石油與天然氣地質，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掃描重建低滲氣藏微觀孔隙結構——以新場氣田上沙溪廟組儲層為例[J].石油與天然氣地質,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] 倪春中，劉春學，張世濤.從巖石露頭裂隙跡線估算裂隙三維空間方向[J].石油與天然氣地質，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—)，男，博士、講師，構造地質學。

國家自然科學基金項目(40902058);云南省應用基礎研究計劃項目(2013FZ026);云南省教育廳科學研究基金項目(2011Y363)。

0253-9985(2014)01-0148-05

10.11743/ogg20140119

TE19

：A