全因子試驗設(shè)計的土壩變形參數(shù)全局敏感性分析

孫一清,杜一峰,楊慶慶,王桂智,甘 磊,周東昊,陳官運,馮先偉

(1.河海大學(xué) 水利水電學(xué)院，江蘇 南京 210098;2.上海市松江區(qū)水務(wù)管理所，上海 201699;3.揚州市勘測設(shè)計研究院有限公司，江蘇 揚州 225007;4.江蘇省水利勘測設(shè)計研究院有限公司，江蘇 揚州 225127)

1 研究背景

大壩是重要的水利工程建筑物，由于存在各種外部和內(nèi)部載荷，必然會產(chǎn)生相應(yīng)的變形，大壩變形達(dá)到臨界極限時，可能會發(fā)生嚴(yán)重破壞甚至完全坍塌[1]，所以有必要對大壩變形進(jìn)行實時監(jiān)測及數(shù)值模擬研究。在眾多壩型中，土壩是世界上最古老的壩型之一[2]，其中均質(zhì)土壩最為經(jīng)濟(jì)實用，因此關(guān)于均質(zhì)土壩的變形研究較多[1,3-6]。其中，Guo等[6]將克里金替代模型、蒙特卡洛模擬、sobol敏感性分析法與一階可靠性方法相結(jié)合，調(diào)查了輸入變量在土壩穩(wěn)定性研究中的重要性；Dong等[4]以中國某高土石壩為例，采用鄧肯-張E-B模型進(jìn)行靜力分析，并通過與監(jiān)測數(shù)據(jù)對比發(fā)現(xiàn)該模型能夠較好地描述土石壩的變形?？梢?，對于土壩的研究已經(jīng)不再局限于有限元分析，越來越多的學(xué)者開始就方法和本構(gòu)模型對土壩的穩(wěn)定性進(jìn)行研究。

在土石壩變形計算較為常用的本構(gòu)模型中[7]，彈塑性模型的參數(shù)很難獲得，計算過程也非常復(fù)雜。在非線性彈性模型中，鄧肯-張模型(Duncan-Chang E-v模型或Duncan-Chang E-B模型)常被用于分析大壩的結(jié)構(gòu)性態(tài)[8]，其中，鄧肯-張E-B模型中每個參數(shù)均具有清晰的物理和幾何意義，不僅可以通過常規(guī)三軸剪切試驗簡單地獲得[9]，還可以通過智能算法進(jìn)行反分析。但該模型參數(shù)的反演需要建立參數(shù)與變形結(jié)果之間的目標(biāo)函數(shù)，最后轉(zhuǎn)化為多目標(biāo)優(yōu)化問題[10]，并且需要大量的模擬數(shù)據(jù)作為支撐，因而研究者們逐漸開始采用敏感性分析方法對模型參數(shù)進(jìn)行研究。在考慮輸入?yún)?shù)的不確定性對模型響應(yīng)輸出的影響程度時，通過敏感性分析可以給出系統(tǒng)中各參數(shù)的重要程度[6]，即可忽略對結(jié)果影響不顯著的參數(shù)，僅考慮對結(jié)果影響顯著的重要參數(shù)，從而提高研究效率[11]。早期更多地采用單因素進(jìn)行敏感性分析，該方法屬于局部敏感性分析的一種[12]，現(xiàn)有的局部敏感性分析方法主要包括有限差分敏感性分析法、場景分解和廣義龍卷風(fēng)圖、蜘蛛圖和單向靈敏度函數(shù)以及基于微分的方法等[13]。但上述方法在不同程度上具有線性、正態(tài)性和局部變化的局限性[14]。事實上，廣義的敏感性分析方法還包括全局敏感性分析方法，與局部敏感性分析相比，全局敏感性分析不僅可以反映出每個因素的影響，還可以反映所有因素相互作用的影響[15]。全局敏感性分析方法可以分為基于回歸的方法、基于試驗設(shè)計和篩選設(shè)計的方法、基于方差的方法和基于元模型的方法等[16]。

