摘要：土的先期固結壓力（p_c）是描述土體受力歷史的重要指標，用卡薩格蘭德（Casagrande）作圖法求解土的先期固結壓力通常受人為因素影響且計算較為復雜。為此，通過研究卡薩格蘭德經驗作圖法已有成果，提出了一種改進方法，并用計算機編程加以實現。該方法將采集的數據點分段定義為貝塞爾曲線的起點和終點，求解出貝塞爾曲線定義式，然后再計算每一段曲線曲率點，最后相互比較即可得出曲線最大曲率點。該方法通過編制相應VBA及EXCEL程序便于工程應用。工程算例表明，該方法降低了人為因素影響，對試驗數據具有較好的適應性，可在巖土工程的土體先期固結壓力確定中推廣應用。

關 鍵 詞：先期固結壓力；貝塞爾插值；最小曲率半徑；Casagrande作圖法

中圖法分類號：TU43

文獻標志碼：ADOI：10.16232/j.cnki.1001-4179.2024.12.023

0 引 言

先期固結壓力p_c是指土層在地質形成歷史上所經受的最大豎向固結壓力，是判斷土體應力歷史的一個重要指標。在具有不同應力歷史的土層變形分析中，先期固結壓力是一個重要的計算參數^［1］。p_c的推求方法大致有3種^［2］：卡薩格蘭德（Casagrande）法、布麥斯脫方法、薛邁脫曼方法。為了解決卡薩格蘭德法中受人為和比例因素影響的問題，一些學者提出了使用數學模型法來確定最大曲率點^［3-5］。為了彌補數學模型法計算復雜的不足，劉林等^［6］建立了適用于正常固結黏土和超固結黏土的壓縮線表達式，然后基于該表達式，推導出先期固結壓力。隨著機器學習（ML）算法的廣泛應用，李超等^［7］采用集成學習算法來建立先期固結壓力預測模型，使用貝葉斯優(yōu)化方法來確定模型的先期固結壓力。Butterfield^［8］提出將e～lgp半對數坐標轉換為ln（1+e）～lgp雙對數坐標來描述土的壓縮性，由雙對數坐標可得到兩條直線，對應兩直線交點的應力即為先期固結壓力p_c。并且Onitsuka ^［9］和Hong ^［10］等通過大量的試驗驗證了雙對數法的有效性。Alibrahim等^［11］根據壓縮系數與壓縮指數推導出一種非線性數學公式，用于擬合先期固結壓力。Zhai等^［12］采用圖解法，通過FX模型確定先期固結壓力。

本文以武漢地區(qū)不同地層土樣為研究對象，選取卡薩格蘭德法，在分析了既有e～lgp曲線數學模型不足的基礎上提出了一種用計算機確定先期固結壓力p_c值的方法。利用VBA語言編寫計算程序，將試驗數據導入計算機即可求解得到土體先期固結壓力p_c值，并自動完成圖解畫圖，可以消除目測定點的人為性，提高求解精度。

1 卡薩格蘭德經驗作圖法回顧

目前，規(guī)范^［13］普遍采用卡薩格蘭德法確定先期固結壓力，其原理是以實驗室固結壓力試驗數據所描繪的離散點為基礎，連接各離散點繪制圓滑e～lgp曲線，利用e～lgp曲線曲率突變點來推求p_c。有很多學者使用計算機對于如何在曲線上找到曲率突變點開展研究，其方法主要分為三類：① 多項式擬合或插值法，如蔡清池^［14］引入e=0.42e₀點，運用Harris模型及多項式擬合求取p_c值；彭靜等^［15］采用等跨度圓弧求圓半徑及拉格朗日多項式插值法求取最小曲率半徑，并通過Matlab語言編制了相應的計算程序求取p_c值；劉用海等^［16］采用3次多項式和最小二乘法回歸擬合壓縮曲線直接求其曲率函數的方法推求p_c值。② 斜率變化法，如劉衍成^［17］用相鄰三點得兩條直線，通過段比較相鄰直線斜率變化幅度大小的方法確定曲率最大點，得到p_c值。③ 三點定圓法，如張書憲^［18］提出相鄰三點逐一比較圓半徑的大小確定曲率最大點；周國云^［19］提出先對試驗數據點插值后以三點求圓半徑比大小的方法。

