粉質黏土未凍水含量預測模型及參數變化規律研究

鄧世磊, 萬旭升, 路建國, 李雙洋, 晏忠瑞

(1.西南石油大學 土木工程與測繪學院,四川 成都 610500; 2. 中國科學院西北生態環境資源研究院 凍土工程國家重點實驗室, 甘肅 蘭州 730000)

0 引言

我國是世界上的第3凍土國，多年凍土面積約占國土面積的22.3%，主要分布于西部的青藏高原和東北的大小興安嶺等地[1]。凍土是一種特殊性質的土，主要由土顆粒和冰、液態水和氣體等對溫度極其敏感的成分組成，它在結構、物理和化學性質上與一般土種均有較大差異。在土體凍結之后，并非所有液態水都會轉換成冰，土體中始終保留著部分液態水，稱之為未凍水[2]。在多年凍土地區施工時，應充分考慮溫度變化對地基的影響[3]。在四季交替變化過程中溫度會不斷變化，所引起的凍土中未凍水含量的變化會對土壤的力學性質產生巨大的影響。從而對當地的建筑、交通工程及基礎造成諸多安全問題，如房屋塌陷、路面隆起和下沉、地基承載力和穩定性明顯變化等。所以，未凍水含量是凍土相關研究的關鍵指標，其隨溫度變化規律在數值模擬中發揮重要作用。

近年來，隨著測量手段和信息技術的發展，未凍水含量的測試方法也變得多種多樣，測試結果越精確，得出的結論更具有科學性。未凍水測量常用方法有量熱法、脈沖核磁共振法(NMR)、頻域反射法(FDR)、頻域傳播法(FDT)及其衍生的多種測量手段[4]，其原理是不同液態水含量對應不同的信號值，可通過不同溫度下的信號值建立未凍水含量與溫度的關系。在大量未凍水含量試驗基礎上，Anderson等[5]基于試驗數據提出了描述未凍水含量隨溫度變化的冪函數模型。徐敩祖等[6]在此模型基礎上，提出了采用一點法和二點法來預測未凍水含量的方法。Sheng等[7]考慮未凍水含量指數函數變化規律，建立了土的凍脹預測模型。冷毅飛等[8]測試了俄石油管道工程沿線大量土樣的未凍水含量，并用此模型進行擬合，得到了良好的效果。

土體的基本物理參數及特征研究為未凍水含量預測提供了新的思路。Dall’Amico等[9]結合土水特征曲線和土體凍結曲線，建立了未凍水含量和溫度的關系式。Topp等[10]建立了一個用介電常數計算土中未凍水含量的三次多項式，并用試驗值確定了常數項的取值。之后，Liu等[11]在試驗數據的基礎上提出了一種利用土壤表觀介電常數和初始鹽濃度來計算土壤未凍水含量的經驗方法。萬旭升等[12]基于理想成冰模型，推導出了利用土壤孔隙比計算土中未凍水含量的計算公式，并引入有效應力原理，建立了負溫下未凍水含量和溫度的關系式。李順群等[13]基于比熱隨負溫的變化過程和潛熱的宏觀表現形式，建立了未凍水含量的反演算法。張翻等[14]根據Johansen導熱系數計算式，反演得到了不同溫度下細粒土的未凍水含量計算公式。Chai等[15]提出了粉質黏土中考慮冰點的未凍水含量物理計算模型，并計算了土中體積水、毛細水和結合水的冰點。Mu等[16]考慮了毛細作用和吸附作用，提出了此模型參數化的計算方法。自然界中凍土的壓力環境復雜，基于此，Zhou等[17]建立了不同壓力下的未凍水含量計算模型。

綜上所述，已有未凍水含量模型種類繁多，但是在未凍水含量計算時，大多需要借助儀器設備測量土體的物理參數，且模型中參數隨土的初始物理狀態變化較大，變化規律沒有統一表述，模型參數難以確定。故本研究選擇3種典型的模型進行研究，通過對粉質黏土的大量試驗數據進行擬合，研究模型參數變化規律，并討論模型的計算精度。最后，給出各模型使用的建議溫度范圍。

1 經典未凍水含量模型

1.1 Anderson和 Tice模型

Anderson和Tice根據未凍水含量隨溫度變化的關系，建立了一個冪函數模型來計算未凍水含量，計算公式如下：

