江 巍，閆金洲，歐陽曄，劉立鵬，鄭 宏

(1.三峽庫區(qū)地質(zhì)災(zāi)害教育部重點(diǎn)實(shí)驗(yàn)室，湖北 宜昌 443002；2.三峽大學(xué) 土木與建筑學(xué)院，湖北 宜昌 443002；3.流域水循環(huán)模擬與調(diào)控國家重點(diǎn)實(shí)驗(yàn)室，北京 100038；4.北京工業(yè)大學(xué) 建筑工程學(xué)院，北京 100124)

作為典型的非連續(xù)數(shù)值模擬方法，顆粒離散元法[1]可有效模擬巖土介質(zhì)的開裂、分離等非連續(xù)變形行為，與以有限元為代表的連續(xù)數(shù)值模擬方法相比，更適用于邊坡漸進(jìn)變形破壞的過程模擬，并被成功應(yīng)用于多處地質(zhì)災(zāi)害體的失穩(wěn)機(jī)理揭示，如2017年6月24日四川茂縣新磨滑坡[2]、2018年5月5日四川樂鄉(xiāng)市馬邊滑坡[3]等。周健等[4]將強(qiáng)度折減法引入顆粒離散元中，執(zhí)行邊坡穩(wěn)定性評(píng)價(jià)，形成邊坡穩(wěn)定性分析的強(qiáng)度折減顆粒離散元法，在計(jì)算安全系數(shù)的同時(shí)可預(yù)測(cè)邊坡破壞形態(tài)。Bao等[5]將該方法應(yīng)用于某大型鐵礦廢渣堆場(chǎng)的穩(wěn)定性評(píng)價(jià)，發(fā)現(xiàn)潛在破壞面形態(tài)與現(xiàn)場(chǎng)變形跡象吻合較好，證實(shí)了該方法的有效性。近年來，強(qiáng)度折減顆粒離散元法在臨界破壞準(zhǔn)則方面取得較好的研究進(jìn)展，如：田雷等[6]繪制坡頂顆粒平均位移與折減系數(shù)關(guān)系曲線，以曲線位移突變?yōu)榕R界破壞；Lu等[7]從能量演化角度出發(fā)，認(rèn)為邊坡總體勢(shì)能的突變可作為整體失穩(wěn)判據(jù)；王培濤等[8]選取黃金分割點(diǎn)作為強(qiáng)度折減顆粒離散元法的折減比例，加速了邊坡強(qiáng)度折減的計(jì)算過程。

通過上述研究，應(yīng)用強(qiáng)度折減顆粒離散元法分析邊坡穩(wěn)定性的能力有所提升，但其在實(shí)際工程中應(yīng)用仍較少。原因在于，顆粒離散元賦予顆粒一系列的細(xì)觀參數(shù)模擬巖土體，應(yīng)用強(qiáng)度折減顆粒離散元法評(píng)價(jià)邊坡穩(wěn)定性時(shí)，如何快速調(diào)整顆粒細(xì)觀參數(shù)取值以對(duì)應(yīng)巖土體抗剪強(qiáng)度折減的問題，仍有待解決[9]。針對(duì)給定的抗剪強(qiáng)度值，采用雙軸壓縮等數(shù)值試驗(yàn)?zāi)Ｐ停貜?fù)調(diào)整顆粒細(xì)觀參數(shù)取值，使得數(shù)值試驗(yàn)測(cè)定的抗剪強(qiáng)度與給定值相符，是運(yùn)用較廣的顆粒離散元細(xì)觀參數(shù)標(biāo)定手段，稱為試錯(cuò)法[10]，但其存在費(fèi)時(shí)費(fèi)力、盲目性大和精確度低等局限[11]。執(zhí)行強(qiáng)度折減法時(shí)，抗剪強(qiáng)度將經(jīng)歷多次折減，相應(yīng)地需要多次標(biāo)定細(xì)觀參數(shù)，試錯(cuò)法的能力無法滿足這一需求。為解決此問題，部分研究將顆粒之間的法向、切向黏結(jié)強(qiáng)度和摩擦系數(shù)同步折減[12]，以表征抗剪強(qiáng)度的折減。由于顆粒細(xì)觀參數(shù)與抗剪強(qiáng)度之間為復(fù)雜的非線性關(guān)系[13-14]，故該手段的合理性存疑。

BР神經(jīng)網(wǎng)絡(luò)能很好地映射輸入?yún)?shù)與輸出參數(shù)之間的非線性關(guān)系，在開采沉降預(yù)測(cè)[15]、煤體瓦斯?jié)B透率預(yù)測(cè)[16]、巖體力學(xué)參數(shù)反演[17]、滑坡易發(fā)性評(píng)價(jià)[18]和地質(zhì)災(zāi)害早期預(yù)警[19]等工程實(shí)踐中的各類問題上得到成功應(yīng)用。本文擬借助BР神經(jīng)網(wǎng)絡(luò)表達(dá)顆粒離散元細(xì)觀參數(shù)與巖土體抗剪強(qiáng)度之間的數(shù)學(xué)關(guān)系，建立根據(jù)抗剪強(qiáng)度目標(biāo)值快速標(biāo)定顆粒細(xì)觀參數(shù)的有效方法，以滿足強(qiáng)度折減時(shí)細(xì)觀參數(shù)快速調(diào)整的要求。

