雙重變異遺傳算法及其在臨界滑動面搜索中的應用

覃 偉

（重慶工程職業技術學院，重慶 402260）

在公路工程、建筑工程、礦山工程、水利工程等領域中，常會遇到大量的邊坡穩定性問題。確定邊坡臨界滑動面對邊坡穩定性分析及坍岸寬度預測等研究具有重要意義[1]。

邊坡臨界滑動面的搜索，本質上是一個優化問題，其具有變量多、約束條件復雜、局部極值點多、非凸性等特點，使得傳統的優化算法（如模式搜索法）在搜索臨界滑動面時，常常陷入局部極值點，給邊坡的穩定性分析帶來許多困難。

遺傳算法（genetic algorithm，GA）由 Holland 教授于 20世紀 70 年代提出，是一種模擬達爾文遺傳選擇和自然淘汰的生物進化過程的搜索尋優算法，具有思想簡單、易于編程實現、算法健壯等優點，同時還具有隱含并行性和全局搜索等顯著特性[2]。因此，一些學者將遺傳算法運用到邊坡臨界滑動面搜索研究中，取得了一定成效[3?8]。如，賀子光等[3]將單純形法和回溯機制引入GEP 中，提出了混合GEP 方法搜索邊坡的臨界滑動面；梁冠亭等[4]采用自適應遺傳優化算法搜索采用抗滑樁支護邊坡的臨界滑動面；朱劍鋒等[5]提出自適應禁忌變異遺傳算法搜索土釘墻臨界滑動面。但是目前的標準遺傳算法及其改進算法仍在不同程度上存在收斂速度較慢、容易進入早熟等缺陷，有待進一步改進和完善。

本文在傳統遺傳算法的基礎上，提出雙重變異策略對遺傳算法進行改進，形成雙重變異遺傳算法（DMGA），并結合簡化Bishop 法，編寫Python 程序，通過對土坡算例臨界滑動面的計算，進一步為邊坡臨界滑動面搜索提供有效方法。

1 臨界滑動面搜索問題的目標函數

求解邊坡的臨界滑動面，即是搜索使邊坡安全系數最小的滑動面。邊坡穩定性計算方法有很多，主要包括剛體極限平衡方法和數值模擬方法[9?10]。對于滑動面呈圓弧形的土質邊坡，由于簡化Bishop 法計算簡單、物理意義明確，在工程實踐中應用廣泛。因此本文采用簡化Bishop 法計算土質邊坡的穩定性。

1.1 簡化Bishop 法

簡化Bishop 法是一種適用于滑動面呈圓弧形的邊坡的穩定性分析方法。其中邊坡安全系數F的計算公式[11]為：

式中：ci—第i條塊滑動面上巖土體的黏聚力/kPa；

φi—第i條塊滑動面上巖土體的內摩擦角/（°）；

li—第i條塊滑動面的弧長/m；

Wi—第i條塊的單寬重量/（kN· m?1）；

αi—第i條塊滑動面傾角/（°），當滑動面傾向與滑動方向一致時，αi取正；當滑動面傾向與滑動方向相反時，αi取負。

當變量mi很小時，計算結果與實際情況不符。根據一些學者的意見，當mi?0.2時，不能采用簡化B ishop 法計算邊坡安全系數[11]。

1.2 目標函數的確定

邊坡臨界滑動面的位置與圓弧形滑動面的圓心半徑R、圓心橫坐標xo、圓心縱坐標yo有關，因此將這3 個量作為目標函數的求解變量。

并不是任意半徑及圓心位置的圓弧都能計算出邊坡的安全系數，如發生圓弧與邊坡不存在交點，或所采用的邊坡安全系數計算方法不適用等情況。因此，臨界滑動面搜索問題可表示為在一定約束條件下求解邊坡安全系數最小值，即求解下面約束問題：

運用外點罰函數法[12]將該約束問題轉化為無約束問題：

其中：

其中，當圓弧滑動面與坡體無交點時，F（x）不存在，則令其為一較大值100；σ為很大的正數，取值1 000。

