融合合成孔径雷达和光学影像的植被覆盖区土壤湿度反演

谢 媛， 严颂华,2*， 马秦生， 陈能成

(1.武汉大学电子信息学院， 武汉 430072； 2.武汉大学宇航科学与技术研究院， 武汉 430072；3.武汉大学测绘遥感信息工程国家重点实验室， 武汉 430072)

在农业生产中，土壤水分对生态环境和水平衡起着至关重要的作用，也是研究土地和大气之间物理过程的重要变量[1-4]。然而由于土壤湿度在时间、空间上具有时空易变性[5]，对土壤水分进行大尺度、高精度且连续的时空监测成为各大学者争相研究的热点。合成孔径雷达(synthetic aperture radar,SAR)在估计土壤参数方面展现了良好的优势[6-11]。对于裸土下土壤湿度反演模型，到目前为止，已形成了多种反演算法[12]，大致分为经验方法、物理方法以及半经验半物理方法几个主要算法。最常用的是IEM(integral equation model)、CIEM(calibrate integral equation model)、AIEM(advanced integral equation model)、Oh和Dubois模型。这些模型的共同之处在于建立SAR数据和土壤参数之间的函数。同时，结合光学图像等多源数据，有助于提高SAR反演土壤湿度的精度。

目前的研究方法中，包括主被动微波遥感相结合共同改进土壤湿度数据测量方法[13]，以及基于被动遥感与光学遥感相结合的手段进行降尺度获取更高分辨率的土壤水分数据[14]取得了一定的进展，但是这些方法大多依赖于额外植被指数产品，缺乏光学遥感数据本身在植被方面突出表现的考虑。实际上，在地表几何结构的影响中，植被的遮挡使得SAR接收到的垂直极化后向散射系数不仅包含了土壤表层后向散射系数，还包含了植被层散射系数。而对于半干旱，或有农田遮挡区域消除植被的影响能更准确地获取土壤表面水分情况。

针对这一问题，依据主动微波遥感与光学遥感数据相结合，利用光学遥感图像中反映的植被层信息，校正土壤水分反演中的植被层散射贡献。其创新性在于既发挥了光学遥感高分辨率特性以及反映地表参数的优势，又充分结合了主动微波遥感全方位、长时间序列的观测特性，实现对土壤水分准确、大面积反演。与此同时采用经验模型来作为土壤湿度反演算法更具有理论与试验依据。

首先给出了总体技术路线，其次是数据处理部分，然后进行算法研究，其中包括用于求取植被覆盖区地表后向散射系数的植被影响消除算法以及用于得到土壤湿度反演结果的裸土下土壤湿度反演算法，之后是基于算法研究给出的试验结果及分析，最后给出结论。

1 总体技术路线

使用Sentinel-1A卫星作为主要数据源，Sentinel-2多光谱数据作为次要数据源，对武汉市豹澥区域的土壤湿度进行反演研究。如图1所示为总体技术路线。

图1 技术路线流程图Fig.1 Flow chart of technical route in this paper

首先对SAR数据以及光学影像进行处理，然后进行植被影响的消除。即从光学影像中提取各个区域植被覆盖度以及各项植被指数。利用考虑了植被覆盖度的修正水云模型，对每个像素点的纯粹土壤表面后向散射系数进行求解。相比于传统的水云模型去除植被影响办法[15-18]，考虑植被覆盖度因素，使得对于每个像元，植被覆盖部分的影响可以达到定量分析，将垂直极化后向散射系数由原始的植被层后向散射系数与经过衰减之后的裸土后向散射系数2个组成部分变成了3部分的组合——植被覆盖部分植被层散射系数、植被覆盖部分经过衰减的纯粹地表散射系数、非植被覆盖部分纯粹地表散射系数。最后得出纯粹地表后向散射系数。

将纯粹地表后向散射系数代入由Dubois模型和Oh模型结合的裸土土壤湿度反演模型。因此对校正了植被层散射贡献的研究区域进行土壤湿度反演，从而得到时间序列上与实测土壤湿度更为接近、更为准确的土壤湿度反演趋势。

2 数据处理

2.1 SAR数据的预处理

Sentinel-1是欧洲航天局哥白尼计划中的地球观测卫星，由两颗卫星组成，载有C波段合成孔径雷达，本文中下载的数据模式为IW SLC模式，空间分辨率为5 m×20 m。图2所示为2020年3月3号的Sentinel-1A数据干涉宽幅模式的斜距单视复数产品图。

对Sentinel-1号数据预处理过程如下：在ESA的专用软件SNAP里面对Sentinel-1号数据进行辐射定标、多视、滤波以及几何校正，经过预处理后从SAR影像中可以得到每个像元的两个参数信息：垂直极化后向散射系数σ和入射角λ。

