曹娟^,1, 张朝^,1, 张亮亮¹, 骆玉川¹, 李子悦¹, 陶福禄²1. 环境演变与自然灾害教育部重点实验室/地表过程与资源生态国家重点实验室 北京师范大学地理科学学部,北京 100875
2. 中国科学院地理科学与资源研究所 中国科学院陆地表层格局与模拟重点实验室,北京 100101

Damage evaluation on soybean chilling injury based on Google Earth Engine (GEE) and crop growth model

CAO Juan^,1, ZHANG Zhao^,1, ZHANG Liangliang¹, LUO Yuchuan¹, LI Ziyue¹, TAO Fulu²1. State Key Laboratory of Earth Surface Processes and Resource Ecology/MEM & MoE Key Laboratory of Environmental Change and Natural Hazards, Faculty of Geographical Science, Beijing Normal University, Beijing 100875, China
2. Key Laboratory of Land Surface Pattern and Simulation, Institute of Geographic Sciences and Resources Research, CAS, Beijing 100101, China

通讯作者: 张朝(1971-), 女, 教授, 主要从事农业灾害、粮食安全以及环境风险方面的研究。E-mail: zhangzhao@bnu.edu.cn

收稿日期:2019-06-15修回日期:2020-04-8网络出版日期:2020-09-25

基金资助:

国家自然科学基金项目.41977405 国家自然科学基金项目.31561143003 国家自然科学基金项目.41571493 国家自然科学基金项目.41621061

Received:2019-06-15Revised:2020-04-8Online:2020-09-25

Fund supported:

National Natural Science Foundation of China.41977405 National Natural Science Foundation of China.31561143003 National Natural Science Foundation of China.41571493 National Natural Science Foundation of China.41621061

作者简介 About authors
曹娟(1994-), 女, 博士生, 主要从事农业系统对全球变化的响应研究。E-mail: caojuan@mail.bnu.edu.cn








摘要
作为中国商品粮的主要生产基地,东北地区频发的低温冷害给中国粮食安全带来了严重的影响,及时、准确和大范围评估低温冷害灾损是降低损失、快速恢复生产的重要前提。本文以鄂伦春为例,提出了一种快速评估低温冷害对大豆生产影响的新方法。首先诊断出该地区典型冷害事件发生的年份为1989年、1995年、2003年、2009年和2018年;然后基于本地化后CROPGRO-Soybean模型设置512组低温冷害和田间管理组合模拟情景;其次构建了1600组包括3个变量（有效积温距平值（CDD）、模拟的叶面积植被指数（LAI）和产量）的冷害脆弱性模型;最后依托Google Earth Engine（GEE）平台逐像元提取大豆关键生育期早晚窗口内最大的宽动态植被指数（WDRVI）及对应的日期（DOY）,将WDRVI转化为大豆种植格点的实际LAI,结合构建的冷害脆弱性模型逐像元计算出产量和减产率。主要发现如下：① 校准后的CROPGRO-Soybean模型能较为准确地模拟不同冷害情景下的大豆生长发育过程;② 大豆遭受全生育期的降温情景（1~3 ℃）的减产幅度比局部降温情景（4个生育期随机生成连续5日温度为0 ℃）的减产幅度大;③ 历史冷害年1989年、1995年、2003年、2009年平均减产率约为9.6%、29.8%、50.5%和15.7%,与实际6.4%、39.2%、47.7%和13.2%的减产率相比,冷害灾损评估结果具有较好的精度且误差均在一倍方差以内;④ 运用该方法评估2018年冷害田间尺度的产量损失,并利用Sentinel-2A影像进行10 m高精度制图。结果显示,该方法能够快速、准确地评估不同尺度的低温冷害灾损情况,为作物估产以及农作物灾害损失评估的业务化运行提供了新的思路。
关键词： 低温冷害;Google Earth Engine（GEE）;CROPGRO-Soybean;大豆;减产率;有效积温距平值（CDD）

Abstract
The frequent occurrence of chilling injury has serious impacts on national-level food security, and it mainly affects the grain yield in Northeast China. Timely and accurate measures are desirable to assess the large-scale impacts, which are the prerequisites for disaster reduction and production recovery. Therefore, we propose a novel method to efficiently assess the impact of chilling injury on soybean. Inner Mongolia is taken as a case study. The specific chilling injury events was diagnosed to occur in 1989, 1995, 2003, 2009, and 2018. The 512 combinations of cold and field management simulation scenarios were established based on the localized CROPGRO-Soybean model. Furthermore, we constructed 1600 cold vulnerability models including CDD (Cold Growth Days), simulated LAI (Leaf Area Index) and yields from the CROPGRO-Soybean model. Finally, we extracted the maximum wide dynamic vegetation index (WDRVI) and corresponding date of the critical windows of early and late growing seasons in the GEE (Google Earth Engine) platform, converted the WDRVI into actual LAI of soybean pixel, and calculated the pixel yield and losses according to the corresponding vulnerability models. The findings show that the localized CROPGRO-Soybean model can accurately simulate the growth and development processes of soybean under different cold scenarios. The soybean yields were reduced due to changes in cold stress during the whole growth period (a decrease of 1-3 ℃), which were greater than those from the local cooling treatments (the temperature of 0 ℃ for 5 consecutive days which are randomly generated during four growth periods). Moreover, simulated historical yield losses in 1989, 1995, 2003, and 2009 were 9.6%, 29.8%, 50.5%, and 15.7%, respectively, which were very close (all errors were within one standard deviation) to the actual losses (6.4%, 39.2%, 47.7%, and 13.2%, respectively). The above proposed method was applied to evaluate the yield loss of 2018 at a pixel scale. Specifically, sentinel-2A image was used for 10 m high-precision yield mapping, and the estimated losses well characterized the actual yield losses from 2018 cold event. The results highlighted that our proposed method can efficiently and accurately assess the chilling injury impact on soybean at different spatial scales. The novel method is also effective for efficient assessment of the impacts of different disasters on other crops.
Keywords：chilling injury;Google Earth Engine (GEE);CROPGRO-Soybean;soybean;yield loss;cold degree days (CDD)


