索朗塔杰 ,杜 军 ,次旺顿珠 ,平措桑旦

（1.西藏高原大气环境科学研究所,拉萨 850001；2.西藏高原大气环境研究重点实验室,拉萨 850001；3.西藏自治区气候中心,拉萨 850001）

引言

2019 年，全球平均温度较工业化前水平高出约1.1℃，是有完整气象观测记录以来的第二暖年。1951～2019 年，亚洲陆地表面年平均气温呈显著上升趋势，速率为0.23℃/10a；我国地表年平均气温也呈显著上升趋势，增温速率为0.24℃/10a[1]。素有“世界屋脊”和“亚洲水塔”之称的青藏高原[2−3]，增温趋势更为显著，1961～2019 年地表年平均气温增温速率达0.37℃/10a[1]。在全球变暖背景下，作为冰冻圈中重要部分之一的冻土也发生了显著的变化。青藏高原多年冻土约占陆地面积的53%，是地球上分布最广泛的高海拔多年冻土[4−5]，冻土温度上升、活动层变厚和冻土面积退缩已对青藏高原气候、水资源、生态环境等产生了明显的影响[6−10]。

冻融指数作为冻土研究的重要参数之一，对冻土变化研究具有非常重要的意义，同时也是气候变化研究的重要指标[11]。近年来，学者们在青藏高原已开展了有关冻融指数的研究。姜逢清等[12]研究指出1966～2004 年青藏铁路沿线7 个气象站点的冻结指数均呈减少趋势，融化指数呈增加趋势。赵红岩等[13]发现1980～2005 年青藏铁路沿线地表冻/融指数存在缓慢的、波动的上升趋势。Wu 等[14]认为自1998 年以来，青藏高原融化指数（TI）和冻结指数（FI）分别呈现出明显的上升和下降趋势，永冻土地区FI的下降幅度比季节性冻土地区更明显；季节性冻土地区TI的上升幅度大于永冻土地区。刘磊等[15]分析得出，1977～2017 年雅鲁藏布江中下游大气冻结指数、地面冻结指数、大气融化指数、地面融化指数的气候倾向率分别为−36.6、−48.7、90.7、115.8℃·d/10a。

羌塘高原约占青藏高原总面积的20%，是世界上海拔最高、气候条件最恶劣的高原，生态环境十分脆弱。作为仅次于格陵兰国家公园的世界第二大陆地自然保护区和我国第二大自然保护区的羌塘国家级自然保护区（以下简称自然保护区）就位于羌塘高原，它也是平均海拔最高的自然保护区，这里气候寒冷而干燥，属于高原寒带季风干旱的气候，生态系统非常脆弱和敏感，生态环境安全的战略地位显得日趋重要[16]。自然保护区大部分地方为多年冻土区，季节性冻土主要分布在其南部边缘地区[17]。活动层厚度由盆地向四周增大，多年冻土区边缘是活动层厚度最大的区域[18]。然而，目前基于冻融指数对自然保护区冻土时空变化的研究未见报道，加之保护区暂无长时间序列观测资料的气象站，给相关研究带来极大困难。因此，本文利用其周边同属羌塘高原湖盆高寒草原区的5 个气象站1971～2019 年逐日平均气温、地表平均地温资料，统计分析了保护区大气和地面冻融指数特征及其时空变化，并预估在RCP4.5 和RCP8.5 两种排放情景下，未来80a（2021～2100 年）大气和地面冻融指数的演变情况，力求为青藏高原冻土研究、生态环境保护和应对气候变化提供科技支撑。

1 资料与方法

1.1 资料

自然保护区内无长时间序列观测资料的气象观测站，本文选取其周边同属羌塘高原湖盆高寒草原区[3]的安多、班戈、申扎、改则和狮泉河5个气象站（图1）1971～2019 年逐日平均气温、地表平均温度、积雪日数等资料。其中，申扎、班戈逐日地表平均温度只有1980～2019 年的资料。自然保护区冻融指数为5 个气象站的平均值；多年平均值为1981～2010 年平均。

图1 自然保护区周边气象站点分布

未来预估数据为国家气候中心发布的《中国地区气候变化预估数据集》Version 3.0 数据中的全球模式数据（CMIP5），CMIP5 数据采用国际耦合模式比较计划第5 阶段中21 个全球大气与海洋环流耦合模式数值模拟集合平均结果[19]。本文利用国家气候中心制作生成的一套包括历史模拟数据（1901～2005 年）和RCP4.5、RCP8.5 两种排放情景下预估数据（2006～2100 年）的月平均资料。

1.2 方法

1.2.1 冻融指数计算方法

冻融指数可分为大气冻融指数和地面冻融指数，分别通过逐日气温和地表温度（地温）数据计算得到。本文中冻结指数统计时间为7 月1 日至翌年6 月30 日，融化指数统计时间为1 月1 日～12 月31 日。冻融指数计算公式[11,13−15]如下：

