1901—2018年三江源地区大气冻融指数变化特征

陈方方， 罗栋梁， 刘 磊， 金会军，3， 李 超

（1.中国科学院西北生态环境资源研究院冻土工程国家重点实验室，甘肃 兰州 730000； 2.中国科学院大学，北京 100049；3.东北林业大学土木工程学院寒区科学与工程研究院，黑龙江 哈尔滨 150040； 4.兰州交通大学测绘与地理信息学院，甘肃 兰州 730070）

0 引言

冻土是高寒生态系统的重要组成部分，其变化深刻影响地表水热平衡、生物地球化学循环乃至整个高寒生态系统的稳定［1］。在气候变暖和人类活动增强叠加作用下，冻土发生了显著退化，表现为地温升高、季节冻深变浅、活动层增厚、融区贯通或扩展、多年冻土厚度减薄乃至消失［2-5］，由此引起生态系统稳定性减弱和资源环境承载压力增大等问题［6］。冻结融化指数是冻土研究相关的重要参数之一，也是判断冻土存在状态与长期变化的重要指标之一［7］。在季节冻土区，年冻结指数是计算年最大冻结深度的重要依据，以评估寒区工程施工的冻胀量等［8］。在多年冻土区特别是资料匮乏的地区，年冻融指数是计算多年冻土和活动层厚度时空分布变化的关键因子之一［9］。

冻土区因交通不便和气候严酷造成实测地气温资料匮乏，导致冻融指数研究以前多基于国家基准台站或自动气象站的观测资料计算。赵红岩等［10］基于地面温度观测资料计算了青藏铁路沿线冻融指数，其分析表明地面冻融指数主要受纬度和海拔的影响。姜逢清等［11］研究了青藏铁路沿线冻融指数变化趋势，结果表明1966—2004年冻结指数减小，而融化指数增加。Wu 等［12］研究发现青藏高原多年冻土区的冻结指数减少速率高于季节冻土区，而融化指数增加速率则小于季节冻土区。刘磊和罗栋梁［13］基于国家基准台站数据分析了雅鲁藏布江流域大气/地面冻融指数时空变化，结果表明1977—2017年大气/地面冻结指数均减小，融化指数均上升。Luo等［14］计算了我国东北北部1972—2005年年平均气温和地面温度、大气/地面冻融指数，结果显示在过去40年气候显著变暖背景下，冻结指数减小，而融化指数显著增加。综上可知，以往研究大多基于国家基准台站观测资料，其空间分布极不均匀，计算时段也偏短，针对整个研究区长时间序列的冻融指数变化研究较少；且大都集中于资料较丰富地区，难以充分反映某一地区冻融指数的整体时空变化情况。

三江源地区位于青藏高原中东部，是长江、黄河、澜沧江的发源地，其特殊的地理位置、丰富的自然资源、重要的生态服务功能使其成为我国重要生态安全屏障。目前三江源地区的植被和水源涵养等高寒生态水文环境要素及经济社会可持续发展不仅得到地方政府和行业部门重视，也引起相关科学研究工作者的广泛关注［15-17］。不少研究指出，高寒植被、生物生产力和土壤化学元素循环等相关生态环境要素的变化与冻土环境及其变化有着紧密联系［18］。因此，研究三江源地区冻融指数的时空变化特征，对于了解该区冻土变化及其与生态环境的协同变化具有较重要意义。本文基于英国东英吉利大学（University of East Anglia）大学气候研究中心（Climatic Research Unit，CRU）的逐月气温数据集，并利用国家基准台站数据验证，计算并分析了三江源地区大气冻融指数的时空变化特征，研究结果可在一定程度上反映其气候和冻土变化特征，同时可为该地区高寒生态环境长时期变化研究提供科学依据。

1 数据和方法

1.1 研究区概况

