中国东北兴安岭地区年冻融频次的分布规律

魏 豪， 张 泽， Andrey MELNIKOV， 金豆豆， 高思如， 冯文杰

（1.中国科学院西北生态环境资源研究院冻土工程国家重点实验室，甘肃兰州730000； 2.东北林业大学土木工程学院/寒区科学与工程研究院，黑龙江哈尔滨150040； 3.东北林业大学东北多年冻土区地质环境系统教育部野外科学观测研究站/东北多年冻土区环境、道路建设与养护协同创新中心，黑龙江哈尔滨150040； 4.俄罗斯科学院西伯利亚分院梅尔尼科夫冻土研究所，俄罗斯联邦雅库茨克117997； 5.中国科学院大学，北京100049）

0 引言

兴安岭地区位于中国东北北部，包括从内蒙古赤峰市到黑龙江漠河市的大兴安岭，以及北起黑龙江岸南至松花江岸的小兴安岭，有中国面积最大、保存相对完好的原始森林，是重要的林业基地和东北边疆的重要门户。兴安岭地区包含大部分多年冻土区和部分季节冻土区，与俄罗斯贝加尔湖和远东地区的不连续多年冻土以及蒙古国山区多年冻土相邻，是中国高纬度多年冻土的主要分布区和寒区工程建设的重要地区。

近年来全球气候变暖日趋严重，IPCC 报告指出在1880—2012年间，全球的近地层年平均气温升高了0.85 ℃［1-2］，中国东北地区的多年冻土和季节冻土都有不同程度的退化，多年冻土区和季节冻土区的南缘都有明显的北移现象［3］。20 世纪50 年代至21世纪初，东北地区多年冻土总面积从4.8×105km2减少为3.1×105km2，冻土南缘向北移动0.1°～1.1°，平均海拔升高160.5 m［4］。同样地，东北地区地表冻结指数的下降和融化指数的上升与全球变暖、冻土退化有着相同的趋势［5-6］。因此，对兴安岭地区的冻土变化及寒区工程建设进行调查研究尤为重要。

俄罗斯冻土学家Ершов［7］在《寒区岩土成因论》中提到，“在自然条件下，一年当中寒区表层的岩土会发生100 余次冻结-融化过程（温度正负变化100余次）。”在一年内发生冻融循环的次数即为年冻融频次，在寒区工程建设和使用过程中，冻融循环是对各类材料的力学性质影响最大的因素，因此年冻融频次的变化对寒区工程的建设、发展以及长期稳定性有着极其重要的影响。

在寒区研究方面，多集中在冻融指数、冻土分布等方面。冻融指数是给定时期内的温度累计值，是对冻土季节性变化进行量化的气候参数。Zhang等［8］结合多年冻土区的年平均地表温度和冻结指数进行研究，发现东北多年冻土区在20 世纪80 年代至90 年代冻土退化严重，普遍由低海拔向高海拔、由连续冻土向不连续冻土退化，在此之后多年冻土区面积略有增加；任景全等［9］通过对吉林省季节冻土区冻融指数的研究发现，在空间上吉林省的冻结指数呈现由北向南逐渐减小的趋势，在时间上冻结指数呈现逐渐减小的趋势。

在室内试验方面，多着眼于工程材料，如岩石、混凝土等在多次冻融循环下的强度等力学性质的变化情况。Jamshidi 等［10］、方丽莉等［11］研究发现冻融循环会通过改变土的结构而导致土的工程性质的改变，进而影响基础设施的耐久性，甚至会导致工程失稳等问题；有学者指出，长期的冻融循环对岩石产生冻融风化、开裂等作用，对许多岩石工程的稳定性产生显著影响［12-15］。这使得表述冻土分布和寒区气候变化的冻融指数在工程实践中存在一定的局限性。而年冻融频次是定量表述寒区冻融循环强弱程度的物理量，能够使实验数据直接应用于寒区工程中，同时也可以和冻融指数进行对比研究，从而更进一步完善寒区研究。

因此，本文基于中国东北兴安岭地区18个气象站点的地表温度数据，对兴安岭地区年冻融频次在空间上的分布情况和在时间上的变化情况进行调查研究，以期为之后在兴安岭地区的气候变化、寒区工程、农林业生产发展等方面的研究提供可靠的数据支撑。

1 数据与方法

1.1 研究区概况

