北麓河流域不同植被覆盖条件下土壤冻融过程中地表土壤含水量变化

胡越中， 王广军， 杜海波， 梁四海， 罗银飞， 董高峰， 彭红明

（1.中国地质大学（北京）土地科学技术学院，北京 100083；2.内蒙古煤田地质局勘测队，内蒙古 呼和浩特 010010；3.中国地质大学（北京）水资源与环境学院，北京 100083；4.青海省环境地质勘查局，青海 西宁 810007；5.青海省环境地质重点实验室，青海 西宁 810007；6.青海省地质环境保护与灾害防治工程技术研究中心，青海 西宁 810007）

0 引言

地表土壤水分是大气降水、植物含水和地下水相互转化的纽带，是植物生长和生存的决定性要素之一［1］，已成为农业、水文、气候、生态等研究领域的重要基础研究内容［2-4］。世界上三分之一的国家和地区，包括我国青藏高原的绝大部分［5］，处于干旱或半干旱地带，水分的匮乏直接制约了当地的生态环境和经济发展。因此世界各国十分关注土壤水分的研究，并为其投入颇多［6］，如2009年欧空局（ESA）发射的SMOS与2015年美国宇航局（NASA）发射的SMAP均为专注于全球土壤水分监测的卫星。

自20 世纪80 年代，微波遥感技术开始广泛应用于区域土壤含水量监测研究［7-8］。由于土壤介电常数和土壤粗糙度影响土壤的后向散射系数，研究者根据地表电磁波传播理论，先后建立起描述后向散射系数与地表物理与几何参数之间的理论模型，为土壤含水量的反演提供了理论基础，如Kirchhoff Approximation（基尔霍夫近似），其又根据土壤的粗糙程度分为POM（物理光学模型），GOM（几何光学模型）两种［9］、SPM（小扰动）［10］、IEM（积分方程）模型［11］等。同时一些学者基于理论模型，结合陆基后向散射计的实测数据，发展出各种适用于主动微波遥感的经验、半经验模型，如Oh 模型［12］，Dubois 模型［13］等。这类模型大多针对裸露地表而言，除土壤纹理结构外，制约其反演精度的重要因素是地表粗糙度［14-16］，通过正确选择雷达参数或多频率、多极化、多入射角数据，可以在一定程度上降低或削弱地表粗糙度的影响；另外，当地表被植被覆盖时，其反演的土壤含水量往往会偏低，使之无法直接应用于有植被覆盖的复杂地表［17-18］。因此在植被覆盖区反演土壤含水量往往需要借助光学遥感提供的参数去除植被层对后向散射系数的影响，如MIMICS（植被微波散射模型）模型［19］，水云模型等［20］。

但总体而言，相较于传统的点尺度物理观测、面尺度的光学或被动微波遥感数据土壤含水量反演方法，主动微波遥感技术，特别是合成孔径雷达（SAR）技术，以其宏观、高效、全天候、高分辨率和穿透性强等特点已经被广泛地应用于农业、水文和生态等领域［21-25］。

北麓河流域位于青藏高原腹地多年冻土区，陆面生态系统以高寒植被为主。已有的研究结果表明，青藏高原土壤冻融对土壤水热变化影响较大［26］，冻结过程有利于土壤维持其水分，而土壤含水量反过来又影响土壤冻融过程及热量分布［27］，进而影响植被的发育程度［28］。在北麓河左冒西孔曲小流域内的定点观测结果表明，土壤冻融过程、土壤含水量和植被生长呈现出复杂的关系［27］。多年冻土活动层开始冻结和消融时间随着植被盖度的减少不断提前，且随着植被盖度减小，活动层含水量变化速率增大，植被起到抑制土壤含水量变化速率的作用［29］。李元寿等［30］的研究结果进一步表明，北麓河左冒西孔曲小流域内多年冻土活动层浅层土壤冻结和融化对植被覆盖度的响应程度较强，接近深层土壤冻结和融化对植被覆盖度的响应程度降低，且伴随草地退化、植被覆盖度降低，土壤含水量对植被覆盖响应越加强烈，植被覆盖度越高，土壤含水量响应越是滞后且土壤含水量变化越加平缓。Wang 等［31］在忽略植被作用的前提下，利用2014—2016 年TerraSAR-X 多角度时间序列数据去除地表粗糙度的影响，基于双组分土壤含水量检索模型完成了北麓河流域内小范围（＜100 km2）土壤含水量反演，但文中没有对年内土壤含水量变化情况进行进一步分析。北麓河多年冻土区活动层冻融过程中，植被覆盖和土壤含水量的响应仍有待于进一步深入研究。

