淺層地溫剖面曲線擬合法估算熱擴散系數

何報寅，劉 杰，徐貴來，丁 超，劉紅衛，楊 柳

（1.中國科學院 測量與地球物理研究所，武漢 430077；2.湖北省地質環境總站，武漢，430034；3.武漢地質工程勘察院，武漢 430051）

淺層巖土的熱物性參數，特別是熱擴散系數（thermal diffusivity，國內文獻也叫導溫系數，或熱擴散率）是淺層地溫能計算和地源熱泵系統設計所必需的重要參數.在現有文獻報道中，巖土熱物性參數的獲得可歸結為3種，即實驗室測量法、現場測試法和文獻參考法［1－4］.實驗室測量法就是采集樣品在實驗室內用巖土熱物性參數測試儀測定.現場測試法是用野外熱物性參數儀在現場直接測定.文獻參考法是根據巖土類型及其含水量等情況，通過查找規范、手冊或相關文獻，直接引用相同或類似巖土的經驗值.

本文從地表淺層巖土熱傳導和熱擴散理論出發，提出一種新的估算方法.即利用理論模型擬合地溫實測剖面曲線求取淺層巖土的熱擴散系數.利用該方法也可同時估算出地表年平均溫度、溫度變幅、初相等大地溫度特性參數，這些參數也是地源熱泵系統設計所必需的重要參數.

1 巖土溫度場理論模型

在不考慮地表短暫的非氣候干擾或系統的擾動等違反一維傳熱假設的因素（如產熱、地形、地下水流、熱性能變化等）條件下，任一時刻在深度z處的淺層地溫T（z）由地溫梯度R（z）和地表時變溫度引起的擾動θ（z）所給定［5］：

其中，T0為地表面平均溫度（℃）；q0為準穩定狀態下地表熱流密度.

巖土通常可以近似看作是一個半無限大介質，其熱量傳遞過程相當于半無限大物體在周期性邊界條件作用下的非穩態導熱過程.根據傳熱學理論，對于勻質的半無限大物體周期性變化邊界條件下的溫度場，可用一維導熱微分方程近似描述為［3，5－6］：

該方程的邊界條件就是z＝0處的地表溫度.隨著太陽輻射的周期性變化，地表溫度變化也可以看作一個簡諧波.對于一個振幅為A0，諧振角頻率為f，初始相位為φ0的地表溫度變化的諧振成分，可表達為余弦函數形式：

式（3）中，θ（0，t）表示半無限大物體在表面z＝0處任何時刻的溫度與平均溫度的差，即θ（0，t）＝T－T0；

應用分離變量法求解（2）和（3）式，得到地表溫度諧振信號向地下傳播引起的任意深度任一時刻的巖土溫度變化［7－8］：

這里，θ（z，t）為地面下深度z（m）處、t時刻巖土的溫度與地表平均溫度T0的差值（℃）；k為巖土綜合參數（m－1），f為簡諧波動頻率；而α為巖土熱擴散系數（也稱導溫系數，單位是m2·d－1），α＝λ／cρ；λ為巖土導熱系數（W·m－1·℃－1）；c為巖土比熱（J·kg－1·℃－1）；ρ為巖土密度（kg·m－3）；A0為地表溫度波動的振幅（℃）；t0是表面溫度波的初始相位出現的時間，單位與諧波周期一致.

于是，對每個溫度簡諧波，其振幅A隨深度的衰減和相位φ隨深度的偏移分別為：

由于太陽輻射的變化有年和日兩個基本的周期變化，也有長周期的變化（如11年左右的太陽黑子周期等），所以地表溫度變化也是多個簡諧周期的疊加.根據熱傳導理論，地表邊界處同時有不同周期τi（i＝1，2，3，…，n.τi＝1／fi）的多個溫度波動，其對溫度場的作用等價于單一波動的疊加：

于是，在不考慮大地熱流的情況下，周期性變化的邊界條件下任意深度z處巖土的溫度隨時間的變化可用下式描述：

而在考慮大地熱流的情況下，可以描述為［5］：

由式（10）和（11）式可以看出，在地層內任意深度的z處，巖土的溫度的變化與表面z＝0處的溫度變化規律相似，都是周期相同的余弦函數規律.但深度為z處的溫度簡諧波的振幅隨著z的增大，逐漸衰減.當溫度振幅衰減到一定程度時，溫度波動很小，幾乎為0，可以忽略不計.在此深度以下，巖土溫度可以看作是常年保持不變，稱為等溫層.根據各城市地表面平均溫度、年地表面溫度波幅和當地巖土熱物性能參數，可以計算出每一時刻、各個深度下的巖土溫度值.

