东北地区大豆冷害指标构建及变化特征

李海楠 朱丽杰 李明倩 姜丽霞 任传友 高西宁3)*

1)(沈阳农业大学农学院， 沈阳 110866) 2)(黑龙江省气象科学研究所， 哈尔滨 150030) 3)(辽宁省农业气象灾害重点实验室， 沈阳 110166)

引 言

大豆是我国主要油料作物之一，播种面积仅次于玉米、水稻、小麦，种植区分布广泛[1]。东北地区作为中国大豆的主产区，大豆品种资源丰富，年平均播种面积及产量占全国的50%以上[2-3]。东北地区大豆一般在4月下旬—5月中旬播种，多在9月成熟，由于大豆性喜温，在生育期内对温度变化敏感，并且东北地区处于中高纬度，春、秋两季气温变化较大，使得大豆在生产过程中易受低温影响，遭遇阶段性低温冷害现象频发，导致大豆产量波动较大[4-7]。国内外对大豆低温冷害的研究多在机理方面[8-11]，对东北大豆冷害指标构建及灾害时空变化特征的研究报道较少。

冷害指标是冷害监测预警和评估的重要依据，是冷害研究的基础[12-13]。国内学者将延迟型冷害指标分为5类：温度距平指标、积温指标、发育期距平指标、热量指数指标和综合指标[14-16]。温度距平指标是利用5—9月平均气温总和的负距平判断冷害发生及灾害等级[17-18]，该指标计算相对简单，准确度较高，可较好地判识典型冷害年，但判识严重冷害年比例偏高，未考虑各发育阶段所需积温，缺少作物各发育阶段之间的联系[19-20]。积温指标主要通过计算作物某一发育阶段的活动积温距平或积温差判断热量条件是否能够满足生长需求[21]，该指标反映的是全生育期的冷害程度，在积温相等情况下，不能反映作物在不同发育阶段受到低温冷害的不同影响[22]。发育期距平指标基于作物发育期进程与低温冷害的较高相关性，采取作物发育期延迟日数表征冷害影响程度[23-24]。低温冷害综合指标主要包含主导指标、辅助指标和参考指标，其核心是主导指标[25]，研究中发现该指标可与作物紧密结合，计算简单，但多用于小范围小尺度研究[26-27]。相比以上指标，热量指数指标根据作物各发育阶段三基点温度，分析不同发育阶段冷害发生情况，充分考虑作物不同发育阶段对热量条件需求程度，具有更明确的生物学意义[28]。王培娟等[29]以热量指数为指导因子，构建东北地区春玉米不同发育阶段冷害指标，为玉米冷害监测预报、风险评估和防灾减灾提供了重要的科技支撑。李新建等[30]通过分析不同播期棉花产量、霜前花率与棉花不同发育阶段的热量指数之间的关系，制定了以热量指数表示的延迟型冷害指标。薛志丹等[31]通过热量指数方法计算大豆种植气候适宜度，对黑龙江省大豆种植进行区划。

利用东北地区大豆发育期观测数据，结合历史气象数据，并收集该区域大豆低温冷害历史灾情数据，从大豆生长所需的热量条件出发，基于热量指数构建大豆低温冷害指标，确定灾害等级阈值，定量分析大豆冷害时空变化特征。通过构建大豆低温指标及其等级阈值，为提高大豆低温监测预警技术、准确评估低温对大豆生长发育的影响，开展大豆防灾减灾气象服务、灾害保险等提供基础支撑。

1 数据与方法

1.1 研究区数据来源与预处理

数据包括气象数据、大豆发育期数据和大豆历史灾情数据。气象数据、大豆发育期数据来源于中国气象数据网(http:∥data.cma.cn/)，其中选取1971—2020年东北地区98个气象站的气象数据，来源于《中国地面气候资料日值数据集(V3.0)》；选取1992—2020年东北地区42个农业气象站的大豆发育期数据(包括大豆播种、出苗、三真叶、开花、结荚、成熟日期)，其中1992—2013年来源于《中国农作物生长发育和农田土壤湿度旬值数据集》，2014—2020年来自于国家气象中心。大豆历史灾情数据来源于《中国气象灾害大典》(黑龙江卷、吉林卷、辽宁卷)[32-34]、《中国气象灾害年鉴》[35]。研究区域内气象站和农业气象站分布如图1所示。

