基于免疫遺傳算法的堆石壩流變反演分析

(上海勘測設計研究院有限公司，上海 200434)

流變是高土石壩變形中的一個重要因素，在設計時正確合理地計入流變變形的影響，并通過施工措施降低堆石體的變形量，是保證大壩安全運行的關鍵。因此在高土石壩設計與施工中，運用數值分析方法計算預測壩體的變形及其趨勢，進行堆石壩的設計優化，是高土石壩持續發展的重要因素[1]。

利用原位監測資料，通過計算和反饋分析土石壩實際工作狀況[2]，對已建成大壩可更好地評價其安全狀況；對在建的大壩，通過原型監測資料進行反分析，對優化堆石壩設計以及改進施工也有重大意義[3]。在土石壩應力變形計算中，參數的選取準確與否關系到數值計算結果是否能符合土石壩的實際情況。反演分析應用現場監測資料反算土石壩參數，能較精確地確定土石壩參數，把數值分析和監測技術有機地結合在一起[4]。近年來，引入了人工神經網絡等智能算法，進一步優化了參數的反分析效率。2005年，張丙印、袁會娜運用人工神經網絡和演化算法進行位移反分析[5]。2006年，朱晟、張美英將基于免疫學原理的遺傳算法應用于堆石壩材料力學參數的位移反分析中，該方法克服了標準遺傳算法收斂方向無序性且容易局部早熟的缺點，在反分析中取得了很好的應用效果[6-7]。

在瀑布溝心墻堆石壩材料室內三軸試驗成果基礎上，本文通過瀑布溝心墻堆石壩的原型監測資料，基于免疫遺傳算法對其流變參數進行反演，并對反演所得參數的合理性進行了驗證分析。

1 計算模型及免疫遺傳反演算法

1.1 增量流變本構模型

國內對堆石流變的研究，最早始于由沈珠江根據堆石流變試驗提出的三參數流變模型[8]，該模型在土石壩的流變分析中得到了廣泛的應用。三參數流變模型以應力全量形式推導，未考慮應力加載路徑。當土石壩有限元分析計算采用增量荷載時，無法考慮該級增量荷載之前流變的影響，堆石體流變的遺傳特性難以體現。因而，本文采用由朱晟在此基礎上提出的可反映應力路徑影響的7 參數增量流變模型[9]，該模型已應用于水布埡土石壩計算分析中[10]。

結合堆石料三軸流變試驗資料，增量流變模型從中得到流變體積模量和剪切模量等模型參數，模型參數分別為ks，ns，kv，Rsf，c，α，其數值可根據三軸流變試驗成果或大壩原型監測資料反演確定。

(1)

(2)

式中，σvi，σsi分別為體積應力和廣義剪應力；Ki，Gi為控制滯后變形部分的參數，稱為流變體積模量和流變剪切模量，是相應加載應力路徑下應力狀態的函數；c，α為試驗常數。

不同應力狀態下，體積應力和最終體積流變近似為線性關系，可以認為流變體積模量一般隨試驗圍壓的變化而變化，以式(3)近似表示：

(3)

式中，Pa為大氣壓力，kv，nv為模型參數，可根據試驗結果整理得到。

不同應力狀態下廣義剪應力和最終剪切流變近似為雙曲線關系，流變剪切模量為

(4)

式中，S為應力水平，Rsf為破壞比，ks，ns為模型參數，可根據試驗資料整理得到。

1.2 目標函數及反演參數

堆石料占大壩體積的絕大部分，是影響壩體沉降的最主要因素，同時參考監測儀器的布置位置，在反分析時主要反演大壩心墻、下游主堆石區和次堆石區的壩料參數[11]。

通過監測點的實測位移量與計算位移值的差異，構造目標函數，形成IGA 抗原。在監測儀器埋設初期時段，各種因素的干擾較大，目標函數需要采用增量時段計算；同時，為了考慮不同測點測值權重和消除不同物理量的量綱，建立目標函數[12-13]如下：

(5)

增量流變模型待反演流變參數為X={ks，ns，kv，nv，Rsf，c，α}。

1.3 基于IGA的流變參數反演算法

1.3.1抗體編碼

