莫林利，趙秀紹，王 旭，鄭 偉

(1.華東交通大學(xué) 軟件學(xué)院，江西 南昌 330013;2.華東交通大學(xué) 土木建筑學(xué)院，江西 南昌 330013)

基于Excel的直剪試驗數(shù)據(jù)處理方法

莫林利1，趙秀紹2，王 旭2，鄭 偉2

(1.華東交通大學(xué) 軟件學(xué)院，江西 南昌 330013;2.華東交通大學(xué) 土木建筑學(xué)院，江西 南昌 330013)

回歸分析中的最小二乘法是處理直剪試驗數(shù)據(jù)誤差最小的方法，但計算工作量大且容易出錯。應(yīng)用Excel軟件給出可用于直剪試驗計算的5種方法，可顯著減少計算工作量并降低出錯率。實際工程數(shù)據(jù)計算對比表明，除基于點(diǎn)線距的規(guī)劃求解法外，其余4種方法均為基于擬合誤差的計算方法，兩種原理的計算結(jié)果差距可忽略不計。對于無黏性土情況，除INTERCEPT函數(shù)與SLOPE函數(shù)聯(lián)合求解法無法實現(xiàn)強(qiáng)制黏聚力為0外，其余方法均可實現(xiàn)。實際應(yīng)用表明Excel圖解法、TREND函數(shù)法及截距和斜率函數(shù)聯(lián)合求解法操作簡單，易于應(yīng)用。

直剪試驗 Excel內(nèi)部函數(shù) 最小二乘法 規(guī)劃求解

邊坡土的抗剪強(qiáng)度是鐵路路基邊坡穩(wěn)定性定量評價和計算滑坡推力的重要參數(shù)，有時內(nèi)摩擦角相差1°～2°，滑坡推力就可能成倍增加［1］。馬顯春等［2］采用正交分析法對滑坡穩(wěn)定性的影響因子進(jìn)行了分析，得出滑坡影響因素敏感性順序為內(nèi)摩擦角＞黏聚力＞地下水位＞重度＞坡高?？梢姾侠泶_定巖土材料的抗剪強(qiáng)度指標(biāo)對于工程設(shè)計、施工等十分重要［3-5］。

直剪試驗的數(shù)據(jù)處理實際上是求出土的抗剪強(qiáng)度參數(shù)，即黏聚力 c和內(nèi)摩擦角 φ。斜截式計算法是在試驗中測得N，M兩點(diǎn)的垂直壓應(yīng)力σN，σM，其對應(yīng)的抗剪強(qiáng)度為τN，τM。通過 N，M的直線方程可以獲得土的抗剪強(qiáng)度指標(biāo)［6］。此法的不足之處是非常明顯的，計算時僅能使用2組數(shù)據(jù)。《土工試驗技術(shù)手冊》［7］建議的方法為繪制抗剪強(qiáng)度與垂直壓力關(guān)系曲線，直線的傾角為內(nèi)摩擦角φ，縱坐標(biāo)的截距為黏聚力c。但是在繪庫侖破壞直線時，一般是繪圖者根據(jù)經(jīng)驗和直觀判斷，盡量讓落在直線兩邊的點(diǎn)數(shù)大致相等，然后求出截距和斜率。由于它是人為地根據(jù)經(jīng)驗和直觀判斷得出的，c，φ值往往帶有一定的人為性和不確定性。采用最小二乘法對試驗數(shù)據(jù)進(jìn)行擬合求解［8-11］，計算量大且容易出錯。本文探討如何巧用Excel的各項功能簡化計算過程，并考慮直剪試驗中c值為負(fù)或強(qiáng)制c=0(無黏性土)的情況，使數(shù)據(jù)處理達(dá)到全自動化。

1 最小二乘法求解原理

為了便于說明各種求解方法的原理與實現(xiàn)方法，以江西撫州八景鎮(zhèn)火車站路基工程一組土樣實測數(shù)據(jù)(見表1)為例進(jìn)行分析。

表1 撫州八景鎮(zhèn)路基土樣測試數(shù)據(jù)