PDF (6246KB)元数据多维度评价相关文章导出EndNote|Ris|Bibtex收藏本文
本文引用格式
曹娟, 张朝, 张亮亮, 骆玉川, 李子悦, 陶福禄. 基于Google Earth Engine和作物模型快速评估低温冷害对大豆生产的影响. 地理学报[J], 2020, 75(9): 1879-1892 doi:10.11821/dlxb202009005
CAO Juan, ZHANG Zhao, ZHANG Liangliang, LUO Yuchuan, LI Ziyue, TAO Fulu. Damage evaluation on soybean chilling injury based on Google Earth Engine (GEE) and crop growth model. Acta Geographica Sinice[J], 2020, 75(9): 1879-1892 doi:10.11821/dlxb202009005

1 引言

大豆是世界上最大的油料作物和重要的蛋白质来源,在世界经济发展中具有重要的地位^[1]。世界上约有2/3的人蛋白质摄量不足,为了解决人口增长与蛋白质紧缺的矛盾,目前许多国家都在积极提高大豆产量^[2]。此外,大豆在中国的播种面积仅次于玉米、水稻、小麦和马铃薯,是最具经济效益和进口需求量最高的作物,因此大豆是关系国计民生的重要基础性、战略性的物资。2019年中央一号文件将“实施大豆振兴计划和奶业振兴行动”作为十大重要点目标之一,明确提出要多途径扩大大豆的种植面积,推进新品种技术示范推广和大豆种植全程机械化,稳定粮食产量,完善大豆生产者的补贴政策。东北地区是中国大豆主产区之一,东北大豆的种植面积和产量位居全国前列,占全国大豆的种植面积和产量的1/3,在中国的大豆产业发展中具有举足轻重的地位^[3]。

全球虽有变暖的趋势,但是作物生长季的增温趋势并不明显,加之近些年来农作物种植带盲目北移和作物晚熟品种的种植面积不断扩大^[4],低温冷害依然是中国北方发生频率最高的农业灾害,也是导致该地区农业灾损最为严重的原因^[5,6]。1951—1980年曾发生8次较强冷害,其中1969年、1972年和1976年3次低温冷害,导致作物减产约57.8亿kg^[7]。东北大豆种植在高纬地带,低温冷害的胁迫会更为严峻^[8]。大豆在生长期间由于热量条件供应不足且温度变幅较大,容易造成生育期延迟和生长不良,出现空壳和秕粒,这些都会降低大豆的产量和品质。因此,快速监测冷害发生的时间并评估低温冷害导致的产量损失显得尤为重要^[9]。

国内外****针对作物低温冷害的监测和产量的损失评估进行了广泛研究,主要分为作物模型方法和遥感方法。虽然校准后作物模型能够实现站点尺度的冷害情景模拟研究,但区域作物模型进行冷害研究需要大区域范围且高精细的地面数据输入信息（品种特性、气象要素、土壤条件、管理措施和田间栽培等）。这些数据在站点尺度易获取,但由于地表环境的非匀质性,将站点尺度数据应用到区域估产上难免会产生误差^[10]。遥感冷害致损监测研究主要是通过遥感手段观测地面植被指数的突变情况,或是通过遥感数据反演气温数据,辅以地面观测数据,构建温度指标识别是否出现冷害及冷害的程度,最后利用冷害指标与受灾年份的产量进行建模分析。遥感方法虽然精度较高,但是受研究区域和年份影响,空间泛化能力较弱,不易于推广。此外,目前针对东北地区玉米和水稻的低温冷害研究较多^[11,12],但是对大豆的冷害研究相对较少。针对上述问题,快速准确地评估东北地区低温冷害对大豆产量的影响显得尤为迫切。

综上,结合作物模型和遥感两种方法的优点进行作物估产和灾损评估是当前研究的热点问题之一。Google Earth Engine（GEE）提供了海量免费共享影像数据并具有强大的数据存储与云计算能力^[13],作物模型CROPGRO-Soybean能动态直观地刻画大豆生长发育和产量形成过程,将两者结合为大豆长势监测和灾害损失评估提供了可能^[14]。本文基于GEE平台并结合CROPGRO-Soybean作物模型尝试提出一种快速的、跨区域的低温冷害对大豆生产影响的评估方法。通过对大豆生长模型进行本地化调参,模拟不同冷害情景并建立大豆冷害脆弱性模型,得到高精度的格点产量和减产率,最后在县级单元上进行精度验证,以期为大豆估产以及低温冷害损失评估提供新的技术和方法。

2 研究区概况与数据来源

2.1 研究区概况

鄂伦春地处大兴安岭山地西北坡,位于大兴安岭山地向呼伦贝尔高平原过渡地段,介于48°50′N~51°25′N和121°55′E~126°10′E之间,总面积约6×10⁴ km²,其中耕地面积约2.1×10³ km²。地形以高原、中山、低山、丘陵为主,属于典型的高原型地貌区。气候属于寒温带半湿润大陆性季风气候,四季明显。年均气温-2.7~0.8 ℃之间,自西向东递增,1月份气温最低可达-30 ℃,七月份气温最高可达37.5 ℃。无霜期平均95 d,年降水量459.3~493.4 mm。风速较小,年均风速1.8~2.9 m/s。鄂伦春是内蒙古自治区重点粮食生产基地,平均年产粮食2亿kg,大豆是该地区种植面积最广的农作物（图1）。

图1

新窗口打开|下载原图ZIP|生成PPT
图1研究区概况

Fig. 1Overview of the study area

2.2 数据与预处理

