基于遺傳算法的巖體結構面產狀的自動聚類方法研究

崔學杰， 陳江軍， 向進洋， 關義濤， 胡小慶

(1.湖北省地質局 水文地質工程地質大隊，湖北 荊州 434020； 2.資源與生態環境地質湖北省重點實驗室，湖北 武漢 430034)

地下水封儲油(氣)洞庫是目前國際上石油、液化氣等能源的主要儲存方式[1]。結構面對巖體結構起控制作用，對地下水封洞庫圍巖的穩定性有重要影響[2-4]。由于地下水封洞庫埋深大，結構面信息的獲取目前主要依靠鉆孔攝像(或電視)技術[5-6]。

結構面數據分析的首要任務是對結構面進行分組，傳統方法是根據結構面走向玫瑰花圖或極點密度等值線圖來進行分組，該方法簡單、直觀，但結果依賴于分析者的個人經驗，缺乏客觀性，同時在數據量較大時不可用[7]。Shanley et al.[8]和Mahtab et al.[9]率先提出了用于結構面產狀分組的聚類方法。Hammah et al.[10-12]采用模糊C均值(FuzzyC-means，FCM)方法進行結構面組數的劃分，并討論了常用的距離量度和聚類有效性指標。FCM是K-means算法結合模糊理論的改進版本，但二者都使用Picard迭代算法來最小化目標函數[13]，因而易收斂于局部最優解，且聚類結果受初始聚類中心影響。Jimenez-Rodriguez et al.[14]嘗試通過對相似度矩陣進行特征分解，從而將原始數據轉換到區分度更顯著的空間，提出了用于結構面分組的譜聚類方法,該方法避免了選擇初始聚類中心的困難，但仍不能完全保證收斂于全局最優解。Xu et al.[15]將問題變量從解空間轉換至混沌空間，開發了基于多尺度混沌優化算法的FCM算法，利用混沌變量的隨機性、有序性和遍歷性搜索從而獲得全局最優解。Li et al.[16]提出了基于粒子群優化算法的K-means算法(KPSO)，該方法的性能略優于Xu et al.[15]的方法。類似地，Song et al.[17]提出了基于量子粒子群的FCM算法(QPSO)。上述方法盡管改善了一些缺陷，但均需要指定結構面組數作為輸入值，否則需要進行大量試算。

遺傳算法(Genetic Algorithm，GA)是一種高效、穩健的隨機優化技術，理論和經驗都表明該方法將提供全局最優解[18-21]。遺傳算法在土木工程和地質工程領域獲得了大量的應用，其有效性獲得廣泛地驗證[22-26]。基于上述認識，本文對傳統遺傳算法做了改進，結合FCM算法，提出了一種用于結構面產狀分組的聚類算法，并將其用于分析某地下水封洞庫結構面產狀數據，從而驗證該方法的有效性。

1 結構面產狀的表示方法

結構面通常被假設為空間中一個幾何平面，其產狀使用傾向α(0°≤α≤360°)和傾角β(0°≤β≤90°)表示，也可用結構面單位法向量或單位矩陣表示。如圖1所示建立空間直角坐標系：x1軸正向指向正北，x2軸正向指向正東，x3軸正向指向上。圖1中的結構面j可以由公式(1)表示。同時，公式(1)也可以表示定義在單位球體上半球面上的點的集合，因此參考文獻[10]也將產狀數據歸類為球狀數據。

圖1 結構面產狀表征示意圖

(1)

2 FCM算法

FCM算法是常用的動態聚類分析算法之一，這里作簡要介紹。假設實測了n個結構面產狀，用集合X表示，即X={x1,x2,…,xn}，其中每一個元素都可以用一個3維向量表示，如公式(1)中的xj。FCM算法的目標是依據數據的相似性或非相似性，根據隸屬度將數據集合X劃分為K個子集(或組)，如{C1,C2,…,CK}。FCM算法通常包含以下步驟[10,13]：

(1) 從數據集合X={x1,x2,…,xn}中隨機選擇K個數據作為初始聚類中心，即z1、z2、…、zK；

(2) 計算模糊隸屬度矩陣U(X)=[ukj]K×n，其中:

(2)

式中:i,k=1,2,…,K；j=1,2,…,n；m為模糊指數，一般取m=2[13]；D(zk,xj)表示xj與Ck組的聚類中心zk非相似性(或距離)。對于結構面產狀等球狀數據，參考文獻[11]建議使用結構面法向量夾角或夾角的函數表示結構面產狀的距離，即

(3)

式中:T表示矩陣轉置。

(3) 根據模糊隸屬度矩陣，計算新的聚類中心z′1、z′2、…、z′K。對于結構面數據等球狀數據，參考文獻[12]建議使用特征分析方法計算聚類中心，對于第k組結構面數據，需要先計算模糊產狀矩陣Sk=[Sk,uv]3×3，其中:

(4)

式中:u,v=1,2,3。計算矩陣Sk的特征值λk1、λk2、λk3和歸一化特征向量ηk1、ηk2、ηk3，如果λk1<λk2<λk3,則最大特征值λk3對應的特征向量ηk3為第k組的聚類中心，即:

