0 引言

冻土通常是指处于0 ℃或0 ℃以下且含有冰的土体和岩体，冻结状态持续2 a及以上的称为多年冻土^[1]。我国多年冻土面积约占国土面积的22.4%，主要分布在青藏高原和东北高纬度地区^[2]。青藏高原由于其海拔高、气候严寒的特点决定了多年冻土的存在和广泛的分布，占我国多年冻土区总面积的70%以上^[3]。214国道，起点为青海省西宁，终点为云南景洪，纵跨青海、西藏、云南，紧邻四川，青海境内长1 084 km。 由于其位于青藏高原东部边缘地带，多年冻土发育，国道214线多年冻土年平均地温高于-1.5 ℃的地区，占多年冻土地区路段的70%～80%^[4]，多年冻土是地气热交换的产物，气候的变化显著地影响着多年冻土的发展过程。多年冻土在全球气候变暖和公路工程作用下正处于退化状态，其工程地质特性十分脆弱^[5-6]。地温场变化是影响冻土的最直接因素，由于地温监测点的布设有限，因此依据现有监测资料，建立多年冻土地温预测模型，对国道214沿线任意位置的地温情况进行预测，将对揭示冻土动态发展变化规律起到积极的推动作用，同时为冻土工程建设提供基础资料。

1 国道214线多年冻土地温分区

为研究国道214沿线(鄂拉山至清水河段)地温空间分布情况，先在道路沿线埋设地温监测设备，与前期设立的监测场地构成214沿线具有连续性的完整监测网络，采集地温监测数据，绘制监测点地温分布剖面图，以便为之后的地温分区做好数据储备。

本研究在冻土发育明显的3个典型路段选取3个路段进行局部地温观测线路布设，分别是花石峡段、龙根查玛依、清水河路段。然后结合不同区段冻土类型、公路病害、土质情况对国道214线进行综合地温分区，明确国道214沿线地温分布情况^[7-9]。

G214沿线冻土按照冻土地温，段落特征(冻土类型、病害类型、土质类型)相似原则，划分为以下7个区域：

Ⅰ区：地温常年在0 ℃左右；基本无冻土分布；位于野马滩段；桩号在K520+800—K532+900之间，该区基本无严重可见病害，段落总长12.1 km。

Ⅱ区：冻土地温常年在0 ℃左右；该区多见季节性冻土；位于黄河沿南段和棉纱岭段北段；桩号在K467+850—K513+600之间，该区道路病害以局部沉陷，网裂最为常见；该段上部土质以砂土为主，下部多见风化砂岩，段落总长45.75 km。

Ⅲ区：冻土地温常年在-0.5 ℃与0 ℃之间；冻土类型以少冰、多冰冻土为主；桩号分布，K318+900—K321+420、K321+420—K327+600、K327+600—K347+300(鄂拉山南半段)、K387+700—K411+500(红土坡段、花石峡北段)、K532+900—K539+000(小野马岭中部和北部)、K578+600—K603+000(查拉坪、查拉根、巴颜喀拉山中部和北部)、K630+200—K632+800、K632+800—K650+420、K650+420—K654+870(查隆穷南段，清水河北段)；该区公路主要病害类型有破碎、横缝、局部沉陷、波浪；该区土质以亚砂土夹碎砾石、碎石土以及砂砾较常见；段落总长68.97 km。

Ⅳ区：冻土地温常年基本保持在-1 ℃；冻土类型以富冰、饱冰为主、含土冰层为辅；桩号分布，K300+960—K318+900(鄂拉山北段)、K347+300—K348+600、K348+600—K369+000(姜路岭段)、K369+000—K372+000、K372+000—K387+700(醉马滩段)、K457+000—K461+500(多次安克郎南端)、K567+000—K578+600(龙根查依玛南端)；道路病害多为网裂、沉陷、断裂板、破碎；土质多为亚砂土和砾石；段落总长64.44 km。

Ⅴ区：冻土地温常年基本保持在-0.5 ℃；冻土类型为多冰、富冰、饱冰冻土；桩号分布，K411+500—K433+000(花石峡南段，长石头山段)、K461+500—K467+850(黄河沿北端)、K539+000—K557+000、K557+000—K567+000(小野马岭北端，野牛沟段，龙根靠段)；该段公路病害以纵裂、沉陷、波浪、断裂板为主；上部土质多见砂、砾砂，碎石土、块石土，下部基本以风化岩石为主；段落总长55.85 km。