兴安岭地区有中国保存比较完好、面积最大的原始森林，森林覆盖率达80.95%，区域范围跨越黑龙江省和内蒙古自治区，属于温带大陆季风性气候。其中，大兴安岭平均海拔为573 m，冬季寒冷干燥且漫长，夏季温暖湿润而短暂，区域内河流较多，以黑龙江水系和嫩江水系为主。兴安岭地区是中国的北部边境地区，包含极北部存在的不连续多年冻土区和一部分岛状多年冻土区以及南部的中-深季节冻土区，其冻土分布情况复杂，年冻融频次变化差异较大，地理环境多样。

本文选取兴安岭及其周边部分地区作为研究区域。采用的气象数据为兴安岭地区1990—2017年18 个气象站0 cm 处逐日最高地表温度、逐日最低地表温度以及日平均地表温度，图1 为各站点分布情况，图2为调查区冻土分布情况［16］。

图1 兴安岭地区地形图及气象站点分布Fig. 1 Topographic map of Xing’anling region and distribution of meteorological stations

图2 兴安岭地区冻土分布［16］Fig. 2 Distribution of seasonally frozen soil，island permafrost and discontinuous permafrost in Xing’anling region［16］

由于兴安岭地区各站点所处地理位置不同，气候环境不同，各站点间年冻融频次差异较大，全部站点在各年份的平均冻融频次不能够详细地描述不同气候环境下的年冻融频次变化情况。研究区由南到北跨越中温带和寒温带两个温度带，由西向东跨越半干旱区、半湿润区和湿润区三个气候区。根据中国气候区划新方案（2010 年）［17］的东北地区气候区划结果对大兴安岭山脉沿线进行分区研究，如图3 所示。寒温带湿润区包含漠河、塔河、新林、图里河等4个站点；中温带湿润区包括呼玛、加格达奇、孙吴、嫩江、黑河等5个站点；中温带半湿润区包括额尔古纳右旗、小二沟、博克图、扎兰屯、阿尔山、索伦等6个站点；中温带半干旱区包括满洲里、海拉尔、新巴尔虎右旗等3个站点。

图3 兴安岭地区气候区划［17］Fig. 3 Division of climate zone in Xing’anling region［17］

1.2 NDVI数据

本文采用的NDVI 数据为Land Term Data Record（LTDR）NDVI（AVH13C1）和MODIS NDVI（MOD13C2）两种NDVI数据集的结合。其中1981—1999 年的NDVI 数据来自LTDR AVH13C1，该数据属于日NDVI数据，其空间分辨率为0.05°，获取数据的网址为http：//ltdr.nascom.nasa.gov/cgi-bin/ltdr/ltdrPage. cgi。2000—2014 年的NDVI 数据来自MODIS NDVI，MOD13C2属于半月合成数据，其空间分辨率为0.05°，获取数据的网址为https：//ladswe.nascom. nasa. gov。这两种数据都已经经过几何校正、大气校正、辐射校正等预处理［18］，并且经过验证，由LTDR 和MODIS 生成的1981—2014 年逐月NDVI数据可靠，可以用到本文的研究中［19］。

1.3 研究方法

1.3.1 年冻融频次的定义

根据实验室中对混凝土、岩石等寒区工程中涉及到的材料进行冻融循环实验时所用到的冻融循环方法，为方便后续研究者们将实验所得数据更好地应用到实际工程建设中，本文规定以一个测量日内日最高地表温度大于0 ℃，同时日最低地表温度小于0 ℃为一次冻融循环。一天内土壤由夜间最低温度升至白天最高温度时，经历一次融化过程，再由白天最高温度降至夜间最低温度时，经历一次冻结过程，视为一次完整的冻融循环。统计一个测量年内冻融循环的次数即为该年的年冻融频次，以次/年表示。冻融循环次数由逐日观测的日最高地表温度和日最低地表温度数据计算得出。

1.3.2 时空变化趋势分析

在年冻融频次时间变化趋势的分析中，采用线性倾向估计的方法研究年冻融频次在长时间序列上的变化趋势及其倾向率。在年冻融频次空间变化趋势的分析中，基于ArcGIS 10.2，将多年平均冻融频次和DEM 数据采用IDW（inverse distance weighted）插值方法进行插值，并采用线性拟合分析方法进一步研究年冻融频次与经度、维度和海拔之间的变化关系。

1.3.3 突变分析

Pettitt 方法是通过检验序列均值变化的时间来确定序列突变点，是一种基于Mann-Whitney的非参数检验方法，可以对存在趋势变化的序列进行突变检验。本文采用Pettitt检验法对兴安岭地区18个站点在长时间序列上的变化情况进行突变分析。在突变分析中需要用到的非参数统计量为

