基于GF-1卫星数据监测灌区灌溉面积方法研究——以东雷二期抽黄灌区为例

宋文龙，李 萌，3，路京选，卢奕竹，史杨军，贺海川

（1.中国水利水电科学研究院，北京 100038；2.水利部防洪抗旱减灾工程技术研究中心，北京 100038；3.首都师范大学 资源环境与旅游学院，北京 100048；4.渭南市东雷二期抽黄工程管理局，陕西 渭南 714000）

1 研究背景

准确监测灌溉面积对于预防干旱、优化灌区水资源利用和水资源调控配置具有重要意义。据统计，世界上约有农田1.5 亿～1.7 亿hm2，约占地球陆地总面积的17%［1-2］，农业用水量约占人类用水的92%［3］。作为农业大国，我国灌溉面积约占总耕地面积的50%，在全球气候变暖、干旱等极端事件增加、经济发展用水需求猛增情况下，水资源短缺形式日趋严峻。准确监测灌溉面积时空分布及其变化，是掌握灌溉用水量、灌溉用水去向的关键，对监测灌溉系统运转、预防干旱、优化作物种植结构、提高灌溉用水效率、优化灌区水资源配置均具有重要意义。

国内外农业灌溉面积监测研究在决策树算法［4-5］、时空螺旋曲线与变化矢量分析［6］、物候特征［7］、支持向量机［8］、自动农田算法［9］、光谱匹配［10］、k-means 和ISOCLASS 算法［11］、自动机器学习算法［12-13］、植被指数与干旱指数反演［14］等方法方面取得诸多进展。其中，光谱匹配技术（SMTs）能够提取长时间序列多期遥感影像的光谱反射率或特征指数变化特征，在周期性农业灌溉面积监测方面具有独特优势。Thenkabai 最先应用这种方法开展土地类型识别与灌溉面积识别研究［10］。2006年，世界水资源管理研究所（IWMI）通过使用多种数据源组成的具有巨大数据量的数据集，根据降雨量与高程信息对全球范围区域进行分类，并使用非监督分类对区域进行聚类，而后使用光谱匹配技术对全球范围内的土地利用类型、灌溉面积分布进行了识别，发布了世界第一份1 km 尺度全球灌溉区域图（GIAM）［15］。2012年，IWMI又针对非洲和亚洲进行了重点研究，绘制了250 m 空间分辨率的灌溉与雨养区分布图［16］。2016年，Ambika 将这种方法应用到了印度，使用了250 m 空间分辨率的MODIS NDVI 数据（MOD13Q1），先通过匹配计算分类，再使用决策树通过NDVI 阈值划分灌溉区域，获得了印度2000—2015年的灌溉区域分布图［17］，相比IWMI 250 m 空间分辨率的灌溉产品提升了精度。Teluguntla 结合定量光谱匹配技术（QSMTs）与农田自动分类算法，使用连续14年的250 m 空间分辨率MODIS NDVI 产品数据对澳大利亚的灌溉信息与种植强度进行了提取［18］。

我国具有气候多样、农业种植结构复杂、田块破碎、灌区信息化水平不高等国情，农业统计缺乏空间信息，可公开获取的灌溉遥感产品缺乏且空间分辨率与精度难以满足我国实际需求［19］，我国灌溉面积遥感监测对数据空间分辨率要求更高，灌溉面积遥感监测技术需要取得进一步突破。本文将针对我国国情，应用较高分辨率的国产卫星遥感数据，通过光谱匹配方法像元尺度应用，并引入自适应阈值算法，构建高分辨率灌溉面积遥感监测新方法，对我国西北干旱半干旱区典型渠灌灌区开展灌溉面积遥感监测研究。

2 研究方法与数据来源

