基于VSWI 和TVDI 差异的河套灌区沈乌灌域耕地灌溉面积遥感监测

2020-09-05李瑞平王思楠范雷雷于明辉

灌溉排水学报 2020年8期

田鑫，李瑞平，王思楠，范雷雷，于明辉

（1.内蒙古农业大学水利与土木建筑工程学院，呼和浩特 010018；2.威海市城市规划技术服务中心有限公司，山东威海 264200；3.内蒙古自治区农牧业大数据研究与应用重点实验室，呼和浩特 010018；4.威海市国土资源地理信息中心，山东威海 264200）

0 引言

【研究意义】我国是个农业发达但水资源匮乏的国家，伴随着我国农业的发展，水资源合理的分配已成为我国农业得以健康发展的基础，而灌溉面积监测作为水资源合理配置的重要依据也越来越受到重视[1-2]。沈乌灌域作为河套灌区的重要组成部分，区域内干旱少雨，蒸发量远大于降水量，水资源紧缺、有效灌溉面积增速减缓缓慢的现状日益严重，因此通过对沈乌灌域的灌溉面积进行监测，为沈乌灌域水资源合理配置提供理论依据十分必要。

【研究进展】20 世纪的90 年代初期，杜文才等[3]利用卫星遥感技术对我国宁夏中宁地区的灌溉情况展开了调查和分析，并绘制出了中宁地区的水资源、农业灌溉设施分布图等。随后，水利部将河南省的灌区作为研究对象，结合美国陆地卫星4～5 号专题制图仪获取的多波段扫描影像以及河南省现有的数据资料等，进行实地的勘察，基本实现了利用遥感技术实现对有效灌溉面积的测量[4]。许迪等[5]对我国山东省簸箕李灌区的农业灌溉系统展开了分析和调查，建立了识别系统，对于灌溉面积测量方法的开发起到了积极的作用。21 世纪初期，我国开展了宁夏引黄灌溉面积及作物种植结构遥感调查项目研究，结合遥感技术、地理信息系统和全球定位系统，实现了对宁夏引黄灌区的灌溉面积、农作物分布、土地结构等数据的收集和分析。易珍言等[6]、沈静[7]分别采用对PDI、MPDI 差值的阈值进行确定的方法，结合少量实测数据，对灌溉面积进行了提取，取得了较好的成果，但均需提前知晓灌溉行为发生的准确时间。

【切入点】基于以上可以发现，遥感技术在灌区灌溉面积监测上具有广泛的应用，基于指数差异阈值的灌溉面积监测方法也取得了一定的进展，但植被供水指数VSWI、温度植被干旱指数TVDI 等在灌溉面积提取方面的应用较少，多用于研究区的旱情监测等，因此，本文尝试采用基于Landsat-8 卫星数据的植被供水指数VSWI、温度植被干旱指数TVDI 结合地表温度LST 变化对研究区的灌溉面积进行提取。

【拟解决的关键问题】探索VSWI、TVDI 在灌溉面积遥感监测模型创建中作为主要参数的可行性。

1 研究区概况与数据来源

1.1 研究区概况

沈乌灌域位于巴彦淖尔市的西南部，三盛公枢纽西北部，106°08′E—107°10′E，40°08′N—40°57′N，沈乌灌域是内蒙古河套灌区的重要组成部分。灌域东南至黄河与鄂尔多斯市隔河相望，西南与阿拉善盟毗邻，东北部与解放闸灌域接壤，西北靠阴山山脉，东西长92 km，南北宽65 km。承担着中国林业科学研究院沙林中心，巴彦淖尔市6 个国营农场，磴口县5个苏木、镇、办事处，杭锦后旗1 个镇，阿拉善盟1个国营农场及鄂尔多斯市杭锦旗巴拉亥镇部分耕地的灌排任务，灌域内土壤类型以壤质砂土和砂质壤土居多，并含有少部分的砂土，沈乌灌域具有明显的大陆性气候特点，冬季严寒，夏季短促。气候干旱少雨，蒸发强烈，年平均降雨量在130～215 mm 之间，年平均蒸发量在2 100～2 400 mm 之间。蒸发量与降水量的严重失衡造成沈乌灌域对引黄灌溉的绝对依赖。

