青藏高原作为地球第三极，由于其独特的高海拔气候条件，成为中低纬度地带面积最广的冻土分布区^[1,2,3]。当土壤温度低于0 ℃，土壤会发生冻结，但由于受到土壤水表面张力以及土壤水中溶质的影响，冻土中的水分始终会保持一定数量的液态形态，称为未冻水。未冻水和冰共同存在于冻土中，两者的相互转化即冻融变化对地表水分循环和能量收支有着重要影响。地表冻融伴随着大量的能水转化，直接影响土壤中的水热分配，进而影响地气之间的水热交换^[4,5]。此外，冻结土壤的导水性能大幅下降，从而形成特殊的产汇流模式，改变区域水循环特征^[6,7,8]。因此，开展地表土壤未冻水及冻融变化的监测对于解析气候变化情景下的冻土水文变化特征具有重要科学意义。

地面土壤湿度观测能够提供较为精确的土壤未冻水以及冻融变化信息，因此近年来Su^[9]和Yang^[10]分别在青藏高原不同气候和下垫面区域建立了卫星尺度的土壤湿度观测网络。但由于强烈的空间异质性，当前针对特定区域的地面观测网络远远不能满足整个高原土壤未冻水以及冻融变化分布特征的监测需求。微波遥感技术因其能够反映地表土壤水分以及冻融变化引起的土壤介电性质差异，因而具备了在全球或区域尺度上监测土壤水分以及冻融变化的潜力^[11,12,13]。

随着SMOS（Soil Moisture and Ocean Salinity）^[14]和SMAP（Soil Moisture Active Passive）^[15]卫星的相继发射成功，全球尺度下的土壤水分产品逐渐推出，同时也引发了利用L波段被动微波遥感监测地表土壤水分和冻融变化的研究热潮^[16]。然而，当前基于L波段被动微波遥感发展的土壤水分产品（如SMOS和SMAP）缺乏冻结土壤的水分变化信息，而冻融监测产品仅能表征地表土壤冻结或融化状态，不能反映冻土中未冻水的动态变化特征。

冻土中未冻水微波遥感监测产品的缺失主要是由于当前广泛使用的土壤介电常数模型仅适用于非冻结土壤，比如Dobson^[17]和Mironov^[18]等发展的模型。考虑到适用于冻结土壤的介电常数模型的发展能够更好地帮助构建更具物理机制的冻融判别标准^[19]，近年来关于冻土介电常数模型的研究也逐渐增强。在冻土介电常数的计算中，Zhang等^[20]考虑了土壤质地、比表面积和土壤未冻水含量之间的联系，Schwank等^[21,22]证实了Birchak等^[23]初步建立的四相混合介电常数模型在冻土中的适用性。近期，Mironov^[19]发展了一种适用于L波段（1.4 GHz）冻结和非冻结土壤的介电常数模型。Zheng等^[24]评估了四相混合介电常数模型和Mironov模型在高原冻结土壤的适用性，发现四相混合介电常数模型更适合高原土壤特征。

对于青藏高原地区，Zheng等^[24]研究发现当前SMAP卫星采用的微波辐射参数化方案高估了植被效应，低估了地表粗糙度影响，由此发展了新的植被和地表粗糙度参数化方案。此外，Zheng等^[24]进一步结合四相混合介电常数模型基于SMAP垂直极化亮温观测反演了冻结和非冻结土壤中的未冻水（θliq）变化特征。然而，这些研究仅使用单通道SCA（Single-Chanel Algorithm）算法对土壤未冻水进行反演，且数据类型局限于SMAP卫星降轨时刻的亮温观测。因此，本文将同时利用SMAP升轨和降轨时刻的亮温观测数据来反演青藏高原玛曲区域的土壤未冻水变化特征，以确定更适用于该地区土壤未冻水反演的数据类型；并通过比较SCA和双通道DCA（Dual-Chanel Algorithm）算法所得结果的反演精度，得到更好的土壤未冻水反演方法。

