土工擊實(shí)試驗(yàn)數(shù)據(jù)處理的解決方案

許小健

（蕪湖市勘察測(cè)繪設(shè)計(jì)研究院有限責(zé)任公司，安徽 蕪湖 241001）

0 引 言

土木工程中的建筑地基處理、基槽回填、公路路基修筑、試驗(yàn)檢測(cè)等，常遇到填土壓實(shí)問(wèn)題，而擊實(shí)試驗(yàn)用規(guī)范標(biāo)準(zhǔn)化的方法模擬了工程現(xiàn)場(chǎng)夯實(shí)原理，實(shí)現(xiàn)了對(duì)土的含水率與干密度關(guān)系的確定，從而為設(shè)計(jì)和施工提供依據(jù)[1]。因此，如何快速、便捷地確定擊實(shí)試驗(yàn)成果，具有重要實(shí)際工程意義。

目前，《土工試驗(yàn)方法標(biāo)準(zhǔn)》（GB/T 50123—2019）和《公路土工試驗(yàn)規(guī)程》（JTG 3430—2020）中均采用圖解法確定擊實(shí)試驗(yàn)成果。眾所周知，圖解法的優(yōu)點(diǎn)是工程技術(shù)人員可以在圖上直觀了解試驗(yàn)數(shù)據(jù)情況，是“規(guī)范”規(guī)定的方法，容易被工程技術(shù)人員所接受；缺點(diǎn)是作圖過(guò)程具有任意性，不同的人員繪制曲線和利用關(guān)系曲線圖讀取最大干密度和最優(yōu)含水率的結(jié)果受人為影響而具有不確定性。因此，尋求一種方便快捷、客觀可靠的試驗(yàn)數(shù)據(jù)處理方法是非常必要的。當(dāng)前，除作圖法外，廣大學(xué)者與工程技術(shù)人員對(duì)擊實(shí)試驗(yàn)數(shù)據(jù)處理方法的研究主要有兩類，一是擬合法[2-4]，二是數(shù)值分析法[5-8]。前者的研究實(shí)質(zhì)是以多項(xiàng)式為基礎(chǔ)利用擬合原理提供擬合解；后者是利用不同插值函數(shù)如Lagrange、Hermite、樣條曲線插值等數(shù)值分析方法處理試驗(yàn)成果。在數(shù)據(jù)處理的實(shí)現(xiàn)手段上，主要采用了Excel[9]、Mathcad[4]、Matlab[10]、Origin[11]、SPSS[12]等商業(yè)軟件進(jìn)行數(shù)據(jù)處理。以上方法不同程度地使得試驗(yàn)數(shù)據(jù)處理過(guò)程得到了改善。為提供一種新的解決方案，本文給出了一種利用VBA（Visual Basic for Applications，簡(jiǎn)稱VBA）結(jié)合Excel內(nèi)部函數(shù)和規(guī)劃求解工具快速自動(dòng)化處理?yè)魧?shí)試驗(yàn)數(shù)據(jù)的解決方案，該解決方案方便、快捷、準(zhǔn)確，且能實(shí)現(xiàn)繪圖成果可視化。

1 基于多項(xiàng)式的擬合模型及求解方案

《土工試驗(yàn)方法標(biāo)準(zhǔn)》（GB/T 50123—2019）和《公路土工試驗(yàn)規(guī)程》（JTG 3430—2020）中規(guī)定標(biāo)準(zhǔn)擊實(shí)試驗(yàn)采用作圖法，即以干密度為縱坐標(biāo)，含水率為橫坐標(biāo)，繪制干密度與含水率的關(guān)系曲線，再獲得曲線上峰值點(diǎn)的縱、橫坐標(biāo)分別為最大干密度ρdmax和最優(yōu)含水率ωop。

為避免人工處理試驗(yàn)數(shù)據(jù)，這里將上述規(guī)范作圖法轉(zhuǎn)換為數(shù)學(xué)求解問(wèn)題。一般情況下，對(duì)于每組試驗(yàn)的m=5個(gè)土樣的含水率和干密度擊實(shí)試驗(yàn)數(shù)據(jù)(ωi, ρi)(i=1, 2, …, m)，可用n次(n=2, 3, 4)多項(xiàng)式曲線模型進(jìn)行描述，也即：