(1)

式中，Wu為質量未凍水含量；W0為初始質量含水率；a和b為模型參數，無實際物理意義；T為土壤溫度，Tf為土壤凍結溫度，無其他因素影響時，一般認為Tf=0 ℃。由于冪函數的數學特性，溫度接近冰點時出現奇異值。

1.2 Michalowski模型

Michalowski考慮了土的殘余水含量，并根據試驗數據提出了以下未凍水含量計算模型：

(2)

式中，Wr為低溫下的殘余質量未凍水含量；μ為模型參數。該模型與冪函數模型相比，解決了當溫度靠近凍結溫度時未凍水含量出現無限大的問題。

1.3 Dall’Amico模型

Dall’Amico在土水特征曲線的基礎上，考慮土體凍結特征，結合Clapeyron方程，給出了低溫下的未凍水含量計算模型：

(3)

(4)

(5)

式中，θu為體積未凍水含量；θ0為初始體積含水率；θr為殘余體積未凍水含量；φ(T)為不同溫度下的土壤基質勢；φw0為總含水率(液體和冰)相對應的基質勢(當土為飽和土時φw0=0)；Lf為熔化潛熱(333.7 kJ/kg)；Tm為大氣壓下水的融化溫度(273.15 K)；g為重力加速度(9.81 m/s2)；α，n，m均為模型參數。在本研究的計算中，將土體近似于飽和狀態，且不考慮凍結溫度的影響。

為了將單位統一化，將所有試驗數據的含水率均轉化為體積含水率，質量含水率與體積含水率的轉換關系為：

(6)

式中，θ為體積含水率；W為質量含水率；ρd為土的干密度；ρw為水的密度。

2 研究方法

2.1 粉質黏土基本物理參數

為了讓分析結果更具有科學性，排除個別特殊數據組對結果產生影響，本研究采集了不同區域的粉質黏土試驗數據進行擬合研究[4,12-13,15,18-23]。組別前面的字母代表不同的文獻來源，后面的數字為數據的組內編號。其中A，B，C，D各為1組數據，H組內為2組數據，F組內為3組數據，其余的E，G，I，J組內各為5組數據，共29組試驗數據。其中，E組數據的初始含水率使用土壤水分儀測量，其他數據組采用烘干法測量，干密度均采用環刀法測量。不同區域的粉質黏土的礦物成分、液塑限和顆粒級配會有所差異。本研究中引用的粉質黏土的液塑限和粒徑分布如表1所示。

表1 液塑限及粒徑分布表Tab.1 Table of liquid limit, plastic limit and particle size distribution

從表1的數據可知，本研究選取的粉質黏土塑限位于12.5～18.6之間，液限位于23.40～36.70之間，塑性指數位于10.9～17.94之間，變化范圍較小。從粒徑分布可以看出，粉質黏土的粒徑主要集中在0.1～0.001 mm之間，顆粒級配較為均勻。

2.2 數據擬合及誤差分析

為了分析模型的精度，繪出試驗數據的T-θu散點圖，采用origin軟件中的非線性擬合的方法，將以上試驗數據組分別用3種模型進行擬合，確定各模型參數的取值。并引入均方根誤差(RMSE)和均差(AD)作為模型精度的評價手段，均方根誤差將直接反映模型的整體預測效果，而均差將反映模型是否會過低或過高地預測未凍水含量。均方根誤差和均值的計算公式為：

(7)

(8)

式中，n為每組數據中數據點數；θe為試驗所測得的體積未凍水含量；θc為未凍水含量模型計算所得的體積未凍水含量。

3 結果與分析

3.1 模型邊界條件

參考各參數的取值，且取Tf=0 ℃，在計算過程中，當溫度接近0 ℃(凍結溫度)時，各模型的邊界條件如下：

(9)

(10)

(11)

從各模型的邊界條件可以看出，當溫度趨近于0 ℃ 時，Anderson 和 Tice模型的預測值接近無限大，而Michalowski模型和Dall’Amico模型的值為土樣的初始含水率。所以在計算均方根誤差和均差時，去除了0 ℃的函數值。Michalowski模型和Dall’Amico模型規避了Anderson和Tice模型在溫度接近0 ℃ 時值出現無限大的現象。

3.2 模型參數對凍結曲線影響

