基于偏差校正和三重组合分析的主、被动微波土壤湿度数据融合

江 渊,王 文,边增淦

(河海大学 水文水资源与水利工程科学国家重点实验室，江苏 南京 210098)

1 研究背景

土壤湿度是地球生态系统中非常重要的组成部分，与地表-大气的热量交换过程有着非常密切的关系，对陆地表面水分蒸散和循环有很强的调节作用，是许多水文模型、气候模型中重要的输入参数[1]。随着遥感技术的发展，微波遥感现已被广泛用于获取大尺度的地表土壤湿度信息。主动微波遥感的重访周期短且不受云雨影响，但受地表粗糙度和植被覆盖的影响严重[2]；被动微波遥感对近地表土壤湿度的敏感度较高，受到大气和云层的影响小，但受地表温度变化的影响严重，时间分辨率高，空间分辨率低[3]。目前已经有一系列的主、被动微波卫星传感器应用于全球土壤湿度的监测。例如：气象业务卫星计划MetOp系列卫星搭载的ASCAT传感器(2006-)[4]，地球观测系统Aqua卫星搭载的AMSR-E传感器(2002-2011)[5]，地球水文变化观测计划GCOM-W1卫星搭载的AMSR2传感器(2012-)[6]，欧洲空间局土壤湿度和海洋盐度卫星SMOS (2009-)[7]，美国宇航局土壤湿度主动被动卫星SMAP等[8]。主、被微波遥感在不同的植被密度区域反演土壤湿度各有优劣。已有研究表明：被动微波土壤湿度产品在裸土或稀疏植被覆盖下具有较高的反演精度，随着植被密度的增加，土壤湿度的精度会降低[9]；主动微波土壤湿度产品在中等植被覆盖下，尤其在植被变化有明显季节性的地区反演精度较高[10]；在高植被密度下，由于植被冠层使得雷达后向散射或地表微波散射大幅衰减，主、被动微波土壤湿度产品在高植被密度下的精度较低[11]。除了遥感卫星观测外，土壤湿度还可以用陆面模型来模拟。陆面模型通过同化大量的地面实测数据，模拟出的地表土壤湿度在绝对数值上较为合理，但陆面模式模拟土壤湿度数据受模型参数化方案、同化方法、土壤参数、大气强迫、地面观测数据等不确定因素的影响[12]，在缺少观测资料的地区不能很好地模拟出土壤湿度的动态变化。微波遥感的优势则是对土壤水分的变化较为敏感，但微波遥感会受到地表粗糙度和植被覆盖等因素的影响，存在不稳定性[1]。因此以模型模拟土壤湿度数据为基准，对遥感土壤湿度数据进行偏差校正，再将多源遥感数据进行融合，能够得到精度更高、时空连续的遥感土壤湿度融合数据。许多学者已经进行了相关研究：Reichle等[13]以全球陆面数据同化系统(GLDAS)模拟土壤湿度数据为基准，利用累积分布函数(cumulative distribution function,CDF)匹配方法对SMMR卫星土壤湿度数据进行偏差校正，然后将校正后的卫星土壤湿度数据与陆面同化数据进行融合。2009年Liu等[14]以AMSR-E为基准，将AMSR-E和TMI、SMMR、SSM数据做CDF匹配，将匹配后的数据融合得到长时间序列的被动微波土壤湿度融合数据。2011年Liu[15]以GLDAS土壤湿度数据为基准，对AMSR-E和ASCAT进行CDF匹配，并根据AMSR-E反演的植被光学厚度(VOD)，在全球不同的植被覆盖度区域对主、被动微波土壤湿度数据赋以不同的融合权重，并进行加权平均。

