基于格点资料的1961−2018年中国霜冻灾害时空变化规律*

白 磊，张 帆,2，文元桥，师春香，吴 晶，尚 明，朱 智，孟俊耀

（1．武汉理工大学智能交通系统研究中心，武汉 430070；2.武汉理工大学航运学院，武汉 430070；3.国家气象信息中心，北京 100081；4.兰州中心气象台，兰州 730020；5.河北工程大学地球科学与工程学院，邯郸 056038；6.北京驰远科技有限公司,北京 101108）

霜冻是由于日最低气温下降，使植物茎、叶处温度下降到0℃以下，导致生长的植物受到冻伤的现象。霜冻日期的变化会显著影响作物产量和生育期，从而通过肥料使用和土壤状态影响农业碳排放[1−3]。1976−2015年中国农业生产受到霜冻灾害损失率呈上升趋势[4]，因此，还需要在气候变化背景下深入研究霜冻灾害精细化时空变化规律，降低霜冻灾害所带来的损失。

1901−2018年全球平均升温0.85℃，而中国升温幅度显著高于全球水平[5]。剧烈的升温会影响霜冻灾害的时空分布，造成霜冻事件大范围减少，无霜期增加[6−7]。在区域尺度上，中国东部地区[8]、东北地区[9]和西北地区[10−11]也表现出初霜日延后，终霜日提前，无霜期延长。上述研究均在站点数据基础上使用反距离插值等算法将分析结果外推到空间上进行分析，在站点分布稀疏和地形复杂地区，存在较大的空间不确定性。

在霜冻灾害的辨识中，通常使用地面0cm 处日最低温度0℃和日最低气温0℃作为阈值，低于这些阈值温度将判断为霜冻灾害发生[10,12−13]。除此，由于百叶箱中气温比地面0cm 温度高约2℃，又有研究将日最低气温2℃作为阈值温度进行霜冻灾害辨识标准。在新疆和东北地区，气温标准辨识精度高，而在华北和陕北地区地温标准辨识精度较高[14]。由此可见，辨识标准存在较强的空间异质性。由于气温相对于地温较容易实现空间外推插值，因此，在大范围区域研究霜冻灾害时可以方便使用，但不同日最低气温辨识阈值的空间辨识效果尚有差异。

目前，已有多套高时空分辨率的格点化气象数据，使用这些数据可以减少站点分析结果外推插值带来的空间不确定性问题，但仍需分析霜冻灾害在格点上的辨识阈值。本研究使用1961−2018年长时间序列的格点化气象数据，研究适合中国区域的格点数据的霜冻辨识阈值、霜冻发生时间及霜冻强度的时空变化趋势，以期为在气候变化背景下合理利用气候资源，规避潜在的气候风险提供数据基础和科学依据。

1 资料与方法

1.1 数据

格点化气象数据为中国气象局国家气候中心制作的CN05.1 数据。该数据集基于中国2400 余个国家级地面气象台站的逐日观测资料，利用薄样板样条插值方法外推插值生成一套中国区域长时间序列格点化观测数据集[15]。格点化数据在制作过程中，一般会对站点数据进行严格的质量控制（数据缺测、均一化等），同时空间外推能够基于合理的假设，避免了传统外推插值中由于站点空间分布稀疏导致外推插值的空间格局分布不合理的现象。该数据集包括日平均、最高和最低地面2m 高度气温等变量，空间分辨率为0.25°，空间覆盖范围为中国大陆地区，时间分辨率为天，覆盖时间长度为1961−2018年。CN05.1 数据已经过细致评估，能够相对准确地表示各气象变量的实际变化情况[16]，可作为中国区域格点化近地面气象场资料。

1.2 方法

1.2.1 霜冻灾害的辨识

相对于地表温度，气温容易在空间上进行外推插值，因此，选择2m日最低气温作为霜冻辨识所需格点化气象数据。分别以日最低气温0℃和2℃作为阈值温度计算初（终）日和无霜日，比较不同日最低气温阈值温度辨识差异，最终选取合适的阈值作为辨识标准。