2 地溫剖面擬合模型

2.1 周期確定

大部分鉆孔測量時間都比較短，不到1h，測量時間內溫度變化不大；另外1天的溫度變化影響的深度較淺，根據公式計算和實測數據驗證，1m深處的溫度可以認為基本上不受地表溫度影響.因此實測數據可以代表當日不同深度處地溫的分布，并且可以忽略1日地表溫度變化對測量數據的影響.

對淺層地表地溫分布影響的主要因子是太陽對地表輻射，太陽輻射存在日夜變化、季節性變化、22年波動和其他周期性變化.其中1年的太陽輻射變化對淺層地表影響比較顯著，地表以下可以達到5℃的溫度波動；另外考慮到淺層地溫能的年內平衡，確定淺層地表溫度分布的主要周期為1年，即f＝1／365d.

2.2 擬合模型和參數求取

在只考慮年周期影響的情況下，即f＝1／365d，地溫的波動方程為：

但對于比較深的孔（孔深大于10m），大地熱流造成的溫度變化已經大于0.06℃，可能影響觀測數據，并且隨著深度增加，影響將進一步加大.另外長周期的太陽輻射溫度對10m以下的地溫波動也有貢獻，但振幅很小.在假定巖土為均一介質的條件下，可以認為大地熱流引起的地溫梯度為一常數.綜合考慮，可在（12）式中加入一個線性項，進行修正，修正式為：

因此，只要測得某一時刻的巖土剖面溫度隨深度變化的曲線，利用（12）或（13）式對曲線進行逆擬合，就可以得到T0、A0、k和β等參數，進而通過（5）式就可以求得巖土熱擴散系數.下文實例中，曲線擬合采用OriginPro7.5非線性擬合工具完成［9］.需要特別指出的是，由于上述模型是在勻質半無限大傳熱介質的假設條件下得到的，因此這里求得巖土熱擴散系數是代表所觀測深度內整個巖土層的平均熱擴散系數.對于巖性差異較大的多層巖土，可以根據地層情況，分段擬合，分層計算各巖層的參數.

2.3 野外地溫剖面數據測量

本文野外溫度剖面測量采用的儀器是TP2000溫度記錄儀，是北京安伏電子技術有限公司的產品.該儀器的主要性能指標見表1.

表1 TP2000溫度記錄儀主要性能指標［10］Tab.1 Main technical parameters of temperature recorder TP2000

測量方法：在野外測量時，把TP2000溫度記錄儀探頭固定在皮尺頂端，垂直吊下井孔，把探頭完全浸入水中.溫度記錄儀將自動記錄井孔中地下水的溫度，測量深度可由皮尺讀得，時間由秒表讀出.測量時，儀器的記錄時間間隔設為10s，每個深度測量時間1～2min，這樣每個深度處就可以獲得6～12個溫度讀數，數據整理時取平均值作為該深度處的地下水溫度.測量深度間隔一般設為10cm，在溫度變化緩慢的深部，間隔適當加寬至20～30cm.

3 實例計算結果和討論

3.1 實例1——洪湖市小港站

該組鉆孔由一眼承壓水孔和一眼潛水孔組成，位于湖北省洪湖市小港監測站內，經緯度為29°55′22.07"N／113°29′31.72"E，1998年 7 月 成 孔，設有套管和過濾管，井管內徑為110mm，屬地下水長期觀測孔.該孔揭示的地層為江漢平原沖洪積相，地層結構呈典型的二元結構［11］.上部0～11.5m為全新統灰綠色～灰黑色淤泥質粘土，中間11.5～17.5m為中更新統灰黑色粘土，下部17.5～27.5m為中更新灰綠～灰色粉砂層.該組孔地下水溫度測量時間為2009年8月11日，即一年中的第223d（t＝223），測量時承壓水位埋深1.05m，潛水位埋深1.05m.

曲線擬合采用（11）式進行.擬合效果如圖1.

圖1 洪湖小港觀測站承壓水孔（a）和潛水孔（b）地溫剖面曲線和擬合結果Fig.1 Temperature profile curves and fitting results for the boreholes（（a）confined water and（b）unconfined water）at Xiaogang Monitoring Station in Honghu

擬合得到的承壓水孔地溫公式為：

