基于離散型二項式系數組合模型的黃土濕陷性評估

任文博， 劉云龍*， 李佳佳， 李碩磊, 李文慶

(1.鄭州大學土木工程學院， 鄭州 450001； 2.中國建筑第七工程局有限公司， 鄭州 450004)

黃土是一種在干旱和半干旱環境中形成的特殊沉積物，在中國分布廣泛,約為6.3×105km2，其中60%為濕陷性黃土[1-2]。濕陷性黃土一般天然含水量較低且結構松散擁有大孔隙，當土體處于低含水量狀態下強度高壓縮性低。然而一旦受到雨水的滲入，土體結構被破壞，強度降低，壓縮性提高，在土體自重或上覆荷載作用下，土體被壓縮，產生濕陷變形[3]。由于以上特性，濕陷性黃土地區的建筑在受到生活用水，管道泄漏和降雨等狀況時極易出現地基沉降，結構開裂和建筑物傾斜等工程問題[4]，因此對黃土濕陷性進行評估至關重要。

黃土的濕陷現象是一個多因素耦合作用的物理化學過程。目前中外解釋黃土濕陷的假說眾多，如可溶性鹽假說、毛細血管假說、欠壓密假說和微結構不平衡吸力假說等，但均不能很好的解釋黃土的濕陷性，因此一般采用試驗的方法對黃土濕陷性進行評估[5-6]。黃土的濕陷性一般用濕陷性系數δs來評估，濕陷性系數為在規定壓力下浸水黃土單位厚度土層的濕陷量，當δs≥0.015時為濕陷性黃土。此外，根據濕陷性系數的大小可將濕陷性黃土分為輕微濕陷性黃土(0.015<δs≤0.030)、中等濕陷性黃土(0.030<δs≤0.070)和強烈濕陷性黃土(δs>0.070)[7]。目前，黃土濕陷性系數的測定方法分為兩種，一種為現場試驗，另一種為室內試驗。現場試驗結果較準確，但周期長投資大，故一般采用室內試驗進行濕陷性系數測定。室內試驗方法又分為雙線法和單線法，雖室內試驗相對于現場試驗較易實現，但由于測定工作量龐大，步驟煩瑣，仍需投入大量的人力、物力和時間且測量精度有限[3]。相對于濕陷性系數，土體的常規物性指標，如孔隙比、干密度、含水量、塑限、液限、塑性指數和液性指數等，易于測量、誤差小且能批量測量。因此一些學者嘗試通過建立土體常規物性指標與黃土濕陷性系數的關系，快速批量的對特定地區黃土濕陷性系數進行預測，取得了豐碩的成果。蘇強[8]、舒志樂等[9]通過對黃土濕陷性系數和物性指標進行相關性分析，分別建立了黃土濕陷性系數與物性指標的一元和多元線性回歸式；李萍等[6]采用多因素分析法建立了濕陷性系數與物性指標間的多元非線性回歸式；邵生俊等[10]、徐志軍等[11]采用因子分析法消除物性指標之間的共線性，建立了相關回歸方程；井彥林等[12]采用最小二乘法支持向量機建立了黃土濕陷性系數預測模型；馮小東[13]采用反向傳播(back propagation，BP)神經網絡對黃土濕陷性進行預測；李瑞娥等[14]利用模糊信息優化技術建立了黃土濕陷性系數與物性指標間的模糊關系；馬閆等[15]建立了黃土濕陷性系數徑向基函數(radial basis function，RBF)神經網絡預測模型；韓曉萌[16]采用自適應模糊神經網絡利用物性指標對黃土濕陷性進行預測。

綜上，雖然不少學者利用黃土基本物性指標建立了不同的黃土濕陷性系數單一預測模型，但單一模型都基于一定的假定，存在自身缺陷，所以精度難以提高。為此，基于4種不同單一預測模型，提出了離散型二項式系數組合模型，使單一模型間能夠充分互補，預測精度大幅度提高。

1 預測模型原理

1.1 單一預測模型

1.1.1 多元線性回歸

在科學研究中，當因變量y受到多個自變量xj(j=1,2,…，m)的影響時,可以采用多元線性回歸對y進行預測[17]。多元線性回歸基本模型為

(1)

(2)

令

則式(2)可寫為

(3)

式(1)表示m+1維空間中的平面，當m=2時為真正的平面，如圖1所示。

β=[b0,b1,b2, …，bm]T=(XTX)-1XTY

(4)

ei為第i組數據的回歸誤差；x1、x2為自變量圖1 真實值與回歸值示意圖Fig.1 Schematic diagram of true value and regression value

1.1.2 BP神經網絡