2014年4月3日，ESA 哥白尼计划（Global Monitoring for Environment and Security）中的地球观测卫星Sentinel-1A 卫星发射升空，并于2017 年4 月对外共享数据，为主动微波反演土壤含水量提供了新的数据支持。针对Sentinel-1A 数据源，El Hajj 等［32］结合IEM 与水云模型建立带有噪声干扰的模拟C波段SAR 数据库，并利用其进行神经网络的训练与验证。文中尝试多种神经网络的输入数据组合，充分研究了土壤含水量与C 波段SAR 数据和植被指数的相关性。结果表明单独使用VV 极化或协同使用VV 与VH 极化反演土壤含水量均能取得较好的效果；Bao 等［33］基于水云模型发展出一种土壤含水量与雷达后向散射系数和植被指数关系的半经验模型，并通过在英国与西班牙地区的实测数据进行拟合来获得模型所需的参数。结果表明该模型可应用于复杂的植被覆盖地表。

本文结合主动微波遥感与光学遥感两种数据源，在改进的“水云模型”去除植被覆盖影响的前提下，协同利用Dubois 模型与Oh 模型以消除难以获取的土壤参数，利用模型正向模拟建立数据库并通过代价函数最小化的方式反演土壤含水量。通过分析时序反演结果，从面尺度研究北麓河流域2018年年内2 037.94 km2土壤冻融过程、土壤含水量和植被覆盖的关系。

1 研究区域及数据

1.1 研究区域概况

北麓河流域位于青藏高原腹地长江源区，地理位置如图1 所示，坐标范围为：92°54′～92°56′ E，34°48′～35°00′ N。研究区海拔为4 510～4 680 m，多年冻土极为发育，活动层厚度一般为2～3 m，地表土层以细粒土（黏土、亚砂土）为主，局部分布少量粗颗粒土（角砾土、碎石土）［34］，属青藏高原可可西里区冲、洪积高平原地貌［35］。

图1 北麓河流域位置示意Fig.1 Overview of Beiluhe basin

研究区植被生长发育情况较好，主要的植被类型为高寒草甸，主要植物种为矮嵩草、短穗兔儿草、藏嵩草等。在高寒荒漠区间或分布着火绒草、垫状点地梅等。该地区属于亚寒带半干旱气候，年均降水量在290～300 mm 之间，集中在6—8 月，占全年80% 以上；年最高气温为19.2 ℃，最低气温为-27.9 ℃，年均气温为-5.0～-3.8 ℃［35］，土壤冻结期为9月—次年4月［34］。

1.2 数据源和预处理

Sentinel-1A卫星处于高度为693 km的太阳同步轨道，轨道倾角为98.18°，轨道周期约为96分钟，重访周期为12天。Sentinel-1A卫星载有C波段合成孔径雷达（SAR），使用TOPSAR的扫描方式，可较好的解决ScanSAR 幅度不匀调制的问题，保证了扫描区内的影像质量均匀［36］。其中心频率为5.405 GHz，极化方式为VV+VH 与HH+HV，入射角度为20°～45°，辐射分辨率为1 dB·（3σ）-1，最大噪声等效散射系数（Noise Equivalent Sigma Zero）为-22 dB。包括4 种工作模式：条带模式（Strip Map Mode，SM），超宽幅模式（Extra Wide Swath，EW），宽幅干涉模式（Interferometric Wide Swath，IW）和波模式（WaveMode，WV）［37］。

