基于多时相SAR的IEM和Oh模型地表土壤含水率反演与融合

黄 对,刘 九 夫，张 建 云,王 文,崔 巍

(1.南京水利科学研究院水文水资源与水利工程科学国家重点实验室，江苏 南京 210029；2.河海大学水文水资源与水利工程科学国家重点实验室，江苏 南京 210098)

0 引言

合成孔径雷达(SAR)空间分辨率高、不受云雨影响，其后向散射系数与地表粗糙度、土壤含水率、植被等因素密切相关，故其L、C、X波段影像适用于大范围地表土壤含水率反演[1，2]。基于SAR影像的土壤含水率反演实际是寻找后向散射系数与土壤含水率关系的过程。相关研究采用经验模型[3-6]、理论模型[7-15]和半经验模型[16-19]进行裸露地表土壤含水率反演，一些学者对理论模型进行了改进，如Baghdadi等[20]对地表粗糙度进行定标以获得更真实的积分方程模型(Intergral Equation Model,IEM)模拟数据，后续有学者开展了基于粗糙度定标的IEM[21]、AIEM[22]土壤含水率反演。对于植被覆盖地表，则需利用水云模型或MIMICS模型去除植被的影响，然后采用裸露地表反演模型进行反演[23-27]。

受模型条件、反演方法与雷达成像特征限制，目前研究主要有基于单期影像[6，16]、两期影像[17，23]或多期影像[18，26]单一模型以及基于单期影像多模型[20，22，27，28]的反演方法，反演过程常需依赖实测土壤含水率数据。不同模型的反演结果必然存在差异，而只利用单期或两期影像的反演也具有一定局限性，由于影像特性(极化方式、入射角)不同，即使基于同一模型其反演细节也会有差异，由此导致不同时相反演结果的精度差异与时空连续性问题，有必要探索基于多时相SAR影像的土壤含水率多模型/多反演方法及反演精度评估，并对时空连续性较差的数据进行融合研究，以获取多时相高精度土壤含水率数据。鉴于此，本文首先利用水云模型分析植被因素的影响，分别基于地表粗糙度定标的IEM、Oh模型，在不利用实测土壤含水率的条件下建立多时相SAR影像的土壤含水率反演方法，利用实测数据分析反演精度，进一步对两种模型的反演结果进行加权平均与算术平均融合，探讨融合方法的适用性，为基于多时相SAR的多模型/方法土壤含水率反演以及多时相高精度土壤含水率数据获取提供参考。

1 研究区域及数据

1.1 研究区域