霜冻日辨识方法：（1）8月1日−12月31日，首个日最低气温低于阈值温度的日期为初霜日；（2）1月1日−7月31日，最后一个日最低气温低于阈值温度的日期为终霜日；（3）在同一年度，终霜日与初霜日间隔日数为无霜期。若上年度未出现，则上年度初、终霜日为空值。

1.2.2 霜冻强度的辨识

不同作物有不同的霜冻强度判断标准[17−19]。参照中国气象行业标准《作物霜冻害等级》[20]和中国气象局《气象灾害预警信号即防御指南》[21]标准，以日最低气温作为判断变量，在每年终霜日与初霜日之间的生长季将霜冻灾害划分为轻度、中度和重度三个强度等级（表1）。

表1 基于2m日最低气温（Tmin）的霜冻强度等级划分标准Table 1 The thresholds of frost intensity based on daily minimum 2m air temperature(Tmin)

1.2.3 霜冻指标趋势分析

以1961−1990年为基准期，计算初霜冻日期、终霜冻日期、无霜期长度和生长季内轻度霜冻、中度霜冻和重度霜冻时间序列的距平变化，并进行7点平滑处理，获得1967−2018年的时间序列，最后使用线性回归方法计算线性变化趋势。

1.2.4 霜冻灾害区划方法

在格点化气象数据的基础上，利用空间聚类分析方法研究中国范围内的霜冻灾害区划。聚类算法采用Python 的sklearn 库中Kmeans 算法，聚类特征采用霜冻日期信息（初霜冻日期、终霜冻日期、无霜期长度）和霜冻强度等级气候态信息（生长季内轻度霜冻、中度霜冻和重度霜冻）。聚类过程中，在比较霜冻日期信息、强度信息和综合两者信息作为特征聚类效果后，选择空间分布特征合理的霜冻灾害区划结果。

2 结果与分析

2.1 霜冻辨识阈值的确定

初霜日、终霜日和无霜期的空间分布具有明显的纬度和高程梯度分布特征。图1a1 和图1a2 分别是0℃标准和2℃标准辨识的初霜日空间分布，由图可见，两者空间分布情况基本一致。新疆天山、阿勒泰山、青藏高原和东北北部地区最早9月前出现初霜冻；9月中旬，新疆北部、内蒙古、东北地区和青藏高原的柴达木盆地等地出现初霜冻；9月中旬−10月中旬，新疆南疆、河西走廊、黄土高原和华北平原等地出现初霜冻；10月下旬−11月下旬，黄淮、江淮和长江中下游流域等地出现初霜冻；四川盆地、云贵高原和华南的大部地区在12月出现初霜冻。对于终霜日日期（图1b1、图1b2)，两个标准辨识的空间格局也基本相似。四川盆地、云贵高原和华南大部地区的终霜日在1月；新疆南疆、河西走廊、黄土高原和华北平原的终霜日在3月；新疆北部、内蒙古、东北地区和青藏高原的柴达木盆地终霜日在4月；新疆天山、阿勒泰山、青藏高原和东北北部地区终霜日在6月。初霜日和终霜日在新疆北疆、青藏高原南部、长江中下游等地区变化剧烈，标准差在10d 以上。而在北方大部分地区标准差在5～10d，变化相对较小。由于0℃标准（图1c1）相对于2℃标准（图1c2）初霜日相对滞后，终霜日相对提前，使得0℃标准无霜期天数普遍多于2℃标准的无霜期。无霜期最长地区主要分布于33°N 以南地区，超过240d。东北和西北北部地区无霜期低于120d。受到地形影响，在天山区、青藏高原地区的无霜期不足60d。

