林玉祥, 林浩東, 莫品強, 褚 鋒, 莊培芝

(1.中交第三航務(wù)工程勘察設(shè)計院有限公司, 上海 200032; 2.中國礦業(yè)大學(xué) 深部巖土力學(xué)與地下工程國家重點實驗室, 江蘇 徐州 221116; 3.山東高速建設(shè)管理集團有限公司, 山東 濟南 250098;4.山東大學(xué) 齊魯交通學(xué)院, 山東 濟南 250002)

隨著連云港沿海開發(fā)規(guī)劃政策的落實,對連云港海積軟土工程特性的研究逐漸引起業(yè)內(nèi)人士的關(guān)注,其中土體滲透特性的研究為當(dāng)前熱點之一。滲透系數(shù)k是能定量反映土體滲透性強弱的重要指標,其數(shù)值的準確確定對滲透特性研究至關(guān)重要。滲透系數(shù)可由土工試驗或現(xiàn)場井孔抽注水試驗獲得,但試驗過程均比較繁瑣且耗時耗力。因此,業(yè)界學(xué)者長期致力于提出一種簡單可靠的滲透系數(shù)確定方法。

國內(nèi)外學(xué)者對粗粒土的滲透特性研究較多[1-2],基于天然土物性指標、粒徑級配曲線等,提出了大量滲透系數(shù)的經(jīng)驗公式。然而,相關(guān)公式通常不適用于細粒土,限制了其在實際工程中的應(yīng)用。為了得到適用于廣泛土類的滲透系數(shù)計算方法,部分學(xué)者利用靜力觸探可以連續(xù)反映土層性質(zhì)的特點,基于靜力觸探數(shù)據(jù)對滲透系數(shù)進行預(yù)測計算。Robertson[3]提出了滲透系數(shù)與土類指數(shù)Ic的平均關(guān)系,并給出了通過靜探數(shù)據(jù)計算Ic的公式。Elsworth等[4-5]和Chai等[6]分別結(jié)合球面流和半球面流模型,采用位錯理論進行推導(dǎo),得到了靜探測試數(shù)據(jù)與滲透系數(shù)的理論關(guān)系式。李鏡培等[7]基于前人研究成果提出了圓柱面徑向滲流模型,推導(dǎo)出了修正后的經(jīng)驗公式。此外,有學(xué)者在滲透特性研究中引入機器學(xué)習(xí)算法,取得了一定成果。許增光等[8]基于實測電阻率,采用BP神經(jīng)網(wǎng)絡(luò)對土體滲透系數(shù)進行預(yù)測,通過計算實例證明了相比經(jīng)驗公式,神經(jīng)網(wǎng)絡(luò)法的平均誤差更小。徐麗等[9]借助極限學(xué)習(xí)機建立滲透系數(shù)與水頭之間的映射關(guān)系,并結(jié)合遺傳算法構(gòu)建了滲透系數(shù)反演分析模型。Zhao等[10]使用測井?dāng)?shù)據(jù)作為輸入,訓(xùn)練了7種機器學(xué)習(xí)模型以預(yù)測珠江口盆地低滲透砂巖的滲透率,其中,XGBoost模型性能最優(yōu)。Tian等[11]通過對3種特征選擇方法和6種機器學(xué)習(xí)模型的比較分析,提出了一種將數(shù)字巖石物理模型與機器學(xué)習(xí)模型相結(jié)合的方法來改進地下多孔介質(zhì)的滲透率預(yù)測。Pham等[12]基于越南峴港項目勘察數(shù)據(jù)構(gòu)建土體滲透系數(shù)預(yù)測模型,并對采用的3種機器學(xué)習(xí)算法進行了比較分析。綜上,前人多致力于對已有經(jīng)驗公式的修正或只與機器學(xué)習(xí)算法進行簡單結(jié)合,而缺少對原始數(shù)據(jù)的合理利用和預(yù)處理分析;此外,已有經(jīng)驗公式或反演模型鮮有考慮滲透系數(shù)在水平向和垂直向的差異,而方向問題在實際應(yīng)用中往往不可忽略。

