基于NDVI时序数据的不同水源灌溉面积分类研究—以人民胜利渠灌区为例

韩宇平 ，冯 吉，陈 莹，黄晓东，代小平 ，潘韶春，殷欢庆

（1.华北水利水电大学，郑州 450045； 2.水资源高效利用与保障工程 河南省协同创新中心，郑州450046；3.水利部 水资源管理中心，北京 100038； 4.河南理工大学，河南 焦作 454000；5.河南省节水农业重点实验室，郑州450045； 6.河南省人民胜利渠管理局， 河南 新乡453000）

0 引 言

我国70%的粮食由灌溉耕地生产，灌溉面积与粮食产量、农民收入、地下水开采等联系紧密。灌溉面积是反映区域水土资源利用的重要指标，耕地灌溉面积本身是灌溉面积的主要构成部分[1]。研究灌溉面积的变化是探索区域粮食产量、农民收入和地下水开采变化的基础。实时监控灌区灌溉面积的变化是提高灌溉管理水平、保障粮食安全、促进灌区可持续发展的基础。

基于遥感数据的灌溉区域分类研究主要集中在灌溉面积和非灌溉面积的区分上。NILTON等[2]基于多时相的Landsat TM数据，采用监督分类的方法提取了巴西东南部的水稻灌溉面积。ZHU等[3]基于NDVI的时间变化和降雨数据构建指标，提取了2000年中国的灌溉面积。沈静[4]、易珍言[5]等基于HJ1A/1B CCD数据，分别采用垂直干旱指数PDI和修正后的垂直干旱指数MPDI提取了内蒙古河套灌区的实际灌溉面积。王啸天[6]构建了基于垂直干旱指数（PDI）的灌区实际灌溉面积监测模型，计算出秦汉灌区4―8月各阶段的灌溉面积与分布。高瑞睿[7]基于气象数据和MODIS 1 km产品数据，计算了河套灌区义长灌区的蒸散发量，用修正的垂直干旱指数模型(MPDI)获取了研究时段始末的土壤含水率，又引入了农作物像元丰度，将1 km尺度灌溉面积转换到田间尺度，并提出了1 km尺度的灌溉面积计算公式，获取了义长灌区的灌溉面积。焦旭[8]、何娇娇等[9]利用Landsat8影像数据，基于遥感地表温度LST及植被供水指数VSWI对石津灌区的灌溉面积进行了提取。宋文龙等[10]基于GF-1较高空间分辨率卫星数据，通过光谱匹配像元尺度应用，并引入OTSU自适应阈值算法，对东雷二期抽黄灌区2018年的主要粮食作物种植强度及其灌溉面积开展了遥感识别提取研究。国内外关于不同水源灌溉面积的遥感分类研究较少。Velpuri等[11]基于NOAA AVHRR10 km、Terra MODIS 500 m、Terra MODIS 250 m和Landsat TM30 m 4种不同分辨率的遥感影像，用非监督分类的方法将印度克里希纳流域区分为地表水灌溉区和地下水灌溉区，发现遥感数据的空间精度越高，灌溉面积提取的精度则越高。Biggs等[12]，Kvisma等[13]基于 MODIS 500 m和 MODIS 250 m遥感影像，采用NDVI时间序列和非监督分类方法提取了印度克里希纳流域的地表水和地下水灌溉区域。国外研究表明基于NDVI时序数据的分类方法在流域尺度可以区分不同的灌溉水源，但该方法在灌区尺度上的适用性还缺乏检验。

MODIS数据的时间分辨率高，但空间分辨率低；而Landsat8数据具有空间分辨率高，时间分辨率低的特点。本文以人民胜利渠灌区为例，基于 Landsat8和MODIS来源的NDVI时序数据，采用监督分类和非监督分类方法对灌区不同水源的灌溉面积进行分类，探讨不同精度的NDVI时序数据和不同分类方法在区分不同水源的灌溉面积中的有效性，以推动灌溉面积遥感解译研究，为灌区管理、水土资源规划等提供技术支持。

1 研究区概况与数据来源

1.1 研究区概况

人民胜利渠灌区是黄河下游兴建的第一个大型引黄自流灌溉灌区。灌区于1952年建成，位于河南省新乡市，介于 34°58′—35°50′N，113°30′—114°27′E之间，灌区范围涉及封丘、滑县、辉县、获嘉、淇县、卫辉、新乡市郊、新乡县、延津、原阳（见图1）[14]。

图1 研究区域图Fig.1 Study area map