z′K=ηk3

(5)

(4) 計算目標函數。計算公式如下：

(6)

(5) 重復步驟(2)—(4)，直到目標函數Jm達到最小值，或U(X)不再發生改變。

3 基于遺傳算法的FCM算法

遺傳算法是由Holland受生物進化和自然遺傳理論啟發提出[20]，適用于空間大、復雜且多峰的搜索和優化問題。該算法是一種全局優化搜索算法，具有簡單通用、魯棒性強等優勢，其核心思想是將問題可能的解編碼為字符串格式(稱為染色體或個體)，并將多個個體組成一個種群；種群遵循自然選擇、適者生存的進化原則，在遺傳算子引導下繁衍后代，產生新的種群；經過若干代后，算法收斂于最優個體。該算法主要包括編碼方法與種群初始化、適應度函數和遺傳算子三個部分。

3.1 編碼方法與種群初始化

傳統的編碼方法通常以聚類中心作為遺傳信息，將其編碼為染色體。對于結構面產狀數據的分組問題，結構面組的平均產狀即為聚類中心。以種群中第i個染色體為例，如果通過現有資料能夠確定有K組結構面，則染色體中將包含K個聚類中心，長度為3K(圖2-a)。在無法預先確定結構面組數時，需要將K設置為不同的整數值進行試算。

為了克服上述缺陷，這里介紹一種新的編碼方法。首先需要預估一個分組的上限值Km(圖2-b)。對于結構面產狀數據，為保守起見，可令Km=10。與傳統方法不同，新方法為每一個聚類中心z附加一個隨機實數變量p(0 ≤p≤ 1)。為了便于說明，這里借助于遺傳學術語，將實數變量p稱為啟動子，而對應的聚類中心稱為外顯子，二者組成一個基因組，如pi1和zi1為染色體i的第1個基因組。啟動子的作用在于控制對應的外顯子是否表達，即聚類中心是否參與聚類過程。為了便于實現，這里采用如下簡單規則：當啟動子>0.5時，對應的外顯子表達。使用該編碼方法時，盡管種群中染色體長度相同，但由于染色體中啟動子為隨機變量，有效的聚類中心數量將是動態變化的，這樣便可以實現結構面組數為動態的聚類分析。

圖2 兩種不同編碼方法示意圖

種群初始化如下所述：對于任一條染色體，首先從均勻分布U[0,1]中隨機產生Km個實數初始化啟動子，再從結構面產狀數據中隨機選擇Km個產狀數據初始化外顯子；如果種群中包含P個染色體，上述步驟將重復P次以完成種群初始化。注意，初始化完成后，應檢查初始種群中的每一條染色體是否合格,若染色體中>0.5的啟動子數量<2，則該染色體不合格，需重新初始化。

3.2 適應度函數

根據自然選擇、適者生存的原則，個體的適應度越高，其繁衍后代并構成新種群的可能性越大。對于聚類問題，常用的做法是將聚類有效性指標作為適應度函數，定量評價染色體所對應的聚類結果的優劣，即適應度值。適用于FCM算法的模糊聚類有效性指標包括目標函數Jm、劃分系數PC、劃分熵指標PE、Xie-Beni指標、Fukuyama-Sugeno指標VFS等[27-28]。這里仍選擇應用較為廣泛的Xie-Beni指標，其定義如下[29]：

(7)

式中：K′表示有效聚類中心數量；σ表征組內緊湊度，值越小表示聚類效果越好；sep為組間距離的最小值，表征組間分離度，值越大表示聚類效果越好；其他符號與前述一致。

由公式(7)可知，當聚類效果最優時，XB指標取最小值。因此，應將適應度函數定義為XB指標的倒數，即

(8)

對于當前種群中的任一染色體，需要按以下步驟計算其適應度值：①提取染色體中有效聚類中心，并確定K′值；②按照FCM算法對結構面產狀數據進行聚類；③將最終的隸屬度矩陣和聚類中心代入公式(8)計算適應度值。

3.3 遺傳算子

選擇、交叉和變異是三種基本遺傳算子，其共同引導種群的進化。

3.3.1選擇算子

選擇算子模擬達爾文進化論中的自然選擇和適者生存過程。適應度越高的個體，有更高的概率被選擇，從而繁殖后代，將遺傳信息傳遞給后代。換言之，只有被選擇的個體才能參與交叉和變異過程。常用的方法有賭輪盤選擇法、錦標賽選擇法、隨機遍及取樣法等[30]。本文選擇易于實現的賭輪盤選擇法。

3.3.2交叉算子

交叉算子交換父代個體的遺傳信息，產生與父代不同的后代。單點交叉、兩點交叉、多點交叉和均勻交叉是常用的交叉方法[20,30]。交叉雖然能夠增加種群的多樣性，但同時也存在破壞優勢基因組合的風險。權衡利弊，這里選擇兩點交叉法，且設置交叉概率為μc(0<μc<1)。在交叉過程中，每一個聚類中心都被視為一個不可分割的三維向量，即交叉點只能位于啟動子和外顯子之間。舉例說明(圖3)，如果個體1、個體2被選為父代個體，且交叉點位置如圖3-a中箭頭所示，則執行交叉后的子代個體如圖3-b所示。交叉完成后，個體3與父代個體1相比較，其獲得了兩個新的外顯子z21、z22，但外顯子z21是否表達取決于啟動子p11，而原始外顯子z13是否表達取決于新的啟動子p23。子代個體4情況相似。兩點交叉法在提高基因重組概率的同時，使得對基因組的破壞仍相對較小。