研究区为美国亚利桑那州Walnut Gulch流域(31°40′～31°47′N，109°50′～110°10′W)，属半干旱区，面积约150 km2。流域年均气温17.7 ℃，年均降雨量330 mm，海拔1 200～1 900 m，地表覆盖为稀疏植被，主要为草地与灌丛[29]。该流域是美国农业研究中心(USDA-ARS，http://www.tucson.ars.ag.gov/dap/)和NASA陆地水文研究的长期研究区域，可提供丰富的观测资料，为基于主动微波遥感反演土壤含水率提供了先验知识与验证资料。

1.2 研究数据

(1)土壤含水率数据。研究区共有19个土壤含水率观测站点(图1)，每20 min记录0～5 cm的土壤体积含水率等参数，已经过盐分等误差校正，采样时间为2004年7月1日-9月30日[30]。

图1 研究区与地表覆盖Fig.1 Study area and land cover

(2)SAR影像。采用5期ENVISAT-ASAR(C波段)交叉极化模式成像数据Level 1级产品，为双极化、地距数字影像[31](表1)，下载网址为http://earth.esa.int。利用NEST软件对影像进行预处理：以5×5像元窗口，采用GAMMA MAP滤波进行影像去噪；进行地理编码与辐射定标，将影像DN值转化为雷达后向散射系数值σ0(m2/m2)；结合DEM(SRTM 1sec grid数据，通过NEST软件获取)进行亮度订正[32]，纠正地形导致的辐射与几何畸变，并将σ0(m2/m2)转化成分贝值σ0(dB)，同时将雷达影像重采样为30 m，并利用研究区边界矢量裁切影像。

表1 研究区合成孔径雷达影像信息Table 1 Information of synthetic aperture radar images in the study area

(3)植被含水量数据。植被含水量(mveg)数据由归一化差分红外指数(NDII)和实测植被含水量数据拟合得到(式(1))[33]，NDII由TM影像第4、5波段光谱反射率计算得到(式(2))。研究区包含0805、0818、0824共 3期植被含水量数据，0714期、0720期数据缺失，以0729期的数据代替。

mveg=0.557692308×NDII+0.13923

(1)

NDII=(B4-B5)/(B4+B5)

(2)

2 研究方法

2.1 水云模型

(3)

(4)

γ(θ)2=exp(-2Bmvegsecθ)

(5)

(6)

2.2 IEM模型

基于IEM模型的土壤含水率反演可分为经验[35]和半经验方法[18]，二者均基于模型模拟数据的散射特征分析构建反演公式，区别在于构建时是否利用实测土壤含水率[23]。

(7)

(8)

式中：p为VH极化方式；A、B、C为对应入射角与极化方式下的系数。

(9)

2.3 Oh模型

(10)

(11)

本文基于Lambert′s law假设入射角与单位面积散射量的关系遵循余弦定律[37，38]，通过对邻近期VH极化影像归一化近似获得VH影像，计算公式为：

(12)

2.4 模型精度评价与融合方法

利用均方根误差(RMSE)、偏差(Bias)评价土壤含水率反演精度，计算公式为：

(13)

(14)

式中：Obsi、Simi分别为i样本的观测值和模型反演值；N为样本数。

(15)

式中：i表示期数，本研究选取0720期、0805期、0824期反演结果进行融合。

3 结果与分析

3.1 去除植被影响的效果分析

计算各期影像各地类去除植被影响前后的后向散射系数差值(图2)和土壤含水率反演结果差值(表2)，可以看出，各地类后向散射系数差值多在0.07 dB以下，土壤含水率差值多在0.002 cm3/cm3以下，说明研究区植被影响不大。结合植被覆盖看，研究时段内研究区NDVI<0.23的面积平均占比约98%，以7月19日为例，沿河树林与灌丛、灌丛与稀疏草地以及灌丛与草地3种混合地类的NDVI均值分别为0.192、0.191、0.186，虽然沿河树林与灌丛NDVI最大值为0.24，但该地类像元数量少，且NDVI高值占比少(NDVI>0.2的像元占比约20%)。考虑研究区植被覆盖度整体较低，导致去除植被影响后反演精度并未显著提升。

图2 各地类去除植被前后的后向散射系数差值Fig.2 Difference between backscatter coefficients after and before removing vegetation for various land cover

表2 各地类去除植被前后的土壤含水率反演结果差值Table 2 Difference between soil moisture inversion results after and before removing vegetation for various land cover cm3/cm3

3.2 0～5 cm土壤含水率反演结果

基于去除植被影响后反演的土壤含水率(图3)进行分析，可以看出，IEM模型的反演值主要在0.1 cm3/cm3以下，时空分布较一致，其中0714期与0818期的土壤含水率较低；Oh模型的反演值多为0.1～0.2 cm3/cm3，0720期、0805期和0824期土壤含水率的空间分布总体相似，其中0805期的土壤含水率略高于其他两期，0714期与0818期的空间分布不连贯。从反演方法看，基于IEM模型反演的0714期、0818期地表粗糙度反演公式有别于其他3期，基于Oh模型反演的0714期、0818期VH极化数据由相邻极化数据进行归一化处理近似获得，且土壤含水率的反演公式也有别于其他3期，由此导致即使利用单一反演模型也存在反演方法差异；而基于同一模型同一反演方法的土壤含水率时空分布一致性强(0720期、0805期、0824期)，反演精度差异小。

图3 基于IEM和Oh模型反演的土壤含水率分布Fig.3 Soil moisture retrieved from ASAR images using IEM and Oh model

由同一模型相邻两期的土壤含水率差值(图4)可知，0824-0818期差值主要表现为土壤含水率减少，基于IEM、Oh模型反演的土壤含水率减少区域占比分别约为62%、45%，增加区域占比分别约为37%、35%；基于IEM、Oh模型的0818-0805期差值减少区域占比均约为40%，增加区域占比分别为59%、37%；0805-0720期土壤含水率增加区域占比分别约为44%、61%，0720-0714期土壤含水率减少区域占比分别约为63%、56%。可见IEM和Oh模型反演的相邻两期土壤含水率的主要变化一致。

图4 基于IEM和Oh模型反演的土壤含水率差值Fig.4 Difference of soil moisture (Δmv) between adjacent dates retrieved by IEM and Oh model

从空间分布上看，Oh模型反演异常值主要位于西南和东北区域，IEM模型反演异常值主要位于西南区域。结合流域植被看，土壤含水率异常值区域与流域植被覆盖度相对较高区域并不对应；结合地形因素看，IEM、Oh模型反演异常值区域对应影像后向散射系数极亮值区，这主要源于地形(坡度>15°)导致的成像高亮，虽然对影像进行了地形校正，但在海拔高且坡度大的区域效果有限。

3.3 实测数据验证与融合分析

利用与影像获取时间一致的19个站点的实测土壤含水率数据进行验证(0714期实测数据缺失)(表3)，可以看出，IEM模型反演的土壤含水率RMSE和Bias值均低于Oh模型，R值均高于Oh模型，IEM模型反演的0805期土壤含水率存在一定程度的低估，而Oh模型反演的土壤含水率则普遍存在高估。

表3 基于IEM、Oh模型反演的土壤含水率精度对比Table 3 Comparison of accuracies of soil moisture retrieved by IEM and Oh model cm3/cm3

由于0714期缺实测验证数据，0818期IEM模型反演结果精度高但Oh模型反演精度低，因此只对其他3期结果进行融合。由融合后的统计结果(表4)可知，0720期、0805期和0824期土壤含水率RMSE相比融合前减少了0.003～0.065 cm3/cm3；融合后Bias相对IEM增加、相对Oh减少；融合后R值均有不同程度的提高。其中，加权平均法的RMSE、Bias均小于算术平均法，R值均大于算术平均法，总体上结合实测数据的加权平均法优于算术平均法。以均值融合结果为例(图5)，融合后的土壤含水率在不同时间上的空间连续性较好。受影像期数、间隔时间、特征限制，基于单个模型不同反演方法获取土壤含水率的能力依然有限(如0805期)，融合方法可视为获取更高精度土壤含水率的一种选择。

表4 融合后的土壤含水率精度对比Table 4 Comparison of accuracies of retrieved soil moisture after fusion cm3/cm3

图5 融合后的土壤含水率分布Fig.5 Distribution of retrieved soil moisture after fusion

4 结论

本文基于水云模型分析植被后向散射系数的影响，并分别利用IEM和 Oh模型建立适合不同时相SAR影像的多反演方法，进一步基于两种模型反演结果建立土壤含水率融合方法，结合实测数据分析融合前后的土壤含水率精度，可为获取高精度土壤含水率的多模型/方法反演、融合研究提供参考。主要结论如下：1)研究期内Walnut Gulch流域NDVI<0.23的面积占比达98%，植被覆盖度整体较低，各地类植被去除前后的土壤含水率反演差值多期均值小于0.002 cm3/cm3，去除植被影响对反演结果的作用不明显。2)基于IEM、Oh模型反演的0720期、0805期、0824期土壤含水率的时空分布较一致，0714期、0818期由于地表粗糙度和土壤含水率反演公式不同，导致其时空分布不同；基于同模型同一反演方法的土壤含水率时空分布的一致性强，反演精度差异小，Oh模型在影像数据未满足其设定条件下的反演效果较差。3)实测数据验证表明，基于IEM模型反演的土壤含水率的RMSE与Bias低于Oh模型，其中 0805期的土壤含水率存在低估，而Oh模型普遍存在高估；均值融合、加权融合方法能一定程度克服使用单模型的缺点，融合后土壤含水率在不同时间上的空间连续性较好，融合后的土壤含水率RMSE减小0.003～0.065 cm3/cm3，与实测站点的相关系数均有不同程度提高。本文所用融合方法可视为获取多时相高精度土壤含水率的有效方式。

本文研究区为植被覆盖度低的半干旱区，植被对反演结果的影响并不明显；基于IEM模型的土壤含水率经验反演方法和地表粗糙度等参数的反演方式在其他区域的适用性有待验证；满足Oh模型条件的最佳极化影像方法有待改进；SAR影像在地形陡峭地区的成像异常问题也有待解决。Sentinel SAR作为 ENVISAT-ASAR后续影像，时空分辨率高，下一步将利用该数据在国内高植被覆盖区深入开展多时相SAR的多模型反演方法与融合方法的适用性研究。