杨 磊 韩丽娟 宋金玲 李 森

1)(北京师范大学遥感科学国家重点实验室，北京 100875) 2)(北京师范大学地理科学学部遥感科学与工程研究院，北京 100875) 3)(国家气象中心，北京 100081)

引 言

在全球气候变化背景下，气候数据统计表明：中国气候自20世纪80年代开始出现明显变暖的趋势，且北方增温更甚于南方[1-3]，北方夏季出现极端高温的现象也在加剧，导致夏玉米产量形成关键期遭遇高温热害的风险不断加大[4]。因此，对夏玉米高温热害监测评估方法的研究尤为重要。

目前对夏玉米高温热害监测评估研究多采用地面气象站点数据，并结合夏玉米发育期数据，分析玉米抽雄吐丝期出现高温热害的次数及空间分布[5-6]。由于地形以及气象站点分布不均，学者们多采用空间插值方法获取无测站分布地区的温度数据[7-9]。刘哲等[10-11]利用空间插值的气象数据，由概率统计方法得到黄淮海地区夏玉米出现高温热害的概率及其空间分布规律。采用遥感数据进行高温热害监测评估研究工作的优势在于实时性和区域性，对研究大范围区域具有较好的代表性与准确性[7]。

利用遥感数据对水稻等农作物的高温监测评估已有较多研究[12-17]，这些研究工作证明地表温度与气温之间具有较强的线性相关关系[18]。本文主要利用MODIS(moderate-resolution imaging spectroradiometer)地表温度产品，建立地表温度与气温的回归关系，采用插值后的气象站观测气温数据对遥感数据缺失区域进行填补，得到卫星数据与地面实测数据融合的长时间气温数据序列，结合夏玉米高温等级热害指标，对2008—2018年黄淮海夏玉米主产区的高温热害进行分析评估。

1 研究区概况

本研究主要集中在黄淮海地区的京津冀区域、河南省和山东省，这里是夏玉米的主要种植区，年降水量为500～900 mm，7—8月降水量约占全年降水量的45%～65%[19]，耕地面积占全国总耕地面积的25%，主要栽培方式为冬小麦和夏玉米轮作，热量和水分都是制约该区域农业发展的因子[20]。

2 数据及预处理

2.1 遥感数据

2.1.1 MODIS数据

本研究所用到的MODIS数据包含MOD09A1地表反射率产品以及MOD11A1/MYD11A1 MODIS地表温度产品。其中地表反射率产品空间分辨率为500 m×500 m，时间分辨率为8 d；地表温度产品分别由Terra和Aqua两颗卫星搭载的MODIS传感器生成，空间分辨率为1 km×1 km，时间分辨率为1 d，该数据通过分裂窗算法反演得到，可获取白天和夜间的地表温度[21]。由于两种地表温度产品的共同使用，每个像元每日可得到4次地表温度的观测数据。为方便研究需将MODIS地表温度产品的像元值转化为以摄氏度为单位的地表温度[22]。MODIS产品均来源于美国太空总署(NASA)遥感数据共享网站(https:∥ladsweb.modaps.eosdis.nasa.gov)，数据时间为2008—2018年每年7—9月，空间范围覆盖京津冀区域、河南省和山东省。

利用MOD09A1地表反射率产品计算夏玉米生长期内的植被指数，通过植被指数的变化从遥感图像上筛选出夏玉米像元作为种植区域；在此基础上，结合MOD11A1/MYD11A1地表温度产品对夏玉米种植区域的地表温度变化进行分析。对MODIS数据的预处理包括影像拼接、投影转换、波段提取、重采样、影像裁剪和去云处理等，该过程主要通过NASA官方提供的MRT(MODIS Reprojection Tool)工具批量处理实现，并对提取出的数据进行滤波重建，消除噪声的干扰。

2.1.2 Landsat数据