圖3 交叉算子示意圖

3.3.3變異算子

變異算子通過產生新的基因來保持種群多樣性，即引導探索新的搜索空間，防止受限于當前局部空間。對于所使用新的編碼方法，啟動子和外顯子將采用不同的變異方法，但變異概率均設置為μm(0<μm<1)。

對于染色體中的任一啟動子pik，變異后的啟動子p′ik按下式計算：

p′ik=pik±0.5

(9)

式中：±號取正和取負的概率均為0.5。變異后的啟動子若<0，則設置為0；若>1，則設置為1。

聚類中心zik的變異是通過向量的旋轉運算實現的。如圖4所示，將zik假想為一虛擬結構面，αik、βik分別表示其傾向和傾角，則：

圖4 變異算子示意圖

(10)

假設新的聚類中心為z′ik，與zik之間的夾角為δ(圖4)。顯而易見，新的聚類中心z′ik是將zik繞原點o旋轉角度δ而產生，而旋轉方位角θ則是隨機的，將從均勻分布U[0,2π)中隨機產生。變異的數學表達如以下公式所示：

(11)

式中：R2(βik)、R3(αik)分別表示將向量繞x2軸、x3軸旋轉的旋轉矩陣。

(12)

(13)

3.4 收斂準則

本文使用精英模型[18,30]來保證最優個體不會因遺傳算子的隨機性而被遺漏。當種群中最優個體的適應度值無明顯變化時，認為算法收斂。

3.5 算法的實現

采用C++語言實現了上述算法。算法采用通用性較強的文本文件(.txt)作為輸入、輸出格式，結構面產狀作為輸入數據，輸出信息包含結構面組數、隸屬度矩陣、平均產狀等。

4 工程實例應用

該地下水封儲油(氣)洞庫位于浙江省東部沿海地帶，主要巖性為上侏羅統高塢組花崗斑巖，結構面比較發育。根據設計運行壓力已確定洞庫埋深，結構面的發育規律對主洞室軸向有重要影響，進而影響洞庫的平面布置和主洞室截面形狀的選擇。為了探明設計埋深區域內結構面發育情況，進行了初步的勘探工作。本文選用了操作豎井附近的1條垂直鉆孔，根據其聲波電視解譯結果(圖5)獲得了從地表至-172 m高程范圍內的386個結構面產狀數據，其等密度極點投影如圖6所示。

圖5 聲波鉆孔電視成像圖(91～96 m)

圖6 結構面等面積上半球極點投影圖

由圖6可知，陡傾角結構面的極點投影較為集中，構成了3個高密度區，按照經驗應分為2組結構面，一組傾向NE，一組傾向SE；中—緩傾角結構面數量少，極點投影比較分散，將這些結構面分為1組或2組都是合理的。因此，若采用目前常見的聚類算法，需要試算來確定結構面組數。為了驗證本文提出的基于遺傳算法的FCM算法的有效性，使用該方法對圖5所示結構面產狀數據進行分組，所用的參數如表1所示，聚類結果如圖7和表2所示。新方法將所獲得的結構面按產狀分為3組，其中中—緩傾角結構面被分為1組(J3)。為了驗證該方法的合理性，使用Li et al.[16]提出的KPSO法計算結構面組數不同情況下的Ray-Turi指標值，結果見表3。由結果可知，Ray-Turi指標在結構面組數為3時取得最小值，即最優組數應取3，與本文提出的方法結果一致。結構面的分組成果將為洞室的設計和圍巖穩定性分析提供參考。

表1 基于GA的FCM算法參數表

表2 使用基于GA的FCM算法的結構面聚類結果

圖7 聚類結果等面積上半球極點投影圖

表3 不同結構面組數情況下KPSO方法聚類結果的Ray-Turi指標值

5 結論

根據結構面產狀進行優勢分組是開展巖體穩定性分析工作的基礎。與傳統方法相比，聚類方法更加客觀且分組結果更加明確，FCM算法是常用的動態聚類分析算法之一，但也存在明顯缺陷。鑒于遺傳算法較強的穩健性和全局最優性，將其與FCM算法結合，提出了一種用于分析產狀數據的聚類方法，該方法的聚類結果不受初始聚類中心影響，總是能夠收斂于全局最優解；而且由于使用了新的染色體編碼方法，該方法無需輸入結構面組數，在難以確定結構面組數的情況下，能夠節約大量時間。與現有聚類方法相比，該方法參數較少、客觀性更強。將該方法用于某地下水封儲油(氣)洞庫結構面產狀數據的聚類分析，獲得了明確、合理的結構面分組，可見此方法具有良好的工程實用性。