基于Sentinel多源遥感数据的河北省景县农田土壤水分协同反演

李伯祥,陈晓勇,2,3①

(1.东华理工大学测绘工程学院,江西 南昌 330013;2.流域生态与地理环境监测国家测绘地理信息局重点实验室,江西 南昌 330013;3.江西省数字国土重点实验室,江西 南昌 330013)

土壤水分(soil moisture，SM)是全球水循环、碳平衡和能量交换的关键变量，在水文管理、气候预测、农业生产等方面发挥着重要作用[1]。土壤水分具有较高的时空异质性，传统基于监测站点进行土壤水分观测，或者烘干称重法能够获取高精度的土壤水分信息，但是不能胜任大范围、时序性的监测任务[2]。随着遥感技术的发展和基础理论的完善，全天候、全天时、穿透力强的微波遥感在土壤水分反演方面得到广泛的应用[3]。相比于光学遥感只能获取表层土壤墒情，微波遥感能够有效地穿透植被覆盖层，并灵敏地探测地表土壤水分动态变化[4]。微波传感器接收的地表后向散射系数与土壤介电特性密切相关[5]，而土壤介电特性主要由土壤含水量决定[6]，因此基于微波遥感估算土壤水分具备良好的物理基础[5-7]。

因目标区域不同，微波遥感反演土壤水分可分为裸露地表区和植被覆盖区[8]。针对裸露地表区，前人提出了Oh模型[9]、Dubois模型[10]和Shi模型[11]等经验-半经验模型，以及几何光学模型[12-13]、物理光学模型[13-14]、小扰动模型[13,15]、积分方程(IEM)模型[16]和高级积分方程(AIEM)模型[17]等理论模型。针对植被覆盖区，则提出了基于植被散射理论的密歇根微波冠层散射模型(Michigan microwave canopy scattering，MIMICS)[18]和Karam模型[19]，与基于理论模型发展而来的Roo模型[20]和水云模型[21]，这些模型通过估算地面总后向散射中的植被后向散射贡献量，从而消除植被的影响，有效获取下垫面的土壤水分信息。随着计算机技术和机器学习的发展，人工神经网络[22]、遗传BP神经网络(GA-BP)模型[23]、支持向量机(support vector machine，SVM)模型[24]、深度学习[25]等机器学习方法也被广泛地应用到土壤水分反演。

微波后向散射系数不仅受到雷达传感器频率、极化方式、入射角的影响，还受到地表植被生长状况和地表粗糙度等观测区域环境因素的影响，其中植被含水量对雷达后向散射系数影响最大[26]。研究表明光学遥感数据可以通过计算植被的各种生理参数更好地补充植被的影响[27]。因此，在通过改进原有的土壤水分反演模型的基础上，充分利用微波和光学遥感的各自优势，减少对地面实测数据的依赖性，实现主被动遥感协同反演成为近年来的研究热点[2]。曾旭婧等[27]以Sentinel-1 SAR数据为基础，并辅以Landsat-8光学遥感数据，对黑河-北安高速沿线地区不同植被覆盖度程度下地表土壤水分进行反演研究。MENG等[28]利用Radarsat-2和Landsat-8遥感数据，结合改进Dubois模型和最优冠层含水量反演模型，反演全生育周期的玉米植被覆盖区地表土壤水分。汪倩倩等[29]利用GF-3和Landsat-8遥感数据，分别计算土壤后向散射系数和干旱指数，并分析两者和土壤水分的关系，通过统计回归分析建立地面资料稀缺农田区域的土壤水分反演模型。

目前在光学遥感数据选择上使用较多的是Landsat-8遥感数据，但Landsat-8重返周期为16 d，空间分辨率为30 m，难以满足小区域农田土壤水分反演高时空分辨率的要求。Sentinel-1和Sentinel-2同属于欧空局Sentinel系列卫星(后面叙述中简称S1和S2)，S2空间分辨率为10 m，S2A/B组成的双星对地观测系统过境重返周期仅为5 d。因此，采用S2搭配S1进行主被动遥感协同反演土壤水分研究在空间、时间和数据配准方面更具有优势。

该研究针对植被覆盖区农田土壤水分反演精度不高、空间分辨率低、植被参数难以准确获取等问题，基于S1主动微波和S2光学遥感数据，采用改进水云模型和Oh模型的组合方法，对植被覆盖地表土壤水分进行定量反演研究。

