基于FY-4 AGRI与Himawari-8 AHI的干旱灾害监测及旱情分析<br/>——以2018年河北省秋季干旱为例*

基于FY-4 AGRI与Himawari-8 AHI的干旱灾害监测及旱情分析
——以2018年河北省秋季干旱为例*

2019-10-14李亚男于文金周鸿渐

灾害学 2019年4期

李亚男，于文金，谢涛，周鸿渐

(1.南京信息工程大学地理科学学院，江苏南京 210044；2.南京信息工程大学遥感与测绘工程学院，江苏南京 210044；3.南京信息工程大学教务处，江苏南京 210044)

在我国，由于季风气候的特性，降水时空分布差异大，大规模的干旱时有发生，给农业生产、人民生活和经济发展造成极大不便。河北省地处我国华北地区，属于温带季风气候，近年来干旱频发，作为我国重要的粮食产地，这严重制约了当地的农业生产。因此及时、精准、有效的干旱监测，有利于准确了解大范围的干旱状况，以便及时采取应对措施减轻灾害对粮食安全产生的危害，这对农业生产，经济建设具有重要的意义。风云四号卫星是我国发射的第二代静止轨道气象卫星，其搭载了四种国际先进的有效载荷，能全面提高我国对天气、灾害和环境的监测能力[1]。多通道扫描成像辐射计AGRI (Advanced Geostationary Radiation Imager)是FY_4A卫星上搭载的传感器之一，用于对地球气象、环境的二维成像观测。其设有14个探测波段，探测波长在0.45～13.8 μm之间，可探测到大气、陆表和海洋的可见光至热红外的发射辐射，为大面积的遥感干旱监测提供了一种新的遥感数据源。Himawari-8卫星是日本的静止轨道气象卫星，于2014年发射，其搭载的16通道的葵花卫星高级成像仪AHI (Advanced Himawari Imager)有着与风云四号卫星相似的通道设置[2]，具有高精度、高光谱、高空间分辨率等特点，对于风云四号卫星的应用研究具有重要的参考价值。

风云四作为新型的卫星遥感数据，目前应用研究不多。本文将利用FY-4 AGRI数据对2018年秋季河北省的干旱状况进行监测，并对比同期的Himawari-8 AHI数据，验证其监测精度，对比分析两种传感器在干旱监测中的表现，分析2018年秋季河北省的干旱状况。本文为探索新型国产气象卫星数据在地表环境监测中的应用进行了有益的尝试，对于新一代国产卫星数据的推广使用有积极的意义。

1 研究区概况与数据

1.1 研究区概况

河北省地处我国华北地区，位于113°04′～119°50′E，36°03′～42°40′N之间，东临渤海，西部与北部分别与山西、内蒙、辽宁相连，南部与山东、河南接壤。全省面积187 693 km2，地势西北高，东南低。西北部为燕山、太行山山地，海拔多在2 000 m以下；北部为坝上高原，属于蒙古高原的一部分，平均海拔在1 200～1 500 m，东南部平原属于华北平原的一部分，平均海拔在50 m左右，面积占全省的41.2%，冬小麦主要种植于此。气候上属于温带大陆性气候，四季分明，夏季高温多雨，冬季寒冷少雨。年均温在-3℃～14℃之间，年平均降水为537.1 mm，降水量夏季最多，冬季最少，秋季高于春季[3]。全省耕地面积约6.2万km2，农作物类型以小麦、玉米为主，主要分布在东南部平原地区(图1)。

图1 研究区及实测样点分布图

2018年秋季，河北省发生了持续性的气象干旱，平均干旱日数为64.8 d，较常年偏多27 d，为2007年以来最多[4]。本文统计了自2000年以来的河北省历年秋季平均降水量与多年平均降水量(图2)，发现2018年秋季，河北省全省平均降水量为33.3 mm，为2000年以来的历史最小值，降水量较常年平均偏少66.4%，属于显著偏少年份。因此就历史同期数据来看，2018年河北省秋季处于干旱状态，因此选取了河北省作为本文的研究区。

