土壤水分是地表水循环、能量循环和生物地球化学循环的关键控制变量，在全球气候系统中扮演着极其重要的角色。作为短期气象预报模型的重要边界条件，它的精确估计能够明显提高天气预报模型的预测能力^[1,2]。诸如土地退化及荒漠化、极端气候事件（如洪涝和干旱）的发生都与土壤水分有着密切的联系。因此，准确获取大范围的土壤水分信息，在灾害预警和监测、生态环境保护以及农业生产等方面均具有重要意义。例如位于中国西北干旱半干旱地区的黑河流域，作为中国第二大内陆河流域，养育着数百万人口，农业和畜牧业是该区域的主要产业。而土壤水分恰恰是农作物和牧草生长的关键因素，也是维系该区域脆弱生态环境的核心因子之一，因此土壤水分的准确估计对该区域显得尤为重要。

土壤水分具有极强的时空异质性，传统的单点测量方法，如重量测量法和时域反射仪法（TDR）无法在短时间内获取大范围的空间分布信息。遥感技术以其独特的大面积观测能力在土壤水分监测方面取得了巨大成功。微波遥感具备全天时、全天候观测的优点，且具有一定的穿透能力，能够探测植被覆盖下的土壤信息^[3]，已成为土壤水分反演的主要数据源，如被动微波的AMSR-E/2、SMOS；主动微波的ERS、 ENVISA/ASAR、ASCAT、ALOS/PALSAR、RADARSAT、Sentinel-1等^[4,5]。其中，Sentinel-1是迄今为止首个向公众免费提供高分辨率SAR数据的雷达卫星计划，该计划由A、B两颗卫星组成，分别于2014年4月和2016年4月发射，两颗卫星组网可提供时间分辨率高达6天的卫星观测。做为ENVISA/ASAR的后继星，Sentinel-1在表层土壤水分反演方面继续受到广泛关注^[6]。因此本文的目标之一就是开展基于Sentinel-1的土壤水分反演研究，探讨其在农田植被覆盖地表土壤水分反演中的潜力。

土壤在微波波段的发射和散射特性均受到其介电特性的影响，而土壤介电常数又与土壤水分的大小直接相关，为微波遥感反演土壤水分奠定了坚实的物理基础^[7]。相对于主动微波遥感，被动微波遥感土壤水分的反演研究起步较早，历史更长，算法也相对成熟。但被动微波遥感的空间分辨率较粗，不适用于较小尺度研究。而主动微波遥感中的合成孔径雷达（SAR）则具有高空间分辨率，能够捕捉细微的地表空间信息特征。但除了受土壤介电特性控制外，主动微波观测的后向散射信号也与地表几何结构（如粗糙度）和植被覆盖密切相关。地物与雷达散射信号之间的相互作用复杂，散射机理混杂（面散射、体散射以及二面角散射相互混合）等，因此基于SAR观测的土壤水分反演方法的发展一直是研究的热点和难点。

在过去的几十年间，基于SAR数据的土壤水分反演取得了诸多进展，形成了多种反演算法，大致可以分为经验方法、物理方法以及半经验半物理几大类。Kornelsen等^[5]对这些方法进行了系统的总结，如变化检测法、优化代价函数法、线性-非线性插值法、查找表法、贝叶斯后验估计法以及人工神经网络等。变化检测方法假设地表粗糙度在短时间内为常量，通过时序接近的影像消除了粗糙度参数对雷达后向散射的影响，建立土壤水分与雷达观测的经验关系^[8,9,10,11]。代价函数优化方法是建立一个代价函数来获取系统误差，然后通过迭代算法最小化该代价函数以获取地表参数的最优估计值。由于优化代价函数法一般基于具有一定物理意义的正向模型和全局优化算法，因此该方法被认为具有较高的普适性和可靠性^{[12,13,14,15]}。机器学习方法主要包括贝叶斯后验估计^{[16,17,18,19,20,21]}、人工神经网络^[6,22,23]和支持向量机^[24]等。此外，植被和地表粗糙度是雷达观测信号重要的贡献因子，也是影响土壤水分反演精度的主要来源，对这两方面的相应参数进行准确估计（或反演）的方法仍需进一步探讨。基于代价函数的土壤水分反演算法能够同步估计多个未知量，包括土壤水分、粗糙度和植被参数。然而在有限的雷达观测（例如双极化）条件下，反演的不确定性随着未知数的增加而增加。增加多源观测对反演算法的代价函数进行有效约束，将有利于提高反演精度，降低不确定性。因此，本文考虑首先利用光学遥感数据对植被参数进行参数化估计，进而利用雷达观测进行土壤水分和粗糙度同步估计，有望降低反演的不确定性，提高反演精度。