1.2 数据来源

本次研究选用的数据为Landsat-8 卫星影像数据，共选用了3 景影像，根据往年灌溉时间，2016 年7—9月，选取每月云量较少的影像，用于研究区8 月底第5 次灌溉面积的提取，详细情况见表1。

表1 研究所用Landsat OLI/TIRS 影像数据 Table1 Landsat OLI/TIRS image data for research institute

收集到耕地野外调查样本379 份，用于模型构建及土地利用类型图的绘制，其中玉米145 个，葵花177 个，瓜类57 个。非耕地采样点采用随机采样的方式，共采集沙地样本点34 个，建设用地采样点25个，林地采样点26 个，同时于8 月初对研究区已灌溉区域进行采样，共计182 个采样点，用于灌溉面积反演结果的检验。

2 模型构建

本文以内蒙古河套灌区沈乌灌域为研究区，基于Landsat-8 数据产品对灌溉面积进行监测，构建灌溉面积遥感监测模型。通过对研究区归一化植被指数NDVI、地表温度LST、植被供水指数VSWI、温度植被干旱指数TVDI 的反演，计算相邻影像的VSWI、TVDI 差值，进而计算耕地实测样本点的VSWI、TVDI差值的均值，通过差值的均值来对每幅影像的灌溉行为进行分析，构建基于VSWI 差异和TVDI 差异的灌溉面积监测模型，结合基于SEE5.0 决策树分类法获取的耕地边界矢量图，对耕地的灌溉面积进行提取。

2.1 归一化植被指数NDVI 反演

遥感光谱数据在计算分析生物量、植被生长状况上都需要借助植被指数，是一项无量纲参数，在光谱数据分析中起到了至关重要的作用。没有任何条件的辅助以及约束，植被状态都可以通过光谱信息表达出来，同时能够对植被的生长活力、生物量、覆盖面积等等这些信息进行定性、定量的评价和分析[8]。1978年，Deering[9]首次提出了RVI（简化比值植被指数），通过非线性的归一化处理后，能够得到NDVI（归一化植被指数），并将其约束在[-1，1]范围内，那么可以得到公式：

式中：ρNIR近红外波段的真实地表反射率；ρR为红光波段的真实地表反射率；NDVI 就是RVI 通过非线性归一化处理后的结果。

2.2 地表温度反演

地表温度是当前我们评估大气层与陆地表层之间的水分以及能量平衡的关键性指标之一，同时也能够反映出土壤墒值、评估作物的产量等，因此对于遥感模型的建立是非常关键的，此外对于土壤的盐碱化、保水能力、气候变化等研究都提供了重要的参考价值。遥感技术的发展推动了地表温度遥感反演算法的精确度和准确度。目前，大范围的地表温度进行监控常常会用到热红外遥感。地表温度反演的算法主要有3种：单窗算法、单通道算法及劈窗算法[10-12]，本研究采用劈窗算法，如式（2）所示：

根据研究区往年的灌溉时间，8 月初的第5 次灌溉和10 月底的秋灌，灌溉作物种类基本涵盖研究区主要农作物，但10 月底前后卫星影像云量普遍较高，因此选取7 月24 日，8 月25 日，9 月16 日3 景影像对地表温度分布情况进行反演，反演结果如图1 所示根据耕地、林地、建设用地、沙地的实测点坐标，计算7—8 月、8—9 月的地表温差LST7-LST8、LST8-LST9，根据表2 可以看出，林地、建设用地、沙地的地表温度7—9 月始终呈下降趋势，而耕地7—8 月地表温度骤降，降幅较其余3 种土地类型更为明显，8—9 月耕地地表温度小幅回升。