1 研究区和数据收集

1.1 研究区概况

选取河北省衡水市景县作为研究区(图1)，景县总面积为1 183 km2，平均海拔30 m以下，地形以平原为主。景县被列为全国商品粮基地县，主要地类类型为耕地，主要种植作物为玉米、小麦和棉花。农业灌溉用水主要依靠开采地下水和自然降水补给。该地区属半湿润温带大陆性季风气候区，年平均气温为12.5 ℃，年平均降水量为554 mm。景县实地考察和采样区范围为37°30′～37°40′ N，116°10′～116°20′ E，地表植被覆盖类型为玉米。

图1 研究区地理位置及采样点分布

1.2 遥感数据处理

研究使用的S1 SAR和S2光学遥感数据可从欧洲航天局( European Space Agency，ESA)免费下载。S1搭载了中心频率为5.405 GHz的C波段合成孔径雷达，具有条带、干涉宽幅、超幅宽和波模式4种对地观测成像模式[30]。主动微波数据选择S1A干涉宽幅模式下的地距(ground range detected，GRD)产品，空间分辨率为5 m×20 m，包括垂直发射水平接收(VH)和垂直发射垂直接收(VV)2种极化模式。GRD产品经过了多视处理，并采用 WGS-84椭球基准投影至地距的聚焦数据，数据获取日期为2018年9月12日，使用ESA提供的SNAP(sentinel application platform)软件进行预处理，包括滤波、辐射定标、几何校正和地理编码，预处理完成后即可获得研究区的后向散射系数(图2)。地表植被覆盖情况使用ESA提供的2018年9月7日S2A的L1C级数据进行分析，并计算各种植被指数(vegetation indices，VI)。L1C数据经过了几何精校正，需要借助ESA提供的Sen2 Cor工具实现辐射定标和大气校正处理，生成L2A级数据。

图2 研究区S1 VH/VV极化的后向散射系数

1.3 实测数据处理

野外观测日期是2018年9月8日。实测采样点数量为21个，采样点坐标信息使用手持式GARMIN-eTrex301定位仪(GPS+GLONASS双系统定位，定位精度1～3 m)记录。使用卷尺测量玉米植株高度，并统计1 m2农田区域内植株的种植密度。每个采样点收集1株玉米植株样本，带回实验室后称量植株鲜重(WF)，然后用电热恒温鼓风干燥箱对植株进行烘干处理，获得植株干重(WD)。

植被层含水量(vegetation water content，CVW)计算公式为

CVW=(WF-WD)·ρ。

(1)

式(1)中，WF为玉米植株鲜重，kg；WD为玉米植株干重，kg；ρ为玉米植株种植密度，株·m-2。

土壤样品采样深度为0～5 cm，样本用铝盒封装带回实验室，将其放进干燥箱中加热48 h直至恒重。采用烘干称重法测定土壤质量含水量，结合采样区的土壤容重将土壤质量含水量(θm)转换为土壤体积含水量(θv)。

(2)

θv=θm×sc。

(3)

式(2)～(3)中，θm为土壤质量含水量，%；m1为原始采样湿土质量，g；m2为烘干后的干土质量，g;sc为土壤容重，由于缺乏相应的土壤容重测定仪器，参照文献[28]，取值1.33 g·cm-3。

经计算，实测样点土壤体积含水量的范围为0.04～0.29 cm3·cm-3，平均值为0.13 cm3·cm-3，标准偏差为0.06 cm3·cm-3。

2 研究方法

2.1 水云模型

植被层对微波信号具有散射和吸收作用，为了获取准确下垫面裸土后向散射系数，首先需要利用水云模型去除植被层的影响。ATTEMA等[21]基于辐射传输方程提出了经典的水云模型方程，该模型假设：(1)植被层由大量形状、大小相同且均匀分布的散射微粒组成；(2)忽略植被与土壤表面的多次散射，仅考虑来自裸土的表面散射和来自植被的体散射；(3)模型中的变量仅为微波入射角、植被含水量和土壤湿度。水云模型的基本表达式为

(4)

(5)

T2=exp(-2τ·secθ)，

(6)

τ=B·CVW。

(7)

2.2 引入植被覆盖度和植株高度改进水云模型

