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

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

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

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

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

1.1 Fisher 分布

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

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

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

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

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

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

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

首先根據(jù)巖石面的產(chǎn)狀(α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 的點乘確定:

式中:θ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表示在巖墻面中從水平面到跡線的順時針角，它可作為跡線角的唯一度量值。

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

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

2.1 模擬退火法

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

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

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

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

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

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

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

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

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

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

3)計算轉(zhuǎn)移概論P:

式中:Xm+1為迭代計算產(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ù)據(jù)估算裂隙三維空間分布的有效性［14-15］，利用個舊錫礦區(qū)高松礦田某一開挖面上得到的裂隙跡線數(shù)據(jù)進行了實例分析，將計算得到的節(jié)理產(chǎn)狀與實測數(shù)據(jù)進行了比較。圖2a 為個舊錫礦區(qū)高松礦田某臨礦山公路的開挖面，長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)相對應(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)，計算出K 值分別等于20，50，100，150 和200 時Fisher 分布。圖4 表示符合Fisher 分布，傾向和傾角均值分別為120°和20°的節(jié)理組，圖4a 中k=50，圖4b 中k=200。本文通過與現(xiàn)場比較，取K 值等于200，計算中所用的參數(shù)A，B，C 所占的權(quán)重分別取70%，30%和0，求得誤差5%范圍內(nèi)θrake和SD(標(biāo)準(zhǔn)方差)值(表1)。

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

4 結(jié)論

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

［1］陳田勇，毛鑫，劉仕銀，等. 利用分形理論計算相對滲透率曲線——以南襄盆地雙河油田核桃園組六油組為例［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］.石油實驗地質(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］馮陣東，戴俊生，鄧航，等. 利用分形幾何定量評價克拉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］余一欣，周心懷，徐長貴，等. 渤海海域新生代斷裂發(fā)育特征及形成機制［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］.中國水利水電科學(xué)研究院學(xué)報，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)組碳酸鹽巖儲層特征研究［J］. 石油實驗地質(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.