直剪試驗求解中，當(dāng)繪制測點(diǎn)的破壞包線時，由于土樣不均勻，試驗讀數(shù)誤差等原因，往往繪制的直線不能通過所有的測點(diǎn)，總是存在一些誤差，根據(jù)這些誤差的平方和為最小，可以確定一條近似的破壞包線。設(shè)直剪試驗中每組n個土樣的垂直應(yīng)力、抗剪強(qiáng)度分別為σi，τi，利用庫侖公式 τ=σtanφ+c進(jìn)行擬合，其中c，φ被視為待定系數(shù)。以偏差+c)］2最小為條件來確定c，φ值。運(yùn)用高等數(shù)學(xué)的多元函數(shù)求極值的方法即可得到當(dāng) M為最小時的 c，φ值。分別對c，φ求一階偏導(dǎo)并令其等于0得［8］

根據(jù)式(1)，式(2)解得

如果用上述數(shù)學(xué)法求解，計算工作量大且易出錯。Excel軟件中有豐富的回歸分析函數(shù)和規(guī)劃求解方法，只要挖掘其功能，就可以簡單有效地完成上述計算。

2 直剪試驗的Excel求解方法

2.1 Excel圖解法

Excel軟件有強(qiáng)大的繪圖工具，先把試驗數(shù)據(jù)填入數(shù)據(jù)表格中，然后插入圖表，并為數(shù)據(jù)添加趨勢線和顯示公式，就可以得到直剪試驗圖。具體實現(xiàn)方法為:選擇正應(yīng)力和剪應(yīng)力數(shù)據(jù)，選擇“插入”→“圖表”→“散點(diǎn)圖”，然后選擇圖形中的點(diǎn)后“添加趨勢線”，在“選項”欄中選擇“顯示公式”，可直接得到擬合公式 y= 0.430 6x+3.260 9，進(jìn)而得 c=3.260 9 kPa，tanφ= 0.430 6，即φ=23.297 1°。

2.2 TREND函數(shù)求解法

2.2.1 TREND函數(shù)求解的基本原理

由于直剪試驗各測點(diǎn)存在線性關(guān)系，這樣就可以方便地利用 Excel中的線性插值函數(shù) TREND。TREND函數(shù)用途:返回一條線性回歸擬合線的一組縱坐標(biāo)值(y值)。即找到適合給定的數(shù)組known_y's和known_x's的直線，并返回指定數(shù)組new_x's值在直線上對應(yīng)的y值。TREND函數(shù)的語法為TREND(known _y's，known_x's，new_x's，const)。其中known_y's為已知的τ=σtanφ+c=σk+c中的τi值集合，known_x's為已知關(guān)系τ=σk+c中可選的σi值的集合，new_x's為需要函數(shù)TREND返回對應(yīng)新σ值的新τ值，const參數(shù)為邏輯值指明是否強(qiáng)制常數(shù)項c為0。當(dāng)擬合發(fā)現(xiàn)c值為負(fù)時，由于黏聚力不可能為負(fù)值，因此可以強(qiáng)制c值為0。

TREND函數(shù)求解原理如圖1所示，很顯然，對于抗剪強(qiáng)度包線τ=σtanφ+c，當(dāng)new_x's=σ=0時，其返回的函數(shù)值為τ0=c，即為黏聚力值。當(dāng)new_x's= σ=1時，函數(shù)返回值 τ1=1×tanφ+c，根據(jù)圖1中的幾何關(guān)系可得tanφ=(τ1－τ0)/1=τ1－τ0。

2.2.2 TREND函數(shù)求解法的實現(xiàn)

將表1中的原始數(shù)據(jù)填入Excel表中，如圖2(屏幕截圖)中所示。

圖1 TREND函數(shù)求解原理

圖2 TREND函數(shù)求解抗剪強(qiáng)度指標(biāo)示例

在C2和D2單元格中，分別輸入“=TREND(B2: B5，A2:A5，0)”，“=ATAN(TREND(B2:B5，A2:A5，1)－TREND(B2:B5，A2:A5，0))*180/PI()”，即可求得該土樣的黏聚力c和內(nèi)摩擦角φ。

TREND(B2:B5，A2:A5，0)在σ=0時可得τ0=c; TREND(B2:B5，A2:A5，1)在 σ=1 kPa時可得 τ1=c +k;TREND(B2:B5，A2:A5，1)－TREND(B2:B5，A2: A5，0)為系數(shù)k。

