青藏工程走廊凍融土熱擴散系數特性與預測模型研究

劉志云， 鐘振濤， 崔福慶， 陳建兵， 彭 惠

（1.長安大學地質工程與測繪學院，陜西西安710054； 2.中交第一公路勘察設計研究院有限公司高寒高海拔地區道路工程安全與健康國家重點實驗室，陜西西安710075）

0 引言

青藏工程走廊作為最重要的陸路進出藏通道，穿越了約550 km 多年凍土區，特殊的工程地質條件引發了凍脹、融沉、開裂、擁包、塌陷及不均勻變形等病害［1-2］。在未來的青藏高速公路工程建設中，高速公路黑色瀝青路面特有的強烈吸熱、聚熱效應將給路基下伏凍土帶來更大的熱沖擊［3］，誘發更為嚴重的工程病害［4-5］。如何保障這一特殊工程場景下公路工程的安全性和穩定性，是寒區凍土工程研究重要課題之一［6］。熱擴散系數是表征土體內溫度擴散速率的物理量，反映了傳熱過程中導熱系數和比熱容兩大基礎參數協同作用的效果［7］。同時，凍土的熱擴散系數一定程度上也是凍土本身對外部熱量輸入響應敏感性的重要標尺，是描述多年凍土特征的關鍵參數之一［8］。因此，揭示多年凍土熱擴散系數的變化規律、提出適用的預測模型對于未來青藏高速的設計施工和病害防治具有重要意義。

熱擴散系數主要獲取方法有理論計算和試驗測試，其中理論計算是基于土壤為常熱擴散系數的半無界介質假設，以一維熱傳導方程為理論基礎獲得計算公式，主要包括振幅法、相位法、對數法、反正切法、數值法和諧波法等［9-10］。許多研究者針對不同計算方法進行了對比分析，Horton 等［11］從計算所需數據條件和計算結果方面對比評價了六種熱擴散系數計算方法，發現具有顯式方程的振幅、相位、反正切和對數方法所需地溫數據少，但計算結果偏差大，而以隱式求解的數值法和諧波法通常以大量實測地溫為計算基礎，結果更為可靠。繆育聰等［12］根據現場觀測地溫數據計算了地表淺層5～20 cm 土壤層的熱擴散系數，對比發現充分利用實測地溫信息的諧波法為最優估算方法。韓炳宏等［13］對青海南部高寒草地土壤的熱擴散系數進行了計算，研究表明熱傳導對流法計算結果及其擬合效果最好，同時發現除干土層外隨土壤深度加深，熱擴散系數逐漸遞減。同時，不少學者根據土壤的傳熱特征對計算方法進行了修正，Gao 等［14-16］綜合考慮了土壤熱傳輸過程中的熱傳導和多孔熱對流方式，給出了一維熱傳導對流方程，采用諧波法和拉普拉斯變換方法推導了熱傳導對流方程的解析解及熱擴散系數的計算公式，并基于實測地溫數據驗證了改進方法的準確性和優越性。原黎明等［17］基于改進的熱傳導對流方法對青藏高原中部活動層表層下5～20 cm 深度的土壤熱擴散系數進行了估算，結果發現5～10 cm 處土壤熱擴散系數顯著大于10～20 cm 深度，且融化季對應熱擴散系數顯著高于凍結季。章永輝等［18］基于青藏高原理塘區野外觀測地溫數據，研究發現利用耦合熱傳導-對流法計算所得熱擴散系數模擬地溫時的準確性最高。在試驗方面，甄作林等［19］利用瞬態平面熱源法測試了原狀和重塑砂土的熱擴散系數，分析發現熱擴散系數與干密度成正比關系，而含水率低時熱擴散系數與含水率為正相關性，達到一定含水率時熱擴散系數趨于穩定或與其呈負相關性。