近年來，敏感性分析越來越多地被應(yīng)用于大壩工程研究中。Li等[17]利用基于方差的全局靈敏度分析方法，從3個不同的層面對大壩的靜力性能進(jìn)行了分析；Liang等[18]采用具有近似矩估計的拉丁超立方體抽樣(Latin hypercube sampling，LHS)方法研究高拱壩地震穩(wěn)定性能的參數(shù)敏感性；Chen等[19]用修正的莫里斯法初步分析了耦合E-v和修正Burgers模型中參數(shù)的靈敏度；Yao等[10]用響應(yīng)面法代替人工神經(jīng)網(wǎng)絡(luò)，描述了E-B模型參數(shù)與大壩沉降之間的敏感性關(guān)系；Ren等[20]以陜西省某土石壩為例，采用Morris方法研究了熱液耦合模型參數(shù)對土石壩溫度場的敏感性；Lakehal等[21]以均質(zhì)土邊坡為例，采用中心組合設(shè)計試驗方法，研究了E-B模型參數(shù)與安全系數(shù)之間的敏感性。與上述敏感性分析方法一樣，全因子試驗設(shè)計同樣屬于試驗設(shè)計和篩選設(shè)計方法中的一種[22]，其雖然需要較多的試驗次數(shù)，但優(yōu)點是可以得到準(zhǔn)確可靠的結(jié)果，并計算出所有因素的主效應(yīng)和交互效應(yīng)，主效應(yīng)即一個因素各水平變化造成結(jié)果變化的程度，且該變化程度不受其他因素各水平變化的影響，交互效應(yīng)即一個因素各水平變化造成結(jié)果變化的程度，且該變化程度受其他因素各水平變化的影響。Golshani等[23]以尾礦壩中的煤精礦為研究對象，采用全因子設(shè)計分析酸堿度、粒度等因素對于煤精礦硫還原的敏感性影響；Kadeethum等[24]采用全因子試驗設(shè)計，針對多孔介質(zhì)裂縫中流體力學(xué)對于油井產(chǎn)能的影響進(jìn)行了統(tǒng)計敏感性分析。

目前，全因子試驗設(shè)計鮮有被用于水工程領(lǐng)域的案例，且多數(shù)水工程中的敏感性分析方法未考慮因素間的交互效應(yīng)，因而缺少對因素的全局敏感性分析。本文以甘肅省某黃土均質(zhì)壩為研究對象，通過有限元法分析程序構(gòu)建反映土石壩結(jié)構(gòu)的靜力分析模型，求解大壩位移變化，再借助Plackett-Burman試驗設(shè)計對均質(zhì)土壩變形參數(shù)進(jìn)行顯著性篩選，并利用全因子試驗設(shè)計對影響均質(zhì)土壩變形的顯著參數(shù)進(jìn)行全局敏感性分析，最后建立了相應(yīng)的回歸模型，以期為黃土均質(zhì)壩變形模擬計算參數(shù)的選取提供參考。

2 敏感性分析方法

以黃土壩殼全計算區(qū)域為研究對象，采用全局敏感性分析方法研究壩料密度和鄧肯-張E-B模型參數(shù)對黃土均質(zhì)壩的敏感性。在全局敏感性分析中，由于全因子試驗設(shè)計只能考慮至多5個因素的敏感性以及因素間的交互效應(yīng)，所以有必要先對考察因素進(jìn)行科學(xué)合理地顯著性篩選，而Plackett-Burman試驗設(shè)計可以實現(xiàn)這一過程，經(jīng)篩選后再進(jìn)行全因子試驗，以得到全局敏感性分析結(jié)果。全局敏感性分析方法具體過程如下：

(1)Plackett-Burman試驗設(shè)計。利用Plackett-Burman試驗設(shè)計對8個計算參數(shù)進(jìn)行試驗設(shè)計，運用方差分析法分析結(jié)果，篩選掉其中不顯著的參數(shù)，保留顯著的參數(shù)。

(2)全因子試驗設(shè)計。對通過Plackett-Burman試驗設(shè)計保留下來的參數(shù)進(jìn)行全因子試驗設(shè)計，對這些參數(shù)進(jìn)行主效應(yīng)分析和交互效應(yīng)分析，繪制出等值線圖及響應(yīng)曲面圖，最后使用方差分析，尋找敏感性較強(qiáng)的幾個參數(shù)及其組合。

