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

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

(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ù)計算對比表明，除基于點線距的規(guī)劃求解法外，其余4種方法均為基于擬合誤差的計算方法，兩種原理的計算結(jié)果差距可忽略不計。對于無黏性土情況，除INTERCEPT函數(shù)與SLOPE函數(shù)聯(lián)合求解法無法實現(xiàn)強制黏聚力為0外，其余方法均可實現(xiàn)。實際應(yīng)用表明Excel圖解法、TREND函數(shù)法及截距和斜率函數(shù)聯(lián)合求解法操作簡單，易于應(yīng)用。

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

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

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

1 最小二乘法求解原理

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

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

直剪試驗求解中，當(dāng)繪制測點的破壞包線時，由于土樣不均勻，試驗讀數(shù)誤差等原因，往往繪制的直線不能通過所有的測點，總是存在一些誤差，根據(jù)這些誤差的平方和為最小，可以確定一條近似的破壞包線。設(shè)直剪試驗中每組n個土樣的垂直應(yīng)力、抗剪強度分別為σi，τi，利用庫侖公式 τ=σtanφ+c進(jìn)行擬合，其中c，φ被視為待定系數(shù)。以偏差+c)］2最小為條件來確定c，φ值。運用高等數(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軟件有強大的繪圖工具，先把試驗數(shù)據(jù)填入數(shù)據(jù)表格中，然后插入圖表，并為數(shù)據(jù)添加趨勢線和顯示公式，就可以得到直剪試驗圖。具體實現(xiàn)方法為:選擇正應(yīng)力和剪應(yīng)力數(shù)據(jù)，選擇“插入”→“圖表”→“散點圖”，然后選擇圖形中的點后“添加趨勢線”，在“選項”欄中選擇“顯示公式”，可直接得到擬合公式 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ù)求解的基本原理

由于直剪試驗各測點存在線性關(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ù)為邏輯值指明是否強制常數(shù)項c為0。當(dāng)擬合發(fā)現(xiàn)c值為負(fù)時，由于黏聚力不可能為負(fù)值，因此可以強制c值為0。

TREND函數(shù)求解原理如圖1所示，很顯然，對于抗剪強度包線τ=σ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ù)求解抗剪強度指標(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ù)值時，可以強制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)有的σ值和τ值計算直線與τ軸的截距，即抗剪強度指標(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 基于點線距的規(guī)劃求解法

2.4.1 基本原理

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

圖3 點線距求解基本原理示意

根據(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 基于點線距的規(guī)劃求解示例

1)打開Microsoft Excel軟件后，按圖4基于點線距的規(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)需要強制 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”。

運行規(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é)果對比。雖然各種方法均采用了最小二乘法，在不強制黏聚力為0的情況下，除基于點線距的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的情況，此時需要強制黏聚力為0。對于強制c=0的情況，INTERCEPT與SLOPE函數(shù)聯(lián)合求解法無法完成強制c=0，其余方法均可強制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ù)法、基于點線距的規(guī)劃求解法和基于擬合誤差的規(guī)劃求解法，重點分析了TREND函數(shù)與基于點線距規(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種求解方法基本原理均采用基于最小二乘法的回歸分析，其擬合誤差保證其與實測值差值的平方和最小，其精度較采用兩點計算的斜截求解法和手工繪圖更能反映各測點對試驗結(jié)果的影響，更能代表土的實際抗剪強度指標(biāo)。

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

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

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

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

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

［5］陳建勝，陳從新，魯祖德，等.強風(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］余東明，姚海林，吳少鋒.三軸試驗抗剪強度參數(shù)值回歸分析法的區(qū)別與修正［J］.巖土力學(xué)，2012，33(10):3037-3042.

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

［12］阮波，張向京，彭意.Excel規(guī)劃求解三軸試驗抗剪強度指標(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