式中：FI（Freezing indices）为大气（地面）冻结指数（℃·d）；TI（Thawing indices）为大气（地面）融化指数（℃·d）；Ti为逐日气温（地温）（℃）；NF为年内气温（地温）<0℃的日数（d）；NT为年内气温（地温）>0℃的日数（d）。

1.2.2 气候统计分析方法

本文采用线性回归和Mann-Kendall（M-K）突变检测等方法[20]，分别统计分析自然保护区年冻融指数的变化趋势、气候突变特征。

2 结果分析

2.1 冻融指数的空间分布特征

表1 给出了1981～2010 年平均冻融指数的空间分布。从表1 可知，自然保护区各站年大气冻结指数（AFI）、地面冻结指数（GFI）分别为1146.6～1801.3℃·d、930.6～1720.1℃·d，年大气融化指数（ATI）、地面融化指数（GTI）依次为728.8～1143.2℃·d、2043.7～3591.1℃·d。ATI和GTI总体上呈自西向东递减的分布，并随海拔升高而减小；AFI和GFI的分布规律不明显，但最大值均出现在安多，最小值出现地不同。

表1 1981～2010 年自然保护区冻融指数多年平均值

就自然保护区而言，AFI平均值为1388.0℃·d，最高值为1804.6℃·d，出现在1998 年；最低值出现在2018 年，仅为979.2℃·d。ATI平均值为1326.6℃·d，最高值为1542.1℃·d，发生在2016 年；最低值为1079.2℃·d，出现在1976 年。GFI平均值为1076.9℃·d，最高值在1979 年，为1507.4℃·d；最低值在2018 年，为743.0℃·d。GTI平均值为2427.8℃·d，最高、最低值分别为2880.8℃·d（2010 年）、1962.9℃·d（1979 年）。从冻融指数的大小来看，GTI>AFI>ATI>GFI。与青藏高原及周边区域比较，自然保护区的FI（AFI、GFI）低于青藏高原永久冻土区[14]，高于青藏高原季节性冻土区[14]、黄河源区[21]和雅鲁藏布江中下游[15]；TI（ATI、GTI）正好相反，比青藏高原永久冻土区偏高，较青藏高原季节性冻土区、黄河源区和雅鲁藏布江中下游偏低。

2.2 冻融指数的时间变化特征

2.2.1 年际变化

如图2a 所示，近49a（1971～2019 年）自然保护区AFI、GFI均表现为显著的减少趋势，平均每年分别减少8.97℃·d 和10.45℃·d，通过显著性检验水平（P<0.001），尤其是1991 年以来，AFI减幅更明显，为−13.85℃·d·a−1（P<0.001）。ATI、GTI均呈现出显著增加趋势（图2b），增幅分别为7.05℃·d·a−1（P<0.001）和11.38℃·d·a−1（P<0.001），以GTI增幅最为明显。地面冻融指数的变化率比大气冻融指数的变化率大，是因为近39a 地表温度的上升速度（0.70℃·10a−1，P<0.001）比气温的上升速度（0.53℃·10a−1，P<0.001）快。

图2 1971～2019 年自然保护区冻融指数年际变化

通过对比与青藏高原及周边地区年冻融指数的变化（表2），结果显示，研究时段内呈现出相同的变化趋势，即AFI、GFI减少，ATI、GTI增加。自然保护区AFI减幅明显大于雅鲁藏布江中下游地区、青藏高原季节性冻土区和黄河源区，低于整个青藏高原及其永久冻土区；GFI减幅大于雅鲁藏布江中下游地区和青藏高原季节性冻土区，但不及于整个青藏高原及其永久冻土区与黄河源区；ATI增幅小于雅鲁藏布江中下游、青藏高原季节性冻土区和黄河源区，接近于整个青藏高原；总体来看，GTI增幅要大于青藏高原及周边地区。

表2 冻融指数变化的对比

因班戈、申扎2 站逐日地表平均地温从1980 年开始有完整资料，便于比较，本文分析了近39a（1981～2019）冻融指数变化趋势的空间分布（表3）。从表3可知，自然保护区各站AFI均呈现为减少趋势，平均每年减少7.75～14.81℃·d（P<0.001），改则减幅最大；特别是近29a（1991～2019 年）减幅更明显，绝大部分站点AFI的减幅在10.0℃·d·a−1以 上，其中安多减幅达17.59℃·d·a−1。ATI在各站点上都表现为显著增加趋势，为5.38～10.66℃·d·a−1（P<0.001），狮泉河最大，申扎最小。所有站点的GFI均为显著减少趋势，为−5.11～−14.26℃·d·a−1（P<0.01），班戈减幅最大，其次是改则，为−12.46℃·d·a−1（P<0.001），申扎最小。GTI在各站点上都表现为增加趋势，增幅为8.90～23.83℃·d·a−1（P<0.001），狮泉河最大，申扎最小，两者相差2.7 倍。