人民胜利渠灌区耕地面积15.08万hm2，灌区设计灌溉面积12.32万hm2，有效灌溉面积9.06万hm2[15]。灌区内主要种植冬小麦和夏玉米，还有水稻、大豆和花生等。灌区的主要灌溉水源为黄河水和地下水，还有极少部分耕地使用共产主义渠水和孟姜女河水进行灌溉。近年来，随着灌区经济与社会的发展，水资源短缺问题日益严峻，引黄水量越来越少[16]。由于引水条件恶化等原因，灌区渠灌面积逐渐减少，井灌区域增多，部分区域地下超采严重。探明灌区渠灌和井灌分布对优化灌区水资源配置，扼制地下水超采具有重要意义。本文中渠灌指纯渠水灌溉，井灌指纯井水灌溉。

1.2 数据来源及处理

1.2.1 Landsat8数据与处理

人民胜利渠灌区的 Landsat8数据来源于中国地理空间数据云平台。该数据为经过辐射校正和几何校正的Level 1T地形矫正影像，空间分辨率为30 m。为了避免不同作物对灌溉面积分类的影响，选取人民胜利渠灌区2016年11月—2017年5月小麦生长阶段的11期Landsat8数据进行分析。

下载的 Landsat8数据无须再做辐射校正和几何校正处理，采用ENVI 5.3对其进行辐射定标、大气校正、图像镶嵌和矢量裁剪等预处理。预处理后的影像数据的NDVI计算式为：

式中：NDVI为归一化植被指数；NIR为近红外波段；R为红光波段。

对NDVI数据进行波段组合，可得到研究区的NDVI时序数据集，Landsat8数据处理流程见图 2。因NDVI时序数据受到太阳高度角、传感器噪声、气溶胶和云污染等自然因素的干扰，致使数据呈现异常点波动变化，本研究使用ENVI的 Savitzky-Golay Filter工具对NDVI时序数据进行了去噪。

1.2.2 MODIS数据与处理

采用的MODIS数据为Terra传感器接收的2011年11月1日—2012年6月16日的中国500 MNDVI5d合成产品MODND1F影像，共46期数据，时间分辨率为5 d。该产品由MODND1D的5 dNDVI的最大值计算得到。

MODISNDVI数据的预处理只需在ENVI中对影像进行几何校正，然后对影像进行矢量裁剪，得到预处理后的卫星影像。之后，将处理后的 MODND1F数据产品波段组合，得到研究区MODISNDVI时序数据集，MODIS数据处理流程见图2。为了剔除其中的异常数据值，提高数据的有效使用性，同样在数据使用前对MODISNDVI时序数据集进行了去噪。

1.2.3 地面实测数据集

采用Oregon739北斗测绘采集器（手持GPS）对研究区不同水源灌溉的耕地进行实地调查，定位精度为3 m，设置测量仪器坐标系统为WGS-84坐标系，测量值采用地理坐标。

2017年8月7―15日，项目组调查了灌区内获嘉、新乡和延津3个县16个村的灌溉水源、种植结构、灌溉日期等情况，采集了38个样本点的坐标信息。其中渠水灌溉样本点11个，井水灌溉样本点11个，城镇和居民点样本16个。

为进一步验证灌溉面积遥感分类结果，2019年3月20―22日，项目组进行了二次调研。通过询问样本村了解该村近几年的灌溉水源、种植结构、灌溉日期等情况，并利用手持 GPS采集不同灌溉水源和种植结构的耕地以及果园、居民点的坐标信息。样本点根据空间分布均匀且重点关注井渠结合灌溉区域的原则确定（样本分布见图 3），共采集样本点 37个，其中，渠水灌溉样本6个，井水灌溉样本23个，井渠结合灌溉样本8个。采样时，在靠近渠道且灌溉水源仅为渠水的大面积耕地内选择渠水灌溉样本；在远离渠道而且灌溉水源仅为井水的耕地内选择井水灌溉样本；井渠结合灌溉样本根据实际调查确定的耕地范围定位。调查发现2017年人民胜利渠灌区的灌溉水源和种植作物与2019年相比差异较小，因此，可以用2019年的调查点来验证2017年的灌溉情况。

此外，项目组在2013年调查了灌区内获嘉县部分村镇的种植结构和灌溉情况，获得渠水灌溉样本点8个，井水灌溉样本点2个，井渠结合灌溉样本点4个。

图2 遥感数据处理流程图Fig.2 Flow chart of remote sensing data processing

图3 样本分布图Fig.3 Sample distribution

2 研究方法