本研究選擇I-1數據樣本進行模型參數對溫度-未凍水含量曲線影響分析，探究模型參數變化時對3種模型擬合曲線的影響。I-1土壤樣本初始體積含水率25.43%，干密度1 560 kg/m3，凍結溫度Tf=0 ℃。在I-1組數據中，模型參數a=15.600，b=-0.537，μ=0.518，α=0.009，n=2.6，m=0.615。改變各模型參數的數值，觀察T-θu曲線形態的變化。

通過改變3種模型中參數的數值，可以對比出在擬合過程中各參數對擬合曲線形態的影響。從圖1(a)中可以看出，在Anderson和 Tice模型中，當a值不變時，b值變化對凍結初始階段(T>-2 ℃)無明顯影響，對后段凍結階段(T<-2 ℃)曲線形態產生影響，b值減小時后段凍結曲線整體下移，b值增大使后段凍結曲線整體上移；當b值不變時，a值會影響曲線的整體上下移動，a值增大會讓曲線整體上移，減小時整體下移。從圖1(b)中可以看出，在Michalowski模型中，參數μ的變化會引起曲線下凹部分上下移動，μ值增大時凍結曲線下凹部分下移，減小時下凹部分上移。從圖1(c)可以看出，在Dall’Amico模型中，參數α，n，m對凍結曲線有類似的影響，其值變化時曲線下凹部分發生變化，當α值不變時，n值增大會使曲線下凹部分上移，減小時下移；相反，當n，m值不變時，α值增大使凍結曲線下凹部分下移，減小時上移。

圖1 參數變化對各模型曲線的影響Fig.1 Influence of parameter change on each model curve

3.3 模型精度分析

通過對試驗數據的擬合，可以得出各數據組的參數取值。不同模型都有不同的擬合特點和效果，為了對比3種模型的計算精度和擬合效果，計算出在各數據組下的RMSE和AD值。畫出了均方根誤差的柱狀圖(圖2)和均差的散點圖(圖4)，比較和分析3種模型的預測準確性。

圖2 RMSE柱狀分布圖Fig.2 Columnar distribution of RMSE

從圖2中可看出，Anderson和Tice模型的均方根誤差計算值明顯低于其他2種預測模型。為了讓結果更清晰化，計算了29組數據的RMSE的平均值，結果如表2所示。

從表2中可以看出，Anderson 和Tice模型的RMSE均值為0.874%，顯著低于Michalowski模型的1.676%和Dall’Amico模型的1.356%。綜上所述，Anderson 和Tice模型的整體預測效果優于其他2種模型，且其計算精度已經非常高，計算數值在實際數值上下浮動僅為0.874%，而Dall’Amico模型略微優于Michalowski模型。

表2 RMSE計算平均值(單位：%)Tab.2 Calculated average values of RMSE (unit: %)

通過原始數據的擬合曲線(圖3)可以發現，C-1組數據的RMSE值明顯高于其他數據組，在該土樣中在凍結溫度附近大量水分凍結，在-0.5 ℃時未凍水含量發生驟降，在-7 ℃時凍結基本完成[14]。凍結溫度附近水分快速凍結及凍結過早完成導致擬合效果較差，數據組的RMSE值偏高。

圖3 C-1組數據擬合曲線Fig.3 Data fitting curves of C-1

從AD散點分布圖(圖4)可以看出，在Anderson和Tice模型中，AD負值僅出現7個，其余22個均為正值，說明Anderson 和Tice模型會高估土中的未凍水含量；而在Michalowski模型和Dall’Amico模型中，AD計算值的負值個數分別為23個和20個，直觀地表明了Michalowski模型和Dall’Amico模型在實際使用中會低估土體的未凍水含量。從圖4中折線的形態可以看出，Anderson 和Tice模型的AD值在0點附近僅有微小的波動，整體預測效果優于Michalowski模型和Dall’Amico模型。相比之下，Michalowski模型和Dall’Amico模型在0點處起伏較大，估值誤差較大。

圖4 AD散點分布圖Fig.4 Scattergram of AD

3.4 模型的整體預測效果評價

