土壤湿度作为陆面水资源形成、转化及消耗过程中的基本参数, 对植被生长、农业生产和区域生态系统健康运转具有重要影响^[1]。遥感技术的发展, 特别是近年来不同传感器、多时相、多波段、高光谱技术的发展, 为大范围、及时、准确土壤湿度遥感监测提供了可能^[2]。自20世纪80年代起, 就有很多学者提出了基于不同遥感数据源及波段信息的土壤湿度反演模型和方法, 其中最具代表性的有热惯量法、微波遥感监测法、热红外遥感监测法、温度-植被综合指数法和光谱特征空间法等^[3-9]。1977年Richardson和Wiegand建立Nir-Red光谱特征空间提出基于土壤背景线的垂直植被指数后, 从遥感光谱特征空间提取土壤水分信息建立土壤湿度遥感反演模型的方法被学者们持续关注, 学者们先后提出了归一化植被指数(Normalized difference vegetation index, NDVI)、增强型植被指数(Extreme value index, EVI)、温度植被干旱指数(Temperature vegetation dryness index, TVDI)、垂直干旱植被指数(Modified perpendicular drought index, MPDI)等监测作物生长期的土壤湿度变化状况, 这种通过遥感影像计算各类指数间接反映土壤含水量的模拟方法应用最为广泛^[10-16]。前述研究虽然取得了不错的成果, 但所采用的遥感数据源多为TM、ETM+甚至中分辨率的MODIS影像, 而基于高分辨率遥感数据的研究却很少, 更没有研究对比不同指数/传感器在干旱半干旱农业区的反演效果和适用性。因此, 本文结合高分辨率遥感数据(GF-1 WFV和Landsat8 OLI)高空间、高时间分辨率和宽覆盖的特点, 以新疆维吾尔族自治区的阿克苏河流域为研究区, 重点探讨不同高分辨率遥感数据源、不同指数在流域尺度土壤湿度反演中的可行性和适用性, 定量监测阿克苏河流域土壤湿度的时空变化信息, 以协调灌区各分闸口的灌溉定额, 促进流域水资源的合理分配和高效利用, 拓展国产高分系列遥感数据在精准农业、农业信息定量提取等方面的应用范围, 为“天地网一体化”的现代农业信息获取和农情信息遥感监测提供理论基础。

1 研究区概况

阿克苏河流域位于新疆维吾尔自治区西部、塔里木盆地西北部, 地处塔里木河上游平原地区, 地理位置为40°—41°35′ N, 78°47′—82°43′ E。具有明显的大陆性季风气候, 干燥少雨, 多年平均降水量约45 mm, 20 m²水面蒸发量1500 mm, 年太阳辐射总量达到6000 MJ/m², 农业全部依靠地表水渠系灌溉和地下水井抽水灌溉。阿克苏流域绿洲是新疆主要的灌溉绿洲农业区和重要的粮食、瓜果生产基地之一, 也是中国重要的棉花生产区。近年来, 随着阿克苏地区社会经济的发展, 流域内地下水位不断降低, 农业生产受到一定影响, 胡杨林等重要的生态系统也逐渐退化。因此, 借助高分辨率遥感数据实现大面积高精度的土壤湿度的实时监测对流域水资源利用及生态安全评价具有极大的现实意义。

2 数据来源与研究方法 2.1 数据来源

选取阿克苏流域的阿瓦提县(灌溉农业为主)、阿拉尔市(灌溉农业为主)和温宿县(林果业为主)3个实验区(图 1), 实测数据为0—10、10—20、20—30 cm深度层的土壤体积含水率(相对值), 采用托普云农TZS-2X-G土壤水分温度速测仪测量(理论相对误差小于3%), 采样时间为2016年7月14日—7月20日, 采样时点前后一个星期天气状况稳定。采样时根据不同土地覆盖类型, 选取面积大于30 m×30 m的规则地块, 在地块中间位置选择2 m×2 m的样方, 样方内随机测量3次取平均值, 同步记录GPS样点坐标、土地覆盖类型和植被密集程度等信息；共采集102个样点(土壤湿度实测数据306个), 其中阿瓦提采样区、阿拉尔采样区和温宿采样区的实测样点分别为45、30个和27个；样点覆盖类型包括小麦、棉花、水稻、果园、林地和裸地, 以农田植被为主(图 1)。