2 雙重變異遺傳算法的設計

2.1 個體編碼與初始種群產生

本文對變量采用一種特殊的實數編碼，在編碼中設置符號基因位與數字基因位。符號基因位存放實數的符號信息，當實數為正或0 時，該基因位的值取1；實數為負時，該基因位的值取?1。數字基因位存放實數的數字部分，實數各位上的數字占一個基因位，且該位上的數字即為對應基因位的值。對每個變量xi均采用上述實數編碼，并將它們的編碼依次連接構成一個染色體。例如，若目標函數的解由3 個變量構成，其中變量x1為36.12，變量x2為?5.099，變量x3為72.562，則對其進行實數編碼如圖1所示。

圖1 實數編碼Fig.1 Real number coding

初始種群內的個體數為N，每個個體通過隨機數產生，個體染色體中的每個數字基因位隨機生成[0,9]區間的整數，每個符號基因位隨機生成1 或?1，若該變量不為負，則對應符號基因位僅產生1。

2.2 適應度值函數

遺傳算法在進化時，以個體的適應度值作為搜索依據。求解優化問題時，需要將目標函數轉化為適應度函數。由于目標函數G（x）為最小值問題，構造適應度值函數如下：

適應度函數值f（x）越大，邊坡的安全系數就越小。

2.3 雙重變異遺傳算法的進化過程

2.3.1 選擇操作

父代個體的適應度值確定后，按照適應度值進行個體的選擇。本文種群的選擇采用精英保留策略[13?14]，每一代種群進化完成后，首先從中找出適應度值最大的一個個體，直接讓其進入下一代，不再對其進行交叉、變異操作。剩余的個體，則采用轉盤式選擇方式進行選擇操作，轉盤式選擇的基本原理是根據每個個體的適應度值的比例來確定該個體的選擇概率[15]。

轉盤式選擇的計算方法：先根據個體的適應度值fi按式（8）計算個體的選擇概率pi，并令p0=0，然后生成[0,1]區間的隨機數r，如果r滿足式（7），則選擇個體j。

其中：

2.3.2 自適應交叉操作