通過對3種模型進行粉質黏土數據的擬合分析，得出了3組模型在不同數據組內的參數取值，并計算了均方根誤差和均差,對擬合結果進行了分析評價。最后計算結果一致表明，Anderson 和Tice模型的計算精度明顯高于Dall’Amico模型和Michalowski模型，在使用過程中會高估土中的未凍水含量，而Dall’Amico模型計算精度略微優于Michalowski模型，且二者都會低估未凍水含量。

4 模型參數變化規律

在擬合過程中，Anderson 和 Tice模型未凍水含量預測模型具有形式簡單、擬合方便和應用范圍廣泛等優點，自從提出以來，被廣大學者推崇使用。Michalowski模型只有1個擬合參數μ，相較于其他2個模型擬合更為簡單。Dall’Amico模型雖然有3個擬合參數，但是m和n相關，且其擬合參數的變化較小，易于確定。

4.1 模型參數的取值范圍

在未凍水含量計算模型中，參數影響函數曲線的形態，確定了凍結過程中未凍水含量隨溫度變化的關系。但是隨著粉質黏土的初始含水率、干密度、顆粒級配和孔隙比的變化，模型參數的值也會改變。將上述29組T-θu試驗數據進行匯總，在數據點集集中區域的上下方分別確定不同模型中的上下邊界，求出上下界曲線的模型參數。與實際擬合得出的參數取值進行對比，并計算參數取值范圍的涵蓋率(涵蓋率為位于上下邊界取值區間內的數據組數與總數據組數的百分比)，結果如表3所示。

表3 上下邊界參數取值Tab.3 Parameter values of upper and lower boundaries

從表3可以看出，在Anderson和Tice模型中，模型參數a的取值范圍為7～45，b的取值范圍為-0.185～-0.65，在29組數據中，參數擬合取值中a最大值為43.589，最小值為7.520。b最大值為-0.185，最小值為-0.733，除了b的最小值外，其他數據均介于上下邊界的b值之間，參數a取值涵蓋率為100%，參數b的涵蓋率為96.6%。在Michalowski模型中，擬合參數μ的取值范圍為0.208～1.200，除D-1組的參數取值2.934外，其他數據樣本取值均位于取值范圍之內，參數μ取值涵蓋率為96.6%。在Dall’Amico模型中，參數α取值范圍為0.002 5～0.012，參數n為2.0～2.8(參數m與n具有一致性)，參數α有6組數據取值位于范圍之外，除E組的5組數據外，參數n的取值全位于范圍之內，參數α的取值涵蓋率為79.3%，n和m的涵蓋率為82.6%。從各模型參數的取值涵蓋率可以看出，由上下邊界確定的參數取值范圍都具有較佳的指導意義。

4.2 初始含水率對模型參數的影響

在不同地區的粉質黏土中，含水率有著明顯的差異，所以，探究模型參數取值與初始含水率的關系尤為重要。下面選取6組來自同一篇文章的數據進行對比分析，保證基本物理性質相同，排除其他干擾因素的影響，分析初始含水率對參數取值的影響。

從圖5(a)～(b)可以看出，在Anderson 和Tice模型中，在E，F，G組土樣數據內，隨著初始含水率增加，參數a的值也在增加，b值減小。在H組的2組土樣數據內，隨著初始含水率的增加，參數a的值增大，b值增大。在I組和J組土樣數據內，隨著初始含水率的增加，參數a在增大，但是b值沒有出現固定增大或減小的現象，其值出現波動的現象。從圖5(c)可以看出，在Michalowski模型中，除了J組數據中隨著初始含水率的增大出現了μ取值減小外，其他5組數據均出現μ取值隨著初始含水率增大而增大的現象。而圖5(d)～(e)反映出在Dall’Amico模型中，參數α和n的取值與初始含水率的關系規律性較差，整體上，α的取值隨初始含水率增大而出現減小的趨勢，n的取值隨著初始含水率的增大而出現減小的趨勢，但二者都有反常的增大和減小。由此可以總結出，在相同的土體中，初始含水率的增大會直接導致a值的增大，μ值有增大的趨勢，但是對b，α，n，m的取值沒有固定的影響，b，α，n，m值的變化可能還與土的凍結過程緊密相關。

圖5 參數取值與初始含水率關系Fig.5 Relationship between parameter value and initial water content

5 討論