關于熱擴散系數影響因素方面，研究主要集中于土壤質地、內部結構、含水率、溫度及深度等。王可里等［20］通過對青藏高原那曲地區非均質土壤熱擴散系數計算發現，熱擴散系數分布具有明顯的深度和季節性變化特征，且冷暖季采用不同熱擴散系數可獲得較好地溫模擬效果。周亞等［21］根據青藏高原觀測站長期監測所得0.8 m 和3.2 m 深度土壤溫度，給出了青藏高原不同地區深層土壤熱擴散系數的年際和季節性變化特征。劉經星等［22］利用數值差分求解方法對不同質地條件土壤的熱擴散系數進行了估算，結果發現土壤質地輕且土熵條件好對應熱擴散系數小。Roxy［23］、安可棟等［24］根據觀測數據研究了土壤水分和熱擴散系數之間關系，研究表明熱擴散系數隨土壤水分增加先增大后減小。邸佳穎等［25］通過室內試驗對比研究了原狀土和裝填土的熱特性差異，結果表明裝填土在處于中等含水率時熱擴散系數較原狀土有所增大，而接近飽和時熱擴散系數趨于一致。馬欣等［26］研究發現土壤熱擴散系數與深度為非線性正相關性，且降水增多會降低熱擴散系數。在計算模型方面，Arkhangelskaya 等［27］以土壤質地、容重和有機碳百分含量作為影響變量，建立了不同土類任意含水率狀態下的熱擴散系數回歸模型。

但應看到的是，目前凍土熱擴散系數的研究報道還相對較少，缺乏針對青藏工程走廊帶典型類別土樣熱擴散系數的大規模測試與預測模型研究。為了進一步揭示青藏工程走廊帶沿線多年凍土的熱物理性質，本文基于青藏公路改擴建（格爾木—拉薩段）工程地質勘察項目912 組土樣導熱系數測試結果和質量加權法計算獲取的比熱容理論值，計算獲得西大灘—唐古拉山沿線典型類別土樣熱擴散系數，分析了凍融土熱擴散系數的分布特征和參數影響規律，提出了基于經驗擬合公式法和RBF 神經網絡方法的凍融土熱擴散系數預測模型，并對不同預測模型的預測效果進行了對比分析，以期為青藏高速公路的熱工設計提供數據參考。

1 熱擴散系數計算方法

熱擴散系數定義為導熱系數與體積比熱的比值，是反映物體增溫快慢的綜合物性參數，計算公式如下：

式中：α為土壤熱擴散系數；λ為土壤導熱系數，W·m-1·K-1；C為體積比熱，kJ·m-3·K-1。

本文青藏工程走廊重塑凍融土導熱系數采用Hot Disk TPS1500 熱常數分析儀以瞬態平面熱源法進行測試，其中測試區段北起西大灘（即走廊帶多年凍土北界）南至唐古拉山共437 km，對應青藏公路里程樁號為K2870～K3307。鉆孔取樣深度自地表向下最深達40 m，共計獲取土樣912組，其中測試融土樣633組、凍土樣858組。

對鉆探獲取的土樣進行土性統計，將土樣數超過15 組的列為典型土樣，按照土的工程分類，可將走廊帶內所取土樣分為黏性土、粉土、全風化巖類、砂土和碎石土共計5 大類12 個亞類，各土類統計如圖1所示。由圖可知，走廊帶內主要以黏性土為主，砂土和碎石土次之，且二者數量相近。

圖1 青藏工程走廊導熱系數測試區段及土性分布Fig.1 Thermal conductivity test section and soil property distribution of Qinghai-Tibet engineering corridor

導熱系數測試主要操作包括：（1）試樣制備：土樣烘干篩分、增濕至天然含水率、靜置悶樣及加壓制樣；（2）測試預處理：平整試樣表面和預先凍結凍土試樣；（3）物性測試：確定測試時間及功率、重復測試和數據誤差對比。其中試樣制備依據《土工試驗方法標準》（GB/T50123—2019）規范內容進行，具體測試系統如圖2所示。

圖2 導熱系數測試系統Fig.2 Thermal conductivity test system

體積比熱按照文獻［28］所給公式進行計算，該研究認為土比熱是多相組分的質量加權平均結果，凍融土比熱因固態冰的存在具有明顯差異，具體計算式如下：