图2 河北省历年秋季平均降水量变化(2000-2018年)

1.2 研究数据

FY-4 AGRI数据采用HDF5作为存储格式，数据集包含了所有的光谱通道和定标通道。本研究采用1km的L1级数据，包含第1，2，3通道的光谱通道和定标通道，数据来源于国家卫星气象中心。Himawari-8 AHI数据采用分辨率为2km的全圆盘L1级数据，包括16个波段的辐射通道和经纬度文件，数据来源于日本宇宙航空研发机构(Japan Aerospace Exploration Agency)。

数据时间为2018-09-08、2018-09-24、2018-10-04、2018-10-24、2018-11-02、2018-11-25这6期数据，数据时间的选取依据是：①为反映干旱的时间变化状态，选取2018年秋季9-11月每月月初及月末的影像；②晴朗少云的天气状况，静止轨道卫星数据的高时间分辨率可在1 d的同一地区内获取多幅影像数据，数据选取的时间为该日内云量较少的影像时间。对所有遥感数据进行定位、定标和大气校正等预处理，生成研究区的反射率数据集。

本文的验证数据采用研究区土壤相对湿度的地基观测数据。本文收集了河北省18个国家级农气站的土壤相对湿度站点数据，取样深度为10cm和20cm。站点分布见图1。

2 研究方法

2.1 垂直干旱指数PDI

植被对红光强烈吸收，对近红外强烈反射，水体对红光和近红外波段均强烈吸收，而裸地反射率从红光到近红外变化很小，因此土壤含水量是影响土壤反射率的重要因素。土壤含水量越高，其在红光和近红外的反射率越低，含水量越低，反之反射率越高。基于土壤的这种光谱特性，阿布都瓦斯提·吾拉木、詹志明等[5-6]利用Landsat ETM+影像构建了近红外和红光的反射率特征空间，表征了不同水分状况下的土壤光谱变化规律，并由此建立了土壤水分监测模型——垂直干旱指数(Perpendicular Drought Index, PDI)。

在NIR-Red光谱特征空间中，二维坐标的每一个点对应影像中的一个像元，植被和水体像元分布于土壤线的上下两侧[7-8]，植被像元和裸土像元呈典型的三角形分布，裸土像元多分布于三角形的下方，这些裸土像元可提取一条土壤线，然后作土壤线的垂线，点到土壤线的垂直距离表征对应像元土壤水分的含量，二者呈负相关关系。其数学表达式为：

(1)

式中：RRED和RNIR分别为红光和近红外波段的反射率，M为NIR-Red光谱特征空间中由土壤点拟合得到的土壤线斜率。PDI的值在0～1之间，值越大，表示干旱越严重。

2.2 土壤线提取

土壤线是指在二维光谱特征空间中土壤像元按照光谱反射率的大小排列而成的直线，它反映了不同水分条件下的土壤含水状况[9]。本文参考秦其明的土壤线自动提取算法来拟合河北省的土壤线，该算法可保证模型精度的同时提高土壤线计算效率，是土壤线提取的常用解决办法。提取步骤如下。

(1)绘制二维光谱特征空间。以红光反射率为横坐标，近红外反射率为纵坐标，将遥感影像的所有像元绘制到二维散点图中。考虑到计算机的性能，直接绘制所有的像元点效率太低，因此可对影像进行适当采样，确定土壤点集。根据经验，当子集数为200～300时，土壤线的波动较小[10]。根据AGRI影像和AHI影像红光反射率值的范围，本文设定的步长为：

StepAGRI= 0.0004；StepAHI= 0.0005。

(2)

(2)确定初始土壤点集。裸土像元一般分布于二维光谱特征空间中散点图的下方。按照步骤(1)中设定的步长，将横坐标红光反射率分成若干子集，找到每个土壤子集中纵坐标的最小点，即每个横坐标分组中近红外反射率的最小值，构成(X, Y)点集，作为土壤线的初始点集。