由試驗所得的離散點繪制而成的e～lgp曲線，理論上呈順時針彎曲狀，且原始的試驗數據點與整條曲線的特性存在關聯性。采用上述多項式擬合或插值法得到的曲線會丟失原始試驗數據點或使曲線發(fā)生順時針和逆時針交替彎曲的狀態(tài)，而且在數據間隔較大的曲線段擬合曲線和插值曲線明顯較卡薩格蘭德法所得曲線彎曲程度要大，這些都改變了卡薩格蘭德的方法所得e～lgp曲線的特征，如圖1所示。

采用上述斜率變化法及三點定圓法，也會由于試驗數據點分布的間距不一而受到干擾。間距很近的3個數據點，中間點位置的變化對所確定的圓弧半徑或兩直線的斜率的影響程度明顯要比間距較大的情形大得多。如圖2所示左邊由間距大的點構成的圓弧，當中間點偏移形成新圓弧時，新舊兩圓弧的曲率變化幅度不大；右邊相距較近的三點構成的圓弧，當中間點偏移形成新圓弧時，新舊兩圓弧的曲率變化幅度很大。由于試驗數據存在誤差，上述現象總是存在。因此斜率變化法、三點定圓法也存在較大缺陷。

2 貝塞爾插值求解以及自動作圖實現

本文著眼于用計算機解決卡薩格蘭德經典圖解法求解土體先期固結壓力的問題，因此需要保持卡薩格蘭德經典圖解法的圖形幾何特性，從而遵從其圖解原理設計出一種新算法。

2.1 貝塞爾曲線分段插值

貝塞爾曲線最早應用于二維圖形相關設計中，其由起始點、終止點和控制點來定義。1962年，法國數學家Pierre Bézier第一個研究了這種矢量繪制曲線的方法，并給出了詳細的計算公式，按照這樣的公式繪制出來的曲線就用他的姓氏來命名，故稱為貝塞爾曲線。該曲線的最大特點是由離散點構成的曲線在每一段連接處是平滑且連續(xù)的，連接處兩端的曲線曲率不會出現較大的變化，因此該方法既保留了e～lgp曲線的特征，又可避免上述三種方法的缺陷。

設P₀、P²₀、P₂是一條拋物線上三個順序不同的點，如圖3所示。過P₀和P₂點的兩切線交于P₁點，在P²₀點的切線交P₀P₁和P₂P₁于P¹₀和P¹₁，則有公式（1）比例成立：

P₀P¹₀/P¹₀P₁=P₁P¹₁/P¹₁P₂=P¹₀P²₀/P²₀P¹₁（1）

當P₀，P₂固定，引入參數t，令上述比值為t：（1-t），即有：

P¹₀=1-tP₀+tP₁（2）

P¹₁=1-tP₁+tP₂（3）

P²₀=1-tP¹₀+tP¹₁（4）

將式（2）和式（3）代入式（4）得：

P²₀=1+t²+2t1-tP₁+t²P₂（5）

當t從0變到1時，它表示了由三頂點P₀、P₁、P₂定義的一條二次貝塞爾曲線。并且點P₂的運動軌跡可以看作前兩個點（P₀、P₁）和后兩個點（P₁、P₂）的線性組合。

依次類推，由4個頂點（P₀、P₁、P₂、P₃）定義的曲線為三次貝塞爾曲線，如圖4所示。

圖4 三次貝塞爾曲線

Fig.4 The cubic Bessel curve

由（n+1）個控制點P_i（i=0，1，…，n）定義的n次貝塞爾曲線P₀ⁿ可被定義為分別由前后n個控制點定義的兩條（n-1）次Bezier曲線P₀^n-1與P₁^n-1的線性組合：

Pⁿ_i=1-tP^n-1₀+tP^n-1₁ t∈［0，1］（6）

n階貝塞爾曲線的通用公式為：

P^k_i=/P_i

（k=0）（1-t）P^k-1_i+tP^k-1_i+1 （k=1，2，…，n）/i=0，1，…，n-k

（7）

