基于格點資料的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）格點數呈減弱狀態，年發生次數較少的地區呈增加狀態。