三軸試驗抗剪強度參數值回歸分析法的區別與修正

余東明，姚海林，吳少鋒

（1. 中國科學院武漢巖土力學研究所 巖土力學與工程國家重點實驗室，武漢 430071；2. 廣東電網公司佛山南海供電局，廣東 佛山 528200）

1 引 言

巖土的三軸壓縮試驗是確定巖土抗剪強度參數黏聚力c和內摩擦角φ值的重要方法。試驗直接獲得軸向破壞強度σ1和對應的圍壓σ3，分別稱為大、小主應力，規范給出采用擬合出 σ1-σ3最佳曲線或基于強度準則擬合破壞包線，再推導出c、φ值的方法[1－3]。以往擬合曲線多采用手工作圖形式，存在較大隨意性。為了得到較為準確的 c、φ值，現在的巖土工作者傾向在試驗數據的基礎上采用較嚴謹的數學理論進行推導或利用軟件分析，以避免由于手工作圖產生的隨意誤差。劉善均等[4]依據破壞應力圓與強度包線盡可能相切，基于最小二乘法推導了求解 c、φ值的方程組，并編制計算機程序來求解方程組。Zambrano-Mendoza等[5]提出用應力圓圓心到強度包線距離進行最小二乘法回歸擬合時，應同時考慮破壞面正應力和剪應力的誤差。陳立宏等[6－7]推導了兩種使用線性回歸理論求解c、φ值的方法，指出了這兩種方法不同之處并初步分析了其原因。阮波等[8]根據非線性規劃理論，應用EXCEL軟件對三軸試驗的強度包線直接進行規劃求解獲取了抗剪強度參數值。本文在Mohr-Coulomb強度準則的假設下，深入最小二乘法的線性回歸原理，分析了兩種最常用的三軸試驗 c、φ值確定方法。揭示了這兩種方法結果差異的根本原因，提出了一種修正方法并從理論與實際兩方面對其進行了驗證。為合理確定巖土三軸試驗抗剪強度參數值提供依據。

2 經典最小二乘法的一元線性回歸

對樣本數據(xi，yi)，若存在yi=a+bxi+εi，并假定：

（1）誤差 εi是數學期望值為 0、方差不變的隨機變量，即 E(εi)=0，Var(εi)不變；

（2）xi為非隨機變量。

則使

滿足式（1）的a，b值是回歸方程y = a + bx的回歸系數的最佳線性無偏估計量，此時：

其中： Lxx、 Lyy與Lxy含義相同。

這就是經典最小二乘法的一元線性回歸計算方法[9－10]。這種方法的最小二乘估算結果最為理想，計算簡單、且易于理解，已成為線性回歸最常用的方法。

3 三軸試驗抗剪強度參數的兩種回歸方法

在進行巖土三軸試驗數據分析時，假設巖土的剪切破壞強度滿足Mohr-Coulomb準則：

利用三軸試驗測量數據，基于線性回歸的原理求解出抗剪強度參數c、φ值，常采用兩種方法。

3.1 σ1-σ3法

利用破壞應力圓由式（5）可推導出

其中，

巖土三軸試驗直接測得試樣的大小主應力σ1i和σ3i，所以可以依據線性回歸原理直接對(σ3i，σ1i)數據進行最小二乘法線性回歸，得到A、B值，再利用式（7）可以計算出c、φ值。這種方法在文獻[6]中被稱為“σ1-σ3法”。

3.2 p-q法

注意到式（6）可以變換為

則式（8）變為

則可以建立p-q坐標系統，將三軸試驗測量數據對（σ3i，σ1i）轉化為（pi，qi）數據，再依據線性回歸原理對得到的（pi，qi）數據進行最小二乘法線性回歸，可以得到C、D值，利用式（9）也可以計算出c、φ值。這里的p、q實際上就是σ3-σ1坐標系下破壞應力圓的圓心橫坐標和半徑。這種方法在文獻[6]中被稱為“p-q 法”。

3.3 “σ1-σ3法”和“p-q法”的區別