式中：j=1, 2, …, n。為確定多項(xiàng)式曲線模型，需要對(duì)多項(xiàng)式系數(shù)進(jìn)行求解，常用的方法是采用最小二乘原理[6-7]，即當(dāng)殘差平方和：

達(dá)到最小值時(shí)，bj值即為所求多項(xiàng)式系數(shù)值。根據(jù)求極值方法，此時(shí)由式（2）分別對(duì)bj求偏導(dǎo)數(shù)并令其等于0，然后得到關(guān)于bj的線性方程組（也稱為法方程），再進(jìn)行線性方程組求解，從而求得多項(xiàng)式系數(shù)bj的值。

“法方程”的求解是確定多項(xiàng)式模型的關(guān)鍵，其求解也稍復(fù)雜，可以通過(guò)自行編程實(shí)現(xiàn)法方程的求解，但對(duì)于非計(jì)算機(jī)專業(yè)的工程技術(shù)人員來(lái)說(shuō)，用C、C++、Fortran等語(yǔ)言進(jìn)行編程將面臨較大的困難，開(kāi)發(fā)不易。文獻(xiàn)[13]雖然采用了VBA編程，但其通過(guò)采用矩陣計(jì)算求解“法方程”，本質(zhì)上其編程與采用其他編程語(yǔ)言實(shí)現(xiàn)難度無(wú)異。鑒于Excel為用戶提供了豐富的工作表函數(shù)，為簡(jiǎn)化編程，本文采用Excel內(nèi)部工作表函數(shù)LinEst函數(shù)直接實(shí)現(xiàn)多項(xiàng)式曲線模型系數(shù)的提取，這樣可以避免對(duì)最小二乘法線性方程組編程求解，從而最大限度地減輕編程的復(fù)雜度。LinEst工作表函數(shù)通過(guò)使用最小二乘法計(jì)算與現(xiàn)有數(shù)據(jù)實(shí)現(xiàn)最佳擬合，并返回描述擬合問(wèn)題的數(shù)組，其函數(shù)調(diào)用語(yǔ)法為L(zhǎng)inEst ([yValues],[xValues], [const], [stats])，其調(diào)用參數(shù)為：（1）yValues為必需參數(shù)，xValues為可選參數(shù)，對(duì)于擊實(shí)試驗(yàn)來(lái)說(shuō)，分別相當(dāng)于干密度和含水率的集合；（2）const和stats為可選參數(shù)。其中stats用于指定是否返回?cái)M合回歸統(tǒng)計(jì)值，如判定系數(shù)r2、殘差平方和SSE等。具體的使用方法可通過(guò)參閱Excel幫助文檔，也可參考本文后續(xù)編程實(shí)現(xiàn)部分。

對(duì)于多項(xiàng)式擬合模型，一旦其系數(shù)確定后，其模型也就確定，但為獲得最大干密度ρdmax和最優(yōu)含水率ωop的值，還需要對(duì)模型進(jìn)行進(jìn)一步求解。由于多項(xiàng)式模型為高次非線性方程，求解不易，而Excel為非線性問(wèn)題的求解提供了規(guī)劃求解工具Solver，使得用戶較易獲得問(wèn)題解。非線性規(guī)劃求解工具Solver采用“Generalized Reduced Gradient (GRG2)”代碼，該代碼由Leon Lasdon（德克薩斯大學(xué)奧斯汀分校）和Alan Waren（克利夫蘭州立大學(xué)）開(kāi)發(fā)，并由Frontline Systems公司封裝后提供。用戶無(wú)需關(guān)心具體的代碼，只需編寫少量的VBA程序進(jìn)行簡(jiǎn)單的調(diào)用即可在Excel中實(shí)現(xiàn)復(fù)雜問(wèn)題的求解。

綜上所述，本文提出的解決方案是通過(guò)在Excel中運(yùn)行VBA程序，先利用Excel內(nèi)置工作表函數(shù)LinEst函數(shù)直接提取多項(xiàng)式擬合模型的系數(shù)，再調(diào)用Excel規(guī)劃求解工具Solver實(shí)現(xiàn)多項(xiàng)式擬合模型的極值求解，從而確定最大干密度和最優(yōu)含水率。