人工神經網絡是基于擁有大量神經元結構的生物大腦建立的數據分析預測模型。其中由于BP神經網絡具有容錯率高、非線性映射和自學習適應等優點，是目前人工神經網絡應用最廣泛的一種[18]。BP神經網絡由一個輸入層、一個或多個隱含層和一個輸出層組成，屬于多層前饋神經網絡，且誤差反向傳遞即計算某一節點誤差時需先求得與之相連的下一節點的誤差值，因此誤差必須自輸出層開始計算。由于不同變量數據的數值大小存在一定差異性，為避免數值較大變量覆蓋較小變量，需先將數據進行歸一化處理[19]。

BP神經網絡基本模型如圖2所示，m個自變量xj(j=1,2,…,m)作為輸入節點，因變量y作為輸出節點。運用程序調整每個節點權值ω，進而得到輸出值，將輸出值與真實值進行對比分析，若精度未達到要求，需再次調整各節點權重ω，直至精度滿足要求。

由Kolmogrov定理可知，一個由m個輸入節點，s個輸出節點和2m+1個隱含節點的3層BP神經網絡可以精確的反映任何映射關系[19]，采用上述3層BP神經網絡開展預測。

p為第一隱含層編號為p的節點；q為第二隱含層編號為q的節點；ωmp為編號為m的自變量與第一隱含層編號為p的節點之間的權重；ωpq為第一隱含層編號為p的節點與第二隱含層編號為q的節點之間的權重；ωq1為第二隱含層編號為q的節點與輸出變量間的權重圖2 BP神經網絡流程圖Fig.2 BP neural network flow chart

1.1.3 支持向量機回歸

支持向量機(support vector machine,SVM)是基于結構風險最小化原則提出的統計學理論，其具有良好的泛化能力，可將問題轉化為凸優化問題，使局部最優解即為全局最優解。如圖3(a)所示，三角形和圓形分別代表了兩類樣本，支持向量機的原理為尋求一個既能保證兩類樣本無誤的分隔開又能使分類間隔最大的最優超平面。對于樣本非線性可分的情況，如圖3(b)所示，可以運用非線性變換將樣本映射到更高的維度中，再尋求最優超平面。

L為最優超平面；L1、L2分別為平行于L且經過各分類樣本中距L最近的樣本的平面圖3 支持向量機示意圖Fig.3 Support vector machine diagram

(4)

(5)

(6)

式中：υ為線性回歸的斜率；x為自變量；b為線性回歸的截距；min為求代數式的最小值。

(7)

(8)

引入Lagrange乘子αi和α′i，以及利用KKT(Karush-Kuhn-Tucker conditions)條件可得到回歸函數為[20]

(9)

現實中大多數據是非線性可分的，因此需要進行非線性支持向量機回歸。非線性支持向量機回歸即先通過非線性映射Φ將數據映射到高維特征空間后再進行線性回歸。通過核函數k(xi,xj)=Φ(xi)Φ(xj)來處理數據，不需要計算映射到高維空間數據坐標，避免了煩瑣的點積運算。目前關于核函數的選取沒有統一標準，一般認為高斯徑向基核函數具有較好的穩定性，故選用高斯徑向基核函數。非線性支持向量機回歸的表達式為

(10)

1.1.4 隨機森林回歸

隨機森林(random forest，RF)是一種基于決策樹(分類樹和回歸樹)的集成學習算法。由于其不易產生過擬合，泛化能力較強，對異常數據具有較好的魯棒性，常用于處理分類和回歸問題。隨機森林回歸原理如圖4所示，運用bootstrap放回隨機抽樣法從原始訓練集中抽取與原始訓練集數據數量相等的訓練子集，共抽取k組，用于訓練k棵回歸樹模型。將預測集數據分別帶入訓練好的回歸樹模型得到k組回歸值，再將k組預測集回歸值取均值即為隨機森林回歸值[21]。

圖4 隨機森林回歸原理Fig.4 Random forest regression principle

1.2 組合預測模型

組合預測模型就是將多種單一模型的預測結果，通過一定規則進行組合，以達到合理利用各單一模型預測結果，克服單一模型缺陷與誤差的目的。若共有T個單一預測模型，n組訓練樣本，H組預測樣本，則組合模型預測結果可表示為

(11)

1.2.1 方差倒數法組合模型

方差倒數法組合模型原理為計算各單一預測模型訓練集預測值與真實值的方差倒數，然后按方差倒數賦予各單一模型權重，方差越小賦予權重越高，權重的計算公式為

(12)

式(12)中：Dt為第t個單一預測模型訓練集方差。

1.2.2 二項式系數法組合模型

