陈鲁皖,韩 玲,王文娟,秦小宝
(长安大学 地质工程与测绘学院,陕西 西安 710064)
土壤水分,是全球气候系统的核心变量之一[1]。对于地球表层的能量和水分平衡具有重要作用,土壤水分的时空分布与动态变化对陆地-大气间热量平衡、大气环流产生显著影响[2]。土壤水分的空间分布对于干旱检测、农业生产和田间灌溉管理也具有极为重要的意义,尤其是在干旱半干旱的中国西北地区。对土壤水分进行动态、实时、大范围和较高精度的估算,利用空间分辨率较高的微波遥感反演是较为有效的手段。
土壤水分反演过程中的不确定性分析是当前研究的热点,E. De Keysera等[3]提出了一种基于线性回归的粗糙度参数化方法,并使用概率分布对其不确定性进行估计,通过反演模型进一步的传播不确定性,得到反演土壤水分的概率分布。J. FernA′ ndez-GA′ lvez等[4]认为大多数土壤水分反演方法依赖于土壤水分与介电特性之间的关系,比较了包含这种关系的9种反演模型,并对土壤水分和介电特性之间的不确定性引起的误差进行量化。Ma C等[5]提出了一种基于贝叶斯定理和马尔可夫链蒙特卡罗方法的土壤水分概率反演(PI)算法,能够定量描述土壤水分反演的不确定性,通过最大似然估计得到高精度的土壤水分估计。李大治[6]通过集合反演来研究土壤水分反演过程中由观测误差、模型参数误差以及不同的反演策略所引起的不确定性。verhoest等[7]研究发现,利用SAR进行土壤水分反演,地表粗糙度参数的不确定性是导致土壤水分反演误差的重要来源。现有的不确定性研究大都集中在如何定量描述模型参数的不确定性,以及研究参数的不确定性对土壤水分反演结果不确定性的影响方面,而对于利用模型参数的不确定性来控制土壤水分反演结果的精度方面的研究较少。
本文首先对AIEM模型的主要参数进行敏感性分析,然后使用Rahman的LUT(look-up tables)方法[8]来反演土壤水分,研究AIEM模型主要影响参数的不确定性,然后基于参数不确定性计算得到不同入射角时满足一定反演精度的主要影响参数的误差控制范围。在反演中使用LUT法计算有效相关长度(lmod)来代替实际测量值,可避免引入粗糙度参数测量过程中的不确定性[9]。根据不同程度的高斯噪声对主要影响参数的取值进行随机扰动,得到大量各参数采样值,输入AIEM模型中得到带噪声的后向散射系数集合;通过LUT方法反演得到Mv,计算得到反演结果的RMSE和标准方差,得到满足一定反演精度的各影响参数的误差量级控制范围;研究得到不同入射角时满足一定反演精度的各影响参数的误差量级控制范围。通过使用ENVISAT ASAR影像和同步地面观测数据试验验证,证明本文得到的各影响参数的误差量级控制范围在较大的入射角范围内都适用。
地面试验数据来源于中国科学院寒区旱区科学数据中心所提供的“黑河综合遥感联合试验”。研究样区位于甘肃省黑河中游的临泽样地,样地较为平坦。临泽地区气候类型为大陆性荒漠草原气候,年蒸发量1 830.4 mm,平均降水量118.4 mm,年平均气温7.7 ℃。选取了与2008-05-24、2008-07-11 ASAR影像同期的临泽采样区地面土壤水分观测数据。
选用欧空局的ENVISAT-1卫星上ASAR传感器获取的SAR影像作为土壤含水量反演的数据源。采用了包含临泽地区的2008-05-24、2008-07-11一共2幅ASAR数据,入射角θ范围依次分别为(15.0°,22.9°)、(31.0°,36.3°),经度范围(99°28′E,101°24′E),纬度范围(38°42′N,40°03′N),入射波段为C波段(f=5.331 GHZ),地面分辨率为12.5 m×12.5 m,工作模式为Alternating Polarization,极化方式为VV和VH两种见图1。其中图1(a)、图1(b)分别为研究区的TM5、4、3假彩色合成影像和ASAR的VV极化影像(2008-07-11)。
为了实现敏感性和不确定性的定量分析,通过几个步骤进行研究:
1)对AIEM模型的各输入参数进行全局敏感性分析,得到各主要参数对模型输出后向散射系数的影响,确定各参数在模型中的贡献大小;
2)对临泽样区的实测数据的取值范围进行统计分析,得到AIEM模型输入参数标定值和主要影响参数的取值范围;
3)根据AIEM模型各主要影响参数的标定值和不同标准方差产生不同量级的高斯噪声,对各参数的取值进行随机扰动,每个量级得到1 000个带噪声的采样值,输入AIEM模型中得到带噪声的后向散射系数集合;
4)通过LUT(look-up tables)方法反演得到Mv,计算反演结果的RMSE和标准方差,研究土壤水分反演RMSE对各影响参数的响应特征,满足一定反演精度的各影响参数的误差量级控制范围;
5)不同入射角时满足一定反演精度的各影响参数的误差量级控制范围。
为了实现敏感性和不确定性的定量分析:
1)对AIEM模型的各输入参数进行全局敏感性分析,得到各主要参数对模型输出后向散射系数的影响,确定各参数在模型中的贡献大小;
2)对临泽样区的实测数据的取值范围进行统计分析,得到AIEM模型输入参数标定值和主要影响参数的取值范围;
3)根据AIEM模型各主要影响参数的标定值和不同标准方差产生不同量级的高斯噪声,对各参数的取值进行随机扰动,每个量级得到1000个带噪声的采样值,输入AIEM模型中得到带噪声的后向散射系数集合;
4)通过LUT(look-up tables)方法反演得到Mv,计算得到反演结果的RMSE和标准方差,研究土壤水分反演RMSE对各影响参数的响应特征,得到满足一定反演精度的各影响参数的误差量级控制范围;
5)继续研究不同入射角时满足一定反演精度的各影响参数的误差量级控制范围。
AIEM模型是计算随机粗糙面单次电磁散射问题的近似方法[10],是目前最广泛使用的裸露土壤表面散射模型。它从电场和磁场所满足的积分方程出发,得到了表面上切向电场和磁场的表达式。
AIEM模型单次散射表达式为[11-12]
(1)
(2)
AIEM模型适用于裸土地区的土壤水分反演,而对于农田及低矮植被覆盖的地区,AIEM模型也适用。Lievens[14]等人在2011年使用水云模型进行了大量试验,研究发现低矮植被尤其是谷物冠层造成散射衰减在很大程度上可以由直接的冠层贡献来补偿,这导致低矮植被覆盖区与裸土相比较雷达散射系数相差很小。因此,本研究未考虑作物覆盖对后向散射系数的影响。
AIEM模型的输入参数包括传感器参数和地表参数。其中传感器参数为入射角θ,地表参数包括Mv,S和l,地表土壤有效温度st和土壤质地(砂土比例sv和粘土比例cv)。
由于AIEM模型比较复杂,输入参数较多,通过对AIEM模型进行参数敏感性分析,得到影响后向散射系数输出的各参数的贡献比例,确定主要影响参数,这样在不确定性分析中就可以对主要影响参数的不确定性进行研究。
本研究中使用扩展傅里叶振幅敏感性分析(Extend Fourier Amplitude Sensitivity Testiing,EFAST)方法来定量分析AIEM模型中各主要输入参数对后向散射系数的影响[15]。EFAST方法的基本思想来自于贝叶斯定理[16],即认为模型输出结果的敏感性可以用模型结果的方差来反映。即参数x的敏感性可以表示为[17]
(3)
式中:Y为模型的输入值;x为输入参数;E(Y|X)为当x取某一定值时Y的期望,varx是x遍历x的取值范围时的方差。EFAST方法可通过傅里叶变换求得模型中参数i忽略与其他参数耦合以后的一阶敏感性指数:
(4)
对于某一参数i在考虑参数耦合作用后的总敏感性指数为
(5)
式中:Si为参数i一阶敏感性指数;STi为参数i的总敏感性指数;E(Y|X-i)为E(Y|Xi)的补集。
本文使用Rahman的LUT方法来计算有效相关长度。Lievens[14]等试验了大量不同的(S,l)组合,发现固定住S的取值,通过微波散射物理模型反演得到校准的l,这样得到的(S,l)组合反演Mv效果较好,并发现S=1 cm,l∈(1 cm,120 cm)时,反演精度较高。
对AIEM模型的主要输入参数的取值进行统一标定。为尽量贴近试验实测数据,地表土壤有效温度取实测平均土壤温度的均值st=22 ℃;土壤质地参数取值根据“黑河流域HWSD土壤质地数据集”中临泽地区砂土和粘土比例的平均值来确定,sv=24%,cv=32%;雷达入射角为θ=30°;有效粗糙度参数S=1.4 cm,l=20 cm;土壤水分Mv=30%,见表1。
表1 AIEM模型参数的标定值
通过敏感性分析,AIEM模型的主要影响参数后,为了定量描述模型主要影响参数的不确定性,根据不同程度的高斯噪声对各参数值进行随机扰动,即以表1中标定值为期望,标准差从小到大设置8个不同量级。每个量级生成1 000个加入高斯噪声的参数集合。将该参数集合输入AIEM模型得到后向散射系数的模拟集合,然后利用LUT方法反演土壤水分Mv,计算土壤水分的RMSE和方差,并通过响应曲线研究各主要影响因子对反演土壤水分的影响。
使用基于AIEM模型的LUT[22]法来反演土壤水分,LUT法的核心思想是建立一个简单高效的查找表(look-up tables,LUT)。首先利用AIEM模型获取不同粗糙度和土壤水分下的后向散射系数集合,S=1.4 cm,l∈(1 cm,120 cm),步长为1 cm,Mv∈(1%,60%),步长为1%,将模拟得到的后向散射数据存储于表格中为LUT表。
然后利用成本函数在LUT表中查找符合条件的后向散射系数对应的土壤水分值。成本函数为
(6)
利用开源Matlab工具Simlab软件的FAST模块来进行敏感性分析,Simlab是由Joint Research Centre of the European Commission基于蒙特卡洛方法设计的非商业软件。使用Matlab编程实现AIEM模型主要输入参数的敏感性分析,在入射频率给定的情况下(C波段,入射频率),AIEM模型的输入参数包括入射角θ、土壤水分Mv、均方根高度S、相关长度l、地表土壤有效温度st和土壤质地(砂土比例sv和粘土比例cv)共7个,假设所有参数都为均匀分布。表2为AIEM模型各参数值域范围,各参数在值域空间内的采样次数为5 000。
表3显示了AIEM模型参数的总敏感性指数STi,STi贡献比例和一阶敏感性指数Si。
表3 AIEM模型参数的敏感性指数
由表3可知,AIEM模型各输入参数中,S,l,Mv和θ这4个参数的贡献比例相加约96%,一起主导了AIEM模型的后向散射系数输出。
根据敏感性分析的结果,选取均方根高度S、相关长度l和雷达入射角θ这3个敏感性参数进行不确定性分析。为了定量描述土壤水分反演过程中S,l和θ这3个参数的不确定性,以表1中的各参数的标定值为期望,不同参数使用不同标准差形成噪声量级:S(0.01,0.02,0.03,0.04,0.05,0.07,0.09,0.1),l(0.2,0.4,0.6,0.8,1,2,3,4),θ(0.2,0.4,0.6,0.8,1,2,3,4),每个噪声量级生成1 000个加入高斯噪声的模型参数集合。图2为参数S加入量级为0.03高斯噪声后得到的带噪声集合的直方图。
图2 带高斯噪声(Std=0.03)的S直方图
从图2中可以看出,加入不同量级高斯噪声进行随机扰动后的参数数值的分布符合正态分布。
将带不同量级噪声的参数S集合输入到AIEM模型中,其他参数使用表1中的标定值,确定带参数S噪声的模拟后向散射系数集合,然后利用LUT法反演土壤水分,计算土壤水分反演值的RMSE,图3为带不同量级噪声的S参数对土壤水分反演结果的RMSE响应折线。
使用同样的方法确定带不同量级噪声的l和θ参数对土壤水分反演结果的RMSE响应折线,如图4和图5所示。
基于SAR的土壤水分反演一般是在给定入射角范围的情况下进行,由图3—图5可以看出,在雷达入射角为θ=30°情况下,以SMOS土壤水分产品的精度为标准,均方根高度S的误差量级至少要控制在0.04即标定值的8%以内,才能满足土壤水分反演精度的要求;同样,相关长度l的误差量级至少要控制在0.8即标定值的12%以内;雷达入射角的误差量级至少要控制在0.6即标定值的6%的范围内。计算S,l和θ这3个参数不同噪声量级集合的方差,并计算带不同量级噪声的参数集合反演得到的土壤水分的方差,见表4。
图3 不同噪声量级S对Mv的RMSE响应折线
图4 不同噪声量级l对Mv的RMSE响应折线
图5 不同噪声量级θ对Mv的RMSE响应折线
由表4可以看出,随着各参数噪声量级和参数集合方差的增大,根据各参数集合反演的土壤水分方差也随之增大,其中基于带噪声的参数S反演的土壤水分的方差变动幅度最为明显,说明由参数S的不确定性引起的反演土壤水分的不确定最大,而由l和θ的不确定性引起的反演土壤水分的不确定则相对较小。
在给定入射角θ=30°情况下的研究,由图6可以看出,入射角的变化对土壤水分反演精度有较大的影响,为定量描述不同入射角情形下敏感性参数在反演中的不确定性,入射角分别取θ=10°,20°,30°,40°,50°共5个值,以表1中的各参数的标定值为期望,不同参数使用不同标准差形成噪声量级:S(0.01,0.02,0.03,0.04,0.05,0.07,0.09,0.1),l(0.2,0.4,0.6,0.8,1,2,3,4),θ(0.2,0.4,0.6,0.8,1,2,3,4),每个噪声量级生成1 000个加入高斯噪声的模型参数集合。将带不同量级噪声的参数S集合输入到AIEM模型中,利用LUT法反演土壤水分,计算土壤水分反演值的RMSE,图6为不同入射角情形下带不同量级噪声的敏感性参数对土壤水分反演结果的RMSE响应折线。
表4 反演土壤水分不确定性量化统计
图6 不同入射角情形下带不同误差量级参数时Mv的RMSE响应折线
由图6可以看出,均方根高度S的噪声干扰在小入射角时对土壤水分反演结果的影响比大入射角时弱,以SMOS土壤水分产品的精度RMSE=0.04为标准,θ∈(10°,50°)时,S的误差量级至少要控制在0.03即标定值的2%以内,才能满足土壤水分反演精度的要求;同样,相关长度l的噪声干扰对土壤水分反演也有相似的影响,即在随入射角的增大,相关长度l的噪声干扰对土壤水分反演结果的影响也逐渐增大,l的误差量级至少要控制在0.6即标定值的8%以内。θ∈(10°,50°)时,不同入射角情形均方根高度S和l的误差量级控制范围见表5。
表5 不同Q时S和i的误差量级控制范围 %
为了测试不同入射角情况下敏感参数不确定性对反演结果的影响,本文选取了包含临泽地区的2008-05-24和2008-07-11一共2幅ASAR数据,入射角θ分别归一化为20°和30°,利用临泽样区中采样点的土壤水分实测数据,将均方根高度固定为1.4 cm,结合同期ASAR影像中的后向散射数据,利用LUT法反演得到采样点相应的有效相关长度和土壤水分。使用表5中敏感参数的控制范围为各敏感参数加入干扰噪声,然后利用LUT法反演得到带噪声的土壤水分,计算各幅ASAR影像中样区采样点土壤水分反演值的RMSE,见表6。
表6的结果表明,通过不同入射角情形下的实测数据验证,本文得到的不同入射角情形均方根高度S和l的误差量级控制范围是有效的。
表6 反演土壤水分RMSE
本文通过定量分析AIEM模型各参数的敏感性,选取均方根高度S、相关长度l和雷达入射角θ这3个敏感性参数进行不确定性分析,通过LUT方法研究反演土壤水分过程中的参数不确定性,确定不同入射角情形下S和l的误差量级控制范围。使用ENVISAT ASAR影像结合实测土壤水分数据,计算采样点土壤水分反演值的RMSE,结果表明,本文计算得到的不同θ时S和l的误差量级控制范围可有效控制土壤水分反演精度,并适用于较大的入射角范围。
[1] 李大治,晋锐,车涛,等.联合机载PLMR微波辐射计和MODIS产品反演黑河中游张掖绿洲土壤水分研究[J].地球科学进展,2014,29(2):295-305.
[2] ENGMAN E T, CHAUHAN N. Status of microwave soil moisture measurements with remote sensing[J]. Remote Sensing of Environment, 1995, 51(1):189-198.
[3] KEYSER E D, VERNIEUWE H, LIEVENS H, et al. Assessment of SAR-retrieved soil moisture uncertainty induced by uncertainty on modeled soil surface roughness[J]. International Journal of Applied Earth Observations & Geoinformation, 2012, 18(1):176-182.
[5] MA C, LI X, NOTARNICOLA C, et al. Uncertainty Quantification of Soil Moisture Estimations Based on a Bayesian Probabilistic Inversion[J]. IEEE Transactions on Geoscience & Remote Sensing, PP(99):1-14.
[6] 李大治.L波段土壤水分反演的不确定性分析及其反演策略研究——以HiWATER PLMR数据为例[D].北京:中国科学院寒区旱区环境与工程研究所,2014.
[7] VERHOEST N E C, LIEVENS H, WAGNER W, et al. On the Soil Roughness Parameterization Problem in Soil Moisture Retrieval of Bare Surfaces from Synthetic Aperture Radar[J]. Sensors, 2008, 8(7):4213.
[8] RAHMAN M M, MORAN M S, THOMA D P, et al. A derivation of roughness correlation length for parameterizing radar backscatter models[J]. International Journal of Remote Sensing, 2007, 28(18):3995-4012.
[9] DAVIDSON M, TOAN T L, MATTIA F, et al. On the characterization of agricultural soil roughness for radar remote sensing studies[J]. IEEE Transactions on Geoscience & Remote Sensing, 2000, 38(2):630-640.
[10] 郭立新,王蕊,吴振森.随机粗糙面散射的基本理论与方法[M].北京:科学出版社,2010.32-39.
[11] SHI J, CHEN K S, LI Q, et al. A parameterized surface reflectivity model and estimation of bare-surface soil moisture with L-band radiometer[J]. IEEE Transactions on Geoscience & Remote Sensing, 2010, 40(12):2674-2686.
[12] 余凡, 李海涛, 张承明,等. 利用双极化微波遥感数据反演土壤水分的新方法[J]. 武汉大学学报(信息科学版), 2014, 39(2):225-228.
[13] 陈晶,贾毅,余凡.双极化雷达反演裸露地表土壤水分[J].农业工程学报,2013,29(10):109-116.
[14] LIEVENS H, VERHOSET N E C, KEYSER E D, et al. Effective roughness modelling as a tool for soil moisture retrieval from C- and L-band SAR[J]. Hydrology & Earth System Sciences Discussions, 2010, 7(4):151-162.
[15] SALTELLI A, TARANTOLA S, CHAN K P S, et al. A quantitative model-independent method for global sensitivity analysis of model output[J]. Technometrics,1999, 41(1): 39-56.
[16] CUKIER R, FORTUIN C, SHULER K,et al.Study of the sensitivity of the coupled reaction systems to uncertainties in rate coefficients: I. Theory[J]. The Journal of Chemical Physics,1973, 59:3873
[17] 郝崇清, 王江, 邓斌, 等. 基于稀疏贝叶斯学习的复杂网络拓扑估计[J]. 物理学报, 2012, 61 (14):497-505.
[18] OH Y, HONG J Y. Effect of Surface Profile Length on the Backscattering Coefficients of Bare Surfaces[J]. IEEE Transactions on Geoscience & Remote Sensing, 2007, 45(3):632-638.
[19] SU Z, TROCH P A, DE TROCH F P. Remote sensing of soil moisture using EMAC/ESAR data[C]// Geoscience and Remote Sensing Symposium, 1996. IGARSS ’96.‘Remote Sensing for a Sustainable Future.’, International. IEEE, 1997:1303-1305 vol.2.
[20] BAGHDADI N, HOLAH N, ZRIBI M. Soil moisture estimation using multi‐incidence and multi‐polarization ASAR data[J]. International Journal of Remote Sensing, 2006, 27(10):1907-1920.
[21] ALVAREZ-MOZOS J, CASALI J, GONZALEZ-AUDICANA M, et al. Assessment of the operational applicability of RADARSAT-1 data for surface soil moisture estimation[J]. IEEE Transactions on Geoscience & Remote Sensing, 2006, 44(4):913-924.
[22] 甄珮珮. 基于粗糙度参数的风沙滩地区土壤水分微波遥感反演模型研究[D]. 西安:长安大学, 2016.