(3)自适应区的选取。步骤(2)中提取的土壤点集并非是绝对的裸土像元，其中可能混入了大量的植被像元。自适应区的选取就是将初始点集的范围缩小，减少非裸土像元的数量，来完成土壤点集的初步筛选。本文初始筛选的方法是分别计算横坐标区间范围内0%～50%、0%～75%、0%～100%、25%～75%、25%～100%和50%～100%的土壤点集的相关系数，选取相关性最强的子集来估算土壤线。

(4)迭代筛选土壤点。将步骤(3)中获取的土壤点集进行循环迭代，每次都把垂直偏差最大的点从土壤点集中删除，以获取最终的土壤点集。最后拟合出土壤线，得到土壤线斜率。

2.3 云水检测

自动提取算法提取的土壤线像元分布于二维光谱空间的下方。这个算法虽然可以避开散点图中左上方的植被像元，但无法解决水体像元和云像元的影响，因此在实际应用中，需对影像进行云水检测，剔除云水覆盖的像元，以获取真实的土壤线。

2.3.1 云检测

本文参考MODIS的可见光和近红外云检测阈值法[11]，并将其运用到风云数据和葵花数据中。

MODIS可见光波段反射率云检测：晴空时，单波段的地表反射率无云的会比有云时低，红光的反射率会比近红外的低；而在有云覆盖的情况下，红光的反射率和近红外的相近。此方法是利用云和地表反射率的差异性来检测云，适用于地表和云反射率差异较大的情况下。MODIS 0.65 μm的红光波段反射率检测适用于陆地，0.86 μm的近红外波段反射率适用于海洋和沙漠。

MODIS可见光和近红外比值反射率云检测：红光和近红外在有云覆盖时反射率的差异很小，而在水面或植被区域差别较大，因此红光和近红外的比值反射率可以用来区分云和水体。

综合以上两种方法，将MODIS云检测算法应用于FY-4 AGRI 和Himawari-8 AHI中，找出AGRI和AHI中与MODIS相近的波段，适用于MODIS可见光与比值阈值法的波段分别为AGRI的第2、第3波段和AHI的第3、第4波段，各传感器对应波段信息见表1。

表1 MODIS与AGRI、AHI对应的云检测波段信息

本文以2018年9月8日的AGRI数据为实验数据，构建NIR-Red光谱特征空间，采用监督分类的方法，提取该影像中云覆盖的像元，计算红光波段反射率和红光-近红外波段的比值反射率在云覆盖中的阈值。最后对影像进行有云和无云覆盖的像元分类，去除有云像元。

2.3.2 水体检测

在NIR-Red光谱特征空间中，水体像元分布在特征空间下部，易与裸土像元混在一起，对土壤线自动提取算法的精度有很大影响，因此选取合适的方法提取水体像元尤为重要。常用的水体提取方法有三种：单波段阈值法、谱间关系法和水体指数法[12-14]。参考程磊[13]在黄土高原地区利用归一化水体指数(NDWI)并结合阈值法识别出了黄土高原的水体和其他下垫面，本文对水体的检测也采取归一化地表参数结合阈值法来识别水体。这种方法能在复杂的下垫面中准确识别出水体，且能有效去除山体阴影，提取出的水体更加真实准确。由于AGRI数据缺少归一化水体指数所需的绿光波段，因此本文选取归一化的植被指数(NDVI)来代替，并结合红光波段，给定阈值，最后检测出水体，去除水体像元。NDVI的表达式如下：

NDVI=(RNIR-Rred)/(RNIR+Rred)。

(3)

式中：Rred和RNIR分别为红光和近红外波段的反射率。NDVI值的范围为-1～1，理论上小于0的值表示覆盖地表是云、雪、水等；0值表示为岩土或褐土等土壤下垫面，大于0的值表示植被覆盖。但在实际应用过程中，水体覆盖的NDVI值也可能会大于0，因此需要对水体像元覆盖的NDVI赋给一个新的阈值。