根据以上非参数统计量可以计算出

一般情况下，认为当P≤0.05时被检验序列中存在突变点。

1.3.4 权重分析

采用多元线性回归权重分析的方法对兴安岭地区年冻融频次在地理因素上的分布情况进行权重分析。在进行多元线性拟合之前，为去除经纬度、海拔数据在数值上对权重分析的影响，首先对经纬度和海拔数据进行标准化处理，然后对年冻融频次与经纬度和海拔的标准化数据进行多元线性回归拟合，整理出经纬度和海拔对年冻融频次影响的权重。

2 结果与分析

2.1 年冻融频次的时间变化

对兴安岭地区18 个站点的年冻融频次数据求平均值，得到全部站点在各年份的平均冻融频次。由图4 可知，兴安岭地区年冻融频次在1990—2017年间呈波动下降的趋势，线性拟合曲线斜率为-2.0，即年冻融频次平均每年下降2.0 次/年，通过了0.05 显著性水平检验，年冻融频次最高值达132 次/年，最低值为66次/年。

图4 1990—2017年兴安岭地区年冻融频次的年际变化Fig. 4 Interannual variation of annual freeze-thaw frequency in Xing’anling region from 1990 to 2017

图5为4 个分区的年冻融频次变化情况及其线性拟合线分布。从图中可以看出，4 个分区内年冻融频次的走势与总体的年冻融频次走势相近，随着时间的迁移，呈现明显的下降趋势，但各区年冻融频次的高低和下降趋势的大小不尽相同。图5（a）表明，在兴安岭北部的寒温带湿润区，年冻融频次拟合曲线斜率为-2.0，即寒温带湿润区年冻融频次平均每年下降2.0 次/年；图5（b）表明，在兴安岭东部的中温带湿润区，年冻融频次拟合曲线斜率为-1.8，即中温带湿润区年冻融频次平均每年下降1.8 次/年；图5（c）表明，在兴安岭南部的中温带半湿润区，年冻融频次拟合曲线斜率为-2.3，即中温带半湿润区年冻融频次平均每年下降2.3 次/年；图5（d）表明，在兴安岭西部的中温带半干旱区，年冻融频次拟合曲线斜率为-1.9，即中温带半干旱区年冻融频次平均每年下降1.9 次/年。以上线性拟合均通过了0.05显著性水平检验。

图5 兴安岭地区各气候区年冻融频次的年际变化Fig. 5 Interannual variation of annual freeze-thaw frequency of four climate zones in Xing’anling region：cold temperate humid area（a），middle temperate humid zone（b），middle temperate sub-humid zone（c）and middle temperate semi-arid zone（d）

注意到年冻融频次在时间线上存在跳跃性较大的区域，因此为了进一步分析年冻融频次在时间线变化中是否存在突变现象，本文采用Pettitt 法对兴安岭地区各站点在时间序列上的变化进行突变分析。对气候区划4个分区内各个站点的平均冻融频次进行整合，得到各分区在每年的平均冻融频次，之后对这一平均冻融频次采用Pettitt 法进行整理，得到非参数统计量Ut，Nmax数据。由图6 可知，兴安岭地区内4 个分区的统计量Ut，N的变化情况大致相同，统计量Ut，Nmax均在2004年达到最高值，从寒温带湿润区到中温带半干旱区，各分区内的统计量Ut，Nmax分别为193、182、195、185，起对应的突变检验值P分别为1.08×10-4、3.20×10-4、8.77×10-5、2.39×10-4，所以各分区内检验值P均小于0.05，则可以说明兴安岭地区年冻融频次在时间序列上的变化发生了显著性突变，突变年份为2004年。

图6 兴安岭地区各气候区Ut，N的变化情况Fig. 6 Variation of Ut，N of four climate zones in Xing’anling region：cold temperate humid area（a），middle temperate humid zone（b），middle temperate sub-humid zone（c）and middle temperate semi-arid zone（d）

由于年冻融频次在2004 年发生突变，所以以2004 年为分界线。1990—2004 年，兴安岭地区年冻融频次平均在118 次/年左右，最大为132 次/年，最小为90 次/年；兴安岭北部的寒温带湿润区年冻融频次平均值在115 次/年左右，最大达到129 次/年，最小为96 次/年；兴安岭东部的中温带湿润区年冻融频次平均值在102 次/年左右，最大为114 次/年，最小为78 次/年；兴安岭南部的中温带半湿润区年冻融频次平均值在132 次/年左右，最大为152 次/年，最小为113 次/年；兴安岭西部中温带半干旱区年冻融频次平均值在119 次/年左右，最大为151 次/年，最小为98次/年。