三江源地区包括长江源区、澜沧江源区和黄河源区，是我国重要的固碳、水源涵养地和生态屏障区［19］（图1），其中长江源区以直门达为流域出口、澜沧江源区以昌都为流域出口、黄河源区以唐乃亥为流域出口［20］。区内冻土、冰川、湖泊、高寒湿地等广泛分布，是长江、黄河和澜沧江流域重要的水源地和补给区，被誉为“中华水塔”［21］。区域气候为典型高原大陆性气候，冬季严寒漫长，夏季温暖短暂；温度日较差大而年较差小，雨热同季，干湿季分明［22］。年均气温-5.6～7.8 ℃，年降水量262.2～772.8 mm，年日照时数2 300～2 900 h［6］。因位于大片连续多年冻土区边缘，三江源地区大片连续、不连续及岛状多年冻土和季节冻土交替分布。

图1 三江源地区地理范围及国家基准台站分布Fig.1 Geographical background and the national meteorological stations distributed in the Three-River Source Region

1.2 资料来源

由于该区气象站稀少且分布不均，且大多在1950 年后才建站观测，因而不便利用台站观测数据开展冻融指数长时间序列的时空变化研究。东英吉利大学（University of East Anglia）气候研究中心（Climatic Research Unit）CRU TS 4.03（以下简称CRU）通过整合已有的若干知名数据集，重建了一套覆盖完整、高分辨率且无缺测的月平均地表气候要素数据集，成为一套较为完整的全球百年尺度气候变化数据集。这套资料和中国已有气候数据相比具有如下优点：第一，青藏高原20 世纪前半期器测资料严重缺乏，CRU 资料尽管因复杂地形和插值方法本身问题存在误差，经比较仍可作为有一定信度的参考资料；第二，中国现有百年温度序列的时间分辨率为年或季，而CRU 资料时间分辨率为月，时间分辨率较高；第三，这个温度数据集基于观测结果统计内插得到，故减少了代用资料的不确定性［23］。本文所选取的CRU 逐月气温空间分辨率为0.5°×0.5°，时间范围为1901年1月至2018年12月；共使用174个格点的逐月气温资料分析三江源地区冻结指数和融化指数的时空变化。

为检验CRU 资料在三江源地区的适用性，本研究从中国气象数据网（http：//data.cma.cn）下载了中国地面气候资料月值数据集进行验证。三江源地区共计19 个国家基准台站，考虑到部分台站被撤销、迁移及观测时间不连续等问题，总共选取了12个站点1960 年1 月至2018 年12 月的气温数据分析验证（表1）。评价方法是将国家基准台站的逐月气温与对应CRU 格点1960—2018 年的逐月气温做统计分析。由于CRU空间分辨率为0.5°，单个格点代表空间范围较大，因此考虑海拔差异和气温垂直递减率赋予CRU 数据月均温相应增量［-0.402 ℃·（100m）-1］［24］，在此基础上比较基准台站和CRU 格点逐月气温统计差异。具体赋值过程为以基准台站海拔为基准，根据其与对应CRU 格点的海拔差值，赋予CRU逐月气温相应增量。

表1 三江源地区国家基准台站与CRU格点逐月气温差值Table 1 Monthly air temperature difference between CRU grid and the national meteorological stations in the Three-River Source Region

1.3 研究方法

1.3.1 冻融指数计算

大气冻结/融化指数（以下均简称冻结/融化指数）的定义为：某个给定时段内日平均气温低于/高于0 ℃的持续时间和其数值乘积的总和，以℃·d或℃·mon 表示。本研究基于逐月气温计算冻融指数，其计算公式如下：

式中：FI、TI分别为冻结指数和融化指数（℃·d）；Ti为第i月的月平均气温（℃）；Di为第i月天数。本研究中，根据每个格点1901年1月到2018年12月的逐月气温计算每个格点1901—2018 年的年冻融指数。根据每个格点所有年份的冻融指数求其平均值，得到1901—2018年的三江源地区平均冻融指数。

