，，

(1.黃河勘測規(guī)劃設(shè)計有限公司 地質(zhì)工程院,鄭州 450003；2.中國地質(zhì)大學(xué)(北京) 水資源與環(huán)境學(xué)院,北京 100083)

1 研究背景

巖體由巖塊和分割它們的結(jié)構(gòu)面組成。結(jié)構(gòu)面是發(fā)育于巖體中，具有一定方向、延伸性和一定厚度的各種地質(zhì)界面，對巖體的強度、變形和滲透性有著決定性影響[1]。在野外，常采用地表露頭測量、探硐和鉆孔巖芯裂隙編錄等方法對結(jié)構(gòu)面數(shù)據(jù)進行采集。工程地質(zhì)、水文地質(zhì)工作中，一個基礎(chǔ)但十分重要的環(huán)節(jié)是對結(jié)構(gòu)面產(chǎn)狀數(shù)據(jù)進行統(tǒng)計分組研究，并確定每組結(jié)構(gòu)面的優(yōu)勢產(chǎn)狀。常用的方法主要是繪制傾向和走向玫瑰花圖以及在吳氏、施密特網(wǎng)上作極點等密度圖[2]。這些方法操作簡單，結(jié)果直觀，易于實現(xiàn)，但是由于地質(zhì)工作者的專業(yè)素養(yǎng)不同，對劃分結(jié)果的主觀判斷難免存在差異,當(dāng)結(jié)構(gòu)面的產(chǎn)狀實測數(shù)據(jù)數(shù)量多且界線較為模糊時更是如此。因此，應(yīng)借助數(shù)學(xué)手段對結(jié)構(gòu)面的產(chǎn)狀數(shù)據(jù)進行客觀劃分。

巖體結(jié)構(gòu)面產(chǎn)狀數(shù)據(jù)的計算機自動分組一直是國內(nèi)外工程地質(zhì)、水文地質(zhì)界研究的熱點問題,國外早在20世紀(jì)70年代便已經(jīng)開始了相關(guān)研究。Shanley和Mahtab[3]首次提出了一套完整的巖體結(jié)構(gòu)面產(chǎn)狀數(shù)據(jù)自動劃分的算法及其數(shù)值實現(xiàn)過程：通過對單位球面上的投影點進行密度分析完成優(yōu)勢組劃分。Hammah和Curran[4]使用模糊C均值聚類方法進行了結(jié)構(gòu)面優(yōu)勢組的自動劃分。Klose等[5]基于隨機逼近法，提出了一種無監(jiān)督的機器自學(xué)習(xí)法，快速完成優(yōu)勢組劃分。Jimenez-Rodriguez和Sitar[6]將譜分析方法和K-means聚類法相結(jié)合進行結(jié)構(gòu)面優(yōu)勢組劃分。Tokhmechi等[7]討論了采用結(jié)構(gòu)面的多個參數(shù)進行分組的重要性。我國在20世紀(jì)80年代也已開始了相關(guān)研究，并達到了國際領(lǐng)先水平[2, 8-9]。近些年來，國內(nèi)一些學(xué)者基于不同數(shù)學(xué)原理，先后提出或改進了結(jié)構(gòu)面自動劃分的算法，并在實際工程中進行了應(yīng)用，獲得了不錯的效果[10-14]。

本文以巖體結(jié)構(gòu)面法向向量為樣本，提出了基于K-means聚類方法和IIndex聚類有效性檢驗指標(biāo)的巖體結(jié)構(gòu)面自動分組方法，并開發(fā)了程序RDAP，服務(wù)于巖體變形、破壞和滲透性的分析與計算。

2 巖體結(jié)構(gòu)面自動分組方法

2.1 結(jié)構(gòu)面法向向量

野外測量得到的結(jié)構(gòu)面產(chǎn)狀數(shù)據(jù)通常使用傾向α和傾角β進行表示。根據(jù)空間解析幾何的相關(guān)知識，基于右手法則，可將產(chǎn)狀數(shù)據(jù)轉(zhuǎn)化為單位向量[2],即

(1)

式中nx，ny，nz分別為結(jié)構(gòu)面法向向量在x，y，z軸上的投影。

2.2 K-means聚類方法

K-means聚類方法是一種基于距離的聚類算法，采用距離作為相似性的評價指標(biāo)，認(rèn)為距離靠近的樣本組成聚類組，聚類目標(biāo)是使各聚類組總的距離平方和最小[15]。K-means算法的基本思想是初始隨機給定K個聚類組中心，并根據(jù)距離按照最鄰近原則把待分類的樣本分配給各個聚類組，然后重新計算各個聚類組的中心，一直迭代，直到聚類組中心的移動距離小于某個給定值。相較于其他聚類方法，K-means聚類方法具有快速、簡單的特點，適合大規(guī)模數(shù)據(jù)集的挖掘。