本文根據(jù)土體靜力觸探數(shù)據(jù)與其滲透特性間的關(guān)聯(lián)性,基于XGBoost機器學(xué)習(xí)算法構(gòu)建一種滲透系數(shù)反演模型。該模型通過對原始數(shù)據(jù)的預(yù)分析選取特征參數(shù)和標簽值,根據(jù)工程數(shù)據(jù)集的特點調(diào)整超參數(shù),且考慮滲透系數(shù)在水平向和垂直向的差異。最后根據(jù)連云港開展的大荷載礦石堆場深厚軟土地基加固研究項目的工程勘察數(shù)據(jù),將本文模型與BPNN模型和前人經(jīng)驗公式進行對比分析。

1 XGBoost反演模型

1.1 XGBoost算法

極限梯度提升算法(Extreme Gradient Boosting, XGBoost)(后文簡稱XGB)在普通樹模型基礎(chǔ)上改進算法以提高精度,充分調(diào)用CPU的多線程并行以使提升樹達到自身的計算極限,是近年來性能最強的機器學(xué)習(xí)算法之一[13]。

基于梯度增強決策樹算法,XGB最小化目標函數(shù)至期望范圍,預(yù)測值計算公式為:

(1)

目標函數(shù)計算公式為:

(2)

最小化目標函數(shù),經(jīng)過正則化項對算法學(xué)習(xí)權(quán)重的平滑,最終得到目標函數(shù)的最優(yōu)解:

(3)

XGB允許使用者根據(jù)實際使用場景調(diào)整損失函數(shù),對于滲透系數(shù)預(yù)測這樣的回歸問題而言,采用最小二乘損失函數(shù):

(4)

XGB算法本質(zhì)是基于梯度提升樹實現(xiàn)的集成算法,涉及參數(shù)極多。考慮到使用的數(shù)據(jù)集較為簡單,故僅對屬于集成算法和弱評估器的部分重要參數(shù)進行調(diào)參(見表1),將其它與提升模型性能無關(guān)的復(fù)雜參數(shù)設(shè)為默認值。

表1 XGBoost算法重要參數(shù)表(部分)Tab.1 Important parameters of XGBoost algorithm (part)

1.2 考慮kv和kh的差異

在研究滲透系數(shù)水平向和垂直向的差異時,Bear[14]利用徑向滲透儀裝置將兩種標準的垂直測量值(kv)與水平值(kh)進行了比較,提出了軟黏土kh/kv的可能值范圍:1～1.5(均質(zhì)沉積物),2～4(發(fā)達宏觀結(jié)構(gòu)),3～15(明顯分層)。以往的滲透系數(shù)經(jīng)驗公式往往不考慮滲透系數(shù)在方向上的差異,本模型將在前人研究基礎(chǔ)上,分別對kv和kh進行預(yù)測,以使預(yù)測結(jié)果更符合工程實際。

2 工程應(yīng)用

2.1 工程概況

場地位于連云港港旗臺作業(yè)區(qū)南部,已由30萬t航道二期疏浚土吹填成陸(工程范圍內(nèi)),土體地下深度約4 m內(nèi)均為吹填土?？讐红o力觸探試驗(CPTU)采用荷蘭Geomil靜力觸探試驗設(shè)備,試驗結(jié)果如圖1所示(以JT33靜探孔為例)。

圖1 Geomil靜力觸探試驗結(jié)果(錐尖、側(cè)壁、孔壓)Fig.1 Geomil cone penetration test results (cone tip, side wall, pore pressure)

工程布置有靜力觸探試驗孔65個、取土試樣鉆孔15個,且每個鉆孔都與一個靜探孔相鄰。根據(jù)取土鉆孔的位置及取土土樣的深度匹配其對應(yīng)的靜探孔數(shù)據(jù),建立CPTU數(shù)據(jù)與室內(nèi)試驗結(jié)果的嚴格對應(yīng)關(guān)系。室內(nèi)土工試驗共設(shè)置了11組垂直向滲透系數(shù)(kv)和8組水平向滲透系數(shù)(kh)測定試驗,以“QT鉆孔號-土樣編號”的形式命名土樣,試驗結(jié)果如表2所示。

表2 土工試驗滲透系數(shù)數(shù)據(jù)集Tab.2 Permeability coefficient data set of geotechnical test

本工程鉆孔取得的土試樣等級為Ⅰ～Ⅱ級,即擾動程度為不擾動～輕微擾動。土樣的采取、蠟封、儲存、運輸及室內(nèi)試驗嚴格按照有關(guān)規(guī)定執(zhí)行,因此認為室內(nèi)土工試驗結(jié)果較為可靠。

2.2 XGBoost模型的訓(xùn)練與預(yù)測

2.2.1 輸入變量和標簽值的選取