區(qū)別于已有顆粒離散元法研究多基于Itasca公司開發(fā)的顆粒流程序РFC的做法，本文采用南京大學(xué)自主研發(fā)的國產(chǎn)顆粒離散元軟件MatDEM[20]。Mat-DEM可根據(jù)研究需要在Matlab中便捷地進(jìn)行二次開發(fā)，自2018年5月發(fā)布以來受到國內(nèi)外學(xué)者的高度關(guān)注，在邊坡失穩(wěn)模式預(yù)測(cè)[21]、高速遠(yuǎn)程滑坡運(yùn)動(dòng)過程分析[22]、泥石流成災(zāi)機(jī)理揭示[23]、光纜與土體的協(xié)調(diào)變形模擬[24]等問題上已取得較好的應(yīng)用成果。新時(shí)代科技競爭日益激烈的背景下，國產(chǎn)科學(xué)計(jì)算軟件開發(fā)和應(yīng)用的重要性日益凸顯。基于MatDEM研究強(qiáng)度折減顆粒離散元法分析邊坡穩(wěn)定性時(shí)顆粒細(xì)觀參數(shù)的標(biāo)定策略，可有效助力國產(chǎn)離散元軟件的應(yīng)用推廣。

1 顆粒細(xì)觀參數(shù)和抗剪強(qiáng)度測(cè)定

1.1 MatDEM的顆粒細(xì)觀參數(shù)

MatDEM通過膠接一系列具有特定力學(xué)性質(zhì)的顆粒來構(gòu)建巖土體模型，與РFC等其他顆粒離散元軟件相同，其細(xì)觀參數(shù)包括物理幾何參數(shù)和接觸模型參數(shù)兩類。物理幾何參數(shù)一般由具體模擬工程對(duì)象的基本特性決定[11]，包括數(shù)值模型尺寸、顆粒直徑、顆粒集合孔隙率、顆粒密度等。其中，顆粒直徑的取值一方面影響模型的宏觀力學(xué)行為，另一方面關(guān)系到模型的計(jì)算規(guī)模和效率，因此必須慎重選取。宿輝等[25]統(tǒng)計(jì)研究后指出，當(dāng)特征長度比L/R（模型長邊尺寸L和顆粒平均半徑R的比值）達(dá)到200以上時(shí)，模型的宏觀力學(xué)性質(zhì)與破壞模式將不受顆粒尺寸效應(yīng)的影響。

MatDEM顆粒之間依靠彈簧相互接觸產(chǎn)生力的作用，如圖1所示。接觸模型參數(shù)包括法向剛度Kn、切向剛度Ks、斷裂位移Xb、切向初始抗剪力Fs0和摩擦系數(shù)μp。

圖1 MatDEM接觸模型示意圖Fig.1 Schematic diagram of contact model in MatDEM

顆粒間法向力Fn與法向相對(duì)位移Xn關(guān)系為：

初始狀態(tài)時(shí)顆粒與相鄰顆粒相互連接，受法向彈簧的拉壓作用。當(dāng)兩顆粒之間的Xn超過斷裂位移Xb時(shí)，法向連接斷裂，顆粒間的拉力消失，僅可承受壓力作用。

顆粒間切向力Fs與切向相對(duì)位移Xs關(guān)系為：

切向連接的斷裂滿足摩爾庫倫屈服準(zhǔn)則：

式中：Fs,max為抗剪力；切向初始抗剪力Fs0為沒有施加法向壓力時(shí)顆粒間的抗剪力，類似于巖土體的黏聚力。法向壓力越大，抗剪力越大。當(dāng)切向力Fs超過抗剪力時(shí)，切向連接斷裂，顆粒間僅?；瑒?dòng)摩擦力-μpFn；對(duì)應(yīng)真實(shí)顆粒之間的膠結(jié)斷裂，法向連接斷裂時(shí)，切向連接也相應(yīng)斷開，反之亦然。

邊坡顆粒離散元數(shù)值模型建立完畢后，模型物理幾何參數(shù)將保持不變。應(yīng)用強(qiáng)度折減法評(píng)價(jià)穩(wěn)定性時(shí)，必須對(duì)接觸模型參數(shù)進(jìn)行調(diào)整，以對(duì)應(yīng)宏觀的巖土體抗剪強(qiáng)度折減。因此，本文細(xì)觀參數(shù)標(biāo)定主要針對(duì)接觸模型參數(shù)進(jìn)行。

1.2 基于雙軸壓縮模型的抗剪強(qiáng)度測(cè)定

顆粒離散元方法中，雙軸壓縮模型常用于測(cè)定顆粒膠結(jié)形成的巖土體抗剪強(qiáng)度。圖2（a）所示的MatDEM雙軸壓縮模型數(shù)值試樣的尺寸為150 mm×300 mm，由5 000個(gè)顆粒組成，直徑為2.5～3.6 mm。模型的頂部和底部為剛性簇單元組成的剛性壓力板，側(cè)壁為柔性簇單元組成的柔性壓力板，保持圍壓不變的情況下可輕微移動(dòng)變形。先施加兩側(cè)壓力板上的圍壓，再給頂部壓力板逐步施加垂直向下位移，單步施加量0.05 mm，平衡后記錄頂部壓力板受力情況，直至試樣最終破壞，如圖2（b）所示。

圖2 雙軸壓縮試驗(yàn)示意圖Fig.2 Schematic diagram of biaxial compression numerical test

以細(xì)觀參數(shù)Kn、Ks、Xb、Fs0和μp分別取值2.9×107N/mm、1.1×107N/mm、4.2×10-4m、6.4×107Рa和0.12為例，在25、50和75 kРa這3個(gè)不同圍壓等級(jí)下進(jìn)行雙軸壓縮試驗(yàn)，得到應(yīng)力-應(yīng)變曲線，如圖3（a）所示。由圖3（a）可見，軸向應(yīng)力達(dá)到第1個(gè)峰值后，應(yīng)力不變情況下，曲線存在明顯應(yīng)變?cè)鲩L段。以破壞時(shí)的軸壓和圍壓為依據(jù)，做莫爾圓并繪制包絡(luò)線，如圖3（b）所示，可測(cè)定上述細(xì)觀參數(shù)設(shè)定下形成的巖土體抗剪強(qiáng)度值為：內(nèi)摩擦角φ為31.2°，黏聚力C為31.9 kРa。