方差分析借助統(tǒng)計量F值來判定顯著性，其計算方法如下：

(1)

式中：MSi為各因素的均方；MSe為誤差的均方。

若Fi>F0.01，說明該因素影響明顯顯著，記為“**”；若F0.01>Fi>F0.05，說明該因素影響顯著，記為“*”；若F0.05>Fi>F0.1，說明該因素有影響，記為“☉”；若F0.1>Fi>F0.2，說明該因素有一定影響，記為“△”；若Fi

鄧肯-張E-B模型共有10個參數(shù)，分別為黏聚力C、初始內(nèi)摩擦角φ0、摩擦角Δφ、破壞比Rf、初始彈性模量基數(shù)K、初始體積模量基數(shù)Kb、初始彈性模量指數(shù)n、初始體積模量指數(shù)m、卸荷再加荷時的彈性模量基數(shù)Kur、卸荷再加荷時的彈性模量指數(shù)nur。在參數(shù)選擇時，由于黃土土壩施工填筑過程中，壩殼室中處在荷載狀態(tài)，且黃土顆粒為散粒體材料，故參數(shù)C、Kur和nur不參與討論，同時施工過程中壩料密度ρ也是一個影響黃土壩殼變形的因素，因此選擇ρ、φ0、Δφ、Rf、K、Kb、n、m這8個參數(shù)進(jìn)行敏感性分析。

2.1 Plackett-Burman試驗設(shè)計

Plackett-Burman試驗設(shè)計被廣泛應(yīng)用于實際工作中，其主要作用是篩選試驗因素。分別對多個影響試驗指標(biāo)強(qiáng)弱程度不詳?shù)脑囼炓蛩厝「?+1)低(-1)兩個水平，借助試驗設(shè)計表格研究試驗因素多個水平間變化對于試驗指標(biāo)的影響，篩選出對試驗指標(biāo)影響顯著的試驗因素，以便減少后續(xù)研究的工作量。Plackett-Burman試驗設(shè)計按規(guī)則生成，排列可具有不唯一性，實際研究中需保留3個以上虛擬變量，即對于M次試驗，所研究因子應(yīng)小于等于M-4個。且試驗次數(shù)M應(yīng)為4的倍數(shù)，但不包括2的冪次方。常用的取值為M=12、20、24、28、36等。

2.2 全因子試驗設(shè)計

全因子試驗設(shè)計不僅可以分析考察因素對于試驗指標(biāo)的敏感性，得到敏感性排序，還可以有效地估計出所有因素的主效應(yīng)和各階交互效應(yīng)，生成交互效應(yīng)較強(qiáng)組合的等值線圖和響應(yīng)曲面圖，并對結(jié)果進(jìn)行方差分析。全因子試驗設(shè)計要求對所考察的全部因素的全部水平至少進(jìn)行一次試驗，其分析結(jié)果真實可靠。在需要考慮因素之間交互效應(yīng)以獲得較精確的分析結(jié)論時，常選擇全因子設(shè)計。全因子試驗設(shè)計通常不超過5個因素且水平為2。

3 工程實例

某水庫大壩工程為黃土均質(zhì)壩，最大壩高為125.00 m，壩頂高程為1 653.00 m。正常蓄水位為1 650.00 m，總庫容為1.40×108m3。形成的有限單元網(wǎng)格共有11 324個結(jié)點和10 877個單元。黃土均質(zhì)壩模型范圍及有限元網(wǎng)格劃分見圖1，壩體各料區(qū)的計算參數(shù)取值如表1所示，分級加載及蓄水過程如表2所示。大壩運行期壩體典型斷面的位移分布如圖2。整體來看，該壩的變形分布符合一般均質(zhì)土壩的分布規(guī)律，竣工期壩體水平位移以壩軸線為中線對稱分布，運行期壩體水平位移在水壓力的作用下稍有偏移；豎向位移呈現(xiàn)均勻下降趨勢，位移最大值在壩體高度的2/3處，且運行期最大值大于竣工期最大值。

圖2 實例工程運行期壩體典型斷面位移分布(單位：mm)

表1 實例工程壩體有限元靜力計算參數(shù)

表2 實例工程有限元計算分級加載及蓄水過程