采用两种高分辨遥感数据, 一是GF-1 WFV多光谱影像, 空间分辨率为16 m, 成像时间2016年6月13日、7月6日、7月22日和8月12日, 主要利用红光波段(Band 3)和近红外波段(Band 4)数据；二是Landsat8 OLI遥感影像, 空间分辨率30 m, 成像时间2016年7月18日, 主要利用红光波段(Band 4)和近红外波段(Band 5)数据。在分析之前, 对各遥感影像进行了预处理, 然后借助ArcGIS10.1将实测样点数据与预处理后的影像进行地理坐标配准。

2.2 研究方法 2.2.1 垂直干旱指数

垂直干旱指数(Perpendicular drought index, PDI)是根据植被-土壤二元组合在红光-近红外二维空间光谱分布变化规律而提出的一种土壤水分反演指数, 它在实践应用中取得了较好效果^{[8-10, 17]}。PDI的计算公式为：

(1)

式中, R_nir和R_red为遥感影像中近红外和红光波段反射率, 分别对应本研究中GF-1 WFV影像的第三、第四波段和Landsat8 OLI影像的第四、第五波段反射率；M为土壤线斜率。

借助ENVI 5.1, 从经预处理的研究区GF-1 WFV遥感影像中提取每个土壤湿度实测样点像元在近红外和红光波段的反射率, 将102个实测样点的反射率在Nir-Red构成的二维光谱特征空间中进行离散化, 然后进行趋势线拟合得到研究区的土壤线方程如公式2。

(2)

根据土壤线的定义可以确定研究区土壤线斜率M为1.2381, 土壤线在纵坐标上的截距I为0.0367。由于同一地块的土壤线差异不大, 基于Landsat8 OLI的土壤湿度反演模型构建时采用同样的土壤线。

2.2.2 改进型垂直干旱指数

PDI没有考虑地表植被覆盖对红光和近红外波段的强散射作用, 因此主要适用于低植被覆盖或裸土地区土壤湿度的遥感反演^{[8-10, 18]}。针对此局限性, 引入植被覆盖度f_v对在Nir-Red光谱特征空间的混合像元进行分解, 克服植被对红光和近红外波段的散射影响, 获取与土壤湿度有关的纯土壤像元反射率, 得到MPDI^[18], 计算公式如下：

(3)

式中, R_{red, v}、R_nir,v为植被在Red和Nir波段的反射率；f_v为植被覆盖度。f_v是指植被(包括叶、茎、枝)在地表的垂直投影面积占统计区总面积的百分比, 本研究主要利用f_v来克服遥感影像中混合像元对土壤湿度光谱信息的影响, 公式如下：

(4)

式中, NDVI_v、NDVI_s分别代表纯植被和裸土的归一化植被指数。借助ENVI 5.1中的Band math工具, 通过红光波段和近红外波段的反射率, 可以计算获取各时期GF-1 WFV和Landsat8 OLI遥感影像的NDVI值；由于研究区地表覆盖复杂, 计算得到的NDVI最大/最小值可能存在误差, 拟取累积概率为5%和95%的NDVI值作为最小值和最大值。

2.2.3 植被调整垂直干旱指数

考虑到覆盖饱和的影响, 引入垂直植被指数(PVI)代替f_v, 作为植被覆盖表征量, 在PVI-PDI二维空间对PDI模型进行修正, 提出了适用于高植被覆盖区土壤湿度反演的植被调整垂直干旱指数(Vegetation adjusted perpendicular drought index, VAPDI)^[8-10]。