當(dāng)土樣為無黏性土或計算出的黏聚力c值為負(fù)值時，可以強(qiáng)制TREND()函數(shù)中const項參數(shù)為0。對表1的測試數(shù)據(jù)求解結(jié)果為:c=3.260 9 kPa，φ= 23.297 1°。

2.3 INTERCEPT函數(shù)與SLOPE函數(shù)聯(lián)合求解法

2.3.1 INTERCEPT函數(shù)與SLOPE函數(shù)聯(lián)合求解原理

在Excel軟件中，函數(shù)INTERCEPT和SLOPE利用現(xiàn)有的x值和y值計算直線與y軸的截距和斜率。其中INTERCEPT(known_y's，known_x's)函數(shù)可以利用現(xiàn)有的σ值和τ值計算直線與τ軸的截距，即抗剪強(qiáng)度指標(biāo)黏聚力 c;SLOPE(known_y's，known_x's)函數(shù)可以利用現(xiàn)有的σ值和τ值計算直線的斜率 k，即抗剪包線斜率tanφ。known_y's表示因變的觀察值或數(shù)據(jù)集合;known_x's表示自變的觀察值或者數(shù)據(jù)集合。

2.3.2 INTERCEPT函數(shù)與 SLOPE函數(shù)聯(lián)合求解法的實現(xiàn)

設(shè)數(shù)據(jù)在 Excel表中的排列仍如圖2所示，則在C2和D2單元格中分別輸入“=INTERCEPT(B2:B5，A2:A5)”，“=ATAN(SLOPE(B2:B5，A2:A5))* 180/PI()”。求解出的結(jié)果為 c=3.260 9 kPa，φ= 23.297 1°。

2.4 基于點(diǎn)線距的規(guī)劃求解法

2.4.1 基本原理

如圖3所示，對抗剪強(qiáng)度庫侖公式τ=σtanφ+c，令tanφ=k，則任意測點(diǎn) A(σi，τi)到直線的距離為di［12］，見式(5)。

圖3 點(diǎn)線距求解基本原理示意

根據(jù)最小二乘法的原理，回歸分析時就是令 d2i和達(dá)到最小，即

式(6)中有k和c兩個變量，實際為線性規(guī)劃求解問題，求解目標(biāo)函數(shù)的最小值。Excel規(guī)劃求解一般步驟為建立目標(biāo)函數(shù)f(x)min，確定約束條件將決策變量用數(shù)學(xué)公式表達(dá)出來，使需要優(yōu)化的目標(biāo)及評判標(biāo)準(zhǔn)量化為這組變量的函數(shù)。

2.4.2 直剪試驗的規(guī)劃求解實現(xiàn)

Excel規(guī)劃求解首先要設(shè)計好一個工作表(如圖4所示)，根據(jù)建立的數(shù)學(xué)模型，將決策變量、約束條件、目標(biāo)函數(shù)等依次排列，并在工作表中將相關(guān)固定信息初始化，具體實現(xiàn)步驟如下。

圖4 基于點(diǎn)線距的規(guī)劃求解示例

1)打開Microsoft Excel軟件后，按圖4基于點(diǎn)線距的規(guī)劃求解所示的格式建立工作表，將有關(guān)數(shù)據(jù)輸入工作表中(其中C2和D2為空白，有待計算)，E2單元格輸入“=((C2*A2－B2+D2)/SQRT(C2^2+1)) ^2+((C2*A3－B3+D2)/SQRT(C2^2+1))^2+ ((C2*A4－B4+D2)/SQRT(C2^2+1))^2+((C2* A5－B5+D2)/SQRT(C2^2+1))^2”。

2)在“工具”欄下拉菜單中找到“規(guī)劃求解”，如果沒有此選項，可以選擇“加載宏”將“規(guī)劃求解”加載進(jìn)來。單擊“規(guī)劃求解”，出現(xiàn)“規(guī)劃求解參數(shù)”對話框，設(shè)置目標(biāo)單元格:$E$2，在“等于”選項中選“最小值”，在“可變單元格”中輸入$F$2:$G$2，在“約束”中單擊“添加”后，逐個輸入下式:“$C$2＞=0”，“$D$2＞=0”，最后按“確定”。當(dāng)需要強(qiáng)制 c值為0時，則條件變?yōu)椤?C$2＞=0”，“$D$2=0”。求解結(jié)果為c=3.254 6 kPa，k=0.430 6，在C3單元格中求得φ=23.298 8°，此時目標(biāo)函數(shù)最小，f(x)min=5.618 9。