圖3 基于雙軸壓縮試驗(yàn)的抗剪強(qiáng)度測(cè)定Fig.3 Determination of the shear strength based on biaxial compression numerical test

試錯(cuò)法重復(fù)調(diào)整顆粒細(xì)觀參數(shù)使得雙軸壓縮試驗(yàn)測(cè)定的抗剪強(qiáng)度符合預(yù)期，但在效率和精度方面均存在局限。強(qiáng)度折減法評(píng)價(jià)邊坡穩(wěn)定性，巖土體抗剪強(qiáng)度值根據(jù)折減系數(shù)不斷調(diào)整。試錯(cuò)法標(biāo)定細(xì)觀參數(shù)的能力不足，成為限制強(qiáng)度折減顆粒離散元法在邊坡穩(wěn)定性評(píng)價(jià)中廣泛應(yīng)用的技術(shù)瓶頸。

必須說明的是，真實(shí)巖土體的變形及強(qiáng)度指標(biāo)往往具備較強(qiáng)的非線性，邊坡不同部位的應(yīng)力水平不同，其抗剪強(qiáng)度指標(biāo)可能并不一致。采用數(shù)值模擬方法分析邊坡穩(wěn)定性時(shí)，引入隨應(yīng)力水平相應(yīng)變化的抗剪強(qiáng)度無疑更為合理，這方面已有學(xué)者開展了探索性的工作[26]?？紤]本文主要擬解決強(qiáng)度折減顆粒離散元法分析邊坡穩(wěn)定性時(shí)試錯(cuò)法標(biāo)定細(xì)觀參數(shù)能力不足的問題，在細(xì)觀參數(shù)標(biāo)定時(shí)，暫未考慮邊坡不同部位應(yīng)力水平差異對(duì)抗剪強(qiáng)度指標(biāo)及其他參數(shù)的影響。

2 基于BP神經(jīng)網(wǎng)絡(luò)的細(xì)觀參數(shù)逆向迭代修正標(biāo)定

2.1 細(xì)觀參數(shù)標(biāo)定的BP神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)設(shè)計(jì)

BР神經(jīng)網(wǎng)絡(luò)是應(yīng)用最廣泛的神經(jīng)網(wǎng)絡(luò)框架之一，通過模擬大腦對(duì)信息處理的機(jī)制尋找輸入信息與輸出信息之間的非線性關(guān)系，由輸入層、隱藏層、輸出層3部分組成，如圖4所示。

圖4 BP神經(jīng)網(wǎng)絡(luò)基本結(jié)構(gòu)Fig.4 Typical structure of BP neural network

BР神經(jīng)網(wǎng)絡(luò)前向傳播時(shí)神經(jīng)元的輸出值為：

式中：l為隱藏層層數(shù)；為第l層、第i個(gè)神經(jīng)元的輸出值；為第l層、第i個(gè)神經(jīng)元的激活函數(shù)；kl-1為第l-1層神經(jīng)元個(gè)數(shù)；為計(jì)算第l層、第i個(gè)神經(jīng)元時(shí)，第l-1層第j個(gè)神經(jīng)元的權(quán)重；為第l層、第i個(gè)神經(jīng)元的偏置項(xiàng)?；跇颖編煊?xùn)練完畢的BР神經(jīng)網(wǎng)絡(luò)，其前向傳播過程事實(shí)上可視為關(guān)于輸入信息x和輸出信息y的非線性函數(shù)：

式中，w和b分別代表所有權(quán)重和偏置信息的集合。

周瑜等[27]嘗試使用BР神經(jīng)網(wǎng)絡(luò)解決顆粒離散元軟件РFC的細(xì)觀參數(shù)標(biāo)定問題，神經(jīng)網(wǎng)絡(luò)構(gòu)建時(shí)以彈性模量、單軸抗壓強(qiáng)度和泊松比等3個(gè)力學(xué)指標(biāo)為輸入，以РFC顆粒的黏結(jié)半徑系數(shù)等4個(gè)細(xì)觀參數(shù)為輸出，樣本庫為400組顆粒的細(xì)觀參數(shù)和數(shù)值試驗(yàn)獲取的巖土體力學(xué)指標(biāo)。然而，其在執(zhí)行“給定巖土體力學(xué)指標(biāo)→神經(jīng)網(wǎng)絡(luò)輸出顆粒細(xì)觀參數(shù)→執(zhí)行數(shù)值試驗(yàn)獲取力學(xué)指標(biāo)→獲取值與給定值比較”的精度測(cè)試后發(fā)現(xiàn)， 10組測(cè)試方案中部分方案精度明顯不符要求。原因在于，BР神經(jīng)網(wǎng)絡(luò)中以巖土體力學(xué)參數(shù)為輸入信息x，以顆粒細(xì)觀參數(shù)為輸出信息y，此時(shí)輸入信息維度少于輸出信息，形成了數(shù)學(xué)上的欠定問題[28]；對(duì)于欠定問題，BР神經(jīng)網(wǎng)絡(luò)必須引入額外約束條件將其轉(zhuǎn)化為適定或超定問題后才能保障求解精度[29]。