圖1 實例工程黃土壩模型范圍及有限元網(wǎng)格劃分

4 敏感性分析

4.1 Plackett-Burman試驗

對ρ、φ0、Δφ、Rf、K、Kb、n、m這8個參數(shù)進(jìn)行Plackett-Burman試驗設(shè)計，每個參數(shù)考察兩個水平，高水平(+1)取試驗設(shè)計參數(shù)的1.2倍，低水平(-1)取試驗設(shè)計參數(shù)的0.8倍，如表3所示。具體試驗方案和試驗結(jié)果(豎向位移V、向下游水平位移H1和向上游水平位移H2)如表4所示，分析結(jié)果如表5所示。

表3 Plackett-Burman試驗各因素水平

表4 Plackett-Burman試驗方案及各方案試驗結(jié)果

由表5中試驗分析結(jié)果可以看出，通過Plackett-Burman試驗設(shè)計對8個變形計算參數(shù)進(jìn)行篩選，各參數(shù)關(guān)于豎向位移V的敏感性排序為Rf>φ0>ρ>Δφ>n>Kb>m>K，其中對于豎向位移影響明顯顯著的參數(shù)為Rf、φ0和ρ，影響顯著的因素為Δφ；各參數(shù)關(guān)于向下游水平位移H1的敏感性排序為Rf>φ0>ρ>Δφ>n>m>K>Kb，其中對于向下游水平位移影響明顯顯著的參數(shù)為Rf和φ0，影響顯著的參數(shù)為ρ和Δφ；各參數(shù)關(guān)于向上游水平位移H2的敏感性排序為Rf>φ0>Δφ>ρ>n>Kb>m>K，其中對于向上游水平位移影響明顯顯著的參數(shù)為Rf和φ0，影響顯著的參數(shù)為ρ和Δφ。

表5 各參數(shù)關(guān)于壩體不同方向位移的試驗敏感性分析結(jié)果

綜合來看，Rf、φ0、ρ和Δφ這4個參數(shù)對于黃土均質(zhì)壩變形的影響更加顯著，因而在計算參數(shù)中需要重點考慮，而n、Kb、m和K這4個參數(shù)可以視為對大壩變形無影響。

4.2 全因子試驗

對上述Plackett-Burman試驗設(shè)計篩選出的對黃土均質(zhì)壩變形影響顯著的參數(shù)Rf、φ0、ρ和Δφ進(jìn)行全因子試驗設(shè)計，研究該4個參數(shù)兩兩之間的交互效應(yīng)。每個參數(shù)考察兩個水平，高水平(+1)取試驗設(shè)計參數(shù)的1.2倍，低水平(-1)取試驗設(shè)計參數(shù)的0.8倍。具體試驗方案及結(jié)果如表6所示。

表6 全因子試驗方案

4個參數(shù)對大壩各方向變形的主效應(yīng)如圖3所示。由圖3(a)和3 (c)可知，造成豎向位移和向上游水平位移變化差異明顯的參數(shù)為ρ、φ0和Rf，且隨著ρ和Rf的增大，豎向位移和向上游水平位移的數(shù)值顯著減小，隨著φ0的增大，豎向位移和向上游水平位移的數(shù)值顯著增大；而無論其他3個參數(shù)如何改變，造成豎向位移和向上游水平位移的數(shù)值變化差異不明顯的參數(shù)為Δφ。由圖3(b)可知，造成向下游水平位移的數(shù)值變化差異明顯的參數(shù)為ρ、φ0和Rf，且隨著ρ和Rf的增大，向下游水平位移的數(shù)值顯著增大，隨著φ0的增大，向下游水平位移的數(shù)值顯著減小，而無論其他3個參數(shù)如何改變，造成向下游水平位移的數(shù)值變化差異不明顯的參數(shù)為Δφ[25]。

圖3 4個參數(shù)對大壩各方向變形的主效應(yīng)圖