研究区位于黄河源区东南部玛曲县（如图1所示），该区域的冻土类型为季节冻土，气候类型为高寒湿润型，气候的特点是冬季寒冷干燥，夏季多雨，年平均气温为1.2 ℃，最冷的月份（1月）平均气温为-10 ℃，最暖的月份（7月）平均气温为11.7 ℃，年均降水量615.5 mm，年均蒸发量1 347.3 mm^[25]。该地区土壤类型主要为粉砂壤土，植被类型主要为高山草甸。

2008年，Su等^[9]在玛曲区域建立了由20个观测剖面构成的土壤温、湿度观测网络。每个观测站点分别在土壤深度5、10、20、40和80 cm处安装了5TM ECH₂O探头用于监测土壤温湿度，并且每15 min记录一次观测读数。当前，该土壤温湿度观测网络已被SMAP卫星任务选定为其土壤水分产品的地面验证点^[9]。SMAP卫星团队选定的验证格网如图1中红框所示，本文选取该验证格网内2016年8月至2017年7月间所有站点5 cm深度的观测平均值用于反演产品的评估和验证。近期，Zheng等^[33]对比评估了冻结时期5TM土壤未冻水观测与Noah陆面模式模拟结果，发现模式模拟与观测的土壤未冻水动态变化一致，从侧面上验证了5TM观测在冻结土壤的有效性。因此，本文直接采用5TM观测用于土壤未冻水反演结果的评估。

本文采用的是SMAP辐射计产品L1C_TB，该产品可以通过美国国家冰雪数据中心（NSIDC,https://nsidc.org/data/smap/smap-data.html）下载获得。选取了2016年8月至2017年7月SMAP卫星降轨（上午6点过境）和升轨（下午6点过境）时刻的360个亮温观测数据，亮温观测值取自覆盖于研究区土壤温湿度观测网络的SMAP像元（即图1红色框）。

本文所采用的前向模拟模型为τ-ω模型，使用的反演算法包括单通道SCA和双通道DCA算法，方法详细描述如下。

τ-ω模型是使用最为广泛的零阶正向微波辐射传输模型^[26]，该模型假设植被为均匀且温度为常量的整体层结发射极化微波。τ-ω模型中的亮温TBp由3部分组成：第一部分为植被直接发射的微波辐射;第二部分为植被发射经地面反射的微波辐射;第三部分为土壤发射经植被衰减后的微波。模型表达如下^[27]：

其中：上标p代表极化方式（水平极化H；垂直极化V），τp为植被光学厚度，ωp为有效散射反照率，TC和TG分别表示植被和土壤的有效温度，γp为植被透射率，rp为粗糙土壤反射率，ψ是入射角，与SMAP的观测角度是相同的，即ψ = 40°。

当前SMAP算法中τp与植被含水量（VWC）线性相关，其中植被含水量由归一化植被指数NDVI计算，而ωp的值由特定植被类型决定^[28]（如本文研究区为草地，ωp取值为0.05）。Zheng等^[24]研究发现SMAP的τp计算方案高估了植被的影响，造成青藏高原草地亮温模拟的高估，由此发展了新的参数化方案，其表达式如下：

假设SMAP过境期间空气、植被和近地表土壤处于热平衡状态，则TC近似等于TG，两种温度可以用单一有效温度Teff代替^[29]，因此式（1）可以重新表达为：

本文采用Choudhury等^[30]给出的理论公式来估计Teff：

其中：Ts(z)是在深度z的温度，衰减系数α(z)与该深度土壤介电常数实部εs'(z)和虚部εs''(z)及观测波长λ相关，其表达式如下：

粗糙地表反射率rp的计算方法如下^[31]：

其中：rsp和rsq分别表示光滑表面对不同极化微波的反射率（其中p=H时，q = V；p=V时，q = H），该反射率可以通过菲涅尔方程计算获得^[26]，h表示有效粗糙度因子，Q表示交叉极化因子，N表示角度校正常数。当前SMAP算法假设h与均方根高度s线性相关，s根据地表类型赋值，Q取值为0，N取值为2。Zheng等^[24]研究发现SMAP的h和Q的计算方案低估了地表粗糙的影响，并引进了新的参数化方案，其表述方式如下^[32]：

其中：s取值为15.6 mm^[24]。