擬合得到的潛水孔地溫公式為：

擬合得到的參數如表2.

表2 小港站觀測孔地溫剖面曲線擬合所得參數Tab.2 Parameters obtained from fitting of temperature profile curves of boreholes at Xiaogang Station

3.2 實例2——武漢市武昌區洪山廣場

該鉆孔為工程地質勘察孔，成孔時間為2009年4月20日，位置為30°32′38.62"N，114°20′1.58"E，孔徑為108mm，不設套管，孔口絕對高程為22.5m，該孔地下水屬承壓水.地層結構為，表層0～1m為雜填土，地面為薄層水泥路面；上部埋深1～10m為黃褐、灰褐色粘土和粉質粘土，軟塑～可塑狀；中部10～43m為灰色粉砂、粉細砂和細中砂，局部夾粉土；下部43～55m為古近系（E）～白堊系（K）強風化砂礫巖和中風化泥質砂巖，泥質膠結.該鉆孔地下水溫度測量時間為2009年7月12日，測量時該孔承壓水水位埋深為1.20m.

考慮到該孔測量深度遠超出10m，需要考慮大地熱流的影響，故采用（12）式進行擬合.擬合效果見圖2.

得到該測點地溫擬合公式：

圖2 武漢市洪山廣場地溫剖面和擬合結果Fig.2 Fitting results and temperature profile curves for the boreholes at Hongshan Square，Wahan

擬合得到的參數如表3所示.

3.3 討論

從這幾例看，所用擬合方程都能很好地對淺層巖土地溫剖面實測數據進行擬合，擬合相關系數可以達到0.99以上.在江漢平原洪湖市小港站，擬合得到的熱擴散系數為0.047～0.050m2·d－1，而有關文獻［1］［2］給出的潮濕的重土（粘土、沙土）的熱擴散系數為0.055m2·d－1，潮濕的輕土（松散的沙子，淤泥）為0.045m2·d－1，二者十分接近.洪湖市小港站是江漢平原地勢最低洼之處，地下水位高，土壤常年多保持為潮濕狀態，該處的淺部地層主要為粘土和粉質粘土，相對較為均一，地面為農業耕作自然狀態.可見，擬合得到的熱擴散系數值的大小與該處巖土類型和含水狀態是相符合的.當然，參考文獻中的數據只能對擬合結果進行大致的檢驗，下一步還需要結合現場測試和實驗室樣品測試進一步檢驗.在武漢市洪山廣場測孔，擬合得到的熱擴散系數為0.090m2·d－1，大于參考文獻中潮濕粘土和沙的值，比砂礫石層的熱擴散系數（約0.085m2·d－1）略大，而接近與巖石和混凝土的熱擴散系數值.總體上比該處巖土層理論上的值偏大.這可能有兩個因素的影響，一是本文擬合得到的熱擴散系數值是反映測量深度內整個巖土層宏觀的熱擴散特性，與某單一巖土的值有一定差異，二是該處為城市中心，附近地面大部分被建筑物、水泥路面和瀝青路面所覆蓋，這些材料本身的熱擴散系數就大，同時它們分布在表面，滲透性能很差，削弱了熱傳質作用，使整個剖面的熱擴散系數偏大.該處擬合得到的線性項修正系數β為－0.077，該系數很小，對整個擬合方程的貢獻也很小.武漢市所在的江漢平原鉆孔地溫梯度約為0.028～0.050［12］，可見，大地熱流對修正系數的貢獻為一半左右，另一半的貢獻可能來源于長期氣候變化和其他因素，如全球變暖趨勢［13－14］以及城市熱島效應［15］.

表3 武漢市洪山廣場測孔地溫剖面曲線擬合所得參數Tab.3 Parameters obtained from fitting of temperature profile curves of boreholes at Hongshan square in Wuhan.

擬合得到的其它大地溫度參數也與所在地區的實際情況相符.根據作者的研究結果，在洪湖和武漢地區，年最低地表溫度一般都出現在一月中旬，年地表溫度變幅在10～14℃ 之間.洪湖地區的年平均地表溫度為18.5℃，武漢地區由于存在很強的城市熱島效應，該值達到21℃ 左右［16－18］.

