王述紅， 朱寶強(qiáng)， 王鵬宇

(東北大學(xué) 資源與土木工程學(xué)院, 遼寧 沈陽 110819)

巖體是由結(jié)構(gòu)面和巖石組成的復(fù)雜塊體，它作為一種非均質(zhì)材料，破壞時往往是沿著結(jié)構(gòu)面斷裂，因此，結(jié)構(gòu)面決定了巖體的強(qiáng)度及穩(wěn)定性[1].然而，由于野外環(huán)境的復(fù)雜性，結(jié)構(gòu)面的很多特征參數(shù)獲取較為困難，并且在利用數(shù)學(xué)手段對結(jié)構(gòu)面進(jìn)行分組時，分組參數(shù)的增加則使得計算的維度大大增加，數(shù)據(jù)之間的差異性使得其相似性也難以度量，因此目前最常用的仍然是依據(jù)結(jié)構(gòu)面的產(chǎn)狀數(shù)據(jù)(即傾向和傾角)進(jìn)行分組.傳統(tǒng)的產(chǎn)狀分組方法多是人為的由玫瑰花圖和極點(diǎn)圖進(jìn)行直觀判斷，主觀性很大，無法準(zhǔn)確定量地給出客觀性分組結(jié)果.

為了避免這種主觀性，在1976年，Shanley等[2]首次提出了一種客觀的聚類分組方法，該方法有著嚴(yán)格的數(shù)學(xué)理論依據(jù)，可以得到較為理想的結(jié)構(gòu)面分組結(jié)果，然而在尋找密度點(diǎn)時小球半徑的合理確定問題一直未能得到解決.后來，Harrison等[3]和Hammah等[4]在總結(jié)前人的研究成果后又提出并發(fā)展了模糊C均值(FCM)聚類算法在結(jié)構(gòu)面分組中的應(yīng)用，均取得了不錯的進(jìn)展.但上述方法本質(zhì)上都屬于局部尋優(yōu)的聚類方法，當(dāng)分組邊界不明確時，極易陷入局部極小值.隨著計算機(jī)技術(shù)的發(fā)展，智能算法等機(jī)器學(xué)習(xí)方法迅速興起，越來越多的學(xué)者將智能算法應(yīng)用于結(jié)構(gòu)面的分組中.李寧等[5]將改進(jìn)遺傳算法引入支持向量機(jī)的分類中，以此對結(jié)構(gòu)面進(jìn)行分組，分組結(jié)果較為客觀，但由于遺傳算法過程較為復(fù)雜，因此該方法應(yīng)用時存在一定的缺陷.后來，Li等[6]、王述紅等[7]、Li等[8]分別將蟻群算法、魚群算法、粒子群算法引入巖體結(jié)構(gòu)面分組的K-means聚類算法中，均得到較為滿意的分組結(jié)果，但是上述算法的復(fù)雜性使得這些方法的分組效率較低，不利于實(shí)際工程應(yīng)用.

鑒于此，本文提出了一種新型的融合模擬退火算法與K-means聚類(simulated annealing algorithm and K-means clustering，SAK)的結(jié)構(gòu)面產(chǎn)狀優(yōu)勢分組方法.該方法通過對K-means聚類算法進(jìn)行優(yōu)化，克服了K-means算法對初始值敏感，從而影響聚類結(jié)果的缺陷[9]，該算法簡單易實(shí)現(xiàn)，最終可搜索到全局最優(yōu)的結(jié)構(gòu)面分組結(jié)果.通過對計算機(jī)模擬生成的結(jié)構(gòu)面數(shù)據(jù)及現(xiàn)場實(shí)測結(jié)構(gòu)面數(shù)據(jù)進(jìn)行分析，并與已有方法進(jìn)行對比，驗(yàn)證了該方法的合理性、高效性和工程實(shí)用性.

1 巖體結(jié)構(gòu)面分組數(shù)學(xué)模型的建立

1.1 結(jié)構(gòu)面產(chǎn)狀的空間表示法

對結(jié)構(gòu)面產(chǎn)狀進(jìn)行分組，先要對結(jié)構(gòu)面傾向α和傾角β數(shù)據(jù)進(jìn)行歸一化處理，常用的方法是將結(jié)構(gòu)面視作無厚度的無限延伸平面，然后用其單位法向量來表示[10].由數(shù)學(xué)知識可知，這可以由空間中的單位球體表示，其在各坐標(biāo)軸上的分量如下：

(1)

則結(jié)構(gòu)面的單位法向量坐標(biāo)可表示為

p=(n1,n2,n3).

(2)

1.2 結(jié)構(gòu)面產(chǎn)狀數(shù)據(jù)之間的相似性度量