式中：Cu和Cf分別為融土和凍土的體積熱容量，kJ·m-3·K-1；Csu、Csf、Cw和Ci分別為融土骨架、凍土骨架、水及冰的質量比熱，kJ·kg-1·K-1；W和Wu分別為總含水量和未凍水含量；ρu和ρf分別為融土和凍土的天然密度，kg·m-3。其中不同土類凍融土骨架及水分質量比熱取值見表1所示。

表1 不同土類土骨架及水質量比熱取值Table 1 Mass specific heat of soil skeleton and water for different soil types

2 凍融土熱擴散系數分布特征及影響因素分析

2.1 凍融土熱擴散系數分布特征

對五類土熱擴散系數進行統計，得到其概率分布如圖3 所示。由圖3 可知，融土熱擴散系數主要分布區間整體依黏性土、粉土、全風化巖類、砂土和碎石土偏右分布，表明融土熱擴散系數主要分布范圍值依次增大。可以看出，土的成分對融土的熱擴散系數有顯著影響，粒徑越小，土的持水性越強，則結合水的比例越高，降低了水分在礦物骨架之間的聯系作用，從而導致其熱擴散系數降低。為排除由取樣隨機性帶來的統計誤差，取累積概率分布區間20%～80%熱擴散系數為對比值，統計得黏性土、粉土、全風化巖類、砂土和碎石土融土熱擴散系數分布區間為0.441×10-6～0.595×10-6、0.485×10-6～0.694×10-6、 0.486×10-6～0.677×10-6、0.589×10-6～0.829×10-6和0.653×10-6～0.858×10-6m2·s-1，均值分別為0.507、0.565、0.586、0.699 和0.759×10-6m2·s-1。可以看到，走廊帶內各土類融土熱擴散系數與粒徑呈正相關性，分布基本排序為黏性土、粉土、全風化巖類、砂土及碎石土依次增大。

圖3 土類熱擴散系數概率分布Fig.3 Probability distribution of thermal diffusivity of soils

同理，圖3 中各土類凍土熱擴散系數主要分布區間整體按黏性土、全風化巖類、粉土、碎石土和砂土依次增大，其累積概率分布區間20%～80%對應的 熱 擴 散 系 數 區 間 為0.742×10-6～1.057×10-6、0.826×10-6～1.134×10-6、0.810×10-6～1.213×10-6、0.996×10-6～1.374×10-6和0.998×10-6～1.410×10-6m2·s-1，均值 分 別 為0.882×10-6、0.976×10-6、0.985×10-6、1.172×10-6和1.194×10-6m2·s-1。對比發現，各土類凍融土熱擴散系數分布規律不盡相同，可能與水分凍結引發土體結構變化有關，總體來說，走廊帶內土類凍土熱擴散系數分布基本排序為黏性土、全風化巖類、粉土、碎石土及砂土依次增大。

對比各土類凍土與融土熱擴散系數發現，凍土熱擴散系數均顯著大于融土，其中黏性土、粉土、全風化巖類、砂土和碎石土均值相對增幅分別為73.96%、74.34%、66.55%、70.82%和54.41%，這是因為凍結狀態下土體內水變為冰后具有更強的導熱能力和較小的比熱。同時，還可發現粒徑較小土類整體具有較大相對增幅，可能與土體凍結過程土骨架結構連接性增強相關，粒徑較小土類表現出更大增幅的導溫性。

2.2 土體天然含水率、干密度對熱擴散系數的影響

對走廊帶內凍融土熱擴散系數與天然含水率、干密度之間的偏相關性進行統計分析，結果如表2所示。由表可知，黏性土、粉土和碎石土凍融土熱擴散系數與干密度呈顯著的正相關關系，這是因為含水率不變的情況下，其質量比熱也不變，而導熱系數隨干密度變大而增大，故隨著干密度變大其凍融土熱擴散系數相應增大。但也存在全風化巖類凍融土及砂土凍土熱擴散系數與干密度為負相關性情況，由圖4可知，在自然環境狀況下不同干密度條件對應熱擴散系數分布隨機且離散，并不具有較強的正相關性。