图1 研究区域气象站与农业气象站分布Fig.1 Distribution of meteorological stations and agro-meteorological stations in the taget region

对于日平均气温数据中个别缺测值，采用相邻两日的气温平均值代替。农业气象站中大豆发育期年份记录为1992—2020年，由于站点间数据记载起始年份不同，且不同年份各站点发育期数据记录也不完整，1992年以来，研究区各站点6个发育期均有完整记录的站年数为585，占所有发育期站年数的48%。分析各站已有数据发现各发育期在不同年份之间的变化差异不大，为了研究方便，选取播种、出苗、三真叶、开花、结荚、成熟6个发育期，以各农业气象站每个发育期历年普遍期平均日期代表该发育期日期。为了更好地基于气象数据进行分析，本文通过反距离权重插值方法对各农业气象站每个大豆发育期日期进行处理，根据各发育期分布图，得到各气象站点尺度的不同发育阶段。根据受灾情况划分灾害等级，分别为轻度、中度、重度3个等级。

根据研究区域内大豆灾情数据，提取出大豆发生冷害时对应发育阶段热量指数，得到各发育阶段不同冷害等级的样本203个，从中选出171个用于指标构建，随机预留出32个用于指标验证(表1)。

表1 大豆各发育阶段不同冷害等级下样本量Table 1 Sample size in different growth stages of soybean with different levels of chilling damage

1.2 研究方法

1.2.1 热量指数

热量指数方法能够避免其他非热量因素干扰，结合了作物的生长发育特性，同时考虑作物不同发育阶段的三基点温度，其计算公式[29]如下：

F(T)=

(1)

(2)

其中，F(T)为逐日热量指数；T为某发育阶段内日平均气温，单位：℃；T0，T1，T2分别为大豆在某发育阶段的适宜温度、下限温度和上限温度，单位：℃。参照东北地区大豆温度适宜指标体系[36-37]并结合东北地区大豆生产实践，得到大豆各发育阶段的三基点温度(如表2所示)。以大豆不同发育阶段内逐日热量指数平均值作为该阶段的热量指数，每个发育阶段的F(T)值越大，表明热量条件越好，越能满足>该阶段大豆生长需要；反之亦然。

表2 大豆各发育阶段的三基点温度Table 2 Three physiological temperatures in different growth stages of soybean

1.2.2 理论概率分布型拟合与检验

Kolmogorov-Smirnov(K-S)分布拟合检验[38]是比较一个频率分布f(x)与理论分布g(x)或两个观测值分布是否一致的非参数检验方法，多用于样本量较少的情况。在灾害样本序列拟合中，正态分布、指数分布、均匀分布是较为常用的概率密度分布函数[39]。样本序列分布型函数拟合的K-S检验结果在0.05显著性水平上差异不显著(P>0.05)，即可接受原假设，表明样本序列服从该分布类型。

1.2.3 冷害指标等级阈值率定

t分布区间估计方法是在样本量较小，总体服从正态分布且总体方差未知的情况下，通过选取某一置信度，即可得到样本平均值的置信区间。该方法已在建立早稻洪涝和高温热害、夏玉米冷害和涝渍灾害指标时得到广泛应用[40-42]。

1.2.4 Mann-Kendall突变检验

Mann-Kendall(M-K)突变检验是一种无分布检验，既可用于检验序列的变化趋势，也可进行突变点检验，该方法在气象领域应用十分广泛[43]。通过构建顺(UF)、逆(UB)时间序列的秩序列，计算顺、逆向样本序列的统计量，利用两者的变化特点分析样本序列的变化趋势。UF值大于0，则表示序列呈上升趋势，小于0则表示下降趋势，当UF值超过0.05显著性水平临界线时，表明上升或下降趋势显著，两条曲线出现交点，且交点在临界值之间，那么交点对应时间便是突变开始时间。

2 结果与分析

2.1 指标构建