2.5 基于測試值與擬合值誤差的規(guī)劃求解

直剪試驗數(shù)據(jù)σi，τi利用庫侖公式擬合直剪試驗中每組n個土樣的試驗數(shù)據(jù)時，c，φ被視為待定系數(shù)。以偏差最小為條件從而確定c，k值。這仍為線性規(guī)劃求解問題，求解目標(biāo)函數(shù)的最小值。其實現(xiàn)方法同2.4中方法類似，僅需要根據(jù)公式更改其目標(biāo)函數(shù)為“=(B2－A2*C2－D2)^2+(B3－A3*C2－D2)^2+ (B4－A4*C2－D2)^2+(B5－A5*C2－D2)^2”。

運(yùn)行規(guī)劃求解后，選擇“保存極限值”即可同時求得c值和k值，然后根據(jù)k值可計算出內(nèi)摩擦角φ值。根據(jù)上述方法，表1數(shù)據(jù)的求解結(jié)果為c=3.261 0 kPa，φ=23.297 1°，此時目標(biāo)函數(shù)最小，M=6.660 9。

3 各種直剪試驗求解方法的對比

為了對比5種求解方法的有效性，以表1中測試數(shù)據(jù)作為第1組數(shù)據(jù)，以福建省福清市江陰鎮(zhèn)海邊固結(jié)淤泥測試數(shù)據(jù)作為第2組數(shù)據(jù)，第2組數(shù)據(jù)如表2所示。

表2 福清江陰鎮(zhèn)淤泥土樣測試數(shù)據(jù)

通過在Excel中求解，5種方法求解結(jié)果如表3所示。根據(jù)表3對2組試驗的計算結(jié)果，分析如下:

1)常規(guī)數(shù)據(jù)計算結(jié)果對比。雖然各種方法均采用了最小二乘法，在不強(qiáng)制黏聚力為0的情況下，除基于點(diǎn)線距的Excel規(guī)劃求解法外，其余方法求解結(jié)果均相同，說明其余4種方法均是基于擬合誤差的求解方法。從數(shù)據(jù)對比可知5種方法計算結(jié)果非常接近，

均可作為直剪試驗的求解結(jié)果。

表3 Excel中5種直剪試驗求解方法計算結(jié)果對比

2)黏聚力為負(fù)值情況處理對比。表3中福州省福清市江陰鎮(zhèn)由于取樣撓動等原因，多組測試數(shù)據(jù)測試結(jié)果擬合后，黏聚力均為負(fù)值，這與理論并不相符。因此對于c值為負(fù)的情況或無黏性土的黏聚力為0的情況，此時需要強(qiáng)制黏聚力為0。對于強(qiáng)制c=0的情況，INTERCEPT與SLOPE函數(shù)聯(lián)合求解法無法完成強(qiáng)制c=0，其余方法均可強(qiáng)制c=0。通過對比可得，對于無黏性土或黏聚力出現(xiàn)負(fù)值情況，不宜采用INTERCEPT與SLOPE函數(shù)聯(lián)合求解法。

3)操作簡易程度對比。從5種方法在數(shù)據(jù)更換后Excel操作繁簡比較，規(guī)劃求解需要重新求解，特別是試驗組數(shù)更改后，目標(biāo)函數(shù)應(yīng)隨著更改，而且會生成中間報告，故對計算機(jī)操作不太熟悉人員有些難度。其它方法均操作簡便，且計算結(jié)果會自動完成更新。

4 結(jié)論

根據(jù)最小二乘法的基本原理及前人的研究成果，提出了基于Excel圖解法、INTERCEPT與SLOPE函數(shù)聯(lián)合求解法、TREND函數(shù)法、基于點(diǎn)線距的規(guī)劃求解法和基于擬合誤差的規(guī)劃求解法，重點(diǎn)分析了TREND函數(shù)與基于點(diǎn)線距規(guī)劃求解法2種方法的求解原理，并針對每種求解方法在Excel軟件中設(shè)計了相應(yīng)公式或?qū)崿F(xiàn)方法。通過幾種直剪試驗的求解方法對比分析，可以得出以下結(jié)論:

