許小健

（蕪湖市勘察測繪設計研究院有限責任公司，安徽 蕪湖 241001）

0 引 言

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

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

1 基于多項式的擬合模型及求解方案

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

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

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

達到最小值時，bj值即為所求多項式系數值。根據求極值方法，此時由式（2）分別對bj求偏導數并令其等于0，然后得到關于bj的線性方程組（也稱為法方程），再進行線性方程組求解，從而求得多項式系數bj的值。

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

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

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

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

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

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

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

上述菜單項的各功能響應均通過VBA程序實現(xiàn)。

3 應用算例

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

表1 各數據組原始試驗數據Table 1 Original test data of each data group set

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

表2 test1數據組的計算結果Table 2 Calculation results of data group set test1

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

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

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

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

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

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

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

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

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

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

4 結 論

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

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

（3）基于Excel VBA的土工擊實試驗數據處理解決方案是一種精確、高效、直觀的數據處理和圖表可視化制解決方案，對于提高試驗和工程技術人員的工作嚴謹性、準確性和效率都有較好的作用；VBA編程門檻較低，值得工程技術人員掌握和學習使用。