圖4 全風化巖類融土熱擴散系數與干密度及天然含水率分布關系Fig.4 Relationship between thermal diffusivity and dry density and natural water content distribution of unfrozen soil of fully weathered rocks

表2 熱擴散系數與影響因素偏相關分析Table 2 Partial correlation analysis between thermal diffusivity and influencing factors

同時，可發現各土類凍融土熱擴散系數與含水率間關系差異顯著，大部分土類熱擴散系數與含水率呈負相關性，但粉土、砂土及碎石土凍土熱擴散系數與含水率呈正相關性，同樣如圖3所示，這可能與天然條件下土類處于不同賦水狀態時，熱擴散系數隨含水率變化規律不同相關。此外，還可發現土體天然含水率與其干密度有強烈的負相關性，凍融土熱擴散系數之間具有顯著正線性相關性。

3 凍融土熱擴散系數預測

3.1 凍融土熱擴散系數預測模型

3.1.1 二元擬合

對走廊帶內凍融土熱擴散系數與干密度、天然含水率之間的擬合關系進行曲線估計后發現，基本都呈為多項式函數形式，得到走廊帶內凍融土熱擴散系數一般擬合公式如下：

式中：α為土體熱擴散系數，10-6m2·s-1；ω為天然含水率，%；ρd為干密度，g·m-3；a1、b1、c1、d1、e1和f1均為待定的擬合系數，各土類擬合值見表3所示。

表3 熱擴散系數二元擬合相關參數Table 3 Related parameters of binary fitting for thermal diffusivity

3.1.2 基于融土熱擴散系數的凍土三元擬合

青藏高原走廊帶內土類分布多樣且分散，各土類礦物成分及組分粒徑等分布界限并不一定明顯，土類間存在相互混雜現象，所以簡單將其分類后采用干密度和含水率對其熱擴散系數進行預測會產生一定偏差。此外，相對于凍土而言，融土試樣制備及測試更為簡單，且融土熱擴散系數結果隱含有決定土性成分和粒徑等信息，前述偏相關分析表明二者具有極強的正線性相關性，同時考慮干密度、含水率對凍土熱擴散系數預測結果進行修正，得到走廊帶內凍土熱擴散系數的三元擬合公式如下：

式中：αf和αu分別為凍土和融土熱擴散系數，10-6m2·s-1；a2、b2、c2和d2均為待定的擬合系數，各土類擬合結果見表4所示。

表4 熱擴散系數三元擬合相關參數Table 4 Related parameters of ternary fitting for thermal diffusivity

3.1.3 徑向基函數（RBF）神經網絡

徑向基函數（radical basis function）神經網絡具有計算量小、學習收斂速度快及泛化能力強的特點，其較強的非線性函數逼近能力在諸多領域得以廣泛應用［29-30］。將干密度、天然含水率及融土熱擴散系數（凍土熱擴散系數三元回歸）作為RBF 神經網絡的輸入層，熱擴散系數作為輸出層，建立以高斯函數為隱含層激活函數的前饋型網絡模型。同時，將走廊帶內凍融土熱擴散系數測試結果按9：1隨機分塊，其中90%樣本數據用于神經網絡模型的訓練樣本，剩余10%用作驗證神經網絡的預測能力。根據預測結果的反饋，對偏差較大的樣本進行剔除，從而有效提高預測精度。圖5 為全風化巖類融土熱擴散系數的RBF 神經網絡模型預測結果，可以看到預測值與實測值較為一致（R2=0.78，數據剔除比例為9%）。

圖5 全風化巖類融土熱擴散系數RBF神經網絡預測結果Fig.5 Prediction results of thermal diffusivity based on RBF neural network of unfrozen soil of fully weathered rocks

3.2 融土熱擴散系數預測模型對比

圖6為各土類融土熱擴散系數兩種預測模型的預測結果，可以看到兩種預測模型所得預測值與實測結果整體較為吻合，大多數樣本點分布在相對誤差為10%范圍內，證明了預測方法的有效性和工程應用價值。此外，對比兩種方法預測值與實測值間差異分布，可發現RBF神經網絡模型整體預測效果更佳。