應該指出的是，本文的模型是建立在一個半無限均質介質的假設基礎上的.對于各層巖性差異大的多層地層，以及地下水位低，包氣帶較厚，或人工填土較厚等情況，該模型不適用.下一步，應在本文的基礎上，發展適用于上述情況的多層介質擬合模型，更好地模擬反演各地層的熱物性參數.另外，需要注意的是，在本文的方法中，是利用測量鉆孔中地下水的溫度代替巖土剖面的溫度的，因此是假定鉆孔中地下水與周圍巖土己達到熱平衡的.這個條件在老鉆孔中是基本滿足的，對于新鉆孔，一般成孔后3～4天就基本達到平衡.為保證充分達到熱平衡，對于新鉆孔，一般要求成孔1周以后才能進行溫度測量.

4 結論

本文研究結果表明，基于淺層巖土熱傳導和熱擴散理論以及一維半無限均勻介質假設推導出的淺層巖土地溫變化數學方程，可以很好地擬合地溫實測剖面曲線，通過曲線擬合能準確地估算淺層巖土整體熱擴散系數，同時可以求取地表年平均溫度、溫度變幅、初相等大地溫度特性參數.該方法建立在可靠的理論和清楚物理機制之上，簡單有效，可以在實際工作中推廣應用.需要注意的是，該方法的適用范圍是一維半無限均勻介質，故所應用的地層條件應符合或大致符合該假設.在進一步的研究中，可在本文的基礎上，發展出針對多層地層的擬合模型，以便更好的擬合實測剖面曲線，更準確地估算大地溫度特性參數和各地層的熱擴散系數.

［1］中華人民共和國建設部.地源熱泵系統工程技術規范（GB 50366—2005）［M］.北京：中國建筑工業出版社，2005.

［2］徐 偉.地源熱泵工程技術指南［M］.北京：中國建筑工業出版社，2001.

［3］陳友明，王宇航，莫志姣.土壤初始溫度模型［J］.湖南大學學報：自然科學版，2007，34（07）：27－29.

［4］胡平放，雷 飛，孫啟明，等.巖土熱物性測試影響因素的研究［J］.暖通空調，2009，39（3）：123－127.

［5］Gonz＇alez－Rouco J F，Beltrami H，Zorita E，and Stevens M B.Borehole climatology：a discussion based on contributions from climate modeling［J］.Climate of the Past，2009，5：97－127.

［6］Carslaw H S and Jaeger J C.Conduction of heat in solids［M］.Oxford Univ.Press，New York，2nd ed.，1959.

［7］Labs K.Underground building climate［J］.Solar Age，1979，4（10）：44－50.

［8］AL－Temeemi A A，Harris D J.The generation of subsurface temperature profiles for Kuwait［J］.Energy and Buildings，2001，33：837－841.

［9］方安平，葉衛平.Origin7.5科技繪圖及數據分析［M］.北京：機械工業出版社，2006.

［10］北京安伏電子技術有限公司.TP－2000－WS土壤溫度記錄 器 ［EB／OL］.http：／／www.a－ volt.com／products／tp2000ws.htm.

［11］何報寅，蔡述明.三峽工程與長江中游淺層承壓水動態［J］.長江流域資源與環境，1999，8（1）：94－99.

［12］徐 明，趙 平，朱傳慶，等.江漢盆地鉆井地溫測量和大地熱流分布［J］.地質科學，2001，45（1）：317－323.

［13］Majorowicz J，Safanda J，Harris R，Skinner W.Large ground surface temperature changes of the last three centuries inferred from borehole temperatures in the Southern Canadian Prairies，Saskatchewan［J］.Global and Planetary Change，1999，20：227－241.

［14］Henry N.Pollack，Jason E.Smerdon，and Peter E.van Keken.Variable seasonal coupling between air and ground temperatures：A simple representation in terms of subsurface thermal diffusivity［J］.Geophysical Research Letters，2005，32，doi：10.1029／2005GL023869.

［15］Yamano M，Goto S，Miyakoshi A，et al.Reconstruction of the thermal environment evolution in urban areas from underground temperature distribution［J］.Science of The Total Environment，2009，407（9）：3120－3128.

［16］陳 峰，何報寅.基于小波變換的城市熱場多尺度分析——以武漢市為例［J］.長江流域資源與環境，2007，16（Z1）：108－113.

［17］陳 峰，何報寅，龍占勇，等.利用Landsat ETM＋分析城市熱島與下墊面的空間分布關系［J］.國土資源遙感，2008，76（2）：56－61.

［18］何報寅，徐貴來，丁 超，等.基于MODIS的武漢城市圈地表溫度場特征［J］.長江流域資源與環境 ，2010，19（12）：1379－1385.