2005—2017 年，兴安岭地区整体多年平均冻融频次在83 次/年左右，最大为95 次/年，最小为66 次/年；兴安岭北部的寒温带湿润区年冻融频次平均值在77 次/年左右，最大达到111 次/年，最小为59 次/年；兴安岭东部的中温带湿润区年冻融频次平均值在73 次/年左右，最大为97 次/年，最小为61 次/年；兴安岭南部的中温带半湿润区年冻融频次平均值在93 次/年左右，最高大为108 次/年，最小为75 次/年；兴安岭西部中温带半干旱区年冻融频次平均值在88次/年左右，最大为102次/年，最小为68次/年。

2.2 年冻融频次的空间变化

由于各站点年冻融频次为长时间序列的数据，为了讨论兴安岭地区各站点间年冻融频次在空间上的变化情况，对各个站点的多年冻融频次数据求平均值，得到各个站点1990—2017 年近30 年年冻融频次的平均值。通过ArcGIS 10.2 对年冻融频次数据进行处理，采用IDW 插值方法将年冻融频次在空间上的变化情况展现在二维等值线地图之上，与等高线插值图进行对比研究，观察年冻融频次与海拔、经纬度之间的关系。图7 为兴安岭地区平均冻融频次插值图，可以看出，兴安岭地区年冻融频次整体上呈现西高东低、南高北低的趋势。冻融频次最大值出现在索伦站，平均冻融频次为139 次/年；最小值出现在黑河站，平均冻融频次为83次/年。

图7 兴安岭地区多年平均冻融频次的插值分布图Fig. 7 Interpolation distribution of averaged annual freeze-thaw frequency in Xing’anling region

同时通过线性回归分析的方法，对多年平均冻融频次与海拔、经纬度进行线性拟合，更加精确地分析年冻融频次在空间因素变化下的变化趋势以及线性倾向率。如图8（a）所示，对兴安岭地区各站点的平均冻融频次与纬度进行线性拟合可以看出，兴安岭地区各站点间的平均冻融频次在由南到北呈现不断递减的线性分布，确定系数R2为0.434，随着纬度的增加，年冻融频次大体上以5.5 次/度的速率递减，拟合线通过了0.05 显著性水平检验。图8（b）为各站点的平均冻融频次随经度变化情况以及线性拟合，可以看出兴安岭地区各站点的平均冻融频次由西向东呈现不断递减的线性分布，确定系数R2为0.284，随着经度的增加，年冻融频次大体上以2.4 次/度的速率递减，拟合线通过了0.05 的显著性水平检验。

如图8（c）所示，兴安岭地区各站点的平均冻融频次在垂直高度上随海拔的升高呈现不断递增的线性分布，确定系数R2为0.140。随着海拔的升高，年冻融频次大体上以23 次/千米的速率递增，拟合曲线未通过0.05 显著性水平检验。斜率接近于0，说明年冻融频次受海拔影响不大。

图8 兴安岭地区多年平均冻融频次随经纬度及海拔的变化Fig. 8 Variation of averaged annual freeze-thaw frequency in Xing’anling region with latitude（a），longitude（b）and altitude（c）

由上述可知年冻融频次在空间变化上会受到经纬度以及海拔的影响，并且随纬度升高而降低，随经度增加而降低，随海拔升高而升高。为了进一步探究影响年冻融频次在空间上变化的因素中经纬度和海拔的影响权重是多少，本文对经度、纬度和海拔在年冻融频次的影响程度上利用多元线性回归分析的方法进行权重分析。首先对经度、纬度和海拔的数据进行标准化处理，使原始数据去量纲化，从而消除原始数据在单位上的限制，使不同单位或量级的指标进行拟合，更好地进行权重分析。然后通过最小二乘法，以多年平均冻融次数为因变量，以经度、纬度和海拔的标准化数据为自变量，进行多元线性拟合，得到表1中的数据。

从表1 可以看出，在满足0.05 显著性水平检验的情况下，经度、纬度和海拔对兴安岭地区年冻融频次多元线性拟合的拟合系数分别为-20.020、-28.580 和-6.332。如果仅考虑地理信息对兴安岭地区年冻融频次的影响，对拟合系数进行标准化处理，得到经度、纬度和海拔对兴安岭地区年冻融频次影响的权重分别为36.5%、52.0%和11.5%，由此可知在影响年冻融频次空间变化的因素中，纬度影响最大，其次是经度，海拔影响最小。

