CN103149220B - 一种单频率微波辐射计土壤水分的反演方法 - Google Patents

Abstract

本发明提出一种单频率微波辐射计土壤水分的反演方法，属于微波遥感领域，该方法不是多参数同时反演方法，也不是经验回归方法，而是利用最少的辅助数据，单参数反演得到土壤水分，该方法克服了多参数反演存在的多解及计算量大的问题，也避免了LPRM算法中对粗糙度设定为定值而造成的反演误差，同时相比SCA算法对辅助数据的过分依赖，本方法需要利用的辅助数据仅为全球免费共享的土壤质地数据，极大的降低了反演土壤水分的业务化应用难度。

Description

一种单频率微波辐射计土壤水分的反演方法

技术领域

本发明涉及微波遥感领域，主要涉及的是一种单频率微波辐射计土壤水分的反演方法。

背景技术

土壤水分是地球生态系统非常重要的组成部分，它对全球水圈，大气圈和生物圈水分和能量的交换起着关键的作用。频繁的全球土壤水分观测对很多环境学科，例如洪水监测，改良土壤侵蚀模型，天气和气候预报以及农业应用都至关重要。因此，获取大范围地表土壤水分时空分布信息是一个迫切需要解决的问题。

传统的土壤水分测量方法是基于点的测量，不仅费时费力，而且很难实现大范围持续监测。而被动微波传感器的空间分辨率通常为几十千米，因此非常适合大尺度的全球性的土壤水分观测。且相比其他遥感手段获取土壤水分，被动微波遥感有其独特的优势：(1)能够全天时工作，(2)能够穿透云层并较少受到大气的影响，(3)通过土壤介电常数与土壤水分直接相关，(4)相对而言对地表粗糙度和植被覆盖不那么敏感。

发明内容

本发明提出了一个新的土壤水分反演方法。该方法基于零阶辐射传输模型和植被光学厚度不受极化依赖的假设。方法中将植被和粗糙度的影响综合为一个影响参数，并利用同频率下H/V两个极化通道将这个综合影响参数表达为只受到土壤水分影响的表达式。最后将这个表达式代回到辐射传输方程中，并利用微波极化差异指数(microwavepolarizationdifferenceindex，MPDI)消去地表温度的影响，进而通过利用非线性迭代算法使得观测的MPDI与模拟的MPDI之间的差的绝对值最小来得到土壤水分。

该反演方法包括如下步骤：

(1)首先设定待反演参数土壤水分mv在一个合理的物理范围区间，然后设定它的迭代步长，所述迭代步长将所述物理范围区间划分成多个网格；

(2)将土壤质地数据与步骤(1)中的土壤水分一同输入Hallikainen经验模型，得到土壤的复介电常数；

(3)将步骤(2)中得到的土壤复介电常数结合卫星的传感器入射角参数输入菲涅尔方程式，得到同频率下H/V两个极化的光滑地表反射率；

(4)将步骤(3)中得到的光滑地表反射率结合粗糙度参数Q、h，输入粗糙度模型Q/H模型，得到粗糙地表的有效反射率；

(5)将步骤(4)中的粗糙地表反射率与植被单次散射反照率ω、植被光学厚度τ_p、土壤和植被冠层的温度T_s和T_c一同输入零阶辐射传输模型，得到给定频率下前向模型模拟的H/V极化的亮温MPDI_sim；

(6)在相同频率下，利用卫星真实观测的H/V极化亮温值，计算得到观测的MPDI_obs，并给出循环收敛的判别准则|MPDI_obs-MPDI_sim|＝min，即在mv不断迭代的过程中，遍历所有的mv迭代值，取所有值中能使得观测的MPDI_obs与模拟的MPDI_sim之间的差的绝对值最小的那个mv值，此时便得到了一个格网对应的mv反演值；

(7)重复步骤(1)～(6)便可以得到所述物理范围内所有格网的my反演值。

优选的是，在步骤(1)中，所述物理范围区间为0～0.6，所述迭代步长为0.001，即每求解一个格网的mv的迭代过程均为：my循环从0开始至0.6结束，每次迭代的步长为0.001，最终mv的取值为满足条件的那个值。

优选的是，在步骤(4)中，其中Q参数选择在6.9GHz时的标定值0.174。