“σ1-σ3法”和“p-q法”都可以應用最小二乘法進行線性回歸求出回歸系數，再利用回歸系數與抗剪強度參數的關系式求解出c、φ。還可以看出“p-q法”的回歸方程形式是直接由“σ1-σ3法”的回歸方程形式推導出來的。從這個意義上說，“p-q法”和“σ1-σ3法”求出的抗剪強度參數值應該是一樣的。而文獻[6]無論從理論推導上還是從試驗結果的分析上都明確得出二者的抗剪強度參數值并不相同，文獻對此進行了分析，認為是σ1的誤差導致了二者的差別，并認為“σ1-σ3法”的結果更為準確，但沒有揭示導致二者差別的根本原因。

σ1和σ3是巖土三軸試驗的測量值。其中σ3是試件所受圍壓，是試驗前設定好的值，在試驗過程中σ3的水平是受到控制的，可忽略誤差認為是固定值。而σ1是在試驗過程中測量出的破壞強度值，是無法控制的，認為其誤差不可忽略，假定是數學期望值為0、且方差不變的隨機變量，即：

用“σ1-σ3法”進行線性回歸，則

可以看出，“σ1-σ3法”完全符合經典最小二乘法線性回歸的假設和形式，故“σ1-σ3法”所得到的A、B估值是A、B真值的最佳解，因而由此推算出的c、φ值也應是c、φ真值的最佳解。

由式（12）可得

即，

此時需要回歸的數據中不僅因變量的測量值qi具有隨機誤差ξi，而且自變量的測量值pi也具有隨機誤差ξi，即pi也是隨機變量。這違反了經典最小二乘法進行線性回歸的前提假設，即自變量觀測值為非隨機變量。這時基于經典的最小二乘法進行線性回歸已經不可用。

但“p-q法”仍然采用經典最小二乘法進行線性回歸，實際上是強行消除自變量的誤差，從式（15）可以看出，自變量的誤差與因變量的誤差是等量的，若因變量誤差不能忽略也就等價于自變量誤差同樣不能忽略。強行忽略自變量誤差的做法必然導致所得到的回歸系數C、D不是C、D真值的最佳解，而由此推算出的c、φ值也不是c、φ真值的最佳解，所以會與“σ1-σ3法”得到的c、φ值不同。

“p-q法”回歸方程的自變量含有 σ1，從而將 σ1的測量誤差引入了自變量，違反了經典最小二乘法線性回歸的基本假定，必然導致回歸結果出現偏差。這就是兩種方法的根本區別，也是文獻[6]中用兩種方法所得到的抗剪強度參數值不同的根本原因。而“σ1-σ3法”則完全滿足經典最小二乘法線性回歸的假定和形式。可以認為“σ1-σ3法”的回歸結果是最接近真實值的，得出的c、φ值也是最接近真實值的。所以對三軸試驗的σ1、σ3直接擬合得到抗剪強度參數值是最合理的。

4 對“p-q 法”的修正

由前述的分析可以看出，“p-q法”違反了經典最小二乘法線性回歸的基本假定，使結果出現偏差。若σ1的誤差可以忽略的話，無論“σ1-σ3法”還是“p-q法”得到的回歸系數都將是真值，兩種方法的回歸方程將完全等價，此時這兩種方法得到抗剪強度參數值完全一樣。但實際上σ1的誤差一般不能忽略，故從根本上說“p-q法”利用經典最小二乘法進行線性回歸是錯誤的。應只考慮用“σ1-σ3法”計算c、φ值。然而在巖土三軸試驗規范[1－2]中規定可利用p-q坐標系統分析抗剪強度參數值，故人們也習慣于“p-q 法”。

“p-q法”雖然符合三軸試驗結果分析的習慣，但得到的結果有偏差，并非最佳解，甚至根本上說是錯誤的，因此，必須要對其修正。結合前面的分析可以認為，“σ1-σ3法”得出的抗剪強度參數值是最合理的，并且是唯一的。而“p-q法”違反了經典最小二乘法線性回歸的基本假定，因此，若要修正“p-q法”就是要修正（pi，qi）線性回歸的方法使得從中得到的抗剪強度參數值與“σ1-σ3法”完全一致。

4.1 正交最小二乘法的基本原理

當xi，yi均為隨機變量，即：

對這種情況下的(xi，yi)進行一元線性回歸，可以采用Golub和Van Loan對經典最小二乘法的推廣，認為其基本依據應是使x方向及y方向的誤差總平方和最?。?/p>

即，