巖土體抗剪強(qiáng)度標(biāo)定MatDEM顆粒細(xì)觀參數(shù)時(shí)，已知信息為內(nèi)摩擦角和黏聚力兩個(gè)指標(biāo)，待確定的顆粒細(xì)觀參數(shù)為法向剛度Kn等5個(gè)指標(biāo)。若按照常規(guī)方式將已知信息作為輸入，以待定信息作為輸出同樣會(huì)導(dǎo)致欠定問題。為避免引入額外約束條件且保障神經(jīng)網(wǎng)絡(luò)前向傳播數(shù)值精度，建立BР神經(jīng)網(wǎng)絡(luò)時(shí)，輸入信息x為顆粒細(xì)觀參數(shù)，輸出信息y為內(nèi)摩擦角和黏聚力，其結(jié)構(gòu)如圖5所示。神經(jīng)網(wǎng)絡(luò)的輸入層神經(jīng)元為5個(gè)，輸出層神經(jīng)元為兩個(gè)；為保障數(shù)值精度，隱藏層設(shè)定為3層，每層有8個(gè)神經(jīng)元。

圖5 細(xì)觀參數(shù)和抗剪強(qiáng)度聯(lián)系的神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)Fig.5 Neural network from micro parameters to shear strength

BР神經(jīng)網(wǎng)絡(luò)的樣本為MatDEM顆粒細(xì)觀參數(shù)和雙軸壓縮試驗(yàn)獲取的相應(yīng)巖土體抗剪強(qiáng)度。參考細(xì)觀參數(shù)的取值范圍[20]，F(xiàn)s0設(shè)置5個(gè)水平，Xb設(shè)置3個(gè)水平，μp、Kn和Ks統(tǒng)一設(shè)置4個(gè)水平，各水平取值見表1。全面正交組合共計(jì)形成960組細(xì)觀參數(shù)進(jìn)行雙軸壓縮數(shù)值試驗(yàn)。顆粒密度ρ取值1 900 kg/m3，顆粒直徑取2.5～3.6 mm。

表1 細(xì)觀參數(shù)水平設(shè)定Tab.1 Level setting of micro parameters

將960組MatDEM顆粒細(xì)觀參數(shù)和雙軸壓縮數(shù)值試驗(yàn)獲取的抗剪強(qiáng)度值作為神經(jīng)網(wǎng)絡(luò)樣本數(shù)據(jù)庫，執(zhí)行BР神經(jīng)網(wǎng)絡(luò)的學(xué)習(xí)訓(xùn)練。首先，對(duì)樣本數(shù)據(jù)進(jìn)行歸一化處理，避免各項(xiàng)數(shù)據(jù)量綱/量級(jí)的不一致對(duì)結(jié)果產(chǎn)生影響。歸一化計(jì)算公式為：

式中，X、Xnorm分別為樣本歸一化之前和歸一化之后的數(shù)據(jù)，Xmin和Xmax分別為數(shù)據(jù)最小值和最大值。

MatDEM是在Matlab框架下編制的，其可便捷地在Matlab中進(jìn)行二次開發(fā)；同時(shí)，Matlab神經(jīng)網(wǎng)絡(luò)工具箱又為BР神經(jīng)網(wǎng)絡(luò)的學(xué)習(xí)訓(xùn)練提供了大量的預(yù)置函數(shù)，因此，所有程序部分可在Matlab中實(shí)現(xiàn)有效集成。BР神經(jīng)網(wǎng)絡(luò)學(xué)習(xí)訓(xùn)練時(shí)，訓(xùn)練函數(shù)選用trainlm函數(shù)，輸出函數(shù)選用pureline函數(shù)，激活函數(shù)選擇tansig函數(shù)，誤差函數(shù)選用均方誤差RMSE。

2.2 細(xì)觀力學(xué)參數(shù)的逆向迭代修正標(biāo)定

BР神經(jīng)網(wǎng)絡(luò)訓(xùn)練完畢后，式（5）中的所有權(quán)重和偏置信息均確定，其事實(shí)上也定義了黏聚力、內(nèi)摩擦角與細(xì)觀參數(shù)之間的復(fù)雜函數(shù)關(guān)系：

式中，f1和f2代表訓(xùn)練好的BР神經(jīng)網(wǎng)絡(luò)中由細(xì)觀參數(shù)逐步運(yùn)算至黏聚力和內(nèi)摩擦角的過程，其直觀表達(dá)式由于權(quán)重、偏置項(xiàng)和激活函數(shù)等的引入很難直接給出，但其真實(shí)存在。

在BР神經(jīng)網(wǎng)絡(luò)框架下，當(dāng)給定輸入信息時(shí)，可進(jìn)行該輸入信息的梯度計(jì)算，因此，仿照非線性方程數(shù)學(xué)求解，可基于式（7）和（8）進(jìn)行給定抗剪強(qiáng)度目標(biāo)值下細(xì)觀參數(shù)取值的數(shù)學(xué)推演。首先，給予一組細(xì)觀參數(shù)初值，使用前向傳播獲取抗剪強(qiáng)度值，將該值與目標(biāo)值的誤差定義為：

式中，C0和φ0分別為黏聚力和內(nèi)摩擦角的目標(biāo)值。設(shè)定誤差限為RMSE0，若誤差不滿足要求，沿誤差梯度下降方向調(diào)整細(xì)觀參數(shù)，直至誤差滿足要求，如圖6所示。圖6中：n為循環(huán)次數(shù)；lr為學(xué)習(xí)率，其引進(jìn)目的為執(zhí)行沿誤差梯度下降方向的步長搜索逼近，本文取0.01。區(qū)別于周瑜等[27]的做法，本文的BР神經(jīng)網(wǎng)絡(luò)以顆粒細(xì)觀參數(shù)為輸入，以巖土體抗剪強(qiáng)度為輸出，通過反推的方式獲取巖土體抗剪強(qiáng)度給定后的顆粒細(xì)觀參數(shù)，故稱之為逆向標(biāo)定。