1.3.2 滑动平均法

滑动平均法又称移动平均法，是对时间序列数据进行平滑和滤波的一种处理方法。通过计算一定滑动长度的平均值，目的是消除统计序列中的随机波动，用平均值显示时间序列的变化趋势［25］。本文利用滑动平均法对冻结/融化指数1901—2018 年时间序列数据进行平滑和滤波。对于样本量为n的序列y，其滑动平均序列可表示为：

式中：k为滑动长度。本文利用Python 3.7编程实现三江源地区1901—2018 年冻融指数的滑动平均处理。根据1901—2018 年所有格点冻结/融化指数平均值，得到三江源地区1901—2018 年冻结/融化指数年际变化。同时选择k=11 年对三江源地区冻结/融化指数的年际变化做滑动平均分析，并根据分析结果对冻结/融化指数滑动平均曲线分段线性拟合。

1.3.3 线性倾向估计

利用一次线性倾向分析冻结/融化指数的变化特征。其公式为：

式中：y为冻结/融化指数；x为年份；a为回归常数；b为反映冻结/融化指数变化速率的回归系数（显著水平P＜0.05）。本文利用Python语言的numpy程序包对三江源地区1901—2018 年冻结/融化指数进行一次线性倾向估计，获得每个格点冻结/融化指数变化的线性倾向率，进而通过空间分析获得整个三江源地区冻融指数线性变化的区域差异。考虑到冻结/融化指数自身的空间差异性，本文通过相对变化量来反映每个格点冻融指数相对其自身的变幅，将1901—2018 年的差值与平均值的比值作为相对变化量。

1.3.4 冻结融化指数空间分布

基于各格点冻结/融化指数平均值与其地理信息的多元统计关系，利用ArcMap 空间插值获得其空间分布特征。对每个格点1901—2018 年的冻结/融化指数求平均值，得到其多年平均值。利用三江源地区500 m 分辨率高程数据重采样获得对应格点平均海拔。基于每个格点的经纬度、海拔与冻结融化指数建立其与经纬度、海拔的多元统计关系。同时，为消除量纲影响，对海拔数值乘以0.01 之后再进行回归分析。最终获得如下关系：

式中：FI、TI分别为冻结指数和融化指数；x、y、h分别为经度、纬度和海拔。在ArcMap 中基于多元线性关系对冻结和融化指数进行空间插值。

2 结果与分析

2.1 CRU数据与国家基准台站数据的比较

CRU 数据集已在北半球及中国部分地区的气候变化和多年冻土研究中广泛使用，具有较高的精度和适用性［26-27］。在青藏高原地区，CRU 资料也有一定的可信度。本文利用相关系数（CRU 数据和国家基准台站数据相关程度）、解释方差R2（CRU 数据可解释国家基准台站数据的方差百分比）、均方根误差RMSE（CRU 数据与国家基准台站数据偏差的平方和与观测次数比值的平方根）及平均绝对误差MAE 来分析CRU 数据在三江源地区的适用性。结果表明，国家基准台站的逐月气温与对应CRU 格点逐月气温数据具有较好的一致性（图2）。两者变化范围均为-30～20 ℃，RMSE=3.29，MAE=1.89。相关系数0.92 及解释方差R2=0.93 表明CRU 逐月气温与国家基准台站数据具有很强的正相关性，适用于三江源地区气温及冻结融化指数的时空变化分析。

图2 CRU逐月气温与国家基准台站逐月气温相关分析Fig.2 Correlation of monthly air temperatures between the CRU and the national meteorological stations in the Three-River Source Region

2.2 冻结指数时空变化特征