Landsat数据空间分辨率(30 m×30 m)较MODIS数据空间分辨率更高,本研究利用与MODIS数据同期的Landsat数据进行辅助研究，主要用于验证夏玉米提取结果的精度。该数据来源于美国USGS网站(https:∥earthexplorer.usgs.gov)，数据时间和空间覆盖范围与MODIS数据相同，对Landsat数据处理包括几何校正、大气校正等。

2.1.3 Google Earth高清影像数据

Google Earth的卫星影像是卫星影像与航拍数据的整合。其中卫星影像来自于快鸟(QuickBird)商业卫星、陆地卫星Landsat7，IKONOS和SPOT5等。本研究利用Google Earth高清影像数据辅助验证夏玉米种植区域提取的位置精度，并帮助识别具体的空间细节信息。

2.1.4 DEM高程数据

本研究借助SRTM(Shuttle Radar Topography Mission，航天飞机雷达地形测绘使命)地形产品数据对研究区域的平原与山区进行划分。该产品由美国NASA和美国国防部国家测绘局(NIMA)联合测量，空间分辨率为30 m×30 m和90 m×90 m两种。本文选用90 m×90 m空间分辨率的DEM数据进行地形划分，通过像元聚合将空间尺度上升到与MODIS数据相同，并根据地理学相关定义，将高程不低于200 m的区域定义为山区，低于200 m的区域定义为平原。该数据来源于地理空间数据云遥感数据平台(http:∥www.gscloud.cn/search)。

2.2 气象数据

气象数据为研究区域中气象站的气温观测数据，主要用于与遥感数据结合，分析夏玉米主产区历年温度的变化情况。该数据来源于国家气象中心，包括2008—2018年每年7—9月逐日最高气温和平均气温。在京津冀区域、河南省和山东省共有353个气象站，其中平原297个站，山区56个站。

2.3 夏玉米发育期及面积验证数据

夏玉米发育期数据来源于国家气象中心，时间为2008—2018年7—9月；面积验证数据包括2008—2018年逐年夏玉米种植面积及受灾面积，种植面积数据来源于国家统计局发布的统计年鉴；受灾面积数据来自国家气象中心农业气象情报产品和期刊新闻等。这些统计数据用于对模拟结果的定量验证分析。

3 研究方法

利用MOD09A1地表反射率产品进行土地覆盖分类，提取出遥感图像上的夏玉米像元，并结合发育期数据，逐像元判别其是否正处于高温热害的敏感时期。然后利用MOD11A1/MYD11A1地表温度产品反演逐日最高气温和平均气温[18]。考虑到温度反演受云层遮挡和气溶胶的影响，本文采用实测气温数据进行插值填补，以弥补遥感数据的缺失[5-6]。最后基于融合后的气温数据，结合夏玉米高温热害指标，对高温热害进行等级划分和对受灾面积进行统计，实现对夏玉米高温热害的监测评估[23]。

3.1 夏玉米种植区域提取

3.1.1 提取方法

结合土地利用分类筛选出耕地类别的像元，然后根据夏玉米生长发育的物候特征，选择7月上中旬的MODIS地表反射率影像，黄淮海平原此时正处于夏玉米播种初期，具有较小的归一化差值植被指数(normalized difference vegetation index，NDVI)值，故将NDVI值大于0.5的区域剔除；选择9月上中旬的影像，提取NDVI值大于0.7的区域，两者进行合并作为夏玉米的种植区域[24-29]。由于两期数据所提供信息的有限性，为充分利用已有数据，需在上述基础上进行二次提取。本研究利用Google Earth高清影像，通过目视选择夏玉米像元作为参考像元，提取其NDVI时间序列并进行Savitzky-Golay滤波(简称S-G滤波)处理；再将第1次分类得到的夏玉米像元进行相同的提取NDVI时间序列操作和滤波处理，与夏玉米参考像元的NDVI时间序列进行相似度比较，相似度判定指标为欧氏距离[29](阈值定为0.7)，将欧氏距离值在0～0.7范围内的像元提取出来作为最终的种植区提取结果。

3.1.2 精度验证

