倪春中，劉春學(xué)，張世濤

(1.昆明理工大學(xué)國(guó)土資源工程學(xué)院，云南昆明650093; 2. 云南財(cái)經(jīng)大學(xué)城市與環(huán)境學(xué)院，云南昆明650221)

裂隙是地球科學(xué)中的一類重要構(gòu)造，在地殼中廣泛分布，與裂隙網(wǎng)絡(luò)滲流關(guān)系密切，因此查清裂隙的空間分布顯得尤為重要。但目前獲取裂隙空間分布特征的途徑，除了在實(shí)驗(yàn)室內(nèi)進(jìn)行小規(guī)模的分形研究外［1-4］，大多數(shù)是通過(guò)鉆孔，坑道開(kāi)挖面以及一些巖石露頭獲取。野外裂隙調(diào)查的方法主要有測(cè)線法［5-7］和統(tǒng)計(jì)窗法［8］。但在實(shí)際工作中存在兩個(gè)問(wèn)題:首先是不同的取樣方法，人為測(cè)量偏差以及儀器誤差都有可能影響結(jié)果精度;其次在野外由于巖石松動(dòng)以及巖塊掉落等原因，運(yùn)用上述兩種方法直接對(duì)陡峭的巖墻進(jìn)行測(cè)量存在很大的安全隱患。本文通過(guò)對(duì)數(shù)字影像提取裂隙痕跡，克服了上述方法的缺點(diǎn)。利用數(shù)字影像處理軟件，提取裂隙痕跡屬性，包括走向、傾角、裂隙間距、大小以及形狀等。裂隙三維空間分布與裂隙痕跡的這些屬性有關(guān)，尤其是走向、傾角。每組裂隙的走向、傾角都有一個(gè)均值，圍繞著這個(gè)均值散布著眾多裂隙。一般認(rèn)為裂隙面法線方向的概率密度可以用Fisher 分布進(jìn)行描述，運(yùn)用Fisher分布密度函數(shù)的前提條件是對(duì)裂隙進(jìn)行分組，根據(jù)參數(shù)K，模擬出裂隙面產(chǎn)狀，同時(shí)由于該函數(shù)是可積函數(shù)，產(chǎn)生隨機(jī)數(shù)較方便。本文先考慮已知巖墻面的產(chǎn)狀，裂隙平均走向、傾角以及Fisher 分布參數(shù)K，計(jì)算巖墻面裂隙痕跡與水平面的夾角，估算出裂隙網(wǎng)絡(luò)面的屬性參數(shù)，在此基礎(chǔ)上，利用模擬退火法對(duì)原問(wèn)題求逆，并將該方法應(yīng)用于個(gè)舊錫礦區(qū)高松礦田，考察其預(yù)測(cè)裂隙網(wǎng)絡(luò)分布的有效性。

1 裂隙面產(chǎn)狀確定其跡線的射線角

如前所述，裂隙方向?qū)α严兜钠渌鼘傩杂绊懼陵P(guān)重要。每組裂隙都有一個(gè)均值，裂隙方向圍繞著這個(gè)均值分布，F(xiàn)isher 分布參數(shù)K 可以用來(lái)生成這些離散值。

1.1 Fisher 分布

巖體裂隙面在空間的傾角用其法線的方向表示。應(yīng)用球面坐標(biāo)系可以給出Fisher 分布模型的表達(dá)式［9-10］為:

式中:θ＇為經(jīng)過(guò)變換后的裂隙面傾角，(°);K 表示裂隙組極點(diǎn)集的集中程度參數(shù)，無(wú)量綱。對(duì)巖石裂隙來(lái)說(shuō)，K 值取值一般介于20～300，其取值方法參見(jiàn)文獻(xiàn)［5］。根據(jù)該公式可以推導(dǎo)出Fisher分布函數(shù)在極點(diǎn)圖上用極坐標(biāo)表示的面積密度，即以真實(shí)平均值θ 為中心的角面積dA 中出現(xiàn)裂隙的概率，其函數(shù)表示為:

式中:θ 為裂隙面傾角真實(shí)平均值，(°)。當(dāng)θ =0時(shí)，f(θ)值最大。該式反映圍繞裂隙均值的對(duì)稱分布。K 值越大，方向變量越集中于真實(shí)平均值。

可以用下述公式近似估計(jì)計(jì)算機(jī)產(chǎn)生的Fisher分布［11］:

式中:Δθ 為裂隙與均值的差，(°);Random(0，1)為介于0 到1 的均勻隨機(jī)數(shù)。為了確定新的傾角和傾向，Δθ 圍繞平均值旋轉(zhuǎn)一個(gè)介于0～360°的任意角。

1.2 裂隙方向和跡線角之間的關(guān)系