三江源地区多年平均冻结指数的空间分布较好地体现了其与海拔的相关性：海拔自西向东逐渐递低，冻结指数也随之减少（图3）。同时多年冻土的连续性自西向东降低，由大片连续多年冻土过渡到不连续及岛状多年冻土甚至季节冻土，冻结指数在空间分布上也呈现了明显的自西向东递减趋势：由最西端的3 400 ℃·d减小到最东端的不到600 ℃·d，平均值为1 940 ℃·d，大部分地区的冻结指数值界于800～2 600 ℃·d。其中冻结指数最大值出现于长江源区西部，最小值则分布在黄河源区东南部。

图3 三江源地区1901—2018年多年平均冻结指数空间分布Fig.3 Spatial distribution of average air freezing index from 1901 to 2018 in the Three-River Source Region

三江源地区1901—2018 年冻结指数整体呈下降趋势，其递减斜率为-1.1 ℃·d·a-1，但并非呈一致性下降趋势，而是随气候波动经历了下降-上升-下降等3 个波动变化阶段（图4）。其中，1901—1943年为波动下降阶段，其变化趋势为-3.4 ℃·d·a-1；1943—1966 年呈波动上升趋势，其变化趋势为8.8 ℃·d·a-1；1966—2018 年又呈波动下降趋势，其变化趋势为-4.3 ℃·d·a-1。

图4 三江源地区1901—2018年冻结指数时间序列变化Fig.4 Time series of air freezing index from 1901 to 2018 in the Three-River Source Region

三江源地区1901—2018 年冻结指数线性倾向率变化范围-2.6～-0.11 ℃·d·a-1［图5（a）］，在空间上自西向东呈减少趋势，以长江源区的变化速率最大，黄河源区的变化速率最小。除长江和黄河源区北部冻结指数变化不显著外，其他区域均显著变化。由图5（b）可知，冻结指数相对变化量具有自南向北减少趋势，黄河源区冻结指数的相对变化量较小。总体而言，三江源地区冻结指数的相对变化量均较小，最小为-0.1，表明三江源地区冻结指数的变化幅度不大。

图5 三江源地区1901—2018年冻结指数线性倾向（a）（实心圆表示显著，P＜0.05）及其相对变化量（b）区域差异Fig.5 Spatial difference of the changing rates（a）（solid circle denotes the change is significant，P＜0.05）and the relative change（b）of the air freezing index from 1901 to 2018 in the Three-River Source Region

2.3 融化指数时空变化特征

融化指数的空间分布与冻结指数呈相反趋势，即随海拔升高而减小，随海拔降低而增大，故三江源地区西部融化指数较小，而东部融化指数较大，在空间分布上自西向东递增（图6）。三江源地区融化指数变化范围为0～1 800 ℃·d，平均值为787 ℃·d。其中长江源区融化指数最小，而黄河源区融化指数最大。

图6 三江源地区1901—2018年多年平均融化指数空间分布Fig.6 Spatial distribution of average air thawing index from 1901 to 2018 in the Three-River Source Region

三江源地区融化指数在1901—2018 年整体呈波动上升趋势，其递增斜率为0.34 ℃·d·a-1，具体而言经历了上升-下降-上升等3 个阶段（图7）。其中，1901—1943 年为波动上升阶段（3.3 ℃·d·a-1），而1943—1981 年为波动下降阶段（-3.1 ℃·d·a-1），1981—2018 年为波动上升阶段（2.9 ℃·d·a-1）。融化指数上升和下降的斜率大致相等，表明变化速率大致相同。

图7 三江源地区1901—2018年融化指数时间序列变化Fig.7 Time series of air thawing indices from 1901 to 2018 in the Three-River Source Region

三江源融化指数1901—2018 年间的线性倾向率变化范围为-0.29～0.95 ℃·d·a-1，且大部分地区融化指数变化斜率大于零，表明融化指数呈上升趋势，只有小部分地区呈下降趋势［图8（a）］。在空间分布上，融化指数的变化率自西向东呈减小趋势，以长江源区融化指数递增速率最大，黄河源区的递增速率最小。融化指数相对变化量的变化范围较大，为0.07～0.86，大多集中在0.1～0.3［图8（b）］。从空间分布来看，三江源地区融化指数1901—2018年的相对变化量以西部长江源头和中部最大，黄河源区东部融化指数的相对变化量较小。