由于靜探孔的探測精度為0.1 m,即隨深度變化每0.1 m可得一組CPTU數(shù)據(jù),為將其與滲透系數(shù)測定值對應(yīng),假定土樣的取土深度范圍內(nèi)每0.1 m的土層滲透系數(shù)都等于該段土樣的滲透系數(shù)值。對CPTU數(shù)據(jù)的異常值進行剔除后,最終總數(shù)據(jù)集共有47組CPTU-kv數(shù)據(jù)和41組CPTU-kh數(shù)據(jù)。模型的輸入變量從CPTU測試指標中選取,根據(jù)經(jīng)驗,考慮使用靜探數(shù)據(jù)轉(zhuǎn)換后得到的歸一化錐尖阻力Qtn和歸一化摩阻比Fr的組合,或直接使用土類指數(shù)Ic作為輸入變量對標簽值進行預(yù)測,相關(guān)計算公式[15]為:

Ic=[(3.47-lgQtn)2+(lgFr+1.22)2]0.5

(5)

Qtn=((qt-σv0)/pa)×(pa/(σ′v0))n

(6)

Fr=fs/(qt-σv0)×100%

(7)

式中:n為應(yīng)力指數(shù),n=0.381×Ic(RW)+0.05×(σ′v0/pa)-0.15;σv0、σ′v0分別為總上覆土壓力和上覆土有效應(yīng)力;qt為孔壓修正錐尖阻力,qt=qc+(1-a)u2,a為探頭面積比,取0.8;pa為參考壓力,pa=100 kPa=0.1 MPa。

對于標簽值的選取,考慮到表2中滲透系數(shù)的數(shù)量級均在10-2以下,故對滲透系數(shù)取對數(shù),將其作為標簽值,并以原數(shù)值設(shè)置對照組進行對比。

本文算法均通過Python的XGBoost模塊來實現(xiàn)。采用默認超參數(shù)值建立預(yù)分析模型,對兩類輸入變量和兩類標簽值(以kv為例)形成的4種預(yù)測組合分別進行訓(xùn)練預(yù)測。以R2為評價指標,取模型運行5次結(jié)果的均值,如表3所示。

表3 預(yù)分析結(jié)果表Tab.3 Results of pre-analysis

從表3中可見,以Qtn+Fr作為輸入變量、lg(kv)作為標簽值時,R2最高?？梢奞tn+Fr比Ic單值能更全面地反映土體的性質(zhì),更適合模型學(xué)習(xí);滲透系數(shù)原數(shù)值過小且接近,將其作為標簽值不利于模型分辨,在取對數(shù)后模型預(yù)測效果有較大提升。所以最后選取Qtn+Fr作為輸入變量、滲透系數(shù)的對數(shù)值作為標簽值建立預(yù)測模型。

2.2.2 模型的訓(xùn)練

1)kv預(yù)測模型

本模型的調(diào)參策略:先通過學(xué)習(xí)曲線調(diào)整弱評估器數(shù)量,再通過網(wǎng)格搜索調(diào)整其余超參數(shù)。將總數(shù)據(jù)集按7∶3的比例隨機劃分為訓(xùn)練集和測試集,為保證嚴謹性,只對訓(xùn)練集而非全集使用K折交叉驗證來估計模型的泛化誤差?？紤]到樣本數(shù)據(jù)量小,取K=5。模型的學(xué)習(xí)曲線圖(以R2為評價指標)如圖2所示。

圖2 訓(xùn)練集學(xué)習(xí)曲線圖Fig.2 Learning curve of training set

圖2中,R2、方差、總泛化誤差分別在弱評估器數(shù)量為29、14、19時達到最佳取值,故最后取num_round=19。經(jīng)過網(wǎng)格搜索對其余參數(shù)進行調(diào)整后,確定了一組最佳的超參數(shù)組合如表4所示。

表4 kv預(yù)測模型的最佳超參數(shù)Tab.4 Optimal hyperparameters of kv prediction model

設(shè)置該組超參數(shù),采用經(jīng)訓(xùn)練集訓(xùn)練的模型作為預(yù)測模型對測試集進行預(yù)測(見圖3)。圖3中,該模型對測試集的預(yù)測值與實測值的決定系數(shù)R2達到0.932,意味著模型的輸入變量即歸一化錐尖阻力Qtn和歸一化摩阻比Fr對滲透系數(shù)的對數(shù)標簽值有較強的解釋性。注意到模型對k=10-3cm/s左右的砂土和k=10-7cm/s左右的軟土均有不錯的預(yù)測表現(xiàn),說明模型對砂土和軟土都能適用,泛化性較強。