采用Landsat数据以及Google Earth高清影像数据进行对比验证，检验利用MODIS数据提取结果在不同地区的适用性，验证内容包括种植面积精度及位置精度的验证。

面积精度的验证借助Landsat数据，利用上述相同的分类方法，将分类结果与MODIS数据提取结果进行对比，计算二者相对误差并判断能否满足研究需求；位置精度验证则是借助Google Earth高清遥感影像数据，将MODIS数据分类结果与Google Earth中高清影像数据进行图层的空间叠置，根据高分辨率影像判读其位置误差。

3.2 温度转换模型

本文建立MODIS地表温度与气温之间的转化关系，将气象站观测的气温数据进行空间插值，并基于地表温度和气温之间的线性回归关系，填补MODIS地表温度数据中的缺失区域，得到最终的气温融合数据。

针对MODIS的数据缺失情况，通常依据地表温度具有空间连续性特征，根据气象站实测气温与像元温度之间的函数关系，估算云覆盖下的地表温度[8]。该方法操作简单且精度较高[30-31]，具有较好的应用前景，因此本研究主要采用该方法对数据缺失区域进行填补。

3.2.1 变量相关性检验

将气象站数据与MODIS地表温度产品进行空间位置叠加，从每日4个时次的地表温度数据中提取出气象站对应的像元值，将有效数据按4:1的比例划分为训练集和测试集，从而保证每年都有数据用于模型构建，即利用2008—2018年80%的数据进行建模；结合多元逐步回归的方法，将日平均气温(Yave)和日最高气温(Ymax)作为因变量，TOD，TON，TYD和TYN4个变量(对应每日4次地表温度观测值)作为自变量建立线性模型，其中，TOD为MOD11A1地表温度产品中白天温度的观测值，TON为MOD11A1地表温度产品中夜间温度的观测值；TYD为MYD11A1地表温度产品中白天温度的观测值，TYN为MYD11A1地表温度产品中夜间温度的观测值。对自变量之间相关性进行检验发现，TON和TYN(每日观测的两次夜间温度)间具有显著相关性，相关系数为0.965，达到0.01显著性水平，因此在后续的建模处理过程中对自变量进行降维处理，以消除模型的多重共线性。

3.2.2 回归模型建立

逐步回归方法的优势为模型自变量是逐步引入模型中，每次只引入1个新的变量，同时对原变量进行显著性检验，如果原变量不再显著，就将其剔除，这样保证每个引入的新变量对模型因变量的影响都是显著的。

首先利用主成分分析方法对上述TOD，TON，TYD和TYN4个变量进行降维处理，将提取到的主成分作为模型中新的自变量进行逐步回归，主成分中各个变量是原始变量经过以下标准化处理得到：

Ti″=(Ti-Tmin)/(Tmax-Tmin)。

(1)

其中，Ti″表示经过标准化处理后的原始变量，Ti表示原始变量，Tmin表示变量在研究区内对应的最低温度累计频率0.05%处的温度值，Tmax表示变量在研究区内对应的最高温度累计频率95.5%处的温度值。

将主成分还原为初始变量，得到最终线性回归模型。在处理多重共线性问题时，由于建立的模型仅用于预测，若拟合程度已达到较好效果，轻微的多重共线则不会严重影响预测的结果，因此采用主成分分析法进行1次降维即可。

3.2.3 模型精度验证

利用2008—2018年7—9月剩余的20%数据进行模型的精度验证，通过比较模型反演结果与气象站实测结果的均方根误差和决定系数描述模型精度。同时，考虑到研究区域的地形以平原和山区为主，模型误差在山区和平原会表现不同，本研究结合DEM高程数据，分别统计这两种地形下，对应像元的模型模拟值与实测值之间的误差。

3.3 夏玉米高温热害等级指标

通过研究温度对玉米产量与品质的影响，结合常用的业务服务等级指标，建立夏玉米高温热害等级指标[32-34](表1)，根据不同等级下日平均气温、日最高气温临界阈值及持续日数，可判断夏玉米高温热害等级。

表1 夏玉米高温热害等级指标