植被覆盖地表的后向散射往往比较复杂，且对于混合像元区域，裸土后向散射是地表总后向散射的重要组成部分[31]。由于真实地表的复杂性，传统水云模型难以满足实际需要。因此，首先需要引入植被覆盖度(vegetation fractional coverage，CVF)对原始水云模型加以改进：

(8)

CVF=(INDV-INDV,min)/(INDV,max-INDV,min)。

(9)

式(9)中，INDV为归一化差值植被指数(normalized difference vegetation index)。根据研究区INDV的累积概率分布情况，取累计分布概率5%和95%对应的NDVI值作为INDV,min和INDV,max，分别取值为0.239 7和0.848 2，对应裸地和完全植被覆盖情况。

真实作物种植区地表情况较为复杂，除了需要考虑植被层含水量，还需重视植被覆盖程度和植被茎秆散射的影响。原始水云模型将整个植被层假设为水平均匀的云层，并没有考虑复杂的植被层多次散射作用。在一些如玉米、小麦等高大农作物覆盖区域和特定波长条件下，植被和土壤间的双次散射在总后向散射中占有一定比例，忽略其影响会对地表土壤水分的反演产生较大误差[30]。

(10)

T2=exp(-2τ·secθ)。

(11)

式(10)～(11)中，σpq1为每单位体积植被冠层的雷达后向散射截面，m2·m-3；kepq为冠层的消光系数。τ由消光系数kepq和植被高度(h)决定：

τ=kepq·h。

(12)

由式(7)和式(12)可推出：

(13)

将式(10)和式(13)代进式(8)，便可得到引入CVF和h改进的水云模型表达式：

(14)

式(14)中，σpq1为常量，以系数A代替。

BINDLISH等[33]提出了水云模型中对应植被、牧场、冬小麦和草地这4种应用场景时的A、B取值，但是A、B属于经验常数，需要根据具体实验环境进行确认。利用 Levenberg-Marquardt非线性最小二乘法和野外采集的实测数据分别拟合VH和VV极化条件下的A、B数值，如表1所示。

表1A、B参数拟合值

Table 1 The fitting values ofAandB

极化模式ABVH-0.022 00.080 8VV-0.020 80.076 9

将A和B取值代入式(14)，得到VH和VV极化的裸土后向散射系数的表达式：

(15)

(16)

2.3 Oh模型解算土壤墒情

改进水云模型从地表总后向散射中去除了植被层的后向散射贡献量，便获得地表裸土直接后向散射部分。根据裸土后向散射系数反演土壤墒情，需要考虑地表粗糙度的影响。Oh模型是基于同极化比(p)、交叉极化比(q)与地表粗糙度和土壤介电常数(ε)之间的关系建立的经验模型，地表粗糙度参数考虑了均方根高度(s)和相关长度(l)。Oh模型是OH等[9]于1992年提出的，经过后续的改进和完善[34-35]，具体表达式如下：

(17)

(18)

(19)

(20)

由于该次实验没有对地表粗糙度进行实地测量，根据研究区野外观测实际情况，均方根高度(s)和相关长度(l)参考文献[28]，合理设置Oh模型VH、VV极化下各输入参数范围和步长(表2)，建立裸土后向散射系数数据库。

表2 后向散射系数数据库输入参数值

Table 2 Input parameters of backscattering simulation database

输入参数入射角(θ)/(°)土壤体积含水量(mv)/(cm3·cm-3)均方根高度(s)/cm相关长度(l)/cm范围30～500～0.300.4～2.21.0～18.0步长50.010.33.0

综上，植被覆盖区农田土壤水分协同反演算法的具体技术流程如图3所示。

图3 技术路线图

3 结果与分析

3.1 植被参数反演

改进水云模型涉及植被层含水量(CVW)、植被覆盖度(CVF)和植被高度(h)，这3个植被影响因子需要建立反演模型。CVW用植被指数法进行反演，通过分析S2A计算的6种光学植被指数和实测植被含水量的相关性，选择最优光学植被指数构建植被层含水量反演模型。植被指数在PIE(pixel information expert)遥感图像处理软件中进行计算，计算的植被指数包括差值植被指数(difference vegetation index，IDV)[36]、改进红边比值植被指数(modified red edge simple ratio index，ImSR705)[37]、归一化水体指数(normalized difference water index，INDW,1610/INDW,2190)[38]、水分胁迫指数(moisture stress index，IMS)[39]、归一化多波段干旱指数(normalized multi-band drought index，INMD)[40]。