本文使用的SAR 数据为IW 模式的Level-1 级地距多视GRD 产品。极化方式为VV+VH，空间分辨率为5 m×20 m，幅宽为240 km，重访周期12 天，成像时间为2018年1月3日—2018年12月29日，共29 景。为抑制SAR 影像斑点噪声，利用GammaMAP［38］滤波对影像进行去噪处理，滤波窗口大小为3×3；再使用SARspace 软件进行辐射定标，地理编码，重采样与裁剪，最终获得研究区空间分辨率为10 m的后向散射系数影像。

此外，采用Landsat-8卫星上搭载的陆地成像仪（OLI）数据反演植被覆盖度，以消除植被冠层散射对土壤含水量反演结果的影响。Landsat-8 OLI 数据包括9 个波段，波长范围为0.43～12.51 μm，重访周期为16天（成像时间与Sentinel-1A最近影像最多相差8 天），空间分辨率为30 m，包括一个15 m 的全色波段。所使用的L1TP 产品已经经过系统辐射校正、几何校正和地形校正，因此只需要进行辐射定标，将原始DN 值转换为大气外层表面反射率，之后利用ENVI软件提供的FLAASH 大气校正模块消除大气影响，并重采样为10 m 分辨率的影像与Sentinel-1A 影像对应。使用Landsat-8 红光（band4）和近红外（band5）波段计算得到归一化植被指数（NDVI）并应用像元二分模型反演植被覆盖度［39］。此外，考虑到没有实测叶面积指数（LAI）数据，无法由NDVI通过线性回归求解LAI，以及采用其它经验模型由于应用区域不同会引入新的误差等问题，本文在前述植被覆盖度反演的基础上，借鉴陈丽等［40］的研究方法，采用多次散射过程冠层模型求解迭代的方法来逐步逼近精确的LAI 数值，以作为后续水云模型的输入参数之一。

为获取研究区表层土壤含水量数据，于2019年8月27日对北麓河流域进行了星地近似同步实地勘测（图1），当日天气阴，气温约为4～19 ℃，表层土壤未冻结。Sentinel-1A 卫星过境时间为北京时间2019 年8 月27 日6 时52 分，实地采样时间为当日08：00—16：00。考虑到地形与交通条件，沿109 国道选择地势平坦，植被长势均匀的地方布设土壤含水量采样点（图1）。采用土壤含水量记录仪（Uni1000土壤含水量记录仪）获取了地表0～8 cm深度的土壤含水量实测数据。采集样方4个角点和中心点共5 个土壤含水量数据，通过距离加权平均计算得到与临近影像像元地理位置相同的土壤含水量。依照上述方法，共采集16个样方的土壤含水量均值数据。

2 研究方法

2.1 改进的水云模型

水云模型是一种适用于低矮植被的半经验后向散射模型，已被证明在C、X 波段和20°～40°入射角的范围内模拟精度较好［41］。本文在考虑单次散射的前提下，将给定像元（像元后向散射系数为植被和裸土混合，即混合像元）的植被覆盖度引入水云模型以进行改进，使混合像元内总后向散射系数分解为植被覆盖区地表贡献和无植被覆盖地表裸土地表贡献，则雷达后向散射分量可以写成：

式中：为给定混合像元总后向散射系数，上标pp为极化方式（本文为同向极化VV 或交叉极化VH）；f为植被覆盖度。为更好描述冠层的散射特性，Prevot等［42］建议使用叶面积指数LAI代替植被含水量，则为植被冠层（水云层）后向散射系数，τ2为双向植被透射率或衰减因子，分别由下式［42-43］给出：

式中：a，b为经验系数，因研究区内土壤种类较多，难以获取空间分辨率足够高的土壤类型数据，因此采用Bindlish 等［44］的研究成果，分别取a，b为allland类型0.0012和0.091；θ为雷达入射角。

则根据式1，裸土地表或植被下垫面后向散射系数可由下式得到：