4 结果与分析

4.1 夏玉米种植区

4.1.1 夏玉米种植区域分布

采用MODIS数据进行夏玉米种植区提取时，考虑到MODIS数据的时间跨度较大，2008—2018年夏玉米种植区域可能发生变化，因此将每年的分类结果进行合并作为夏玉米种植分布的平均状况。结果如图1所示，图中黄绿色部分代表夏玉米种植区域。

图1 2008—2018年研究区域夏玉米分布平均状况

4.1.2 夏玉米种植面积精度验证

采用MODIS和Landsat数据提取的夏玉米样区(图2)，验证夏玉米种植面积提取精度，结果表明：尽管二者的提取结果在细节上有差别，但总体分布情况吻合程度较高。从MODIS数据中提取的夏玉米种植面积为2.27×105hm2，从Landsat数据中提取的夏玉米种植面积为2.46×105hm2，二者的相对偏差为8.27%，满足研究精度需求。

图2 MODIS与Landsat数据种植区提取结果

用河南省和山东省的夏玉米种植面积进行定量验证，通过查询国家统计年鉴数据，2017年山东省玉米的种植面积为4.00×106hm2，河南省玉米的种植面积为3.99×106hm2；利用遥感数据提取的结果分别为4.50×106hm2和4.36×106hm2，相对偏差分别为12.5%和9.0%。

4.1.3 夏玉米分布位置精度验证

利用Google Earth高清影像数据与MODIS数据提取结果进行对比，结果显示二者提取结果的边界轮廓十分吻合，仅在西部山区存在一些错分的像元和东部存在一些漏分的像元，总体分类结果具有较高的位置精度。

4.2 气温与地表温度转换模型与验证

(2)

(3)

(4)

将X1，X2和X3作为新的自变量进行逐步回归，并还原为初始变量后，得到最终的线性回归模型为

Yave=0.085TOD+0.444TON+0.176TYD+

0.323TYN+0.299。

(5)

同理，地面实测的最高温度的线性模型为

Ymax=0.177TOD+0.152TON+0.223TYD+

0.281TYN+8.380。

(6)

对该模型进行变量显著性检验以及模型共线性检验，发现模型均达到极显著水平，方差膨胀因子值均为1.000，已不存在多重共线性，其决定系数为0.8171，均方根误差为2.08℃。当分别考虑该模型在平原和山区的表现时，平原区域模型拟合的决定系数为0.8202，均方根误差为1.28℃；山区模型拟合的决定系数为0.8109，均方根误差为2.32℃，二者均达到0.001显著性水平(图3)。综上所述，该模型具有较高精度，能够准确反映当前温度状态，在平原区域表现相对更好，且误差精度与已有研究结果吻合[6]，整个研究区域的模型结果见表2。

图3 平原和山区模拟气温与实测气温散点分布

表2 气温拟合模型

4.3 2008—2018年夏玉米高温热害受灾面积

通过对2008—2018年有效数据进行模型拟合得到夏玉米主产区高温热害的面积变化(表3)，其中2011年和2012年夏玉米主产区未遭受明显的高温热害。

表3展示了2008—2018年历史遭受高温热害的面积变化，2009年、2013年、2014年、2017年和2018年均遭受较大范围的高温热害，其中以2017年和2018年的受灾情况最为明显，在这些受灾区域中，河南省的受灾程度最为严重，京津冀区域次之，山东省遭受高温热害程度较轻。此外，在这两年间亦出现了历史上少有的三级高温热害情况，面积分别为1.86×105hm2和8.98×105hm2，受灾区域集中在河南省南部和西北部，以及山东省西部等。

由于夏玉米实际遭受高温热害的面积较难统计，主要通过相关统计数据对本研究结果进行验证[35-38]。2017年7月下旬至8月上旬河南省出现大范围连续10～15 d超过35℃的高温天气；2018年7月下旬至8月上旬，山东、河北、河南3省高温日数和平均最高气温均达1981—2018年历史同期最高值，高温时段与夏玉米抽雄-吐丝期吻合，夏玉米遭受不同程度的高温热害，导致结实率下降。从统计得到的高温热害发生强度和减产率，跟受灾面积互相印证了本文结果的合理性。