2 基于Excel VBA的解決方案實(shí)現(xiàn)

VBA是一種具有簡(jiǎn)單易用特點(diǎn)的編程語(yǔ)言，其編程系統(tǒng)中加入了面向?qū)ο蟮臋C(jī)制，把Windows編程的復(fù)雜性程序進(jìn)行了封裝，并提供了一種可視界面的設(shè)計(jì)方法。同時(shí)，VBA程序賴以運(yùn)行的應(yīng)用程序開(kāi)發(fā)工具VBE（Visual Basic Editor，簡(jiǎn)稱VBE）寄生于Office環(huán)境。由于Office軟件的安裝與使用較為普遍，故VBA的運(yùn)行和存放也無(wú)需特殊的環(huán)境設(shè)置，這進(jìn)一步降低了編程開(kāi)發(fā)的門檻。

為方便使用，本文基于Excel設(shè)計(jì)了程序用戶界面，如圖1所示，該界面保留了與Excel界面的風(fēng)格統(tǒng)一和美觀。Excel運(yùn)行時(shí)，程序會(huì)在Excel的菜單欄形成一個(gè)“土工擊實(shí)試驗(yàn)”菜單項(xiàng)，該菜單項(xiàng)包含相應(yīng)功能選項(xiàng)。

圖1 程序用戶界面Fig. 1 Application user interface

上述菜單項(xiàng)的各功能響應(yīng)均通過(guò)VBA程序?qū)崿F(xiàn)。

3 應(yīng)用算例

為展示本文解決方案的便捷性和有效性，現(xiàn)以表1中擊實(shí)試驗(yàn)數(shù)據(jù)組舉例說(shuō)明。其中test1試驗(yàn)數(shù)據(jù)來(lái)源于文獻(xiàn)[7]，test2和test3試驗(yàn)數(shù)據(jù)來(lái)源于文獻(xiàn)[10]，test4試驗(yàn)數(shù)據(jù)來(lái)源于文獻(xiàn)[4]，test5來(lái)源于文獻(xiàn)[15]。

表1 各數(shù)據(jù)組原始試驗(yàn)數(shù)據(jù)Table 1 Original test data of each data group set

啟動(dòng)Excel后，默認(rèn)自動(dòng)打開(kāi)圖1所示“土工擊實(shí)試驗(yàn)數(shù)據(jù)處理”選項(xiàng)卡。以表1中test1為例，在Excel表格的B2:B6單元格內(nèi)依次輸入含水率試驗(yàn)數(shù)據(jù)ωi，在C2:C6單元格內(nèi)依次輸入干密實(shí)試驗(yàn)數(shù)據(jù)ρi；輸入完成后，點(diǎn)擊“模型擬合計(jì)算與極值求解”，程序會(huì)自動(dòng)批量依次完成n=2、n=3、n=4階多項(xiàng)式曲線模型的擬合與極值求解工作，并顯示結(jié)果如表2所示。表2中，Excel表格F～J列中自動(dòng)填充多項(xiàng)式系數(shù)結(jié)果；在K和L列分別顯示多項(xiàng)式擬合附加回歸統(tǒng)計(jì)值r2和SSE；在M和N列分別顯示最優(yōu)含水率ωop和最大干密度ρdmax。經(jīng)過(guò)比對(duì)，本文處理test1試驗(yàn)數(shù)據(jù)的結(jié)果與文獻(xiàn)[7]的4次多項(xiàng)式曲線擬合模型結(jié)果完全一致，且與ωop=20.98%，ρdmax=1.781 g/cm3也完全相同。但文獻(xiàn)[7]的做法需要在Matlab軟件的命令窗口中根據(jù)軟件反饋信息不斷與軟件交互輸入，便捷性稍差，而本文則實(shí)現(xiàn)了在Excel環(huán)境中通過(guò)一次簡(jiǎn)單點(diǎn)擊操作就能立即顯示數(shù)據(jù)自動(dòng)化處理結(jié)果，操作簡(jiǎn)單便捷。