图8 三江源地区1901—2018年融化指数线性倾向（a）（实心圆圈表示显著水平P＜0.05）和相对变化量（b）区域差异图Fig.8 Spatial difference of the change rate（a）（solid circle denotes the change is significant，P＜0.05）and the relative change（b）of the air thawing index from 1901 to 2018 in the Three-River Source Region

3 讨论

3.1 CRU数据的适用性

气候变化研究通常需使用百年以上长时间序列的气候资料，而青藏高原器测资料时间短、观测站稀缺且分布极不均匀，因此长时间序列气候和环境变化以及多年冻土的变化常基于ERA、CRU、MERRA2 等再分析资料进行研究［28-30］。目前CRU逐月气温再分析资料在冻土变化研究中仍存在不足。一方面，CRU 逐月气温是月尺度资料，相较于日均气温在冻融指数计算方面存在一定偏差，特别是冻结和融化过程频繁交替的夏初或秋冬月份偏差尤其大。当月内逐日气温均高于或低于0 ℃时，基于逐月气温的计算自然可真实反映冻融指数变化；但在冻融过程频繁交替时，即逐日气温在月内围绕0 ℃反复波动时，基于逐月气温则较难获得真实冻融指数。青藏高原作为典型的高海拔冻土区，日内温差变化大，其在春季和秋季的冻融交替过程中，温度在0 ℃附近频繁波动，因此基于CRU 等再分析资料获取准确的冻融指数比高纬度地区面临更严峻的挑战［31］。另一方面，由于CRU 数据空间分辨率为0.5°×0.5°，每个格点所包含的实地面积约为3 000 km2，而青藏高原地形异常复杂，在较小空间范围内海拔变化可达数百米乃至上千米。考虑到气温的垂直递减率，基于再分析资料计算冻融指数势必存在低估或高估的现象；与青藏高原相比，CRU 数据在北半球其他地区的适用性更高［31］。若要精准获得青藏高原地区冻融指数的时空变化规律，需在考虑复杂地形和局地因子影响的基础上对现有再分析资料降尺度以发展更高时空分辨率的气候数据。在今后研究工作中，拟将CRU 逐月气温与已有逐日气温再分析资料或空间分辨率更高的再分析资料［25，32］与地形因子校正、实地观测资料相结合，以尽量减小误差。

3.2 冻融指数的时空变化及其影响

对比冻融指数的空间分布，可以发现冻结和融化指数的空间分布大致呈相反趋势，即自西向东，随着海拔逐渐下降及多年冻土连续性逐渐降低，冻结指数递减，融化指数递增；表明从长江源区到黄河源区年平均气温大致在上升。由冻结指数与经纬度和海拔的多元线性关系可知，冻结指数与纬度和海拔正相关，而与经度负相关（图1和图3）。三江源地区海拔自西向东在不断降低，冻结指数与经度的负相关从一定程度上可理解为冻结指数与海拔正相关。因此，三江源地区冻结指数的空间分布呈现出海拔越高、冻结指数就越大。由图3 可知，由南向北，冻结指数由小于800 ℃·d逐渐增加到2 600 ℃·d，体现了冻结指数在空间分布上的纬度效应。这与曹斌等［33］在研究黑河流域冻结指数的空间分布，何彬彬等［34］研究北疆地区冻结指数的影响因素的结果一致，即冻结指数主要受纬度和海拔的双重影响。

