基于多时相高分一号影像的土壤湿度反演

杨丹阳, 严颂华, 杨永立, 田 茂

(1.武汉科技大学信息科学与工程学院，武汉 430080；2.武汉晴川学院北斗学院，武汉 430204；3.武汉大学宇航科学与技术研究院，武汉 430072)

土壤湿度又称土壤含水量，是量化陆地及大气能量交换的关键参数，精准地监测土壤湿度的变化对农业、气候、水文、陆地生态环境等领域的应用具有重要意义[1-2]。相对于传统的土壤湿度监测方法，光学遥感技术具有高效率、低成本、大范围、高分辨率、快速获取、动态监测能力强等优点[3]，借助光学遥感数据反演土壤湿度引起了众多学者的关注与研究，同时也是目前遥感技术应用的重要发展方向之一[4-5]。

基于可见光-近红外波段的干旱指数反演土壤湿度的方法被陆续提出，先后有：归一化植被指数(normalized difference vegetation index,NDVI)、基于两波段增强型植被指数(two-band enhanced vegetation index,EVI2)[6]、垂直干旱指数(perpendicular drought index,PDI)[7]、修正型垂直干旱指数(modified perpendicular drought index,MPDI)[8]、植被调整垂直干旱指数(vegetation adjusted perpendicular drought index,VAPDI)[9]等。近年来张月等[10]基于NDVI、EVI指数监测藏北牧区土壤湿度，发现二者均呈现较强的相关性；李喆等[11]以湖北漳河灌区为研究区域，验证了PDI模型在湿润地区监测的可行性；杨学斌等[12]在内蒙古干旱监测中发现MPDI模型具有更高的监测精度；吴春雷等[9]提出VAPDI指数，并验证了该模型在农田植被区土壤水分监测方面的适用性。

但是先前的土壤湿度研究大多使用单一时间点的Landsat TM、ETM+和MODIS等中低分辨率光学遥感影像作为数据源。高分一号(GF-1)光学遥感卫星作为中国首颗自主研制的高分辨率对地观测卫星[13]，具有重访周期短(4 d)、空间分辨率高(16 m)、免费获取数据等优势，在农作物信息提取等方面取得了显著的成果[14-15]，但用于土壤湿度反演方面的研究正处于起步阶段，尤其是对不同季节、不同地域的土壤湿度监测研究甚少。

因此，利用GF-1的成像优势，选择武汉市江夏区豹澥镇野外观测场作为研究区域，从时间角度出发计算2019—2020年不同时相数据的PDI、基于NDVI的MPDI(MPDIN)、基于EVI2的MPDI(MPDIE)和VAPDI干旱指数，并基于四种指数分别建立土壤湿度反演模型，对其进行精度评估，在此基础上研究该区域的空间格局分布，进一步分析VAPDI指数对土壤湿度的反演能力。

1 研究区域与数据源

1.1 研究区域

研究试验区为豹澥野外观测场，如图1所示。该地区气候温和，雨量充沛，年平均气温16 ℃，冬季(1月)平均气温2～4 ℃，夏季(7月)平均气温 30 ℃ 左右，年平均降雨量1 240 mm，全年无霜期约为260 d。

图1 豹澥野外观测场区Fig.1 Baoxie field observation field

1.2 数据源

遥感数据采用空间分辨率为16 m的GF-1多光谱宽幅覆盖(wide field of view,WFV)影像，有效载荷技术指标如表1所示。主要用到的产品级别为1A级(相对辐射校正产品)TIFF格式影像，主要波段数据为红光波段和近红外波段。为了减少大气和云层对实验结果的干扰，同时考虑季节因素对土壤湿度空间格局的影响，经过筛选得到高分数据成像时间为：2019年12月—2020年6月，共10景。此遥感数据可以从中国资源卫星应用中心下载。

表1 GF-1 WFV有效载荷技术指标Table 1 Technical specifications of GF-1 WFV