得到的a、b值，是同時考慮x和y兩個方向誤差的最佳估值。這里所推廣的最小二乘法被稱為整體最小二乘法[11]，也被稱之為正交最小二乘法[12]。

4.2 修正的“p-q法”的線性回歸系數

當用x代替σ3，y代替σ1，由式（6）、（10）得：

由式（15）可知，pi和qi都是隨機變量，首先違反了經典最小二乘法關于自變量為非隨機變量的前提假設；且pi、qi這兩個隨機變量的誤差項都是εi，也違反了正交最小二乘法關于自變量和因變量的誤差項要相互獨立的前提假設。故既不能采用經典最小二乘法也不能采用正交最小二乘法對（pi，qi）進行線性回歸。

但注意到此時情況與正交最小二乘法的回歸比較接近，即自變量和因變量都是隨機變量，故決定在正交最小二乘法的基礎上進行修正，提出：

當xi、yi均為隨機變量，即

前提條件為：εi、ξi是相關聯的，即二者非獨立；E(εi)= E(ξi)= 0；Var(εi)= Var(ξi)為不變值。

對這種情況下的(xi，yi)進行一元線性回歸，認為其基本依據應是使下式成立：

在式（22）中引入常數λ是用來表示因 xi、yi的誤差不獨立所帶來的影響，暫且稱λ為誤差非獨立影響系數。若能通過 xi、yi的誤差關系求解出常數λ的話，代入式（22）則可得出(xi，yi)的一元線性回歸系數。稱這種方法為修正的正交最小二乘法。

可以看出，（pi，qi）剛好可以滿足上述提出的修正的正交最小二乘法的前提條件，由式（15）還可以得到 pi、qi這兩個隨機變量的非獨立關系是二者的誤差項完全相同。基于本文提出的修正的最小二乘法進行（pi，qi）的線性回歸的方法，本文稱之為修正的“p-q法”。

由式（22）可得

即，

若能根據式（24）得到的C、D值的話，則修正的“p-q法”將可用于巖土抗剪強度指標的分析。

由式（24）得：

其中， i =1, 2, …,n。

由式（15）得：

由式（25）、（26）并注意到式（28），得：

由式（27）、（28）、（30）得：

由式（31）可得：

由式（31）～（33）及（28）得：

可得

由式（30）、（31）可得

由式（15）知 pi、qi誤差項相等，結合式（36）～（38）可得

則

式（30）即為修正的“p-q法”的線性回歸系數，再利用式（9）就可得到c、φ的最佳估值。

4.3 修正的“p-q法”與“σ1-σ3法”等效

修正的“p-q法”得出的c、φ如果是真值的最佳估值的話，那么它們應與“σ1-σ3法”得出的c、φ估值完全一致，即無論采用何種方法，c、φ最佳估值只有一組。若修正的“p-q法”與“σ1-σ3法”等效，則由式（7）、（9）可得：

故若能證明“σ1-σ3法”的回歸系數A、B值與修正的“p-q 法”的回歸系數 C、D 值滿足式（41）、（42），則即證明了修正的“p-q法”和“σ1-σ3法”等效。

由式（39）可得

注意到式（29），由式（20）可得：

需要注意的是，這里x =σ1，y =σ1。

將式（44）代入（43）可得：

對式（19）利用經典的最小二乘法進行回歸，由式（2）可得

由式（45）、（46）可知，式（42）成立，同理可證式（41）也成立。故修正的“p-q法”與“σ1-σ3法”是等效的得證。修正的“p-q法”得出和“σ1-σ3法”同樣的c、φ值，它們是最合理的估值。

5 計算實例

四川省某擬建隧道圍巖主要為千枚巖，對其進行室內三軸強度試驗，取天然含水狀態下的一組三軸試驗數據，如表1所示。

分別采用“σ1-σ3法”和修正的“p-q法”確定抗剪強度參數c、φ值。

表1 天然狀態下千枚巖三軸強度試驗數據Table 1 Triaxial test data of phylliteunder natural state

（1）σ1-σ3法

由式（2）、（3）可得 B=4.61，A=49.06，則回歸方程為

將A、B值代入式（7），求解可得：c = 11.42 MPa，φ = 40.06o。

（2）修正的“p-q 法”

由式（40）可得D = 0.64，C = 8.74，則回歸方程為