通过K-S分布拟合检验(0.05的显著性水平)对大豆5个发育阶段、3种不同灾害等级下的15组热量指数序列，分别进行正态分布、均匀分布和指数分布模型概率密度函数拟合与检验(表3)，结果表明：有15组冷害样本总体服从正态分布，有11组总体服从均匀分布，3组服从指数分布，在以上3种分布函数对任意一个样本序列的检验中，正态分布P值均为最高。因此本文可以选用通过显著性检验的正态分布作为热量指数序列的最优拟合函数。

表3 3种函数对大豆热量指数样本序列拟合的K-S检验结果Table 3 K-S test for results of 3 functions for fitting heat index samples of soybean

基于正态分布，对各热量指数样本序列进行t分布区间估计，基于置信度为95%的热量指数区间的上限值确定各等级冷害的阈值，得到冷害等级指标(表4)。以第2发育阶段轻度冷害为例说明，它表示当我国东北大豆在出苗-三真叶阶段的热量指数为(67.5，73.5]时，有95%的概率发生轻度冷害。

由表4可以看到，在同一发育阶段内，热量指数越小，大豆发生冷害的等级越高，受灾情影响越严重；在相同的冷害等级下，大豆三真叶-开花阶段、开花-结荚阶段、结荚-成熟阶段的冷害等级区间阈值较高，播种-出苗阶段、出苗-三真叶阶段相比较低，说明大豆在生长发育的中后期热量需求较高，生育前期所需热量条件较低。

表4 基于热量指数的东北地区大豆不同发育阶段冷害等级指标Table 4 Chilling damage level indicators in different growth stages of soybean in Northeast China based on heat index

2.2 指标验证

利用随机预选出的不同发育阶段、不同冷害等级下的32个样本，根据冷害发生的时间、地点、受灾等级，计算受灾地区临近站点热量指数，并与新建立的冷害指标进行一致性检验(表5)。由表5可知，验证的32个样本中，与真实灾害等级完全匹配的共27个，验证准确率达到84.4%；在不同冷害等级内统计，样本验证准确率超过81.8%；各冷害等级下的等级指标在不同发育阶段内分析，样本验证准确率最低出现在播种-出苗阶段，为77.8%，出苗-三真叶阶段准确率为87.5%，三真叶-开花阶段的准确率为80%，开花-结荚阶段准确率最高为100%，结荚-成熟阶段准确率为85.7%。对未通过验证的5个灾害样本再次分析，发现真实的冷害等级均比指标等级低1级。

表5 大豆冷害指标的验证准确率(单位：%)Table 5 Verification accuracy for chilling damage indicator of soybean(unit：%)

2.3 东北地区大豆低温冷害的时间变化特征

2.3.1 变化趋势

图2为东北地区全生育期内冷害次数时间变化曲线，可以看到，1972，1974，1997，2009年的冷害次数相对较多，每个站点的平均冷害次数在1.6次·a-1以上，次数偏多年份主要集中在20世纪70年代，与文献[44]给出的东北地区较强的3年低温冷害年相一致。《中国气象灾害大典(黑龙江卷)》也记载1972年和1976年低温造成的粮食减产幅度超过20%，1983年、1997年也造成大田作物严重受害[32]。

图2 1971—2020年东北地区大豆 全生育期内冷害次数Fig.2 Chilling damage frequency during the whole growth period of soybean in Northeast China from 1971 to 2020

图3为1971—2020年东北地区大豆各发育阶段内不同等级冷害及总冷害次数随时间变化。由图3可知，大豆在结荚-成熟阶段发生总冷害次数最高(平均每站0.36次·a-1)，播种-出苗阶段、三真叶-开花阶段总冷害次数相对较低，出苗-三真叶阶段、开花-结荚阶段最低(平均每站0.17次·a-1)。在同一发育阶段内，冷害次数随冷害等级的增高而减少。

三真叶-开花阶段总冷害次数最高值出现在1983年，平均每站0.72次；结荚-成熟阶段总冷害次数最高值出现在1997年，为平均每站0.99次；两年中重度冷害分别为平均每站0.24次、0.30次。1983年6月黑龙江省出现异常低温，月平均气温较常年偏低2.9℃，为1949年以来最低，吉林和辽宁省内大部分地区作物也受到冷害影响；1997年9月中旬受极地冷空气影响，东北大部分地区出现降温天气，大田作物受害严重。

