地表参数遥感反演的前提就是建立正向模型, 即根据物理基础建立起地表参数与卫星获得的亮温之间的函数关系: 该模型通过输入地表参数, 如土壤水分、植被含水量和地表温度, 根据亮温与地表和大气相关变量间的关系, 得到相应的亮温数据。一般的正向模型将卫星与地表间的介质分为大气、植被和土壤3层, 植被层可看作是一个位于粗糙土壤之上的单次散射层, 因而需要从地表、植被和大气3方面考虑。地表土壤的发射率通过介电常数和观测角与土壤水分相关。除此之外, 植被覆盖度、地表粗糙度、土壤温度以及地形等因素都对土壤发射率有影响。

在粗糙地表情况下, 传统模拟裸露地表微波辐射特性的方法是采用Q/H模型以及一系列由此衍生出来的模型。本研究引入Shi等在2006年研究并提出的Q_p参数化模型。该模型以AIEM模型提供的模拟数据库为参照, 对比传统的Q/H模型和Q_p参数化模型在H、V两个极化和极化比值下的模拟精度, 发现Q_p参数化模型更适合AMSR-E等大入射角辐射计, 并且给出了从6~37 GHz的Q_p计算公式。该公式建立了Q_p参数与粗糙度参数s(均方根高度)和l(相关长度)的商(即s/l)之间的函数关系^[5]。

在有植被覆盖的情况下, 选择τ -ω 模型来模拟植被顶层的亮温 Tbp7。该模型适用于低矮植被覆盖地表(如农田), 且对低频下植被覆盖地表的微波辐射描述逼近真实情况。植被的覆盖虽然减弱了来自土壤的辐射, 但同时增加了植被本身的辐射, 使得土壤水分变化对植被顶层亮温的变化贡献减小, 即在有植被覆盖的时候, 土壤水分的敏感性降低而植被含水量的敏感性升高, 从而土壤水分反演变得困难。此时的土壤水分反演思路主要有两种: 一是估算并去除植被的影响, 转化为在裸露地表下反演土壤水分; 二是将土壤水分和植被含水量同时反演。本研究采用第二种思路。

卫星观测的亮温实际上是大气顶层的亮温, 在微波辐射信号从地表到达传感器的传输过程中, 还受到大气的影响, 一般采用经验关系来消除^[7]。

参数敏感性分析是为了辅助反演时选择合适的通道组合以及分析反演结果, 通常表示为在给定的参数基准值x₀和变化范围X_j下, 归一化的参数敏感性系数值^{[1, 2]}。由于引入了Q_p模型代替Q/H模型作为正向模拟的模型, 同时增加了一个待反演参数— — 粗糙度, 并且扩大了植被含水量的范围, 因此, 以下将要重新分析所有通道(除了89 GHz以外)对各个待反演参数的敏感性, 并将其与Njoku等在6~18 GHz 3个通道的工作进行对照。

表3列出了裸露地表情况下AMSR-E 10个通道(除89 GHz以外)的亮温对地表粗糙度(m)、土壤含水量(s_w)、地表温度(t_s)的敏感性。当实际地表为裸露土壤或接近裸露土壤(植被覆盖率很低)时, 相比植被含水量, 更受关注的地表参数是m、s_w和t_s。为了简化正向模型, 这里固定植被含水量v_wc=0 kg/m²; 为了方便比较不同频率和极化的地面发射率对这3个参数的敏感性, 表中的亮温是尚未考虑大气影响的近地面亮温, 与表4和表5中针对植被覆盖地表的大气顶层的亮温是不同的。这样的前提并不影响对反演算法的检验。

表3

Tab.3

表3(Tab.3)

表3 裸露地表情况下的各参数敏感性 Tab.3 Parameters’ sensitivities in case of bare soil surface 参数 基准值 x0 变化范围 Xj 敏感性/K 6H 6V 10H 10V 18H 18V 23H 23V 37H 37V m=s/l 0.4 0.3 12.6 2.3 12.3 2.1 11.1 1.9 10.7 1.8 10.2 1.7 sw/(m3· m-3) 0.02 0.42 259.5 127.9 209.9 97.1 142.0 59.0 117.6 47.1 86.0 32.9 ts/℃ 0 20 12.3 16.8 11.2 16.5 12.3 17.1 13.4 17.6 15.2 18.4