优选的是，在步骤(5)中，采取单次散射反照率为零，即ω＝0，植被光学厚度不受极化影响，即τ_H＝τ_V，植被、土壤温度相等，即T_s＝T_c＝T，将粗糙度和植被的影响合并成一个综合的影响因子，即exp(-2τ-h)，然后通过利用同频率下H/V极化双通道的辐射传输方程将这个综合影响因子表达成只受到土壤水分影响的表达式，即exp(-2τ-h)＝f(mv)，因为my在整个算法过程中一直都是一个迭代的已知值，因此粗糙度和植被的影响exp(-2τ-h)＝f(my)也可以计算得到，也就是说，除了温度未知，其他的参数都“已知”，接着，利用微波极化差异指数MPDI消去温度的影响，进而可以得到模拟的MPDI_sim。

本研究提出的土壤水分反演算法，既不是多参数同时反演算法，也不是经验回归算法，而是利用最少的辅助数据(仅土壤质地数据)单参数反演得到土壤水分。该算法克服了多参数反演存在的多解及计算量大的问题，也避免了LPRM算法中对粗糙度设定为定值而造成的反演误差，同时相比SCA算法对辅助数据的过分依赖，本算法需要利用的辅助数据仅为全球免费共享的土壤质地数据，极大的降低了反演土壤水分的业务化应用难度。

本算法最大优势在于同时考虑了植被和粗糙度对地表辐射亮温的影响，在土壤水分的反演过程中不需要利用除土壤质地数据外的任何辅助数据(如地表粗糙度，植被含水量等)。同时，该方法是完全基于物理的，没有区域依赖性，且只需要利用一个频率的亮温观测值即可进行土壤水分的估算。因此，该算法同样能应用于SMOS卫星和SMAP计划来进行全球土壤水分的反演。

附图说明

图1为利用站点实测的土壤水分对算法反演的土壤水分及NASA官方土壤水分产品进行长时间序列的对比和验证，(a)2008.7.1～2008.10.31；(b)2009.4.1～2009.10.31；(c)2010.4.1～2010.7.31。

图2为利用站点实测的土壤水分对算法反演的土壤水分及NASA官方土壤水分产品进行散点图的对比和验证。

图3为土壤水分反演方法流程图。

具体实施方式

本算法基于零阶段辐射传输模型，即通常所说的τ-ω模型。在这个模型中，电磁波在大气及植被层间多次散射的影响被忽略。模型的公式如(1)所示：

T_Bp＝(1-R_sp)·T_s·exp(-τ_p)+(1-ω)·T_c·[1-exp(-τ_p)]·[1+R_sp·exp(-τ_p)](1)

其中下标p表示水平或垂直极化，R_sp表示土壤有效反射率，T_s和T_c分别表示土壤和植被冠层的温度，τ_p表示观测方向的植被光学厚度，ω表示植被的单次散射反照率。

为了简化反演过程，在算法中假设ω等于零，并假设τ_p不受极化的影响，这种假设在卫星尺度下被广泛应用。同时，进一步假设植被温度和土壤温度相等，即T_s＝T_c＝T。则公式(1)可以简化为如下公式(2)所示：

T_Bp＝T·[1-R_sp·exp(-2τ)](2)

粗糙地表的有效反射率可以通过Q/H模型计算得到，Q/H模型的计算公式为：

R_sp＝[(1-Q)·R_op+Q·R_oq]·exp(-h)(3)

其中下标p和q表示不同的极化方式，Q和h分别表示Q/H模型中的两个粗糙度参数，将公式(3)代入公式(2)，并进行部分系数的合并，可以得到如下结果：

T_Bp＝T·[1-R′_sp·exp(-2τ-h)](4)

其中R′_sp＝(1-Q)·R_op+Q·R_oq。可将Q参数固定为一个全球的定标参数，利用h参数去刻画地表粗糙度的空间变异性。在本算法中，Q参数选择在6.9GHz时已经标定的值，即0.174。因此，根据菲涅尔方程，在特定极化方式的频率和入射角下，R′_sp只与土壤的介电特性有关。算法选择利用Hallikainen经验模型来实现土壤水分和土壤介电常数之间的转换。

所述Hallikainen经验模型公式为：

${\mathrm{\ϵ}}^{\′}=(a_0+a_{1}S+a_{2}C)+(b_0+b_{1}S+b_{2}C)m_v+(C_0+c_{1}S+c_{2}C)m_v^2$