为了构建土壤湿度反演模型并验证土壤湿度反演结果，还需要用到与地面同步的表层土壤湿度实测数据。因此，采用NB-IoT低功耗、高精度土壤湿度传感器(图2)在观测场建立5个土壤湿度地面观测站点，其中站点D用于建立反演模型，站点A、站点B、站点C、站点E用于验证评估(站点A为无植被覆盖区，站点B、C、E为有植被覆盖区)，以实现豹澥试验区土壤湿度的长期动态连续监测，采集时间与遥感数据成像时间相对应。具体的土壤湿度实测数据来自环境监测数据平台。

图2 土壤湿度传感器Fig.2 Soil moisture sensor

2 研究原理

2.1 土壤线

由于土壤在红光波段(Red)和近红外波段(NIR)的二维光谱特征空间中呈一种线性关系，因此根据这种关系拟合出来的直线称之为土壤线，如图3中直线BC。表达式为

Rnir=aRred+b

(1)

式(1)中：a为土壤线斜率；b为土壤线在纵坐标上的截距；Rnir和Rred分别为近红外波段和红光波段光谱反射率。

图3 垂直干旱指数原理图Fig.3 Schematic diagram of perpendicular drought index

2.2 垂直干旱指数

图3中，通过坐标原点作一条垂直于土壤线BC的法线L，此时任意一点E(Rred，Rnir)到直线L的距离越小代表土壤越湿润，反之亦然。Ghulam等[7]以此提出了一种便捷有效的土壤湿度反演模型——垂直干旱指数(PDI)，计算公式为

(2)

式(2)中：PDI的范围为0～1，PDI越大，对应的土壤湿度越干旱，反之亦然。

2.3 修正型垂直干旱指数

图3中植被覆盖情况不同的A、E、D三点到法线L的距离相等，可以说明这三点的PDI值相等，依据上述的PDI设计原理可以得到这三处土壤湿度值相等。但是A、E、D三处有着不同的植被覆盖度，因此它们的土壤湿度值也不相同，这就说明了PDI模型不具有普遍性。因此引入植被覆盖度因子去除PDI指数没有考虑到的植被部分的影响，建立土壤湿度快速监测模型——修正型垂直干旱指数MPDI[8]提取混合像元中土壤的反射率，计算公式为

(3)

式(3)中：Rv,red和Rv,nir分别为植被端元的红光和近红外波段光谱反射率；fv为植被覆盖度，计算公式为

(4)

式(4)中：VI表示植被指数；VIv和VIs为纯植被和纯裸土的归一化植被指数，分别取累积概率为5%和95%的VI作为VIs和VIv。文中用到两种植被指数：归一化植被指数NDVI、基于两波段增强型植被指数EVI2，计算公式分别为

(5)

(6)

2.4 植被调整垂直干旱指数

为进一步减小混合像元中植被因素的影响，本文采用吴春雷等[9]在PDI模型的基础上引入表征植被覆盖程度的垂直植被指数(perpendicular vegetation index,PVI)，构建PVI-PDI二维特征空间(图4)，并提出了适用于有植被覆盖情况下的植被调整垂直干旱指数VAPDI，计算公式为

(7)

式(7)中：PDI(A)、PVI(A)分别代表图4中PVI最大值A点对应的PDI、PVI指数，其中PVI的计算公式为

(8)

图4 PVI-PDI特征空间内像元散点分布示意图Fig.4 Distribution schematic diagram of pixels in PVI-PDI feature space

2.5 基于干旱指数的土壤湿度模型

基于PDI、MPDIN、MPDIE、VAPDI的计算公式，计算成像时间内试验区各像元的指数值，并对观测站点的土壤湿度实测值与相应的PDI、MPDIN、MPDIE、VAPDI指数值进行线性拟合，建立土壤湿度反演模型，其表达式为

SM=mI+n

(9)

式(9)中：SM为土壤湿度反演值；I为4种干旱指数之一；m为模型斜率；n为模型截距。

2.6 土壤湿度遥感监测技术流程

首先对多幅GF-1原始影像进行预处理，得到红光和近红外波段的光谱反射率；其次，利用反射率值提取土壤线斜率和截距，并计算研究区域的PDI、MPDIN、MPDIE、VAPDI指数值；接着将10 d原位土壤湿度实测值与对应的4种干旱指数分别构建4种土壤湿度反演模型，并进行定量评估；最后选择精度最高的反演模型对豹澥区域进行土壤湿度空间格局分析。技术流程如图5所示。