依照式（2）、（3）、（4），利用Sentinel-1A 影像处理得到的VV，VH 极化后向散射系数值以及与其时间最接近的Landsat-8 影像处理得到的植被覆盖度和叶面积指数，计算得到土壤表层的后向散射系数值。由此，土壤含水量的反演转化为去除植被覆盖影响后土壤表层后向散射系数刻画的问题，可以使用针对裸露地表反演土壤含水量的模型代替水云模型中的土壤表层后向散射系数函数。目前常用的经验、半经验模型有Dubois 模型，Oh 模型，Zribi模型等［45］。本文在前述考虑单次散射和混合像元求得的裸土地表或植被下垫面后向散射系数的基础上，综合Dubois 模型和Oh 模型，进行土壤含水量的反演。

2.2 模型的综合使用

Dubois 模型是基于地面实测数据建立的经验模型。在1.5～11 GHz 频率范围，30°～65°入射角范围内，Dubois模型可以较好地模拟VV与HH后向散射。但是在地表较为粗糙［ks＞2.5，k为自由空间波束，s为均方根高度（cm）］，土壤水分较高时（ms＞35%，ms为土壤体积含水量）Dubois 模型将不再适用。由Dubois经验模型［13］：

式中：ε为土壤介电常数；λ为波长（cm）。则可推导出：

Oh 提出了同向极化和交叉极化的后向散射系数经验模型［12］。模型的适用的条件为：4%＜ms＜29.1%，0.13＜ks＜6.98，10°≤θ≤70°。

将式（6）代入Oh模型：

则可去除地表粗糙度的影响，式（7）中，ms为土壤体积含水量。进一步结合TOPP公式［46］：

可以建立起土壤体积含水量ms和土壤同向极化VV和交叉极化VH的关系模型，因公式过于复杂，为简化起见，本文表示为：

本文利用查找表法［47］解决式（9）ms反演问题。对于某一像元，利用式（4）去除植被影响后，根据，ms与的关系模型对进行正向模拟（为真实值，ms范围设定为0.01～0.99 cm3·cm-3），从而建立起查找表。并利用与的平方误差代价函数S查找结果，即使S最小化时相对应的ms最优解作为最终的反演结果。S表达如下：

式中：为观测VH 极化后向散射系数；为使用关系模型模拟的VH极化后向散射系数。

2.3 反演流程

本文反演土壤含水量的流程如图2所示：

图2 反演流程Fig.2 The inversion process

（1）首先选择Sentinel-1A 影像及与其时间最为接近的Landsat-8 影像，并对其进行预处理。利用Landsat-8 影像计算得到研究区的植被覆盖度及叶面积指数。

（2）逐像元进行判断，若其植被覆盖度大于0，则利用改进水云模型进行处理，最终获得VV 与VH极化土壤表层的后向散射系数。

（3）综合Dubois、Oh、Topp模型，消除中间量ks，建立ms、和三者间的关系式。并利用查找表法对每一像元进行反演，最终获得研究区土壤含水量影像。

3 结果与分析

3.1 精度验证

图3为实测土壤含水量与反演土壤含水量的散点图：实线为回归方程，虚线为1∶1 直线，反演土壤水含量误差棒为1.96 倍标准差（95%置信区间），实测土壤含水量误差棒为仪器误差；回归方程调整决定系数（AdjustedR-square，Adjusted-R2）为0.6848，均方根误差（Root Mean Square Error，RMSE）为0.039 cm3·cm-3。当土壤含水量较高时，会出现明显低估的现象，从而造成拟合曲线斜率小于1。Dubois模型与Oh模型的最佳适用范围分别为：土壤含水量在0～0.35 cm3·cm-3与0.04～0.291 cm3·cm-3以内。因此可归因于土壤含水量超出最佳适用范围，从而造成的低估现象。此外，在本文中多种模型运用过程中的累积误差也是造成实测土壤含水量与反演土壤含水量存在差异的原因。