圖3 kv預(yù)測模型結(jié)果可視化對比Fig.3 Visual comparison of results of kv prediction model

2)kh預(yù)測模型

與kv預(yù)測模型類似,經(jīng)過學(xué)習(xí)曲線和網(wǎng)格調(diào)參后,確定了一組適用于kh預(yù)測的最佳超參數(shù),如表5所示。

表5 kh預(yù)測模型的最佳超參數(shù)Tab.5 Optimal hyperparameters of kh prediction model

該模型在測試集上的表現(xiàn)見圖4。圖4中,因砂土水平向滲透系數(shù)數(shù)據(jù)較少,故模型對砂土的滲透系數(shù)預(yù)測效果略有下降,對軟土預(yù)測效果仍較好。總體實測值與預(yù)測值的R2為0.891,仍為較高水平。

2.2.3 模型預(yù)測效果

將上述kv和kh預(yù)測模型應(yīng)用于項目CPTU數(shù)據(jù)集,將土工試驗結(jié)果的實測值與模型預(yù)測值直接進行對比,如圖5所示。

圖5中,考慮到工程實踐中滲透系數(shù)的計算可接受一個數(shù)量級范圍內(nèi)的誤差[16],故加入距實測值一個數(shù)量級的兩誤差線,該范圍內(nèi)的預(yù)測值都認為達到精度要求。為每個取土孔及土樣編號指定特定符號,將該土樣深度范圍內(nèi)每0.1 m對應(yīng)的滲透系數(shù)預(yù)測值顯示在對比圖中,該深度范圍內(nèi)滲透系數(shù)的平均值以較大的符號顯示。從圖5可見,除了-lg(k)=3附近的兩個數(shù)據(jù)在誤差范圍外,其余數(shù)據(jù)的預(yù)測值均與實測值較為接近。兩個不準數(shù)據(jù)來源于QT06-101、QT14-02土樣,均為砂土,這可能與土工試驗結(jié)果中砂土數(shù)據(jù)較少有關(guān),但該土樣深度范圍內(nèi)滲透系數(shù)預(yù)測值的平均值與實測值相當(dāng)接近。在測試集上,kv和kh的實測值與預(yù)測值取對數(shù)后的R2分別為0.932和0.891,均達較高水平。

3 各預(yù)測方法結(jié)果對比

3.1 Robertson經(jīng)驗公式

Robertson[3]提出了一種土壤滲透系數(shù)k與土類指數(shù)Ic的平均關(guān)系式,實測值在該經(jīng)驗公式圖中的分布情況如圖6所示。

圖6 實測值k-Ic分布圖Fig.6 k-Ic distribution of measured values

圖6中,由于數(shù)據(jù)量較小,實測值與公式曲線離散度較大。將Robertson公式的計算值分別與實測值的垂直向和水平向滲透系數(shù)(取對數(shù)后)進行比較,垂直向的R2為0.330,水平向的R2為-0.022(R2為負意味著預(yù)測值的預(yù)測效果還不如實測值的均值),總體R2為0.282。注意到垂直向的R2遠高于水平向,說明Robertson公式可能更適合預(yù)測kv而非kh。

3.2 (半)球面流模型經(jīng)驗公式

Elsworth等[4-5]、Chai等[6]、李鏡培等[7]分別結(jié)合球面流、半球面流、圓柱面徑向滲流模型推導(dǎo)出了CPTU測試指標與水平向滲透系數(shù)之間的理論關(guān)系式,其適用于正?；蜉p微超固結(jié)的黏性土和松散的無黏性土。將CPTU測試指標和探頭參數(shù)代入上述各經(jīng)驗公式,將計算值與土工試驗結(jié)果的水平向滲透系數(shù)實測值進行對比,如圖7所示。

圖7 預(yù)測值與實測值對比圖Fig.7 Comparison of predicted and measured values

圖7中,三個經(jīng)驗公式的計算值與實測值取對數(shù)后的R2分別為-3.290、0.571、0.505,可見,經(jīng)典Chai公式表現(xiàn)最好。

3.3 BP神經(jīng)網(wǎng)絡(luò)模型