二項式系數法組合模型是基于預測集預測結果賦予各模型權重。先求出各單一預測模型的H個預測樣本的預測平均值,然后按均值從小到大對模型增序排列，再利用二項式系數賦予權重，使均值排序越靠近中間位置的模型權重越大，權重的計算公式為[22]

(13)

1.2.3 離散型二次項系數法組合模型

二項式系數法組合模型雖能使預測集的預測值靠近整體預測均值排序居中的單一預測模型的預測值，但預測集各模型整體預測均值排序與單個預測樣本預測值排序并不完全相同，甚至相反。這就可能使某一預測樣本的組合預測值更靠近最大值或最小值，而不是排序居中的預測值。因此，針對這一問題提出了離散型二項式系數法組合模型，該方法即對每一個預測樣本的各單一模型預測值進行增序排列，并用二項式系數賦予權重，以保證每一個預測樣本的組合預測值都能靠近該預測樣本單一模型預測值增序排列的中間值，第h個預測樣本的預測值增序排列第t位的單一預測模型權重計算公式為

(14)

式(14)中：ωt,h為第h個預測樣本的預測值增序排列第t位的單一預測模型權重，h=1,2,…，H。

2 建立數據庫及預測模型的算法實現

2.1 數據庫的創建

利用在陜西省“引漢濟渭”工程中收集到的60組濕陷性黃土相關數據[23]，建立了濕陷性黃土物性指標數據庫，包括濕陷性系數、干密度、孔隙比、天然含水量、飽和度、塑限、液限、塑性指數和液性指數等。數據庫前50組作為訓練集，通過訓練得到該地區利用黃土物性指標預測濕陷性系數的預測模型。后10組數據作為預測集，用來驗證預測模型的可行性與準確性。

2.2 離散型二項式系數組合模型的算法實現

黃土濕陷性系數預測模型利用MATLAB編程實現，如圖5所示。具體步驟如下。

圖5 預測模型MATLAB算法實現流程Fig.5 MATLAB algorithm implementation process of prediction model

步驟1導入數據庫前50組數據計算相關系數表，并結合隨機森林重要性指數排序綜合確定擬用物性指標。

δs為濕陷性系數；ρd為干密度；e為孔隙比；ω為天然含水量；Sr為飽和度；wL為液限；wP為塑限；IP為塑性指數；IL為液性指數圖6 各物性指標間相關系數及散點圖Fig.6 Correlation coefficients and scatter plots among various physical properties

步驟2導入已確定物性指標相關數據，運用多元線性回歸得到單一預測模型1，將導入數據進行歸一化處理，先后進行支持向量機回歸(SVR)、BP神經網絡預測和RF回歸，得到單一預測模型2、3和4。

步驟3將后10組數據導入各單一預測模型，得到4種單一模型預測結果，再分別運用方差倒數法、二項式系數法和離散型二項式系數法計算模型權重，得到3種組合模型預測結果。

3 黃土物性指標的選取

3.1 相關系數法

相關系數R表示兩變量的相關程度，計算如式(15)所示， 且|R|≤1，|R|越接近1，相關性越大[17]。

(15)

為了選擇合適的物性指標，首先計算各物性指標間的相關系數進行相關性分析，如圖6所示。上三角為各物性指標間的相關系數，下三角為對應的散點圖，對角線為各物性指標的頻率分布直方圖。易見濕陷性系數δs與各物性指標相關性排序為:飽和度Sr>干密度ρd>孔隙比e>液性指數IL>天然含水量ω>塑性指數IP>液限wL>塑限wP，又由于干密度與孔隙比相關性接近1(0.986)，故兩者取其中之一(干密度)分析即可。而塑性指數、液限和塑限與濕陷性系數的相關系數均小于0.15，相關性很差，因此不參與預測。按照相關系數法選取飽和度、干密度、液性指數和天然含水量作為預測物性指標。

3.2 隨機森林重要性指數法

隨機森林不僅能夠用于分類與回歸，還可進行變量重要性分析，得到各因變量的重要性指數。如前文所述，隨機森林回歸在采用bootstrap抽樣方法時，若原始訓練集有n個樣本，則在某一子訓練集的抽取中，每個樣本未被抽中的概率均為(1-1/n)n，當n為50時，概率約為36.4%，這部分未被抽中的樣本稱為袋外(out of bag, OOB)數據。將每棵樹的袋外數據帶入與之相對應的回歸樹可得到袋外誤差errOOB1。然后對袋外數據的某一因變量xi施加噪聲干擾(即改變其數值)，得到干擾后袋外誤差errOOB2。依次計算出k棵回歸樹的errOOB1和errOOB2，代入隨機森林重要性指數(feature importance measures, FIM)計算公式[式(16)]，計算結果如圖7所示。

(16)