图3 实测与反演土壤含水量散点图Fig.3 Scatter plot of soil moisture content from measurements and inversions

3.2 土壤含水量年度变化分析

整体上看，北麓河流域1 月上旬（1 月3 日）—3月下旬（3 月28 日）土壤含水量基本无变化；从4 月初（4 月9 日）开始缓慢增加，5 月初（5 月3 日）—6 月中旬（6 月20 日），土壤含水量开始迅速增加并达到最高；6 月中旬后（6 月20 日），除个别地区略有波动外，北麓河流域土壤含水量在在整个夏季保持较稳定的状态；从9 月上旬（9 月12 日）开始，土壤含水量开始下降，经过两个月的下降期至11 月中上旬（11月11日）达到稳定，土壤含水量不再有明显的变化，并一直持续到年末（12月29日）（图4）。

图4 2018时序土壤含水量分布Fig.4 Time series of soil moisture content in 2018

北麓河地区5 月中旬开始进入夏季融化期［26］，且6 月中旬以后，北麓河地区的降水类型主要以降雨为主［30］，因此，降雨和升温导致的活动层融化是6月土壤含水量激增的重要原因；夏季土壤含水量的波动变化则主要是由于降雨量变化造成的；而9 月—11 月的土壤含水量下降则是由于降水减少且主要为固态降水造成。

但需要注意的是，北麓河流域部分地区的土壤含水量时序变化不尽相同，这可能与植被覆盖或土壤质地有关，因此有必要根据植被覆盖和土壤质地等选取典型区进行对比分析。

为弱化地形对雷达后向散射系数的影响，且考虑到研究区植被覆盖普遍较低（夏季最高植被盖度普遍在30%以下，个别水热及土壤养分较好区域夏季最高植被盖度在90%左右）选取地形平坦、植被覆盖度能进行显著区分及土壤质地差异明显的4个区域作为典型研究区（图1），分析研究区土壤冻融过程、土壤含水量和植被覆盖的关系。为减弱土壤含水量空间变异性的影响，每个区域大小为3×3像元。

图5为各典型区的年内土壤含水量变化。对比各典型区土壤含水量变化折线，可以发现：

图5 典型区土壤含水量变化折线图Fig.5 Temporal variations in soil moisture content in the four typical areas

（1）相较于其它典型区，90%植被盖度区多年冻土活动层在11 月—4 月期间含水量基本无变化，保持稳定，但随着植被覆盖度的降低，在此期间，其它典型区（沙化地表除外）的土壤含水量则波动较大。

分析原因，这可能与地表均质性或土壤质地对雷达后向散射系数的影响有关。90%植被盖度区（主要是高寒草甸）土壤地表土质为亚砂土、亚黏土，且地表腐殖质较厚［5］，地表均质性较好，进而减少了对雷达后向散射系数的干扰，且在此期间，降水主要以固态水形式，活动层冻结后未冻水分很小。

活动层与大气的水汽交换主要在土壤表层以升华等形式进行，水汽交换量很小［26］，而植被层又具有保温作用［25］，从而使高植被覆盖区土壤含水量保持较好的稳定性。但低植被覆盖区砾石较多且分布不均，或裸土地表景观破碎化较大（图6），而对雷达后向散射系数影响较大，进而使反演的多年冻土活动层土壤含水量波动性较大。沙化地表由于主要是细砂、粉砂组成，均质性较好，故沙化地表土壤含水量在此期间也保持了较好的稳定性。

图6 低植被覆盖区地表实地照片Fig.6 Photos of the low vegetation cover areas

90%植被盖度的典型区土壤含水量从4 月9 日开始稳步缓慢增加至5 月3 日，推测该时间段内土壤升温较慢［26］。考虑到卫星重访周期，该区多年冻土活动层消融时间应在4 月9 日之后，这与沱沱河埋深4 cm 活动层消融时间为4 月15 日的野外观测结果［27］吻合（沱沱河距北麓河直线距离＜90 km），进一步验证了本文观测结果的准确性。同理，根据雷达反演的土壤含水量年内曲线，推测该区冻结期应在10月18日后。