同样以9月8日河北省的AGRI数据为例，采用监督分类的方法，提取该影像中水体覆盖的像元，计算水体像元在红光波段反射率和NDVI上的阈值。最后对影像进行水体和非水体的像元分类，剔除水体覆盖的像元。

3 结果分析

3.1 云水检测结果

3.1.1 云检测结果

以2018年9月8日的河北省AGRI数据为实验数据，该日河北省西北部有明显的云覆盖区域。图3为2018年9月8日的FY-4 AGRI云检测结果图，图3a为云检测前，图3b为云检测后，图3c为云像元在二维光谱空间中的散点图，红色部分代表云。对比云检测前后的结果，可以看出基于红光反射率和红光-近红外的比值反射率的阈值法可以较好地提取出云像元。全部影像云检测的阈值结果如表2所示。

表2 AGRI和AHI云检测阈值结果

图3 河北省2018年9月8日FY-4 AGRI云检测结果

3.1.2 水体检测结果

同样以9月8日河北省的AGRI数据为例，计算该日水体像元的阈值范围，分别为Rred<0.1和NDVI<0.26，结果见图4。图4a中可以看到，该阈值可以将大部分的水体像元去除，但仍然有少部分的靠近裸土像元的像元未被去除；图4b显示了被去除掉的水体像元，如黑色框中是河北省与渤海交界处的深水区，红色框中是河北省内的水库、河流等水域。此方法可以去除掉大部分水域尤其对深水区的去除效果明显，但在水面较浅的地区，水体像元与裸土区的像元在NDVI的值上有部分重合而无法全部去除。

图4 河北省2018-09-08 AGRI影像水体检测结果

为了保留NIR-Red光谱特征空间中相对完整的裸土像元，给定NDVI的阈值不应过高。利用NDVI结合红光波段反射率的阈值法，将全部传感器影像进行水体检测，阈值结果见表3。

表3 AGRI和AHI传感器水体检测的阈值结果

比较AGRI和AHI传感器水体检测的计算结果，发现单波段反射率的阈值AHI整体上比AGRI要大，但在NDVI比值反射率上二者间阈值结果相差不大，这是由于两种传感器对光谱响应的差异造成的，在单波段反射率上AHI传感器比AGRI传感器整体偏大。通过给定合理的阈值，AGRI数据和AHI数据的水体检测结果基本一致。

3.2 土壤线提取结果

基于土壤线自动提取算法对全部影像进行土壤线提取(表4)。可以看到随时间的不同，土壤线略有波动，但大体上是稳定的，土壤线斜率稳定在0.8～1.2左右。两个传感器提取的土壤线的最大值均出现在2018年10月24日，土壤线斜率的波动大体相同。可以认为基于该方法提取的土壤线是可靠的。

表4 AGRI和AHI的土壤线提取结果

3.3 PDI旱情分析

根据3.2中拟合的土壤线斜率，利用公式(1)对河北省的AGRI和AHI遥感影像计算垂直干旱指数PDI，将其结果导入ArcGIS中进行渲染(图5)。PDI的值在0～1之间波动，值越大，代表该地旱情越严重。结果显示，河北省2018年秋季整体呈干旱的趋势。9月，PDI指数相对较低，全省土壤墒情状况较好，9月8日AGRI计算的PDI均值为0.189，9月24日为0.2；10月干旱进一步发展，10月4日AGRI计算的PDI均值为0.22，10月24日为0.27；11月干旱程度最为严重，PDI值整体较高，11月2日AGRI计算的PDI均值为0.31，11月25日为0.24。河北省的干旱主要发生在东南大部分地区和西北、东北部少部分地区，东北西南山地植被区的PDI指数变化不明显，中南部华北平原对PDI的响应更明显。