表1 经纬度及海拔对兴安岭地区年冻融频次影响的权重Table 1 Weight of influence of longitude，latitude and altitude on annual freeze-thaw frequency in Xing’anling region

2.3 年冻融频次与NDVI变化情况对比

研究发现，北半球约70%的多年冻土地区在生长季（4—10 月）和72%的地区在秋季期间归一化植被指数增加，冻土变暖与绿化呈正相关关系［20］。中国东北多年冻土区植被在生长季的NDVI 变化情况，在空间上大致表现为由西向东逐渐增加的趋势，在时间上随着时间的增加呈现显著的增加趋势。如图9 所示，中国东北多年冻土区NDVI 由1990年的0.56增加到2014年的0.65。

图9 1990—2014年东北地区NDVI的年际变化Fig. 9 Interannual variation of NDVI in Northeast China from 1990 to 2014

冻土区地表年冻融频次除了会受到经纬度和海拔的影响之外，同时也受到植被覆盖度、山坡坡向等多种因素的综合影响。从图10 可以看出，中国东北多年冻土区归一化植被指数NDVI 随时间序列呈显著的波动上升趋势，而年冻融频次呈波动下降趋势。采用Pettitt 法对东北地区NDVI 数据进行突变分析，得到NDVI 变化的统计量Ut，N在2000年取得最大值-150，其对应的突变检验值P=0.00049＜0.05，说明东北地区NDVI 在2000 年发生了突变。

图10 兴安岭地区NDVI与年冻融频次的年际变化对比Fig. 10 Comparison of interannual variation between NDVI and annual freeze-thaw frequency in Xing’anling region

图11将大兴安岭山脉周边地区的年冻融频次与NDVI 变化情况进行相关性拟合，在拟合过程中令年冻融频次滞后于对应NDVI 年份4 年，发现随着植被覆盖率的增加，兴安岭地区的年冻融频次呈显著下降趋势，经线性拟合的确定系数达到0.460，由此可知兴安岭地区年冻融频次与NDVI变化情况具有良好的负相关性。兴安岭地区年冻融频次在2004 年发生突变，与NDVI 数据对比有明显的滞后性，说明年冻融频次变化会受到当地植被覆盖率的影响，且这种影响有一定的滞后性。植被覆盖率的增长会逐渐减缓地表对热能的接收和散发，这个过程需要一定的时间，从而导致年冻融频次的降低，这个过程需要一定的时间，所以导致年冻融频次变化滞后于NDVI的变化。

图11 兴安岭地区年冻融频次与NDVI散点图及线性拟合Fig. 11 Annual freeze-thaw frequency and NDVI scatter diagram and linear fitting in Xing’anling region

3 结论

本文通过对中国东北兴安岭地区及周边地区18个站点1990—2017年的年冻融频次进行整理，讨论了兴安岭地区年冻融频次在时间和空间上的分布情况，并利用权重分析、突变分析等方法，分析了地理因素对年冻融频次变化的影响以及年冻融频次的突变情况，并结合东北地区NDVI 变化情况进行了对比。得出以下结论：

（1）兴安岭地区年冻融频次在空间上大体上呈现随经度、纬度升高而降低，随海拔升高而升高的趋势，进一步分析得到，经度、纬度和海拔在影响年冻融频次变化中以经纬度影响为主，海拔影响较小。

（2）受到气候变化的影响，兴安岭地区年冻融频次随着时间的推移，整体呈现明显下降的趋势，且在2004 年前后发生突变。同时进行年冻融频次与NDVI数据对比分析后发现，年冻融频次与NDVI间有较好的负相关性，且东北地区NDVI 在2000 年左右发生突变，年冻融频次的突变与之相比，存在一定滞后性，说明年冻融频次在一定程度上受到当地植被覆盖率变化情况的影响。

调研年冻融频次分布情况可以使实验室中得到的岩石冻融循环试验数据直接对接寒区工程，对寒区工程的建设有着重要的指导意义，也可弥补冻融指数在工程实践中的不足之处，同时对寒区的气候变化和农林业发展等有着重要的借鉴意义。但本文数据存在空间分辨率和精度上的不足，在讨论年冻融频次的影响因素时，也只考虑了经纬度和海拔的影响，对其他影响因素暂未考虑，期待未来加深研究。