本文對傳統(tǒng)K-means聚類方法稍加改進，沿襲田開銘和萬力的方法[2]，使用結(jié)構(gòu)面法向向量與聚類組最大概率方向(最大概率方向等價于聚類組結(jié)構(gòu)面法向向量合矢量的方向)的夾角定義聚類距離，而不是K-means聚類方法通常使用的歐幾里德距離。這樣的改進突出了結(jié)構(gòu)面產(chǎn)狀數(shù)據(jù)的矢量性，具有更加顯著的工程、水文地質(zhì)意義。改進后的K-means聚類方法包括以下幾個步驟：

(1)根據(jù)預(yù)定義聚類組數(shù)，隨機指定若干聚類組中心，計算每個結(jié)構(gòu)面產(chǎn)狀樣本數(shù)據(jù)與聚類組中心之間的角度。

(2)將每個結(jié)構(gòu)面產(chǎn)狀樣本數(shù)據(jù)分配給與它夾角最小的聚類組中心，形成新的聚類組。

(3)將各組結(jié)構(gòu)面產(chǎn)狀樣本數(shù)據(jù)的最大概率方向作為新的聚類組中心。

(4)重復(fù)步驟(2)和步驟(3)，直至所有聚類組的中心固定，結(jié)構(gòu)面產(chǎn)狀樣本數(shù)據(jù)分配固定。

2.3 聚類有效性檢驗指標(biāo)

K-means聚類方法均屬于動態(tài)聚類法。動態(tài)聚類法通常無法直接確定聚類優(yōu)勢組數(shù)。實際應(yīng)用中，除了人為根據(jù)經(jīng)驗預(yù)先設(shè)定聚類組數(shù)外，常根據(jù)聚類有效性指標(biāo)確定最優(yōu)聚類組數(shù)。聚類有效性指標(biāo)較為常用的有：Davies-Bouldin (DB) Index，Dunn’s Index，Calinski Harabasz(CH) Index，Xie Beni(XB) Index[13]。近期，Maulik和Bandyopadhyay[16]提出了一種新型聚類有效性檢驗指標(biāo)Iindex,即

(2)

(3)

式中：N為結(jié)構(gòu)面聚類組數(shù)；v表示屬于第i聚類組的結(jié)構(gòu)面的法向單位向量；cNi表示第i聚類組結(jié)構(gòu)面中心向量，即最大概率方向;d(v,cNi)表示v和cNi兩向量之間的歐幾里得距離。

IIndex聚類有效性檢驗指標(biāo)使用各聚類組中心距離的最大值度量組間分離度，使用聚類組中各成員與組中心的距離和度量組內(nèi)分離度。IIndex聚類有效性檢驗指標(biāo)認(rèn)為：IIndex取值愈大，其所對應(yīng)的分組效果愈佳。此外，實例驗證，當(dāng)聚類距離定義為向量夾角時，相較于上述常用的聚類有效性檢驗指標(biāo)，使用IIndex聚類有效性檢驗指標(biāo)選擇優(yōu)勢聚類組數(shù)的效果更好[13]。

2.4 巖體結(jié)構(gòu)面自動分組程序RDAP

RDAP是巖體結(jié)構(gòu)面自動分組程序英文名稱Rock Discontinuities Auto-classification Program的簡稱。RDAP使用MatLab科學(xué)計算語言編制計算程序和GUI界面。原理上，RDAP首先根據(jù)Iindex進行聚類有效性檢驗，選取Iindex最大值所對應(yīng)的聚類組數(shù)作為結(jié)構(gòu)面優(yōu)勢組數(shù)，然后將其作為預(yù)定義聚類組數(shù)代入改進的K-means聚類算法實現(xiàn)巖體結(jié)構(gòu)面產(chǎn)狀數(shù)據(jù)的計算機批量自動分組(圖1)。RDAP的自動分組結(jié)果可作為Dips,Unwedge,Steronet等商業(yè)軟件的輸入數(shù)據(jù)，服務(wù)于巖體變形、破壞和滲透性的分析與計算。

圖1 I index聚類有效性檢驗指標(biāo)隨聚類優(yōu)勢組數(shù)的變化Fig.1 RDAP interface showing the plot of cluster validity index I versus cluster numbers

3 算法驗證

Shanley和Mahtab[3]于1976年首次提出了一種結(jié)構(gòu)面優(yōu)勢組自動劃分的方法及其數(shù)值實現(xiàn)過程，并以算例的形式給出了286條結(jié)構(gòu)面數(shù)據(jù)及自動分組結(jié)果。學(xué)術(shù)界對這組數(shù)據(jù)的分組形式已經(jīng)形成了較為統(tǒng)一的認(rèn)識，并常將其作為檢驗其他分組方法正確與否的試金石[5]。Shanley和Mahtab使用其算法，將這286條結(jié)構(gòu)面數(shù)據(jù)分為了3個聚類組(圖2)。

圖2 Shanley和Mahtab[3] 分組結(jié)果與RDAP分組結(jié)果的對比Fig.2 Classing results by Shanley and Mahtab[3] and by RDAP