因此，本文的目标是开展基于SAR和光学遥感数据的农田土壤水分联合反演研究。在反演方法方面，强调主动微波与光学遥感结合，同步估计土壤水分、地表粗糙度和植被参数。数据方面，采用多时相的Sentinel-1和Landsat 8数据，一方面寄望于获取高精度高分辨率的土壤水分反演结果；另一方面，通过这两种数据协同的土壤水分反演研究，发挥两种观测数据各自的优势，也可为应用国产卫星数据（例如高分3号和高分6号）联合反演土壤水分提供方法上的参考。

以黑河流域（图1（a））中游张掖绿洲（图1(b)）为主要研究区，开展土壤水分反演算法发展和实验研究。研究区主要农作物包括制种玉米、蔬菜和小麦等（图1（c））。绿洲外围被戈壁、荒漠和湿地等包围。该区域气候干旱、年降水量少，灌溉水是土壤水分补给的主要来源。然而，夏季日照强烈、土壤蒸发量大，土壤水分的动态变化又十分剧烈。2012年黑河流域生态水文过程综合遥感联合观测试验在该区域开展，期间布设了大量的生态、水文和气象要素观测系统，积累了丰富的观测数据^[25]。本文选择2012年设置的核心试验区（图1（b）中红框范围）开展土壤水分反演研究，主要基于以下两点考虑，一是该区域的地面观测数据能为算法发展和验证提供数据；二是本区域内地表覆被类型多样，非常具有代表性，可作为异质地表土壤水分反演研究的示范。

Sentinel-1卫星是欧洲航天局（ESA）哥白尼计划中的地球观测卫星，由两颗卫星组成，载有C波段合成孔径雷达，中心频率5.4 GHz。Sentinel-1A拥有多种成像模式，包括Stripe Map（SM）、Interferometric Wide Swath （IW）、Extra Wide Swath （EM）等，本文采用IW模式的地距多视产品（Ground Range Detected, GRD）。该模式数据空间分辨率为5 m×20 m。采用ESA提供的软件SNAP进行了数据处理，包括辐射定标、辐射校正、光斑滤波、地形多普勒校正以及地理编码等。为了和光学数据Landsat 8的空间尺度匹配，将处理后的SAR数据重采样为空间分辨率30 m。需要说明的是，为了尽可能使Sentinel-1和Landsat 8过境时间匹配，同时考虑地面同步观测数据的可获取性（无线传感器网络为5月中旬到9月下旬，自动气象站全年观测），选择Sentinel-1数据的标准为：①假设植被含水量和结构在短期内变化不大，认为Sentinel-1与Landsat 8的过境时间差不超过5天的数据均可匹配；②土壤水分受降水、灌溉、植物蒸腾和土壤蒸发等因素影响，夏季变化剧烈，因此要求地面观测数据与Sentinel-1数据在日尺度上同步。基于此准则，共选用了18对Sentinel-1与Landsat 8影像，获取影像信息详见表1。

Landsat 8是NASA陆地卫星系列中最新一代的光学和热红外遥感卫星，于2013年2月11日发射升空。本文采用其上搭载的陆地成像仪（Operational Land Imager, OLI）传感器的多波段数据估算相关指数，进而反演植被含水量，为土壤水分反演提供辅助。由于获取的L1T级别文件已经通过了几何校正，并已投影到UTM/WGS84系统下，因此对其进行的预处理包括辐射定标和大气校正。本文采用ENVI/IDL的批处理工具ENVITask一次性完成辐射定标和大气校正两个环节。鉴于研究区大气透过率高，大气校正采用快速大气校正方法ENVIQUACTask，可满足研究需要。Landsat 8影像过境日期详见表1。

本文采用黑河流域生态水文过程综合遥感联合观测试验期间建设的地面自动气象站（超级站）和无线传感器网络（Wireless Sensor Network, WSN）节点观测数据对反演结果进行验证。超级站常年提供多层土壤、大气通量观测。地上气象要素包括气温、风速、风向和辐射等，地下观测变量包括土壤水分、温度和热通量等。采用表层（~5 cm）土壤水分观测数据集，其观测频率为10 min，对观测数据进行了日平均，用以验证反演结果。此外，2012年5月中旬到9月下旬，以该站为中心，在周围1 km×1 km范围内架设9套WSN节点，提供两层（4和10 cm）土壤水分、温度观测。数据观测频率也为10 min，本文采用日平均的4 cm土壤水分开展反演算法验证，有效验证样本数信息详见表1。