如图 2所示, 三角形ABC内所有土壤湿度等值线近似为直线, 且交于A点。在PVI=0的裸土区, 三角形ABC内任一点E的PDI可以近似用土壤湿度等值线AE与横坐标轴的交点F的PDI代替, OF的长度即为E点修正后的PDI, 即PVI趋于0时, VAPDI等于PDI。

根据三角形相似原理得到任意X点的VAPDI计算公式如下：

(5)

其中, 垂直植被指数PVI的计算公式如下：

(6)

式中, I是土壤线表达式的截距。借助ENVI 5.1中的Band math工具, 获取各时期GF-1 WFV和Landsat8 OLI遥感影像的PVI值。

3 结果与分析 3.1 土壤湿度反演模型构建 3.1.1 土壤实测湿度的描述性统计

以研究区102个表层(0—10 cm)土壤湿度实测数据为例进行描述性统计分析, 结果见表 1。各植被类型的表层土壤湿度均值和中值接近, 说明研究区土壤湿度整体分布较为均匀；裸土的中值和均值存在差异, 且变异系数达到43.8%, 说明研究区裸土的干、湿度存在一定差异；土壤湿度最大值和最小值分别为0.485和0.088, 对应的覆盖类型分别为水稻和裸土；102个实测样点的土壤湿度平均值和标准差分别为0.225和0.081。

3.1.2 土壤湿度反演模型构建

依据各采样点的土地利用类型、植被覆盖度、地形、坡度坡向等具体指标, 从102个采样点中选取68个作为建模样本集, 剩余的34个样点作为反演模型精度验证和评价的验证样本集。通过ENVI 5.1中的Band math工具, 计算获取GF-1 WFV(2016年7月22日)与Landsat8 OLI(2016年7月18日)影像的f_v、PVI两个参数和PDI、MPDI和VAPDI, 借助SPSS 19.0对68个建模样本3个深度层的土壤湿度实测数据和对应的PDI、MPDI和VAPDI分别进行拟合, 发现线性拟合效果最好, 汇总回归模型见表 2。

从表 2可以看出, 拟合模型均为负相关, 除了基于GF-1 WFV影像的20—30 cm深度层的PDI、MPDI拟合模型和基于Landsat8 OLI影像的20—30 cm深度层的PDI拟合模型通过0.05的显著性检验外, 其余模型均通过0.01的显著性检验, 具备统计学意义；尤其是各指数均与0—10 cm深度层的土壤湿度相关性最强, 平均R²达到了0.68, 说明基于光学遥感影像近红外和红光波段构建的遥感反演指数对近地表层土壤湿度信息具有较强的敏感性, 能够模拟和监测更大范围的土壤湿度的空间变化, 但对地下较深层次的土壤湿度反演精度略低。

从各回归模型的决定系数R²来看, MPDI、VAPDI的拟合效果要明显优于PDI, 这是因为在Nir-Red二维光谱特征空间中, 各像元的反射率由土壤、植被甚至其他地物信息共同决定, PDI指数没有考虑到植被覆盖对土壤光谱信息的影响, 因此无法完全反映出表层土壤湿度的实际水平。而MPDI和VAPDI分别利用了不同的植被覆盖表征量对混合像元的光谱信息进行修正, 因而所表达的土壤湿度信息更精确, 与土壤湿度实测值之间的相关性高于PDI。

从表 2来看, 对于WFV, PDI与土壤湿度间的相关性随深度增加而明显减小, 而OLI的PDI与3个深度层土壤湿度间的相关性差异不大, 但同样都是在20—30 cm深度层的相关性最小, 标准误差也相对较大。说明对于WFV和OLI两种遥感影像, 当深度超过20 cm, PDI对土壤湿度的敏感度不高, 难以准确地反映该深度层的实际土壤湿度状况。MPDI与土壤湿度的相关性也随着土壤深度增加而减小, 基于WFV和OLI的MPDI均与0—10 cm土壤湿度相关性最强, R²达到了0.7, 说明基于近红外和红光波段信息构建的反演指数对土壤湿度的敏感性随深度增加而降低。在3个深度土层上, 基于WFV与OLI的MPDI和土壤湿度实测值的相关系数水平基本相当。两种遥感影像的VAPDI与3个深度层土壤湿度之间都存在线性负相关；其中在0—10 cm深度层的相关性最强, R²分别达到了0.7469和0.7869, 在10—20 cm和20—30 cm深度层相关性相对较低, 但高于其他两个指数；就传感器而言, 基于WFV影像的回归拟合模型的各项系数均优于OLI影像。