表3 1981～2019 年自然保护区各站点冻融指数的变化趋势

2.2.2 年代际变化

图3 给出了1971～2019 年自然保护区冻融指数的10 年际距平变化。分析表明，FI在1970s～1990s 为正距平，进入21 世纪后转为负距平，表现为逐年代递减的变化特征。TI则呈现出逐年代递增的年代际变化特征，即在1970s～1980s 为负距平，1990s～2010s 年代为正距平，尤其是2010s的GTI较常年值偏高244.3℃·d，比1970s 偏高了407.5℃·d。在30 年际的变化上，1980～2010 年与1971～2000年比较，AFI、GFI分别偏低94.3℃·d、120.7℃·d，ATI、GTI分别偏高67.6℃·d、103.7℃·d，以地面冻融指数变幅较大。

图3 1971～2019 年自然保护区冻融指数年代际变化

2.2.3 突变分析

图4 给出了近49a 自然保护区冻融指数的M-K 突变检验。如图4a所示，AFI的UF（正序特征曲线）在1971～1983 年呈振荡上升态势，多数年份为正值，1988～2019年UF 为负值且呈明显下降趋势，UF 在2003～2019 年突破了−1.96，这表明2003 年以后AFI持续减少；UF和UB（逆序特征曲线）在2001 年出现交叉点，且在±1.96 之间，确定突变点发生在2001 年，突变后AFI偏低了277.3℃·d。同理可知，GFI、ATI、GTI分别在1999 年、1993 年和1998 年发生了气候突变，突变后较突变前分别偏低286.8 ℃·d、偏高183.5 ℃·d 和偏高327.3℃·d。以上分析表明，自然保护区冻融指数发生了明显变化，地面冻融指数时间转折点发生在20 世纪90 年代末，ATI突变时间最早，较AFI突变点偏早8 年。突变后FI降幅明显，TI升幅更大，这也佐证了气候变暖的事实。

图4 1971～2019 年自然保护区冻融指数的M-K 检验（a.AFI，b.GFI，c.ATI，d.GTI)

丁一汇等[22]研究表明，我国年平均气温在1998年发生了突变，年平均气温升高趋势从1998 年以后开始停滞，甚至下降。1991～2000 年是我国气温显著上升的10 年，该年代也是我国和全球的大部分地区气候有显著变化的时段。杜军等[23]分析认为，1971～2017年自然保护区年平均气温突变点发生在1995 年，较我国气温突变年份略早。以上分析表明，近49a 自然保护区多数冻融指数的突变年份比平均气温偏晚，这也说明冻融指数的突变确实是受气候变暖的影响。

2.3 冻融指数变化的气候成因分析

冻融指数是根据逐日平均气温（地温）计算得到，因而与温度具有直接的关系，考虑到冻土作为一种特殊的介质显著影响了地表与大气之间的能量交换[11]，冻融指数与气候因子之间的关系也值得探讨。表4给出了冻融指数与主要气候因子的相关系数。如表所示，FI与平均气温、平均最低气温和降水量呈显著的负相关，与日照时数、平均风速、积雪日数和最大积雪深度存在显著的正相关关系。而TI与上述气候要素的相关性正好相反，即与平均气温、平均最低气温和降水量呈正相关，与日照时数、平均风速、积雪日数和最大积雪深度呈负相关。分析表明，近49a 自然保护区FI（TI）减少（增加），与平均气温、平均最低气温显著升高[24]、降水量增加[24]、平均风速减小[25]密切相关，还与积雪日数和最大积雪深度的减少[23]有关。

表4 1971～2019 年自然保护区冻融指数与气候要素的相关系数

2.4 未来气候变化对冻融指数的影响

2.4.1 未来气温变化

CMIP5 试验全球模式对青藏高原气候模拟能力较为稳定[26−27]。本文利用给出的21 个模式结果，计算了未来80a（2021～2100 年）自然保护区气温变化的预估结果。表5 列出了RCP4.5 和RCP8.5 两种排放情景下，近期（2021～2040 年）、中期（2041～2070 年）、远期（2071～2100 年）3 个时期以及未来80a（2021～2100 年）自然保护区年平均气温距平，可以看出：在两种排放情景下不同时段内年平均气温均呈升高趋势，与气候基准年（1991～2020 年）相比，RCP4.5 和RCP8.5 两种排放情景下，未来80a（2021～2100 年）平均气温分别升高1.95℃和3.26℃。

表5 RCP4.5 和RCP8.5 排放情景下未来80 年（2021～2100 年）自然保护区年平均气温变化

2.4.2 未来冻融指数的变化

