浅层地温剖面曲线拟合法估算热扩散系数

何报寅，刘 杰，徐贵来，丁 超，刘红卫，杨 柳

（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.