本研究使用的NDVI时序数据是指不同时段的NDVI数据按时间先后顺序叠加而成的数据集。可以认为每个时间代表1个波段，在时间轴上每个波段的NDVI值就会形成1条NDVI时间序列曲线。理论上，不同水源灌溉的作物NDVI时间变化曲线具有不同的特征。由于渠灌的灌溉水量多于井灌，渠灌作物的NDVI峰值更大，并且能在高值区保持更长时间；井灌的灌溉水量较少，井灌作物的NDVI峰值较小，NDVI下降快[13]。渠灌的供水时间不灵活，用水不方便，通常井灌区域先灌溉，渠灌区域后灌溉，井灌区域的NDVI比渠灌区域增长地更早[12]。由于井灌区灌水更频繁，井灌作物的NDVI值下降的幅度小于渠灌，井灌作物的NDVI曲线的波动程度较小。基于不同水源灌溉区的NDVI时序曲线的特征，可以采用监督分类和非监督分类方法对灌溉水源进行分类。

2.1 监督分类

监督分类原理为：将多个波段的NDVI叠加图像作为1个多维空间，每个样区的光谱在此空间中转换成1个由其在各个波段上像元值决定的单位向量。将其他样区向量与典型样区向量进行比较，如果二者夹角在一定范围内，则认为是同一种类型[17]。这种分类方法不仅能有效避免邻近像元散射的影响，而且对于光谱曲线形状的匹配较为理想。监督分类在本研究中的思路为基于实际调查获取典型区域的灌溉水源、种植作物情况，根据实际调查点的NDVI时间序列曲线，提取不同灌溉水源灌溉的耕地的典型NDVI曲线，根据典型曲线对全灌区的灌溉水源进行分类。采用ENVI的“Spectral Library Builder”命令将不同类别的典型样本的NDVI时序数据曲线存储为ENVI光谱库，采用波谱角填图（Spectral Angle Mapper，SAM)工具对灌区其他区域的光谱曲线和典型样本的光谱曲线进行匹配和分类。经过多次试验，本研究将0.4°定为最大向量夹角。

2.2 非监督分类

非监督分类，也称“聚类分析”或“点群分类”，是在多光谱图像中搜寻、定义其自然相似光谱集群的过程。非监督分类不必对地物获取先验知识，仅靠图像上不同类地物光谱或纹理信息进行特征提取，再统计特征的差别以达到分类的目的，最后对已经分别的属性进行确认。非监督分类在本研究中的思路为根据曲线的差异利用非监督分类方法对NDVI时序数据进行自动识别分类。根据分类结果提取各类别的曲线，分析曲线变化特征合并类似曲线，得到最终的分类结果。本文采用K-means非监督分类方法比较灌区不同栅格的NDVI时间序列曲线特征进行分类。K-means利用各聚类中对象的均值获得一个“中心对象”（引力中心）来进行计算，然后迭代地重新配置中心对象，直至完成分类过程[18]。采用ENVI中的 K-Means Classification工具进行分类。通常分类数量越多，得到的分类结果越理想；迭代次数越大，得到的结果越精确，运算时间也越长。本研究选择的分类数量为40，最大迭代次数为15。

3 结果与分析

3.1 基于Landsat8数据分类结果

3.1.1 监督分类结果

2017和2019年所有调查样本的光谱曲线见图4。分析不同灌溉类型的NDVI时序数据曲线可知：不同水源灌溉作物的NDVI曲线具有相似的波动特征，但NDVI的数值不同。井渠结合灌溉类型的NDVI数值最高，曲线的波峰值最大；井水灌溉类型的DNVI数值最低，曲线的波峰值最小；渠水灌溉类型的NDVI数值以及曲线的波峰值介于井渠结合灌溉和井水灌溉之间。

为了选取准确的典型样本曲线，去除同类型中明显不符合变化趋势的曲线，选择变化趋势相同且大小相似的多条曲线中较平均的曲线，选取了8个样本点的NDVI时序曲线作为典型样本光谱曲线（图4右）。

图4 基于Landsat8数据的调查样本和典型样本光谱曲线图Fig.4 Spectrum curve of survey sample and typical sample based on Landsat 8 data

图5 基于Landsat8 数据的监督分类结果Supervision classification result based on Landsat 8 data

图6 基于Landsat8数据的非监督分类结果Fig.6 Unsupervised classification result based on Landsat 8 data

图7 基于Landsat8数据的非监督分类各类别曲线Fig.7 Unsupervised classification curves based on Landsat 8 data