2.1 研究区域概况东雷二期抽黄灌区（见图1），位于东经109°10′—110°10′，北纬34°41′—35°00′，横跨陕西省富平、大荔、蒲城三县，属黄土波状台原区，海拔介于385～600 m，地势东南低，西北高。灌区属温带大陆性季风气候，春季温暖干旱，夏季多暴雨，并有短暂大风，冬季严寒干燥，雨雪稀少，年平均降水量519～552 mm，干旱年份仅为360～390 mm，且降水时空分配不均衡，50%以上雨量集中在7—9月，年蒸发量约为1700～2000 mm。灌区灌溉方式为抽取黄河水进行渠道引水漫灌。设计灌溉面积8.5 万hm2，分春灌、夏灌和冬灌。主要作物为小麦和玉米，占灌区种植总面积的90%以上。主要作物生长周期及灌溉时间详见表1。

图1 东雷二期抽黄灌区示意图

表1 主要作物生长周期及灌溉时间

2.2 数据

2.2.1 卫星遥感数据 数据源为2018年GF-1 号16 m 多光谱数据，来源为中国资源卫星公共网站平台（http：//www.cresda.com/CN/），详细参数见表2。使用2018年研究区云覆盖量低于15%的影像，并编译成由92 个波段组成的数据集（共有23 景数据，每景4 个波段）。此数据集用于准确获取不同类别的时序光谱特征信息［20-21］。

对遥感数据进行辐射定标、几何校正、大气校正等预处理，并计算NDVI 得到研究区各像元的NDVI 时间序列曲线。NDVI 能够较好地反映植被生长状态，进行归一化处理后可以更好地限制数据范围，有利于后续的光谱匹配计算，用于提取灌溉信息。

表2 GF-1 多光谱传感器信息

2.2.2 实测样点数据 像元尺度的光谱匹配算法需要端元光谱特征信息，即灌溉区域小麦和玉米的NDVI 时间序列曲线，因此在研究区开展野外调研与采样工作，获得研究区主要地物的真实地面样点，用于提取灌溉面积及精度验证。

结合历史灌溉情况与卫星遥感影像，选择具有典型代表意义的灌溉区域与非灌溉区域，在选取灌溉区域时应注意避免附近有建筑物与道路的干扰，避开农田的边界处，以防在取点时样本点中同时包含农田与其他干扰地物。最终在研究区内共选取225 个真实样点数据（见图2），主要记录样点位置、作物类型及灌溉情况等信息。

图2 样点分布

2.2.3 IWMI 灌溉产品 使用国际水管理研究所（International Water Management Institute，IWMI）的数据产品IWMI-Irrigated Area Map Asia（2000—2010）and Africa（2010）作为验证参考数据（http：//waterdata.iwmi.org/）。该产品使用250m 空间分辨率MODIS NDVI 合成数据（MOD13Q1），通过ISODATA 算法进行非监督分类，之后依据各种类的NDVI 时间序列曲线确定土地利用类型，得出亚洲和非洲的灌溉区域分布图。IWMI 灌溉产品将农业区域分为灌溉区和雨养区，并且根据作物种植强度将灌溉区具体划分成单季作物灌溉区域、双季作物灌溉区域和连续作物灌溉区域。

2.3 研究方法有别于IWMI 采用的聚类光谱匹配，本文针对高分辨率遥感影像，将光谱匹配方法应用于像元尺度，保证所有像元的时序光谱曲线参与匹配计算，减少聚类过程造成的误差，并引入OTSU 自适应阈值算法自动确定灌溉面积提取阈值。OTSU 算法的引入可保证每个像元计算相似度信息的同时，自动确定相似度合理变化范围，对于不同年份与不同数据情况有良好变动适应性，无需人为干预即可提取灌溉区域。基于像元尺度光谱匹配计算的灌溉面积遥感监测方法适用于高分辨率卫星遥感数据，满足提取小地块灌溉信息的需求，能够提高灌溉面积监测结果的精度。

2.3.1 光谱匹配方法 光谱匹配法是一种量化端元光谱（样本光谱）与目标光谱（待测光谱）相似度的方法。根据光谱匹配原理，光谱匹配计算方法可以分为统计算法、光谱波形特征算法、光谱编码匹配算法以及特征空间算法。本文使用统计算法和光谱波形特征算法计算光谱相似度，通过三个指标对研究区主要作物端元光谱与目标光谱匹配程度进行定量分析［10］，具体算法如下：