表3 裸露地表情况下的各参数敏感性 Tab.3 Parameters’ sensitivities in case of bare soil surface

从表3可以看出, m和s_w在H极化时的敏感性明显强于V极化时的敏感性; t_s在两种极化时的敏感性相当。m和t_s在各个频率的敏感性相当, s_w在低频的敏感性明显强于在高频的敏感性。所以, 同样的通道数目, 低频组合比高频组合更容易获得较好的反演结果。此外, 由于正向模型包括了Q_p模型, 则在H极化的粗糙地表反射率中V极化的光滑地表反射率的贡献不为0, 同时使用双极化数据可以消除极化差, 较之单一的H极化的频率组合更容易获得较好的反演结果。

在有植被覆盖的地表, 随着植被覆盖度的增大, v_wc会增加, 同时冠层对地表辐射的衰减效应增强。表4和表5分别列出了v_wc≤ 1.5 kg/cm²和v_wc> 1.5 kg/cm²两种情形下大气顶层亮温对s_w、v_wc及t_s的敏感性。

表4

Tab.4

表4(Tab.4)

表4 低植被覆盖情况下的各参数敏感性 Tab.4 Parameters’ sensitivities in case of lower vegetation cover surface 参数 基准值 x0 变化范围 Xj 敏感性/K 6H 6V 10H 10V 18H 18V sw/(m3· m-3) 0.02 0.42 282.1 121.9 230.4 91.7 156.8 54.8 vwc/(kg· cm-2) 0.5 1.0 31.7 43.8 52.7 71.5 95.8 122.8 ts/℃ 0 20 11.5 17.1 10.2 16.7 11.4 17.4

表4 低植被覆盖情况下的各参数敏感性 Tab.4 Parameters’ sensitivities in case of lower vegetation cover surface

表5

Tab.5

表5(Tab.5)

表5 高植被覆盖情况下的各参数敏感性 Tab.5 Parameters’ sensitivities in case of higher vegetation cover surface 参数 基准值 x0 变化范围 Xj 敏感性/K 6H 6V 10H 10V 18H 18V sw/(m3· m-3) 0.02 0.42 182.7 79.0 122.3 48.6 65.2 22.8 vwc/(kg· cm-2) 1.5 4.0 114.5 137.8 169.7 192.2 229.5 245.0 ts/℃ 0 20 9.8 13.4 7.9 11.3 6.2 8.7

表5 高植被覆盖情况下的各参数敏感性 Tab.5 Parameters’ sensitivities in case of higher vegetation cover surface

从表4、5可以看出, s_w的敏感性在H极化时显著强于在V极化时的; v_wc和t_s的敏感性在两种极化时基本相当; 敏感系数最小的是t_s, 预测3个参数反演结果误差最大的是t_s。随着v_wc值的变大, s_w和v_wc的敏感性反而大大变弱, t_s的敏感性也变弱, 反演的难度增加; s_w敏感性强弱变化趋势为: 6GHz> 10GHz> 18GHz, 而v_wc敏感性强弱的变化趋势为: 6GHz< 10GHz< 18GHz, t_s的敏感性随频率改变的变化幅度很小。

为了与Njoku算法进行对比, 本文所采用的迭代算法仍然是基于Levernberg-Marquardt原理^[2]。

这里的反演问题归结为一个结合了微波遥感定量反演理论的非线性系统的最优化问题, 涉及到一些数学原理, 故考虑使用MATLAB来辅助实现, 以便于调用其工具箱的矩阵运算和求解方程组等函数。MATLAB的最优化工具箱提供的lsqnonlin函数虽然可以选择Levernberg-Marquardt方式来实现优化, 但层层调用, 且封装了其他的优化方式, 不利于修改和调用。本研究在MATLAB帮助菜单下的web source中找到Miroslav Balda提供的LMFnlsq函数^[6], 参考并自定义算法。