由于需要求解e～lgp曲線最大曲率點，采用三次貝塞爾曲線可保證曲線的光滑性和連續(xù)性。假設室內固結試驗采集到的數據點為P₀，P₁，…，P_n，將這些數據點分為n-1段，即L_i=P_iP_i+1（i=0，1，…，n-1）。每一分段的前后兩點（P_i，P_i+1）定義為三次貝塞爾曲線的起點和終點，求解三次貝塞爾曲線的定義式，即可求得每一段最大曲率點，最后比較每一段的最大曲率點，求解出整個e～lgp曲線的最大曲率點。

上述可知，三次貝塞爾曲線是由前后4個控制點定義的。因此在實際試驗中，當某個采集數據點由于誤差偏離真實值時，三次貝塞爾曲線受到誤差的影響要小于斜率變化法及三點定圓法，如圖5所示。

2.2 基于VBA程序實現先期固結壓力求解和自動作圖

（1）根據以上理論，編寫VBA程序，從EXCEL中啟動VBA，將EXCEL對應的單元中數據導入VBA程序中，并用這些數據完成相應計算。將計算結果自動輸送到EXCEL預先設定功能應用的單元格中。

（2）VBA中的貝塞爾插值函數通過相互嵌套調用來求解e～lgp曲線的最大曲率點和該點切線斜率，并輸出到EXCEL單元格。

（3）在EXCEL中調用前期固結試驗數據的統(tǒng)計模板和卡薩格蘭德經典圖解繪制應用模版。在VBA程序中完成相應運算后將結果輸回EXCEL時，單元格預設功能函數將自動完成作圖。

3 計算實例

試驗樣品取自武漢某堤段勘察項目鉆孔原狀樣，地層代號如下：3-4淤泥質土，3-5粉質黏土、粉土、粉砂互層土（簡稱互層土），7-1黏土，每層3組，呈灰色、灰褐色。按照GB/T 50123—2019《土工試驗方法標準》^［13］用環(huán)刀切取試樣，試樣尺寸直徑61.8 mm，高度20 mm。試驗設備為單杠桿固結儀，加壓荷載分別為6.25，12.5，25，50，100，200，300，400，600，800，1 200，1 600，2 400，3 200 kPa。試樣在每級荷載下固結不少于24 h，待變形穩(wěn)定后再施加下一級荷載，直至最大荷載固結穩(wěn)定，試驗過程中室內無振動。

根據室內固結試驗成果，利用e～lgp貝塞爾插值求解法求取先期固結壓力，并采用作圖法及l(fā)n（1+e）～lgp雙對數法進行了驗證^［8］，其土層物理指標及計算成果見表1，軟件自動計算成圖見圖6～9。

由表1試驗成果可知：3-4層淤泥質土，3-5層粉質黏土、粉土、粉砂互層土，7-1層黏土，三種不同深度土樣均處于欠固結狀態(tài)，和實際情況較為吻合。同時，由貝塞爾插值求解法與作圖法得到的先期固結壓力p_c值相對誤差為-2.3%～1.5%，貝塞爾插值求解法與雙對數法得到的先期固結壓力P_c值相對誤差為-2.7%～0.1%，這說明本文采用貝塞爾插值求解先期固結壓力的方法是可行的，能夠滿足現行規(guī)范推薦方法要求。

4 結 論

本文在對國內外各學者提出的先期固結壓力確定方法開展研究的基礎之上，給出了求解e～lgp曲線最大曲率點的計算方法，通過VBA語言編制了相應的計算程序。該方法有以下優(yōu)點：

（1）求得的曲線均通過試驗數據點，避免了多項式擬合或插值法得到的曲線會丟失原始試驗數據點的情況。

（2）本文方法的每一段曲線均與下一段曲線相切，保證整個曲線的光滑性，采用三次曲線即可滿足求解曲率函數的要求，避免了高階多項式插值法使曲線發(fā)生順時針和逆時針交替彎曲突變的情況。

（3）算例表明本文方法對試驗數據具有較好的適應性，且與EXCEL兼容，計算速度快，操作簡便，但該方法還有待進一步的完善和試驗數據的檢驗。

參考文獻：

［1］MATLOCK H S.Correlation for design of laterally loaded pile in soft clay［C］∥Proceedings of 2nd Offshore Technology Conference，Houston：［s.n.］，1970.