為了避免傾向相差約180°的兩組高陡傾角結(jié)構(gòu)面分組時出錯的情況，本文采用任意兩結(jié)構(gòu)面p1和p2之間所夾的銳角γ的正弦值作為兩結(jié)構(gòu)面之間產(chǎn)狀的相似性度量[11]，即

γ=arccos|p1·p2T|.

(3)

則兩結(jié)構(gòu)面單位法向量之間的距離為

(4)

1.3 結(jié)構(gòu)面分組的目標(biāo)函數(shù)

假設(shè)有n個結(jié)構(gòu)面Fi(i=1,2,…,n)，其單位法向量分別為pi(i=1,2,…,n)，可劃分為k組，每組聚類中心為cj(j=1,2,…,k)，定義uij為第i個結(jié)構(gòu)面屬于第j個分組的隸屬度，定義目標(biāo)函數(shù)J為所有結(jié)構(gòu)面單位法向量p與各分組中心c之間的距離和(即結(jié)構(gòu)面分組的總類間離散度)，即

(5)

(6)

式中：pi分別為第i個結(jié)構(gòu)面的單位法向量坐標(biāo);cj和cs為各分組中心的法向量坐標(biāo);m為權(quán)值分配系數(shù)，一般取1～2，本文取m=2;d(pi,cj)和d(pi,cs)分別為第i個結(jié)構(gòu)面到第j個和第s個分組中心之間的距離，采用式(4)計算.由式(6)可知：結(jié)構(gòu)面參數(shù)到各結(jié)構(gòu)面分組中心的距離越小，聚類的誤差越小，分組的結(jié)果也更精確.

2 基于SAK的巖體結(jié)構(gòu)面產(chǎn)狀分組

2.1 算法基本原理

模擬退火算法[12](SA)在20世紀(jì)80年代由Metropolis首次提出的一種啟發(fā)式隨機(jī)搜索算法，其基本原理來源于固體的退火過程.在SA算法中，目標(biāo)函數(shù)由內(nèi)能模擬，控制參數(shù)為溫度和溫度冷卻參數(shù)，通過對初始解重復(fù)執(zhí)行“擾動產(chǎn)生新解—計算目標(biāo)函數(shù)差—由Metropolis準(zhǔn)則判斷是否接受新解”的過程，逐步進(jìn)行優(yōu)化，算法迭代終止時的當(dāng)前解即為近似的最優(yōu)解.該算法具有漸近收斂性和并行性，已在理論上被證明是一種以概率1收斂于全局最優(yōu)解的全局優(yōu)化算法[13].

本文將其引入結(jié)構(gòu)面分組中，組成全局尋優(yōu)能力強(qiáng)的SAK算法，通過對K-means算法聚類結(jié)果進(jìn)行優(yōu)化，從而有效改進(jìn)K-means算法易受初始聚類中心影響的缺陷.

2.2 SAK算法的實(shí)現(xiàn)過程

在SAK算法中，初始溫度t0和溫度冷卻參數(shù)q的選取是非常重要的，它們對算法的收斂速度和全局最優(yōu)性有很大的影響.當(dāng)q取值較大時，算法的收斂速度大大降低;當(dāng)q取值較小時，溫度下降過快，此時則不易獲得全局最優(yōu)解.為了使最初產(chǎn)生的初始解能夠按照Metropolis準(zhǔn)則被接受，算法一開始就應(yīng)達(dá)到準(zhǔn)平衡狀態(tài)，因此選取初始溫度為t0時的聚類結(jié)果t0=J0作為算法的初始解，溫度冷卻參數(shù)由優(yōu)化效果適當(dāng)選取，降溫過程如下:

t(b+1)=t(b)q

(7)

式中：t(b)為當(dāng)前循環(huán)次數(shù)下的溫度值;t(b+1)為進(jìn)一步循環(huán)后的溫度值;q為冷卻參數(shù)，略小于1.

另外，SAK算法中新解的產(chǎn)生是通過對當(dāng)前解隨機(jī)擾動得到，擾動公式如下:

r=fix[rand()×n+1]

(8)

h=fix[rand()×k+1]

(9)

(10)

基于SAK算法的結(jié)構(gòu)面優(yōu)勢分組流程如下：

1) 首先對待分組的結(jié)構(gòu)面數(shù)據(jù)集進(jìn)行歸一化處理并執(zhí)行K-means聚類，將聚類分組的結(jié)果作為初始解s，由式(6)計算出初始目標(biāo)函數(shù)值Js；

2) 初始化溫度t0=Js、最大迭代次數(shù)L、溫度冷卻參數(shù)q和每個溫度t下的循環(huán)次數(shù)l(即Metropolis鏈長)；

3) 通過隨機(jī)擾動產(chǎn)生新解s′，這代表隨機(jī)產(chǎn)生了一個新的結(jié)構(gòu)面分組方式，計算此時對應(yīng)的目標(biāo)函數(shù)值Js′；