根据以上分析，造成耕地温度变化的原因是灌溉，进一步分析，耕地7—8 月的温度降低若是降雨引起，首先其温度降幅应低于或持平于其他土地利用类型，其次其地表温度变化趋势应与其余3 种土地利用类型基本相同。

表2 不同土地类型地表温度时间变化 Table 2 Time variation of surface temperature in different land types ℃

图1 沈乌灌域地表温度反演结果 Fig.1 Inversion results of surface temperature in Shenwu irrigation area

2.3 VSWI、TVDI 差异模型构建

Carlson 等[13]首次提出了植被供水指数（Vegetation Supply Water Index，VSWI）的概念，其是在地表温度和植被指数的基础上建立的一种监控模型，能够实现对植被的干旱情况进行观察。具体的表达式为：

式中：地表温度以TS表示，单位为k；归一化植被指数为NDVI；图像增强系数为B。

Sandholt 等[14]首次提出了温度植被干旱指数(Temperature Vegetation Dryness Index，TVDI)的概念，其是在Ts-NDVI 特征空间基础上发展而来的。这种干旱指数能够计算出植被指数以及地表温度，其数据来源于遥感影像，为此可以有效地避免出现土壤水分状态监控中的问题，大大降低了植被覆盖度的影响程度，进而提高监测数据的准确度，为此对于地形复杂、面积较大区域的干旱监测具有很好的适应性。计算式（4）如示：

式中：湿边方程以Ts-min表示；干边方程以Ts-max表示；地表温度以Ts表示。TVDI 的值域为[0，1]。

植被供水指数和温度植被干旱指数可以描述植被的供水情况和土壤的含水率，植被供水指数越大，植被的供水情况越好，温度植被干旱指数越小，土壤含水率越高，这2 种情况都能证明灌溉行为的发生，因此可以利用以上2 种指数对灌溉面积进行监测，对于灌溉面积的提取，关键在于确定灌溉后的影像，基于以上分析可以发现，灌溉后影像的植被供水指数VSWI 会高于灌溉前的影像。而TVDI 会低于灌溉前的影像，以此可以得出灌溉行为与上述2 种指数的关系，如式（5）和式（6）：

式中：Iv、IT分别表示研究区某个像元处的灌溉程度；VSWIt1、TVDIt1为前一幅影像的植被供水指数和温度植被干旱指数；VSWIt2、TVDIt2为后一幅影像的植被供水指数和温度植被干旱指数。

在排除降水影响的情况下，对于植被供水指数，如果后一幅影像某一处像元的植被供水指数大于前一幅影像，即Iv小于0，可以认为在该像元处，发生了灌溉行为，且Iv值越小，灌溉行为越明显。对于温度植被干旱指数来说，如果后一幅影像某一处像元的温度植被干旱指数小于前一幅影像，即IT大于零，则可以认为在该像元处发生了灌溉行为。且IT值越大，灌溉行为越明显。

按照所收获到的地表温度的反演结果以及归一化植被指数的运算数据，以VSWI 运算公式为基准，对B 进行取值，即B=10 000，在运算之后得出2016年7 月24 日、8 月25 日、9 月16 日的植被供水指数，根据式（3），对各幅影像的VSWI 进行计算，并根据式（5），计算时间相邻影像的VSWI 差值，根据耕地采样点经纬度数据，得到耕地采样点VSWI 差值的均值，结果见表3。从表3 可以看出，8 月耕地采样点的植被供水指数比7 月同比高12.5，相比9 月同比高8.7，基于以上分析，结合本文构建的基于VSWI差异的灌溉面积监测模型可以看出，8 月的灌溉行为十分明显，可以作为灌溉面积提取的依据。

表3 耕地采样点VSWI 差值均值 Table3 Mean value of VSWI difference of cultivated land measured points