算法的输入参数为最大迭代次数、阻尼系数初始值、待反演参数值的边界范围、代价函数的容差及待反演参数的容差, 具体实现步骤如下:

(1) 定义代价函数Fun和雅克比矩阵J。Fun即各个通道的亮温观测值(用模拟值加上高斯随机噪声表示)与用估计值x_k(第k次迭代后的x)计算的亮温估计值之差的加权平方和, 用各个通道观测噪声标准差平方的倒数作为权重; J可以用解析函数的一阶泰勒导数或者有限数值差分求导。由于模型公式的复杂性, 若直接计算一阶偏导数作为雅克比矩阵, 不但计算复杂且可能会导致无法收敛, 所以这里采用后一种方式, 但其缺陷是计算量大, 速度慢。

(2) 初始值设定。初始值包括待反演参数x₀(可以任意选取, 但为了尽快收敛, 一般取它的中间值)、迭代次数最大值MaxIter、阻尼系数l(缺省值是0。接近0时, 切换到GN方法来搜索解; 接近无穷大时, 切换到SD方法来搜索解。如果代价函数的值逐渐变小, 1/scale变小, 加快收敛速度, 反之, scale变大, 增加解的搜索范围, scale的取值一般在2~10之间)、代价函数容差FunTol及待反演参数容差XTol。

(3) 调用函数。在MATLAB下自定义。

(4) 限制待反演参数值的边界范围。边界范围参见表2。如果迭代反演的结果超出边界范围, 则调整算法的输入参数, 重复步骤(3)。

(5) 满足结束条件时(达到最大迭代次数、小于代价函数容差FunTol及小于待反演参数容差XTol), 结束迭代。

在裸露地表条件下, 待反演参数有综合粗糙度(m)、土壤含水量(s_w)及地表温度(t_s)。此时在反演误差的贡献中少了植被层的影响, 土壤水分反演的精度较高。实验发现, 在6GHz和10GHz的双极化观测组合下, 地面温度的反演误差较大。这里假设已经获得了6GHz、10GHz、18GHz、23GHz及37GHz双极化的大气参数。大气校正后, 由大气顶层的亮温得到近地面的亮温, 然后再进行反演。实验发现, 加入37GHz双极化观测后, 地表温度反演的精度明显提高, 这与前面的参数敏感性分析结果是一致的。4通道算法的最佳频率组合是6GHz和37GHz的双极化, 共4个通道的近地面亮温, 并且, 10个通道的观测组合得到的反演精度最高, 图1显示了分别添加0 K和0.3 K的高斯白噪声后3个地表参数的反演结果。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 裸露地表情况下m、sw、ts的反演结果Fig.1 Retrieval results of m、sw and ts in case of bare soil surface

从图1可见, 高斯随机噪声从0 K增加到0.3 K时, 反演精度下降, 说明该算法对噪声比较敏感。表6列出了不同噪声条件下, 3个地表参数的反演精度。

表6

Tab.6

表6(Tab.6)

表6 不同高斯噪声条件下的地表参数反演精度 Tab.6 Retrieval RMSE of m、sw and st under various guass noise conditions RMSE 噪声/K 0 0.3 RMSE(sw)/(m3· m-3) < 0.000 1 0.002 RMSE(m) < 0.000 1 0.006 RMSE(ts)/℃ < 0.000 1 0.773

表6 不同高斯噪声条件下的地表参数反演精度 Tab.6 Retrieval RMSE of m、sw and st under various guass noise conditions

任意改变初始值, 反演结果基本不变。将反演结果输入到正向模型中, 各个通道亮温的方差也均低于添加的噪声, 表明多通道反演具有降噪效果。