4個參數(shù)兩兩交互對大壩變形的效應(yīng)如圖4所示。由圖4可以看出，對于豎向位移、向下游水平位移和向上游水平位移來說，僅有ρ×Rf和φ0×Rf兩組交互效應(yīng)的效應(yīng)線相對角度明顯較大，說明ρ×Rf和φ0×Rf的交互效應(yīng)顯著[26]，即一個因素各水平變化造成結(jié)果的差異會受另一因素各水平變化的影響，而ρ×φ0、ρ×Δφ、φ0×Δφ、Rf×Δφ4組交互效應(yīng)的效應(yīng)線相對角度較小，說明交互效應(yīng)不顯著。

圖4 4個參數(shù)兩兩交互對大壩變形的效應(yīng)圖

由交互效應(yīng)分析可知，ρ×Rf和φ0×Rf兩組效應(yīng)需要重點關(guān)注。設(shè)定另外兩個變量在最佳參數(shù)時，繪制兩組效應(yīng)下大壩豎向位移、向下游水平位移和向上游水平位移的等值線，如圖5所示。由圖5可看出，在保持φ0為26°和Δφ為5°不變時，隨著ρ和Rf的增大，豎向位移和向上游水平位移的數(shù)值逐漸減小(位移量增大)，而向下游水平位移的數(shù)值逐漸增大(位移量增大)；保持ρ為17 g/cm3和Δφ為5°不變，隨著φ0的減小和Rf的增大，豎向位移和向上游水平位移的數(shù)值逐漸減小(位移量增大)，而向下游水平位移的數(shù)值逐漸增大(位移量增大)。再對ρ×Rf和φ0×Rf兩組交互效應(yīng)下繪制豎向位移、向下游水平位移和向上游水平位移相應(yīng)的響應(yīng)曲面圖，如圖6所示。由圖6可見，對于3個方向的位移而言，ρ和Rf越大，則位移越大，當(dāng)ρ為20.4 g/cm3、Rf為0.9時，位移達(dá)到最大值；φ0越小、Rf越大，則位移越大，當(dāng)φ0為20.8°、Rf為0.9時，位移到達(dá)最大值。

圖5 ρ×Rf和φ0×Rf兩組交互效應(yīng)下大壩各方向位移等值線 圖6 ρ×Rf和φ0×Rf兩組交互效應(yīng)下大壩各方向位移響應(yīng)曲面

方差分析結(jié)果見表7，表7中方差分析結(jié)果為不同參數(shù)及其組合的敏感性統(tǒng)計量Fi值，F(xiàn)i值越大則敏感性越強(qiáng)。由于Plackett-Burman試驗篩選出的顯著參數(shù)Rf、φ0、ρ和Δφ相互之間的三階及四階交互效應(yīng)顯著性太小，在此不予分析。

表7 4個參數(shù)及其交互效應(yīng)敏感性的方差分析結(jié)果

由表7可知，對于豎向位移，Rf、φ0、ρ和Δφ4個參數(shù)及其交互效應(yīng)的敏感性排序為Rf>φ0>φ0×Rf>ρ>ρ×Rf>ρ×φ0>Δφ>Rf×Δφ>φ0×Δφ>ρ×Δφ，其中，Rf、φ0、φ0×Rf和ρ的影響明顯顯著，ρ×Rf和ρ×φ0的影響顯著；對于向下游水平位移，Rf、φ0、ρ和Δφ4個參數(shù)及其交互效應(yīng)的敏感性排序為Rf>φ0>φ0×Rf>ρ>ρ×Rf>ρ×φ0>Δφ>Rf×Δφ>φ0×Δφ>ρ×φ0。其中，Rf、φ0、φ0×Rf和ρ×Rf的影響明顯顯著，ρ×φ0有影響；對于向上游水平位移，Rf、φ0、ρ和Δφ4個參數(shù)及其交互效應(yīng)的敏感性排序為Rf>φ0>φ0×Rf>ρ>ρ×Rf>ρ×φ0>Δφ>Rf×Δφ>φ0×Δφ>ρ×φ0。其中，Rf、φ0、φ0×Rf和ρ×Rf的影響明顯顯著，ρ×φ0有影響。

通過方差分析得到4個參數(shù)及其交互效應(yīng)的敏感性排序為Rf>φ0>φ0×Rf>ρ>ρ×Rf>ρ×φ0>Δφ，與前文Plackett-Burman試驗分析得到的敏感性排序Rf>φ0>ρ>Δφ相符合。而φ0×Rf和ρ×Rf的敏感性是需要重點關(guān)注的。