4) 判斷此時的目標(biāo)函數(shù)值Js′是否為最優(yōu)解，如果是則保存此時的結(jié)構(gòu)面分組方式為最優(yōu)分組、Js′為最優(yōu)目標(biāo)函數(shù)值，否則執(zhí)行步驟6)；

5) 計算目標(biāo)函數(shù)差△J=Js′-Js，并判斷其是否小于0，若小于0則接受新解s′作為下一個當(dāng)前解，否則按照Metropolis準(zhǔn)則(即以概率exp(-△J/t))接受新解s′；

6) 判斷當(dāng)前迭代次數(shù)是否達(dá)到最大迭代次數(shù),若是，則執(zhí)行步驟7)，若否，則繼續(xù)執(zhí)行步驟3)～5)；

7) 判斷是否達(dá)到終止溫度，若是，則算法結(jié)束，輸出當(dāng)前聚類劃分結(jié)果作為結(jié)構(gòu)面優(yōu)勢組分類結(jié)果；若否，則降低溫度，繼續(xù)執(zhí)行步驟3)～6).

2.3 聚類有效性評價

為了避免對最優(yōu)分組結(jié)果判斷的單一性，本文采用模糊分類系數(shù)F和分類熵指標(biāo)H進(jìn)行聚類有效性評價，這兩類指標(biāo)也是最常用的衡量聚類有效性的函數(shù)[5]，其計算公式分別為

(11)

(12)

式中：n為結(jié)構(gòu)面數(shù)據(jù)集的個數(shù);uij為第i個結(jié)構(gòu)面屬于第j個分組的隸屬度，由式(5)計算所得;a為對數(shù)的底數(shù)，a∈(1，+∞)，規(guī)定當(dāng)uij=0時uijloga(uij)=0，本文取a=10.當(dāng)F越大，H越小，表明分類的模糊度越小，聚類效果越好.因此，對于同一方法，當(dāng)F較大、H較小時，結(jié)構(gòu)面分組結(jié)果較好;而對于不同方法，F(xiàn)相對較大，H相對較小的方法為較好的方法.

3 算法準(zhǔn)確性驗(yàn)證

為了驗(yàn)證SAK算法應(yīng)用于結(jié)構(gòu)面產(chǎn)狀優(yōu)勢分組中的準(zhǔn)確性，采用計算機(jī)隨機(jī)模擬生成了3組界限并不明顯的120個結(jié)構(gòu)面數(shù)據(jù)(服從正態(tài)分布).表1為這三組結(jié)構(gòu)面數(shù)據(jù)的詳細(xì)參數(shù)及兩種算法聚類后分組中心的對比結(jié)果;表2為聚類有效性評價結(jié)果;圖1為這3組結(jié)構(gòu)面數(shù)據(jù)的極點(diǎn)及密度等值線圖;圖2a和圖2b分別為采用K-means算法和SAK算法分組的結(jié)果.SAK算法相關(guān)參數(shù)設(shè)置為：初始溫度tbegin=20 ℃，終止溫度tend=0.1 ℃，溫度冷卻參數(shù)q=0.92，每個溫度t下的循環(huán)次數(shù)(即Metropolis鏈長)l=400，最大迭代次數(shù)L=500.

表1 結(jié)構(gòu)面數(shù)據(jù)參數(shù)及分組中心Table 1 Discontinuity data parameters and grouping centers

表2 隨機(jī)數(shù)據(jù)聚類有效性檢驗(yàn)Table 2 Clustering validity test of random data

由表1中分組結(jié)果與已知結(jié)果的對比中可以看出，SAK算法的分組中心與已知中心更為接近；表2顯示了兩種算法的聚類有效性結(jié)果，當(dāng)分組數(shù)為2和3時，結(jié)果較為接近，但綜合比較F和H兩指標(biāo)后可知最優(yōu)分組數(shù)為3.此時，從圖1、圖2a和圖2b的對比中，也可以很明顯地看出SAK算法準(zhǔn)確性更高，分組結(jié)果更符合圖1中實(shí)際的情況.

4 工程實(shí)例

重慶市三環(huán)高速公路合川至長壽段興隆隧道位于重慶市渝北區(qū)木耳鎮(zhèn)，隧址區(qū)屬構(gòu)造侵蝕丘陵地貌，隧道大體沿垂直構(gòu)造線方向布設(shè)，與巖層走向呈大角度相交，穿越地層主要為侏羅系上沙溪廟組地層，地層分布連續(xù).圍巖巖性主要為侏羅系中統(tǒng)上沙溪廟組泥巖和砂巖，中風(fēng)化巖體較完整，發(fā)育高傾角的構(gòu)造裂隙.隧址區(qū)無活動性斷裂、泥石流、滑坡等不良地質(zhì)現(xiàn)象，地下水類型主要為第四系松散土層孔隙水及基巖裂隙水.