2.2 光学数据处理

Sentinel-2号是高分辨率多光谱成像卫星，携带一枚多光谱成像仪，可覆盖13个光谱波段。对光学数据进行处理目的是对原始影像进行大气校正，再利用归一化水体指数(normalized difference water index, NDWI)得到研究区域非水体影像。因此，在常规预处理基础上，添加对光学数据进行水体裁剪的处理部分。光学数据处理分为两步：常规预处理、裁剪水体。

2.2.1 常规预处理

光学数据的常规预处理主要对影像进行大气校正，校正光谱分辨率。利用Sen2Cor处理库对Sentinel-2 L1C级数据进行大气校正，得到L2A级数据，采取10 m分辨率进行处理，导出ENVI文件。

2.2.2 裁剪水体

由于对于土壤湿度的反演工作针对的是陆地部分，若不进行水体的去除，则在反演过程中，水体会被认为是具有植被覆盖部分，因此，水体的存在会对反演土壤湿度产生干扰。利用光学数据进行水体的去除手段是利用归一化水体指数NDWI进行决策树的分类。首先在ENVI(environment for visualizing images)软件中，利用Layer Stacking模块进行红、绿、蓝和近红外四个波段的融合，然后计算NDWI，计算公式为

(1)

式(1)中：Rgreen和Rnir为绿波段和近红外波段的光谱反射率；σNDWI为归一化水体指数。

根据光谱图像NDWI，将其进行分类，即在光学影像中可以认为σNDWI<0的部分为非水体部分，σNDWI≥0的部分为水体部分，再进行决策树的分类，得到非水体部分shp图，对原始图像进行剪裁，即可得到去除水体之后的图像，也就是光学影像处理完最后结果。

剔除水体之后的影像即为下一步处理的基础影像，其中包含了四个波段的反射率Rred、Rgreen、Rblue、Rnir。

3 算法研究

实验中所涉及的算法主要分为两个步骤：第一步旨在通过修正的水云模型，从垂直极化后向散射系数中消除植被对于后向散射系数的影响；第二步则通过Dubois模型和Oh模型结合，利用介电常数与土壤湿度之间经验关系，结合第一步求解的纯粹土壤表面散射系数反演植被区域土壤湿度。

3.1 植被影响的消除

经过对Sentinel-1A SAR图像进行预处理得到垂直极化后向散射系数σ与入射角θ，以及对光学影像数据处理得到新的光学影像之后，可以利用3.1.1节进行纯粹土壤表面散射系数的求取，3.1.2节是实验结果示例，3.1.3节则利用实验结果分析了植被对于后向散射系数的影响与植被覆盖度的相关性。

3.1.1 植被影响消除步骤

将预处理之后的光学影像作为输入，求解纯粹土壤表面散射系数步骤如下：

(1)利用波段反射率数据求解归一化植被指数(normailized differential vegetation index,NDVI)。作为植被影响因素的重要参数，求解公式为

(2)

式(2)中：Rred和Rnir分别表示红波段和近红外波段的光谱反射率；σNDVI代表归一化植被指数。

(2)利用NDVI指数，根据植被覆盖度的定义求出每个像素点的植被覆盖程度σVFC，计算方式为

(3)

式(3)中：σVFC为每个像素点的植被覆盖程度参数；σNDVI,s为没有植被覆盖的像素点纯粹裸土表面NDVI；σNDVI,veg为完全由植被覆盖的像素点NDVI。本文定义σNDVI,s是图像中累积概率为5%的NDVI，而σNDVI,veg是累积概率为95%的NDVI。

对于2020年3月20日的武汉小部分区域植被覆盖度如图3所示。

图3 2020年3月20日武汉小部分区域(包含豹澥)植被覆盖度分布图Fig.3 Distribution map of vegetation coverage in a small area of Wuhan (including Baoxie area) on March 20, 2020

(3)利用NDVI以及VFC求解植被层散射系数以及去除植被之后纯粹地表散射系数。

水云模型定义垂直极化后向散射系数由植被后向散射系数和植被下地表后向散射系数组成。可通过光学图像求解每个像元植被覆盖程度，考虑植被覆盖度VFC，将一个像元看作是一部分植被覆盖区域和另一部分裸土区域。因此本文中修正的水云模型定义垂直极化后向散射系数σ由以下三个部分组成：植被覆盖地区的植被散射系数和经过衰减之后的植被覆盖下纯粹地表散射系数以及裸土区域的纯粹地表散射系数，表现形式如图4所示。

图4 单个像素区域后向散射系数分布示意图Fig.4 A pixel grid back scattering coefficient distribution model of modified water cloud model

由上述分析以及模型示意图，可以得出修正水云模型形式为