圖6 BP神經(jīng)網(wǎng)絡(luò)逆向標(biāo)定細(xì)觀參數(shù)Fig.6 Micro parameters determined by BP neural network

對(duì)于逆向標(biāo)定得到的細(xì)觀參數(shù)，執(zhí)行雙軸壓縮數(shù)值試驗(yàn)，檢查抗剪強(qiáng)度與給定值的吻合程度，即精度檢查。精度檢查的結(jié)果可能并不理想，原因在于BР神經(jīng)網(wǎng)絡(luò)樣本庫質(zhì)量未能達(dá)到要求。訓(xùn)練完畢的BР神經(jīng)網(wǎng)絡(luò)前向傳播過程可視為雙軸壓縮數(shù)值試驗(yàn)的一種替代，樣本庫的樣本數(shù)量越多，其替代的可靠程度越高。為保障效率，神經(jīng)網(wǎng)絡(luò)的樣本庫的樣本數(shù)量有限，故前向傳播過程替代雙軸壓縮數(shù)值試驗(yàn)的可靠度未必足夠，導(dǎo)致逆向標(biāo)定所得細(xì)觀參數(shù)精度不足。

對(duì)于非線性方程求解而言，已訓(xùn)練完畢的BР神經(jīng)網(wǎng)絡(luò)雖然精度不足，但提供了粗略的搜索方向，因此，逆向標(biāo)定獲取的顆粒細(xì)觀參數(shù)將比給定的細(xì)觀參數(shù)初值更加趨近于最終解。鑒于此，當(dāng)精度檢查結(jié)果不理想時(shí)，將逆向標(biāo)定得到的細(xì)觀參數(shù)與其雙軸壓縮數(shù)值試驗(yàn)獲取的抗剪強(qiáng)度作為一組新樣本添加到已有樣本庫中，重新進(jìn)行BР神經(jīng)網(wǎng)絡(luò)的學(xué)習(xí)訓(xùn)練，此過程可稱為“樣本修正”。BР神經(jīng)網(wǎng)絡(luò)經(jīng)過重新學(xué)習(xí)訓(xùn)練后，式（7）和（8）中的函數(shù)f1和f2將相應(yīng)地調(diào)整變化，以前次標(biāo)定獲得的細(xì)觀參數(shù)為初值，再次進(jìn)行逆向標(biāo)定和精度檢查。重復(fù)執(zhí)行逆向標(biāo)定—精度檢查—樣本修正的過程，直至精度檢查順利通過，則顆粒細(xì)觀參數(shù)的標(biāo)定過程完成，整個(gè)流程如圖7所示。

圖7 細(xì)觀參數(shù)標(biāo)定的完整流程Fig.7 Whole calibration process of micro parameters

該流程中，每一輪的逆向標(biāo)定均以上一輪標(biāo)定獲得的細(xì)觀參數(shù)為初值，故顆粒細(xì)觀參數(shù)的最終取值獲取類似于數(shù)值求解中的迭代方式，因此，將這一總體流程稱為顆粒細(xì)觀參數(shù)的逆向迭代修正標(biāo)定。

2.3 細(xì)觀參數(shù)標(biāo)定策略的有效性測(cè)試

為測(cè)試逆向迭代修正標(biāo)定策略的有效性，給定抗剪強(qiáng)度目標(biāo)值C=3.0 kРa、φ=19.6°，執(zhí)行標(biāo)定流程以獲取顆粒細(xì)觀參數(shù)。經(jīng)過4輪逆向標(biāo)定后，細(xì)觀參數(shù)精度檢查通過，每輪逆向標(biāo)定獲取的細(xì)觀參數(shù)和雙軸壓縮數(shù)值試驗(yàn)獲取的抗剪強(qiáng)度值見表2，抗剪強(qiáng)度與目標(biāo)值的誤差如圖8所示。

表2 每輪標(biāo)定得到的細(xì)觀參數(shù)及測(cè)定的抗剪強(qiáng)度Tab.2 Micro parameters obtained in each round and its resulted shear strength

圖8 兩種神經(jīng)網(wǎng)絡(luò)誤差對(duì)比Fig.8 Comparison on results of the two neural networks

以抗剪強(qiáng)度給定值為標(biāo)準(zhǔn)，細(xì)觀參數(shù)逆向迭代修正標(biāo)定執(zhí)行過程中，第1輪標(biāo)定的黏聚力相對(duì)誤差為27.52%，摩擦系數(shù)相對(duì)誤差為59.70%；迭代1次后黏聚力相對(duì)誤差降至4.19%，摩擦系數(shù)相對(duì)誤差降至15.21%；4次迭代后，黏聚力相對(duì)誤差降至0.70%，摩擦系數(shù)相對(duì)誤差降至0.36%。因此，該方法在顆粒細(xì)觀參數(shù)的標(biāo)定上具有可靠的數(shù)值精度。同時(shí)，必須注意的是，前面2輪標(biāo)定得到的細(xì)觀參數(shù)精度不足。但是，逐輪的樣本修正后，標(biāo)定結(jié)果的數(shù)值精度將迅速提升。其原因在于，如第2.2節(jié)所述，每輪逆向標(biāo)定獲取的顆粒細(xì)觀參數(shù)將比上輪標(biāo)定獲取的細(xì)觀參數(shù)更趨近于最終解，那么，樣本修正可理解為在最終解附近區(qū)域?qū)ι窠?jīng)網(wǎng)絡(luò)不斷加密。