采用浮點數編碼方式，模型共7個參數，因此就有7個基因。假設初設種群數為N，隨機產生7×N個小于1的實數，αi1,αi2,…,αi7(i=1,2,…,N)取初始種群第i個抗體的7個基因分別為

(6)

式中，aji,bji為第i個抗體中第j個基因的取值范圍，可參考類似工程分析取值。重復N次，就可以取得初始抗體群X1,X2,…,XN。

1.3.2構造適應度函數

首先將獲得的初始抗體群按目標函數計算值進行升序排列，采用Michalewicz提出的非線性排序的選擇概率計算公式建立適應度函數：

f(Xi)=α(1-α)i-1

(7)

式中，i為個體排序序號，α∈(0,1)。

從式(7)可以看出，個體被選中的概率與個體的排序位置有關，排序位置越前適應度越大，因此排序前的個體被選中概率大，合理的適應度函數可避免算法局部收斂或搜索停滯。

1.3.3記憶細胞庫更新

將每代種群中的所有個體根據適應度值的大小進行降序排列，選取前M個個體存入記憶細胞庫，完成記憶細胞庫得換代更新過程。

1.3.4基于濃度的適應度調整

為保證個體的多樣性，提高算法的全局搜索能力，根據抗體間歐氏距離和適應度來計算抗體間的相似度和抗體濃度。記抗體Xi和Xj的歐氏距離為D(Xi,Xj)，適應度分別為f(Xi)和f(Xj)，取適當常數η>0，t>0，如果兩個抗體滿足式(8)：

(8)

則認為抗體Xi和Xj相似，抗體Xi的濃度即為與抗體XJ相似的抗體個數除以抗體總數，記為Ci。則抗體Xi適應度函數根據濃度修正為

(9)

式中，α，β∈(0,1)，fmax(Xi)為所有抗體的最大適應度值，Ci為抗體Xi的濃度。

從式(9)可以得出：對于抗體Xi，采用上述濃度調節后，濃度對選擇幾率影響減小，若抗體濃度高且適應度較高，獲得的適應度補償就較少，可避免使計算陷入局部極優；而抗體濃度不高但適應度高，得到適應度補償就會相對較高，增加了該抗體的選擇概率。

1.3.5選擇操作

個體的選取概率P由根據濃度修正過的適應度函數值計算得來：

(10)

圖1 土石壩流變參數反演分析流程Fig.1 Program of back analysis for earth-rock dam rheological parameters

基于選取概率，用輪盤賭選擇方式來產生新的個體，適應度值較高的個體容易被選擇，同時由于是采用相對值進行判斷，避免過大適應度值的絕對性影響，避免過早收斂。

1.3.6交叉與變異

交叉運算選擇算術交叉算子，由兩個個體的線性組合而產生兩個新的個體，交叉概率取0.8。變異采用的是均勻變異算子，確定作為變異點的各個基因座的位置后，通過變異概率對各個變異點進行隨機變異，變異概率取0.1。

1.3.7群體更新

將最終獲得的個體和記憶細胞庫中本身的M個體合并在一起，得到一個新種群，然后再根據適應度值進行降序排列，將前M個個體更新進入記憶細胞庫中，下一代群體則為剩余N個個體。

1.3.8終止條件

當目標函數值低于某個值或進化代數超過預先設定值，則反演計算終止，此時得到參數即為所需值。

結合免疫遺傳算法和堆石壩流變有限元分析過程，流變參數反演過程如圖1所示。

每個壩體分區堆石料的力學特性參數包括瞬變參數(鄧肯張E-v模型參數：Rf，K，n，G，F，D)和流變參數(ks，n，kv，nv，Rsf，c，α)。在填筑期，占壩體大部分體積的堆石料的瞬時變形占主導，影響位移計算值的主要因素是E-v模型參數，因此在反演獲得流變參數后還需要對室內三軸試驗的E-v模型參數進行反演修正。

2 實際工程應用

2.1 瀑布溝堆石壩概況