除此之外，获取了2016年5月和8月多个站点的植被参数测量数据，包括植被高度、生物量、叶面积指数和植被覆盖度等（其中植被含水量的测量采用采样烘干法），用于植被含水量估算结果的验证。详细的测量日期及有效样本数见表1。

植被的干扰是影响土壤水分反演精度的重要因素。理论模型虽然可以较好地刻画微波辐射传输过程，但所需参数多，应用到反演中较为困难。因此，考虑到本研究所收集的植被地面观测数据情况，采用水云模型对植被效应予以校正。该模型是Attema和Ulab为研究农作物而提出的经验模型，有着广泛的应用^[27]。模型假设：①植被是由许多大小和形状相同的微粒组成，可以被看作是一层“云”;②植被—土壤之间的多次散射被忽略。因此，雷达接收到的后向散射信号可表达为：

其中：σtotal0为雷达接收到的总后向散射系数，σveg0为植被散射部分，σsoil0为土壤散射部分；θ为入射角，vwc(kg/m2)为植被含水量；τ2为植被观测透过率，可由植被含水量参数化获得。A和B为模型的经验参数，对比前人的研究，研究区植被状况与文献[28]接近，故本研究中A，B分别取值0.001 4和0.084。通过对式（1）变形可获得土壤表面散射系数为：

早在上世纪70年代，光学遥感数据，如NOAA/AVHRR，Landsat TM/ETM，就已经应用于植被含水量的估算研究中。逐渐形成了三大类方法，包括光谱反射率方法、光谱指数法和辐射传输模型方法。其中光谱指数法相对于反射率方法包含更多的光谱信息，精度高于反射率方法；相对于辐射传输方法，光谱指数法较为简单，因而应用最为广泛。其中闻熠等^[29] 基于修正的土壤调节植被指数（MSAVI）在黑河中游农田开展了植被含水量估算，该方法通过系统分析研究区完整生长季的植被含水量和MSAVI的相互关系，构建了植被含水量与MSAVI的关系模型（式5），取得了理想的结果。鉴于本文研究区和植被覆盖类型与之相同，故本文首先采用该方法开展植被含水量估算。

其中:

此外，Gao等^[30]利用地面实测数据，针对不同作物类型，如玉米、谷物、小麦和草地等，系统评估了多种光谱指数（如归一化植被指数（NDVI）及归一化水体指数（NDWI））及相关的植被含水量估算方法，包括Jackson方法和Chen方法等^[31,32]，在此基础上发展了针对不同作物类型的植被含水量估算方法。本文为充分探讨多种植被含水量估算方法对土壤水分反演算法的影响，采用Gao等^[30]针对玉米的估计方法（式7）：

其中：

上式中，NIR代表近红外波段，RED为红波段，SWIR为短波红外波段。考虑到Landsat 8数据的波段特点，式（8）和（10）选用文献[30]中SWIR中心波长1 640 nm的形式。

在植被—土壤后向散射分解和植被含水量估算的基础上，土壤后向散射系数（公式4）由地表土壤介电属性（介电常数）和几何属性（粗糙度）决定。以往的研究发现，地表粗糙度对雷达后向散射系数的敏感性不亚于（甚至高于）土壤水分对其的敏感性^[33]。因此，本文采用全局优化算法最小化代价函数的策略进行粗糙度和土壤水分同步估计。具体来讲，就是采用裸土表面散射模型模拟与消除植被散射部分之后的雷达观测构造如式（11）的代价函数J，采用复型洗牌优化算法（SCE-UA）^[34]最小化该代价函数，从而获得全局最优解。

其中：σobs,k0为消除植被散射贡献的雷达后向散射系数，由于研究区在2014~2016年只有单极化VV的雷达数据可用，因此σobs,k0为Sentinel-1观测的VV极化后向散射系数。σsim,k0为正向模型的后向散射系数模拟值。SCE-UA算法需要给定待优化参数的取值范围，为此根据研究区的情况，土壤水分的范围设定为0.05 ~ 0.40 m³/m³，地表均方根高度的范围为0.2 ~ 2.0 cm。为了使反演在满足精度要求的前提下尽量快速高效，本文选择半经验模型Oh模型（2004年模型）^[35]进行正向模拟。其一，该模型是由土壤水分等变量及经验系数直接模拟后向散射系数及极化比，避免了由土壤水分模拟介电常数的中间过程；其二，该模型仅用均方根高度（RMSH）一个参数来表征地表粗糙度，使得未知参数个数减小，从而有效降低反演过程中的欠定性；其三，该模型具备与积分方程模型（IEM）^[36]、高级积分方程模型（AIEM）^[37]等理论模型相近的模拟性能，而且计算过程简单，耗时少。