Ⅵ区：冻土地温在-0.5~1.0 ℃；冻土类型富冰、饱冰、含土冰层；桩号分布，K433+000—K457+000(多次安克郎中部和北部)、K603+000—K615+000、K615+000—K630+200(巴颜喀拉山南端，查隆穷中部和北部)；该段公路病害以纵裂、沉陷、破碎为主；土质多见亚砂土夹砾石，碎石土，块石土；段落总长51.2 km。

Ⅶ区：冻土地温基本保持在0 ℃，冻土类型富冰、饱冰、含土冰层；桩号分布，K513+600—K520+800(棉纱岭南端，野马滩北端)、K654+870—K662+900、K662+900—K664+040、K664+040—K669+330(清水河南段)；该段公路病害严重，以沉陷为主，局部冻胀翻浆；上部土质多见淤泥质土、亚黏土、碎石土，下部为风化泥岩；路段总长21.76 km。

根据国道214沿线不同路段的地温监测数据，以冻土地温为评价目标，以平均地温为主要特征参数，辅以冻土类型、土质类型、公路病害等间接反映冻土地温的特征参数，根据各个特征参数里的不同指标对冻土地温表征程度的差异，分别对国道214线年平均地温、冻土类型、土质类型以及公路常见病害4个特征参数里的指标进行量化，再以里程桩号为横坐标，以4组量化值为纵坐标分别绘制折线图，组合4组折线图，将特征相近的路段合并，最终得到7个国道214线地温分区，并绘制了国道214线地温分布直观图，见图1，为后续的地温预测提供基础数据。

2 地温预测理论模型的建立

模型假设条件：(1)地表无短暂的非气候干扰，如产热，地形，地下水流等；(2)无系统扰动，如热性能变化等^[10-12]。

则一年中任意时刻，在深度z处的浅层地温T(z)受地温梯度R(z)和地表时变温度引起的扰动θ(z)影响，用公式表示为：

(1)

式中，T₀为地表平均温度；q₀为准稳定状态下地表热流密度。

多年冻土区基本可以近似看作是一个半无限大介质，其热量传递过程可类似于一个半无限大物体在周期性边界条件作用下的非稳态导热过程。根据热传导理论，对于均质半无限大物体周期性变化边界条件下的温度场，可用一维导热微分方程近似表示为：

(2)

方程边界条件为z=0处的地表温度。太阳辐射具有周期特性，导致地表温度具有周期特性，地表温度变化也可以近似为一个简谐波。假设地表温度变化简谐波的振幅为A₀，谐振角频率为f，初始相位为φ₀，将冻土地面任意时刻的温度与平均温度做差，即：

(3)

应用分离变量法求解式(3)、(4)，得到地表温度谐振信号向地下传播引起的任意深度任一时刻的岩土温度变化：

(4)

(5)

(6)

式中，θ(z,t) 为地面下深度Z处、t时刻岩土的温度与地表平均温度T₀的差值；k为冻土综合参数；f为简谐振动频率；α为冻土热扩散系数；λ为冻土导热系数；c为冻土比热容；ρ为岩土密度；A₀为地表温度波动的振幅；t₀为地表温度波的初始相位出现的时间，单位与简谐波周期一致。

对于每个温度简谐波，振幅A随深度的衰减和相位φ随深度的偏移分别为：

(7)

(8)

太阳辐射的变化有年和日两个基本周期的变化，也有长周期的变化，地表温度变化是多个简谐周期的叠加，根据热传导理论，地表边界处同时有不同周期 τ_i(i=1,2,3,…,n) 的多个温度波动，将这些波动叠加有：

(9)

(10)

不考虑大地热流，周期性变化的边界条件下任意深度z处冻土的温度随时间的变化有：

(11)

考虑大地热流，有：

(12)

由此可见，地层内任意深度处，冻土温度变化与表面的温度变化规律相似，都可以描述为周期相同的余弦函数规律。且简谐波的振幅随着深度的增大而衰减，当温度振幅衰减到一定程度时，温度几乎不再波动，也就是不同冻土地区地温变化规律分析中[b,+∞]深度内冻土的稳定温度。

根据各个地温分区地表平均温度，年地表温度波幅和岩土热物理性能参数，可以计算处每一时刻，各个深度下的岩土温度值。

3 地温拟合模型

根据热传导理论，结合国道214线冻土地温随深度呈两阶段变化特点，将地温曲线沿深度方向划分为地温波动区与地温稳定区两大特征区域。应用ORIGIN软件，对地温波动区采用非线性拟合算法进行曲线拟合，地温稳定区选用常函数进行拟合，分别得到7个地温分区的两阶段地温曲线拟合模型。