其中，a₀，a₁，a₂，b₀，b₁，b₂，c₀，c₁，c₂为标定系数，S和C分别表示沙土和粘土的百分含量。在不同的频率下，a₀，a₁，a₂优选的值如表1所述。

因此，如果土壤质地数据已知，那么R′_sp可以表示为土壤水分的函数(即R′_sp＝f(mv))。同时，在公式(4)中，由于植被和粗糙度的影响都是以指数的形式出现的，因此可以将两者综合为一个影响参数，即exp(-2τ-h)。进一步令exp(-2τ-h)＝a，那么在相同频率下，公式(4)中H/V极化的亮温可以表示为：

T_Bh＝T·(1-R′_sh·a)(5)

T_Bv＝T·(1-R′_sv·a)(6)

结合公式(5)和公式(6)，植被及粗糙度的综合影响参数a可以表示为：

$a=\frac{T_{Bv}-T_{Bh}}{T_{Bv}\·R_{sh}^{\′}-T_{Bh}\·R_{sv}^{\′}}---\left(7\right)$

从(7)式可以看出，在H/V极化亮温已知的情况下，综合影响参数a可以表示为土壤水分的函数，即a＝exp(-2τ-h)＝f(my)。将公式(7)回代入公式(5)和(6)可以得到：

$T_{Bh}=T\·(1-R_{sh}^{\′}\·\frac{T_{Bv}-T_{Bh}}{T_{Bv}\·R_{sh}^{\′}-T_{Bh}\·R_{sv}^{\′}})---\left(8\right)$

$T_{Bv}=T\·(1-R_{sv}^{\′}\·\frac{T_{Bv}-T_{Bh}}{T_{Bv}\·R_{sh}^{\′}-T_{Bh}\·R_{sv}^{\′}})---\left(9\right)$

从公式(8)和公式(9)可以看出，相同频率下，H，V极化亮温只受到土壤水分和地表温度的影响。另一方面，微波极化差异指数MPDI在许多研究中常被用来消去地表温度的影响。MPDI的计算公式如下所示：

$MPDI=\frac{T_{Bv}-T_{Bh}}{T_{Bv}+T_{Bh}}---\left(10\right)$

将公式(8)和公式(9)代入公式(10)，可以得到前向辐射传输模型模拟的MPDI，即MPDI_sim，可以看出MPDI_sim只与土壤水分有关。最后，可以构造观测的MPDI与模拟的MPDI之间的差的绝对值为前向过程，利用非线性迭代算法使前向过程最小，从而得到土壤水分。

参见附图3所示，本方法的具体实现过程如下：

(1)首先设定待反演参数土壤水分(mv)在一个合理的物理范围区间：如0～0.6，然后设定它的迭代步长，如0.001；(即每求解一个格网的mv的迭代过程均为：mv循环从0开始至0.6结束，每次迭代的步长为0.001，最终mv的取值为满足条件的那个值)

(2)将土壤质地数据(免费共享)与步骤(1)中的土壤水分一同输入Hallikainen经验模型，得到土壤的复介电常数；

(3)将步骤(2)中得到的土壤复介电常数结合卫星的传感器参数(具体指的是入射角)输入菲涅尔方程式，得到同频率下H/V两个极化的光滑地表反射率；

(4)将步骤(3)中得到的光滑地表反射率结合粗糙度参数Q、h，输入粗糙度模型Q/H模型，得到粗糙地表的有效反射率，其中Q参数取标定值0.174(6.9GHz)；

(5)将步骤(4)中的粗糙地表反射率与其他参数如植被单次散射反照率ω、植被光学厚度τ_p，土壤和植被冠层的温度T_s和T_c一同输入零阶辐射传输模型，得到给定频率下前向模型模拟的H/V极化的亮温。为了简化反演过程，采取单次散射反照率为零(ω＝0)，植被光学厚度不受极化影响(τ_H＝τ_V)及植被、土壤温度相等(T_s＝T_v＝T)的假设。在这些合理假设的基础上，通过对辐射传输方程的巧妙变形，将粗糙度和植被的影响合并成一个综合的影响因子(即exp(-2τ-h))。然后通过利用同频率下H/V极化双通道的辐射传输方程将这个综合影响因子表达成只受到土壤水分影响的表达式(即exp(-2τ-h)＝f(my))。因为mv在整个算法过程中一直都是一个迭代的已知值，因此粗糙度和植被的影响exp(-2τ-h)＝f(my)也可以计算得到。也就是说，除了温度未知，其他的参数都“已知”。接着，利用微波极化差异指数MPDI消去温度的影响，进而可以得到模拟的MPDI(MPDI_sim)；