3.2 反演模型的精度评价

将基于两种传感器的PDI、MPDI和VAPDI拟合模型得到的不同深度层的土壤湿度反演值, 分别与对应的34个验证样点的土壤湿度实测值进行验证分析, 并分别计算验证样点土壤湿度反演值与实测值的决定系数(R²)、平均绝对误差(MAE)、平均相对误差(MRE)和均方根误差(RMSE)等指标值, 来验证和定量评价反演模型的精度(表 3)。

由表 3可知, 从各精度评价指标值来看, 3个土壤深度层中, 基于WFV和OLI两种传感器的PDI、MPDI、VAPDI土壤湿度反演模型均在0—10 cm处反演精度最高, 说明光学遥感影像更适合于表层土壤湿度的反演, 而对较深层次土壤湿度信息的敏感性较弱, 原因是近红外和红光的穿透能力较弱, 以反映近地表的光谱信息为主。

对比3种指数的反演结果来看, 基于WFV和OLI两种传感器的VAPDI土壤湿度反演模型在0—10 cm深度层反演结果的精度均明显高于PDI和MPDI；在敏感性较差的10—20 cm和20—30 cm处, 3种指数的精度评价指标虽然相差不大, 但VAPDI和MPDI的反演结果稍优于PDI。从GF-1 WFV和Landsat8 OLI遥感数据源的反演结果精度来看, 在同一深度层上, 对于不同的反演指数, 两种影像的反演精度各有优劣, 其中基于WFV的PDI和VAPDI反演结果的R²、MAE、MRE和RMSE均优于OLI影像；而利用MPDI构建模型时, 二者的反演效果基本一致, 反演精度无明显差别。

3.3 表层土壤湿度反演与动态监测

根据前述分析, 拟选择两种遥感影像的VAPDI土壤湿度反演模型, 对阿克苏河流域0—10 cm表土层的土壤湿度进行反演, 通过对实际反演效果的对比分析, 提出推荐的遥感数据源用于土壤湿度大规模动态监测。

3.3.1 表层土壤湿度反演结果分析

以经过预处理的研究区GF-1 WFV(2016.7.22)、Landsat8 OLI(2016.7.18)遥感影像为基础, 在ENVI/IDL中计算出各像元的VAPDI值, 利用构建的土壤湿度反演模型获取每个像元的土壤湿度值, 得到研究区表层土壤湿度的空间分布格局(图 3和图 4), 图中白色部分为掩膜掉的水体与城区, 红色表示土壤湿度偏低, 由红色到深蓝色土壤湿度逐渐增高。

通过图 3、4可以看出, 两种遥感影像反演得到的阿克苏河流域表层土壤湿度空间格局基本一致, 即靠近水源的区域土壤湿度高, 而远离水源的地区土壤湿度值一般偏低。原因主要是阿克苏河流域自然降水量极小, 流域农业生产需水几乎全部依靠地表渠系和地下水抽水灌溉, 而且靠近河流、水库及湖泊等自然水源的地区一般为农业生产区, 如上游托什干河两侧的乌什县、库玛拉河两侧的温宿县南部、中游的阿克苏河两岸的阿瓦提灌区和阿克苏市中部以及中下游阿克苏河与塔里木河交汇处的阿拉尔灌区, 灌溉渠系密集, 能够长期保障该区域农作物用水需求, 表层土壤湿度相对较高, 均在0.3以上；流域内三大水库胜利水库、上游水库、多浪水库和流域西南部的千鸟湖地区土壤湿度也较高；远离水源的地区多为天然植被或荒漠戈壁区, 水网密度小且缺乏人工灌溉和维护, 受流域干燥气候影响大, 土壤湿度值一般处于0.2以下, 部分区域土壤湿度值甚至接近为0。

