余東明，姚海林，吳少鋒

（1. 中國(guó)科學(xué)院武漢巖土力學(xué)研究所 巖土力學(xué)與工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，武漢 430071；2. 廣東電網(wǎng)公司佛山南海供電局，廣東 佛山 528200）

1 引 言

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

2 經(jīng)典最小二乘法的一元線性回歸

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

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

（2）xi為非隨機(jī)變量。

則使

滿足式（1）的a，b值是回歸方程y = a + bx的回歸系數(shù)的最佳線性無(wú)偏估計(jì)量，此時(shí)：

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

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

3 三軸試驗(yàn)抗剪強(qiáng)度參數(shù)的兩種回歸方法

在進(jìn)行巖土三軸試驗(yàn)數(shù)據(jù)分析時(shí)，假設(shè)巖土的剪切破壞強(qiáng)度滿足Mohr-Coulomb準(zhǔn)則：

利用三軸試驗(yàn)測(cè)量數(shù)據(jù)，基于線性回歸的原理求解出抗剪強(qiáng)度參數(shù)c、φ值，常采用兩種方法。

3.1 σ1-σ3法

利用破壞應(yīng)力圓由式（5）可推導(dǎo)出

其中，

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

3.2 p-q法

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

則式（8）變?yōu)?/p>

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

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

“σ1-σ3法”和“p-q法”都可以應(yīng)用最小二乘法進(jìn)行線性回歸求出回歸系數(shù)，再利用回歸系數(shù)與抗剪強(qiáng)度參數(shù)的關(guān)系式求解出c、φ。還可以看出“p-q法”的回歸方程形式是直接由“σ1-σ3法”的回歸方程形式推導(dǎo)出來(lái)的。從這個(gè)意義上說(shuō)，“p-q法”和“σ1-σ3法”求出的抗剪強(qiáng)度參數(shù)值應(yīng)該是一樣的。而文獻(xiàn)[6]無(wú)論從理論推導(dǎo)上還是從試驗(yàn)結(jié)果的分析上都明確得出二者的抗剪強(qiáng)度參數(shù)值并不相同，文獻(xiàn)對(duì)此進(jìn)行了分析，認(rèn)為是σ1的誤差導(dǎo)致了二者的差別，并認(rèn)為“σ1-σ3法”的結(jié)果更為準(zhǔn)確，但沒(méi)有揭示導(dǎo)致二者差別的根本原因。

σ1和σ3是巖土三軸試驗(yàn)的測(cè)量值。其中σ3是試件所受圍壓，是試驗(yàn)前設(shè)定好的值，在試驗(yàn)過(guò)程中σ3的水平是受到控制的，可忽略誤差認(rèn)為是固定值。而σ1是在試驗(yàn)過(guò)程中測(cè)量出的破壞強(qiáng)度值，是無(wú)法控制的，認(rèn)為其誤差不可忽略，假定是數(shù)學(xué)期望值為0、且方差不變的隨機(jī)變量，即：

用“σ1-σ3法”進(jìn)行線性回歸，則

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

由式（12）可得

即，

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

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

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

4 對(duì)“p-q 法”的修正

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

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

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

當(dāng)xi，yi均為隨機(jī)變量，即：

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

即，

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

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

當(dāng)用x代替σ3，y代替σ1，由式（6）、（10）得：

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

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

當(dāng)xi、yi均為隨機(jī)變量，即

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

對(duì)這種情況下的(xi，yi)進(jìn)行一元線性回歸，認(rèn)為其基本依據(jù)應(yīng)是使下式成立：

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

可以看出，（pi，qi）剛好可以滿足上述提出的修正的正交最小二乘法的前提條件，由式（15）還可以得到 pi、qi這兩個(gè)隨機(jī)變量的非獨(dú)立關(guān)系是二者的誤差項(xiàng)完全相同?；诒疚奶岢龅男拚淖钚《朔ㄟM(jìn)行（pi，qi）的線性回歸的方法，本文稱之為修正的“p-q法”。

由式（22）可得

即，

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

由式（24）得：

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

由式（15）得：

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

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

由式（31）可得：

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

可得

由式（30）、（31）可得

由式（15）知 pi、qi誤差項(xiàng)相等，結(jié)合式（36）～（38）可得

則

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

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

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

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

由式（39）可得

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

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

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

對(duì)式（19）利用經(jīng)典的最小二乘法進(jìn)行回歸，由式（2）可得

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

5 計(jì)算實(shí)例

四川省某擬建隧道圍巖主要為千枚巖，對(duì)其進(jìn)行室內(nèi)三軸強(qiáng)度試驗(yàn)，取天然含水狀態(tài)下的一組三軸試驗(yàn)數(shù)據(jù)，如表1所示。

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

表1 天然狀態(tài)下千枚巖三軸強(qiáng)度試驗(yàn)數(shù)據(jù)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法”計(jì)算出的抗剪強(qiáng)度參數(shù)值完全一致。

6 結(jié) 論

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

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

[1]南京水利科學(xué)研究院. SL237-1999 土工試驗(yàn)規(guī)程[S].北京: 中國(guó)水利水電出版社, 1999.

[2]交通部公路科學(xué)研究院. JTG E40-2007 公路土工試驗(yàn)規(guī)程[S]. 北京: 人民交通出版社, 2007.

[3]中國(guó)水電顧問(wèn)集團(tuán)成都勘測(cè)設(shè)計(jì)研究院. DL/T 5368-2007 水電水利工程巖石試驗(yàn)規(guī)程[S]. 北京: 中國(guó)電力出版社, 2007.

[4]劉善均, 王韋, 許唯臨. 莫爾圓求解 c、φ值的最佳擬合[J]. 四川大學(xué)學(xué)報(bào)（工程科學(xué)版）, 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]陳立宏, 陳祖煜, 李廣信. 三軸試驗(yàn)抗剪強(qiáng)度指標(biāo)線性回歸方法的討論[J]. 巖土力學(xué), 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]陳立宏, 陳祖煜, 李廣信. 線性回歸抗剪強(qiáng)度指標(biāo)方法的改進(jìn)[J]. 巖土力學(xué), 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規(guī)劃求解三軸試驗(yàn)抗剪強(qiáng)度指標(biāo) [J]. 鐵道科學(xué)與工程學(xué)報(bào), 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]張小蒂. 應(yīng)用回歸分析[M]. 杭州: 浙江大學(xué)出版社,1991: 78－82.

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

[11]黃開斌, 俞錦成. 整體最小二乘問(wèn)題的解集與極小范數(shù)解[J]. 南京師大學(xué)報(bào)（自然科學(xué)版）, 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]李玉武, 劉咸德. 正交最小二乘法及其在分析化學(xué)中的應(yīng)用[J]. 礦巖測(cè)試, 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.