在北方，大部分地区0℃标准辨识的初霜日（图1a1）晚于2℃标准（图1a2）10～20d。0℃标准能够辨识出新疆准噶尔盆地西侧低海拔地区和中天山巨大的高程差异引起初霜日梯度变化，能辨识出新疆东部伊吾淖毛湖地区的初霜日晚于邻近的哈密巴里坤地区10～20d；东北辽河平原晚于松嫩平原10d 的空间格局，而2℃标准下无明显差异；川渝地区和湖南、江西等地相差10d 的空间格局；湖南、江西和福建为非均匀的空间结构，而2℃标准下为均匀空间分布。对于终霜日，两个标准辨识的主要差异在新疆准噶尔盆地、黄土高原中部和长江以南地区。0℃标准能捕捉到江西和湖南地区交接的丘陵地带终霜日日期晚于周边河流平原地区的空间分布格局。对于2℃阈值标准，细节则不明显。可见，0℃标准辨识霜冻现象相比2℃空间上更加合理，故本研究以0℃阈值作为标准辨识霜冻现象。

图1 以日最低气温0℃（1）和2℃（2）为辨识标准分析的多年（1961−2018年）气候态下霜冻日期及其标准差的空间分布Fig.1 Spatial pattern of multi-year climatological mean frost date in 1961−2018 and their standard deviations detected by the threshold of daily minimum temperatures of 0℃（1） and 2℃（2）,respectively

2.2 霜冻强度变化分析

图2 为不同强度霜冻多年平均发生频数的空间分布和不同时期差值。由图2a1 可见，对于轻度霜冻，东北、西北、青藏高原西部和华南北部地区每年发生3～6 次；黄土高原、内蒙古高原、大/小兴安岭和长白山地区中轻度霜冻每年发生6～8 次；黄淮地区、山东半岛，轻度霜冻每年发生8～12 次；青藏高原东部和新疆天山每年发生10～16 次。对于中度霜冻（图2b1），黄土高原、蒙古高原、大兴安岭北部、长白山、青藏高原东部和山东半岛地区，每年发生3～4.5 次；西北地区（除新疆南疆）和长三角地区每年发生1～2 次。对于重度霜冻（图2c1），内蒙古中部地区每年发生1.5～2.5 次；黄土高原大部分地区、山东西部、河南东部、东北北部和东部地区发生0.5～2.5 次。

图2 轻度（a）、中度（b）和重度（c）等级霜冻1961−2018年平均发生次数（1）及其与1961−1990年差值（2）的空间分布Fig.2 Spatial pattern of the multi-year mean annual frequency(1)of light frost(a),moderate frost(b)and severe frost(c)in 1961−2018 and the difference(2)between the occurrences times of different frosts intensity between 1991−2018 and 1961−1990(former−latter)

对比1991−2018年和1961−1990年可以发现，两个时期霜冻强度变化的空间格局存在明显的差异。在图2a2 中，与1961−1990年相比，1991−2018年轻度霜冻在中国北方大部分地区年平均发生次数增加1～2 次，而在长江中下游等南方地区减少2～3 次。对于中度霜冻（图2b2），除东北山区、黄土高原和秦岭等地区呈现增加态势外，1991−2018年中国大部分地区中度霜冻年平均发生次数减少0.5～1.5 次。对于重度霜冻（图2c2），除内蒙古中部地区在1991−2018年年平均发生次数减少了0.5～1 次，全国大部分地区重度霜冻增加0～0.5 次。整体上，1991−2018年相对于1961−1990年不同强度的霜冻在中国北方呈微弱增加状态，在南方呈现明显减少状态。

图3 为不同强度早霜冻和晚霜冻在1961−2018年的空间分布。由图可见，总体上，不同强度晚霜冻的发生频次高于早霜冻。对于轻度霜冻，在中国东部两者的空间分布不同。早霜冻在北方较多，晚霜冻在淮河以南的地区相对较多。对于中度霜冻，早霜冻和晚霜冻整体上的空间格局基本一致。在内蒙古中东部、吉林和辽宁地区，分别为1.8～2 次·a−1和3～4 次·a−1。对于重度霜冻，早霜冻和晚霜冻空间分布大体一致。在内蒙古中东部、吉林和辽宁地区早霜冻和内蒙古中部的晚霜冻是发生次数均为0.4～1 次·a−1。