进一步对比两种遥感数据源的土壤湿度反演结果可以发现, 图 4中不同区域的土壤湿度的高低层次更加分明, 土壤湿度的空间异质性更加明显, 这主要是GF-1 WFV遥感影像空间分辨率为16 m, 稍高于Landsat8 OLI影像30 m的空间分辨率, 在土壤湿度反演模型精度相当的情况下, 较高的空间分辨率能够表现出更加精细的地表信息, 能够更详细的表达土壤湿度的空间异质性, 能够为实现精确到地块的农作物灌溉和区域生态安全评价提供参考。

3.3.2 表层土壤湿度动态监测

为进一步探讨流域土壤湿度的时空演变规律, 选择2016年6月13日、7月6日、7月22日和8月12日4个时点的GF-1 WFV遥感影像, 对流域0—10 cm深度层土壤湿度的时空动态变化进行反演监测和分析。按前述同样的方法, 得到4个时点的0—10 cm土壤湿度空间分布图(图 5)。

由图 5可知, 研究区土壤湿度较高的区域主要分布在沿河流两侧的乌什县、温宿县南部、阿瓦提县北部、阿克苏市东部及阿拉尔市中部等流域中上游地区, 该区域是阿克苏河流域主要的农业生产区, 灌溉渠系密集, 且6—8月份为农作物的主要生长期, 农田植被覆盖度高, 土壤持水能力强, 土壤湿度高且相对稳定, 对这一类土壤湿度相对稳定的流域中上游农田植被区, 为提高水资源利用率, 减少沟渠水系漫灌带来的水资源损失, 农业需水灌溉可采用滴灌措施；而乌什县北部、温宿县和阿克苏市东部、阿瓦提县西部和南部、柯坪县、沙雅县等流域边缘及下游地区土壤湿度相对较低, 且该地区土地利用类型大多为灌木林、荒草地和裸地, 植被覆盖度低, 土壤湿度受地表蒸散发影响较大, 时空变化相对强烈；对这类流域边缘及下游表层土壤湿度时空变化强烈的非农业生产区, 应加强人工防护林、生态林的种植与维护力度, 防止进一步的土地荒漠化及天然林地、草地的退化, 从而保障阿克苏河流域的农业可持续发展和绿洲生态系统的健康运行。