定義節(jié)理平面法向量的傾向和傾角為(αj，βj)，巖石面法向量的傾向和傾角為(αf，βf)(圖1)。

首先根據(jù)巖石面的產(chǎn)狀(αj，βj)可以分別計(jì)算垂直于節(jié)理平面(j 向量)和巖石面(f 向量)的法線坐標(biāo)，跡線矢量t 為矢量j 與矢量f 的叉乘t =j×f，裂隙走向的單位矢量可以表示為s =(sx，sy，sz)=(cosαf，sinαf，0)。

裂隙跡線角θt可以用s 與t 的點(diǎn)乘確定:

式中:θt為跡線與水平面的銳夾角，(°);sx，sy，sz分別為裂隙走向單位矢量在x，y，z 軸上的分量;tx，ty，tz分別為跡線方向單位矢量在x，y，z 軸上的分量。

圖1 裂隙方向與跡線角關(guān)系Fig.1 Relationship between fissure direction and trace dip

由于θt可出現(xiàn)在不同的象限內(nèi)，不是巖石面上跡線角的唯一值，為了確保跡線角的唯一性，本文采用射線角θrake表示在巖墻面中從水平面到跡線的順時(shí)針角，它可作為跡線角的唯一度量值。

2 二維裂隙痕跡反演其三維空間分布

基于以上數(shù)據(jù)，本文利用模擬退火法，對(duì)逆問(wèn)題求解，就可以得到裂隙的傾向和傾角。

2.1 模擬退火法

模擬退火優(yōu)化方法最早由Kirkpatrick 和Cerny 等人提出和研究［12-13］，分析了優(yōu)化領(lǐng)域里的許多優(yōu)化組合和運(yùn)用統(tǒng)計(jì)物理求解物理系統(tǒng)最低能量狀態(tài)問(wèn)題，并基于此導(dǎo)出了模擬退火尋優(yōu)法的基本框架。

考慮一個(gè)物理系統(tǒng)s，設(shè)s 共有N 個(gè)狀態(tài)。當(dāng)達(dá)到平衡后，系統(tǒng)將隨機(jī)地處于一種狀態(tài)，并且處于狀態(tài)xi的概率滿足Boltzmann 分布，即:

式中:ui為s 處于狀態(tài)xi時(shí)具有的能量，J;Tk為系統(tǒng)處于熱平衡狀態(tài)下對(duì)應(yīng)的溫度，℃。

從Boltzmann 分布可以看出，溫度越低，系統(tǒng)處于具有小能量狀態(tài)的概率也越大，當(dāng)溫度趨于絕對(duì)零度時(shí)系統(tǒng)處于具有最小能量狀態(tài)的概率趨于1。所以，為了使物體s 具有理想的規(guī)則晶格結(jié)構(gòu)(對(duì)應(yīng)于能量最小的狀態(tài))，通常的作法是給s加足夠高的溫度，使其達(dá)到液化狀態(tài)，然后緩慢地逐步降低溫度，直到溫度接近于絕對(duì)溫度。這樣，s 將隨著溫度的下降，逐漸地從高能量狀態(tài)過(guò)渡到能量越來(lái)越低的狀態(tài)，直至能量小或接近最小。這個(gè)過(guò)程，物理上稱為退火過(guò)程。

顯然，如果把優(yōu)化問(wèn)題的狀態(tài)看成物理系統(tǒng)的狀態(tài)，目標(biāo)函數(shù)值看成物理系統(tǒng)的能量，則優(yōu)化問(wèn)題和物理系統(tǒng)求得最低或接近最低能量狀態(tài)問(wèn)題是一致的。模擬退火尋優(yōu)法就是基于此而導(dǎo)出的。

2.2 求解裂隙面產(chǎn)狀

由統(tǒng)計(jì)物理的退火過(guò)程，可以構(gòu)造如下模擬退火仿真優(yōu)化算法。

1)設(shè)定初始狀態(tài)X(0)，定義節(jié)理平面法向量的傾向?yàn)棣羓和傾角為βj，并假定αj0=360°和傾角βj0=90°，C=C0，C0為設(shè)置的控制參數(shù)。

2)隨機(jī)產(chǎn)生新?tīng)顟B(tài)X(m)，傾向αjm和傾角為βjm，估計(jì)優(yōu)化目標(biāo)值:

式中:μ，σ，θ 分別表示平均值、標(biāo)準(zhǔn)偏差以及偏斜度;A，B，C 為常數(shù);I 表示二維圖像中痕跡方向角分布的參數(shù);T 表示根據(jù)公式(2)-公式(9)由模擬退火法產(chǎn)生的參數(shù)。

3)計(jì)算轉(zhuǎn)移概論P(yáng):