將方差分析結(jié)果繪制成大壩變形標(biāo)準(zhǔn)化效應(yīng)的Pareto圖，如圖7所示。由Pareto圖可以更加直觀地觀察判斷出各因素的顯著性。其判斷方法是：如果某因素對應(yīng)的矩形條超過標(biāo)準(zhǔn)化效應(yīng)閾值，則該因素為顯著影響因素[25]。

圖7 大壩變形標(biāo)準(zhǔn)化效應(yīng)的Pareto圖

由圖7可知，Rf、φ0、ρ3個參數(shù)的主效應(yīng)顯著，φ0×Rf、ρ×Rf的二階交互效應(yīng)顯著，Δφ及其余二階效應(yīng)均不顯著，與前文中主效應(yīng)及交互效應(yīng)分析結(jié)果一直。將Δφ及其余二階項剔除，對大壩變形標(biāo)準(zhǔn)化效應(yīng)的Pareto圖中顯著的因素建立回歸模型。

豎向位移的回歸模型為：

V=7.89+0.269ρ-0.441φ0-15.05Rf-

0.570ρ·Rf+0.766φ0·Rf

(2)

向下游水平位移的回歸模型為：

H1=-7.61+0.455ρ+0.437φ0+14.62Rf-

0.817ρ·Rf+0.784φ0·Rf

(3)

向上游水平位移的回歸模型為：

H2=8.94+0.633ρ-0.579φ0-16.91Rf-

1.107ρ·Rf+1.026φ0·Rf

(4)

對于豎向位移回歸模型，R2=0.8905，R2(預(yù)測)=0.7198，R2(調(diào)整)=0.8358；對于向下游水平位移回歸模型，R2=0.9493，R2(預(yù)測)=0.8702，R2(調(diào)整)= 0.9239；對于向上游水平位移回歸模型，R2=0.9171，R2(預(yù)測)= 0.7877，R2(調(diào)整)= 0.8756，表明各方向位移的回歸模型的擬合性較優(yōu)。

綜合來看，Rf、φ0、φ0×Rf、ρ、ρ×Rf對黃土均質(zhì)壩變形敏感性更強(qiáng)，可以在考慮交互效應(yīng)的正交試驗中重點考察。

5 結(jié) 論

為了全面探究大壩變形參數(shù)對于均質(zhì)土壩變形的影響，本文選取甘肅省某黃土均質(zhì)土壩作為案例，通過有限元分析程序建立了均質(zhì)土壩靜力分析模型，分析了該壩在竣工期和運行期的位移變化規(guī)律；使用Plackett-Burman試驗設(shè)計對鄧肯-張E-B模型參數(shù)及壩料密度進(jìn)行對于大壩變形的顯著性篩選；對篩選出的顯著參數(shù)進(jìn)行全因子試驗設(shè)計，研究對于均質(zhì)土壩變形影響顯著的參數(shù)及其組合，得到以下主要結(jié)論：

(1)該土壩的變形分布符合一般均質(zhì)土壩的分布特點，竣工期壩體水平位移呈對稱分布，運行期壩體水平位移在水壓力的作用下由對稱分布逐漸向下游變形。豎向位移呈現(xiàn)均勻下降趨勢，最大值在壩體高度的2/3處。

(2)對于Plackett-Burman試驗設(shè)計中的ρ、φ0、Δφ、Rf、K、Kb、n、m這8 個參數(shù)，需要重點考慮的是Rf、φ0、ρ和Δφ，而n、Kb、m和K4個參數(shù)可以視為對大壩變形無影響。

(3)在全因子試驗設(shè)計中分別進(jìn)行了主效應(yīng)分析、交互效應(yīng)分析、方差分析，并且對交互效應(yīng)明顯顯著的參數(shù)組合進(jìn)行等值線圖和響應(yīng)曲面圖分析，并建立相應(yīng)的回歸模型。得到的結(jié)論是：被篩選出的4個參數(shù)中，Rf、φ0和ρ3個參數(shù)的主效應(yīng)明顯顯著，且ρ和Rf、φ0和Rf兩個參數(shù)的交互效應(yīng)敏感性不能被忽視。