图5 土壤湿度反演技术流程图Fig.5 Flow chart of soil moisture retrieval technology

3 实验过程与结果分析

3.1 影像预处理

基于ENVI 5.3软件实现遥感卫星影像的预处理，包括绝对辐射定标、大气校正、几何校正、投影转换和图像裁剪等，最终得到光谱反射率。首先，将中国资源卫星应用中心下载到的绝对辐射定标系数带入定标公式中，并在ENVI的bandmath中进行辐射定标计算；其次，利用ENVI的FLAASH Atmospheric Correction工具完成大气校正，其中本文主要用到的大气模型为中纬度夏季，气溶胶模型为城市气溶胶；然后以三次多项式为计算模型对大气校正后的影像进行几何校正(将经过系统几何校正的2A级产品作为基准影像)；接着将影像的投影转换为经纬度投影“Geographic Lat/Lon”，基准为“WGS-84”；最后，依照豹澥试验区的经纬度进行影像裁剪。

3.2 土壤线提取

获取土壤线的方法通常有两种：常规方法和自动算法[16]。其中常规方法是在ENVI 5.3中按照分类图像选出纯净裸土像元，并在Nir-Red二维光谱空间中拟合出土壤线；自动算法通过自适应区间选择来自动获取土壤线参数，具体操作步骤如下：首先，利用归一化水指数(normalized difference water index,NDWI)阈值法和监督分类法对预处理后的GF-1影像中的水体部分进行掩膜剔除。其次，在ENVI5.3的Scatter Plot工具中，将影像的Band3作为横坐标、Band4作为纵坐标画出所有像元反射率的二维平面散点图，对图中最下方的散点以ROI(region of interest)方式保存。由于所得到的像元数量过多，大大降低了寻找所有土壤点(横坐标所对应最小纵坐标值的点)的效率，因此我们将散点中Band3反射率的最大值与最小值的差值作为组长，并分为100组，找到每组Band3中对应的最小Band4反射率值，组成初始土壤点集。接着以Band3总区间范围的0～50%，0～75%，0%～100%，25%～75%，25%～100%，50%～100%作为区间子集，计算每个子集的决定系数，决定系数最大的子集作为最佳土壤点集。最后，利用最优区间的红光和近红外波段的反射率拟合出土壤线方程。

以2020年1月29日为例，利用常规方法和自动算法进行土壤线提取并比较，两种方法的拟合情况如图6所示，并反映了豹澥试验区红光和近红外波段反射率的关系。从拟合结果中可以看出，二者得到的土壤线斜率相差不大分别为1.049 6和1.000 7，截距分别为0.064 6和0.092 2，但自动算法提取土壤线的R2较高，为0.995 5，常规方法的R2次之，为0.9809。

图6 基于两种方法拟合土壤线Fig.6 Fitting the soil lines based on two methods

为进一步比较两种提取方法的准确性，利用2019—2020年的GF-1数据提取土壤线，图7所示为各参数的变化情况。图7(a)、(b)中，通过自动算法和常规方法提取到的土壤线斜率和截距无显著差异，图7(c)中，常规方法的决定系数(R2)的均值为0.973 5，自动提取算法的决定系数(R2)的均值为0.993 6，后者的拟合效果较好，因此采用自动算法提取豹澥区域土壤线更具有优势。

图7 两种提取土壤线方法的相关参数变化情况Fig.7 Changes of relevant parameters of two methods of soil line extraction

3.3 土壤湿度模型的建立与验证

3.3.1 土壤湿度模型的建立

收集了2019—2020年共10 d站点A、站点B、站点C、站点D、站点E的GF-1遥感影像数据，并将提取到的红光、近红外波段的光谱反射率和土壤线斜率a、截距b值分别代入式(2)、式(3)、式(7)中，获取不同站点每天的PDI、MPDIN、MPDIE和VAPDI值。由于站点C处传感器出现电量耗尽的情况，因此造成后3 d的数据缺失。

图8 4种指数值随实测土壤湿度的变化趋势Fig.8 Variation trend of 4 types of index values with measured soil moisture