本文以隧道洞口段現(xiàn)場實(shí)測的118個基巖露頭結(jié)構(gòu)面產(chǎn)狀數(shù)據(jù)為例，進(jìn)一步驗(yàn)證所提SAK模型在結(jié)構(gòu)面分組中的合理性和工程實(shí)用性.赤平投影下結(jié)構(gòu)面產(chǎn)狀的極點(diǎn)及密度等值線圖如圖3所示.

由圖3可看出，各組結(jié)構(gòu)面之間的界限并不明顯，根據(jù)該圖可大致判斷出分組數(shù)為2～4，因此對結(jié)構(gòu)面數(shù)據(jù)分別采用K-means和SAK算法進(jìn)行聚類分析，并與文獻(xiàn)[8]中的KPSO算法進(jìn)行對比，利用F和H兩指標(biāo)進(jìn)行聚類有效性評價，對比結(jié)果如表3所示.由表中結(jié)果可以看出，當(dāng)分組數(shù)為2時，聚類效果最好，此時三種算法聚類結(jié)果分別如圖4a、圖4b和圖4c所示.K-means算法、SAK算法和文獻(xiàn)[8]中KPSO算法迭代速度對比結(jié)果如圖5所示；結(jié)構(gòu)面優(yōu)勢產(chǎn)狀統(tǒng)計結(jié)果見表4.

表4 結(jié)構(gòu)面優(yōu)勢產(chǎn)狀分組結(jié)果Table 4 Grouping results of discontinuity dominant orientation

表3 實(shí)測數(shù)據(jù)聚類有效性檢驗(yàn)Table 3 Clustering validity test of measured data

對比圖3、圖4a和圖4b，不難看出，SAK算法的結(jié)構(gòu)面優(yōu)勢分組結(jié)果更符合密度等值線圖所反映的分組區(qū)域，且與圖4c中KPSO算法的分組結(jié)果幾乎一致.原因主要是K-means算法通過隨機(jī)生成初始聚類中心進(jìn)行聚類，因此其聚類結(jié)果容易不理想，而采用全局尋優(yōu)能力強(qiáng)的模擬退火算法可以對K-means算法聚類結(jié)果進(jìn)行優(yōu)化，從而改善原始K-means聚類算法由于初始聚類中心選擇不當(dāng)而使聚類效果不理想的缺陷.根據(jù)表3中聚類有效性檢驗(yàn)結(jié)果的對比也可以看出，無論分組數(shù)為幾組，SAK算法及KPSO算法的聚類的結(jié)果都要大大優(yōu)于K-means算法的聚類結(jié)果，并且本文所提方法略優(yōu)于文獻(xiàn)[8]中的KPSO算法，進(jìn)一步驗(yàn)證了所提算法具有較高的精度.

另外，從圖5中可以很明顯地看出，同樣都是對K-means算法進(jìn)行優(yōu)化，文獻(xiàn)[8]中算法整體迭代速度大大降低，而本文提出的SAK算法的迭代速度則降低較少，本文方法的迭代速度要大大優(yōu)于文獻(xiàn)[8]中KPSO算法(大約提高50%)，表明了算法的高效性.因此，基于模擬退火算法K-means聚類(SAK)的結(jié)構(gòu)面優(yōu)勢分組方法更適合應(yīng)用于結(jié)構(gòu)面的分組中，具有較強(qiáng)的合理性和工程實(shí)用性.

5 結(jié) 語

本文將模擬退火算法引入了結(jié)構(gòu)面的優(yōu)勢分組中，提出了一種新型的基于模擬退火算法K-means聚類(SAK)的結(jié)構(gòu)面優(yōu)勢分組算法.該算法通過逐次迭代后進(jìn)行最優(yōu)解的精確搜索，最終搜索到全局最優(yōu)的結(jié)構(gòu)面分組結(jié)果，有效克服了K-means算法對初始值敏感，從而影響聚類效果的缺陷，避免了人為劃定分組方式的主觀性.結(jié)合計算機(jī)模擬生成的結(jié)構(gòu)面數(shù)據(jù)及重慶市三環(huán)高速公路興隆隧道洞口段現(xiàn)場實(shí)測的結(jié)構(gòu)面數(shù)據(jù)，將該算法與文獻(xiàn)[8]中提出的KPSO算法進(jìn)行對比，證明了SAK算法在聚類準(zhǔn)確性、聚類精度及迭代速度上均較優(yōu)，可以得到較為合理的結(jié)構(gòu)面分組結(jié)果，有一定的推廣應(yīng)用價值.