4.4 夏玉米高温热害空间分布

上述研究表明：夏玉米主产区遭受高温热害的面积增加，结合表3与近些年数据模型拟合的结果，分析得出受灾范围主要集中在研究区域的中部，即京津冀区域南部，河南省北部以及山东省西部地区。此外，2017年和2018年夏玉米主产区遭受大范围不同程度的高温热害。鉴于此，本研究以2017年和2018年为重点研究年份，对夏玉米主产区高温热害及其等级的空间分布进行分析(图4)，并对受灾面积进行统计，结果表明：2017年夏玉米主产区遭受高温热害面积为6.68×106hm2，主要集中在河北省南部，河南省大部及山东省西部部分区域，其中一级受灾面积为4.25×106hm2，二级受灾面积为2.25×106hm2，三级受灾面积为1.86×105hm2，主要集中在河南省。2018年受灾面积有所增加，达到1.01×107hm2，主要集中在河北省东部和南部，河南省东部以及山东省中部地区，其中一级受灾面积为4.08×106hm2，二级受灾面积为5.07×106hm2，三级受灾面积为8.98×105hm2。

图4 2017年和2018年夏玉米主产区高温热害等级分布

5 结论与讨论

本文基于遥感数据和气象站数据，建立了夏玉米高温热害监测评估模型。利用该模型对2008—2018年夏玉米主产区的高温热害进行评估，得到以下主要结论：

1)地表温度与气温具有显著的线性相关关系，平原的日平均气温模拟的决定系数在0.8以上，日最高气温模拟的决定系数在0.7以上；山区精度略低，日平均气温模拟的决定系数在0.7以上，日最高气温模拟的决定系数在0.6以上，二者均达到0.001 显著性水平，且日平均气温的反演精度高于日最高气温的反演精度。

2)2008—2018年黄淮海夏玉米主产区遭受高温热害的面积呈增加趋势，且空间分布具有相似性，受灾范围主要集中在京津冀区域南部、河南省北部以及山东省西部。

3)2017年和2018年黄淮海夏玉米主产区遭受较严重的大范围高温热害。2017年夏玉米主产区高温热害主要集中在河北省南部，河南省大部以及山东省西部的部分区域；2018年的受灾面积有所增加，集中在河北省东部和南部、河南省东部以及山东省中部地区。

本研究的局限性体现在以下几个方面：①在遥感数据及图像分类方法方面，本研究基于500 m×500 m空间分辨率的MODIS地表反射率数据，根据夏玉米生长发育的物候特征提取夏玉米的种植区域，提取结果的精度还需进一步提升。如使用Sentinel数据、国产高分卫星数据等具有更高空间分辨率的遥感影像，结合更加准确的分类算法能够有效提高分类结果的精度。②在建模方法方面，本研究通过多元逐步回归并结合主成分分析，模拟出2008—2018年7—9月黄淮海夏玉米主产区的日平均气温和日最高气温，气温模拟的精度直接影响最终结果。针对建模过程中出现的多重共线性问题，可采用更加简易且准确的方法进行消除，如利用随机森林算法进行建模，该方法是一种有监督的学习算法，以决策树为“基学习器”，通过随机且有放回的采样构建训练集，按照多棵树预测结果的平均值作为最终结果。由于基于决策树模型的分类机制，使用随机森林算法建模不会出现多重共线性的问题，且可以解决过拟合问题，该方法的建模过程更简单且具有较高精度。③在影响夏玉米生长发育的气候条件方面，本研究仅考虑了气温对夏玉米生长发育的影响，而在实际情况中，多种气象因素共同影响夏玉米的长势及产量，如降水可在一定程度上缓解因温度过高导致的玉米植株损伤，因此在模型中引入降水等其他气象因素将有助于提高热害风险区域评估的准确度。