Iindex聚類有效性檢驗指標(biāo)隨聚類組數(shù)的變化圖表明：Iindex聚類有效性檢驗指標(biāo)在聚類組數(shù)為3時等于1.72，為最大值(圖1)。因此，根據(jù)本文所提分組方法，286條結(jié)構(gòu)面數(shù)據(jù)最優(yōu)可以分為3個聚類組，這與Shanley和Mahtab所識別的最優(yōu)聚類組數(shù)目一致。具體來講，使用RDAP，第1組識別得到了114條結(jié)構(gòu)面，第2組59條結(jié)構(gòu)面，第3組113條結(jié)構(gòu)面(圖1)。為了直觀地進行比較，我們將Shanley和Mahtab的分組結(jié)果和RDAP的分組結(jié)果分別繪制在極點圖上(圖2)。經(jīng)比較，兩者的分組結(jié)果基本一致，證明了本文所提分組方法的可靠性。

需要強調(diào)的是，除了可以給出分組結(jié)果，RDAP還能根據(jù)聚類組的最大概率方向給出每組結(jié)構(gòu)面的優(yōu)勢、代表性產(chǎn)狀(圖1)，便于進行后續(xù)工程地質(zhì)穩(wěn)定性分析、水文地質(zhì)的計算和防滲灌漿鉆孔優(yōu)勢方位的選定。

4 實例應(yīng)用

某地下工程位于廣東江門市附近，坐落于花崗巖巖體內(nèi)部。工程地下建筑物主要為斜井、豎井及實驗廳3部分。目前，工程面臨的主要問題是施工涌水。

4.1 涌水裂隙分組

豎井在施工過程中觀測并記錄了42條主要涌水裂隙(表1)。將這42條涌水裂隙的傾向和傾角數(shù)據(jù)作為輸入數(shù)據(jù)，選擇“程序自動分組方法”。RDAP將這42條主要涌水裂隙自動識別為5個裂隙組(表1)，并可給出這5個裂隙組的最大概率方向，即優(yōu)勢產(chǎn)狀(表2)。由于第2和第3組涌水裂隙數(shù)據(jù)少(只有2條)，代表性差，可將其舍去，認(rèn)為花崗巖巖體中主要發(fā)育3個裂隙組。在自動分組結(jié)果的基礎(chǔ)上，根據(jù)組內(nèi)每條涌水裂隙的出露高程(剔除異常點)、隙寬描述和每組裂隙的優(yōu)勢產(chǎn)狀，可以計算得到每組裂隙的真張開度和真裂隙間距(表2)。

4.2 灌漿鉆孔最佳方位的選取

涌水量的預(yù)測與防治是該工程施工中的重點與難點。灌漿防滲處理技術(shù)是工程涌水防治的有效、常規(guī)方法。其中，灌漿鉆孔最佳方位的選取是影響灌漿效果的關(guān)鍵因素之一。

假設(shè)灌漿鉆孔的傾向為a，傾角為b(與水平面之間的夾角)，灌漿鉆孔的空間向量可以表示為

(4)

式中m，n，p為灌漿鉆孔在x，y，z軸上的投影。

那么，灌漿鉆孔揭露的視裂隙間距與真裂隙間距的關(guān)系可以表示為

(5)

表1 豎井涌水裂隙基本情況Table 1 Basic information of gushing fractures

表2 涌水裂隙統(tǒng)計分組結(jié)果Table 2 Classing result of gushing fractures

式中：sa為視裂隙間距；s為真裂隙間距；θ為灌漿鉆孔與裂隙面法向向量之間的夾角。

灌漿鉆孔的最佳方位是在同等壓力下吃漿量最大的方位。換言之，應(yīng)是灌漿鉆孔穿過最多、最寬裂隙的方位[17]。若忽略裂隙寬度、粗糙度、充填情況等因素對灌漿效果的影響，灌漿鉆孔最佳方位的選取等價于單位長度的灌漿鉆孔盡可能多地與各組裂隙相交，即獲得最大的灌漿鉆孔綜合視裂隙密度的方位。轉(zhuǎn)化成數(shù)學(xué)語言，即為求解以下目標(biāo)函數(shù)的最大值，其表達式為

(6)

圖3 綜合視裂隙密度三維渲染圖Fig.3 Three-dimensional plot of the apparent fracture density

5 結(jié) 語

(1)以巖體結(jié)構(gòu)面法向向量為樣本，提出了基于K-means聚類方法和IIndex聚類有效性檢驗指標(biāo)的巖體結(jié)構(gòu)面自動分組方法，并開發(fā)了巖體結(jié)構(gòu)面自動分組程序RDAP。

(2)通過與經(jīng)典文獻進行對比，本文所提結(jié)構(gòu)面自動分組方法的可靠性得到了論證。

(3)使用本文提出的結(jié)構(gòu)面自動分組方法對某地下工程的實測涌水裂隙資料進行了預(yù)處理，初步計算了灌漿鉆孔的最佳方位，為工程涌水的防治提供了依據(jù)。