3.1 周期的确定

由于大部分钻孔测量时间均不到1 h，测量时间内温度变化不大，1 d的温度变化影响深度较浅，根据公式计算和实测数据验证，实测数据能够代表当日不同深度处地温分布，忽略1 d地表温度变化对测量数据的影响。

对浅层地表地温分布影响的主要因素是太阳辐射作用。太阳辐射周期性变化，主要包括日夜交替变化，季节交替变化，22 a太阳黑子波动，以及其他一些周期性的变化，但其中1 a的太阳辐射变化对浅层地表温度影响比较显著，地表以下可能达到5 ℃ 左右的地温波动；另外考虑浅层地温能年度内的平衡，最终确定浅层地表温度分布的主要周期为1 a，即f=1/365 d。

3.2 拟合模型建立和参数确定

在只考虑年周期影响的情况下，地温的波动方程为：

(13)

对于比较深的孔，需要考虑大地热流造成的温度变化影响，假定冻土为均一介质的条件下，可以认为大地热流引起的温度梯度为一常数，综合考虑，在式(3)中加入线性项，进行地温修正有：

(14)

因此，只要测得某一点位某一时刻冻土地温随深度变化的曲线，利用式(13)、(14)对曲线进行拟合，就可以得到该点位的地表平均温度T₀、地表温度简谐波振幅A₀、地温简谐波振动频率f、冻土综合参数k，以及大地热流梯度β。

观察地温变化规律分析时，国道214沿线地温剖面曲线，所有曲线的共性在于，从[0,b]深度内，曲线形状近似为一个简谐波，在[b,+∞]深度，曲线近似为一条垂直于x轴的直线，该直线与X轴的交点为该分区的地温稳定值，记为C_i，地温分区不同，直线与X轴的交点不同。

将曲线做分段处理，结合式(12)，建立分段拟合曲线模型：

(15)

式中，i=1,2,3,4,5，6,7； T_i(z,t) 为i区某一深度，某一时刻地温值；z为深度；t为地温曲线测试的时间，例如在9月2号测量的地温曲线，9月2号为一年中的第247 d，因此t=247；T₀ⁱ为i区地表平均温度；A₀ⁱ为i区地表温度简谐波振幅；f_i为i区地温简谐波振动频率；t₀ⁱ为i区地表温度波的初始相位出现的时间；C_i为i区稳定地温值。

下述所有地温剖面曲线的浅层(z∈［0,b))曲线拟合均采用OriginPro8.0中的非线性拟合工具完成；对于 z∈［b,+∞ )，地温值为该区稳定地温值。

4 基于地温分区的G214线地温曲线拟合

基于式(15)，结合地温区间划分，应用OriginPro8.0中的非线性拟合工具，分别对不同地温分区的地温曲线进行数据拟合。本文只针对Ⅰ区和Ⅲ区给出拟合过程。

4.1 Ⅰ区地温曲线拟合

选取具有代表性的I 区测温点位，由于该区无冻土，地温监测点位布设较少，仅在K522+245桩号处布设了一处点位，该点位测量时间为2011年11月20号上午11:15。该区b=5 m，长年稳定地温在0 ℃,见图2。

考虑拟合精度，在监测数据的基础上，用内插法将数据点加密，加密后的数据点如表2所示。

该测点地温拟合公式：

(16)

该区无冻土，稳定地温0 ℃，C_Ⅳ=0，冻土地温从5 m深处开始稳定，b=5。

按照地温分区结果，I 区地温剖面曲线类似，Ⅰ区有着相似的温度简谐波形式，可由该区某一测点浅层地温曲线模拟，与分段地温曲线拟合模型，得到 I 区地温拟合公式(拟合参数见表3)：

(17)

当z=5时，无论t取何值，计算得到T_I(z,t)在取整后均在0 ℃左右，曲线在间断点处连续性较好。

4.2 Ⅲ区地温预测

查表6，Ⅲ区地温年增长率是0.14 ℃/a，一年按365 d计算，平均每天地温增长0.14/365=3.836×10^-4 ℃/d；参照Ⅰ区地温曲线拟合方法，得到Ⅲ区地温曲线拟合公式：

(18)

地温升高，将曲线向地温增大方向平移，则Ⅲ区N年份第t d的地温曲线公式为：

(19)

式中，N为年份；t为第N年的第t天；b^*为地温稳定深度； T_Ⅲ(b^*,t) 为 Ⅲ 区稳定地温深度以下N年份第t天的地温。