各植被指数与实测植被层含水量的相关性分析如表3所示。由S2A的近红外波段ρ842和短波红外波段ρ1610、ρ2190计算得到的归一化水体指数INDW,1610、INDW,2190和水分胁迫指数IMS，这3种植被指数与实测植被含水量具有良好的相关性，R2均大于0.60。其中INDW,1610的拟合效果最好，R2为0.834 3，均方根误差(root mean square error，RMSE)为0.737 7 kg·m-2。因此，采用INDW,1610构建植被层含水量反演模型。植被层含水量反演结果图4(a)所示。

表3 植被指数和实测植被层含水量拟合关系

Table 3 Correlations of the vegetaion indices and thein-suitmeasured vegetation water content data

植被指数公式1)植被层含水量(CVW)反演模型1)R2均方根误差(RMSE)/(kg·m-2)差值植被指数(IDV)IDV=ρ842-ρ665CVW=0.002 6×IDV+0.344 60.306 91.508 8改进红边比值植被指数(ImSR705)ImSR705=ρ740-ρ443ρ705+ρ443CVW=3.204 3×ImSR705+1.925 50.166 81.653 8水分胁迫指数(IMS)IMS=ρ1610ρ842CVW=-32.218 5×IMS+24.343 70.617 31.120 6归一化多波段干旱指数(INMD)INMD=ρ842-(ρ1610-ρ2190)ρ842+(ρ1610-ρ2190)CVW=-57.2193×INMD+43.13720.460 41.330 8归一化水指数(INDW)INDW,1610=ρ842-ρ1610ρ842+ρ1610CVW=33.0177×INDW,1610-2.074 30.834 30.737 7INDW,2190=ρ842-ρ2190ρ842+ρ2190CVW=29.155 8×INDW,2190-8.344 60.619 11.118 3

1)变量右下标数字为Sentinel-2对应波段。

S2遥感数据首先计算归一化差值植被指数，然后根据植被覆盖度公式〔式(9)〕便可得到如图4(b)所示的研究区植被覆盖度反演结果。研究发现，实测玉米植株高度和雷达极化比(VH/VV)具有良好的相关性(R2=0.460 9)〔图4(c)〕，可根据拟合关系式建立研究区玉米植株高度反演模型，反演结果如图4(d)所示。为了有效得到植被覆盖区域，需要利用归一化差值植被指数进行面积提取。结合研究区的实际情况，当INDV阈值设置为大于0.3时能有效剔除裸露土地以及城镇建设用地，获得更为精准的植被种植区域范围。

3.2 去除植被影响的土壤水分反演

从图中可以看到，相同采样点的VH极化后向散射系数低于VV极化，裸土后向散射系数普遍比地表总散射系数低。去除植被影响后的VH极化后向散射系数变化范围为-8.06～-0.94 dB，平均变化-4.41 dB；VV极化后向散射系数变化范围为-3.02～0.95 dB，平均变化-1.49 dB。校正前后VH极化后向散射的变化幅度大于VV极化，说明VH极化信号的传输过程更容易受到植被层的影响。

土壤体积含水量的预测结果和野外实测数据如图6所示。VH极化的决定系数R2为0.399 3， RMSE为0.055 6 cm3·cm-3，平均绝对误差(mean absolute error，MAE)为0.046 7 cm3·cm-3；VV极化的预测结果优于VH极化，VV极化条件下的决定系数R2为0.653 0，RMSE为0.040 1 cm3·cm-3，MAE为0.032 7 cm3·cm-3。VH极化条件下土壤水分出现预测结果偏高现象，散点离散于1∶1线附近；VV极化条件的预测结果和实测值较为接近，散点大部分位于1∶1线附近。

图4 植被含水量反演分布、植被覆盖度分布、玉米植株实测高度和交叉极化比关系、玉米冠层高度反演结果

极化后的祼土后向散射系数；极化后的祼土后向散射系数； 极化后的地表总后向散系数；极化后的地表总后向散射系数。

图6 VH和VV极化方式下土壤水分预测值与实测值对比