表2 test1數(shù)據(jù)組的計(jì)算結(jié)果Table 2 Calculation results of data group set test1

以test2、test3試驗(yàn)組數(shù)據(jù)為例，為便于討論，提供求解最大干密度與最優(yōu)含水率結(jié)果見(jiàn)表3。限于篇幅，本文所求各階多項(xiàng)式曲線模型系數(shù)不再展示，后同。由表3可見(jiàn)，本文所得結(jié)果與文獻(xiàn) [10]利用Matlab商業(yè)軟件處理結(jié)果相比，也幾乎一致，僅test2數(shù)據(jù)組在多項(xiàng)式階數(shù)較高（n=3、n=4）時(shí)所求得含水率精度上略有差異，但最大誤差也不超過(guò)0.15%，可以忽略不計(jì)，其原因是在計(jì)算過(guò)程中對(duì)數(shù)據(jù)精度的設(shè)置差異所致。

表3 各方法所得結(jié)果的比較（test2、test3數(shù)據(jù)組）Table 3 Comparison of results obtained by different methods (data set test2, test3)

再以test4數(shù)據(jù)組為例，文獻(xiàn)[4]中提供了多項(xiàng)式階數(shù)n=3的求解情況，如表4所示?？梢?jiàn)，本文處理test4試驗(yàn)數(shù)據(jù)的結(jié)果（n=3階）與文獻(xiàn)[4]采用Mathcad商業(yè)軟件所得結(jié)果完全一致，與《公路土工試驗(yàn)規(guī)程》（JTJ 051—93）圖解法結(jié)果也較為接近。由此可見(jiàn)，本文解決方案便捷、有效。

表4 各方法所得結(jié)果的比較（test4數(shù)據(jù)組）Table 4 Comparison of results obtained by different methods (data set test4)

最后以test5數(shù)據(jù)組為例，在多項(xiàng)式階數(shù)n=2時(shí)，各方法所得結(jié)果均較為一致。在多項(xiàng)式階數(shù)n=3、n=4時(shí)，文獻(xiàn)[15]的Newton迭代計(jì)算擬合法所求結(jié)果與本文所得結(jié)果差異較大，如表5所示。不難發(fā)現(xiàn)，對(duì)于多項(xiàng)式階數(shù)n=3時(shí)，文獻(xiàn)[15]中擬合法和數(shù)解插值法所求最大干密度1.63 g/cm3和1.59 g/cm3，明顯小于擊實(shí)試驗(yàn)實(shí)測(cè)數(shù)據(jù)的1.67 g/cm3，較為不合理。對(duì)于多項(xiàng)式階數(shù)n=4時(shí)，文獻(xiàn)[15]求得最優(yōu)含水率為14.2%，明顯低于其圖解法結(jié)果18.6%和本文所求結(jié)果18.63%，也不合理。為便于觀察，利用設(shè)計(jì)的“導(dǎo)出擊實(shí)試驗(yàn)成果圖”功能模塊，繪制n=4時(shí)的含水率與干密度關(guān)系曲線圖如圖2中（a）所示，可見(jiàn)文獻(xiàn)[15]的ωop=14.2%不符合關(guān)系曲線趨勢(shì)，同樣可以得出其不合理的結(jié)論。此外，從判定系數(shù)r2值來(lái)看，本文所求擬合模型效果更好。

表5 各方法求得結(jié)果的比較（test5數(shù)據(jù)組）Table 5 Comparison of results obtained by different methods (data set test5)

為避免產(chǎn)生上述test5所得試驗(yàn)結(jié)果偏差較大的情況發(fā)生，利用程序可以自動(dòng)批量繪制、導(dǎo)出含水率與干密度關(guān)系曲線圖的功能，可以起到輔助用戶直觀掌握最大干密度和最優(yōu)含水率情況，同時(shí)也可以起到橫向驗(yàn)證數(shù)據(jù)處理結(jié)果與成果圖所得圖解結(jié)果一致性的作用。繪圖過(guò)程無(wú)需人工進(jìn)行調(diào)整，保留了圖解法的直觀。