2.3.2 M-K突变检验

图4为东北地区全生育期冷害次数时间序列的M-K突变检验曲线。由图4可以看到，冷害次数在20世纪70年代初期呈增加趋势，80年代初期也存在小幅度增加。由图4还可以看到，1993年冷害次数出现由高到低的突变，1993年前每站冷害平均次数为1.24次·a-1，1993年后每站冷害平均次数为0.94 次·a-1，且2004年后冷害次数下降趋势明显。

图3 1971—2020年东北地区大豆各发育 阶段不同等级冷害及总冷害次数Fig.3 Chilling damage frequency in different growth stages of soybean with total chilling damage frequency in Northeast China from 1971 to 2020

图4 1971—2020年东北地区大豆全 生育期冷害次数M-K统计量曲线Fig.4 M-K statistic curves of chilling damage frequency during the whole growth period of soybean in Northeast China from 1971 to 2020

2.4 东北地区大豆低温冷害的空间变化特征

2.4.1 全生育期及不同发育阶段冷害频次空间分布

图5为1971—2020年东北地区大豆全生育期和不同发育阶段内冷害频次空间分布。由图5可以看到，全生育期内冷害频次分布大体上呈由南向北增加，黑龙江省最北部大兴安岭地区的漠河、塔河附近地区为冷害的高发地区，平均频次超过3.6次·a-1；吉林省东南部冷害频次也相对较高，平均频次超过2.7次·a-1；辽宁及吉林省大部分地区冷害频次相对较少，平均频次在0.9次·a-1以下。各发育阶段冷害频次在空间分布上呈一致的变化特征，即均以黑龙江省最北部和吉林省东南部为极高值中心。结荚-成熟阶段高值区域最广，从黑龙江省海伦附近一直延伸到黑龙江省最北部，且冷害频次增减梯度明显；出苗-三真叶阶段在整个研究区域的冷害频次分布差距最小，平均频次均在0.8次·a-1以下；出苗-结荚阶段东北大部分地区处在冷害频次低值区，平均频次低于0.2次·a-1。

2.4.2 不同等级下全生育期冷害频次空间分布

对于全生育期，轻度冷害频次高值区主要位于吉林省东部地区，平均每站0.9～1.2次·a-1，中度冷害和重度冷害的频次空间分布特征大致相同，均以黑龙江省北部和吉林省东部为高值中心向四周递减，中度冷害高频次为每站0.6～0.9次·a-1，重度冷害高频次达每站2.4～3.0次·a-1以上(图略)。辽宁省、吉林省大部分地区各等级下每站冷害频次均为0～0.6次·a-1，黑龙江省冷害频次变化梯度较大。

图5 1971—2020年东北地区大豆全生育期和不同发育阶段冷害频次分布Fig.5 Distribution of chilling damage frequency in the whole growth period and different growth stages of soybean in Northeast China from 1971 to 2020

2.4.3 冷害频次年代际变化

全生育期冷害频次的年代际变化如图6所示，冷害高发区在各年代分布比较一致，高值区多位于黑龙江北部、吉林东南部地区，冷害频次平均每站超过4次·a-1，辽宁省和吉林省大部分地区一直处在冷害频次的低值区，平均每站低于1次·a-1。近年东北地区气温逐渐升高，积温显著增加，低温冷害减少[45]。20世纪70—80年代冷害频次高值区在黑龙江省漠河周边扩大明显，低值区也在向北延伸，20世纪80年代至21世纪初高值区逐渐缩小，低值区继续扩大，这与吕佳佳等[46]、王萍等[20]研究的灾害主要发生区域及变化特征相符合。

中度、重度冷害频次年代际分布大致相同，高值区主要位于大兴安岭及长白山附近地区，随年代变化高值区逐渐缩小。轻度冷害频次年代际分布有所不同，高值区位置和范围变化明显，20世纪70年代位于黑龙江、吉林东部地区，70年代后范围缩小，并向黑龙江北部以及吉林东南部地区移动(图略)。