2.2.1 遥感数据 遥感具有快速、宏观和动态3个重要特点,能够获取大范围地表信息,在作物生长监测和产量预报中发挥巨大作用^[15]。在GEE平台上,Landsat 5地表反射率产品（Landsat 5 Surface Reflectance）和Sentinel-2A地表反射率产品（Sentinel-2A Surface Reflectance）已经过辐射矫正、几何矫正和大气矫正等数据预处理。首先通过JavaScript API在线访问和筛选东北地区大豆关键生育期（5—9月）内Landsat 5（1984—2012年,空间分辨率30 m,重访周期为16 d）和Sentinel-2A（2017—2019年,空间分辨率10 m,重访周期为5 d）所有可用影像,其次利用地表反射率产品的质量控制（QA）波段进行去云处理（Landsat 5和Sentinel-2A影像的去云代码链接分别为https://code.earthengine.google.com/dc10f7a140c2fe51aa30-dfab08和https://code.earthengine.google.com/6c7863acf90a59607fd4ec2cd6beadde）。经过以上预处理之后才能够获得与地面实测曲线拟合度较好的影像光谱曲线,计算准确的植被指数。

2.2.2 气象数据和土壤数据 1981—2018年鄂伦春气象站的日平均、最高、最低气温、降雨和日照时数来自于国家气象信息中心（National Meteorological Information Center）,太阳辐射值参考Angstrom^[16]计算得到。CROPGRO-Soybean模型的土壤输入数据获取自中国土壤科学数据库（Soil Science Database）,包括容重、有机碳、全氮、pH及其粉粒、砂粒和粘粒的百分比。土壤反照率,径流潜力和排水率参照Hoogenboom等^[17]的程序估算获得。影响光合速率的土壤肥力因子（SLPF）设定为1.0,初始土壤含水量为田间持水量。

2.2.3 土地利用数据和田块实测数据 鄂伦春10 m的大豆种植面积图和模型本地化所用到的2014—2017年6个田块实测数据来源于当地保险公司。如表1、表2所示,数据主要包括田块的位置,大豆生长基本数据（品种熟性、产量结构和发育期状况）和田间管理数据（翻耕、除草、灌溉和施肥等）。1981—2015年鄂伦春县级产量数据来源于农业农村部种植业管理司（http://www.zzys.moa.gov.cn）和内蒙古统计年鉴。研究区河流、县级和乡镇行政区划数据来源于中科院资源环境数据云平台（http://www.resdc.cn）的1:25万基础地理信息数据,数据投影基准统一转化为WGS-1984。

Tab. 1
表1
表1鄂伦春6个田块的基本信息
Tab. 1The basic information of six crop fields

基本信息	五岔沟	兴盛	春亭阁	春林	乃木河	小库莫
经度(°)	124.42	124.49	124.49	124.42	124.11	124.17
纬度(°)	50.10	49.97	49.83	49.83	49.4	49.7
高程(m)	486	393	380	365	381	447
品种	黑河38	黑农35	合丰50	合丰25	合丰39	黑河18
熟性	中熟	中熟	中熟	中晚熟	中熟	中熟
栽培方式	平作条播	条播	条播	垄作	条播	平作条播

新窗口打开|下载CSV

Tab. 2
表2
表26个田块的关键生育期、种植密度和产量观测数据
Tab. 2The observation information of key growth period, planting density and yield of crop fields

年份	观测数据	五岔沟	兴盛	春亭阁	春林	乃木河	小库莫
2014	播种日期	5月16日	5月18日	5月1日	5月18日	5月28日	5月11日
	开花普遍期	7月14日	7月12日	6月30日	7月6日	7月18日	6月30日
	开花成熟期	9月26日	9月28日	9月18日	9月28日	9月28日	9月24日
	种植密度(株/m2)	39.69	30.42	47.73	36.75	37.98	38.47
	产量(kg/hm2)	1160	1350	1650	1150	610	1050
2015	播种日期	5月16日	5月16日	5月28日	5月8日	5月17日	5月18日
	开花普遍期	7月16日	7月10日	7月2日	7月12日	7月16日	7月8日
	开花成熟期	9月16日	9月14日	9月28日	9月28日	9月28日	9月8日
	种植密度(株/m2)	47.94	35.6	38.18	26.24	25.84	28.58
	产量(kg/hm2)	980	1100	1410	1030	510	860
2016	播种日期	6月6日	5月16日	5月6日	5月22日	5月18日	5月18日
	开花普遍期	7月20日	7月14日	6月28日	7月28日	7月2日	7月2日
	开花成熟期	9月30日	9月26日	9月18日	9月30日	9月18日	9月14日
	种植密度(株/m2)	32.97	30.55	40.04	42.9	40.8	35.82
	产量(kg/hm2)	1280	1450	1850	1370	670	1110
2017	播种日期	5月24日	5月23日	5月21日	5月18日	5月24日	6月2日
	开花普遍期	7月8日	7月12日	7月22日	7月8日	7月1日	7月14日
	开花成熟期	9月28日	9月24日	9月22日	9月18日	9月18日	9月24日
	种植密度(株/m2)	43.5	35.76	27.56	36.9	40.6	52.36
	产量(kg/hm2)	1350	1500	1240	1500	750	1200

新窗口打开|下载CSV

3 方法

3.1 冷害年的确定

确定历史冷害年份是进行冷害损失评估的第一步。在本文中,参考潘铁夫等^[18]提出的标准：规定大豆关键生育期5—9月有效积温（Growing Degree Days, GDD）的总和低于历史有效积温均值100 ℃且该年的大豆年产量低于前1个无灾年的年份被认定为冷害年。本文不对冷害等级做区分,只是判断是否发生低温冷害。分析图2可见,1988年到2018年,积温距平值小于-100 ℃的年份有1989年、1991年、1995年、2003年、2009年和2018年,这与张胜利2017年识别的鄂伦春自治旗的冷害年基本一致^[19]。但实际1991年鄂伦春大豆产量并没有造成损失,因此本文只选取1989年、1995年、2003年、2009年和2018年为鄂伦春地区低温冷害年。同时选取前一年（1988年、1994年、2002年、2008年和2017年）作为无灾年进行对比分析。