图3 1961−2018年早霜冻（1）和晚霜冻（2）中轻度霜冻（a）、中度霜冻（b）和重度霜冻（c）发生频数的空间分布Fig.3 Spatial pattern of the frequency of light(a),moderate(b)and severe(c)frosts in the early frost period(1)and late frost period(2)frosts in 1961−2018

图4 展现了不同强度多年平均早（晚）霜冻发生次数的空间分布和不同时期差异。由图可见，对于轻度霜冻，早（晚）霜冻在华北和东北地区，1991−2018年比1961−1990年多年平均增加0.4～0.8（0.6～1.2）次，在长江以南地区多年平均减少0.4～0.8（0.6～2.1）次。对于中度霜冻，早霜冻和晚霜冻的空间分布格局基本一致，都是秦岭淮河以南的地区1991−2018年比1961−1990年多年平均呈减少态势，北方（除内蒙古中部）呈现分区的增加态势。对于重度霜冻，早（晚）霜冻在东北（新疆北疆、华北和内蒙古东部）呈现增多态势，在华北、内蒙古中部呈现减少态势。

图4 轻度（a）、中度（b）和重度（c）早霜冻（1）和晚霜冻（2）在1991−2018年与1961−1990年发生次数差值的空间分布Fig.4 Comparison of the multi-year mean annual occurrences of early frosts(1)and late frosts(2)with different intensities[light frost(a),moderate frost(b)and severe frost(c)]between 1991−2018 and 1961−1990

2.3 霜冻趋势变化分析

由图5a1 可见，1961−2018年中国大部分地区初霜日以1～3d·10a−1（P＜0.05）的速率呈显著推迟的趋势，在新疆北疆和青藏高原趋势高达4～6d·10a−1（P＜0.05）。对于终霜日（图5b1），黄土高原、河西走廊、新疆大部分地区以0～1·10a−1的微弱趋势推后，青藏高原、天山、东北北部、淮河流域以3～5d·10a−1（P＜0.05）的速率呈显著提前趋势。对于无霜期（图5c1），中国大部分地区以1～3d·10a−1（P＜0.05）的速率呈显著增加趋势；淮河流域变化速率为6～7d·10a−1（P＜0.05）；西部的青藏高原、天山等地区，显著增加趋势在7d·10a−1（P＜0.05）以上。对于轻度霜冻年发生次数（图5a2），南方地区以0.2～0.6 次·10a−1的速率略微减少，而在北方以0.2～0.4次·10a−1的速率轻微增加。对于中度霜冻年发生次数（图5b2），内蒙古中部、青藏高原东部和淮河流域以0.1～0.4 次·10a−1的速率略微减少，华北和东北地区以0.1～0.3 次·10a−1的速率略微增加。对于重度霜冻年发生次数（图5c2），内蒙古中东部大部分以0.1～0.3 次·10a−1的速率略微增加，淮河流域以0.3～0.5次·10a−1的速率略微减少。

图5 1961−2018年初霜日(a1)、终霜日(b1)和无霜日(c1)及轻(a2)、中(b2)、重(c2)度等级霜冻发生次数变化趋势的空间分布Fig.5 Spatial pattern of early frost,last frost and frost-free days’ linear trend(Figs.a1,b1 and c1)and different frost intensities’ linear trend(Figs.a2,b2 and c2)in 1961−2018

2.4 霜冻气候分区及其变化特点分析

2.4.1 霜冻灾害聚类划区