图8 基于Landsat8数据的非监督分类结果合并图Fig.8 Unsupervised classification results graph based on Landsat 8 data

根据典型样本曲线进行分类，结果见图5。渠水灌溉主要分布在人民胜利渠的渠首，井渠结合灌溉分布在灌区的中上游位置，而井水灌溉分布在整个灌区。各类别面积为：渠水灌溉面积6 685.92 hm2；井水灌溉面积128 193.48 hm2；井渠结合灌溉面积3 601.80 hm2；城镇和居民点面积53 874.36 hm2。

3.1.2 非监督分类结果

基于Landsat8数据的非监督分类数量为 40，分类结果见图6，各类别曲线见图7。分析曲线特征并对比实际土地利用资料，将40种类别进行合并，得到最终的分类结果见图8。由图8可知基于Landsat8数据的非监督分类只分出耕地和非耕地2种类别，并没有区分出不同水源的灌溉面积。耕地面积113 702.85 hm2；非耕地面积79 112.43 hm2。

3.2 基于MODIS数据分类结果

3.2.1 监督分类

2013年调查样本的光谱曲线见图9。分析不同灌溉类型的NDVI时序数据曲线可知：渠水灌溉类型的NDVI时序数据曲线的波峰值最大，波谷值介于井灌和井渠结合灌溉之间，NDVI在高值区的持续时间较长；井渠结合灌溉类型的NDVI曲线的波峰值较大，波谷值最大，NDVI在高值区的持续时间较长；井水灌溉类型的NDVI时序数据曲线的NDVI峰值最小，NDVI在高值区的持续时间最短。该结果和理论假设基本一致，即灌溉水源的灌溉水量越大，NDVI的波峰值越大，NDVI在高值区的持续时间更长。井渠结合灌溉曲线与渠灌曲线的差别不大，可能原因为井渠结合灌溉的区域用渠水较多，井水较少。本研究选取了4个样本点的NDVI时序曲线作为典型样本的光谱曲线（见图9）。

图9 基于MODIS数据光谱曲线Fig.9 Spectral curves of survey samples and typical samples based on MODIS data

根据典型样本曲线进行分类，分类结果见图10。各类别面积为：渠水灌溉面积35 576.19 hm2；井水灌溉面积99 029.77 hm2；井渠结合灌溉面积26 365.75 hm2；城镇和居民点面积33 256.00 hm2。渠水灌溉主要分布在渠首位置，井渠结合灌溉主要分布在灌区中游，井水灌溉分布在中游和下游。分类结果与 2012年的实际灌溉情况相比，渠水灌溉和井水灌溉比较符合实际情况，井渠结合的分类效果不太理想。与Landsat8数据的监督分类结果相比，各类别面积的分布特征相近，但是分类精度较低。这是由于 MODIS数据的空间分辨率较低，使得分类结果成片显示。

图10 基于MODIS数据的监督分类结果图Fig.10 Supervised classification result graph based on MODIS

图11 基于MODIS数据的非监督分类结果图Fig.11 Unsupervised classification result graph based on MODIS

图12 基于MODIS数据的非监督分类各类别曲线图Fig.12 Unsupervised classification curves based on MODIS data

图13 基于MODIS数据的非监督分类合并图Fig.13 Unsupervised classification merging graph based on MODIS data

3.2.2 非监督分类

图例数据的非监督分类数量为40，结果见图11，各类别的NDVI时间序列曲线见图12。分析曲线特征并对比实际土地利用资料，将40种类别进行合并，得到最终的分类结果见图 13。各类别面积为：渠水灌溉面积17 260.78 hm2；井水灌溉面积79 097.26 hm2；井渠结合灌溉面积47 985.67 hm2；城镇和居民点面积49 989.47 hm2。与监督分类结果相比，非监督分类的结果较差，渠灌面积和井灌面积较少，井渠结合灌溉面积、城镇和居民点的面积较多，而且各类别面积的分布与实际情况不符，尤其是井渠结合灌溉面积。基于Landsat8数据的非监督分类只分类出耕地面积与非耕地面积，与之相比，基于图例数据的分类效果较好。

3.3 分类精度分析

3.3.1 基于Landsat8数据的分类结果验证与分析

将2019年和2017年实测GPS验证样本与监督分类结果进行比对，验证结果见表1。井水灌溉面积的分类精确度最高，但渠水灌溉和井渠结合灌溉面积的分类精度均小于60%。井渠结合灌溉类型的分类精度最低，验证井渠结合样本点后，发现错误分类为井水灌溉，可能这些区域用井水灌溉的比例高于渠水。除去城镇和居民点样本，不同水源灌溉面积的总分类精度为73.58%。

