基于Sentinel-2A影像的OPTRAM模型及其改进模型的土壤水分估算研究

朱燕香,潘剑君,白浩然,康翔

(南京农业大学资源与环境科学学院,江苏 南京 210095)

土壤水分在农业生产、气象变化和水文过程中具有重要作用,是干旱监测和管理的关键[1]。土壤水分传统研究耗时耗力,且研究区域受限[2],而遥感技术为大范围、多时空监测土壤水分提供了有力手段,同时还具有快速、及时和便捷等优点[3]。遥感技术发展至今,各国为监测地表土壤水分已发射多种高时空分辨率的卫星,如欧洲土壤湿度与海洋盐度(SMOS)卫星[4]、高级散射仪(ASCAT)[5]以及主被动微波卫星SMAP[6]等。

近年来,国内外学者对土壤水分监测进行了大量的研究[7-10]。条件植被温度指数(VTCI)是利用条件植被指数、条件温度指数以及距平植被指数对研究区年度间相对干旱程度进行监测的指数[11]。詹志明等[12]利用近红外波段和红光波段反射率建立光谱特征空间,根据土壤水分在该特征空间的分布规律,提出遥感干旱监测新指数——垂直干旱指数(PDI),改进的垂直干旱指数(MPDI)[13]是在垂直干旱指数的基础上加入植被覆盖度后提出的干旱监测指标。吴春雷等[14]根据垂直干旱指数和垂直植被指数的空间分布规律提出了植被调整垂直干旱指数(VAPDI)。

传统的梯形模型(thermal-optical trapezoid model,TOTRAM)是一种将热遥感数据和光学遥感数据相结合反演土壤水分的模型。如:由Sandholt基于地表温度(Ts)与归一化植被指数(NDVI)之间关系的空间数据关系,提出的温度植被干旱指数(TVDI)[15]。它存在2个明显的问题:模型所用遥感数据必须同时包含光学波段和热红外波段,故不适用于不提供热波段数据的卫星(例如:Sentinel-2);地表温度不仅受土壤水分的影响,还受大气条件如风速、空气温度和空气湿度等影响[16]。因此,传统的梯形模型需要多种对应的大气参数进行校准。为解决这2个问题,Sadeghi等[17]提出基于物理的OPTICAL TRApezoid模型(OPTRAM)反演地表土壤水分。OPTRAM模型是由NDVI与短波红外变换反射率(STR)组成的,它不需要热波段数据,可用于仅提供光学波段的卫星;并且,短波红外反射率不会随大气条件发生显著变化。因此,OPTRAM模型可以免于大气参数校准。

Sentinel-2是具有时空分辨率高、重返周期短等优势的遥感卫星,但由于该卫星发射时间较短,目前光学遥感的土壤水分反演研究主要集中在基于Landsat 或MODIS等卫星的大中型研究区,应用Sentinel-2影像反演小区域土壤水分的研究较少。本文以Sentinel-2A影像为数据源,应用OPTICAL TRApezoid模型(OPTRAM)对南京市六合区北部地区进行土壤水分反演适用性研究,并利用NDVI-DFI像元三分模型对该模型进行改进,使土壤水分反演精度提高,并分析改进模型反演取得更高精度的原因,从而推动后续的土壤水分反演工作。

1 研究区域与数据

1.1 研究区概况

研究区位于江苏省南京市六合区(118°33′～119°00′E,32°21′～32°37′N),约有750 km2,地貌大部分属宁镇扬山区,北部为丘陵岗地,南部为平原岗地,平均海拔20 m。该地区位于北亚热带湿润的季风气候区,雨热同期,日照充足,四季分明,年平均气温为16 ℃,年平均降雨量为1 157 mm。研究区土地利用类型多样,以农用地为主,林地所占比例较大,全区林地覆盖率约为6%。主要农业种植物是油菜、水稻和小麦,数据采集时为小麦成熟期。

1.2 实测数据

2019年5月22—23日对研究区进行野外数据采集。获得的数据包括TDR(time-domain reflectometry)水分仪测得的表层10 cm土壤水分数据、GPS点位数据以及采样点的土地利用状况数据等。研究区地理位置及采样点分布如图1所示。每个采样点测定5个数据,取平均值作为采样点的实测值。本次采集数据共39个,用于模型的拟合验证。

1.3 遥感数据

Sentinel-2A是2015年6月23日成功发射的第1颗“哥白尼计划”的光学遥感卫星,有13个光谱波段,重返周期为10 d,空间分辨率最高可至10 m,在非商业卫星中具有极大的优势,常被用于全球陆地观测以及自然灾害的监测与管理等。