（2）虽然30%植被盖度的低植被覆盖区11月—3 月期间波动较大，但从图4 中仍然可以看出，3 月28 日—次年4 月21 日，土壤含水量有明显增加趋势。同因卫星重访周期，推测30%植被盖度区多年冻土活动层消融时间应在3 月28 日后。同时根据土壤含水量变化曲线，该区的冻结期应早于90%植被盖度的典型区。

（3）裸土地表从3 月4 日开始，土壤含水量开始呈现波动性上升趋势，因此推测裸土地表活动层的开始融化时间应在3 月4 日后。对比后续土壤含水量推测在3 月16 日后研究区有一次降温过程，致使3 月28 日土壤含水量再次降低。而同期90%植被盖度典型区含水量无变化，可能是由于植被保温效果，导致气温波动对土壤含水量的影响较少所致。

（4）沙化地表在6 月8 日后土壤含水量才开始出现明显攀升，10 月18 日后下降到低值点，推测埋深在8 cm 以上的沙化地表受冻融作用影响很小，其含水量主要受降水影响。

（5）总体上看，植被覆盖对活动层冻融过程具有明显的迟滞作用，即植被覆盖度越高，活动层冻结和消融时间越滞后，这与刘光生等［29］和李元寿等［30］的研究结论一致。

4 结论与展望

本文基于2018 年Sentinel-1A 时序数据，辅以Landsat-8 光学遥感数据，利用改进的“水-云”模型，综合Dubois 模型、Oh 模型和TOPP 公式，完成了北麓河流域2018年土壤含水量的反演，并基于实测数据验证精度。通过对年内土壤含水量时序反演结果的分析，得到了土壤冻融过程、土壤含水量和植被覆盖的关系。主要结论如下：

（1）本文提出的在改进的“水云”模型基础上，综合Dubois 模型、Oh 模型和TOPP 公式进行土壤含水量反演，有效去除了植被覆盖和地表粗糙度的影响。16 个点的同步实测结果表明：回归方程为Adjusted-R2为0.6848，RMSE为0.039 cm3·cm-3。

（2）地表均质性或土壤质地对雷达后向散射系数的影响较大。11 月—次年4 月期间，因高植被覆盖度典型区地表均质性较好，加之降水主要以固态水形式，活动层与大气的水汽交换量小，以及植被的保温作用，含水量基本无变化，保持稳定，但随着植被覆盖度的降低，在此期间，其它典型区（沙化地表除外）的土壤含水量则波动较大。

（3）植被覆盖对活动层冻融过程具有明显的迟滞作用，即植被覆盖度越高，活动层冻结和消融时间越滞后。推测裸地、30%植被盖度和90%植被覆度的典型区2018 年消融开始时间分别在3 月4 日后、3月28日后和4月9日后。90%植被盖度典型区冻结期应在10 月18 日后，而裸地和30%植被盖度典型区则早于90%植被盖度典型区。

（4）本文研究结果中北麓河地区年内土壤含水量变化趋势与前人定点观测所得到的趋势基本吻合，表明基于Sentinel-1A 年度时序数据，从流域面尺度研究土壤冻融过程、土壤含水量和植被覆盖的关系是可行的。

多年冻土活动层冻融期土壤水分以未冻水与冰两种形式存在，散射机理较为复杂；同时因研究区冻融气候环境恶劣，未能采集该时期的研究区土壤含水量数据，因此冻融土壤含水量反演的方法与精度有待进一步研究和验证。另外，10 月30 日Sentinel-1A数据缺失，导致本文无法给出上述典型区的准确冻结日期。

土壤含水量与降水等密切相关，而北麓河流域又地处青藏高原腹地，气候多变，多数情况下研究区内各地夏季降雨不同步，要准确地探究北麓河流域多年冻土区，特别是夏季土壤含水量与植被覆盖的关系，仍是后续要持续探究解决的难题之一。