考虑到Zheng等^[24]的研究同样针对玛曲地区，本研究将采用由此发展的新的植被（即式（3））和地表粗糙度（即式（8）-（9））参数化方案。

光滑表面微波反射率主要采用菲涅尔方程进行计算：

最近，Zheng^[22]等证实了四相混合介电常数模型^[23]对玛曲地区冻土介电常数εs估计的适用性，因此本文采用四相混合介电常数模型计算εs。四相混合介电常数模型将湿润的土壤组分分为空气、液态水、冰和土壤基质：

其中：θs是土壤孔隙度，θ和θliq分别代表土壤总含水量和液态含水量。指数η等于0.5，空气，冰和土壤基质的介电常数分别为εair = 1、εice = 3.2 +i⋅ 0.1及εmatrix = 5.5 +i⋅ 0.2^[21]。考虑到土壤温度的影响，使用Dobson^[17]模型估算液体水的介电常数（εw）。对于未冻结土壤，θ和θliq的值相等，该值可由土壤湿度观测网的5TM探头观测获取。对于冻结的土壤，θliq的值通过5TM探头观测得到，θ的值则通过对冻融前后θliq实测值的线性插值估算获取^[22]。

本文的反演算法流程如图2所示，通过多次迭代计算，输出亮温模拟值（TBsim,ip，i为迭代次数）与SMAP观测值（TBsmapp）之间代价函数（CFi）最小时CFmin对应的土壤未冻水含量（θliq）为反演值，算法具体描述如下。

输入一个SMAP卫星亮温观测值（TBsmapp）并配置好τ-ω模型（3.1节）的相应参数，分别给定θliq和CFmin的初值为0 m³m^-3和300 K²，将θliq代入τ-ω模型计算得到本次迭代模拟的TBsim,ip，并将其与TBsmapp差的平方赋值给CFi。若CFi小于CFmin，将CFi重新赋值给CFmin。然后θliq自加0.001 m³m^-3，进行下一次迭代，直至θliq加到0.5 m³m^-3即迭代500次后跳出循环，最后统计与CFmin对应的θliq即为本次土壤未冻水含量的反演值。

图3（a）和3（b）分别展示了2016年8月至2017年7月SMAP卫星降轨和升轨过境时刻玛曲地区土壤未冻水含量反演值与实测值（Obs）的时间序列。其中反演值分别基于垂直极化亮温观测的单通道算法（SCA-V）、水平极化亮温观测的单通道算法（SCA-H）和双通道算法（DCA）获取。图4进一步以散点图的形式展示了图3中3种反演值与实测值的相关性。表1为土壤未冻水反演值与实测值的误差统计分析，统计指标包括无偏均方根误差（ubRMSE）、偏差（Bias）、均方根误差（RMSE）和相关系数（R）。

对于降轨时刻，结合图3（a）、图4（a）~（c）以及表1发现：3种土壤未冻水的反演值均较为同步地反映了实测值的动态变化，反演值与实测值相关性较强，相关系数R均大于0.9，其中SCA-V反演值与实测值的相关系数最高（R=0.92）。然而，从图3（a）可以看出，基于SMAP降轨时刻的反演结果在10~11月和4~5月等冻融交替季节对观测值存在明显低估。上述情况在图4（a）~（c）中表现为较多实测值在0.2~0.4 m³m^-3范围内的散点明显向右偏离1∶1斜线（图中红色实线）。表1结果显示，3种结果的ubRMSE均大于0.05 m³m^-3，未达到SMAP任务的预期目标^[29]（即ubRMSE≤0.04 m³m^-3）。

对于升轨时刻，结合图3（b）、图4（d）~（f）及表1发现：3种土壤未冻水的反演值均同步反映了实测值的动态变化，但在1~2月的反演结果稍有低估，其中SCA-H低估最多，并且SCA-H的反演结果在3~6月出现4次高估的峰值，说明其曲线的波动起伏较大，SCA-V和DCA的反演结果相对稳定。图4（d）~（f）散点图中各趋势线（黑色虚线）与1∶1斜线几乎重合，体现了反演值与实测值极强的相关性，各相关系数（R）均大于0.9，其中SCA-V和DCA的相关系数最高（R=0.95），且ubRMSE的值分别为0.035 m³m^-3和0.039 m³m^-3，达到了SMAP任务的预期目标。