為說明BР神經(jīng)網(wǎng)絡(luò)解決數(shù)學(xué)欠定問題時(shí)的局限性，圖8補(bǔ)充繪制了欠定BР神經(jīng)網(wǎng)絡(luò)輸出結(jié)果的誤差情況。該欠定BР神經(jīng)網(wǎng)絡(luò)以抗剪強(qiáng)度為輸入，以顆粒細(xì)觀參數(shù)為輸出，隱藏層層數(shù)和每層神經(jīng)元的布置仍如圖5所示。采用同樣的數(shù)據(jù)樣本庫，BР神經(jīng)網(wǎng)絡(luò)訓(xùn)練完畢后，直接輸入抗剪強(qiáng)度給定值得到顆粒細(xì)觀參數(shù)。用欠定BР神經(jīng)網(wǎng)絡(luò)輸出的顆粒細(xì)觀參數(shù)進(jìn)行雙軸壓縮模型試驗(yàn)，黏聚力的相對(duì)誤差為44.67%，摩擦系數(shù)相對(duì)誤差為45.85%。

3 雙軸壓縮數(shù)值模型并行試驗(yàn)技術(shù)

BР神經(jīng)網(wǎng)絡(luò)標(biāo)定細(xì)觀參數(shù)的全過程中，數(shù)據(jù)樣本庫建立是一項(xiàng)基礎(chǔ)而繁瑣的工作。為加快樣本庫建立過程，開發(fā)雙軸壓縮數(shù)值模型并行試驗(yàn)技術(shù)，以充分發(fā)揮MatDEM中GРU計(jì)算模式的能力[30]。Mat-DEM中有GРU和CРU兩種計(jì)算模式，執(zhí)行不同功能時(shí)調(diào)用不同的計(jì)算模式。一般來說，數(shù)值模型建立和材料參數(shù)賦值為CРU計(jì)算模式。受力后的顆粒平衡迭代計(jì)算則會(huì)根據(jù)顆粒數(shù)量在GPU和CPU兩種模式中自動(dòng)切換，單元數(shù)量較少時(shí)選擇CPU模式，單元數(shù)量增加到一定程度則選擇GPU模式。單個(gè)雙軸壓縮數(shù)值模型試驗(yàn)盒中，顆粒數(shù)量一般在5 000左右，無法發(fā)揮GPU模式計(jì)算優(yōu)勢(shì)。

本文對(duì)MatDEM進(jìn)行二次開發(fā)，實(shí)現(xiàn)雙軸壓縮數(shù)值模型的并行試驗(yàn)；同時(shí)，構(gòu)建多個(gè)雙軸壓縮模型，每個(gè)模型的細(xì)觀參數(shù)單獨(dú)賦值，壓力板位移同步施加，顆粒受力后的平衡迭代同步進(jìn)行，壓力板數(shù)據(jù)同步記錄。此時(shí)，參與同步計(jì)算的顆粒數(shù)量將達(dá)數(shù)十萬，利用專業(yè)工作站的高性能GPU，可有效提高計(jì)算效率。圖9為54個(gè)雙軸壓縮數(shù)值模型的排列分布情況和并行試驗(yàn)獲取的試樣破壞時(shí)的位移云圖。

圖9 54個(gè)雙軸壓縮模型的排列分布和位移云圖Fig.9 Distribution and displacement diagrams of 54 biaxial compression models

試驗(yàn)研究發(fā)現(xiàn)，雙軸壓縮數(shù)值模型并行試驗(yàn)所需要的時(shí)間一方面取決于GPU的計(jì)算能力，另一方面取決于雙軸壓縮數(shù)值模型的排列分布。當(dāng)雙軸壓縮數(shù)值模型的排列分布呈現(xiàn)圖9所示的正方形布局時(shí)，并行試驗(yàn)的效率最高。其原因在于，以顆粒離散元為理論基礎(chǔ)的MatDEM，執(zhí)行逐時(shí)間步的增量計(jì)算，必然涉及動(dòng)力系統(tǒng)的能量平衡和耗散問題。默認(rèn)阻尼下，MatDEM數(shù)值模型完成一次標(biāo)準(zhǔn)平衡，即應(yīng)力波從模型一端傳至另一端所需的時(shí)間Ttotal為：

式中：Npack表示模型最長維度上堆積的顆粒數(shù)；ΔT為動(dòng)力計(jì)算采用的時(shí)間增量步步長，取值需遠(yuǎn)小于彈簧簡諧振動(dòng)周期，以精確模擬數(shù)值模型彈性變形過程。多個(gè)雙軸壓縮數(shù)值模型并行試驗(yàn)時(shí)，Mat-DEM會(huì)將同步試驗(yàn)的數(shù)值模型視為整體。雙軸壓縮數(shù)值模型總個(gè)數(shù)一定時(shí)，采用正方形布局時(shí)，整體系統(tǒng)在兩個(gè)維度上的顆粒數(shù)基本接近，效率最高。

硬件配置Intel（R） Xeon（R） E-2224@3.40 GHz CPU，NIVIDIA Quadro RTX 4000顯卡的情況下，單個(gè)雙軸壓縮模型的試驗(yàn)耗時(shí)約50 min。同樣配置條件下，54個(gè)雙軸壓縮模型并行試驗(yàn)總耗時(shí)約480 min，節(jié)約82%的計(jì)算時(shí)間。因此，采用并行試驗(yàn)技術(shù)進(jìn)行雙軸壓縮數(shù)值模型試驗(yàn)，可充分發(fā)揮GPU計(jì)算模式的能力。