本文选取表5 中在RCP4.5 和RCP8.5 两种排放情景下气温增暖值来预估自然保护区冻融指数的变化情况。首先建立冻融指数与年平均气温的回归方程（表6），通过显著性检验水平（P<0.0001）；然后利用回归方程分别计算在RCP4.5、RCP8.5 两种排放情景下，近期、中期、远期3 个时期以及未来80a 自然保护区冻融指数的预估值。

表6 羌塘自然保护区冻融指数与年平均气温的回归方程

与1991～2020 年比较（表7），在RCP4.5、RCP8.5两种排放情景下，近期、中期、远期自然保护区均表现出FI减少、TI增加的变化特征，以远期最为明显；在RCP8.5 排放情景下，AFI、GFI分别减少829.4℃·d、918.6℃·d，ATI、GTI依次增加673.6℃·d、1078.3℃·d。未来80a，在RCP4.5 排放情景下，自然保护区AFI、GFI分别减少322.8℃·d、357.6℃·d，ATI、GTI分别增加262.2℃·d、419.7℃·d；在RCP8.5 排放情景下，冻融指数的变率更大。

表7 RCP4.5 和RCP8.5 情景下未来80a 自然保护区冻融指数的变化值

3 结论与讨论

本文利用1971～2019 年其周边同属羌塘高原湖盆高寒草原区5 个气象站的逐日平均气温、地表平均地温资料，统计分析了保护区大气和地面冻融指数特征及其时空变化，并预估在RCP4.5 和RCP8.5 两种排放情景下，未来80a（2021～2100 年）大气和地面冻融指数的演变情况，得到如下主要结论：

（1）自然保护区ATI和GTI总体上呈自西向东递减的分布，并随海拔升高而减小；AFI和GFI的分布规律不明显。从冻融指数的大小来看，GTI>AFI>ATI>GFI。与青藏高原及周边地区比较，自然保护区的FI（AFI、GFI）低于青藏高原永久冻土区，高于青藏高原季节性冻土区、黄河源区和雅鲁藏布江中下游，TI（ATI、GTI）正好相反。

（2）近49a（1971～2019 年）自然保护区AFI、GFI分别以−8.97℃·d·a−1、−10.45℃·d·a−1的速率显著减少。ATI、GTI则表现为显著增加趋势，平均每年分别增加7.05℃·d 和11.38℃·d。与青藏高原对比，AFI、GFI减幅小，ATI增幅接近，GTI增幅大。

（3）在10 年际尺度变化上，FI在1970s～1990s 为正距平，2000s～2010s 为负距平，表现为逐年代递减的变化特征；而TI相反，呈逐年代递增的变化特征。在时间转折上，近49a 自然保护区AFI、ATI、GFI、GTI的突变点分别出现在2001年、1993年、1999年和1998年；ATI突变时间最早，较AFI偏早8 年。

（4）近49a 自然保护区FI（TI）减少（增加），与平均气温、平均最低气温显著升高、降水量增加、平均风速减小密切相关，还与积雪日数、最大积雪深度的减少有关。

（5）在RCP4.5、RCP8.5 两种排放情景下，近期、中期、远期自然保护区都表现为FI减少、TI增加的变化特征，以远期最为明显。中等排放（RCP4.5）情景下，未来80a 自然保护区AFI、GFI将分别减少322.8℃·d、357.6℃·d，ATI、GTI分别增加262.2℃·d、419.7℃·d；RCP8.5 排放情景下，冻融指数的变化率变得更大。

青藏高原FI显著下降，TI显著上升，这与高原气温明显升高关系密切，本文分析表明冻融指数与平均气温、平均最低气温存在极显著的相关关系。近50a来青藏高原气温上升率高达0.37℃/10a，远高于全国的增温率（0.16℃/10a）[28]，且研究时段距今越近，气温增温率越大，表明近期增暖更为明显[29−32]。在气候变暖背景下，本区冻融指数变化趋势与青藏铁路沿线[12−13]、青藏高原[14]、黄河源区[21]以及雅鲁藏布江中下游地区[15]是一致，也与我国北疆[33]、东北[34]、黑河流域[11]，以及北美、加拿大和阿拉斯加[35]等地相同。这也表明，FI下降、TI上升在北半球已是普遍现象。尽管20 世纪90 年代出现了全球变暖的停滞期[36−38]，但自1998年以来，青藏高原TI明显上升，FI明显下降[14]。这意味着，青藏高原永久冻土区将加速变暖。

本文基于地面气象观测数据，分析了自然保护区冻融指数的时空演变特点，这既有助于掌握季节性冻土的变化特征，也对草地生态系统退化研究有重要的意义。但仅从单相关系数分析了气温、降水、日照时数对冻融指数的影响，未考虑其他影响因素，这些关系背后的机理是什么，有待下一步研究。