图8所示为4个地面观测站的原位土壤湿度与相应的干旱指数的变化趋势，可以看出，4种指数的变化趋势与实测土壤湿度的变化趋势刚好相反，例如图8(b)，实测土壤湿度值在1月29日呈现最大值为0.295 2，当天计算出的PDI、MPDIN、MPDIE、VAPDI为0.181 9、0.124 3、0.140 5、0.139 6，分别对应这段时间内的最小值；实测土壤湿度值在5月3日呈现最小值为0.117 2，当天计算出的PDI、MPDIN、MPDIE、VAPDI为0.225 6、0.223 5、0.232 2、0.238 4，分别对应这段时间内的较大值，图8(a)、图8(c)、图8(d)可以显示同样的结果。因此，可以较好地说明这4种干旱指数的意义，即指数值越大，土壤湿度越小，地表越干旱，反之亦然。

将采集到的10 d站点D原位土壤湿度值依次与对应的PDI、MPDIN、MPDIE、VAPDI指数值进行线性拟合分析，建立土壤湿度反演模型，四种模型的拟合结果如图9所示。从拟合结果可以看出，线性模型下的PDI、MPDIN、MPDIE、VAPDI指数与实测土壤湿度之间均呈现负的相关关系，且决定系数(R2)分别为0.698 5、0.713 3、0.714 3和0.755 2，这表明4种指数与土壤湿度之间存在较强的线性关系；此外，PDI的决定系数小于MPDIN、MPDIE和 VAPDI，说明MPDI、VAPDI的拟合效果较优于PDI，这是因为将研究区域内的水体像素剔除后，剩下的像元是由植被和土壤等端元组成的混合像元，基于PDI的反演模型没有剔除混合像元中的植被信息，具有一定的局限性，而修正后的MPDI和VAPDI反演模型考虑了植被因素的影响，并减少这一部分带来的偏差，因此在土壤湿度反演模型的建立上有着明显的优势。

3.3.2 反演模型的精度评价

提取成像时间内站点A、站点B、站点C、站点E处像元对应的PDI、MPDIN、MPDIE、VAPDI指数值，分别代入上述建立的4种土壤湿度反演模型中，得到4处站点10 d的反演值。利用豹澥试验区收集到的真实数据与计算到的预测数据进行验证分析，绘制出二者的二维散点图，如图10所示；并计算土壤湿度反演值与实测值的决定系数R2、平均相对误差(mean relative error, MRE)、均方根误差(root mean square error,RMSE)等评价指标，对4种反演模型进行精度评估，如表2所示。决定系数衡量反演值与实测值之间的拟合优度，平均相对误差反映反演值与实测值之间平均误差的大小，均方根误差描述反演值与实际值之间的离散程度。

图10中，1:1线代表斜率为1的直线，散点位于1:1线附近表明计算得到的土壤湿度估计值能够较真实地反映土壤湿度情况，如图10(a)，实测值在0.22～0.27区间范围内对应的土壤湿度反演值较为精确，位于1:1线上方区域土壤湿度被高估，1:1 线下方区域土壤湿度被低估，出现这种现象的原因可能是由于严重的红光饱和问题[17]。通过对比四种模型的误差拟合结果与斜率为1的直线偏差情况可以发现，VAPDI偏差更小，优于其他模型，图10(b)～图10(d)可以得到同样的结果。

图9 PDI、MPDIN、MPDIE、VAPDI与实测土壤湿度线性拟合结果Fig.9 Linear fitting results of PDI, MPDIN, MPDIE,VAPDI and measured soil moisture

图10 不同模型下土壤湿度反演值与实测值的二维散点图Fig.10 Two-dimensional scatter plots of soil moisture inversion and measured values under different models

表2 基于GF-1遥感影像数据土壤湿度反演模型的评价指标值Table 2 Evaluation index value of soil moisture retrieval model based on GF-1 remote sensing image data