图6 中基本输入为每个格点的经纬度。图6a 是额外增加霜冻日期信息（初霜冻日期、终霜冻日期和无霜日日数）作为输入的聚类。图6b 是增加多年平均不同霜冻强度发生次数信息（轻度霜冻、中度霜冻和重度霜冻）作为输入的聚类。图6c 为综合全部信息作为输入的聚类结果。前两类聚类与原始的数据输入大体一致。综合考虑的信息与考虑霜冻日期信息大体一致，说明Kmeans 算法给霜冻日期信息的权重高于不同霜冻强度发生次数信息。在图3c 中，蓝色区域为亚热带分区，包含华西、华中、华东和华南北部地区；红色区域为暖温带分区，地理上包含华北、黄土高原南部和新疆天山以南的南疆区域；绿色区域为中温带分区，包含西北大部分地区，东北地区的东北和南部地区；粉红区域为寒温带区域，包括东北西部和北部、青藏高原边缘地区和新疆天山和阿勒泰山与盆地过度的地区；黄色区域为青藏高原分区，包含青藏高原和新疆天山、阿勒泰山地区。该分区结果与中国气候分区大体一致，能够在空间上进行合理解释，故以图6c 中的霜冻区划作为分区依据。

图6 以多年平均霜冻日期（a）、霜冻等级（b）及其综合（c）作为特征的霜冻灾害区划Fig.6 Frost hazards regionalization using the multi-year mean frost date features(a),multi-year mean annual occurrences of frost intensities(b)and mixed features of all above(c)

2.4.2 各区霜冻情况变化

由图7 可见，整体上，1961−2018年中国不同霜冻分区内霜冻日期（图7a、c 和e）格点数（面积）平均距平发生了显著变化。在5 个分区中，青藏高原霜冻日期变化最为剧烈。在1991−2018年，青藏高原的初霜日日期延后趋势为4d·10a−1，终霜日日期提前趋势为 5d·10a−1，无霜期日数增加趋势为9d·10a−1。亚热带变化幅度在这5 个分区中最小。对于轻度霜冻年发生次数（图7b），1991−2018年28a间5 个分区均呈现减少趋势，其中以亚热带地区减少最为剧烈。对于中度霜冻年发生次数（图7d），亚热带在近28a 间呈现稳定减少趋势，而中温带和寒温带呈现剧烈波动变化。对于重度霜冻年发生次数（图7f），仅有中温带在近28a 间呈现显著增加趋势，其他分区无显著变化。

图7 不同霜冻日期和不同等级霜冻年发生次数区域平均距平的年际变化趋势Fig.7 The interannual variation of regional area mean anomalies’ different frost dates and annual occurrences of different frost intensities

2.4.3 各区霜冻情况覆盖面积变化

图8a 中，初霜日日期在8月1−11日和11月1日−12月 1日，1961−1990年平均格点数多于1991−2018年平均格点数。在先前的分析中，青藏高原和新疆天山等高寒山区（图6c、图a1）在8月初发生初霜冻，华北等暖温带地区在11月发生初霜冻。说明青藏高原和部分暖温带的初霜日日期在推后。10月 1−31日，1991−2018年平均格点数多于1961−1990年平均格点数，说明西北和东北等中温带地区初霜日日期整体在推后。在图8c 中，1月11日−2月10日、3月2−22日、4月1−30日、6月1−20日，1991−2018年平均格点数多于1961−1990年平均格点数，说明5 个分区终霜日格点不同程度推后，以暖温带（3月2日−22日）幅度最大。在图8e 中，60～90d、150～200d 和280～350d，近28a 格点数明显增加。整体上，无霜日相邻日期分区从左到右呈现“减少—增加”规律，前期分区减少的格点数，在后续日期分区增加。对于轻度霜冻（图8b），以多年平均发生6 次为界，＜6 次的格点数在增加，≥6次的格点数在近28a 减少。对于中度霜冻（图8d），1991−2018年发生3 次的格点数少于1961−1990年平均格点数。重度霜冻（图8f），多年平均发生3 次的格点数同中度霜冻。从3 个不同霜冻强度发生的格点数变化可以发现，1991−2018年霜冻发生的强度相对1961−1990年整体呈减弱趋势。

图8 1961−1990年与1991−2018年不同霜冻日期和霜冻强度年发生频次，对应的格点面积对比Fig.8 Comparison of the grid area of different frost dates and different frost intensity ’s annual frequency

3 结论与讨论