1)Excel中的 TREND函數(shù)求解法，截距函數(shù)INTERCEPT與斜率函數(shù) SLOPE聯(lián)合求解法均為Excel內(nèi)部函數(shù)法，大大地減小了編輯公式的工作量，是一種操作簡便計算方法，在無黏性土或計算黏聚力＜0的情況下，TREND函數(shù)求解法最為方便。

2)5種求解方法基本原理均采用基于最小二乘法的回歸分析，其擬合誤差保證其與實測值差值的平方和最小，其精度較采用兩點(diǎn)計算的斜截求解法和手工繪圖更能反映各測點(diǎn)對試驗結(jié)果的影響，更能代表土的實際抗剪強(qiáng)度指標(biāo)。

3)在5種求解方法中，當(dāng)計算的黏聚力 c值為正值時，各種方法求解結(jié)果相近，均可作為試驗計算結(jié)果。當(dāng)由于某種原因而引起擬合c值為負(fù)值或土樣為無黏性土?xí)r，c=0才正確。Excel圖解法、TREND函數(shù)法、兩種規(guī)劃求解法均可實現(xiàn)強(qiáng)制黏聚力為0，而INTERCEPT與SLOPE函數(shù)聯(lián)合求解法無法實現(xiàn)強(qiáng)制截距為0，故不能采用。

［1］李向國.滑帶土抗剪強(qiáng)度參數(shù)的確定［J］.西部探礦工程，2010(2):37-38.

［2］馬顯春，王雷，趙法鎖.滑坡穩(wěn)定影響因子敏感性分析及治理方案探討［J］.地質(zhì)力學(xué)學(xué)報，2008，14(4):381-388.

［3］劉遠(yuǎn)明，劉杰，夏才初.不同節(jié)理表面形貌下非貫通節(jié)理巖體強(qiáng)度特性直剪試驗研究［J］.巖土力學(xué)，2014，35(5): 1269-1273.

［4］衡帥，楊春和，曾義金，等.基于直剪試驗的頁巖強(qiáng)度各向異性研究［J］.巖石力學(xué)與工程學(xué)報，2014，33(5):874-883.

［5］陳建勝，陳從新，魯祖德，等.強(qiáng)風(fēng)化角巖力學(xué)—變形特性的直剪試驗研究［J］.巖土力學(xué)，2010，31(9):2869-2874.

［6］王曉斌，張鐵男，董志國.直剪試驗中的斜截式計算方法［J］.水利科技與經(jīng)濟(jì)，1999，5(2):82-83.

［7］南京水利科學(xué)研究院土工研究所.土工試驗技術(shù)手冊［M］.北京:人民交通出版社，2003.

［8］談云志，王世梅.直剪試驗數(shù)據(jù)處理方法的對比分析［J］.三峽大學(xué)學(xué)報(自然科學(xué)版)，2005，27(2):132-134.

［9］林立軍.二次曲線擬合法在無機(jī)結(jié)合料穩(wěn)定材料擊實試驗數(shù)據(jù)處理中的應(yīng)用［J］.福建建筑，2012，30(5):83-85.

［10］余東明，姚海林，吳少鋒.三軸試驗抗剪強(qiáng)度參數(shù)值回歸分析法的區(qū)別與修正［J］.巖土力學(xué)，2012，33(10):3037-3042.

［11］陳輝，顏勝才，孟軍濤.K30試驗數(shù)據(jù)的曲線擬合分析方法［J］.鐵道建筑，2012(11):67-70.

［12］阮波，張向京，彭意.Excel規(guī)劃求解三軸試驗抗剪強(qiáng)度指標(biāo)［J］.鐵道科學(xué)與工程學(xué)報，2009，6(5):57-60.

(責(zé)任審編 李付軍)

TU411

:ADOI:10.3969/j.issn.1003-1995.2015.09.29

2014-11-03;

:2015-02-10

國家自然科學(xué)基金項目(51208198);江西省教育廳項目(JXJG-13-5-15);華東交通大學(xué)科研基金(12RJ02)

莫林利(1977— )，女，河南安陽人，講師，碩士。

1003-1995(2015)09-0102-04