基于K-means聚類方法和I Index聚類有效性檢驗(yàn)指標(biāo)的巖體結(jié)構(gòu)面自動(dòng)分組及應(yīng)用

，，

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

1 研究背景

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

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

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

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

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

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

(1)

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

2.2 K-means聚類方法

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

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

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

(2)將每個(gè)結(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 聚類有效性檢驗(yàn)指標(biāo)

K-means聚類方法均屬于動(dòng)態(tài)聚類法。動(dòng)態(tài)聚類法通常無(wú)法直接確定聚類優(yōu)勢(shì)組數(shù)。實(shí)際應(yīng)用中，除了人為根據(jù)經(jīng)驗(yàn)預(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]提出了一種新型聚類有效性檢驗(yàn)指標(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聚類有效性檢驗(yàn)指標(biāo)使用各聚類組中心距離的最大值度量組間分離度，使用聚類組中各成員與組中心的距離和度量組內(nèi)分離度。IIndex聚類有效性檢驗(yàn)指標(biāo)認(rèn)為：IIndex取值愈大，其所對(duì)應(yīng)的分組效果愈佳。此外，實(shí)例驗(yàn)證，當(dāng)聚類距離定義為向量夾角時(shí)，相較于上述常用的聚類有效性檢驗(yàn)指標(biāo)，使用IIndex聚類有效性檢驗(yàn)指標(biāo)選擇優(yōu)勢(shì)聚類組數(shù)的效果更好[13]。

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

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

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

3 算法驗(yàn)證

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

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

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

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

4 實(shí)例應(yīng)用

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

4.1 涌水裂隙分組

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

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

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

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

(4)

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

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

(5)

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

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

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

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

(6)

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

5 結(jié) 語(yǔ)

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

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

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