σ=(1-σVFC)σs+σVFCσveg+σVFCτ2σs

(4)

式(4)中：双层衰减因子(透过率)τ2和植被层后向散射系数求解公式[19]为

τ2=exp[-σNDVI/cosθ]

(5)

σveg=σNDVIcosθ(1-τ2)

(6)

式(4)中第1项为模型图给出的非植被覆盖地区纯粹地表散射系数；第2项为模型图中植被层散射系数；第3项为经过植被双层衰减后纯粹地表散射系数。σs为需要求解的纯粹土壤表面散射系数。

(4)最后，可以根据SAR雷达影像总体后向散射系数σ和入射角度θ，以及由光学影像参数得到植被覆盖区域纯粹地表散射系数表达式为

(7)

3.1.2 植被消除影响试验示例

为验证3.1.1所述植被影响的消除算法，对试验场地豹澥区域做了相关性研究。在豹澥区域布设的九个原位土壤湿度计实测站点情况如图5所示。

图5 原位土壤湿度计布设情况Fig.5 Layout of in situ soil hygrometer

其中每个站点相对于原始的SAR图像中的垂直极化后向散射系数，去除植被散射系数之后都有一定程度的衰减，将光学图像中每个站点于2020年3月20日的植被覆盖度与在时间序列上植被影响的极化值平均衰减量做相关性验证，所有站点中由于站点3、6、7原位土壤湿度计电池没电故不予考虑。实验数据主要由2019年11月28日至2020年4月20日这一时间段，组成时间序列。以2020年3月3日为例，利用植被消除算法对每个站点垂直极化后向散射系数σ进行植被影响的消除，得到纯粹地表后向散射系数σs，对于每个站点，坐标以及去除植被前后垂直极化值、入射角情况如下表1所示。

3.1.3 植被消除影响实验分析

下面分析植被覆盖度与时间序列各站点垂直极化后向散射系数平均衰减量相关性。其中时间序列垂直极化后向散射系数平均衰减量用Δσ表示，植被覆盖度选取的是2020年3月20日光学影像得到的每个站点的VFC，各站点Δσ计算方法如式(8)所示，计算结果如表2所示。

表1 2020年3月3日植被影响消除算法处理结果Table 1 Processing results of vegetation impact elimination algorithm on March 3, 2020

(8)

式中：求和范围j为每个站点在时间序列2019年11月28日到2020年04月20日上的12个周期；σj为第j个周期站点i的后向散射系数，即植被消除算法输入；σj,s为第j个周期站点i的植被消除算法输出。

绘制植被覆盖度与垂直极化值平均衰减量Δσvv相关曲线如图6所示。

表2 原位土壤湿度计站点植被覆盖度与ΔσTable 2 Vegetation coverage of in-situ soil hygrometer sites and Δσ

图6 植被覆盖度与垂直极化平均衰减量相关性曲线Fig.6 Correlation curve between vegetation coverage and average attenuation of vertical polarization

由表2、图6可知，站点植被覆盖度与时间序列植被对垂直极化值平均衰减量有着非常强的相关性，皮尔森相关系数达到了0.83，因此，对于植被覆盖大的区域，更有必要去考虑植被的影响并消除。

3.2 裸土下土壤湿度反演

对于裸土下土壤湿度的反演，利用Dubois模型和Oh模型作为介电常数和均方根高度的求解方法。其中，Dubois模型给出了垂直极化值σvv和水平极化值σhh分别关于介电常数和土壤均方根高度的表达形式。Oh模型则给出了σvv和σhh之间的关系表达式。因此通过结合这两个模型，首先对介电常数和均方根高度进行求解，其次，利用介电常数和土壤湿度之间的经验方程，对裸土下土壤湿度值进行反演。

3.2.1 求解介电常数和均方根高度

(1)Dubois模型形式为

(9)

(10)

式中：k为波数；λ为波长；θ为入射角；ε为介电常数；s为均方根高度。

为以Dubois模型为基底建立方程组，需要求解水平极化值。对于水平极化值，引入Oh模型，依据Oh模型中垂直极化值与水平极化值比值来求解。

(2)Oh模型的表达形式为

(11)

(12)

式中：

(13)

(3)建立方程组。结合Dubois模型[式(9)、式(10)]以及Oh模型[式(11)]，也就是垂直极化值与水平极化值的比值，可以得到关于介电常数和均方根高度的二元方程，即

(14)

(15)

根据式(14)、式(15)组成的方程组，代入后向散射系数和入射角，即可求解出均方根高度s和介电常数ε。

3.2.2 根据介电常数反演土壤湿度

根据介电常数ε和土壤湿度Mv的经验模型对土壤湿度进行反演，其中经验模型为