以上分析表明，Anderson和 Tice模型計算精度較高，但在凍結溫度附近不收斂，而Michalowski模型和Dall’Amico模型在凍結溫度處的計算值為初始含水率，且在凍結結束階段計算值更接近殘余水含量。根據3個模型的計算誤差，現擬將3個模型根據最佳溫度適用范圍，分溫度區段結合使用。以未凍水含量測試溫度點一致的F組和I組試驗數據為例，對比各模型在各溫度點的計算值與實際值的差值，繪出AD絕對值分布圖，如圖6所示。

圖6 AD絕對值分布Fig.6 Distribution of absolute values of AD

根據差值的分布的特點，將0～-15 ℃溫度區間從-2 ℃和-10 ℃分成3個階段，依次為劇烈相變階段、主要凍結階段和凍結穩定階段。從圖6中可以看出，Dall’Amico模型在0～-2 ℃劇烈相變階段的計算精度較高；Anderson 和Tice模型在主要凍結階段，即-2～-10 ℃的計算值更接近實際值；而Michalowski模型在凍結穩定階段精度高于其他2個模型。基于此，現將3個模型根據以上溫度區間的計算特性，把土中未凍水含量計算分成3個區段，提出未凍水含量分段計算公式：

(12)

為了驗證式(12)的可行性，現用F組與I組的數據進行計算精度校核，對比3個經典未凍水含量模型和式(12)的RMSE值，計算結果如表4所示。

表4 F組和I組RMSE計算平均值(單位：%)Tab.4 Calculated Average Values of RMSE in group F and group I (unit: %)

從表4的F組和I組數據的RMSE計算結果來看，在3個模型中，Anderson 和 Tice模型計算精度較高，RMSE值僅為0.654%，擬合效果較佳，其他2個模型分別為1.140%和1.013%。在使用式(12)對土中的未凍水含量進行分段預測之后，整體的RMSE計算平均值降至0.572%，整體預測效果優于其他3個模型。由此可以得出，式(12)結合了3個模型的優點，能更好地模擬土中未凍水含量隨溫度的變化。

土的液塑限、顆粒級配等會影響未凍水含量變化。已有研究表明，塑性指數對土體凍結特征曲線過冷階段無明顯影響；塑性指數的增大會減緩土中水分凍結的速率，使凍結階段曲線變得平緩,對凍結完成階段未凍水含量幾乎無影響[24]。土中孔徑影響未凍水含量的變化，未凍水含量改變隨著孔徑減小而減小，土中孔徑大小與顆粒級配密切相關，顆粒越小，孔徑越小，水分凍結愈發困難[6]。然而，土中未凍水含量模型參數與土體基本物理量之間還未建立具體關系，大部分參數物理意義不明確，故在后期研究中要考慮土的物理參數建立新的未凍水含量預測模型。

6 結論

本研究對已有粉質黏土試驗數據進行擬合分析，通過計算模型的均方根誤差RMSE和均差AD，對Anderson和Tice提出的冪函數模型、Michalowski考慮初始含水率和殘余水含量的指數函數模型，Dall’Amico結合Clapeyron方程提出的未凍水含量預測模型進行評價，歸納了模型參數的取值范圍，討論了參數取值和初始含水率的關系。結論如下：

(1)對于特定土而言，Anderson 和 Tice模型具有較高的預測精度，但是在冰點處出現奇異值，誤差較大，且隨初始含水率變化，預測值離散性較大。Dall’Amico模型整體性較好，可較好地反映不同初始含水率變化規律，而Michalowski模型模擬效果相比較差。Anderson和Tice模型會高估土中未凍水含量，Michalowski模型和Dall’Amico模型會低估未凍水含量。

(2)Anderson 和 Tice模型中，參數a取值范圍變化最大，為5～45，反映了模型隨土的不同初始物理狀態波動較大；Michalowski模型中，參數μ取值范圍為0.208～1.200；Dall’Amico模型中，參數α變化最小，取值范圍為0.0025～0.012，n為2.0～2.8，反映出模型穩定性較好。土的初始含水率的增大會引起參數a和μ的值增大，對b，α，n，m值沒有明確的影響。

(3)基于各模型的計算特性建立的分段計算公式結合了三者的優點，劇烈相變階段使用Dall’Amico模型計算，主要凍結階段和穩定凍結階段分別使用Anderson和Tice模型和Michalowsk模型，可較好地預測土中的未凍水含量。