式中:Xm+1為迭代計(jì)算產(chǎn)生的新解;f(Xm+1)為式(6)代入新解求得的優(yōu)化目標(biāo)函數(shù)值;c 為初始控制參數(shù)。

4)判斷接受或放棄新解Xm+1。當(dāng)P(Xm+1=Xm)＞r(r 為設(shè)定的允許誤差)，接受Xm為新解;否則，放棄Xm+1為新解。

5)判斷終止條件，滿足條件則停止迭代;否則轉(zhuǎn)向2)。

3 實(shí)例分析

為了考察二維裂隙數(shù)據(jù)估算裂隙三維空間分布的有效性［14-15］，利用個(gè)舊錫礦區(qū)高松礦田某一開(kāi)挖面上得到的裂隙跡線數(shù)據(jù)進(jìn)行了實(shí)例分析，將計(jì)算得到的節(jié)理產(chǎn)狀與實(shí)測(cè)數(shù)據(jù)進(jìn)行了比較。圖2a 為個(gè)舊錫礦區(qū)高松礦田某臨礦山公路的開(kāi)挖面，長(zhǎng)16.2 m，高約3 m。該巖墻面的產(chǎn)狀為255°∠80°，巖性主要為白云質(zhì)灰?guī)r，發(fā)育有多組節(jié)理。經(jīng)數(shù)字影像處理，獲取179 條裂隙痕跡，提取了裂隙痕跡傾角三維分布圖(圖2b)。根據(jù)裂隙痕跡與水平面夾角θrake值分布直方圖，巖墻面裂隙痕跡以90°為界分為兩組(圖3)，均值分別為49.0°和128.2°。

圖2 巖墻面裂隙痕跡分布影像a)相對(duì)應(yīng)的裂隙痕跡三維分布b)Fig.2 Distribution of fissure trace in rock surface image (a)corresponding with the 3D distribution of the fissure trace (b)

圖3 裂隙痕跡θrake值頻度直方圖Fig.3 Histogram showing the value of fissure trace

依據(jù)裂隙叢聚規(guī)律，節(jié)理面傾向取值范圍為［0°，359°］，傾角取值范圍為［0°，90°］。根據(jù)公式(2)，計(jì)算出K 值分別等于20，50，100，150 和200 時(shí)Fisher 分布。圖4 表示符合Fisher 分布，傾向和傾角均值分別為120°和20°的節(jié)理組，圖4a 中k=50，圖4b 中k=200。本文通過(guò)與現(xiàn)場(chǎng)比較，取K 值等于200，計(jì)算中所用的參數(shù)A，B，C 所占的權(quán)重分別取70%，30%和0，求得誤差5%范圍內(nèi)θrake和SD(標(biāo)準(zhǔn)方差)值(表1)。

為了求得裂隙面的三維空間分布，根據(jù)所提取裂隙痕跡數(shù)據(jù)，利用Matlab 語(yǔ)言編程，然后根據(jù)最優(yōu)函數(shù)公式(6)及判別式(7)求解最小誤差(設(shè)定r值≤0.05)。裂隙面產(chǎn)狀估算結(jié)果如表1 所示，可見(jiàn)模擬情況比較接近巖體裂隙的實(shí)際分布情況。

4 結(jié)論

圖4 符合Fisher 分布的節(jié)理組極點(diǎn)圖(傾向120°，傾角20°)Fig.4 Point diagram of joint set in agreement with the Fisher distribution (mean dip direction of 120°and mean dip of 20°)

表1 裂隙痕跡模擬結(jié)果Table 1 Simulation of fissure traces

由巖墻面裂隙痕跡估算裂隙的空間產(chǎn)狀，可以通過(guò)影像數(shù)字處理，運(yùn)用計(jì)算機(jī)自動(dòng)提取裂隙跡線信息，有效克服了在陡峭的巖墻面上無(wú)法運(yùn)用測(cè)線法和統(tǒng)計(jì)窗法進(jìn)行裂隙測(cè)量的缺點(diǎn)，消除了人為運(yùn)用上述兩種方法造成的誤差，同時(shí)也可以彌補(bǔ)在野外工作中遺漏的部分裂隙信息。裂隙跡線進(jìn)行分組后，確定每組裂隙跡線傾角均值，根據(jù)現(xiàn)場(chǎng)裂隙分布的集中程度，可以確定Fisher 分布常數(shù)K，結(jié)合巖墻面方向以及鏡頭方向與巖墻面法線夾角，利用模擬退火優(yōu)化仿真算法，估算得到裂隙的傾向和傾角，實(shí)現(xiàn)從二維裂隙生成三維裂隙網(wǎng)絡(luò)，大大減少了人工現(xiàn)場(chǎng)測(cè)試裂隙工作量。