在植被覆盖条件下, 待反演参数有土壤含水量(s_w)、植被含水量(v_wc)、地表温度(t_s)。Njoku等在算法的模拟研究和实测数据验证中发现, 当v_wc< 1.5 kg/m²时, 数据反演结果误差较小, 并且在模拟研究中, 在给模拟数据库添加了0.3 K的高斯白噪声后, 3个参数的反演误差分别小于0.06 g/cm³、0.1 kg/m²和1.5℃。但该算法不适合反演v_wc> 1.5 kg/m²的数据。笔者认为其原因为: ①反演算法的实现仅从数学原理出发, 没有考虑遥感模拟数据库的参数范围以及遥感模型的适用范围, 迭代可能得到了无意义解, 从而在调用遥感模型代码时提前跳出程序; ②当植被含水量大于1.5 kg/m²时, 辐射亮温对土壤含水量的敏感性大大减小, 而代价函数在逼近最小值的过程中存在多解, 在满足某一个迭代结束条件之前, 反演得到的结果随机性较大。

针对上述情况, 本研究对反演方法进行了改进, 具体措施为:

(1)对地表参数进行抽样, 试用反演算法, 获得最高反演精度所对应的一些参数配置, 并将这些参数配置作为经验值;

(2)在反演算法中增加对反演参数值边界条件的限制(迭代算法的缺陷就在于对超定方程组求最优解时会产生多解, 为了保证最优解落在地表参数的物理范围或者实际范围内, 给出了一个数学角度上的修正方案, 即一旦超过范围就取上下限), 假如迭代点的土壤水分含量小于0, 则令它为0.01, 如果大于0.45, 则令它为0.44;

(3)计算有限差分雅克比矩阵代替。如果将正演的植被模型看作一个关于待反演的3个地表参数的函数, 这个函数是非常复杂的非线性函数, 不方便计算每次迭代后3个地表参数的偏导数, 所以在计算雅克比矩阵时, 用有限差分雅克比矩阵代替, 这样设计有利于保证反演结果的收敛和算法的稳定性, 但是算法的时间和空间复杂度变大。

从2.2节参数敏感性分析中得到用3个低频、双极化、共6个通道数据一起参加反演能获得最好的反演结果的推断, 在实验中证实了这一判断的正确性。当用6个通道同时反演3个参数时, 图2显示了分别在噪声0 K和0.3 K的条件下的反演结果。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 植被覆盖地表情况下sw、ts、vwc的反演结果Fig.2 Retrieval results of swts and vwc in case of vegetation cover surface

从图2可见, 在高斯白噪声为0 K时的2 772个点中, 只有几个点反演结果误差过大, 准确率为99.75%; 高斯白噪声为0.3 K时, 表7列出了3个地表参数在不同的v_wc的变化范围内的反演精度, 并与Njoku等的结果进行了对比。由于篇幅所限, Njoku等的反演结果请参考文献[1, 2], 这里仅列出其反演精度以便进行对照。

分析认为, 出现少数反演误差比较大的点, 主要是因为迭代反演的不稳定性。v_wc和s_w的敏感性强弱相当, 且明显强于t_s, 所以反演精度RMSE(s_w)≈ RMSE(v_wc)> RMSE(s_t)。随着v_wc的增大, v_wc和s_w的敏感性都变弱, 反演误差相对较大的点集中在v_wc值靠近边界最大值的一端, 这一现象在噪声增大为0.3 K时更明显, 同时也说明算法对噪声比较敏感。

表7

Tab.7

表7(Tab.7)

表7 不同vwc变化范围内3个地表参数的反演精度(高斯白噪声为0.3 K) Tab.7 Retrieval RMSE of sw、vwc and ts under various vwc value range (Noise=0.3 K) RMSR 0.5≤ vwc≤ 1.5 1.5< vwc≤ 5.5 0.5≤ vwc≤ 5.5 0≤ vwc≤ 1.5(Njoku) RMSE(sw)/(m3· m-3) 0.012 0.017 0.016 0.06/1.2=0.05 RMSE(vwc)/(kg· m-2) 0.011 0.037 0.036 0.10 RMSE(ts)/℃ 1.170 1.490 1.310 1.50

表7 不同vwc变化范围内3个地表参数的反演精度(高斯白噪声为0.3 K) Tab.7 Retrieval RMSE of sw、vwc and ts under various vwc value range (Noise=0.3 K)

任意改变初始值, 反演结果基本不变。将反演结果输入到正向模型中, 各个通道亮温的方差也均低于添加的噪声, 表明多通道的反演方法取得了降噪效果。