(6)在相同频率下，利用卫星真实观测的H/V极化亮温值，计算得到观测的MPDI(MPDI_obs)，并给出循环收敛的判别准则|MPDI_obs-MPDI_sim|＝min。即在mv不断迭代的过程中，遍历所有的mv迭代值，取所有值中能使得观测的MPDI与模拟的MPDI之间的差的绝对值最小的那个mv值。此时便得到了一个格网对应的mv反演值。

(7)重复步骤(1)～(6)便可以得到研究范围内所有格网的mv反演值。

该方法能够同时考虑植被和粗糙度对地表辐射亮温的影响，只需要利用最少的辅助数据和一个频率的亮温观测值就能对地表土壤水分进行很好的估算。为了进一步展示算法的效果和相比于其他算法的优势。研究选择利用中国青藏高原玛曲土壤水分监测网络的土壤水分实测数据对本算法反演的土壤水分和美国国家航空航天局(NASA)的AMSR-E官方土壤水分产品进行对比和验证。因为卫星观测和地面站点观测在水平尺度(主要对应空间分辨率)和垂直尺度(主要对应观测深度)都不相同，因此对于算法精度的验证往往采用算法反演值与实测数据值进行时间序列的比较，以便观察它们在一个长时间范围内的吻合程度和变化趋势。在本方法中，选择对比的时间段为2008年7月1日至2008年10月31日，2009年4月1日至2009年10月31日，2010年4月1日至2010年7月31日，其中年份中间缺省的时期为实验区的冰冻期，此时土壤已经冻结，传统的土壤介电模型不再适用，因此这些天数被排除在实验之外。同时相比白天，夜晚时土壤、植被冠层和空气的温度更加均衡稳定，因此研究选择对比的时间为ASMR-E降轨时间(凌晨1点30时)。利用站点实测的土壤水分对算法反演的土壤水分及NASA官方土壤水分产品进行长时间序列的对比和验证，如图1所示：(其中SM为soilmisture即土壤水分的缩写)

从图1可以很明显的看出，相比NASAAMSR-E的官方土壤水分产品，利用算法反演得到的土壤水分值无论是在整体的趋势上还是具体的数据值上，都与地面实测数据更加一致。玛曲土壤水分监测网络由于地处潮湿寒冷气候，在每年的4月至10月期间有频繁的降雨，因此导致土壤水分随时间有较为明显的变化。算法反演的土壤水分值能够很好的捕捉到土壤水分随时间的变化趋势。而NASA官方土壤水分产品的数值在绝大多数时都稳定在0.15左右，几乎没有动态变化，这与实际情况明显不符。另一方面，标准偏差(Standarddeviation，STD)经常用来衡量样本的动态变化，表2对实测土壤水分值，算法反演土壤水分值和NASA官方土壤水分产品的标准偏差进行了统计。从表2中可以看出，算法反演值的STD与实测水分值的STD比较接近，而NASA官方土壤水分的STD则很小，仅为0.13，远远低于实测水分值的STD，进一步证明算法反演的土壤水分值在捕捉地表土壤水分的动态变化上要明显优于NASA官方土壤水分产品。

表2算法反演土壤水分值及NASA官方土壤水分产品相对于实测土壤水分值的误差统计：

STD表示标准偏差，Bias表示偏差，MAE表示平均绝对误差，RMSE表示均方根误差

为了对比更加全面，进一步对算法反演的土壤水分值及NASA官方土壤水分产品与实测土壤水分值进行散点图对比，如图2所示，并对两者的误差进行统计，如表2所示。从图2和表2可以看出，算法反演的土壤水分精度要明显优于NASA官方土壤水分产品。前者更靠近1比1直线，在各个时期的均方根误差(RMSE)均低于0.1。而NASA官方算法则明显低估了地表的土壤水分值，整个时期的均方根误差高达0.235，说明了该算法在青藏高原地区并不适用。表2中其他的误差统计量包括偏差(Bias)，平均绝对误差(MAE)则进一步显示了本研究的算法相比NASA官方算法在反演土壤水分的精度上有了全面的改善和提高。