在遺傳算法中，通過交叉操作既可以使種群多樣化，又可以引導搜索方向，促進優秀個體的產生。本文采取單變量交叉方式，即從種群中選出任意2 個個體（個體1、個體2）進行交叉，當產生的[0,1）區間的隨機數小于交叉概率Pc時，隨機選擇一個變量xi，使個體1 與個體2 中變量xi所對應的基因位的值進行交換。

交叉概率是影響遺傳算法性能的關鍵因素之一。盧明奇等[16]為了解決傳統的自適應遺傳算法在進化初期常存在停滯現象、進化后期易陷入局部極值點的問題，提出了一種改進的自適應交叉概率。但該方法參數較多，增加了參數設置的難度。本文在考慮適應度值及進化代數的情況下，提出一種復合自適應交叉概率，使交叉概率Pc隨個體適應度值及進化代數進行調整，Pc的計算公式如下：

式中：fmax—當前種群中的最大適應度值；

f′—待交叉個體中較大的適應度值；

—當前種群的平均適應度值；

T—總的進化代數；

t—當前進化代數。

2.3.3 雙重變異操作

在傳統的遺傳算法中，變異操作本質上是隨機搜索，并不能保證產生改進的后代[17]，進而常導致遺傳算法的收斂速度較慢。為此，本文借鑒模式搜索法中探測移動的思想，提出探測變異操作，并與傳統的直接變異操作相結合，構成雙重變異策略對遺傳算法進行改進，使遺傳算法的收斂速度及全局搜索能力得到提升。雙重變異操作包括探測變異操作與直接變異操作。其中，探測變異操作將提升算法的局部尋優能力，直接變異操作將使算法具備全局尋優能力。

為了敘述方便，將采用雙重變異策略的遺傳算法稱為雙重變異遺傳算法（DMGA），將僅采用探測變異策略的遺傳算法稱為探測變異遺傳算法（DEMGA），將僅采用直接變異策略的遺傳算法稱為直接變異遺傳算法（DIMGA）。

（1）探測變異操作

探測變異的具體操作是在隨機選擇1 個基因位后分5 種情況進行。

情況1：若選擇的是符號位變異，則直接進行變異探測，并計算適應度值，判斷是否發生變異（圖2）。

圖2 探測變異操作情況1 示意圖Fig.2 Schematic diagram of the first case of detecting mutation operation

情況2：若選擇的基因位為變量xi的第1 個數字位，則直接進行變異探測，并計算適應度值，判斷是否發生變異。該變異探測分3 種子情況（圖3）。子情況1，若該位的值為0，則變異為[1,9]區間的隨機整數；子情況2，若該位的值為9，則變異為[0,8]區間的隨機整數；子情況3，若該位的值為[1,8]區間的整數m，則變異分兩個方向進行探測，一個方向是變異為[m+1,9]區間的隨機整數，另一個方向是變異為[0,m?1]區間的隨機整數，然后將適應度值較大者與變異前適應度值進行比較，判斷是否變異。

圖3 探測變異操作情況2 示意圖Fig.3 Schematic diagram of the second case of detecting mutation operation

情況3：若選擇的基因位位于變量xi第1 個數字位之后，且該位的值位于區間[1,8]，則直接進行變異探測（同情況2 的子情況3），并計算適應度值，判斷是否發生變異。

情況4：若選擇的基因位位于變量xi第1 個數字位之后，且該位的值為0，則該變異探測分2 種子情況（圖4）。子情況1，當其前一位的值等于0 時，則該位變異為[1,9]區間的隨機整數；子情況2，當其前一位的值大于0 時，則變異分兩個方向進行探測，一個方向是該位變異為[1,9]區間的隨機整數m，另一個方向是該位值也變異為m，且前一位的值作減1 變異，然后將適應度值較大者與變異前適應度值進行比較，判斷是否變異。

圖4 探測變異操作情況4 示意圖Fig.4 Schematic diagram of the fourth case of detecting mutation operation

情況5：若選擇的基因位位于變量xi第1 個數字位之后，且該位的值為9，則該變異探測分2 種子情況（圖5）。子情況1，當其前一位的值等于9 時，則該位變異為[0,8]區間的隨機整數；子情況2，當其前一位的值小于9 時，則變異分兩個方向進行探測，一個方向是該位變異為[0,8]區間的隨機整數m，另一個方向是該位值也變異為m，且前一位的值作加1 變異，然后將適應度值較大者與變異前適應度值進行比較，判斷是否變異。

圖5 探測變異操作情況5 示意圖Fig.5 Schematic diagram of the fifth case of detecting mutation operation

（2）直接變異操作

直接變異操作是指將待變異個體的每個變量xi所對應的染色體等分為前后兩段，在每一段中隨機選擇1 個基因位進行直接變異，變異后的適應度值不控制變異的發生。直接變異能夠提供種群的多樣性，避免早熟的產生。

“昆南”陰平聲字“他”的唱調（《南西廂·佳期》【十二紅】“愛他兩意和諧”，687），其中的即為兩節型過腔。雖然這個過腔的音樂材料都相同，都來自于本唱調音階的級音，但由此組成的樂匯或句型可以分作兩節，為第一節級音性過腔，為第二節級音性過腔。這個過腔即是由同一種音樂材料組成的兩節型過腔。

（3）雙重變異操作

當待變異個體適應度值f接近當前最佳個體時，希望待變異個體能夠具有局部尋優能力，同時也具有跳出局部極值點，搜索全局最優的能力；當待變異個體適應度值f距離當前最佳個體較遠時，希望個體向著適應度值增加的方向變異。將探測變異與直接變異結合在一起形成雙重變異操作，便可以實現上述目標。

雙重變異操作的實施策略：

（4）變異概率

變異操作是在一定的概率下發生的，當產生的[0,1）區間的隨機數小于變異概率Pm時，將會發生變異。變異概率的大小決定著新個體產生的機會及計算量的大小。較小的變異概率，產生新個體的機會較小；較大的變異概率，遺傳算法的計算量較大，且趨于純粹的隨機搜索算法。

為了能達到好的尋優效果，在進化的早期期望能有較大的變異概率，以增加搜索的廣度；在進化的晚期期望能有較小的變異概率，以提高進化的速度。同時，也希望當待變異個體適應度值接近當前最佳個體時，減小變異概率，使較好的個體得到保留；當待變異個體適應度值距離當前最佳個體較遠時，有較大的變異概率，使較差的個體容易被淘汰。

基于上述原因，本文在綜合考慮待變異個體的適應度值f及進化代數對變異概率影響的情況下，提出一種復合自適應變異概率Pm，其計算公式如下：

其中，a為區間（0,1）內的常數，a越大，變異概率Pm越大；k為區間[0,1]內的常數，k越大，變異概率Pm也越大。反映了待變異個體的適應度值距離當前種群中的最大適應度值的相對距離；當一定時，Pm將隨t的增大而減小，即進化早期變異概率較大，使算法有較大的搜索廣度；當t一定時，越大，Pm就越小，即待變異個體適應度值越接近當前最佳個體，較好的個體就更容易得到保留。

2.4 雙重變異遺傳算法的計算步驟

步驟一：隨機生成實數編碼（圖1）的初始種群，初始種群中的個體數為N，且N為偶數。

步驟二：將個體染色體解碼為實數，根據式（6）計算初始種群中每個個體的適應度值，并令當前進化代數t=0。

步驟三：根據式（8）計算選擇概率，并令j=0。

步驟四：如果j=0，將種群中適應度值最大的一個個體直接選出，不進行后續的交叉及變異操作，直接復制到下一代，然后根據式（7）采用轉盤式選擇方式選擇1 個個體，進行步驟六；如果j>0，根據式（7）采用轉盤式選擇方式選擇2 個個體，進行步驟五。

步驟五：將選出的2 個個體進行單變量交叉操作，進行步驟六。交叉操作方法：當產生的[0,1）區間的隨機數小于交叉概率Pc時，隨機選擇一個變量xi，使個體1 中變量xi所對應的基因位的值與個體2 的進行交換。其中，交叉概率Pc根據式（9）計算得到。

步驟六：對個體進行雙重變異操作，進行步驟七。雙重變異操作方法：當時，產生一個[0,1]的隨機數 γ，若γ>0.5，采用探測變異；若γ ?0.5，采用直接變異。當時，采用探測變異。探測變異操作、直接變異操作按2.3.3 節所述方法進行。變異概率Pm按式（10）計算。

步驟七：如果（j+1）×2

步驟八：將新生成的種群個體的染色體解碼為實數，根據式（6）計算適應度值。如果t==T，計算結束；如果t

說明：本文所述的直接變異遺傳算法、探測變異遺傳算法與雙重變異遺傳算法的計算步驟的不同之處主要體現在步驟六。若將步驟六改為：“對個體進行直接變異操作”，則為直接變異遺傳算法的計算步驟；若將步驟六改為：“對個體進行探測變異操作”，則為探測變異遺傳算法的計算步驟。

3 澳大利亞計算機應用協會考核題分析

本文將雙重變異遺傳算法（DMGA）與簡化Bishop法相結合編寫Python 程序，以1987年澳大利亞計算機應用協會（ACADS）設計的考核題1（a）、考核題1（c）[18]為求解對象，對這兩道題分別進行50 次獨立計算，并以ACADS 提供的裁判程序答案檢驗算法的尋優效果。然后，在分別僅進行直接變異、探測變異的情況下，對上述兩題也進行50 次獨立計算，對比分析DMGA 算法的性能。

3.1 考核題及其裁判程序答案

考核題1（a）為一均質土質邊坡，如圖6所示，坡高10 m，坡度為26.565°，黏聚力為3 kPa，內摩擦角為19.6°，土體重度為20 kN/m3。要求確定邊坡的臨界滑動面及最小安全系數。

圖6 考核題1（a）坡面示意圖（坐標平移后）[18]Fig.6 Profile of slope in EX1（a）（after coordinates are modified）[18]

考核題1（c）為一非均質土質邊坡，如圖7所示，坡高10 m，坡度為26.565°，黏聚力、內摩擦角及土體重度見表1。要求確定邊坡的臨界滑動面及最小安全系數。

圖7 考核題1（c）剖面圖（坐標平移后）[18]Fig.7 Profile of slope in EX1（c）（after coordinates are modified）[18]

表1 考核題1（c）的材料性質[18]Table 1 Material characteristic in EX1（c）[18]

ACADS 給出的考核題1（a）邊坡最小安全系數的裁判程序答案有：0.990，0.991，1.000；考核題1（c）邊坡最小安全系數的裁判程序答案有：1.385，1.390，1.406。

3.2 計算參數

本文給定圓弧滑動面半徑R、圓心橫坐標xo、圓心縱坐標yo的搜索區間分別為：R∈（0,99]，xo∈

[?99,99]，yo∈[0,99]。

計算所涉及的相關參數取值如下：

（a）滑體土條寬度按R/50 進行取值，且單個條塊不跨越地形剖面的轉折點。當含有地形轉折點、滑動面剪出口及滑動面后緣的條塊，其寬度不足R/50 時，上述相應控制點即為條塊端點；

（b）變異概率計算公式（10）中參數a=0.5，k=0.6；

（c）實數編碼的串長為18（每個變量的編碼長度為6，且變量精確到3 位小數），如圖1所示；

（d）種群數N 為50；

（e）初次最大進化代數T為200，當最大進化代數與最佳個體出現代數之差小于等于30 時，最大進化代數增加30 代，最大進化代數的最大值為500。

3.3 計算結果及過程分析

3.3.1 雙重變異遺傳算法計算結果

表2、表3 分別是基于雙重變異遺傳算法對考核題1（a）、1（c）的計算結果（限于篇幅，僅列出了部分計算結果）。

表2 考核題1（a）的計算結果（DMGA）Table 2 Calculated results of EX1（a）（DMGA）

表3 考核題1（c）的計算結果（DMGA）Table 3 Calculated results of EX1（c）（DMGA）

計算結果統計顯示，對于考核題1（a），50 次計算中，安全系數的最小值為0.985 2，最大值為0.993 7，平均值為0.985 7，與該考核題的裁判程序答案基本一致，說明對于均質邊坡該算法可靠有效；對于考核題1（c），50 次計算中，安全系數的最小值為1.394 9，最大值為1.407 0，平均值為1.397 9，與該考核題推薦的裁判答案基本一致，說明對于非均質邊坡該算法也可靠有效。

3.3.2 三種算法計算結果比較

基于雙重變異遺傳算法（DMGA）、直接變異遺傳算法（DIMGA）、探測變異遺傳算法（DEMGA），分別對兩道考核題計算50 次的統計結果（表4）顯示，DMGA算法搜索得到的安全系數平均值小于DIMGA 算法、DEMGA 算法的計算結果，說明DMGA 算法具有更強的全局搜索能力。并且DMGA 算法計算所得的標準差較小，說明該算法具有較好的魯棒性。

表4 考核題1（a）、1（c）的計算結果統計分析Table 4 Statistical analysis of computation results of EX1（a）and EX1（c）

3.3.3 三種算法計算過程對比

（1）單次計算的進化過程對比

圖8 為采用DMGA 算法、DIMGA 算法、DEMGA算法對考核題1（a）、1（c）進行的第1 次計算的進化過程曲線；表5 為采用DMGA 算法、DIMGA 算法、DEMGA算法對考核題1（a）、1（c）進行的第1 次計算的收斂過程對比表。圖8 顯示DEMGA 算法在計算考核題1（a）時，最大適應度值曲線明顯偏低，表明陷入了局部極值點；DIMGA 算法計算的最大適應度值曲線在進化早期明顯偏低，且對考核題1（a）、1（c）分別在第116代和第151 代才收斂（表5），表明其收斂速度較慢；DMGA 算法分別在第66 代和第72 代時收斂，其計算得到的安全系數低于DIMGA 算法和DEMGA 算法，顯示出優于DIMGA 算法和DEMGA 算法的尋優能力。

表5 考核題1（a）、1（c）算法收斂過程對比表Table 5 Comparison of convergence processes of EX1（a）and EX1（c）

圖8 考核題1（a）、1（c）第1 次計算的進化過程曲線Fig.8 Evolution curve of the first calculation of EX1（a）and EX1（c）

（2）50 次計算的平均進化過程對比

為了進一步對比三種算法進化過程的統計變化規律，根據各算法50 次計算的歷代最佳個體的適應度值計算出歷代最佳個體的平均適應度值，繪制出前200 代的平均進化過程曲線（圖9）。曲線顯示，與上述單次計算對比類似，DEMGA 算法的局部尋優能力較強，在進化的早期，能快速搜索到局部極值點附近，但其跳出局部極值點的能力弱，容易出現早熟現象；DIMGA 算法具有全局搜索的特點，局部尋優能力較弱，收斂速度較慢；DMGA 算法將探測變異與直接變異結合，既具有較強的局部搜索能力，又具有較強的跳出局部極值點的能力，能夠在搜索的廣度與深度上達到較好平衡，因此其全局尋優能力強。

圖9 考核題1（a）、1（c）的平均進化過程曲線Fig.9 Average evolution curve of EX1（a）and EX1（c）

4 工程實例

某海堤邊坡[19?20] 為一非均質土質邊坡，其剖面及土層物理力學參數如圖10所示，現采用雙重變異遺傳算法結合簡化Bishop 法搜索該邊坡的臨界滑動面及最小安全系數。

圖10 邊坡剖面圖[19?20]Fig.10 Slope profile[19?20]

圓弧滑動面半徑及圓心坐標的搜索區間分別為：R∈（0,99]，xo∈[?99,99]，yo∈[0,99]。且需滿足約束條件[20]：6≤yo≤30。式（10）中的參數k=0.7，其他參數設置與3.2 節相同。

經過1 次獨立計算，根據歷代的適應度值（最大值、平均值）與進化代數做出進化過程曲線（圖11），該曲線顯示了歷代進化的個體適應度值的變化情況。在進化的初期（第7 代之前），種群的最大適應度值接近于0，說明進化初期種群中的個體均位于非可行域內。第7 代之后，遺傳進化使種群中的一部分個體進入可行域，種群的最大適應度值及平均適應度值總體上隨進化代數的增加而增加。當進化到第83 代時，種群的最大適應度值為0.370 89，隨后直到第200代進化結束時均未更新，算法收斂。

圖11 進化過程曲線Fig.11 Evolutionary process curve

計算結果顯示，邊坡的安全系數F=1.696 2，臨界滑動面的圓心坐標xo=6.859 m，yo=11.910 m，半徑R=14.910 m，臨界滑動面左側與坡面交點的橫坐標為xL=?2.11 m，臨界滑動面右側與坡面交點的橫坐標為xR=20.55 m。與文獻[20]采用簡化Bishop 法計算的滑面位置（xL=?2.29 m,xR=20.17 m,yo=11.65 m）基本一致，且安全系數小于文獻[20]計算的最小值（Fmin=1.727），說明算法可靠有效，且尋優效果更佳。

5 結論

（1）本文提出的雙重變異遺傳算法，能夠在搜索的廣度與深度上達到較好平衡，使待變異個體適應度值距離當前最佳個體較遠時，能夠向著適應度值增加的方向進行搜索；待變異個體適應度值接近當前最佳個體時，算法既具有局部尋優能力，也具有跳出局部極值點，搜索全局最優的能力。

（2）與僅進行直接變異或探測變異的遺傳算法相比，雙重變異遺傳算法具有更強的全局搜索能力及魯棒性，與簡化Bishop 法結合，能夠有效地搜索出邊坡的圓弧形臨界滑動面。