对比降轨和升轨时刻土壤未冻水反演结果的主要特点，发现升轨时刻的反演结果避免了图3（a）在冻融交替的过渡季节时存在的低估现象。与降轨时刻相比，升轨时刻反演值与实测值间的无偏均方根误差（ubRMSE），偏差（Bias）及均方根误差（RMSE）分别降低了22%、86%及34%，相关系数（R）平均高出4%。

综上，基于升轨数据的反演结果优于降轨数据，升轨时刻基于SCA-V和DCA算法的反演结果皆达到SMAP卫星的任务目标，其中SCA-V的效果最优。这对Zheng等^[24]研究中数据类型选取方面有一定的补充和拓展。一方面，本文研究结果进一步验证了SMAP官方产品选择SCA-V算法作为标准算法的可靠性；另一方面，本研究发现升轨时刻DCA土壤未冻水的反演结果同样理想，精度也满足SMAP的预期目标^[29]，这对加强多极化信息的综合利用有一定启发。

Chan等^[35]人通过验证发现，SCA-V算法在SMAP 2级土壤水分产品（L2_SM_P）中效果最优，并将其作为SMAP的标准算法。这里将基于本文SCA-V得到的未冻水反演结果（SCA-V）与SMAP标准产品（SMAP Product）进行对比。图5列出了二者与实测值的时间序列，由于目前SMAP产品不包含冻结期土壤未冻水的计算，因此图中绿色曲线数据有所空缺。表2分别给出了SMAP产品（降轨：89条数据；升轨：103条数据）及对应日期下本文的SCA-V反演值与实测值之间的误差统计分析。

对于降轨时刻，结合图5（a）和表2发现：在2016年10月与2017年4~5月这两个时间段内（图中以橙色方框标出），SMAP产品对实测值存在明显低估，SMAP产品与实测值的ubRMSE大于0.04 m³m^-3，产品精度未达到预期的任务目标。SCA-V反演值有效改善了SMAP产品对实测土壤水分的低估，且与实测值的ubRMSE，Bias及RMSE平均分别降低了18%、53%及28%，R平均提高了8%。此外，本文SCA-V反演值的ubRMSE达到了SMAP卫星产品的任务目标。

对于升轨时刻，结合图5（b）和表2发现：SMAP产品同样在橙色方框内的两个时段存在明显低估，与实测值的ubRMSE大于0.04 m³m^-3，产品精度未达到任务目标。本文SCA-V反演值同样有效降低了SMAP产品对研究区土壤水分的低估，且与实测值的ubRMSE，Bias及RMSE平均分别降低了26%、76%及39%，R平均提高了8%。此外，本文SCA-V反演值的ubRMSE达到SMAP卫星产品的任务目标。

综上，与SMAP标准产品相比，采用改进的植被和地表粗糙度参数化方案后的SCA-V算法在暖季土壤水分的反演精度更高，基于SMAP升轨和降轨时刻亮温观测的反演值ubRMSE皆达到了卫星任务目标。此外，本文的反演算法因采用四相混合介电常数模型，可以进一步反演冻结期的土壤未冻水，因此更适用于青藏高原地区冻融土壤条件下的水分反演。

Zheng等^[24]的研究分别评估了SMAP默认参数化方案及修改后的植被光学厚度、地表粗糙度和介电常数等参数化方案对亮温模拟及土壤未冻水反演的影响，发现当前SMAP任务采用的微波辐射参数化方案高估了青藏高原地区植被效应，低估了地表粗糙度的影响。本文在此基础上进一步讨论耦合了改进后的各参数化方案对土壤未冻水反演结果的综合影响，发现与SMAP标准产品相比，通过四相混合介电常数模型配置更新后的SMAP算法在暖季反演得到的土壤水分精度更高。此外，本研究与Chan等^[35]的研究结果保持一致，均发现SCA-V算法的效果最优。