因为在算法中通过对辐射传输方程的巧妙变形，将植被和粗糙度的影响合并成了一个综合的影响因子，这样就将原来需要求解的两个未知数变成了一个。然后通过利用同频率下H/V极化双通道的辐射传输方程将这个综合影响因子表达成了只受到土壤水分影响的表达式，这样就将此综合参数与待反演参数联系起来，避免了多个参数的同时求解，且无需对植被和粗糙度的影响做较大的假设。最后算法通过利用微波极化差异指数MPDI消去地表温度的影响，进而利用非线性迭代算法使得观测的MPDI与模拟的MPDI之间的差的绝对值最小来得到土壤水分。这样就不需要利用辅助数据或者ka波段进行温度的反演，有效的较少了算法对辅助数据的依赖性和卫星频段的使用，经过细微的调整，该算法便可直接应用于最新的只具备单个波段的SMOS卫星和SMAP计划进行土壤水分的反演。同时，由于该算法是基于辐射传输模型提出的，具有很强的理论基础，所以该方法没有地区依赖性，相比经验回归算法能够更好的应用于全球的土壤水分反演。

Claims (4)

1.一种单频率微波辐射计土壤水分的反演方法，其特征在于，该方法包括如下步骤：

(1)首先设定待反演参数土壤水分mv在一个合理的物理范围区间，然后设定它的迭代步长，所述迭代步长将所述物理范围区间划分成多个网格；

(2)将土壤质地数据与步骤(1)中的土壤水分一同输入Hallikainen经验模型，得到土壤的复介电常数；

(3)将步骤(2)中得到的土壤复介电常数结合卫星的传感器入射角参数输入菲涅尔方程式，得到同频率下H/V两个极化的光滑地表反射率；

(4)将步骤(3)中得到的光滑地表反射率结合粗糙度参数Q、h，输入粗糙度模型Q/H模型，得到粗糙地表的有效反射率；

(5)将步骤(4)中的粗糙地表反射率与植被单次散射反照率ω、植被光学厚度τ_p、土壤和植被冠层的温度T_s和T_c一同输入零阶辐射传输模型，得到给定频率下前向模型模拟的H/V极化的亮温MPDI_sim；

(6)在相同频率下，利用卫星真实观测的H/V极化亮温值，计算得到观测的MPDI_obs，并给出循环收敛的判别准则|MPDI_obs-MPDI_sim|＝min，即在mv不断迭代的过程中，遍历所有的mv迭代值，取所有值中能使得观测的MPDI_obs与模拟的MPDI_sim之间的差的绝对值最小的那个mv值，此时便得到了一个格网对应的mv反演值；

(7)重复步骤(1)～(6)便可以得到所述物理范围内所有格网的mv反演值。

2.基于权利要求1所述的方法，其特征在于，在步骤(1)中，所述物理范围区间为0～0.6，所述迭代步长为0.001，即每求解一个格网的mv的迭代过程均为：mv循环从0开始至0.6结束，每次迭代的步长为0.001，最终mv的取值为满足收敛条件的那个值。

3.基于权利要求1所述的方法，其特征在于，在步骤(4)中，其中Q参数选择在6.9GHz时的标定值0.174。

4.基于权利要求1所述的方法，其特征在于，在步骤(5)中，采取单次散射反照率为零，即ω＝0，植被光学厚度不受极化影响，即τ_H＝τ_V，植被、土壤温度相等，即T_s＝T_c＝T，将粗糙度和植被的影响合并成一个综合的影响因子，即exp(-2τ-h)，然后通过利用同频率下H/V极化双通道的辐射传输方程将这个综合影响因子表达成只受到土壤水分影响的表达式，即exp(-2τ-h)＝f(mv)，因为mv在整个算法过程中一直都是一个迭代的已知值，因此粗糙度和植被的影响exp(-2τ-h)＝f(mv)也可以计算得到，也就是说，除了温度未知，其他的参数都已知，接着，利用微波极化差异指数MPDI消去温度的影响，进而可以得到模拟的MPDI_sim。