瀑布溝水電站位于大渡河中游、四川省漢源縣和甘洛縣境內，最大壩高186 m，總庫容53.90億m3。水庫正常蓄水位850.00 m，死水位790.00 m。在壩體中布置的各種監測儀器中，壩體位移的原型監測數據相對準確，取得的資料也較多，可以作為反演分析的基礎。瀑布溝心墻堆石壩在0+240共布設3組電磁式沉降儀VE1～VE3，其具體分布見圖2。

2.2 流變參數反演分析過程

在對堆石料流變參數進行反演之前，需要首先確定大壩各個分區的瞬變參數。根據室內三軸試驗提供的資料，瀑布溝心墻壩各個分區的參數見表1，并在此基礎上進行流變參數的反演。

大壩現場實際測量的變形觀測值包括瞬間變形分量和流變變形分量。為了盡量排除自重等外荷載所引起的瞬時變形分量的干擾，從堆石體的總變形中獲得合理的流變變形，選擇1～3個月大壩停工階段或者大壩竣工后外荷載基本不變的時期，此時堆石的變形主要為流變變形，反演分析得出的流變參數較為合理。

根據施工情況，瀑布溝心墻堆石壩不存在完全停工階段，但對某些測點來說，在某些時段壩體的填筑方量較小且對監測點影響較小，瞬間變形影響較小，適于進行流變分析。因此，可以認為此時段內監測點相對位移為該時段內對應節點的流變量。根據瀑布溝施工過程，共選取3段時間進行分析：① 2008年7月21日至2008年9月5日，該時段內大壩心墻及下游主次堆石區維持至769 m高程不變，只有上游堆石區局部由769 m填筑至775 m高程，如圖3所示。② 2008年10月21日至2008年12月15日，該時段大壩心墻及下游主次堆石區填筑至795 m高程，只有上游堆石區780 m至786 m高程局部進行堆石體填筑，如圖4所示。③ 2009年4月25日至2009年6月23日，大壩進行壩體下游壓重部分施工，如圖5所示。可基本認為這3段參考水布埡面板堆石壩施工與運行性狀反演研究[10]，流變參數初始取值范圍見表2。

圖2 瀑布溝0+240斷面電磁式沉降儀布置Fig.2 The electromagnetic settling instrument at 0+240 section of Pubugou Core Rockfill Dam

圖3 截止至2008.9.5 大壩斷面(0+240)Fig.3 The dam section up to 2008.9.5(0+240)

時間內大壩各監測點的沉降值主要是流變變形，可以利用這段時間內測得的資料進行流變參數的反演分析[14-16]。

表1 瀑布溝心墻壩壩料E-v模型參數Tab.1 E-v mode parameters of Pubugou core rockfill dam

目標函數的進化過程如圖6所示，到29代時，目標函數基本穩定，到35代時，停止計算，得到反演參數。心墻及主次堆石反演所得流變參數見表3。

2.3 E-v模型參數反演修正

分析時間為2008年1月31日至12月31日，大壩從738.0 m高程填筑到799.0 m高程。參考室內試驗參數，擬定心墻堆石、主堆石和次堆石區鄧肯張E-v模型參數的搜索區間如表4所示，最后得到的反演修正參數如表5所示。

3 反演成果分析驗證

為驗證所獲得的參數的合理性，選取典型監測結點，分析其沉降、水平位移、應力隨時間變化的過程，并與實測值進行比較。0+240斷面典型監測點位置如圖7所示。

表2 流變模型參數的初始取值范圍Table.2 The variation range of rheological parameters

表3 反演所得的流變參數Table.3 The rheological parameters obtained from the inversion

3.1 堆石測點時程曲線對比分析

3.1.1沉 降

圖8為典型測點的沉降時程曲線對比圖。由圖中可見：典型測點的計算沉降與實測沉降數值上相差不大，隨時間的變化規律也基本相同。但在填筑初期，計算沉降與實測沉降較為吻合，填筑后期，計算沉降略大于實測沉降。主要原因為僅以2008年7月21日至9月5日實測沉降反演得到的參數僅能較好地反映此時間段以前的土石壩堆石流變及瞬變參數，對于之后的填筑過程模擬會存在一定的誤差。總體而言，各典型測點的計算沉降與實測沉降吻合度較好，說明由反演所得的土石壩瞬變參數和增量流變參數較為合理，能滿足對土石壩進行數值分析的要求。