故而，按照国家气象局气候司所发布的相应规范对干旱等级进行划分，即重旱、中旱、轻旱、正常、湿润依次为10 cm 土壤相对湿度<40%、40%～50%、50%～60%、60%～80%、>80%；而基于VSWI 则依次为0～5、5～10、10～15、15～25、25～35。获取了研究区7—9 月的旱情分布情况，见图2。从图2 同样可以看出，8 月耕地的植被供水指数情况确实相比7月和9 月明显提高。在获取沈乌灌域研究区地表温度反演结果及NDVI 运算结果之后，在对研究区TVDI进行运算的时候选用了以ENVI 感图像处理软件为基础二次开发的TVDI 功能模块。按照NDVI 对应的地表温度，从起最高、最小值点对干、湿边进行拟合[15-17]，从而得出该区域Landsat-8 数据7 月24 日、8 月25日、9 月16 日的相应拟合曲线，见图3。

图2 基于VSWI 的旱情分布 Fig.2 Drought distribution map based on VSWI

图3 沈乌灌域TS-NDVI 特征空间拟合曲线Fig.3 Ts-NDVI characteristic space fitting curve in Shenwu irrigation area

基于干旱等级划分标准把研究区TVDI 分成了湿润、正常、轻旱、中旱、重旱5 个等级，即：0～0.2、0.2～0.4、0.4～0.6、0.6～0.8、0.8～1.0[18]，图4 为详尽的结果。利用式（4）对3 幅影像的温度植被干旱指数进行计算，并根据野外耕地实测点，计算相邻影像实测点的TVDI 差值的均值。从表4 可以看出，8 月实测点的TVDI 均值小于7 月和9 月，根据基于TVDI差异的灌溉面积监测模型可以看出，8 月的灌溉行为较为明显，因此选取8 月的影像对灌溉面积进行提取。

表4 耕地实测点TVDI 差值均值 Table4 Mean value of TVDI difference between the measured points of cultivated land

3 结果与分析

3.1 土地利用分类

综合Landsat-8 影像的光谱特征，植被指数，ISODATA 分类结果作为样本集的属性值，将研究区5种主要的地物类型作为训练数据集，基于SEE5.0 的决策树分类方法进行土地利用分类，其对于耕地、林地、沙地、建设用地的分类精度分别为90.0%、88.5%、88.2%、88.0%，总精度为88.7%，其土地利用类型分类结果如图5 所示，进而得到沈乌灌域耕地矢量边界。

图5 沈乌灌域土地利用类型 Fig.5 Land use type map of Shenwu irrigation area

图6 基于耕地矢量边界的沈乌灌域灌溉面积分布 Fig.6 Distribution map of irrigation area in Shenwu irrigation area based on cultivated land vector boundary

表5 灌溉面积监测结果精度评定 Table5 Accuracy evaluation of irrigation area monitoring results

3.2 耕地灌溉面积提取

根据耕地矢量边界，将基于VSWI 差异模型、TVDI 差异模型所获取的研究区第5 次灌溉面积进行提取，结果如图6 所示可以看出，二者具有较高的一致性，提取的灌溉面积区域基本相同，同时根据统计数据对2 种模型的精度进行对比，结果如表5。根据反演已灌溉点位占灌溉采样点的百分比，VSWI 和TVDI 的模型精度分别为85.3%和89.7%。

4 讨论

通过指数阈值选取对灌溉面积进行提取，关键在于对于灌溉行为发生时间的确认及阈值的确定，相较之前的成果，本次研究通过地表温度的变化确定灌溉行为发生的时间，针对河套灌区等存在大量沙地的区域可取的较好的成果。今后的TVDI 模型创设，会用到一些参数，譬如地表温度、植被信息等[19-22]，在遥感参数繁杂的推算期间或许会对一定的经验和半经验公式加以运用，关键参数反演结果在地表复杂度较高的区域内会出现较大的偏差，从而使TVDI 监测旱情情况的精准度无法得到有效地确保。因而，在后期探究中，需要将各个模型参数的适应范畴和运算条件理清，结合多年数据，对模型参数进行进一步优化，为定量化遥感分析打好基础。