由图8可知，Landsat8数据的非监督分类结果并不理想，只分类出土地利用类型，并没有分类不同水源的灌溉类型，表明基于Landsat8数据的非监督分类方法不适用于研究区灌溉水源的分类。

3.3.2 基于MODIS数据的分类结果验证与分析

基于2013年验证样本分别与监督分类和非监督分类结果进行比对验证分类精度，结果见表2。采用监督分类，井渠结合灌溉的分类准确度最高，分类精度为 75%，但渠水灌溉和井水灌溉的分类精度很低，并且城镇和居民点的分类精度也不理想。除去城镇和居民点样本，不同水源灌溉面积的总分类精度为57.14%。由于2013年的样本点数据局限在获嘉县，并且样本点太少，影响分类结果的验证效果。城镇和居民点验证结果较差原因可能为MODIS数据的空间精度太低，导致混合像元问题比较严重，使得城镇和农村居民点的曲线特征无法清晰地表现出来。

在非监督分类的结果中，井水灌溉类型的分类结果最好，渠水灌溉类型的分类结果较好，井渠结合灌溉类型的分类最不理想。除去城镇和居民点样本，不同灌溉水源的总分类精度为 64.28%。分类结果不理想的原因可能为：非监督分类方法不适用或图例数据的空间分辨率太低。由于精度较高的Landsat8数据的非监督分类结果较差，可以认为非监督分类方法不适用于灌区尺度的灌溉水源分类。

表1 Landsat8分类结果精度评价表Table 1 Accuracy evaluation of landsat8 classification results

表2 MODIS分类结果精度评价表Table 2 Accuracy evaluation of MODIS classification results

4 讨 论

在 Landsat8数据的监督分类结果中，渠水灌溉和井水灌溉耕地的NDVI时间序列曲线变化与理论假设基本一致，即井灌作物的NDVI峰值低于渠灌作物，渠灌作物的NDVI峰值低于井渠结合灌溉作物。图例数据的监督分类结果中，井灌作物的NDVI值最小，但渠灌和井渠结合灌溉的NDVI值差别不大，并且同种类别的NDVI时间序列曲线的变化差距较大，对分类结果有一定的影响。虽然不同水源灌溉作物的NDVI数值表现出差异，但不同水源的NDVI值的变化速率和变化时间的差异不明显，还需要进一步的研究分析。

本文采用的卫星数据中，图例数据的时间分辨率高，但空间分辨率低；而 Landsat8数据的空间分辨率高，时间分辨率低，数据精度对不同灌溉水源的分类有很大的影响。在灌区尺度，数据的空间分辨率对分类精度的影响大于时间分辨率。采用时间和空间分辨率均较高的数据预期能进一步提高灌区灌溉面积的分类精度。本研究没有排除降水的影响，分类结果存在一定的误差。人民胜利渠灌区的情况复杂，不同水源灌溉的耕地交错分布，对分类精度也有一定的影响。此外，人民胜利渠灌区有少部分耕地使用河水灌溉，使得这部分耕地的分类结果显示为渠水灌溉。不同水源灌溉的耕地的典型NDVI时间序列曲线是进行监督分类的关键。NDVI时间序列曲线特征的影响因素较多，曲线特征与灌溉水源的关系还有待进一步研究。本研究提取的典型NDVI曲线仅为一年的特征曲线，该曲线在其他年份的适用性还需进一步检验。此外，对于与人民胜利渠灌区的灌溉情况和气候差异较大的灌区，本研究采用的研究方法的适用性还需探讨。如何快速获取可适用于不同区域和不同年份的特征明显的典型曲线是进一步研究的重点。

5 结 论

本文以人民胜利渠灌区为例，基于 Landsat8遥感数据和图例数据获取NDVI时序曲线，应用监督分类和非监督分类的方法提取研究区 2017年和 2012年不同水源的灌溉面积。研究发现：①非监督分类方法并不适用于灌区尺度的灌溉水源分类。②图例500 m空间分辨率的影像数据造成的混合像元问题较重，在很大程度上影响了分类的精度，监督分类精度仅为57.14%。③Landsat8影像数据的空间分辨率较高，基于NDVI时序数据采用监督分类的分类效果较好，分类精度可达 73.58%。④典型样本曲线在一定程度上影响分类的结果，获得典型样本曲线是采用监督分类方法的重点，也是提高灌溉水源分类精度的关键。