值得說(shuō)明的是，本文所提解決方案在處理test4和test5試驗(yàn)數(shù)據(jù)時(shí)，發(fā)現(xiàn)當(dāng)n=4時(shí)，Excel規(guī)劃求解工具Solver容易陷入局部解，所求最大干密度偏低，最優(yōu)含水率結(jié)果易產(chǎn)生變動(dòng)；而在n=2、n=3時(shí)，均未發(fā)現(xiàn)此現(xiàn)象。經(jīng)過(guò)多次試驗(yàn)求解發(fā)現(xiàn)，這與規(guī)劃求解工具Solver求解迭代計(jì)算的初始值有關(guān)。而軟件缺乏解決實(shí)際問(wèn)題對(duì)初始值的選擇，這是軟件所采用的計(jì)算方法所致。上述文獻(xiàn)[15]中所求test5試驗(yàn)組所得不合理結(jié)果亦應(yīng)屬此原因所致。為解決這一問(wèn)題，考慮到含水率的取值要存在合理的實(shí)際意義，所以其初值的選擇也并非難事，只要將自變量迭代初始調(diào)整為試驗(yàn)組含水率數(shù)據(jù)的平均值即能較好地解決該問(wèn)題，通過(guò)在“模型擬合計(jì)算與極值求解”源程序調(diào)用規(guī)劃求解工具Solver時(shí)添加cells(4,"M").Formula= "=AVERAGE(B2:B6)"語(yǔ)句即可解決。當(dāng)然，也可以選擇試驗(yàn)數(shù)據(jù)組已經(jīng)獲取的最大干密度所對(duì)應(yīng)的含水率作為自變量迭代初始值，這也是一種可選擇的方案。

此外，值得一提的是，由文中5組試驗(yàn)組數(shù)據(jù)情況發(fā)現(xiàn)，在多數(shù)情況下，當(dāng)多項(xiàng)式階數(shù)為n=2時(shí)，擬合效果稍差，如圖2（b）中所示含水率與干密度關(guān)系圖曲線峰值點(diǎn)，其擬合計(jì)算所得最大干密度較實(shí)測(cè)數(shù)據(jù)偏低。當(dāng)n=3、n=4時(shí)，所求最大干密度與最優(yōu)含水率的結(jié)果更為符合原始數(shù)據(jù)規(guī)律，但n=4時(shí)，含水率與干密度趨勢(shì)曲線的兩側(cè)尾部易產(chǎn)生“翹曲”，圖2（a）中曲線尾部?jī)啥撕腿鐖D2（c）中曲線左側(cè)尾部所示，n=3時(shí)，則容易取得數(shù)據(jù)成果和所繪制曲線趨勢(shì)特征的“均衡”，如圖2（d）所示。

圖2 含水率與干密度關(guān)系曲線圖Fig. 2 Relation curve between moisture content and dry density

4 結(jié) 論

（1）基于多項(xiàng)式曲線擬合原理，將圖解法轉(zhuǎn)換為數(shù)學(xué)求解問(wèn)題，再利用Excel規(guī)劃求解功能獲得擊實(shí)試驗(yàn)成果的解決方案，克服了人為因素的影響，避免了人工圖解法處理試驗(yàn)數(shù)據(jù)較為繁瑣和具有任意性的劣勢(shì)。

（2）實(shí)現(xiàn)了與Excel風(fēng)格統(tǒng)一的可視化程序界面，程序界面美觀、便捷、高效，實(shí)現(xiàn)了通過(guò)簡(jiǎn)單點(diǎn)擊操作就能進(jìn)行數(shù)據(jù)自動(dòng)化處理和成果圖繪制工作。

（3）基于Excel VBA的土工擊實(shí)試驗(yàn)數(shù)據(jù)處理解決方案是一種精確、高效、直觀的數(shù)據(jù)處理和圖表可視化制解決方案，對(duì)于提高試驗(yàn)和工程技術(shù)人員的工作嚴(yán)謹(jǐn)性、準(zhǔn)確性和效率都有較好的作用；VBA編程門檻較低，值得工程技術(shù)人員掌握和學(xué)習(xí)使用。