对比AGRI和AHI结果图发现，整体上基于AGRI计算的PDI指数与基于AHI计算的PDI指数在干旱变化趋势上相对一致，但AGRI的PDI变化层次较AHI的更为清晰，反映出AGRI对土壤水分变化的敏感性。

对PDI指数进行干旱分级。参考前人基于PDI模型的干旱指数分级研究[15-16]，分析PDI与实测土壤相对湿度的关系，并结合河北省的实际情况，建立了PDI干旱等级的分级标准(表5)。

图5 河北省9-11月PDI结果图(上：AGRI；下：AHI)

图6 河北省2018年秋季干旱监测结果

表5 PDI、土壤相对湿度与干旱等级间的对应关系

利用这一标准对PDI指数进行干旱分级。由上述PDI结果可知，基于AGRI传感器和AHI传感器计算得到的PDI分布趋势基本一致，所以这里仅对AGRI影像进行干旱分级(图6)，9月旱情初步发展，河北省大部分地区土壤墒情良好，9月初仅华北平原和西北山地出现少量轻度干旱，9月下旬干旱面积有所发展，但全省仍然是以轻度干旱为主。10月旱情加重，受旱面积和干旱程度较9月明显增加，河北省中南部的华北平原地区，受旱面积最大，以轻度和重度干旱为主，河北省西部和东北部还出现少量重度干旱。10月下旬到11月上旬，旱情进一步发展，以轻度干旱和中度干旱为主，旱情呈现从南到北发展的趋势。到11月下旬，旱情有所缓解，河北省基本以轻度干旱为主。可以看到河北省的主要干旱地区集中在中南部的华北平原，华北平原是我国冬小麦的主要种植区域，华北平原旱情的发展给当地的农业生产造成了影响。

4 PDI精度验证及对比

本文收集了河北省18个国家级农气站站点的土壤相对湿度数据，分别与PDI进行相关性分析，验证结果如表6所示。可以看出，基于AGRI反演的PDI结果与实测数据呈负相关关系，二者大体上有较好的相关性，相关系数稳定在0.6左右，整体上11月25日的相关性最好，10月4日的相关性最差，最大值为0.82，最小值为0.46；基于AHI反演的PDI结果与实测数据呈负相关关系，相关系数在0.5～0.7左右，整体上9月24日的相关性最好，11月2日的相关性最差，最大值为0.74，最小值为0.42。

表6 AGRI和AHI与实测数据的验证精度

注：*通过0.05置信度检验；**通过0.01置信度检验

对比AGRI和AHI的精度结果，发现传感器反演的PDI均与实测数据呈负相关关系，且相关性普遍在0.5以上。AHI传感器在20 cm土壤深度的相关性要更高，但在10 cm土壤深度上，AGRI与实测土壤相对湿度数据的相关性更稳定。总体上，20 cm的土壤湿度数据与PDI之间的相关性均强于10 cm的土壤湿度。这是由于10 cm深度的土层易受近地表风速和气象条件的影响，并不能准确反映地表的反射率特性，而20 cm深度的土壤更接近作物根部，土壤水分对作物的影响更大，更能真实反映出作物的受旱情况，这与现有的研究发现一致。朱琳等[16]利用风云3号卫星MERSI数据对内蒙古东辽市和赤峰市进行PDI干旱监测，20 cm深度的土壤水分相对数据与PDI的相关性明显高于10 cm深度的土壤水分数据，R2在0.41～0.96之间；李峰等[17]同样利用MERSI数据对2010年山东省秋冬季节干旱进行PDI监测，20 cm深度的土壤水分与PDI的相关性更好，相关性在0.437～0.609之间。

以2018年9月8日这天的PDI结果为例，将基于AGRI和AHI反演的站点所在PDI值进行对比(图7)，发现AHI反演的PDI值整体偏高，这是由传感器对光谱响应的差异造成的，总体趋势上，AGRI反演的站点PDI值与AHI反演的值趋势相当。