1978−2017 年西藏高原冻融指数时空变化特征

廖 莹，范继辉，李怡颖，程根伟

（1. 中国科学院、水利部成都山地灾害与环境研究所，四川 成都 610041；2. 中国科学院大学，北京 100049）

冻土是指温度在0 ℃或0 ℃以下并含冰的各种岩石或土，主要分布在对气候变化敏感的高纬度和高海拔区域[1]。2007–2016 年，全球变暖导致全球平均多年冻土温度增加了0.29 ± 0.12 ℃[2]，青藏高原的多年冻土范围自1970 年以来缩减了约10 万km2[3-4]。特别是近地表的正负积温变化将直接影响多年冻土和季节冻土的消长。冻融指数是给定时期内的温度累计值[5]，是近年来较为流行的量化冻土季节性变化的气候参数，主要包括季节冻结/融化指数[6]、近似冻结/融化指数[7]和年冻结/融化指数[8]等。其中，年冻结/融化指数是计算一年内所有低于/高于0 ℃的日均温之和[9]，该方法已广泛运用在北半球[10-11]、加拿大[12]、挪威[13]、蒙古[14]和我国的青藏高原[15-16]、东北[17-18]和北疆[19-20]等地区。由于冻融指数及其变化对地表能量交换以及生态系统过程有重大影响，因此被作为气候变化的敏感指示器[21]，在气候变化评估、冻土制图以及寒区工程环境评估等领域应用广泛[22]。例如，Shi 等[8]利用冻融指数的连续空间分布数据，结合气候预测模型[23]估算了北半球的冻土范围和面积。Peng 等[12]利用空气融化指数和环境参数模拟了北半球1971-2000 年的活动层厚度。

西藏高原作为青藏高原的主体，平均海拔4 000 m以上，是我国高寒冻土的主要分布区，多年冻土面积近70 万km2，约占整个青藏高原多年冻土面积的一半[24]。在气候变化的背景下，作为“全球气候变化敏感区”的西藏高原表现出整体变暖、变湿的特征[25]。而气温的升高对高原冻土影响强烈，加之人类活动引起的负面影响，近30 年活动层厚度呈现整体增大趋势[26-27]。但由于西藏高原的气象站点极度缺乏，部分站点建成较晚，整体存在分布稀疏、数据缺失等问题，因此气候变化背景下西藏高原冻融指数变化的研究还较为欠缺，仅有的几例研究主要聚焦在雅鲁藏布江干流或工程沿线区域[15-16,28]，这些区域靠近人类活动区，受人类活动和工程建设影响强烈，尚不能代表西藏高原整体的冻融特征。因此本文将在ArcGIS 等软件的支持下，利用西藏高原共34 个气象站点的日均气温数据，对1978-2017 年冻融指数的时空变化特征进行分析，探讨其潜在的影响因子，为下一步估算西藏高原冻土范围和面积、活动层厚度等提供数据，并为西藏高原农牧业生产、草地生态环境评估、工程建设和冻害防治等提供参考。

1 材料与方法

1.1 研究区概况