图6 1971—2020年东北地区大豆全生育期冷害频次年代际变化Fig.6 Decadal variation of chilling damage frequency in the whole growth period of soybean in Northeast China from 1971 to 2020

续图6

3 结论与讨论

本文利用1971—2020年东北地区98个气象站的逐日平均气温数据、1992—2020年东北地区42个农业气象站的大豆发育期及历史灾情数据，采用热量指数、K-S分布拟合检验、t分布区间估计、趋势分析、M-K突变检验等方法，构建东北地区大豆低温冷害指标并对灾害变化特征进行分析，主要结论如下：

1) 在同一发育阶段内，热量指数越小，大豆发生冷害的等级越高，受灾情影响越严重；在相同的冷害等级下，大豆三真叶-开花-结荚阶段的冷害等级区间阈值较高，播种-出苗-三真叶发育阶段相对较低，大豆在生长发育的中后期热量需求较高，发育前期所需热量条件较低。

2) 由32个独立灾害样本验证得到灾害等级与冷害指标的匹配率达84.4%。在不同冷害等级下，样本验证准确率超过81.8%；在不同发育阶段，播种-出苗阶段准确率为77.8%，出苗-三真叶阶段准确率为87.5%，三真叶-开花阶段准确率为80%，开花-结荚阶段准确率为100%，结荚-成熟阶段准确率为85.7%。本文构建的东北地区大豆低温冷害指标可以较好反映东北地区大豆冷害发生情况。

3) 东北地区大豆受到低温冷害影响较严重的年份有1972，1974，1976，1983，1995，1996，1997，2009年，冷害次数在1993年前后发生突变，1993年后呈下降趋势，2004年后下降趋势明显。大豆在结荚-成熟阶段发生冷害次数最高，平均每站0.36次·a-1，播种-出苗阶段、三真叶-开花阶段冷害次数相对较低，出苗-三真叶阶段、开花-结荚阶段最低，平均每站0.17次·a-1。在同一发育阶段，不同等级下冷害次数随冷害程度的增高而减少。

4) 东北地区大豆冷害频次分布由南向北增加，黑龙江省最北部大兴安岭地区的漠河、塔河附近地区为冷害高发地区，吉林省的东南部长白山地区为次高发地区；各发育阶段内的冷害频次在空间分布上呈一致变化，即均以黑龙江省最北部和吉林省东南部为极高值中心。20世纪70—80年代冷害频次高值区在黑龙江省漠河周边扩大明显，低值区也在向北延伸，20世纪80年代至21世纪初高值区逐渐缩小，低值区逐渐扩大。

研究发现黑龙江省北部的大兴安岭地区以及吉林省东南部长白山地区冷害频发，虽然从气象上看这些地区低温冷害相对严重，但从大豆种植分布看，这些区域多以林地为主，大豆种植比例相对较低，低温所带来的经济损失低于其他主产区。本文基于热量指数构建大豆冷害指标，从生物学角度出发，考虑农作物的三基点温度，排除非热量因素干扰，将农作物生长发育特征更好地结合在一起。指标使用方便，所需数据获取简单，根据大豆所处发育阶段、三基点温度及日平均气温数据，便可进行大豆冷害监测，具有实际的生物学意义。在对大豆发育期数据处理方面，由于某一农业气象站各发育期气象条件存在年际差异，同时部分年份数据不完整，因此以发育期历年平均值代表各站发育期日期[29,40]；由于农业气象站与气象站数量和空间分布不一致，本文采用空间插值方法获得各站点尺度大豆发育期日期。以上处理方法可能造成气象条件异常年份获得的大豆发育期日期与实际生产中的日期存在误差，而发育期在热量指数计算过程中至关重要，该问题的解决有赖于合理设置农业气象站数量和空间分布，规范农业气象站农田管理及观测技术，结合遥感及大数据技术等农情监测，获取更加完备的农业气象数据。东北地域广阔，种植大豆品种会有差异，本文所采用的三基点温度参数只能反映部分品种温度需求，且指标构建在产量方面考虑欠缺，后续研究需要参考更多大豆品种的三基点温度，并考虑结合作物生长模型开展研究以完善指标体系[47]。