根据地面样点数据提取端元光谱C=（s1，s2，…，sn），从数据集中获取目标光谱Xi，j=（b1，b2，…，bn）。

（1）形状定量。光谱关联相似度（Spectral Correlation Similarity，SCS）计算公式如下：

式中：Xi，j为目标光谱的NDVI 值；C 为端元光谱的NDVI 值；i，j 分别为第i 行、第j 列；n 为数据集层数；μi，j为目标光谱NDVI 均值；μs为端元光谱NDVI 均值；σX为目标光谱的NDVI 标准偏差，σC为端元光谱的NDVI 标准偏差。

SCS 是衡量端元光谱与目标光谱NDVI 时间序列曲线形状相似度的指标。SCS 值越高，端元光谱与目标光谱的NDVI 时间序列曲线就越相似。当SCS 值为0 时，相似度最小；当SCS 为1 时，相似度最大，即两种光谱完全一致。所有SCS 值在0～1 范围之外的像元即可认为与端元光谱相似度差异过大，直接剔除。

（2）距离定量。欧氏距离（Euclidian Distance Similarity，EDS）是用于衡量端元光谱与目标光谱在光谱空间中距离的指标，在光谱空间中距离越近则越相似。EDS 计算公式如下：

通常欧氏距离越大则代表着两种光谱差异性越大，反之则差异性越小。为了方便计算将上述公式获得结果通过归一化处理，将其标准范围处理至0 到1。具体公式计算如下所示：

式中，EDSnormal值的范围为在0～1，它从光谱特征空间距离上测量端元光谱与目标光谱NDVI 时间序列曲线间相关性的大小。M，m 分别代表最大、最小的欧氏距离。通过对EDS 的计算可以把SCS 运算结果中满足条件像元点进行计算，理想状态下，0 表示二者完全一致，1 表示二者完全不相关。此外，EDS 受像元数量的影响，一旦有新的像元参与计算，该式中的M 和m 即可能会发生变化，对不同研究对象有较好的变动适应性。同时欧氏距离的局限性在于不同年份气候情况可能导致NDVI 光谱曲线会发生上下平移，进而会影响判断结果的一致性。因此，仅使用光谱距离的特征描述往往还不足以较好的完成光谱匹配的工作，需要进一步对端元光谱与目标光谱的相似程度进行量化。

（3）形状距离综合定量。光谱相似值（Spectral Similarity Value，SSV）综合了SCS（形状定量）与EDS（距离定量）的特点，从NDVI 时间序列曲线形状和光谱特征空间距离两方面衡量了端元光谱与目标光谱间的形似度，其计算公式如下：

式中SSV 值越小，光谱之间越相似，其范围通常介于0～1.414 之间。不同作物之间，NDVI 时间序列曲线的差异较大，SSV 值高；对于同一种作物，灌溉区域的NDVI 时间序列曲线比非灌溉区域具有更高的一致性，SSV 值更低。

2.3.2 OTSU 算法 以实地采样获取的灌溉区域玉米和小麦的NDVI 时间序列曲线分别作为端元光谱，利用上述衡量光谱相似度的指标计算研究区每个像元的SSV 值。同类型的光谱相似度越高，SSV值越低；不同类的光谱相似度越低，SSV 值越高。需要确定合理分割阈值，当SSV 值小于该阈值时与端元光谱识别为同一类别。因此，引入OTSU 自适应阈值算法计算SSV 分割阈值来判断是否为小麦或玉米的灌溉区域，小于该阈值即为灌溉区域，从而识别研究区的灌溉面积空间分布情况，极大减少了人为干预造成的影响。

OTSU 算法是无参量的一种自适应阈值选取方法［22］，该方法可以给定一个临时阈值，然后计算阈值两侧范围的像元SSV 的方差值，当方差达到最大所对应的像元SSV 为最佳阈值，即目标区域与背景区域差异最大化，而后将相似度大于或等于该阈值的归并成一类，小于该阈值的归并成另一类，得到相应的二值化图像。具体原理算法如下：