根据热传导理论，结合国道214线冻土地温随深度呈两阶段变化特点，将地温曲线沿深度方向划分为地温波动区与地温稳定区两大特征区域，应用ORIGIN软件，对地温波动区采用非线性拟合算法进行曲线拟合，地温稳定区选用常函数进行拟合，分别得到7个地温分区的两阶段地温曲线拟合模型。

5 案例分析

在国道214桩号K535+800处，钻取地质勘查孔，成孔时间为2011年11月28日上午11点10分，位置路左15 m处，钻孔编号BZK14，该孔勘探结果主要为少冰、多冰冻土，地基土主要为角砾土、细砂、块石土，多年冻土人为上限约为5 m，预测该地地温随深度的分布情况。

预测步骤：

Step1：查地温分区表4。

对应该点位桩号，位于地温分区第Ⅲ区。

Step 2：由勘查资料，可知预测年份N=2011；查各区地温曲线拟合年份表5。

可知该地温曲线拟合年份也为2011年9月2日，即N₀=2011；查表6。

可知该点位地温分区地温年均增长率v_Ⅲ=0.14，换算成每天的地温增长率${v_{{\rm{Ⅲ}}}} = {{0.14} \over {365}} = 3.836 \times {10^{-4}}^\circ {\rm{C/}}d$；11月28日是2011年的第362 d，9月2日是2011年的第275 d，则t=362-275=87；根据第Ⅲ区地温曲线拟合公式，代入公式中的浅层地温预测模型:

(20)

得该点位2011年11月28日的浅层地温曲线公式：

(21)

Step 3：绘出全深度曲线，见图3。

Step 4：根据浅层地温预测曲线读出稳定地温深度。

观察预测曲线，该点位2011年11月28号稳定地温深度b^*=5 m。

Step 5：计算T_Ⅲ(b^*,t)，作为z=5 m 深度以下的地温：

(22)

则桩号K535+800断面5 m以下稳定地温为-0.270 ℃。

Step 6：写出K535+800桩号处2011年11月28日地温曲线符合公式：

(23)

曲线绘制如图4所示。

预测验证：

为验证该桩号2011年11月28号预测的全地温曲线的结果是否准确，特在同一天实际钻孔监测地温随深度的变化情况。监测数据如表7所示。

对比K535+800断面地温预测和实测曲线，如图5所示。

二者大致形状相似，地温随深度变化的趋势是一致的，但在浅层预测地温比实际地温偏高一些，尤其在1~2 m之间，预测地温比实际地温偏高，稳定地温预测没有太大变化，稳定地温深度基本一致，用分段函数进行地温断面曲线预测，总体趋势预测良好。

6 结论

(1) 本文根据国道214线不同地点的地温曲线分布形态、年平均地温，辅以所在区域的冻土类型、土质类型以及公路常见病害等特征，将国道214线进行路段划分，最终得到国道214沿线地温7个分布区，依据该7个地温分区，绘制了国道214沿线地温分布直观图。

(2) 根据不同区域冻土地温随深度的变化规律，将地温曲线沿深度方向划分为地温波动区与地温稳定区两大特征区域，应用拟合软件，对于地温波动区，采用非线性拟合算法进行曲线拟合，地温稳定区选用常函数进行拟合，最终分别得到7个地温分区的两阶段地温曲线拟合模型。

(3) 根据地温随时间的变化规律，分别计算7个地温分区的年平均地温增长率，结合地温曲线拟合模型，按照平移理论，得到7个地温分区的地温预测模型及参数；根据地温预测模型，可对国道214沿线地温情况进行预测。

(4) 地温预测案例表明，由地温预测模型计算该测点某一时间点的地温作为该测点的理论地温，绘制出理论地温剖面曲线，将理论地温剖面曲线与该地区实测地温剖面曲线相比较，发现二者大致形状相似，地温随深度变化趋势基本一致，浅层地温理论值与实测值稍有不同，但差别不大，而稳定地温理论值与实测值相同，稳定地温深度也基本一致，用该地温预测分段函数模型进行的地温断面曲线预测，总体趋势预测良好。

(5) 在对国道214线进行地温预测的过程中，由于地温监测数据点比较稀疏，需要先对深度用内插法进行点的加密，建议今后在地温观测过程中，沿深度方向多布设一些测点，这样收集的地温沿深度方向的数据将更加准确，则对地温曲线的拟合以及之后的地温预测理论模型的建立也会更加准确。