4 邊坡穩(wěn)定性算例分析

以澳大利亞計(jì)算機(jī)應(yīng)用協(xié)會(huì)（簡稱ACADS）的兩個(gè)邊坡經(jīng)典考題為例，采用MatDEM執(zhí)行邊坡穩(wěn)定性的強(qiáng)度折減法分析，考察本文提出的顆粒細(xì)觀參數(shù)標(biāo)定策略的實(shí)用性。穩(wěn)定性分析時(shí)，臨界破壞準(zhǔn)則參照田雷等[6]的做法，記邊坡模型中所有顆粒的位移平均值為模型平均位移S，以其發(fā)生突變?yōu)檫吰逻_(dá)到極限狀態(tài)的判據(jù)。

4.1 算例1：均質(zhì)土質(zhì)邊坡

圖10為均質(zhì)土質(zhì)邊坡，土體力學(xué)參數(shù)為C=3.0 kPa，φ=19.6°，重度γ=20.0 kN/m3。ACADS給出的安全系數(shù)推薦值為1.00，對(duì)應(yīng)的滑動(dòng)范圍如圖10橙色區(qū)域所示。

圖10 算例1模型及ACADS推薦解的滑動(dòng)范圍Fig.10 Model of Example 1 and the failure region by ACADS

在MatDEM中建立該邊坡的顆粒離散元模型如圖11所示，模型含82 698個(gè)直徑67～96 mm的顆粒。此時(shí)，根據(jù)特征長度比L/R要求[25]將用于測(cè)定抗剪強(qiáng)度的雙軸壓縮數(shù)值模型尺寸調(diào)整為4 000 mm×8 000 mm，以避免顆粒尺寸效應(yīng)影響細(xì)觀參數(shù)標(biāo)定。

圖11 邊坡的MatDEM模型和臨界失穩(wěn)狀態(tài)位移Fig.11 MatDEM model and its critical displacements

分析穩(wěn)定性時(shí)，初始折減系數(shù)設(shè)為0.70；若邊坡處于穩(wěn)定狀態(tài)，則將折減系數(shù)增加0.10繼續(xù)計(jì)算，折減系數(shù)達(dá)到1.10時(shí)邊坡失穩(wěn)破壞。在折減系數(shù)0.95到1.10范圍內(nèi)進(jìn)一步細(xì)化分析，最終確定邊坡穩(wěn)定系數(shù)為1.02。邊坡模型平均位移S隨著強(qiáng)度折減系數(shù)Fsr增大的變化曲線如圖12所示。由圖12可知：Fsr達(dá)到1.02前，S基本呈水平直線；數(shù)據(jù)檢查發(fā)現(xiàn)，S隨著Fsr增加仍略有增加，表明抗剪強(qiáng)度下降會(huì)影響顆粒平衡后的位置，宏觀上體現(xiàn)為變形輕微增長；Fsr達(dá)到1.02后，S值隨Fsr增加而急劇變大，表明邊坡進(jìn)入臨界失穩(wěn)狀態(tài)。邊坡達(dá)到臨界失穩(wěn)狀態(tài)時(shí)的位移，如圖11所示。此時(shí)邊坡自坡腳到坡頂已出現(xiàn)明顯貫通的滑動(dòng)面，邊坡沿該滑動(dòng)面將開始發(fā)生滑動(dòng)。

圖12 算例1的S-Fsr曲線Fig.12 S-Fsr curves for Example 1

強(qiáng)度折減過程中，抗剪強(qiáng)度不同折減系數(shù)折減后對(duì)應(yīng)的細(xì)觀參數(shù)見表3。由表3可知，由于抗剪強(qiáng)度與細(xì)觀參數(shù)之間的非線性映射關(guān)系非常復(fù)雜，抗剪強(qiáng)度的折減很難通過直接對(duì)細(xì)觀參數(shù)進(jìn)行折減來實(shí)現(xiàn)。

表3 不同折減系數(shù)下土體抗剪強(qiáng)度與細(xì)觀力學(xué)參數(shù)Tab.3 Micro parameters and soil shear strength corresponding to reduction coefficients

4.2 算例2：堆石壩模型

算例2為一由堆石體、過渡層、心墻和反濾層組成的堆石壩，壩高146 m，壩底寬648 m，其形態(tài)如圖13所示。

圖13 堆石壩模型及ACADS推薦解的滑動(dòng)區(qū)域Fig.13 A rock fill dam and the failure region by ACADS

壩體模型各種材料的力學(xué)參數(shù)見表4，由于算例2是用于穩(wěn)定性分析，其參數(shù)與實(shí)際工程中各種材料的參數(shù)取值存在一定差異。ACADS對(duì)該堆石壩模型給出的安全系數(shù)推薦值為1.976，對(duì)應(yīng)的滑動(dòng)面如圖13所示。

表4 堆石壩模型各材料力學(xué)參數(shù)Tab.4 Mechanical properties of materials in rockfill dam model

在MatDEM中建立該模型的顆粒離散元模型如圖14所示，模型含89 768個(gè)直徑為0.33～0.48 m的顆粒。顆粒細(xì)觀參數(shù)標(biāo)定時(shí)，按照特征長度比L/R的要求[25]，將雙軸壓縮數(shù)值模型尺寸調(diào)整為20 m×40 m。

圖14 堆石壩的MatDEM模型Fig.14 MatDEM model for the rock fill dam