3.1 讨论

对于霜冻的时空演变规律，在站点尺度已有若干研究。在以往研究中，由站点计算的初霜日变化趋势在新疆[10]、青海[22]、山西[23]分别为 0 ～1.2d·10a−1、3d·10a−1和2d·10a−1，本研究格点数据计算的初霜日变化趋势分别为0～3d·10a−1、2～3d·10a−1和1～3d·10a−1。可以看出，在小区域内，格点数据计算结果接近已有站点数据计算结果，但有更强的趋势变化。在西北地区[24]和全国范围[7]，本研究初霜日变化趋势均为1.4d·10a−1，与先前研究1.8d·10a−1和1.4d·10a−1接近。对于终霜日变化趋势，本研究中国变化范围为−3～1d·10a−1，以往研究平均为−2 d·10a−1，基本能够反映出终霜日提前的现象。对于无霜期，本研究呈现的延长趋势为1～5d·10a−1，而基于站点的趋势为3.4 d·10a−1。由此可见，无论初、终霜日，还是无霜期天数，格点尺度计算出来的趋势总体上与站点计算的趋势一致，但强度上有显著的差异。在空间分布上，格点数据与站点数据计算结果空间结构相似，但基于格点数据的趋势空间分布有更多的细节。对于初、终霜冻日期，本研究格点化结果，基本接近上述区域研究结果。

现有研究在中国大陆范围基本使用700 余站点，本研究中0.25°格网有15000 以上的格点，这使西部站点分布稀疏的地区数据的空间代表性增强。然而，以中国作为研究范围，0.25°格网基本适用，但对于省或者县域精细化的研究工作中（如霜冻保险等），则需要更高分辨率的数据区刻画局部小气候。目前，中国气象局的陆面同化系统CLDAS[25]，能够实时提供1km 高空间分辨率的气象数据，可以为今后的区域尺度精细化研究提供数据支持。

3.2 结论

（1）日最低气温0℃相比2℃作为阈值辨识初/终霜冻和无霜期的空间分布更合理，更适合辨识霜冻。初/终霜冻日期和无霜期天数空间分布具有随着纬度和高度变化明显的地理分布特征。在中国大部分地区，初霜冻日期呈现1～3d·10a−1显著延后的趋势，在新疆北疆和青藏高原显著的推后趋势可达4～6d·10a−1。对于终霜冻日期，在青藏高原、天山、东北北部、淮河流域呈现3～5d·a−1显著提前趋势。对于无霜期，中国大部分地区呈现1～3d·10a−1显著增加趋势

（2）对于不同霜冻强度的年发生次数，1991−2018年相对于1961−1990年，不同强度的霜冻年累计发生次数在中国北方有微弱的增加态势，南方呈现明显减少的态势，以轻度霜冻最为明显。对于中度/重度霜冻，内蒙古中东部增加趋势最为显著。全年内，晚霜冻不同霜冻强度年发生频次高于早霜冻。

（3）综合考虑霜冻日期信息（初霜冻日期、终霜冻日期和无霜期日数）和多年平均不同霜冻强度发生次数信息（轻度霜冻、中度霜冻和重度霜冻）的霜冻气候分区，可将全国分为青藏高原、寒温带、中温带、暖温带和亚热带5 个分区。对于霜冻日期，青藏高原变化最为剧烈，亚热带变化幅度最小。在青藏高原，初霜日以4d·10a−1的趋势延后，终霜日以5d·10a−1的趋势提前，无霜期日数以9d·10a−1趋势增加。对于霜冻强度年发生次数，5 个分区均呈现减少趋势，在1991−2018年，亚热带地区减少最为明显。

（4）从不同霜冻日期的格点数统计看，青藏高原变化剧烈，其次是中温带和寒温带。从霜冻强度年发生次数对应格点数统计看，中度强度以上霜冻强度面积呈减弱状态；轻度霜冻频发地区（＞6次·a−1）格点数呈减弱状态，年发生次数较少的地区呈增加状态。