冻结指数为一年中的负积温，融化指数则为一年中的正积温。冻融指数的变化不仅反映了土壤冻结状态的变化，也在一定程度上反映了气温变化和多年冻土对气候变化的响应幅度。气温变化是导致冻融指数变化的直接原因，其剧烈程度与冻融指数的变化幅度相一致。本文研究表明三江源地区1901—2018 年间冻结指数呈波动下降、融化指数呈波动上升趋势，这大体与易湘生等［27］、陈晓光等［35］、张士锋等［20］的结果一致，即三江源地区自二十世纪五六十年代以来气温呈波动上升趋势，以冬季增温最显著。本文的研究表明，冻融指数气候倾向率以长江源区和澜沧江源区的变化速率最大，大致呈自西向东递减的趋势，这表明长江源区和澜沧江源区是整个三江源地区气温升幅最大的地区。这与窦睿音［36］对三江源地区近半个世纪的气候变化研究相一致，即澜沧江源区与长江源区增温速率最快，分别为0.35 ℃·d·a-1、0.32 ℃·d·a-1；而黄河源区增温速率最小，仅为0.19 ℃·d·a-1。此外，冻融指数的气候倾向率随海拔升高而增大，反映出地温较低、连续性较高的多年冻土其退化更严重。Ran等［37］发现青藏高原海拔3 600 m 处年均气温增温速率为0.33 ℃·d·a-1，而海拔5 200 m 处的增温速率为0.49 ℃·d·a-1，这与三江源地区自西向东海拔不断降低、高海拔大片连续多年冻土区增温速率最大的结果较为一致。

冻融指数高低反映了土壤的冻融起止时间及持续时长，因此对生态系统、地表水文过程、农业生产、基础设施等都会产生重要影响［8］。冻融指数作为生态系统的一个重要指示因子，受季节冻融循环影响明显的生态系统过程可分为生长季、土壤冻结初期、冻结期和融冻期等四个关键时期［38］。冻融指数通过决定植物生长发育时期的时长，从而进一步影响动植物的生长进程、动物群落的结构和功能以及相关土壤生化特性。融化指数增大、冻结指数减小，意味着融化期延长、冻结期缩短，表明植被生长季将延长，植物生产力可能得到提高，高寒植被的生态服务功能也将提升；反之，则植被生长季缩短，植物生产力减弱，其生态服务功能亦可能削弱。整体上，三江源地区1901—2018 年间冻结指数减小、融化指数增大，表明高寒植被在过去100 年间可能有了更长的生长季，因而高寒生态环境可能有了一定程度的改善。

4 结论

本文基于国际上广泛使用的CRU 再分析资料，对1901—2018 年间三江源地区的冻结和融化指数的时空变化及其差异进行了比较分析，主要结论如下：

（1）三江源地区冻结指数的变化范围为600～3 400 ℃·d，大多数冻结指数集中分布于800～2 600 ℃·d；在空间分布上，纬度和海拔越高，冻结指数就越大。融化指数的变化范围为0～1 800 ℃·d，空间分布呈现自西向东递增的趋势。

（2）三江源地区冻结指数在1901—2018 年整体呈波动减少趋势，递减斜率为-1.1 ℃·d·a-1，大致经历了下降（1901—1943 年）-上升（1943—1966 年）-下降（1966—2018 年）等三个阶段；融化指数整体波动增加，递增斜率为0.34 ℃·d·a-1。

（3）三江源地区冻结指数线性倾向斜率变化范围-2.6～-0.11 ℃·d·a-1。在1901—2018 年间呈降低趋势。融化指数线性倾向斜率变化范围-0.29～0.95 ℃·d·a-1，大部分地区融化指数呈上升趋势。自西向东，冻结指数递减斜率在减小，融化指数递增速率也在减小。冻结指数相对变化量都较小，最小为-0.1，说明三江源地区冻结指数的变化幅度不大。而融化指数相对变化量的变化范围较大，为0.07～0.86 ℃·d·a-1，但大多集中在0.10～0.30 ℃·d·a-1。

（4）冻融指数的时空变化在一定程度上反映了气候变化对高寒植物生长的影响，三江源地区冻融指数过去100年的时空变化特征可能意味着高寒植被有了更长的生长季。