采用強(qiáng)度折減法分析該堆石壩模型的穩(wěn)定性，模型的平均位移S和強(qiáng)度折減系數(shù)Fsr關(guān)系曲線如圖15（a）所示。以S突變?yōu)槟Ｐ团R界失穩(wěn)判據(jù)，確定該堆石壩安全系數(shù)為1.975。模型達(dá)到臨界失穩(wěn)狀態(tài)時(shí)的位移如圖15（b）所示。圖15（a）顯示，即便模型處于穩(wěn)定狀態(tài)，模型平均位移S已超過2 m。其原因在于，顆粒離散元邊坡模型建立完畢之后，由于存在重力勢(shì)能，顆粒自平衡過程中會(huì)導(dǎo)致模型產(chǎn)生一定的變形，類似于重力作用下的沉降變形。事實(shí)上，算例1的邊坡模型處于穩(wěn)定狀態(tài)時(shí)，也存在0.002 m的模型平均位移S。算例2模型平均位移S偏大，其原因主要為，施加于模型上的邊界約束條件，堆石壩模型僅在底部邊界設(shè)置位移約束，而算例1邊坡模型在底部和左右兩側(cè)邊界均設(shè)置位移約束，因此堆石壩模型重力效應(yīng)可充分發(fā)揮，變形相對(duì)較大；此外，堆石壩模型高度達(dá)146 m，大尺寸規(guī)模下，顆粒自平衡過程中的模型變形也相對(duì)較大。

圖15 堆石壩的安全系數(shù)結(jié)果和臨界失穩(wěn)狀態(tài)位移Fig.15 Safety factors and the critical displacements for therock fill dam

強(qiáng)度折減過程中，材料顆粒細(xì)觀參數(shù)隨著折減系數(shù)變化而調(diào)整的情況見表5。

表5 不同折減系數(shù)下堆石壩模型中材料的抗剪強(qiáng)度與細(xì)觀力學(xué)參數(shù)Tab.5 Micro parameters and shear strength of material in rockfill dam model corresponding to reduction coefficients

兩個(gè)算例的執(zhí)行情況顯示，本文提出的細(xì)觀參數(shù)標(biāo)定策略，可滿足MatDEM等顆粒離散元軟件采用強(qiáng)度折減法分析邊坡穩(wěn)定性時(shí)，迅速調(diào)整顆粒細(xì)觀參數(shù)取值，以對(duì)應(yīng)巖土體抗剪強(qiáng)度折減的需求。安全系數(shù)結(jié)果方面：算例1安全系數(shù)為1.02，與ACADS推薦的安全系數(shù)相差0.02；算例2安全系數(shù)為1.975，與ACADS推薦的系數(shù)相差0.001?；瑒?dòng)面結(jié)果方面：MatDEM獲取的臨界失穩(wěn)狀態(tài)下模型的滑動(dòng)區(qū)域與ACADS推薦的滑動(dòng)范圍基本吻合?？傮w上，MatDEM獲取的安全系數(shù)和滑動(dòng)范圍與ACADS推薦解之間具有良好的可比性，表明采用本文策略標(biāo)定的顆粒細(xì)觀參數(shù)，在對(duì)應(yīng)不同程度折減后的巖土體抗剪強(qiáng)度方面具有滿意的數(shù)值精度。因此，在應(yīng)用顆粒離散元強(qiáng)度折減法評(píng)價(jià)邊坡穩(wěn)定性時(shí)，本文提出的顆粒細(xì)觀參數(shù)逆向迭代修正標(biāo)定策略可發(fā)揮積極作用。

5 結(jié) 論

針對(duì)以MatDEM為代表的顆粒離散元法軟件應(yīng)用強(qiáng)度折減技術(shù)評(píng)價(jià)邊坡穩(wěn)定性時(shí)的顆粒細(xì)觀參數(shù)標(biāo)定問題，以顆粒細(xì)觀參數(shù)為輸入，抗剪強(qiáng)度為輸出構(gòu)造BР神經(jīng)網(wǎng)絡(luò)，建立顆粒細(xì)觀參數(shù)的逆向迭代修正標(biāo)定方法，開發(fā)雙軸壓縮數(shù)值模型并行試驗(yàn)技術(shù)以加快神經(jīng)網(wǎng)絡(luò)樣本庫樣本的獲取，并采用ACADS的兩個(gè)邊坡穩(wěn)定性分析經(jīng)典考題驗(yàn)證了顆粒細(xì)觀參數(shù)逆向迭代修正標(biāo)定方法的能力和精度。基于研究過程中的發(fā)現(xiàn)，可得結(jié)論如下：

1）細(xì)觀參數(shù)逆向迭代修正標(biāo)定方法可有效解決顆粒離散元強(qiáng)度折減法分析邊坡穩(wěn)定性時(shí)的顆粒細(xì)觀參數(shù)標(biāo)定問題。標(biāo)定得到的細(xì)觀參數(shù)可良好對(duì)應(yīng)巖土體的抗剪強(qiáng)度，數(shù)值算例計(jì)算結(jié)果與推薦解之間的可比性好。

2）MatDEM的二次開發(fā)在Matlab程序?qū)崿F(xiàn)便捷，使得MatDEM具備很好的研究潛力和應(yīng)用價(jià)值。本文開發(fā)的雙軸壓縮數(shù)值模型并行試驗(yàn)技術(shù)有效加速了神經(jīng)網(wǎng)絡(luò)樣本庫的建立過程，可為其他小型模型的重復(fù)試驗(yàn)工作提供參考。

3）使用BР神經(jīng)網(wǎng)絡(luò)執(zhí)行參數(shù)反演問題研究時(shí)，必須注意神經(jīng)網(wǎng)絡(luò)輸入信息和輸出信息維度的大小關(guān)系，以避免形成數(shù)學(xué)上的欠定問題。