BP神經(jīng)網(wǎng)絡(luò)雖然性能上比XGBoost和隨機森林等樹模型稍遜,但優(yōu)勢在于結(jié)構(gòu)簡單,能夠提取出較為簡潔的公式以供使用。訓(xùn)練過程不再贅述。

將kv和kh預(yù)測模型的權(quán)值和閾值參數(shù)在三層BP神經(jīng)網(wǎng)絡(luò)中反算出來,代入提取公式可得:

(8)

lg(kh)=-0.79·tansig(-0.29·Fr+

2.67·Qtn+1.99)-1.66·tansig(2.47·Fr-

1.62·Qtn+0.98)-0.08·tansig(-0.65·Fr+

2.16·Qtn-2.67)+0.15

(9)

將CPTU數(shù)據(jù)代入上述公式計算,并將結(jié)果反歸一化后可得滲透系數(shù)的BPNN預(yù)測值。 將實測值與預(yù)測值直接進行對比,如圖8所示。

圖8 BPNN模型預(yù)測值與實測值對比圖Fig.8 Comparison of predicted and measured values of BPNN model

總體來說,BPNN預(yù)測模型對kv預(yù)測效果較好,R2為0.824,對kh預(yù)測效果較差,R2為0.612。

3.4 各預(yù)測方法結(jié)果對比

將基于CPTU數(shù)據(jù)的各預(yù)測方法所得預(yù)測值與土工試驗實測值進行對比,如表6所示。

表6 各預(yù)測方法預(yù)測值與實測值對比表Tab.6 Comparison of predicted and measured values by each prediction method 單位: cm/s

表6中,取土孔所取土樣的淤泥、粉質(zhì)黏土、細砂、粉砂作為均質(zhì)沉積物,滲透系數(shù)實測值都基本符合kh/kv=1～1.5。XGB模型的預(yù)測值也符合該相關(guān)關(guān)系,但BPNN模型的kh預(yù)測值普遍低于kv預(yù)測值,所以選擇XGB模型進行滲透系數(shù)預(yù)測更合理。

計算各預(yù)測方法的計算值與實測值(取對數(shù)后)的距離,使用R2、均方誤差MSE、均方根誤差RMSE、平均絕對誤差MAE、平均絕對百分比誤差MAPE 5個常用評估指標進行評價,以R2為主要評價標準降序排序,如表7所示。其中,對于機器學(xué)習(xí)預(yù)測模型,則根據(jù)測試集而非全集的表現(xiàn)進行評估。

表7 各預(yù)測方法計算值與實測值距離的評估指標Tab.7 Evaluation indexes of distance between calculated and measured values by each prediction method

從表7中可見,基于XGB的kv預(yù)測模型表現(xiàn)最好,其次是基于XGB的kh預(yù)測模型?；贐P神經(jīng)網(wǎng)絡(luò)的預(yù)測模型表現(xiàn)雖然不如XGB模型,但也優(yōu)于其它經(jīng)驗公式。前人的經(jīng)驗公式中,Chai公式有著最好的預(yù)測表現(xiàn)。由于本文所用數(shù)據(jù)集的數(shù)據(jù)量有限,在其他地區(qū)的土壤滲透系數(shù)預(yù)測中,XGB預(yù)測模型不一定能有如此表現(xiàn),但對于連云港地區(qū)類似吹填土區(qū)域的滲透特性研究,還是有一定的參考價值。

4 結(jié) 論

1) 本文根據(jù)土體的CPTU數(shù)據(jù)與其滲透系數(shù)之間的關(guān)聯(lián)性,提出了一種基于XGBoost機器學(xué)習(xí)算法的滲透系數(shù)反演模型,該模型適用于廣泛土類,且考慮了水平向和垂直向滲透系數(shù)的差異。

2) 通過對輸入變量和標簽值的合理選取以及對最優(yōu)超參數(shù)組合的確定,模型最終達到了較高的預(yù)測精度。該XGBoost預(yù)測模型的kv預(yù)測值與實測值的R2可達0.932,kh預(yù)測值與實測值的R2可達0.891,遠高于傳統(tǒng)的BPNN模型和前人經(jīng)驗公式。

3) 本文的工程實例表明,基于機器學(xué)習(xí)算法構(gòu)建的XGBoost模型合理可行。XGBoost模型學(xué)習(xí)速度快、預(yù)測精度高、需調(diào)參數(shù)少,適合復(fù)雜原理相關(guān)的土體參數(shù)反演,對于滲透特性研究具有一定的參考價值。