首先分别基于NDVI、MSAVI和NDWI 3个指数进行了植被含水量估算，并利用2016年地面采样数据对结果进行了验证，如图2所示。图中vwc_ndvi、vwc_msavi、vwc_ndwi分别表示基于NDVI、MSAVI、NDWI方案估计的植被含水量。可以看出，3种指数方法的反演结果与实地测量值具有很好的相关性，决定系数（R²）分别是0.898、0.856和0.928，这表明3种指数都能够较好地捕获植被含水量的变化趋势。然而，3种方法都存在总体上的明显低估，其中基于MSAVI方法低估最为明显；相对而言，基于NDWI和NDVI的方法对植被含水量的估计基本能够接近地面观测值。整体上低估的原因可能是由于当植被覆盖度逐渐增大甚至完全覆盖地表时，可见光和红外波段只能观测到植被冠层信息。此时，光谱指数已接近饱和，对植被含水量的变化不再敏感，从而导致对植被含水量的明显低估。另外，需要注意的是，当植被含水量很小时，NDWI方法对植被含水量存在一定程度的高估。这是由于该指数是由受水分胁迫非常敏感的1 640 nm波段计算而来，该波段不仅受植被水分控制，而且受到土壤湿度状况的影响。当植被刚刚出苗时，卫星观测信号不仅包含植被信息，而且还含有土壤表层的信息，例如土壤水分，因此会导致该指数对植被含水量的高估。最后，需要注意到是，本文由于观测数据有限，对植被含水量的反演结果验证可能存在一定的片面性，但从整体来看，NDVI和NDWI方法能够较好地估计植被含水量的变化趋势，获得了较为理想的精度，均方根误差（RMSE）小于0.9 kg/m²。随着地面观测数据和卫星影像的积累，有望提高反演的精度。

在植被含水量估算的基础上，采用迭代优化算法同步估计了土壤水分和地表粗糙度。地表粗糙度的反演结果取值范围在0.5 ~1.5 cm之间，由于没有相应的粗糙度测量数据对其进行验证，这里仅展示土壤水分的反演验证结果。为了系统评价不同植被含水量估算方案对土壤水分反演结果的影响，相应地将3种植被含水量的估计结果代入水云模型中进行植被透过率参数化，进而开展土壤水分反演，并采用地面无线传感器网络和自动气象站观测数据对结果进行验证。如图3所示，总体来看，无论采用何种植被含水量估算方案，土壤水分反演结果与地面观测值呈现了良好的一致性，具有较高的相关性（R²分别为0.564、0.484和0.70）和较低的均方根误差（RMSE分别为0.084、0.076和0.073 m³/m³）。这一方面说明了该反演算法能够较好地进行土壤水分估计，另一方面说明了结合Sentinel-1雷达与Landsat 8多光谱遥感数据进行土壤水分反演的可行性。相比较而言，基于NDWI的植被含水量估算方法能够得到最优的土壤水分估计结果，其具有最小的均方根误差和最大的决定系数。这表明准确估计植被含水量或准确地参数化植被透过率对土壤水分的反演算法至关重要。

在土壤水分反演算法充分验证的基础上，以黑河中游试验核心观测区2016年5月29日观测数据为例，获得的土壤水分空间分布如图4所示。可以看出，该区域内土壤水分呈现强烈的异质性，从0.05 m³/m³的极干状态到0.43 m³/m³的接近饱和状态。空间分布上，在研究区东北偏西方向上土壤水分呈现明显的高值，西边及西北部分区域土壤水分值相对较低。结合灌溉信息定性分析发现高值区为夏季第二轮灌溉区域，而干燥区域尚未进行该轮次灌溉，距上次灌溉过去20多天，强烈的夏季辐射导致土壤水分严重蒸发，故而呈现干燥状态。因此，可以认为该反演获得的土壤水分空间分布合理，体现了该区域由于灌溉导致的土壤干湿状况的空间分布特征。

从上述验证结果可以看出，本文提出的结合Sentinel-1和Landsat 8 遥感数据进行土壤水分反演的方法取得了较好的效果。尽管用于验证的地面观测数据较为有限，但其结果具有一定的代表性，可为相关研究提供一定的参考。基于此，对于本研究中的几个细节问题，值得进一步予以讨论。