对4个时点的表层土壤湿度的时空动态监测结果分析发现, 阿克苏河流域表层土壤湿度空间格局稳定但随时间变化特征较为复杂, 研究区降水稀少, 主要依托区域完善的沟渠体系进行地表灌溉, 因而以阿克苏河为轴线的绿洲农业核心区, 土壤湿度保持在35%以上, 特别是当人工开闸放水进行灌溉时, 农田植被区的土壤湿度会明显增强；另外, 在中国气象数据网(http://data.cma.cn/)查阅2016年6—7月的全国24 h降水量分布图, 发现2016年6月28日前后阿克苏河流域北部、2016年7月6—12日全流域尤其是11—12日阿拉尔、沙雅县等地区有过一次持续降水, 尽管日降水量在几毫米, 但对土壤湿度的影响还是较大的, 7月22日反演的土壤湿度也明显增强, 特别是阿拉尔、沙雅县的大部分区域, 这也与植被指数估算土壤湿度存在5—10天的延迟天数吻合^[10]。

4 讨论

目前关于区域土壤湿度的遥感反演研究中, 采用MODIS数据和Landsat TM/ETM数据偏多, 而Landsat8 OLI/TIRS系列数据的开放, 为区域土壤水分监测提供了新的数据源^[19]；GF-1遥感卫星作为中国高分辨率对地观测系统国家科技重大专项的首发星, 其WFV传感器具有16 m空间分辨率、800 km幅宽、4天覆盖周期的特点, 在农业监测与评估应用方面具有独特的优势^[20], 使得在小尺度层次上进行大范围土壤湿度空间格局研究成为可能。本研究在前人研究基础上尝试以GF-1 WFV和Landsat8 OLI两种高分辨率遥感影像为数据源, 利用PDI、MPDI、VAPDI对阿克苏河流域进行不同深度层土壤湿度反演与动态监测, 验证GF-1 WFV反演土壤湿度的可行性和获取最适反演深度；而在获取区域土壤湿度后, 可结合从水源到田间输配水的复杂过程构建的流域水资源可达性模型和农作物种植格局、农作物灌溉定额等, 提出格网化或灌区化的灌溉用水量的快速测评技术, 以协调灌区各分闸口的灌溉定额, 对促进流域水资源的合理分配和高效利用具有较强的应用前景。

本研究不同植被指数的反演结果表明, MPDI和VAPDI能够在一定程度上克服混合像元对土壤湿度光谱信息的影响, 反演的精度要比PDI高, 尤其是高植被覆盖度区, 采用垂直植被指数(PVI)修正的VAPDI反演精度最佳, 这与杨学斌等学者^[16]的PDI、MPDI与植被覆盖区实测土壤含水量的R²分别为0.37、0.5355以及吴春雷等学者^[9]的PDI、VAPDI与裸土、麦茬、土豆、豇豆4种植被覆盖类型土壤实测含水率的R²分别为0.630、0.504、0.571、0.543和0.599、0.523、0.602、0.585的研究结论相似。这主要是由于PDI指数在构建拟合模型时, 将遥感影像中各像元反射率视为单一裸土反射率, 忽视了混合像元的影响, 土壤湿度信息区分能力较差, 因而从遥感机理上来说其反演结果的误差较大, 主要适合于裸土或者植被稀疏的地区。而经f_v和PVI修正的MPDI和VAPDI建立的土壤湿度反演模型, 通过对遥感影像中的混合像元进行不同程度的分解, 使用植被和裸土像元二者的综合反射率进行拟合, 就能更加精确地反演出表层土壤的湿度信息, 反演结果更接近实测值, 也更适用于植被覆盖度较高的区域。PVI在表征植被覆盖信息时能较好地消除土壤背景的影响, 且相比f_v不易出现植被覆盖饱和现象, 更能精确地表达植被覆盖信息, 因此在高植被覆盖度区域, VAPDI比MPDI更具有优势。

在土壤-植被-大气系统中, 土壤湿度的光谱信息不仅会受到植被覆盖的影响, 还受大气状况、传感器差异、土壤背景等复杂环境因素的影响, 因此在后续的研究中, 还需要进一步对不同遥感传感器、不同气候区的土壤湿度反演效果进行对比分析, 提取合适的因子对现有土壤湿度反演指数进行修正, 特别是融合温度等参数进行协同反演, 以提高模型的反演精度, 从而找到一套适用范围广、精确度高的土壤湿度遥感反演方法。

5 结论

1) 对于GF-1 WFV和Landsat8 OLI遥感影像, PDI、MPDI和VAPDI与研究区0—10、10—20、20—30 cm深度层的土壤湿度实测值之间均存在一定的线性负相关关系, 其中与0—10 cm表层土壤湿度的相关性最强, 平均R²达到了0.68, 表明基于光学遥感影像近红外和红光波段反射率构建的土壤湿度反演指数对近地表层土壤湿度信息更敏感, 但对地下较深层次的土壤湿度反演精度略低。

2) 相比PDI指数模型, MPDI与VAPDI反演精度更高, 同时具有快速、高效的特点；受植被覆盖饱和的影响, 在高植被覆盖区宜采用VAPDI。

3) 两种遥感影像反演得到研究区表层土壤湿度空间分布图的总体格局一致, 但空间分辨率较高的GF-1 WFV遥感影像的反演结果相比Landsat8 OLI更精细, 得到土壤湿度空间分布结果更有层次性, 能够反映出相近地块土壤湿度的异质状况。