土壤体积含水量反演结果如图7(a)所示，从图中可以看到，土壤体积含水量集中在0～0.30 cm3·cm-3，空白部分为剔除的土壤水分反演异常值区域和人工建设地块区域，如城镇建设用地、公路等。图7(b)则是各个区间范围土壤体积含水量的频数分布情况，可以看到研究区土壤水分存在偏低现象，土壤体积含水量主要集中在0～0.10 cm3·cm-3区间。根据土壤墒情的划分标准[41]：5%以下属于干土；8%左右属于灰墒；12%～15%属于黄墒；>15%～20%属于合墒(褐墒)；>20%属于饱墒(黑墒)。研究区农田土壤水分大部分属于干土和灰墒，作物出现水分亏缺现象，需要人工适时开展灌溉作业，保证农作物的良好生长和收成。

图7 土壤体积含水量(θv)反演结果和频数分布

4 讨论

在前人研究的基础之上，结合微波与光学遥感的各自优势与特点，以景县玉米农田野外实验验证了S1与S2遥感数据协同反演植被覆盖区农田土壤水分的可行性。S1 VV极化模式具有更强的穿透性，可更灵敏地探测植被覆盖区下垫面的土壤水分信息。因此，VV极化条件下土壤水分的反演在R2、RMSE和MAE误差指标方面均优于VH极化方式，这与已有的研究成果相一致[4,27,30]。郭交等[22]基于S1和S2遥感数据反演了陕西省杨凌示范区冬小麦种植区农田地表土壤水分。由于其观测时间正处于冬小麦的生长季节，植株种植稀疏、高度较低，植被对微波遥感的吸收和散射作用不强。水云模型+Oh模型的组合反演模型的R2和RMSE分别为0.650和0.025 cm3·cm-3，反演精度和该研究较为接近。

该研究的误差主要有2个方面：(1)S1微波遥感影像和S2光学遥感影像自身配准误差和校正误差，以及这2种影像数据空间分辨率匹配误差，这些误差都会对最后土壤水分反演精度产生影响。(2)土壤水分和植被含水量的空间分布具有较强的空间异质性，采集的样点数据并不能完全客观、真实地反映10 m空间尺度(遥感影像数据重采样空间分辨率)实际情况。土壤水分样品和玉米植株样品从野外运输到实验室的过程中存在水分蒸发现象，这也会对植被含水量反演模型和土壤水分反演模型的建立产生影响。

该研究存在2个方面的不足：(1)遥感数据在实时性获取和多极化建模方面存在不足。一方面，S1A/B和S2A/B的重复过境时间分别是6和5 d，并不能满足农田土壤水分反演实时性的要求。另一方面，S1A数据只包含VH/VV双极化方式，相比于ALOS-2、Radarsat-2、TerraSar-X等全极化SAR影像数据缺少多极化建模的能力。不同极化的组合方式往往包含多方面的信息，这些信息在土壤水分遥感反演的过程中往往是相当重要的[42]。(2)野外数据采集过程只对玉米种植的农田区域进行了观测分析，实测样点数量有所不足，景县地区农田还种植了一定面积的冬小麦。基于单一地表作物建立的土壤水分反演模型具有一定的局限性，为了提高模型的推广应用价值，今后将在收集大量样点实测数据的基础上开展对种植或混种了多种类型农作物的农田地表土壤水分反演做进一步研究。

5 结论

基于Sentinel-1A主动微波遥感数据和Sentinel-2A光学遥感数据，以归一化水体指数INDW,1610计算植被含水量，并结合改进水云模型、Oh模型以及野外实地考察获取的观测数据，建立河北省景县研究区半经验土壤水分反演模型，实现了微波与光学遥感协同反演植被覆盖区农田土壤水分。

(1)Sentinel-2遥感影像计算的归一化水体指数INDW,1610和实测植被层含水量有着良好的相关性，R2为0.834 3，RMSE为0.737 7 kg·m-2，可利用INDW,1610建立植被层含水量反演模型。

(3)Sentinel-1的VV极化模式下，改进水云模型和Oh模型的组合方法具有较高的土壤水分反演精度，决定系数R2为0.653 0， RMSE为0.040 1 cm3·cm-3，MAE为0.032 7 cm3·cm-3，这3项反演精度评价指标均优于VH极化。

致谢：感谢欧洲航天局提供Sentinel系列遥感数据，北京航天宏图信息技术股份有限公司提供PIE系列遥感软件。