本文的研究在SMAP数据类型选取和反演算法方面对Zheng等的研究进行了有效的补充和拓展。 Zheng等^[24]的研究只选择了SMAP卫星早上6点降轨过境时刻的亮温观测数据进行土壤未冻水反演，本研究同时利用了升降轨数据，发现升轨比降轨时刻得到的土壤未冻水反演结果精度更高。此外，该研究只利用了SCA算法进行土壤未冻水反演，本研究发现DCA算法的反演结果同样理想，反演精度达到了SMAP卫星任务要求。

由于采用的Mironov^[18]介电常数模型仅适用于非冻结土壤，当前SMAP产品不包含冻结期土壤未冻水的计算。本研究通过引入考虑冰水共存的四相混合介电常数模型（即式（12）），成功反演了季节冻土冻结和非冻结时期土壤未冻水的动态变化特征。近期，Zheng等^[34]研究发现L波段微波辐射计在冻结和非冻结时期的穿透深度表现一致，冻结和非冻结时期微波辐射特性的差异主要受介电常数变化的影响，因此本文在冻结和非冻结时期采用相同的微波辐射模型和反演算法进行土壤未冻水的反演。

需要注意的是，本研究未对5TM在冻结时期的观测进行标定和验证，但Zheng等^[33]对比评估了冻结时期5TM土壤未冻水观测与Noah陆面模式模拟结果，发现模式模拟与观测的土壤未冻水动态变化一致，从侧面上验证了5TM观测在冻结土壤的有效性。Zheng等^[34]研究发现冻结和非冻结时期L波段微波辐射特性的差异主要受介电常数变化的影响，当前SMAP任务采用的Mironov^[18]介电常数模型不适用于冻结土壤，而适用于冻结土壤的介电常数模型没有综合考虑不同的土壤质地、变化的土壤温度和微波频率等因素^[24]，因此介电常数模型的进一步发展是今后改进土壤未冻水反演精度的关键。

本文主要研究了利用被动微波遥感技术反演负温环境下的土壤未冻水含量，所采用的数据源为SMAP卫星辐射计所记录的亮温观测数据，研究区为青藏高原东部黄河源区的玛曲区域。研究方法为基于零阶辐射传输模型（即τ-ω模型）的SCA和DCA算法，通过四相混合介电常数模型配置更新后的SMAP算法进行研究区土壤未冻水反演。主要结论如下：

（1）基于SMAP卫星不同过境时刻亮温观测及不同算法，对研究区土壤未冻水的反演值均较同步地反映了实测值的动态变化特征（相关系数R均大于0.9），但降轨时刻各算法的反演结果在冻融交替的过渡季节存在明显低估。

（2）无论是用SCA还是DCA反演算法，升轨比降轨时刻得到的反演结果精度更高。这可能是由于上午6点卫星降轨过境时土壤冻融存在异质性，因此不如下午6点卫星升轨过境时的结果。这为将来类似研究的数据源筛选提供了参考。

（3）升轨时刻SCA-V对该研究区土壤未冻水的反演精度最高。无偏均方根误差ubRMSE=0.035 m³m^-3，满足SMAP卫星任务要求，这为将来类似研究的反演算法选择提供了参考。

（4）升轨时刻DCA对该研究区土壤未冻水的反演结果同样理想，ubRMSE=0.039 m³m^-3，同样也达到了SMAP卫星任务要求，这对将来加强多极化信息的综合利用有一定启发。

（5）与SMAP标准产品相比，本文采用改进的植被和地表粗糙度参数化方案以及四相混合介电常数模型的SCA-V算法对研究区未冻水的反演精度更高，更适用青藏高原地区冻融土壤条件下未冻水的反演。由于SMAP的标准算法是对全球尺度而言，且SMAP产品不包含冻结期土壤未冻水的计算，因此必须结合研究区的植被覆盖度和地表粗糙度对参数化方案进行改进，才能得到研究区土壤冻融时期更真实的未冻水含量。

需要注意的是，本文的研究区域仅限于玛曲区域，且观测数据的限制本研究只在一个像元上进行了反演结果验证。今后的工作亟需将发展的算法应用到整个青藏高原，并探讨算法在区域尺度上的适用性和精度，期望能为解析气候变化情景下青藏高原地区的冻土水文变化特征提供科学参考。