［1］陳田勇，毛鑫，劉仕銀，等. 利用分形理論計(jì)算相對(duì)滲透率曲線——以南襄盆地雙河油田核桃園組六油組為例［J］.石油與天然氣地質(zhì)，2012，33(4):578-581.Chen Tianyong，Mao Xin，Liu Shiyin. Fractal theory-based calculation method of relative permeability curves—a case study from the Hetaoyuan Formation in Shuanghe oilfield，Nanxiang Basin［J］.Oil ＆ Gas Geology.2012，33(4):578-581.

［2］王津義，彭金寧，王彥清，等.貴州南部印支晚期——燕山期構(gòu)造變形特征［J］.石油實(shí)驗(yàn)地質(zhì)，2012，34(4):362-367.Wang Jinyi，Peng Jinning，Wang Yanqing，et al.Tectonic deformation features from Late Indosinian to Yanshanian in southern Guizhou［J］.Petroleum Geology ＆ Experiment，2012，34(4):362-367.

［3］馮陣東，戴俊生，鄧航，等. 利用分形幾何定量評(píng)價(jià)克拉2氣田裂縫［J］.石油與天然氣地質(zhì)，2011，32(6):928-933.Feng Zhendong，Dai Junsheng，Deng Hang，et al. Quantitative evaluation of fractures with fractal geometry in Kela-2 gas field［J］.Oil ＆ Gas Geology.2011，32(5):928-933.

［4］郭愷，仝兆岐，滕厚華.地震采集面元尺度與成像橫向分辨能力分析［J］.石油與天然氣地質(zhì)，2012，33(1):141-147.Guo Kai，Tong Zhaoqi，Teng Houhua. Analysis on bin size of seismic survey and lateral resolution of imaging［J］.Oil ＆ Gas Geology.2012，33(1):141-147.

［5］Priest S D，Hudson J A.Estimation of discontinuity spacing and trace length using scanline surveys［J］.International Journal of Rock Mechanics，Mining Sciences and Geomechanics Abstracts，1981，18(3):183-197.

［6］Cruder D M.Describing the size of discontinuities［J］.International Journal of Rock Mechanics and Mining Sciences and GeomechanicsAbstracts，1977，l4(3):133-137.

［7］余一欣，周心懷，徐長(zhǎng)貴，等. 渤海海域新生代斷裂發(fā)育特征及形成機(jī)制［J］.石油與天然氣地質(zhì)，2011，32(2):273-279.Yu Yixin，Zhou Xinhuai，Xu Changgui，et al. Characteristics and formation mechanism of the Cenozoic faults in the Bohai Sea waters［J］.Oil ＆ Gas Geology，2009，30(2):273-279.

［8］Paul P H. Estimating the mean length of discontinuity trace［J］.International Journal of Rock Mechanics and Mining Sciences and GeomechanicsAbstracts，1981，18(3):221-228.

［9］陸峰，王俊奇.Fisher 模型在巖體裂隙面模擬中的應(yīng)用［J］.中國(guó)水利水電科學(xué)研究院學(xué)報(bào)，2011，8(4):309-313.Lu Feng，Wang Junqi.Application of Fisher Model in rock fracture simulation［J］. Journal of China Institute of Water Resources and Hydropower Research，2011，8(4):309-313.

［10］Liu Chunxue，Koike K. Extending multivariate space-time geostatistics for environmental data analysis［J］.Mathematical Geology，2007，39(2):289-305.

［11］Priest S D. Discontinuity analysis for rock engineering［M］.London:Chapman ＆ Hall，1993:473.

［12］Kirkpatrick S，Gelatt C D，Vecchi M P.Optimization by simulated annealing［J］.Science，1983，220(4598):671-680.

［13］Cerny V.Thermodynamical approach to the traveling salesman problem:an efficient simulation algorithm［J］.Journal of Optimization Theory and Applications，1985，45(1):41-51.

［14］Koike K，Liu Chunxue，Sanga T. Incorporation of fracture directions into 3D geostatistical methods for a rock fracture system［J］. Environment Earth Science，2012，66(5):1403-1414.

［15］李旭兵，劉安，曾雄偉，等.雪峰山西側(cè)地區(qū)寒武系婁山關(guān)組碳酸鹽巖儲(chǔ)層特征研究［J］. 石油實(shí)驗(yàn)地質(zhì)，2012，34(2):153-157.Li Xubing，Liu An，Zeng Xiongwei，et al. Characteristics of carbonate reservoirs of Cambrian Loushanguan Formation to the west of Xuefeng Mountain［J］. Petroleum Geology ＆ Experiment，2012，34(2):153-157.