当影像数据量级为L（G=0，1，…，L-1）时，初始阈值为h 将图像分为目标区域T 与背景区域B 两部分，pi代表SSV 为i 的像元出现的概率，ni代表SSV 为i 的像元数量，N 代表像元总数，T、B 两部分概率为：

对应的T、B 两部分的均值为：

则存在阈值h0即为最佳阈值。

3 结果分析

3.1 主要作物端元光谱通过野外实地采样，获得了玉米与小麦灌溉区域的端元光谱（见图3）。图3（a）为玉米灌溉区（6月至10月）9 个实测样本点的NDVI 时间序列曲线，从中筛选5～6 条相似程度最高的曲线作为端元参考光谱，然后再使用端元光谱均值化进行光谱匹配计算。均值化可以保证所获得的均值光谱在代表该种类光谱特征信息的同时减少数据波动对结果的影响，且进行种类光谱廓线匹配计算时缩小同类别差异性。同理，图3（b）为小麦灌溉区（10月至6月）共10 个实测样本点的NDVI 时间序列曲线，提取小麦端元光谱均值并计算。

图3 玉米与小麦端元光谱

3.2 主要作物灌溉区域相似性参数空间分布根据像元尺度光谱匹配算法得到的主要作物灌溉区域识别相似性参数SCS、EDS 和SSV 分布如图4所示，研究区内未能灌溉到的地方主要位于灌区最北（区域1）、中部（区域2）与右下位置（区域3）。根据实地调查了解，最北侧主要是由于没有渠系覆盖；而蒲城附近（区域2）是由于靠近城镇，在进行城镇建设时灌溉渠系损坏导致灌溉情况较差；右下位置为大荔系统，以种植果树为主，主要作物种植面积相对较少。因此提取结果与实际情况相符。

图4 主要作物灌溉区域相似性参数空间分布

3.3 灌溉区域相似性分割阈值通过光谱匹配方法获得了研究区主要作物小麦与玉米相似度量化空间信息。进而对获得的量级分布直方图使用OTSU 算法得到两种作物对应的相似度分割阈值（见图5），结果表明小麦灌溉区最佳分割阈值为0.3985，玉米灌溉区最佳分割阈值为0.4639。

图5 OTSU 法计算直方图

3.4 主要作物灌溉区域分布应用小麦和玉米灌溉区域相似性分割阈值，得到主要作物小麦和玉米的灌溉区域空间分布图（见图6）。根据作物种植强度情况，研究区内的灌溉区域可分为双季轮作灌溉区（小麦和玉米轮作）、单季小麦灌溉区和单季玉米灌溉区，以双季轮作灌溉区为主，单季作物灌溉区主要分布在城镇周围与渠系未灌溉区域，由于地块破碎、种植结构复杂造成。多数非植被覆盖区域（如城镇、村庄、渠道、道路等）识别效果较好，但在一些小的建筑物周围有小部分误识别区域，是由于混合像元所导致。

对灌区各灌溉子系统（各子系统分布参见图2）的主要粮食作物小麦和玉米的灌溉面积进行统计（见表3）。各灌溉子系统灌溉面积由大到小依次是流曲、孙镇、兴镇、荆姚、刘集、蒲城和大荔。除大荔外，其他灌溉子系统均是小麦和玉米轮作灌溉区为主，占比大于50%，荆姚占比达到81.72%。大荔果树种植面积较大，粮食作物灌溉主要是单季小麦和单季玉米，小麦和玉米轮作灌溉较少。

图6 主要作物灌溉区域分布图

表3 各灌溉子系统灌溉面积情况 （单位：hm2）