圖4 截至2008.12.15 大壩斷面(0+240)Fig.4 The dam section up to 2008.12.15(0+240)

圖5 截至2009.06.23 大壩斷面(0+240)Fig.5 The dam section up to 2009.6.23(0+240)

圖7 0+240斷面典型測點位置示意Fig.7 Position diagram of typical monitoring point at 0+240 section

圖6 目標函數隨進化代數變化圖Fig.6 Change process of target function

表4 E-v模型參數的初始取值范圍Table.4 The variation range of E-v mode parameters

3.1.2水平位移

圖9為典型測點水平位移時程曲線對比圖。各典型測點的計算與實測位移在數值上相差不大，變化規律類似，但存在部分數值在填筑初期吻合度較差。這是因為在水平位移較小情況下，監測儀器存在一定的初始誤差。

3.1.3應 力

從圖10中可以看出，土石壩堆石的計算正應力變化規律與實測正應力變化規律一致，并且數值上也相差不大，說明了反演所得參數較為合理。

3.2 防滲墻測點時程曲線對比分析

與堆石位移監測值相比較，防滲墻位移測值較為精確。由圖11可以看出，由反演所得參數計算得到的防滲墻水平位移與實測水平位移較為接近，且變化規律一致，較為準確地模擬了防滲墻的變形情況，說明了反演所得土石壩參數的合理性。

3.3 壩體變形應力分析

大壩有限元計算分析采用鄧肯張E-v模型和增量流變模型，參數均采用上述分析成果。

3.3.1壩體變形

圖12～13為0+240號斷面竣工期壩體上下游方向水平位移和豎向位移等值線分布圖。竣工期，0+240斷面上下游方向位移均指向壩坡方向，向上游的最大水平位移為42.6 cm，位于上游高程730 m處；而向下游移動的最大水平位移為42.0 cm，位于下游高程720 m處。壩體的最大沉降為195.4 cm，位于高程740 m處的心墻部位。

表5 反演所得的E-v模型參數Table.5 The E-v mode parameters obtained from inversion

3.3.2壩體應力

圖14～15分別為0+240斷面竣工期壩體主應力等值線分布圖。此時壩體主應力基本與壩坡平行，并從壩頂部向壩基處呈現加大的趨勢。壩體第一主應力極值為3.60 MPa，第三主應力極值為1.30 MPa。

竣工期壩體最大沉降為195.4 cm，約占壩高的1.1%，向下游水平位移增大至42.6 cm，向上游最大位移為42.0 cm，基本符合瀑布溝心墻壩的變形應力特性，與實測資料吻合良好。

圖8 堆石典型監測點的沉降時程Fig.8 Settling procedure of typical monitoring points of rockfill

圖10 典型監測點的土壓力時程Fig.10 Soil pressure at typical monitoring points of rockfill

圖9 堆石典型監測點的水平位移時程Fig.9 Horizontal displacement at typical monitoring points of rockfill

圖11 防滲墻典型監測點的水平位移Fig.11 Horizontal displacement at typical monitoring points of cut-off wall

圖12 竣工期0+240斷面水平向位移(單位：cm)Fig.12 Horizontal displacement of 0+240 section at completion period (unit: cm)

圖13 竣工期0+240斷面豎直向位移(單位：cm)Fig.13 Vertical displacement of 0+240 section at completion period (unit: cm)

圖14 竣工期0+240斷面第一主應力(單位：MPa)Fig.14 The major principal stress of 0+240 section at completion period (unit: MPa)

圖15 竣工期0+240斷面第三主應力(單位：MPa)Fig.15 The minor principal stress of 0+240 section at completion period (unit: MPa)

4 結 論

(1) 對典型測點的位移應力及壩體變形應力分析表明，運用免疫遺傳算法對瀑布溝心墻堆石壩反演的流變參數較為真實地反映了土石壩堆石料的流變特性。

(2) 在免疫遺傳反分析過程中，每一次求解過程都需反復調用正分析程序，以優化初始參數搜索空間，減少目標函數穩定進化代數，這對減少計算時間，提高反演效率有著重要意義。