首先，除了植被含水量的估算方案对土壤水分反演结果的影响外，水云模型中的经验参数A和B对土壤水分的反演结果也有影响^[38]。本文研究区的植被类型虽然多样，但玉米占有主导性地位，故采用文献[28]中所推荐的参数，而在反演实验中未专门测试不同的参数得到的反演结果。为进一步分析参数A和B的影响，以水云模型为框架，针对两个参数进行了局部敏感性分析，分析了雷达后向散射系数对参数变化的响应程度（图5）。可以看出，随着参数A的变化，交叉极化后向散射系数呈现明显的上升趋势，同极化后向散射系数未见明显变化。而随着参数B的变化，同极化和交叉极化雷达后向散射系数均呈现明显的下降趋势。由此可见，这两个参数对后向散射系数均较为敏感，在采用水云模型进行植被冠层散射建模及土壤水分反演时，需尽可能获得它们的“真值”，以确保模拟或反演的准确性。但这两个参数又并非可直接观测的物理量，因此，可能的途径是进一步发展针对参数A和B的像元尺度参数化或标定方案，特别是考虑到不同的植被种类以及不同生长期的植被状态。

其次，迭代循环设置是利用迭代优化算法进行参数估计时需要考虑的问题，因为样本量的大小一方面决定了估计精度的收敛问题，另一方面也涉及到计算耗时。本文在开展正式反演之前测试了不同迭代次数对SCE-UA算法反演精度和计算耗时的影响，最终发现当迭代次数在100左右，反演结果的均方根误差趋于稳定收敛。当迭代次数小于此值时，反演误差就呈现明显的波动；但当迭代次数大于此值时，误差不会明显减小而计算耗时却又增加。当然，计算耗时也取决于正向模型的复杂程度和收敛标准的确定。本文采用相对简单的Oh模型进行正向模拟，反演的研究区较小，影像的空间分辨率为30 m，因此计算耗时问题并不明显。但我们建议在采用SCE-UA算法进行参数估计或者遥感地表参数反演时，特别是面向大尺度的土壤水分产品生产时，有必要首先进行迭代次数对误差收敛性和计算耗时的敏感性测试。

此外，地表粗糙度的影响一直是SAR反演土壤水分面临的核心难点之一。本文采用粗糙度和土壤水分同步估计的反演策略也是基于粗糙度对雷达后向散射系数高度敏感的现实考虑。但由于本文没有相应的粗糙度观测数据对反演结果进行验证，因此未对粗糙度反演结果的精度加以讨论。事实上，近期的研究发现^[39]，即便有地面粗糙度观测数据，由于地表均方根高度、特别是相关长度的尺度依赖性和尺度效应，反演的粗糙度与地面实测也不一定吻合，可能的解决办法是继续发展与遥感观测尺度匹配的粗糙度参数估计方法。

本文结合Oh模型、水云模型和复型洗牌全局优化算法，发展了一种土壤水分、植被含水量和地表粗糙度协同反演方案，采用Sentinel-1 SAR和Landsat 8 多光谱数据开展了实验研究，并采用黑河中游的自动气象站和无线传感器网络观测数据、地面人工观测数据分别对反演的土壤水分和植被含水量进行了验证。主要取得如下结论：

（1） 在迭代优化算法的框架下，结合主动微波和多光谱遥感数据，可同时估算植被含水量、土壤水分和粗糙度，从而有效避免了由于植被参数和粗糙度导致土壤水分病态反演的问题。

（2） 本文提出的土壤水分反演算法取得了理想的土壤水反演结果（RMSE≤0.073 m³/m³），与地面观测值具有较高的相关性（R²≥0.7）。

（3） 土壤水分反演精度依赖于植被含水量的准确估计，本文的案例发现基于NDWI的方法能够较好地估计植被含水量，与地面实测数据相比，R²=0.928，RMSE为0.885 kg/m²，进而有效地改善了土壤水分的估计精度。

当前诸多的地学应用，如环境保护监测、精准农业等，都需要高精度的地表土壤水分产品。传统上，受限于成像幅宽和重访周期，很难通过SAR观测直接反演的办法得到较高空间分辨率的土壤水分实用产品，而Sentinel-1双星系统的运行为此提供了契机和可能^[40]；甚至，通过协同多源SAR传感器的时序数据，有望获得高空间（< 50 m）和高时间分辨率（<3 d）的业务化土壤水分产品^[41]。除此之外，要获得质量可靠的土壤产品，还需要进一步发展相关的反演方法，特别是对植被和地表粗糙度的处理^[42,43]，进而能将相关的研究成果应用于国产卫星数据（如高分系列），服务于国家社会经济发展。