由圖7可見，重要性排序前五的物性指標與相關系數法相同，故按照隨機森林重要性指數法也選取飽和度、干密度、液性指數和天然含水量作為預測物性指標。

圖7 隨機森林重要性指數Fig.7 Random Forest importance exponent

3.3 選取物性指標數據

綜合相關系數與隨機森林重要性指數，選取飽和度Sr、干密度ρd(單位：g/cm3)、液性指數IL和天然含水量ω(單位：%)4種物性指標構建黃土濕陷性系數δs預測模型。訓練集(1～50)與測試集(51～60)相關數據如表1所示。

4 預測結果與精度分析

4.1 預測結果

將表1數據導入前文編寫的MATLAB程序中，可分別得到多元線性回歸(模型1)、BP神經網絡(模型2)、SVR(模型3)、RF回歸(模型4)、方差倒數法組合模型(模型5)、二項式系數法組合模型(模型6)和離散型二項式系數法組合模型(模型7)的預測集預測值如表2所示。上述4種單一模型和3種組合模型預測曲線如圖8所示。

表1 訓練集與預測集樣本數據

表2 各預測模型預測集濕陷性系數預測值

圖8 各模型預測曲線Fig.8 Each model prediction curve

4.2 模型預測精度分析

為了更清晰的分析圖8的預測結果，分別繪制了各預測模型相對誤差圖及殘差圖，如圖9、圖10所示。由圖9可知，4種單一模型(模型1～模型4)的最大相對誤差均超過15%，甚至將近40%，傳統組合模型(模型5和模型6)最大相對誤差均位于10%～15%，而所提出的組合預測模型(模型7)最大相對誤差小于10%，為8.98%。由圖10可知，四種單一模型(模型1～模型4)以及模型5的最大殘差均大于0.05，模型6和模型7的最大殘差均小于0.05，且模型7最大殘差最小。

圖9 模型預測相對誤差Fig.9 Relative error of model prediction

圖10 各模型預測殘差Fig.10 Model prediction residuals

此外，選用絕對誤差平方和(the sum of squares due to error，SSE)、標準誤差(standard error, SE)、相對誤差平方和(the sum of squares due to relative error, SSRE)、相對標準誤差(standard relative error, SRE)、平均絕對百分誤差(mean absolute percentage error, MAPE)和希爾不等系數(Theil inequality coefficient, Theil IC)精度判定指標開展模型精度分析。精度指標計算公式如式(17)所示，精度指標數值越小，模型精度越高，7種模型精度指標計算結果如表3所示。由表3可知，4種單一預測模型中模型1和3精度指標均優于模型2和4。對比單一模型與組合模型精度指標，可見組合模型精度指標整體優于單一模型精度指標，此外所提出的組合預測模型(模型7)所有精度指標相對于其他模型均為最優。

表3 模型預測精度指標對比

(17)

綜合圖9、圖10和表3的分析結果可知，3種組合預測模型預測精度整體高于4種單一預測模型。這是由于組合模型能夠對各單一模型進行模型互補，在某一樣本點即使某一單一預測模型有較大誤差，組合模型也可以利用其他單一模型給予彌補。本文提出的離散型二項式系數法組合預測模型(模型7)各指標均最小，預測精度最高，是因為該方法能夠先大致判斷在某一樣本處4種單一模型預測偏差大小，然后基于偏差大小賦予模型不同權重，偏離平均值越大賦予權重越小，以使組合預測結果偏差盡可能小，大大提高預測精度。

5 結論

黃土濕陷性系數是濕陷性黃土地區工程設計中的必要參數。通過創建黃土濕陷性系數及其物性指標數據庫，采用MATLAB數據分析軟件建立了4種單一預測模型、2種傳統組合預測模型和1種新型組合預測模型，并進行精度分析，得出如下主要結論。

(1)采用相關系數和隨機森林重要性指數綜合判定該地區物性指標對濕陷性系數的重要性，從大到小依次為飽和度、干密度、孔隙比、液性指數、天然含水量、塑性指數、液限和塑限，并選取飽和度、干密度、液性指數和天然含水量4種物性指標來開展黃土濕陷性系數預測。

(2)提出一套新型組合模型算法(離散型二項式系數法組合模型)，所提模型能夠使預測值盡可能地靠近每個樣本排序居中的單一模型預測值，以達到模型互補提高預測精度的目的。

(3)利用SSE、SE、SSRE、SRE、MAPE和Theil，IC6種精度指標及相對誤差和殘差圖對4種單一預測模型、2種傳統組合預測模型和所提出的新型組合預測模型進行精度分析，發現組合預測模型精度整體高于單一預測模型，且所提出的新型預測模型各精度指標均最小，精度最高，最大相對誤差為8.98%，平均相對誤差為3.43%。