3.5 精度验证使用实测样点数据对提取结果进行混淆矩阵精度验证（见表4），使用高精度数据提取灌溉面积与种植强度结果的总体精度为88.27%，Kappa 系数为0.8308。对于分布最广的双季轮作灌溉区域的提取效果最佳，精度达到90.14%；非灌溉区域的提取精度为89.80%，其主要与单季小麦、单季玉米灌溉区发生混淆，与单小麦灌溉区误差较大，种类错分误差为13.04%；单季玉米灌溉区提取精度为89.47%；单季小麦灌溉区的提取精度相对最差，为78.26%。结果分析发现主要是混合像元对提取结果影响最大，单季小麦与单季玉米灌溉区主要分布在城镇附近，受人类活动影响，导致种植结构分布复杂且细小散碎，提取难度较大。

表4 GF 基于实测样点的结果精度评价 （单位：%）

3.6 合理性分析图7为国际水管理研究所（IWMI）发布的MODIS-250m 全球灌溉区域分布成果的研究区部分，是目前可公开获得精度较高的灌溉区域分布数据，对其重采样成16 m 分辨率作为参考数据。东雷二期研究区总面积为99 008.06 hm2，IWMI 研究区灌溉面积结果统计为82 745.37 hm2，使用GF-1 结果统计灌溉面积为81 571.58 hm2，二者计算灌溉面积结果差异为1.42%。对本研究获得的灌溉区域与IWMI 灌溉区域数据进行空间比较，空间重叠率为87.29%。但二者在不同种植强度灌溉区类型及其分布方面提取结果存在较大差异，据灌区实际调研与本研究结果显示，研究区主要包括双季轮作（小麦和玉米轮作）灌溉区域与单季作物（小麦、玉米）灌溉区域，而IWMI 发布的成果仅分为双季作物灌溉区域与三季/连续作物灌溉区域，本研究成果与实际情况相符，而IWMI 成果有误。

如图8所示，本研究成果与IWMI-Irrigated Area Map Asia（2000—2010）相比，在空间分辨率与地物识别方面具有显著优势，可以更有效区分出建筑物、道路、渠道与田块信息；受数据分辨率限制，IWMI 结果无法获取灌溉田块分布细节，较多的建筑物与道路存在混淆错分情况。 与IWMI 结果比较，本研究成果能更好地反映我国灌区实际灌溉情况，一方面是由于采用了更高空间分辨率的遥感数据，另一方面相较于基于聚类的光谱匹配计算方法，基于像元尺度的光谱匹配方法适用于较高空间分辨率的灌溉区域遥感监测研究。

图7 IWMI 发布的灌溉区域分布

图8 基于GF-1 的灌溉面积监测结果与IWMI 结果细节对比

4 结论

本研究结合野外调研采样，应用16 m 较高空间分辨率GF-1 卫星影像，获取了研究区主要作物的时间序列光谱特征信息，构建了像元尺度光谱匹配方法，并使用OTSU 算法自动获得灌溉区域分割阈值，获得主要作物种植强度及其灌溉区域分布，取得的主要结论如下：

（1）2018年东雷二期抽黄灌区灌溉面积为81 571.58 hm2，其中双季轮作（小麦与玉米轮作）灌溉面积为40 335.88 hm2，单季小麦灌溉面积为15 276.94 hm2，单季玉米灌溉面积为14 059.14 hm2。各灌溉子系统灌溉面积由大到小排序依次是流曲、孙镇、兴镇、荆姚、刘集、蒲城和大荔。

（2）使用地面调研采样数据对灌溉面积提取结果进行验证，总体精度为88.27%，Kappa 系数为0.8308。对双季轮作（小麦与玉米轮作）灌溉区域与单季玉米灌溉区域提取效果较好，精度分别为90.14%和89.47%；单季小麦灌溉区域提取精度较差，精度为78.26%，易与其他类型地物发生混淆。

（3）本文提出的像元尺度的光谱匹配计算方法，更适用于基于高空间分辨率遥感数据的灌溉信息提取，通过引入自适应阈值算法，最大程度减少人为影响，保证了灌溉信息提取的较高精度。对比IWMI 结果，本研究成果可以更有效区分城镇、村庄、道路、渠道、田块等边界信息，可有效识别小田块灌溉分布，在作物种植强度及其灌溉区域分布方面更符合实际情况。