三重组合(triple collocation)分析(下文简称为TC分析)可以用来估计大尺度卫星遥感土壤湿度观测的随机误差方差和信噪比，已成为评价土壤湿度产品的重要方法[16]。Scipal等[17]通过TC分析估算ASCAT、AMSR-E土壤湿度数据、ERA-Interim土壤湿度再分析数据各自的误差方差，证实TC分析结果能够客观反映ASCAT和AMSR-E数据的误差特征。Leroux等[18]采用TC分析对SMOS、AMSR-E和ASCAT土壤湿度数据误差进行估计。吴凯等[19]利用ASCAT、AMSR-E土壤湿度数据以及ERA-Interim土壤湿度再分析数据，通过TC分析得到了3种土壤湿度数据的误差方差和信噪比，并结合MODIS土地覆盖类型数据进行分析，发现TC分析能够在草原、农田和裸地得到较为客观的土壤湿度误差方差。Gruber等[20]基于TC分析计算ASCAT和AMSR2数据的误差方差，通过误差方差计算融合权重，并提出基于相关显著性水平检验将像元分类的融合方案，该方法也作为欧洲空间局气候变化倡议(climate change initiative，CCI)项目中土壤湿度数据集ESA-CCI SM3.x版本融合算法的基础。

目前能够在线发布的土壤湿度数据产品均来自ASCAT、AMSR2、SMOS、SMAP、FY3B/C等遥感卫星(或传感器)。从现有对这些数据产品的评估结果来看，ASCAT产品质量最好，FY3B次之，SMOS产品质量较差[21-22]。由于不同遥感观测数据之间存在系统性误差，在数据融合前要先进行偏差校正。常用的偏差校正基准为GLDAS表层土壤湿度(例如ESA-CCI融合产品和NOAA/NESDIS发布的SMOPS融合产品[23])。然而中国陆面同化系统(CLDAS)在模拟中国区域土壤湿度方面比GLDAS具有更高的精度，CLDAS与实测数据相比，具有较高的相关系数，偏差在-0.04～0.04 m3/m3之间，在青藏高原地区也具有较高的精度[24]。因此本文以CLDAS表层土壤湿度产品为基准，对主、被动微波土壤湿度产品进行偏差校正，同时为了便于与ESA-CCI融合产品进行对比，选用与ESA-CCI融合产品近年来一致的数据源，即主动微波土壤湿度产品(ASCAT-A/B)、被动微波土壤湿度产品(AMSR2，SMOS)采用以TC分析为核心的误差分析方法计算主、被动微波土壤湿度数据的权重，进而采用加权平均的方法实现多源微波土壤湿度数据融合，在一定程度上提高了中国陆地区域微波遥感土壤湿度产品的精确度。

2 数据来源与研究方法

2.1 数据来源

本文中使用的数据有微波遥感土壤湿度产品(ASCAT-A/B，AMSR2，SMOS)、校正基准数据，即CLDAS表层土壤湿度数据、验证数据，即土壤墒情站实测土壤体积含水量数据。另外采用ESA-CCI土壤湿度融合产品和ERA-Interim土壤湿度再分析资料进行精度对比分析。时间范围为2017年3月-2018年12月，空间范围为中国陆地区域。