(1.993+0.002a+0.015b)+(38.086-

0.176a-0.633b)Mv+(10.720+

(16)

式(16)中：a为黏土含量,%；b为沙粒含量,%，a、b在不同频段下取值不同，本文中取频率为5.55 GHz时含量，即a=0.515 1,b=0.134 3。

由此，即给出了利用垂直极化后向散射系数和入射角信息求解裸土区域土壤湿度的反演算法。

4 土壤湿度反演实验及验证分析

为了验证土壤湿度反演趋势的正确性， 选取2019年11月28日至2020年4月20日这一时间段进行分析，总共12景数据，其中，有2月20号1幅影像元数据丢失。土壤湿度反演步骤如第1节中技术路线所述，具体步骤如第2、3节详细说明。

与实测土壤湿度比较结果分为两个部分：时间序列土壤湿度反演趋势结果分析和空间序列土壤湿度反演结果分析。

4.1 时间序列土壤湿度反演结果分析

根据以上处理结果，对所选时间序列土壤湿度进行反演。进行单个站点时间序列方面的分析。表3以站点1为例，列举时间序列上反演的土壤湿度信息，反演模型输入部分包括垂直极化后向散射系数、植被影响消除后纯粹土壤表面散射系数、入射角，代入模型可以得到植被影响去除前后土壤湿度反演值。

其他站点情况如下，绘制其他站点Mv和Mvs分别与SM的相关曲线。

图7所示为6个站点去除植被前后土壤湿度反演曲线与实测土壤湿度变化曲线的相关性，可以得出以下结论：

表3 站点1时间序列土壤湿度反演结果以及实测信息Table 3 Soil moisture inversion results and measured information on time series of site 1

(1)土壤湿度反演趋势与原位湿度计的实测值变化趋势都有很高的相关性。大部分站点皮尔逊相关系数均达到了0.6左右，即强相关。对于站点2 相关系数达到了0.8以上，呈现极强相关趋势。

(2)相比于植被影响消除之前，去除植被影响后反演的土壤湿度值与实测值变化趋势更相关，即证明实验结果对于利用SAR影像反演土壤湿度的反演结果有了很好的提升。也证明利用Sentinel-1 SAR影像和Sentinel-2多光谱数据结合来反演土壤湿度一定时间内变化趋势相比单利用SAR微波后向散射系数反演具有更高可信度。

4.2 空间序列植被去除前后土壤湿度反演图像结果分析

图8所示为2020年3月3日豹澥区域植被去除前后土壤湿度反演结果。

图8 植被去除前后土壤湿度反演结果对比Fig.8 Comparison chat of soil moisture inversion results before and after vegetation removal

由于豹澥区域植被覆盖面积较少，选取附近红色区域发现相比于去除植被影响之前总体后向散射系数反演的土壤湿度，去除植被影响后，对于土壤湿度黄土部分和植被覆盖代表的湿部分边界更为明显，整块区域有明显的几何形状，而未去除前的土壤湿度的区域划分则很模糊。

同时，对于植被覆盖度越大地区，理论上植被对于土壤湿度反演影响则越大，因此，绘制每个站点P2与P1差值曲线，与站点植被覆盖度做相关性分析，结果如图9所示。

图9 植被覆盖度与土壤湿度反演提升效果相关曲线Fig.9 Correlation curve of VFC and soil moisture inversion lifting effect

由图9可以看出，对于植被覆盖度越大的地区，去除植被前后土壤湿度与实测土壤湿度相关性差值越大，且两者相关性达到了0.8以上，及证明试验结果符合理论。

5 结论

利用Sentinel-1 SAR影像数据与Sentinel-2光学影像数据对土壤湿度反演过程中的植被影响进行了研究。通过对SAR图像进行预处理，提取所研究区域的垂直极化后向散射系数以及入射角信息，同时利用Sentinel-2号的光学影像求解归一化植被指数NDVI以及植被覆盖度VFC，得到植被层散射系数和垂直极化衰减量，代入修正水云模型求解出植被覆盖区域纯粹地表后向散射系数，最后利用裸土下土壤湿度求解算法，对研究区域进行土壤湿度的反演。证实所述方法对比直接利用SAR影像垂直极化后向散射系数反演土壤湿度，提升了与实测值的相关系数，对于利用SAR影像反演土壤湿度的研究领域有着更好的参考价值和实际意义。

方法不足之处在于：没有对植被的类型进行更进一步的研究，对于不同植被覆盖产生的影响只是做了统一性处理，其次，在分辨率方面，最终结果采用的是5 m×20 m的空间分辨率，因此对于更加精细的区域进行土壤湿度的反演还存在一定的研究空间。在该论述方案中，对于空间序列方面土壤粗糙度的影响缺乏探讨，克服这些问题，将是下一步研究的重点。