將 C、D 值代入式（9），求解可得：c= 11.42 MPa，φ = 40.06o。

可以看出“σ1-σ3法”和修正的“p-q法”計算出的抗剪強度參數值完全一致。

6 結 論

（1）從本質上揭示了兩種線性回歸方法得出的三軸試驗抗剪強度參數值不一致的原因。指出“p-q法”由于自變量為隨機變量，違背了經典最小二乘法線性回歸的基本假定，因而得到的結果不可靠；而“σ1-σ3法”則完全滿足經典最小二乘線性回歸法的假定和形式，所得出的結果是真值的最佳估值。

（2）參考正交最小二乘法的原理，對“p-q法”提出修正的最小二乘法形式。使“p-q法”的自變量和因變量同為隨機變量、且誤差項不獨立的前提條件在修正的最小二乘法中得到滿足。稱之為修正的“p-q法”并推導了修正后的回歸系數。從理論上證明了修正的“p-q法”與“σ1-σ3法”得到的抗剪強度參數值完全一致，給出的一個實際算例也驗證了這兩種方法得到的c、φ值相同。

[1]南京水利科學研究院. SL237-1999 土工試驗規程[S].北京: 中國水利水電出版社, 1999.

[2]交通部公路科學研究院. JTG E40-2007 公路土工試驗規程[S]. 北京: 人民交通出版社, 2007.

[3]中國水電顧問集團成都勘測設計研究院. DL/T 5368-2007 水電水利工程巖石試驗規程[S]. 北京: 中國電力出版社, 2007.

[4]劉善均, 王韋, 許唯臨. 莫爾圓求解 c、φ值的最佳擬合[J]. 四川大學學報（工程科學版）, 2002, 34(6): 40－42.LIU Shan-jun, Wang Wei, XU Wei-lin. The optimal fitting method for computing c、φ using circle of Mohr[J].Journal of Sichuan University (Engineering Science Edition), 2002, 34(6): 40－42.

[5]ZAMBRANO-MENDOZA O, VALKO P P, RUSSELL J E. Error-in-variables for rock failure envelope[J].International Journal of Rock Mechanics and Mining Sciences, 2003, 40(1): 137－143.

[6]陳立宏, 陳祖煜, 李廣信. 三軸試驗抗剪強度指標線性回歸方法的討論[J]. 巖土力學, 2005, 26(11): 1785－1789.CHEN Li-hong, CHEN Zu-yu, LI Guang-xin. Discussion of linear regression method to estimate shear strength parameters from results of triaxial tests[J]. Rock and Soil Mechanics, 2005, 26(11): 1785－1789.

[7]陳立宏, 陳祖煜, 李廣信. 線性回歸抗剪強度指標方法的改進[J]. 巖土力學, 2007, 28(7): 1421－1426.CHEN Li-hong, CHEN Zu-yu, LI Guang-xin. A modified linear regression method to estimate shear strength parameters[J]. Rock and Soil Mechanics, 2007, 28(7):1421－1426.

[8]阮波, 張向京, 彭意. Excel規劃求解三軸試驗抗剪強度指標 [J]. 鐵道科學與工程學報, 2009, 6(5): 57－60.RUAN Bo, ZHANG Xiang-jing, PENG Yi. Programming solver tools of Excel evaluate shear strength parameters from results of triaxial tests[J]. Journal of Railway Science and Engineering, 2009, 6(5): 57－60.

[9]張小蒂. 應用回歸分析[M]. 杭州: 浙江大學出版社,1991: 78－82.

[10]G. A. F. 塞伯. 線性回歸分析[M]. 北京: 科學出版社,1987: 125－128.

[11]黃開斌, 俞錦成. 整體最小二乘問題的解集與極小范數解[J]. 南京師大學報（自然科學版）, 1997, 20(4): 1－5.HUANG Kai-bin, YU Jin-cheng. Solution sets and minimum norn solution for the basic total least squares problem[J]. Journal of Nanjing Normal University(Natural Science), 1997, 20(4): 1－5.

[12]李玉武, 劉咸德. 正交最小二乘法及其在分析化學中的應用[J]. 礦巖測試, 2001, 20(4): 257－262.LI Yu-wu, LIU Xian-de. Orthogonal least squares method and its application in analytical chemistry[J]. Rock and Mineral Analysis, 20(4): 257－262.