(1)微波遥感土壤湿度产品。ASCAT-A/B L2主动微波土壤湿度产品由欧洲气象卫星应用组织(EUMETSAT)(http://hsaf.meteoam.it/)发布，该产品采用TU Wien water retrieval package (TU WARP)[25]算法反演土壤表层(0～2cm)的土壤湿度饱和度(数值范围0～100%)，空间分辨率为25km；AMSR2 L3被动微波土壤湿度产品是由日本宇宙航空研究开发机构(JAXA)发布[26](https://gportal.jaxa.jp/)，该产品以地表参数反演模型(land parameter retrieval model，LPRM)为反演算法得到地表土壤体积含水量，空间分辨率25 km；SMOS土壤湿度数据产品包括L2和L3产品。SMOS L2产品采用L-band microwave emission of the biosphere (L-MEB)模型反演得到0～5 cm土壤体积含水量，空间分辨率为35～50 km[27]。L3产品是在L2产品的基础上通过加权平均和最优插值得到，其空间分辨率为25 km。本文采用的是SMOS L3产品，由西班牙巴塞罗那专家中心(BEC)发布(http://bec.icm.csic.es/)。

(2)参考数据与验证数据。CLDAS近实时逐日土壤体积含水量分析产品v2.0，由国家气象信息中心发布(http://data.cma.cn)，包含5层土壤体积含水量数据，空间分辨率为0.0625°×0.0625°。该数据主要应用Noah-MP等3种陆面模式，利用大气驱动数据和初始场信息驱动陆面模式得到土壤体积含水量[28]。本文中使用土壤表层(0～5 cm)土壤体积含水量数据。为了与遥感数据匹配，将其每4×4网格求平均得到0.25°×0.25°数据；ESA-CCI土壤湿度融合产品(http://cci.esa.int/)，以GLDAS v2.1陆面模式土壤湿度为基准对两种主动微波土壤湿度产品(ERS和ASCAT)以及5种被动微波遥感土壤湿度产品(SMMR、SSM/I、TMI、WINDSAT及AMSR-E/AMSR2)进行偏差校正，通过融合算法得到1978-2018年6月、空间分辨率为0.25°×0.25°的逐日土壤湿度融合产品[29]；ERA-Interim是欧洲中期天气预报中心(ECMWF)第三代土壤湿度再分析资料(https://apps.ecmwf.int/)。它采用了四维变分分析(4D-Var)实现了再分析资料质量的提升[30]。ERA-Interim土壤湿度再分析资料经过许多研究者检验，在不同季节均可以稳定地重现观测数据变化的空间分布特征[31]。本文中使用了表层0～7cm ERA-Interim土壤湿度再分析数据。

(3)地面实测数据。实测土壤湿度数据来源于云南、河南、浙江、拉萨土壤墒情数据库以及黑河计划数据管理中心。本文使用了36个测站的5～10 cm处土壤湿度。河南、云南、浙江的站点采样时间均为2017-2018年，拉萨采样时间均为2017-2018年4月，黑河流域采样时间为2017年。土壤体积含水量正常值范围为0～50%，将小于0或大于50%的观测值按无效值处理。由于主、被动微波遥感反演的土壤湿度数据为表层土壤湿度，但关于表层的具体深度并没有严格的定义，通常认为微波遥感土壤湿度有效范围为地表0～5 cm[29]，而本文所选的CLDAS数据为地表0～5 cm土壤湿度，ERA-Interim数据为地表0～7 cm土壤湿度，实测数据均选择地表5～10 cm处的土壤体积含水量。由于观测深度的不一致，在偏差校正和验证评估时会带来一定误差，但地表0～10 cm土壤湿度间通常存在密切的关系，垂直变化差异较小[32]，因此本文忽略了由于数据观测深度的不一致所带来的误差。

2.2 研究方法

本文融合多种遥感土壤湿度产品的总体过程是：在基于CLDAS表层土壤湿度数据对多种微波土壤湿度数据进行偏差校正的基础上，对遥感数据与模型数据进行皮尔逊相关分析和TC分析来获取遥感数据的权重，然后对多种遥感土壤湿度产品进行加权平均。所有数据统一采用WGS84地理坐标系统，空间分辨率为0.25°×0.25°，流程如图1所示。

2.2.1 偏差校正 由于主动微波数据表示土壤湿度的饱和度百分比(数值范围0～100%)，被动微波以及模型数据表示土壤体积含水量百分比(数值范围通常为0～50%)。除了数据范围不同，遥感数据也存在着一定的系统性偏差，因此在数据融合前必须将多种来源的遥感数据调整到统一的数值范围，并进行偏差校正。累积分布函数(cumulative distribution function,CDF)匹配方法是一种分段线性匹配技术，可以校正不同数据之间的系统性偏差，已经被广泛地应用于遥感数据融合领域[33]。本文中以经过重采样为0.25°×0.25°的CLDAS表层(0～5 cm)土壤湿度为基准，将ASCAT-A、ASCAT-B、AMSR2、SMOS进行逐网格CDF匹配。

图1 多源微波遥感土壤湿度数据融合流程

2.2.2 三重组合分析 三重组合(triple collocation,TC)分析是估计数据随机误差方差的方法，需要3个误差相互独立的数据。本文中将使用主动微波数据(即ASCAT-A与ASCAT-B的算术平均值，下文中用下标A表示),被动微波数据(即AMSR2、SMOS的算术平均值，下文中用下标P表示)，以及模型模拟数据(即CLDAS表层土壤湿度数据，下文中用下标M表示)组成TC分析的3组数据。TC分析中误差模型表示为[34]：

总的来看，“三师”工作重点解决了大学生在校学习期间学习、生活和思想上中遇到的突出问题。一是解开思想疙瘩，针对一部分文化知识学习基础差的学生，认真做好心理辅导工作，为其鼓劲打气，帮助树立信心。二是解决生活中的实际问题，针对大学生择友、交友、同学间和睦相处、互帮互助等问题，通过谈心谈话发现问题后予以协调解决。三是解开疑惑，及时了解学生思想政治上的各种疑惑，特别是针对反响较为强烈的评优、助困等方面疑惑，严格按照国家相关政策给予解释和答复。

i=αi+βiΘ+εi

(1)

式中：i为土壤湿度的观测值；Θ为土壤湿度真值；αi和βi为观测数据i关于真值Θ的系统误差；εi为数据i均值为0的随机误差。

误差模型有以下4个基本假设[35-36]：(1)观测值与真值之间为线性关系；(2)观测数据的误差具有稳定性；(3)3个数据的随机误差相互独立(不同数据的随机误差之间的协方差为零)；(4)同种数据的随机误差与真值相互独立(随机误差与真值之间的协方差为零)。在TC分析中，3个数据的随机误差方差表示为[36]：

(2)

(3)

主、被动微波数据的权重则采用其均方误差比计算得到，计算公式如下：

(4)

2.2.3 基于相关显著性的融合方案 主、被动微波土壤湿度产品在植被密度不同的地区精确度差异较大，将主、被动微波土壤湿度产品精确度相近的地区进行加权平均，可以增加这些地区的观测次数，减少随机误差；精确度差异较大地区单独使用主动或被动微波数据，可以结合主、被动微波土壤湿度产品优势。根据Gruber等[20]提出的方法，选用模型数据为参考，对主动微波数据和被动微波数据进行逐网格皮尔逊相关分析，判断遥感数据与模型数据以及遥感数据之间是否显著相关(p值是否小于0.05)，是则标记为1，否则标记为0。根据3种数据之间的相关情况，将研究区域划分为5种融合方案：(1)加权平均；(2)算术平均；(3)只使用主动微波数据；(4)只使用被动微波数据；(5)忽略像元。具体分类标准如表1所示。

表1 基于相关显著性的融合方案

注：A为主动微波数据；P为被动微波数据；M为模型模拟数据。

2.2.4 评估指标 本文中对遥感数据的精度分析采用的评估指标包括相关系数(r)、平均偏差(MBE)、均方根误差(RMSE)。各指标的具体计算公式如下：

(5)

(6)

(7)

式中：xi和yi分别为不同的土壤湿度数据；n为统计样本数。

3 结果分析与讨论

3.1 土壤湿度融合结果

对主、被动微波数据进行加权平均得到2017年3月-2018年12月的逐日融合土壤湿度数据。对融合后的逐日数据取平均值得到逐月土壤湿度融合数据并计算中国7大气候分区的面平均值，绘制不同气候区逐月变化过程线。从图2可以看出，7大分区的土壤湿度大致可分为4个梯队。湿润区和热带区土壤湿度最大，全年土壤湿度在0.3 m3/m3左右浮动；东北区和半湿润区土壤湿度数值相近，在0.12～0.28 m3/m3变化且随季节变化最为明显；半干旱区和青藏区位于第3梯队，均在0.1～0.2 m3/m3波动，土壤湿度也呈现季节变化特征；干旱区的土壤湿度最小，全年均在0.11 m3/m3以下且随季节变化幅度很小。

经过融合的遥感土壤湿度数据能够合理的体现出中国陆地土壤湿度的空间分布特征和随季节变化的气候特征，下文将通过多种数据进一步定量分析融合数据的精度。

图2 2017-2018年融合土壤湿度数据不同气候区逐月平均值变化过程线

3.2 融合土壤湿度数据与多种土壤湿度产品的对比分析

3.2.1 融合土壤湿度数据与CLDAS数据对比分析 以CLDAS数据作为基准进行偏差校正和TC分析，为了对比基准数据和融合土壤湿度数据的差异性，将融合数据与CLDAS数据进行对比分析。首先将融合土壤湿度、CLDAS土壤湿度逐日数据求面平均值并绘制逐日变化过程线，见图3。由图3(a)可以看出，融合数据和CLDAS数据的总体变化趋势基本一致，均表现出明显的季节特征，但逐日间的波动有所差别。图3(b)展示了两者之间良好的相关关系，决定系数R2=0.94。为进一步验证融合数据与CLDAS数据的差异，计算融合数据与CLDAS数据之间逐网格的相关系数(r)、平均偏差(MBE)和均方根误差(RMSE)并绘制像元统计直方图，如图4所示。

图3 2017-2018年融合土壤湿度和CLDAS土壤湿度全国面平均值逐日变化过程线和相关性散点图

图4 融合土壤湿度数据与CLDAS的相关系数、平均偏差、均方根误差像元统计直方图

从图4(a)可以看出，大部分地区的相关系数集中在0.3和0.7附近，最大相关系数为0.95，最小为-0.63，平均相关系数0.56；图4(b)表示融合数据与CLDAS数据的平均偏差，由于偏差校正时以CLDAS数据为基准，融合数据与CLDAS数据的平均偏差较小，大部分像元的平均偏差在0附近，平均值仅为0.001 m3/m3。偏差校正消除了遥感土壤湿度产品的系统性误差，再通过加权平均得到的多源遥感融合土壤湿度数据与CLDAS土壤湿度数据的全国平均值随季节变化趋势基本一致，但在日尺度上的小范围波动特征不同。为了进一步验证融合数据的精度，下文将融合数据与其他土壤湿度产品进行对比分析。

3.2.2 融合土壤湿度数据与ESA-CCI融合数据对比分析 本文所用的遥感数据与ESA-CCI融合数据在研究时间内所用数据是一致的，因此将融合土壤湿度数据与ESA-CCI融合数据(数据重叠期2017年3月-2018年6月)进行评估分析，计算融合数据与ESA-CCI融合数据之间的相关系数(r)、平均偏差(MBE)、均方根误差(RMSE)并绘制像元统计直方图，如图5所示。从图5(a)可以看出，相关系数主要分布在0.75左右，最大值为0.98，最小值为-0.97，全国平均值为0.62。由此可以认为融合土壤湿度与ESA-CCI融合数据在全国大部分地区有良好的相关性；由图5(b)可以看出，平均偏差主要集中在-0.1～0.1 m3/m3之间，最大值为0.14 m3/m3，最小为-0.36 m3/m3，全国平均偏差为-0.023 m3/m3，说明融合数据总体略小于ESA-CCI融合数据。融合数据与ESA-CCI融合数据的均方根误差基本接近于0，平均值仅为0.004 m3/m3，说明融合土壤湿度数据与ESA-CCI融合数据之间变化差异较小。

3.2.3 融合土壤湿度数据与ERA-Interim再分析土壤湿度对比分析 为了检验融合数据的精确度，选用不同与CLDAS的模型同化土壤湿度数据进行验证。计算融合数据与ERA-Interim数据之间的逐网格相关系数(r)、平均偏差(MBE)、均方根误差(RMSE)并绘制像元统计直方图，如图6所示。由图6(a)可以看出，正相关像元数量大于负相关像元数量，相关系数最大为0.92，最小为-0.71，平均相关系数为0.42；图6(b)表明大部分地区平均偏差值在-0.2～0.1 m3/m3之间，最大偏差为0.31 m3/m3，最小偏差为-0.40 m3/m3，全国平均偏差为-0.01 m3/m3，说明从全国平均来看融合数据比ERA-Interim数据偏小。图6(c)显示融合数据与ERA-Interim数据均方根误差，最大均方根误差为0.07 m3/m3，平均值为0.012 m3/m3。

3.3 不同产品与实测数据的对比分析

为了进一步验证融合土壤湿度数据的精确度，将融合数据、ESA-CCI融合数据、CLDAS数据与36个站点的实测土壤湿度进行对比分析，表2为每个站点的实测数据与站点所在位置网格的不同产品在相同观测时间内的相关系数(r)，平均偏差(MBE)和均方根误差(RMSE)，并将表中的数据绘制为箱线图，如图7所示。在图7的箱线图中，柱状体表示数据的上、下四分位数范围，柱状体中的横线表示中位数，×表示平均数，两端横线表示数据的上、下限，横线以外的圆点为离群点。

图5 融合土壤湿度数据与ESA-CCI SM融合数据的相关系数、平均偏差、均方根误差像元统计直方图

图6 融合土壤湿度数据与ERA-Interim再分析土壤湿度相关系数、平均偏差、均方根误差像元统计直方图

由图7(a)可以看出，融合数据、CLDAS数据与站点实测土壤湿度都表现出较好的相关关系。除了一个离群点以外，融合数据与实测数据的相关系数为0.30～0.88，平均相关系数大于0.5，相关系数总体上高于ESA-CCI融合数据和CLDAS数据。由图7(b) 可以看出，3种数据与实测数据相比，大部分站点均表现出正偏差，其中融合数据的平均偏差更接近0，数据较为集中，没有出现异常值，上、下四分位数约在-0.01～0.06 m3/m3之间，CLDAS数据和ESA-CCI融合数据的平均偏差较大，出现了异常值。从图7(c)为不同产品与实测数据的均方根误差来看，除了少量异常点外，均方根误差均在0～0.03 m3/m3之间，满足目前全球气候观测系统(GCOS)对遥感观测土壤体积含水量与站点实测值的均方根误差小于0.04 m3/m3的精度要求[29]。融合数据的均方根误差平均值和中位数都更接近于0，且异常值较少，说明融合数据与实测数据的离散程度较小，稳定性更好。综上所述，通过与实测土壤湿度对比可以发现，融合数据比ESA-CCI融合数据和CLDAS数据具有更高的相关系数，较小的平均偏差和均方根误差，总体上精度高于ESA-CCI融合数据和CLDAS数据。

3.4 讨 论

由于主动微波数据与模型数据土壤湿度代表的物理含义不同，即使具有同样的物理含义的数据(如被动微波数据与模型数据)也可能具有不同的数值范围，因此在融合不同产品之前必须进行校正，使不同遥感产品具有相同的物理含义与数值范围。目前最常用的校正方法是选取一个数据基准进行CDF匹配，本文也采用了这种方法。多源微波遥感数据逐日面平均值在偏差校正前后的CDF图如图8所示。

表2 融合数据、ESA-CCI融合数据、CLDAS数据与实测土壤湿度对比分析

图7 融合数据、ESA-CCI融合数据、CLDAS数据与实测数据相关系数、平均偏差、均方根误差箱线图

由图8可看出，对于被动微波数据来说，未校正前AMSR2和SMOS数据整体偏小，尤其AMSR2数据范围分布在0～0.15 m3/m3之间，大幅低于中国区域土壤湿度的实际值。偏差校正消除了遥感数据的整体性偏差，使微波遥感土壤湿度数据与CLDAS数据在时间序列上具有相同的数值范围和相似的累计概率分布。但CDF匹配可能对误差分析结果产生影响，CDF匹配是在时间序列上逐网格地对数据进行分段线性拉升，匹配过程会将模型数据的分布特征加入到遥感数据中[15]，因此可能对TC分析结果造成影响。

TC分析中假设3个数据集的误差互不相关。由于主、被动微波土壤湿度反演模型的反演机理存在很大差异，因此主动微波数据(ASCAT)与被动微波数据(AMSR2、SMOS)之间的误差互不相关。对于模型数据而言，主要误差来自于模型物理机制和驱动数据的误差，因此原始主、被动微波数据和CLDAS数据三者之间误差互不相关，满足TC方法的假设条件。但是CDF匹配会改变原始数据的误差特征，使不同数据之间的误差可能存在相关关系。在这种情况下，TC分析只考虑数据自身的误差方差对分析结果造成的影响。

图8 多种土壤湿度产品的累积分布函数(CDF)图

4 结 论

通过对多源微波遥感土壤湿度数据进行偏差校正、TC误差分析、加权平均可以增加同一像元内的观测频次，减少观测误差。土壤湿度数据融合结合了主、被动微波产品在不同植被覆盖度下的优势，将融合数据与CLDAS数据、ESA-CCI融合数据、ERA-Interim数据、实测土壤湿度进行对比，可以发现虽然多源遥感数据以CLDAS数据为基准进行了偏差校正，使融合数据与CLDAS数据偏差很小，但逐日的土壤湿度波动存在差异。融合数据比ESA-CCI融合数据和ERA-Interim数据略有偏小，与实测数据的相关系数更大，平均偏差和均方根误差更小。本文虽参考ESA-CCI的融合方法，但由于所用的偏差校正基准不同，使得融合土壤湿度数据在中国地区的精确度上优于ESA-CCI融合数据，因此本文的研究可以认为是对ESA-CCI融合数据在中国区域的改进。