图2

新窗口打开|下载原图ZIP|生成PPT
图2东北地区低温冷害年

Fig. 2The cold years in Northeast China

3.2 模型的本地化

模型CROPGRO-Soybean内嵌在DSSAT 4.7中,已被广泛应用于世界各地的大豆生长发育模拟^[20,21]。在CROPGRO-Soybean模型中,15个遗传参数（表3）可以确定1个品种。本文利用田块测量的2014—2017年的资料,首先初始化遗传参数,将迭代次数设定为6000次,基于DSSAT模型自带的GLUE（Generalized Likelihood Uncertainty Estimation）模块计算大豆品种参数;其中2014—2015年实测数据用于模型的校准,2016—2017年作为检验年份。依次模拟研究区6个田块的品种遗传参数,最后采用均方根误差（Root Mean Square Error, RMSE）和相对均方根误差（Relative Root Mean Square Error, RRMSE）对关键生育期（开花期和成熟期）和产量的模拟值与实测值的偏差进行估计。

Tab. 3
表3
表3CROPGRO-Soybean 模型参数和田块矫正前后的遗传参数
Tab. 3CROPGRO-Soybean model parameters and the genetic coefficients of crop fields

参数	定义	默认值	五岔沟	兴盛	春亭阁	春林	乃木河	小库莫
CSDL	临界短日长(h)	12.15	14.03	12.23	13.06	14.58	11.9	14.03
PPSEN	生长发育光周期响应斜率(/h)	0.2	0.235	0.294	0.287	0.229	0.146	0.304
EM-FL	出苗到开花期的时间间隔(d）	21	13.09	18.33	15.46	16.53	26.14	22.74
FL-SH	第一次开花到第一次结荚时间间隔(d)	6	6	6	6	6	6	6
FL-SD	第一次开花到第一次结实时间间隔(d）	12	19.82	21.06	18.42	18.39	11.26	12.06
SD-PM	第一次播种到成熟的时间间隔(d)	26	37.56	24.5	36.63	31.27	22.25	34.33
FL-LF	第一次开花到叶片停止生长时间间隔(d)	20	20	20	20	20	20	20
LFMAX	最大叶片光合速率(mg CO2/(m2 s))	1.03	1.023	1.052	1.034	1.011	1.196	1.257
SLAVR	在标准生长条件下, 品种比叶面积(cm2/g)	385	311.2	303.8	301	317.6	337.9	301.1
SIZLF	整个复叶的最大叶面积(cm2)	137	138.1	145.1	138	141.2	188.5	217.6
XFRT	每日生长中分配给种子和外壳最大比例	1	1	1	1	1	1	1
WTPSD	单粒种子的最大重量(g)	0.155	0.162	0.157	0.161	0.181	0.195	0.186
SFDUR	所有荚内籽粒灌浆持续天数(d)	22	25.42	24.88	24.22	25.36	21.62	21.94
SDPDV	标准生长条件下每荚平均粒数(荚数)	2.2	2.415	2.276	2.09	2.24	2.397	1.794
PODUR	最佳生长条件下, 到结痂期持续天数(d)	13	13	13	13	13	13	13

新窗口打开|下载CSV

（1）RMSE=1n∑i=1n(Si-Oi)2
（2）RRMSE=1n∑i=1n(Si-Oi)2Oavg
式中：O_i和S_i分别为观测值和模拟值;O_avg为观测平均值;n为样本量。

3.3 产量损失估算

大豆产量损失估算主要分为5个步骤：① 基于校准后的CROPGRO-Soybean模型,模拟当地的不同种植管理和低温冷害情景并记录模拟产量,叶面积植被指数（Leaf Area Index, LAI）及其日期（Day of Year, DOY）;② 依据DOY计算每种情景下大豆生育期内≥ 10 ℃的有效积温距平值（Cold Growth Days, CDD）;③ 利用模拟的产量、CDD和关键生育期内的两期LAI组合构建冷害脆弱性模型;④ 依托GEE平台反演大豆栅格点的实际LAI及其获取日期（DOY）,并计算实际的CDD;⑤ 依据DOY逐像元匹配冷害脆弱性模型,逐像元计算产量和减产率。

3.3.1 模拟低温冷害情景 本地化校准后的作物模型能够较好地模拟气候因素对作物生长发育影响的潜能,评估低温冷害对作物的影响^[22]。本文根据当地实际记录的田间管理信息,设置不同播种日期和种植密度的组合生成一系列输入参数。依据实际田块记录数据,播种日期在5月1号至6月6号范围内等间隔（5 d）设置,分别为5月1日、6日、11日、16日、22日、27日和6月1日、6日共8种情景;种植密度在19.42~52.36株/m²间随机产生,分别为19.42株/m²、22.77株/m²、31.82株/m²、35.86株/m²、38.87株/m²、41.06株/m²,47.43株/m²和52.36株/m²共8种情景。此外,本文共设置1种正常天气情景和7种冷害情景：① 生长季内逐日降温1~3 ℃的3种冷害情景;② 分别对大豆的出苗到开花期、开花到结荚期、结荚到鼓粒期和鼓粒到成熟期,随机设置连续5 d最低温度等于0 ℃的4种冷害情景。总的来说,每个田块共有512（8 ×8 ×8）种模拟情景。最后记录512种模拟情景下所有田块生育期内的逐日LAI和对应产量。

3.3.2 计算生育期内有效积温距平值 计算每次模拟情景下大豆生育期内≥ 10 ℃的有效积温,减去1981—2018年相应生育期内的平均有效积温得到此次模拟的≥ 10 ℃积温距平值（℃）,本文选用该指标（CDD）刻画低温冷害对最终产量的影响^[23]。由于每次模拟的播种日期和天气情况不同,因此每次模拟得到的CDD也不相同。CDD计算公式如下：

（3）CDD=∑n=1nDi-GDDˉ
式中：CDD为每次模拟得到的有效积温距平值;D_i为模拟生育期第i天的≥ 10 ℃的日均温（℃）;n为每次模拟的生育期长度; GDDˉ为相同时段内1980—2018年日均温度≥ 10 ℃活动积温的平均值。

3.3.3 多情景下构建冷害脆弱性模型 以研究区大豆的普遍开花期为界限,前后各扩展40 d,形成大豆关键生育期内的早晚窗口。以所有田块多情景下（6 ×512）模型模拟的产量（Yield）为因变量,早晚窗口组合的2期LAI以及全生育期内有效积温距平值（CDD）为自变量构建回归方程,共建立1600（40 ×40）组方程。

（4）Yield=β0+β1×LAI1,d+β2×LAI2,d+β3×CDD
式中：Yield为每次模拟的产量;LAI_{1, d}为早窗口内第d天的LAI;LAI_{2, d}为晚窗口的第d天内的LAI值;CDD为每次模拟计算的有效积温距平值; β0、 β1、 β2和 β3分别为常数项、早窗口LAI_{1, d}回归系数、晚窗口LAI_{2, d}回归系数和CDD回归系数,保存生成的1600套回归系数以便后续计算。

3.3.4 冷害年逐格点反演实际LAI和估算减产率 研究表明宽范围动态植被指数（Wide Dynamic Range Vegetation Index, WDRVI）对大豆的LAI很敏感,能够很好的反演出大豆LAI^[24]。因此本文利用GEE平台逐像元提取早晚时间窗口内的最大WDRVI值及影像获取日期（DOY）,并将WDRVI转换为大豆栅格点的实际LAI^[24]。依据早、晚时间窗口内影像获取DOY的组合匹配1600组方程的回归关系系数进行逐像元产量估算,生成格点尺度的产量分布图,并计算其减产率。由于缺乏其他田块的实测数据,因此本文将格点尺度的产量和减产率取均值得到县级单产和减产率,最后利用实际记录的县级数据对研究区的减产率和灾损情况进行精度评价。减产率通过“（无灾年产量-灾害年产量）/无灾年产量×100%”表示。

（5）WDRVI=αρNIR-ρredαρNIR+ρred
（6）LAI=-ln(1.79-WDRVI)-0.532/0.3
式中： ρNIR近红外波段反射率; ρred为红光波段反射率;α是权重因子（取值0.2）,引入α是为了缩小近红外波段和红光波段对植被指数贡献之间的差距。

4 结果与分析

4.1 CROPGRO-Soybean模型评估

对比2014—2017年大豆的关键生育期和产量的实测值与模拟值,CROPGRO-Soybean可以很好地捕捉大豆的开花（Time to Anthesis as Days after Planting, ADAT）和成熟日期（Physiological Maturity Day, MDAT）（图3）,平均开花和成熟日期的RMSE分别小于5 d和6 d,RRMSE低于9.3%和7.9%。产量的RMSE为93.3 kg/hm²,RRMSE为8.0%。整体来看,成熟日期的模拟优于产量和开花日期模拟效果,但是模拟误差都在较为合理范围内（RRMSE ≤ 10%）,表明本地化后的模型能够体现鄂伦春地区大豆主要生育期和产量的变化情况。研究表明,经过校准后的CROPGRO-Soybean模型能够模拟不同天气情景和管理情景下大豆生长发育过程和产量。

图3

新窗口打开|下载原图ZIP|生成PPT
图3大豆的模拟值与实测值比较

Fig. 3Comparison between simulated and observed variables of soybean

4.2 冷害脆弱性模型结果

图4展示的是校准后CROPGRO-Soybean兴盛田块LAI值模拟结果,发现每次模拟的LAI值差异较大,这是由于在不同的天气条件和农户田间管理方式下,作物生长的过程及最终产量存在较大差异。从图4可以看出,本文选取的早晚窗口基本能够覆盖每次模拟的关键生育期,LAI的值在2~7之间变化,这些模拟结果为多元线性回归模型的训练提供了良好的代表性样本。

图4

新窗口打开|下载原图ZIP|生成PPT
图4512次CROPGRO-Soybean模型输出的LAI

Fig. 4Daily leaf area index (LAI) outputs from 512 simulations of the CROPGRO-Soybean model



早窗口和晚窗口的前部分日期组合构建的模型的决定系数呈黄色,而窗口后半部分组合的决定系数明显要高得多（呈现深红色）,但所有窗口内的日期组合建立模型的可靠性都较高（R²>0.8,图5）,尤其是早窗口的191~210 d与晚窗口所有日期的组合（R²>0.95）。该时段内植被指数更接近于最大值,所以此时期的LAI更能反映出大豆的最终产量。当然,低的回归模型训练误差并不能保证在实际数据上有良好的表现,因为许多输入数据（天气、土壤和田间管理措施等）存在一定误差的,且作物模型并不能详尽地捕捉作物生长发育的所有相关过程^[25]。但高的决定系数确实表明了在大豆关键生育期内对LAI的两次观测值与CDD的组合可以得到较为准确的估算产量。

图5

新窗口打开|下载原图ZIP|生成PPT
图5回归模型的决定系数

Fig. 5The coefficient of determination for the regression models

4.3 估算产量损失

4.3.1 基于GEE反演的实际LAI 基于GEE平台反演的2018年大豆关键生育期内早晚窗口最大LAI（图6a、6b）及其对应的影像获取DOY（图6c、6d）,早晚窗口LAI的取值范围为6~7.2,与实际模拟较为一致,进一步证明CROPGRO-Soybean模型模拟的低温冷害情景能较为准确的刻画当地低温和产量之间的关系。虽然研究区内早晚窗口内大豆种植栅格像元点的可用影像日期不尽一致,但是通过组合关键生育期的2次LAI,能够突破遥感影像时空不连续的限制,进而由点尺度的作物模型模拟推演到面尺度的产量估算,实现高精度格点尺度的灾害损失评估。

图6

新窗口打开|下载原图ZIP|生成PPT
图62018年大豆生育期内早晚窗口Sentinel-2反演的最大LAI值及其对应的日期

注：a、c分为早窗口最大LAI值及其对应的提取日期;b、d分别为晚窗口最大LAI值及其对应的提取日期。
Fig. 6The maximum LAI and its specific observation dates of early and late growing season windows obtained from Sentinel-2 in 2018



4.3.2 不同冷害情景下校准后的CROPGRO-Soybean模型模拟的产量 本文共设置7种低温冷害情景,用来研究不同冷害情景对大豆产量造成的影响。情景1~情景3分别为全生育期温度降低3°、2°和1°的实验,情景5~情景8分别为出苗到开花期、开花到结荚期、结荚到鼓粒期和鼓粒到成熟期,随机设置连续5 d最低温度为0 ℃,情景4代表正常天气。大豆生长过程中遭遇冷害低温会出现普遍的减产,模拟的大豆产量在不同的低温情景下差异较大（图7）。在全生育期温度降低越多,大豆产量越低（情景1~3 ℃）。在大豆的生长季中,逐日最低温度降低1 ℃,生育期内GDD大约减少100 ℃·d,产量分布在910.5~1242.5 kg/hm²;逐日均温降低2 ℃,生育期内GDD大约减少200 ℃·d,产量分布在571.8~838.0 kg/hm²;逐日均温降低3 ℃,生育期内GDD大约减少300 ℃·d,产量分布在273.8~435.7 kg/hm²。大豆属于喜温作物,尤其在进入花芽期后需要较高的温度,低温会抑制后续各个生育期的生长,从而导致后期鼓粒期受阻,出现减产甚至绝收的情况。从 情景5~情景8可以看出,大豆在不同生长发育期内对低温冷害的响应不同,在结荚到鼓粒期降温对产量的影响最大,这是因为大豆生长后期对温度特别敏感,鼓粒期温度低,种子发育会受到影响,进而增加秕粒并延迟成熟,最终对产量造成影响,导致减产,这与前人的研究结果一致^[26]。

图7

新窗口打开|下载原图ZIP|生成PPT
图7不同冷害情景下校准后的CROPGRO-Soybean模型模拟产量

Fig. 7Estimated yields by calibrated CROPGRO-Soybean model under different cold injury scenarios



4.3.3 估算历史冷害年份县级减产率 根据鄂伦春历史气象资料分析获得1989年、1995年、2003年、2009年和2018年5个冷害年,且这5年对应的前一年（1988年、1994年、2002年、2008年和2017年）都是正常年。由于2015年以前没有Sentinel-2A遥感卫星影像数据,只能采用空间分辨率为30 m的Landsat 5卫星影像,生成30 m的田块间产量分布图。2017年和2018年采用Sentinel-2A卫星影像,生成10 m的田块间产量分布,并计算低温冷害年的减产率。由于缺乏田块的实测数据,故将逐像元模拟格点尺度的产量和减产率取均值得到县级产量和减产率。由图8展示的低温冷害年实际的县级减产率（2018年的实际记录县级产量数据缺失）和模拟的县级减产率可以看出,历史减产率（图8）的变化与大豆生育期内有效积温减少量（图2）分布基本一致。如2003年的生育期内GDD最少,遭遇的低温冷害程度最高,实际减产率最高;2003年模型模拟的减产率高达50.5%,实际为47.7%;减产率最低的是1989年,模型模拟减产率为9.6%,而实际减产率是6.4%;除1995年模型模拟的减产率略高于实际减产率外,1989年、2003年和2009年模型模拟的减产率均略高于实际情况。所有历史冷害年减产率模拟的误差都在一倍方差内,表明该方法能够较为准确地评估由于低温冷害造成的大豆产量损失状况。

图8

新窗口打开|下载原图ZIP|生成PPT
图8实际减产率与模拟减产率的对比

Fig. 8The comparison between actual and simulated yield losses



4.3.4 估算2017年和2018年产量的空间分布 基于上述历史冷害年损失的精度验证结果,发现本文提出的方法能够较为准确的评估研究区低温冷害导致的产量损失状况。因此利用该方法生成研究区最近1次冷害年（2018年）和无灾年（2017年）的大豆产量空间分布图（图9）。取像元平均值得到县级产量值,2017年和2018年鄂伦春县级估测产量分别为2004.6 kg/hm²和1540 kg/hm²。由图9可以看出,冷害年（2018）区域产量估测值整体上明显低于无灾年（2017）,这与实际情况相符。从图9①~图9⑧产量模拟空间分布的局部放大图可以看出,遭受相同的冷害,不同田块对低温冷害天气的减产反应不一样,图9①~图9②和图9⑤~图9⑥可以看出2018年遭受冷害后,对应田块的产量明显低于2017年,但也有少量区域（图9③和图9⑦）2年的产量基本持平,甚至有极少的田块图9④和图9⑧可以明显看出产量有轻微增加的现象。通过实地调研考察发现,这是由于不同农户的田间管理措施、灌溉条件和地形因素的差异导致的。

图9

新窗口打开|下载原图ZIP|生成PPT
图9无灾年(2017年)和冷害年(2018年)产量空间估测

注：①~④分别表示2017年局部放大区域的地块产量估测,⑤~⑧表示2018年相应位置的地块产量估测。
Fig. 9Spatial distribution of estimated yields in normal year (2017) and cold year (2018)



4.3.5 2018年的冷害减产率评估 从2018年低温冷害年份减产率的空间分布特征可以看出（图10）,相较于2017无灾年,2018冷害年鄂伦春产量有明显的减少。整体来看,2018年主要表现为大面积轻微减产,基本覆盖整个研究区,重度减产主要分布在研究区的东北地区,无灾区主要零星的分布在研究区南部地区。鄂伦春2018年减产率达到23.1%。通过统计乡镇单元低温冷害导致的减产率（表4）,发现约98%以上的乡镇表现为减产,除了二十里村未受灾,减产率为-29.94%。表明2018年低温冷害事件对鄂伦春的大豆整体上造成较为严重的影响。

图10

新窗口打开|下载原图ZIP|生成PPT
图102018年相对2017年的减产率的空间分布

注：A和B分别表示2018年鄂伦春局部放大区域的田块减产率。
Fig. 10The spatial distribution of yield losses in 2018 relative to 2017



Tab. 4
表4
表4大豆乡镇减产情况统计
Tab. 4The estimated yield losses of soybean at village level

乡镇名	减产率(%)	乡镇名	减产率(%)	乡镇名	减产率(%)
卡日楚村	21.23	卧北村	29.77	向阳村	24.35
岭南村	14.64	扎西村	19.54	乌尔其村	31.14
讷尔克气村	22.61	五岔沟村	15.86	朝阳猎民村	21.19
青松沟村	24.52	乃木河村	25.15	十六栋房村	32.28
团结村	19.38	马尾山村	18.34	诺敏河村	6.54
渔场村	8.76	库维地村	22.58	小二红村	18.08
新兴村	9.55	跃进村	28.83	龙头村	29.25
欧肯河村	18.24	奎勒河村	5.29	新丰村	20.63
兴盛村	28.04	东升村	21.52	卧罗河村	30.79
二十里村	-29.94	小库莫村	17.63	铁东村	30.05
乌鲁布铁猎民村	9.53	都柿沟村	21.30	大库莫村	21.58
红旗村	22.65	朝阳沟村	22.67	二根河村	7.97
马鞍山村	18.52	新发村	26.92	东升村	21.52
春林村	26.02	春亭阁村	16.09	毛家铺村	18.77

新窗口打开|下载CSV

5 结论与讨论

本文基于GEE平台和CROPGRO-Soybean作物模型在10 m像元尺度评估了鄂伦春大豆冷害损失,实现了较高精度的产量制图。发现校准后的CROPGRO-Soybean模型能够模拟大豆对不同气候情景和田间管理措施的响应,大豆在整个生育期内持续遭受冷害（1~3 ℃）造成的减产比单一生育期（4个生育期随机生成连续5 d温度为0 ℃）更加严重。对比大豆县级产量的模拟值与实际值,得出该方法能够有效的模拟大豆产量;且估算的冷害年县级减产率与实际情况基本一致,生成的格点尺度的减产率分布图基本上可以表征田块间产量损失的空间差异性。相比传统的站点作物模型和遥感估测灾害损失,该方法更加高效,限制少,空间泛化能力强,可应用于大尺度以及高精度田块间的大豆产量估算和灾损评估。该研究为大豆产量以及灾害损失评估的业务化运行提供了新思路。

然而,由于缺少田块实际产量损失调研数据,使得我们必须将像元尺度的产量整合到县级进行模型结果评估。今后应收集更多的田块数据进行结果验证。为了实现精细制图,本文以Sentinel-2A影像为数据源,然而Sentinel系列卫星在2015年开始运行,导致2015年以前的灾害只能借用30 m的Landsat 5系列卫星数据制图分析。因此,未来研究可着重于多源遥感数据融合和多种作物模型耦合的产量预报和损失评估,进一步提高该方法估产和损失评估的精度。

参考文献 View Option 原文顺序
文献年度倒序
文中引用次数倒序
被引期刊影响因子

[1] Das A, Parida S K. Advances in biotechnological applications in three important food legumes
Plant Biotechnology Reports, 2014,8(2):83-99.
[本文引用: 1]

[2] Li Wenbing, Song Chunxiao, Chang Xingchao, et al. Drought resistance evaluation of 20 soybean varieties under drought stress
Journal of Northeast Agricultural University, 2019,50(4):1-10.
[本文引用: 1]

[ 李文滨, 宋春晓, 苌兴超, 等. 干旱胁迫下20个大豆品种抗旱性评价
东北农业大学学报, 2019,50(4):1-10.]
[本文引用: 1]

[3] Liu X, Jin J, Wang G, et al. Soybean yield physiology and development of high-yielding practices in Northeast China
Field Crops Research, 2008,105(3):157-171.
[本文引用: 1]

[4] Liu Dan, Yu Chenglong, Du Chunying. Dynamic monitoring rice delayed type chilling damage based on remote sensing in northeast China
Transactions of the Chinese Society of Agricultural Engineering, 2016,32(15):157-164.
[本文引用: 1]

[ 刘丹, 于成龙, 杜春英. 基于遥感的东北地区水稻延迟型冷害动态监测
农业工程学报, 2016,32(15):157-164.]
[本文引用: 1]

[5] Hao T Y, Wang S G, Shang K Z, et al. Research progress and outlook of low temperature chilling injury in Northeast China
Meteorological and Environmental Research, 2010,11(3):85-91.
[本文引用: 1]

[6] Han Rongqing, Li Weijing, Ai Wanxiu, et al. The climatic variability and influence of first frost dates in Northern China
Acta Geographica Sinica, 2010,65(5):525-532.
[本文引用: 1]

[ 韩荣青, 李维京, 艾婉秀, 等. 中国北方初霜冻日期变化及其对农业的影响
地理学报, 2010,65(5):525-532.]
[本文引用: 1]

[7] Zhang Jianping, Wang Chunyi, Zhao Yanxia, et al. Impact evaluation of low temperature to yields of maize in Northeast China based on crop growth model
Acta Ecologica Sinica, 2012,32(13):4132-4138.
DOI:10.5846/stxbURL [本文引用: 1]

[ 张建平, 王春乙, 赵艳霞, 等. 基于作物模型的低温冷害对中国东北三省玉米产量影响评估
生态学报, 2012,32(13):4132-4138.]
[本文引用: 1]

[8] Zhou Liwei. Spatiotemporal characteristics of maize low temperature and cold damage and its effect on maize yield in Heilongjiang Province
[D]. Harbin: Northeast Agricultural University, 2017.
[本文引用: 1]

[ 周立威. 黑龙江省玉米低温冷害时空特征及其对玉米产量的影响
[D]. 哈尔滨: 东北农业大学, 2017.]
[本文引用: 1]

[9] Liu Xiaofei, Zhang Zhao, Shuai Jiabing, et al. Effect of chilling injury on rice yield in Heilongjiang Province
Acta Geographica Sinica, 2012,67(9):1223-1232.
[本文引用: 1]

[ 刘晓菲, 张朝, 帅嘉冰, 等. 黑龙江省冷害对水稻产量的影响
地理学报, 2012,67(9):1223-1232.]
[本文引用: 1]

[10] Sinclair T R, Rawlins S L. Inter-seasonal variation in soybean and maize yields under global environmental change
Agronomy Journal, 1993,85(2):406-409.
[本文引用: 1]

[11] Li Xiufeng, Ma Shuqing, Zhao Huiying, et al. Evaluation on maize chilling injury based on WOFOST model in Hetao irrigation region in Inner Mongolia
Chinese Journal of Agrometeorology, 2016,37(3):352-360.
[本文引用: 1]

[ 李秀芬, 马树庆, 赵慧颖, 等. 基于WOFOST模型的内蒙古河套灌区玉米低温冷害评价
中国农业气象, 2016,37(3):352-360.]
[本文引用: 1]

[12] Chen De. Study on monitoring and evaluating the chilling injury of rice in Northeast China by remote sensing and crop model
[D]. Nanjing: Nanjing University of Information Technology, 2017.
[本文引用: 1]

[ 陈德. 基于遥感和作物模型的东北水稻低温冷害监测评估
[D]. 南京: 南京信息工程学院, 2017.]
[本文引用: 1]

[13] Gorelick N, Hancher M, Dixon M, et al. Google Earth Engine: Planetary-scale geospatial analysis for everyone
Remote Sensing of Environment, 2017,202(3):18-27.
[本文引用: 1]

[14] Boote K J, Jones J W, Batchelor W D, et al. Genetic coefficients in the CROPGRO-soybean model
Agronomy Journal, 2003,95(1):32-51.
[本文引用: 1]

[15] Chen Zhongxin, Ren Jianqiang, Tang Huajun, et al. Progress and perspectives on agricultural remote sensing research and applications in China
Journal of Remote Sensing, 2016,20(5):748-767.
[本文引用: 1]

[ 陈仲新, 任建强, 唐华俊, 等. 农业遥感研究应用进展与展望
遥感学报, 2016,20(5):748-767.]
[本文引用: 1]

[16] Angstrom A. Solar and terrestrial radiation
Quarterly Journal of the Royal Meteorological Society, 1924,50(210):121-126.
[本文引用: 1]

[17] Jones J W, Hoogenboom G, Porter C H, et al. The DSSAT cropping system model
European Journal of Agronomy, 2003,18(34):235-265.
[本文引用: 1]

[18] Pang Tiefu, Fang Zhanshen, Zhao Hongkai, et al. Crop Low Temperature Damage and Its Defense. Beijing: Agriculture Press, 1983: 20-25.
[本文引用: 1]

[ 潘铁夫, 方展森, 赵洪凯, 等. 农作物低温冷害及其防御. 北京: 农业出版社, 1983: 20-25.]
[本文引用: 1]

[19] Zhang Shengli. Analysis on the characteristics of chilling injury in Oroqen
China Southern Agricultural Machinery, 2017,48(4):181-182.
[本文引用: 1]

[ 张胜利. 鄂伦春自治旗低温冷害特点分析
南方农机, 2017,48(4):181-182.]
[本文引用: 1]

[20] Wang Xinli, Wang Futang, Qiu Guowang. Application of optimizing theory to determining genetic parameters involved in CERES-Soybean model
Journal of Applied Meteorology, 1995,6(a01):49-54.
[本文引用: 1]

[ 王信理, 王馥棠, 裘国旺. 系统优化理论在CERES—大豆模式品种参数求解中的应用
应用气象学报, 1995,6(a01):49-54.]
[本文引用: 1]

[21] Zhu Dawei, Jin Zhiqing. Impacts of changes in both climate and its variability on food production in Northeast China
Acta Agronomica Sinica, 2008,34(9):1588-1597.
[本文引用: 1]

[ 朱大威, 金之庆. 气候及其变率变化对东北地区粮食生产的影响
作物学报, 2008,34(9):1588-1597.]
[本文引用: 1]

[22] Pan Zhuokun. Integration of remote sensing and crop growth model for regional low temperature impact monitoring early warning and yield estimation
[D]. Hangzhou: Zhejiang University, 2015.
[本文引用: 1]

[ 潘灼坤. 耦合遥感信息与作物生长模型的区域低温影响监测
[D]. 杭州: 浙江大学, 2015.]
[本文引用: 1]

[23] Ma Shuqing, Xi Zhuxiang, Wang Qi. Risk evaluation of cold damage to corn in Northeast China
Journal of Natural Disasters, 2003,12(3):137-141.
[本文引用: 1]

[ 马树庆, 袭祝香, 王琪. 中国东北地区玉米低温冷害风险评估研究
自然灾害学报, 2003,12(3):137-141.]
[本文引用: 1]

[24] Anthony N R, Anatoly G, Yi P, et al. Green leaf area index estimation in maize and soybean: Combining vegetation indices to achieve maximal sensitivity
Agronomy Journal, 2012,104(5):1336-1337.
[本文引用: 2]

[25] Lobell D B, Thau D, Seifert C, et al. A scalable satellite-based crop yield mapper
Remote Sensing of Environment, 2015,164(11):324-333.
[本文引用: 1]

[26] Sang Shupeng. Study on measures to cope with low temperature damage in different growth stages of soybean
Soybean Science & Technology, 2013,5(1):53-54.
[本文引用: 1]

[ 桑树鹏. 大豆不同生育期内应对低温冷害措施的研究
大豆科技, 2013,5(1):53-54.]
[本文引用: 1]