圖6 各土類融土熱擴散系數預測結果對比Fig.6 Comparison of prediction results of thermal diffusivity of unfrozen soil of different soil types

3.3 凍土熱擴散系數預測模型對比

圖7為各土類凍土熱擴散系數四種預測模型的預測結果，同樣可看到四種預測模型具有良好的預測性能。同時，對比四種方法預測值分布范圍及與實測值間差異可發現，將融土熱擴散系數納入回歸模型的三元預測方法較二元預測方法具有更寬的樣本分布范圍，且整體預測精度更高。

圖7 各土類凍土熱擴散系數預測結果對比Fig.7 Comparison of prediction results of thermal diffusivity of frozen soil of different soil types

3.4 走廊帶內凍融土熱擴散系數預測方法對比及應用總結

表5和表6 分別為凍融土熱擴散系數預測模型所得結果的相關系數R2和相對誤差在10%內占比，由表可知融土熱擴散系數兩種預測方法中RBF 神經網絡在剔除樣本比例少的情況下，各土類預測精度均顯著高于二元擬合方法，其中粉土融土預測效果最好，相關系數R2和相對誤差10%內占比分別為0.90 和82%。同時，對比凍土熱擴散系數預測模型可發現，三元模型預測效果明顯優于二元模型，進一步說明融土熱擴散系數包含有重要的土壤成分信息，可有效彌補不完全參數回歸模型的缺陷。此外，三元RBF 神經網絡方法要優于三元擬合方法，說明RBF 神經網絡方法可高效捕捉熱擴散系數與影響因素間關系特征，進而有效提升預測精度，其中全風化巖類預測效果最好，相關系數R2和相對誤差10%內占比分別為0.96和90.77%。

表5 熱擴散系數預測模型相關系數R2Table 5 Correlation coefficient R2 of thermal diffusivity prediction models

表6 熱擴散系數預測模型相對誤差在±10%內占比Table 6 The relative error of the thermal diffusivity prediction models is within±10%

依據前述分析可知，青藏工程走廊帶各土類凍融土熱擴散系數分布離散，不同預測方法應用總結如下：

（1）二元擬合方法針對于熱擴散系數主體分布區間，其應用范圍較窄且誤差較大，但獲取形式簡單（只需干密度和含水率參數），可滿足一般工程估算的相應要求；

（2）三元擬合方法局限于凍土熱擴散系數，且不可避免高成本和冗長的實驗測試過程，但具有適用樣本分布范圍廣、預測性能好及應用簡單的優點；

（3）RBF 神經網絡需以大量實測樣本數據和數學算法為基礎，其開放性不強且應用相對復雜，但具有應用范圍廣、預測精度高的優勢，總體來看該方法為熱擴散系數最佳預測方法。

4 結論

本文針對青藏工程走廊典型土類的凍融土熱擴散系數進行了計算，分析了熱擴散系數分布特征和參數影響規律，提出了不同預測模型并對比了預測效果，得出主要結論如下：

（1）走廊帶土類以黏性土、砂土和碎石土居多，熱擴散系數與土類粒徑整體呈正相關性，具體為融土熱擴散系數按黏性土、粉土、全風化巖類、砂土及碎石土依次增大，凍土熱擴散系數按黏性土、全風化巖類、粉土、碎石土及砂土依次增大。

（2）熱擴散系數與干密度及天然含水率相關性隨土類及凍融狀態差異明顯，凍、融土熱擴散系數呈顯著正線性關系。

（3）融土熱擴散系數預測模型中，二元RBF 神經網絡預測效果顯著優于二元經驗公式；凍土熱擴散系數預測模型中，基于融土熱擴散系數的三元預測模型的預測精度明顯高于二元模型，且三元RBF神經網絡模型預測效果更佳。

（4）綜合熱擴散系數不同預測模型預測效果和誤差分析可得，二元擬合方法形式簡單，可滿足一般工程估算需求；三元擬合方法需以融土實驗結果為基礎，適用樣本分布范圍廣且預測精度高；RBF神經網絡模型具有最優的預測精度且應用范圍最廣，為最佳預測模型。