通过比较表2，并结合研究区域实际植被覆盖度的分布情况可知，在无植被覆盖的站点A处，基于PDI、MPDIN、MPDIE、VAPDI反演模型得到的反演值与实测值的R2均在0.735 0左右、MRE均在4.50%左右、RMSE均在1.10%左右，可以看出这四种模型的反演精度相差不大；在有植被覆盖处，如站点C，基于PDI、MPDIN、MPDIE、VAPDI反演模型的验证结果中R2分别为0.583 9、0.658 0、0.709 9、0.727 7，MRE分别为5.89%、5.86%、5.21%、4.24%，RMSE分别为1.82%、1.70%、1.69%、1.48%，可见后三者的反演精度均高于PDI，与实测值更为接近，站点B、站点E可以得到类似的结果。

这主要是因为在光谱反射率中，红光和近红外波段的细微变化对其结果会产生较大的影响，其中植被端元对红光谱有较强的吸收特性、对近红外光谱有较强的反射特性，但是PDI没有考虑到这一因素带来的影响，因此只适用于裸土和低植被区域，无法准确地将中、高植被区域的土壤湿度真实情况反映出来。而MPDIN指数中添加了可以反映植被冠层的植被指数NDVI，MPDIE指数中添加了不受大气干扰和土壤背景噪音影响的植被指数EVI2，VAPDI中添加了不易出现植被覆盖饱和的垂直植被指数PVI，这些模型克服了植被信息对土壤湿度的影响，这与前文所述的结论是一致的。

结合图10和表2进一步比较基于MPDI和VAPDI指数的反演模型，发现基于VAPDI指数在高植被覆盖区域估算到的土壤湿度要比MPDI更为准确，这是因为MPDI中fv和植被像元反射率的定量评估较为复杂，通常使用经验阈值，这在极端条件下可能会导致不确定性和差异，因此在高密度植被区域采用VAPDI模型效果最好。

3.4 土壤湿度空间格局分析

为了进一步验证基于VAPDI指数反演土壤湿度的准确性，选用基于VAPDI指数建立的土壤湿度反演模型对研究区域2019年12月27日、2020年1月29日、3月14日、4月16日4个时间点进行空间格局分析，分布图如图11所示。

图11中4个时间点的土壤湿度数据主要集中在10%～32%，这与我们实测到的土壤湿度范围相匹配。结合相关的气象数据可知，图11(a)、图11(b)为气候较干旱的12月—次年1月，1处椭圆区域和2处区域的颜色占比最大的为红色，代表土壤湿润的蓝色像元较少，整体像元颜色分布差异不大；而图11(c)、图11(d)为气候较湿润的3—4月，不仅存在土壤较干旱的黄色和红色像元，同时也存在较多代表土壤湿润的蓝色像元，如2处和3处区域，整体像元颜色分布差异较大也较为离散。这说明在有植被遮挡的情况下，由10 d建立的VAPDI模型可以清晰地将豹澥区域的土壤湿度空间分布情况反映出来，且有较好的区分地表土壤湿度的能力。

图11 豹澥区域4个时间点的土壤湿度空间分布Fig.11 Spatial distribution of soil moisture at four time points in Baoxie area

4 结论

采用高时空分辨率的GF-1 WFV遥感影像作为数据源，以豹澥野外观测场为研究区域，基于2019—2020年6个月的PDI、MPDIN、MPDIE和VAPDI指数值与观测站点实测数据之间的线性关系，分别建立4种土壤湿度监测模型，并在豹澥试验区进行验证。结果表明：

(1)利用常规方法和自动提取算法对多天的土壤线参数进行提取并比较，发现通过两种方法得到的斜率和截距大致相同，而二者10 d决定系数R2的均值分别为0.973 5、0.993 6，因此自动算法效果更好，可以降低土壤线对模型带来的偏差。

(2)站点D计算到的4种干旱指数与相应实测值间的决定系数R2分别为0.698 5、0.713 3、0.714 3、0.755 2；结合站点A、站点B、站点C、站点E处土壤湿度反演模型的评价指标(R2、MRE、RMSE)可以发现，在有植被覆盖的情况下基于VAPDI指数的模型相对于其他三种模型比具有较高的反演精度。

由于光学遥感数据多受云量等天气因素的影响，使得土壤湿度的提取变得困难和复杂，后续将引入对地物穿透能力较强的微波遥感技术，减少大气的影响，进一步提高对土壤湿度反演的精度。