本研究获得2019年5月22日2景Sentinel-2A L1C影像,均来自于欧洲空间局(https://scihub.conpernicus.ed/dhus/#/home)。Sentinel-2A L1C影像是经过正射校正和几何精校正后的大气表观反射率产品,但是并没有经过辐射定标和大气校正。利用欧洲空间局提供的Sen2cor插件,生成能直接使用的Sentinel-2A L2A数据。重采样分辨率至10 m,在ENVI中将所需的6个光谱波段(Band 2、Band 3、Band 4、Band 8、Band 11和Band 12)合成,再进行影像拼接和裁剪,并运用归一化水指数(normalized difference water index,NDWI)和NDVI指数阈值法进行水体掩膜处理,监督分类法进行不透水层掩膜处理。

2 研究方法

2.1 OPTRAM模型

有研究发现短波红外转换反射率(shortwave infrared transformed reflectance,STR)与裸土和植被覆盖区根系土壤水分之间有线性物理关系[16]。Sadeghi等[17]以此为基础提出一种新的光学梯形模型(OPTICAL TRApezoid Model,OPTRAM)。该模型根据STR与NDVI的特征空间分布进行土壤水分估算。

(1)

式(1)中:STR代表短波红外转换反射率;STRd和STRw分别代表土壤的干、湿状态的STR。

(2)

式(2)中:RSWIR表示短波红外的反射率。本文中短波红外波段为Sentinel-2A中的中心波长为2 190 nm的Band 12。

STRd=id+sdNDVI

(3)

STRw=iw+swNDVI

(4)

式(3)、式(4)中:id、sd和iw、sw分别表示干边和湿边拟合方程的截距和斜率,用于求解梯形模型中有固定间隔的NDVI的短波红外转换反射率最大值STRmax和短波红外转换反射率最小值STRmin。目前,已有不少研究将遥感植被指数量化进而分析植物生长状态并将其与根系土壤水分联系起来。土壤水分状况会影响植被的水分状况,从而改变植被的光谱特征。本文假定土壤水分与植被含水量存在线性关系。

由上述公式可以看出该模型与传统的梯形模型类似,不同之处在于以STR代替地表温度(land surface temperature,LST)。OPTRAM模型的STR-NDVI散点图也存在干边和湿边之分。干边代表地表干旱情况下的STR,湿边代表地表水分充足情况下的STR。

2.2 改进的OPTRAM模型

2.2.1 NDVI-DFI像元三分模型非光合植被(non-photosynthetic vegetation,NPV)是指不能将太阳能转化为用于其自身生长和发育的化学能的植物[18],包括衰老或死亡的植被。常见非光合植被有凋落物、作物茬、枯叶等。在传统梯形模型和OPTRAM模型中常用NDVI作为植被量的衡量指标,但是模型中真正需要的信息为光合植被的覆盖量,所以单纯用植被指数去计算会出现一定的误差。国内外学者发现将干枯燃料指数(dead fuel index,DFI)用于估算非光合植被覆盖度时具有很大的潜力,将其与NDVI相结合建立NDVI-DFI模型,通过该模型的二维空间分布可以将混合像元分解为土壤(BS)、光合植被(PV)和非光合植被(NPV)3个端元。

(5)

(6)

式中:RSWIR1、RSWIR2、RRED、RNIR分别为Sentinel-2A影像的Band 12、Band 11、Band 4和Band 8。为避免水体和建筑用地等影响产生异常值,在计算指数前对影像进行指数阈值法和监督分类法掩膜处理。

图2为NDVI-DFI像元三分模型示意图,特征空间近似表现为三角形,其中NPV的NDVI高、DFI低,位于三角形的左上部;BS的NDVI和DFI均很低,位于三角形的左下部;PV的NDVI高,DFI很低,位于右侧中部;混合像元分布在三角形内部[19]。

2.2.2 OPTRAM模型的改进在ENVI中对预处理过的Sentinel-2A影像进行最小噪声分离(MNF)变换降维,选择能代表大部分信息的前5个分量进行投影。然后设定参数,其中迭代次数为2 000,阈值系数为2.5,得到像元纯净指数(PPI)。最后取PPI大于3且又靠近NDVI-DFI模型特征空间定点的像元为纯净像元,分别计算3个顶点纯净像元的平均指数作为相应端元的特征值。

用PPI指数法确定BS、PV、NPV端元值,利用ENVI 5.2软件中的Triangular Linear Spectral Unmixing扩展工具对预处理的Sentinel-2影像进行分解,得到光合植被覆盖度(Fpv)、非光合植被覆盖度(Fnpv)和土壤覆盖度(Fbs)的RGB合成图。用Fpv代替OPTRAM模型中的NDVI,用以改进OPTRAM模型。

3 结果与分析

3.1 NDVI-DFI模型反演光合植被覆盖度

3.1.1 NDVI-DFI模型可分性从图3可知:当NDVI小于0.25时,DFI随NDVI的增大而增大;当NDVI大于0.25时,DFI随NDVI的减小而减小。同时能够看出该特征空间分布表现为三角形。BS、PV、NPV在NDVI和DFI的分布具有明显的差异,可以将BS、PV和NPV 3个端元区分开来,即NDVI-DFI模型用于本研究区混合像元分解可行。因此,本文基于Sentinel-2A影像所构建的NDVI-DFI特征空间满足该模型线性分解的基本假设。

3.1.2 混合像元三分结果研究区的植被覆盖度空间差异较大,植被覆盖达到100%的区域主要位于西北部以及中部地区的山地和林地。像元二分模型假设混合像元为土壤和植被2种端元组成,本文分别取NDVI累积频率置信度5%和95%作为NDVIsoil和NDVIvegetable的近似值对混合像元进行分解。由ENVI对像元二分植被覆盖度(Fv)的6 261 761个像元统计结果(图4-A)可知,零植被覆盖度为研究区的4.89%,全植被(100%)覆盖度为5.15%。在0

像元三分模型得到的光合植被覆盖度(图4-B)与像元二分模型得到的植被覆盖度数值大致一致,但是整体上像元三分模型得到的光合植被覆盖度要大于像元二分植被覆盖度。山湖水库正下方的16号采样点位于恋山坝上草原,该样点无裸露的土壤,像元三分模型得出光合植被覆盖度(Fpv)为50.78%,非光合植被覆盖度(Fnpv)为49.22%,像元二分结果的植被覆盖度(Fv)为62.66%,土壤覆盖度(Fbs)为37.34%;山湖水库上方的28号采样点为一片荒地,土壤、枯草和绿草混合交错,像元三分结果显示:Fbs、Fnpv和Fpv分别为22.22%、49.79% 和27.99%,像元二分结果Fv为35.96%,Fbs为64.04%。Fpv低于Fv的原因可能是像元二分模型直接将非光合植被端元当作土壤端元进行处理,无法判别非光合植被和土壤,还可能因为像元二分模型存在将部分土壤与光合植被混肴的问题。

3.2 OPTRAM及其改进模型特征空间的构建结果

取NDVI=0.01为固定间隔,切割STR与NDVI组成的二维空间散点图,分别找到间隔内STR的最大值(STRmax)和最小值(STRmin),最后线性拟合STR-NDVI值,得到干边和湿边。OPTRAM特征空间分布图以NDVI为横坐标,STR为纵坐标。

由图5和表1可知:与传统的梯形模型相似,OPTRAM和改进的OPTRAM的特征空间分布均表现为梯形,干、湿边拟合效果良好,其中OPTRAM的STR-NDVI空间的干边拟合决定系数(coefficient of determination,R2)为0.56,湿边R2为0.34,改进的OPTRAM的STR-Fpv特征空间的干、湿边拟合度都优于STR-NDVI特征空间,干边R2为0.60,湿边R2为0.42,2个模型的干边均表现为负相关,湿边均表现为正相关。

表1 STR-NDVI和STR-Fpv的干、湿边拟合Table 1 The fitting dry and wet edges of STR-NDVI and STR-Fpv

NDVI能够增加植被在遥感影像中的反应特征,同时减少土壤反射信息的影响。但是NDVI在反映植被覆盖时仍有一定的局限性,当植被覆盖度较低,NDVI对植被量的指示作用下降[20]。土壤和非光合植被的光谱反射曲线相似,仅用NDVI是很难区分2个端元,所以STR-NDVI特征空间中当NDVI较小时,与STR的特征空间里的散点分布比较离散。由于NDVI的饱和性问题,当NDVI达到0.8左右,就无法再表达更加茂盛的植被覆盖度情况。由图5可以看到:STR-NDVI特征空间中当NDVI大于0.8时,与STR的拟合曲线更趋于曲线,而非直线。STR-Fpv特征空间中用光合植被覆盖度Fpv表示植被覆盖状况,在植被量过低或过高时,解决了NDVI对植被覆盖度反应灵敏度下降的问题。2个模型的特征空间的干、湿边拟合精度结果显示,改进的OPTRAM模型更具优势。

3.3 与实测含水量的相关性验证

利用野外获取的GPS经纬度定位信息,提取对应的各采样点的TVDI、OPTRAM、OPTRAM(Fpv)值,与TDR水分仪实测的39个地表10cm土壤水分含量进行最小二乘回归分析。从图6可见:土壤水分与遥感反演指数均表现为良好的负相关线性关系,TVDI、OPTRAM和OPTRAM(Fpv)模型计算的指数值越高,土壤含水量越低,土壤干旱越严重。因此,本文所用指数均能用于土壤水分反演研究。OPTRAM与实测土壤含水量的拟合决定系数(R2=0.38)与较为成熟的TVDI方法的拟合结果(R2=0.39)相近。这一结果表明OPTRAM模型在南京六合区反演土壤水分的可适用性,其反演精度与TVDI反演精度近似。

改进的OPTRAM(Fpv)模型的模拟值与实测数据的拟合效果达显著水平(R2=0.55,RMSE=0.03,P<0.01),反演得到的指数能反映土壤水分情况,相较于OPTRAM模型的拟合效果(R2=0.38,RMSE=0.03,P<0.01)大大提高。去除非光合作用植被(NPV)的影响,改进的OPTRAM(Fpv)模型能更加准确地反映研究区的植被覆盖状况。改进的OPTRAM(Fpv)模型的反演效果要比OPTRAM模型的效果好,有更高的拟合精度和更好的模型稳定性。

4 结论与讨论

本文基于Sentinel-2A影像通过构建STR-NDVI特征空间建立OPTRAM模型对南京市六合区进行土壤水分反演研究。通过NDVI-DFI像元三分模型将研究区混合像元分解为土壤、光合植被和非光合植被3种端元,将得到的光合植被覆盖度(Fpv)引入OPTRAM模型,构建STR-Fpv特征空间,建立改进的OPTRAM模型。

由于研究区为多种地表覆盖混合地区,并非在常应用NDVI-DFI模型的草原等地区,生成的特征空间没有表现为一个理想的三角形。研究区以农用地为主,林地所占比例也较大。5月为植被生长茂盛时期,大量的光合植被存在,土壤和非光合植被覆盖,导致DFI值较低。同时,研究区的纯净光合植被像元要远远多于纯净的土壤像元和非光合植被像元,使空间分布倾向于光合植被。

蜡熟期小麦开始衰亡,到完熟期体内水分含量降低至10%左右,叶绿素含量几乎为0,干物质停止累积,麦秆逐渐失去弹性,光合作用减弱,吸收率下降,红光反射率上升[21]。植被指数是基于植被叶绿素在0.69μm处的强吸收,像元三分模型根据光谱反射率计算得到的NDVI和DFI将完熟期小麦归类为非光合植被。完熟期冬小麦反射率较基底枯落物(长期的NPV)反射率偏高,计算得到的DFI偏低。由于裸土像元与非光合植被光谱相似,像元二分模型直接将完熟期的小麦作土壤处理,这也是像元二分模型将土壤与非光合植被混肴的根本原因。像元三分模型得到的光合植被覆盖度(Fpv)与像元二分模型得到的植被覆盖度Fv分布总体相似,但是Fpv值要略低于Fv。

OPTRAM是基于STR-NDVI特征空间构建的模型,其核心是干、湿边拟合方程的确定[22]。OPTRAM模型的STR-NDVI特征空间分布与改进的OPTRAM(Fpv)模型的STR-Fpv特征空间分布为明显的梯形,干、湿边拟合效果显著,均表现干边为负相关,湿边为正相关。STR-Fpv特征空间的干、湿边拟合方程的决定系数都大于STR-NDVI特征空间的决定系数。

传统的梯形模型(如TVDI)是较为成熟的土壤水分反演研究方法。温度植被干旱指数(TVDI)利用植被指数与地表温度的强相关性和对土壤水分信息的互补,通过构建LST/NDVI特征空间,确定植被指数在水分充足和严重缺水条件时的地表温度之间的关系,达到反演土壤水分的目的。近年来,国内外不少学者用TVDI来监测土壤干旱情况[23-25]。TVDI、OPTRAM和改进的OPTRAM(Fpv)与实测土壤水分数据拟合结果显示为显著负相关(P<0.01),OPTRAM与TVDI的反演效果近似,在南京地区OPTRAM模型具有可适用性。相对于OPTRAM模型,改进的OPTRAM(Fpv)模型与实测土壤水分的拟合精度大大提升。