［2］錢家歡，殷宗澤.土工原理與計算［M］.北京：中國水利水電出版社，1996.

［3］趙明志，羅強，蔣良濰，等.原狀土e-lgp曲線的S形函數表達及簡化加載固結試驗方法［J］.實驗力學，2018，33（5）：769-780.

［4］孫政，何振華，陳宗輝，等.長沙粉質黏土的結構性與先期固結壓力研究［J］.湖南文理學院學報（自然科學版），2019，31（3）：64-70.

［5］安之煥，陳杰，唐占元，等.土體初始結構強度的表征及其數值確定方法［J/OL］.水利水電技術（中英文）.http：∥link.cnki.net/urilid/10.1746.tv.20240128.1434.004.

［6］劉林，魯曉兵，王淑云.確定前期固結壓力的一個簡單數學模型［J］.土木工程學報，2019，52（增2）：30-34.

［7］李超，汪磊，陳洋，等.基于貝葉斯集成學習算法的土體先期固結壓力預測模型［J］.地球科學，2023，48（5）：1780-1792.

［8］BUTTERFIELD R.A nature compression law for soils［J］.Geotechnique，1979，29（4）：469-480.

［9］ONITSUKA K，HONG Z，HARA Y，YOSHITAKE S.Interpretation of oedometer test data for natural clays［J］.Soils and Foundations，1995，35（3）：61-70.

［10］HONG Z，ONITSUKA K.A method of correcting yield stress and compression index of Ariake clays for sample disturbance［J］.Soils and Foundations，1998，38（2）：211-222．

［11］ALIBRAHIM B，UYGAR E.Nonlinear calculation method for one-dimensional compression of soils［J］.Arabian Journal for Science and Engineering，2022，47（4）：4865-4877.

［12］ZHAI Q，ZHU Y，RAHARDJO H，et al.Mathematical model for the determination of the pre-consolidation pressure using the graphical method［J］.Arabian Journal of Geosciences，2023，16（1）：69.

［13］住房和城鄉(xiāng)建設部.土工試驗方法標準：GB/T 50123-2019［S］.北京：中國計劃出版社，2019.

［14］蔡清池，謝漢康.一種確定先期固結壓力的數值計算方法［J］.工業(yè)建筑，2021，5（5）：164-167

［15］彭靜，張勝軍，何嬌，等.基于計算機求解土的先期固結壓力［J］.資源環(huán)境與工程，2015，29（2）：202-205.

［16］劉用海，朱向榮，常林越.基于Casagrande法數學分析確定先期固結壓力［J］.巖土力學，2009（1）：211-214.

［17］劉衍成.前期固結壓力的計算和繪圖程序［J］.工程勘察，1987（3）：10-15.

［18］張書憲.用計算機繪圖確定先期固結壓力的一種方法［J］.巖土工程界，2000，8（3）：45-47.

［19］周國云.一種確定前期固結壓力p_c的電算方法［J］.水文地質工程地質，1987（5）：34-36.

（編輯：鄭 毅）

Solution of pre-consolidation pressure of soil based on Bessel interpolation

QIN Haijun¹，ZHANG Shengjun²，YANG Mingwei²

（1.Gansu Civil Engineering Research Institute Co.，Ltd.，Lanzhou 730000，China； 2.Yangtze River Survey and Technology Research Institute，Ministry of Water Resources，Wuhan 430011，China）

Abstract： The pre-consolidation pressure （p_c） of soil is an important index that describes the force history of soil.Solving the pre-consolidation pressure of soil by Casagrande drawing method is usually affected by human factors and complicated to calculate.Through studying Casagrande empirical drawing method，we put forward an improved calculation method that is realized by computer programming.In this method，the collected data points are defined as the beginning and end points of the Bessel curves segmentally，and the definition equations of the Bessel curves are given.We then calculate the curvature points of each section of curves，and finally we derived the maximum curvature points by comparing each curvature point.The new method is compiled into corresponding calculation programs by VBA and EXCEL programming to facilitate the application.It shows that this method can reduce the influence of human factors and possesses better adaptability to test data，thus it is worthy of promotion in determining pre-consolidation pressure in geotechenique engineering.

Key words： pre-consolidation pressure；Bessel interpolation；minimum curvature radius；Casagrande drawing method