西藏高原(26°50′ - 36°53′ N，78°25′ - 99°06′ E)位于北半球中纬度地带，喜马拉雅山脉及其以北地区，青藏高原西南边陲，面积120.22 万km2，绝大部分区域(92%以上)位于海拔4 000 m 以上的高寒生态区[29](图1)。西藏高原具有独特、多样的气候特征，大致表现为西北地区严寒干燥，东南地区温暖湿润。2017 年，全区年均气温在-1.1～12.1 ℃，由东南向西北递减；干湿分明，年均降水量在117.4～929.4 mm，也呈东南向西北递减的分布规律(西藏自治区气象局http://xz.cma.gov.cn)；植被从东南向西北依次呈现森林、草甸、草原和荒漠的植被景观，与气温和降水的空间分布高度相关。

图1 西藏高原气象站点分布图Figure 1 Study areas and distribution of national meteorological stations on the Tibetan Plateau

1.2 数据来源及预处理

本研究使用的数据包括：①1978-2017 年西藏34 个气象站点日均温数据(距地面2 m 处气温)来源于中国地面气候资料日值数据集(http://data.cma.cn)，由于建站较晚，南木林、贡嘎、墨竹工卡、类乌齐、浪卡子、洛隆和八宿7 站点仅有1992-2017 年的数据；②数字高程模型(digital elevation model, DEM)数据(分辨率90 m × 90 m)和植被归一化指数(normalized difference vegetation index, NDVI)数据(2008-2017)来源于中国科学院资源环境科学数据中心(http://www.resdc.cn)；③人口、地区生产总值(gross domestic product, GDP)和第一产业增加值数据为34 个气象站点所在县(市) 2013 - 2016 年的多年均值，来源于《中国县域统计年鉴(县市卷-2014)》、《中国县域统计年鉴(县市卷-2015)》、《中国县域统计年鉴(县市卷-2016)》和《中国县域统计年鉴(县市卷-2017)》[国家数字电子图书馆(http://www.nlc.cn)]。

考虑到气象数据中存在缺测、漏测等情况，利用缺测站点与周边站点间的数据，按照空间相关性插值，对缺测部分进行插补。对不同的缺测情况采用不同的插补方法，具体为：①单日缺测：取缺测日前后1 d 平均值；②连续2 日缺测：第1 缺测日取其前2 d 平均值，第2 缺测日取其后2 d 平均值；③缺测连续2 日但不超1 年：取前后各1 年同时段的平均值；④缺测日数超1 年：利用有数据年份的数据与周边站点建立相关关系，采用最相关的2～3 个站点数据进行插值。

1.3 研究方法

1.3.1 冻融指数计算

冻融指数包括冻结指数和融化指数，是冻结/融化期内所有气温小于/大于0 ℃的累计值。考虑到西藏高原冷暖季分明，为保证计算的冷季连续，本研究以每年7 月1 日至翌年6 月30 日为周期计算冻结指数，以一个日历年的1 月1 日至12 月31 日为周期计算融化指数[10]。其计算公式分别为：式中：FI(freezing index)、TI(thawing index)分别为冻结指数和融化指数(℃·d)；Ti和Tj分别为负值逐日温度(℃)和正值逐日温度(℃)；NF和NT分别为年内温度低于和高于0 ℃天数。

1.3.2 趋势分析

采用一元线性回归分析计算34 个气象站点气温、冻结指数和融化指数的变化趋势并通过最小二乘法拟合的斜率反映，当斜率 ＞ 0 时，表示增加趋势，当斜率 ＜ 0 时，表示降低趋势。斜率计算方法为[30]：

式中：SlopeX为站点x的变化斜率；i为序列序号；n为序列长度；Xi为站点x第i年的值。

1.3.3 空间插值

通过反距离插值法(inverse distance weighted)对34 个气象站点的冻融指数进行空间插值[31-33]，以获得空间连续的冻融指数均值和变化趋势。反距离插值法是对所有样本点线性加权，以估计未知点的值。权重与未知点与样本点之间的距离成反比，距离越近，权重越大。具体公式为：

本研究中1978-1991 年的冻融指数空间结果由基于拉萨、察隅等27 个站点数据空间插值结果计算得到，1992-2017 年的结果由基于增加站点后的34 个站点数据空间插值计算得到。

1.3.4 相关关系分析

采用Pearson 相关系数反映冻融指数与纬度、海拔、人口和GDP 等因子间的相关程度[34]。计算公式如下：

式中：rij为冻结(融化)指数i和某影响因子j之间的相关系数，n为样本数。

由于各因子间会相互影响，造成相关关系结果的不准确，因此本文利用偏相关分析(partial correlation analysis)，在控制其他影响因子的基础上，计算某一影响因子与冻结(融化)指数的偏相关关系。公式如下[35]：

式中：rij,k为控制某影响因子k的作用后，冻结(融化)指数i和另一影响因子j之间的偏相关系数。控制变量增加时的公式以此类推。

2 结果与分析

2.1 西藏高原气温时空变化特征

基于1978-2017 年34 个气象站点数据，计算西藏高原年平均气温，并分析其时空变化特征(图2)。图2a 表明，以34 个气象站点代表的西藏高原1978-2017 年平均气温在3.04～5.72 ℃，均值4.41 ℃，并呈波动上升趋势，气候倾向率为0.06 ℃·a-1。图2b显示，近40 年西藏高原年平均气温整体呈增加趋势，类乌齐附近区域的气温增速最快，仅有八宿附近的气温呈现下降趋势。

图2 1978−2017 年西藏高原年均气温年际变化趋势(a)和空间变化趋势(b)Figure 2 Inter-annual variation (a) and spatial variation trends (b) of air temperature on the Tibetan Plateau during 1978−2017

2.2 冻融指数年际间变化特征

西藏高原冻结指数的年际变化过程显示，近40 年西藏高原冻结指数介于386.8～886.1 ℃·d，最低值出现在2017 年(386.8 ℃·d)，最高值出现在1982 年(886.1 ℃·d)，均值为591.2 ℃·d (图3a)；冻结指数总体呈现波动下降的趋势(P＜ 0.01)，下降速率9.2 ℃·d·a-1。冻结指数的距平统计结果(图3b)表明西藏高原冻结指数年际变化十分明显，冻结指数距平值在近40 年内逐渐降低(P＜ 0.01)，21 世纪以后，西藏高原各年的冻结指数小于近40 年均值。

图3 1978−2017 年西藏高原冻结指数年际变化趋势(a)和距平变化(b)Figure 3 Inter-annual variation (a) and variations in the departures (b) of the freezing index on the Tibetan Plateau during 1978 − 2017

由西藏高原融化指数年际变化过程可以看出，近40 年西藏融化指数范围介于1 960.3～2 521.5℃·d，最低值出现在1986 年(1 960.3 ℃·d)，最高值出现在2009 年(2 521.5 ℃·d)，均值为2 209.3 ℃·d(图4a)。融化指数整体呈现波动上升的趋势(P＜0.01)，上升速率为11.1 ℃·d·a-1。融化指数距平统计图4b 表明，西藏高原融化指数整体年际差异较大，融化指数距平值在近40 年内逐渐上升(P＜ 0.01)，在1995-2002 年内存在数值的小幅震荡，1995 年以前，西藏高原各年的融化指数小于近40 年平均值，自2002 年以后则大于近40 年平均值。

图4 1978−2017 年西藏高原融化指数年际变化趋势(a)和距平变化(b)Figure 4 Inter-annual variation (a) and variations in the departures (b) of the thawing index on the Tibetan Plateau during 1978 − 2017

2.3 冻融指数空间变化特征

西藏高原34 个气象站点冻结指数与融化指数的统计特征(表1)显示安多的冻结指数多年均值最大(1 759.95 ℃·d)，察隅的冻结指数多年均值最小(0.03 ℃·d)。冻结指数均值的空间分布特征呈现自东南向西北逐渐递增的趋势(图5a)。1978-2017 年间，除八宿和察隅之外，所有站点的冻结指数均呈下降趋势(表1)，其中，类乌齐的冻结指数下降速率最快，倾向率为-17.4 ℃·d·a-1，林芝的冻结指数下降速率最慢，倾向率为-0.5 ℃·d·a-1。西藏高原冻结指数的倾向率整体呈现自东南向西北递减的空间分布特征(图5b)，表明西藏高原西北地区冻结指数下降速率较西藏高原东南地区更快。

图5 1978−2017 年西藏高原冻结指数均值(a)和变化趋势(b)Figure 5 Mean annual values (a) and spatial variation trends (b) of the freezing index (FI) on the Tibetan Plateau during 1978 − 2017

据西藏高原34 个气象站点融化指数统计特征(表1)显示，察隅的融化指数多年均值最大(4 446.66℃·d)，安多的融化指数多年均值最小(940.71 ℃·d)。融化指数均值在西藏高原东部和中南部个别区域出现高值中心(图6a)。1978-2017 年，除八宿附近区域的融化指数呈下降趋势，西藏高原整体经历了融化指数的快速上升，类乌齐和拉萨的融化指数上升最快(图6b)，倾向率分别为22.2 和18.4 ℃·d·a-1(表1)。

图6 1978−2017 年西藏高原融化指数均值(a)和变化趋势(b)Figure 6 Mean annual values (a) and spatial variation trends (b) of the thawing index (TI) on the Tibetan Plateau during 1978 − 2017

表1 西藏高原34 个气象站点冻结指数与融化指数的统计特征Table 1 Summary of mean annual freezing and thawing indices for 34 meteorological stations on the Tibetan Plateau

2.4 冻融指数的影响因子分析

采用Pearson 相关系数[36]来衡量冻融指数与自然因子(如纬度、经度、海拔)、NDVI 和社会经济因子(如人口、GDP、第一产业增加值)共7 个因子间的线性相关性强弱(表2)。表中显示，冻结指数和融化指数与自然因子的相关关系更强。冻结指数与纬度、经度、海拔和NDVI 都显著相关，与纬度和海拔呈极显著正相关关系(P＜ 0.01)，与经度和NDVI 呈显著负相关关系(P＜ 0.05)；融化指数与经度和海拔分别呈现显著正相关(P＜ 0.05)和极显著负相关关系(P＜ 0.01)，与纬度和NDVI 关系不显著(P＞ 0.05)。冻结指数和融化指数与社会因子之间并无显著的相关关系。Pearson 相关分析表明，冻结指数与海拔呈现最高正相关关系，融化指数与海拔呈现最高负相关关系。

冻融指数与7 个因子的偏相关结果(表2)表明冻结指数与纬度和海拔呈极显著正相关关系(P＜0.01)；融化指数与海拔呈极显著负相关关系(P＜0.01)，与NDVI 呈显著负相关(P＜ 0.05)。冻结指数与海拔呈现最高正相关关系，融化指数与海拔呈现最高负相关关系。

表2 西藏高原冻融指数与影响因子的Pearson 相关及偏相关系数Table 2 Pearson correlation and partial correlation coefficients of the freezing and thawing indices and influencing factors on the Tibetan Plateau

综合以上两种分析结果，认为冻融指数与海拔的相关性最强。将冻结指数和融化指数与海拔的偏相关关系以散点图的形式更直观地表达(图7)，可以观察到冻结指数与海拔有很强的正线性关系(P＜0.01) (图7a)，融化指数与海拔具有很强的负线性关系(P＜ 0.01) (图7b)，R2值均超过了0.5。

图7 西藏高原冻结指数(a)和融化指数(b)与海拔的偏相关散点图Figure 7 Scatterplot of the partial correlation between the freezing index (a) and thawing index (b) and altitude on the Tibetan Plateau

3 讨论

在全球变化背景下，青藏高原地区气温持续上升已是得到观测验证的事实[37-38]，本研究中计算的西藏高原气温倾向率0.06 ℃·a-1与SROCC 评估[39]表明的增温速率0.03 ℃·a-1和杨春艳等[40]得到的西藏近50 年气温倾向率每10 年升高0.58 ℃的结果基本一致。西藏高原1978-2017 年冻融指数及其变化趋势均符合环北极地区和基于气象站点观测数据的青藏高原冻融指数范围和倾向率(表3)，由于不同数据来源、区域范围差异和所选气象站点不同，存在具体数值上的偏差。值得注意的是，基于不同数据来源的青藏高原冻融指数及其变化趋势差异较大，可能是因为全球格网气候数据产品的空间分辨率较大(2.5° × 2.5°/0.5° × 0.5°)[8]，更易忽视下垫面的复杂程度，尤其是青藏高原地区地势起伏剧烈，低分辨率的数据产品难以模拟复杂的气温分布情况，导致计算结果的出入较大。因此，针对尺度较小的区域，利用气象站点的观测数据计算冻融指数会更接近真实值。

表3 不同数据来源和不同区域的冻融指数及其变化趋势对比Table 3 Comparison of freezing and thawing indices and the variation trends of different data sources and different regions

本研究基于各站点的冻融指数及其变化趋势进行空间插值以获得区域内冻融指数连续空间分布特征。冻结指数均值呈现由东南向西北逐渐增加的特征与气温的分布特征相关，而对气候变化的敏感度更高也是导致藏北地区冻结指数下降速率较西藏高原其他地区更高[40]的主要原因。融化指数在西藏绝大部分区域都呈上升趋势，并在拉萨、林芝和类乌齐等区域呈现出点状高值中心，可能与近年来农业开垦、城市扩张、产业发展等较强的社会经济活动有关，同时，这些区域也是西藏生态退化区[29]，对人类活动的影响更为敏感。

本研究归纳的4 个自然因子与西藏高原冻融指数都有一定的相关关系，但纳入考虑的3 个社会经济因子与冻融指数并无较显著的相关关系，说明冻融指数作为一项气候参数更易受到自然条件的影响，尤其是与气温关系紧密的自然因子。由于纬度位置决定地球上的热量分配，近地面气温也存在垂直高度上的变化规律，所以纬度和海拔与气温的关系最紧密，并与冻融指数显著相关。由于西藏高原气象站点极为稀少，仅有的34 个标准气象站也大多分布在海拔3 500～4 500 m 的位置，因此本研究利用空间插值获得的冻融指数连续空间分布特征所能代表的范围和精度有限。虽考虑到冻融指数与海拔具有极显著相关关系，但目前成熟的插值方法尚未能将各影响因素的贡献度作为插值过程中的权重，下一步研究将尝试插值方法上的创新，优化冻融指数及其变化趋势的空间结果。

土壤的冻胀、融沉及其年内循环会造成土壤力学性质的变化，严重时会发生工程冻害，影响土壤稳定性和上部建筑的安全[42]。冻融指数是绘制多年冻土和季节性冻土分布图[23]、估算土壤冻结深度[43-44]、测量冻土活动层厚度[45-46]以及指示冻融作用强度和特征的重要指标，近年来在寒区工程环境评估领域应用广泛[47]。随着西藏自治区人口增长、工农业发展、旅游开发、公路铁路网逐步完善与基础设施建设的开展，预防工程冻害是重中之重，本研究依据最新的气象资料计算分析冻融指数并获得连续的空间分布数据，可以为下一步估算西藏高原多年冻土面积、模拟多年冻土空间分布范围、测量各区域土壤冻结深度等提供重要参数数据，为西藏高原农牧业生产、草地生态环境评估、工程布局和冻害治理等提供科学参考。

4 结论

通过上述分析，可以得出以下几点主要结论：

1) 1978-2017 年西藏高原的气温在波动中上升，多年平均值为4.41 ℃，气候倾向率为0.06 ℃·a-1，约是全球同期升温率的2 倍，对全球变化高度敏感。

2) 1978-2017 年西藏高原冻结指数范围在386.8～886.1 ℃·d，平均值为591.2 ℃·d，融化指数范围在1 960.3～2 521.5 ℃·d，平均值为2 209.3 ℃·d。冻结指数和融化指数的变化趋势相反，分别随时间波动下降和波动上升，变化速率为-9.2 和11.1 ℃·d·a-1。

3) 近40 年冻结指数均值和下降速率均呈现出自东南向西北递增的空间分布规律。融化指数均值较高和上升明显的区域集中在以拉萨为主的西藏高原中南部和东部。

4) 海拔是影响西藏高原冻融指数的关键因子，与冻结指数呈现极显著正相关关系，与融化指数呈现极显著负相关关系。