CN106226260A - 一种结合微波和红外遥感影像的土壤水分反演方法 - Google Patents

Abstract

本发明是一种结合微波和红外遥感影像融合的土壤水分反演方法。本发明结合PDI指数和VSWI指数的适用性特点，根据地表植被覆盖度的不同，建立了基于综合干旱指数（CDI）的遥感监测模型。先使用红外遥感数据得到综合干旱指数结果，然后结合微波遥感数据，求得综合干旱指数反演结果与微波遥感数据的相关系数，建立回归方程，影像融合同时对微波数据进行降尺度，结果可将红外数据反演的综合干旱指数结果转化为真实的土壤水分值，并且用微波遥感数据对红外数据的空值区域进行填补，即可得到完整的高分辨率土壤水分产品。该方法简单实用，提高了大范围土壤水分反演的效率。

Description

一种结合微波和红外遥感影像的土壤水分反演方法

技术领域

本发明属于遥感图像处理领域，特别涉及一种结合微波和红外遥感影像融合的土壤水分反演方法。

背景技术

陆表土壤含水量不仅是旱情监测的重要指标，还是气候、水文、生态、农业等领域的重要参数，也是全球气候变化的重要组成部分。遥感技术的发展为获取大范围地表土壤水分信息提供了有效手段。

光学/热红外遥感主要基于土壤及土壤上植被的光谱反射特性及表面发射特性来估算土壤水分。当研究区域较大时，尽管光学/热红外遥感反演地表土壤水分有着较高的空间分辨率，但受云层的影响严重，造成了结果在空间上的不连续，不能有效对土壤水分进行业务化的监测。被动微波遥感有着较好的物理基础，受云层的影响不大，但星载被动微波传感器的空间分辨率相对较低，应用受到一定限制。本发明结合红外和微波遥感数据，旨在得到完整的高分辨率土壤水分产品，能提高大范围区域的陆表土壤水分监测效率，便于土壤水分产品的业务化生产。

发明内容

本发明根据土壤植被覆盖度的不同，综合使用垂直干旱指数(PDI)和植被供水指数(VSWI)建立综合干旱指数(CDI)监测模型。先使用红外遥感数据根据综合干旱指数模型得到反演结果，然后结合微波遥感数据，求得综合干旱指数模型反演结果与微波遥感土壤水分数据的相关系数，建立回归方程，进行影像融合同时对微波数据进行降尺度，将红外遥感数据反演的结果转化为真实的土壤水分值，并且用微波遥感数据对红外数据的空值区域进行填补，从而得到完整的高分辨率土壤水分产品。

采用的具体技术方案如下：

一种结合微波和红外遥感影像的土壤水分反演方法，包括以下步骤：

步骤1：干旱指数PDI适合低植被覆盖区，VSWI适合高植被覆盖区，故基于红外遥感影像，以NDVI指数作为区别植被覆盖程度的阈值，建立综合PDI指数与VSWI指数的综合干旱指数；

步骤2：进行综合干旱指数与微波遥感土壤水分数据的相关性分析，为得到最佳拟合方程，以完成两种数据的融合；

步骤3：根据步骤2得到的最佳拟合方程，融合微波和红外影像的结果。本融合算法融合同时可以进行对微波数据的降尺度，并且用微波数据补全综合干旱指数的空值区，最终得到完整土壤水分产品。

所述步骤1在基于红外遥感影像建立综合干旱指数时，具体包括子步骤：

步骤1.1，计算影像的垂直干旱指数，其计算公式为：

$PDI=\frac{1}{\sqrt{M^2+1}}(R_{red}+{\mathrm{MR}}_{nir}),$

式中，R_red为经过大气校正的红光波段反射率，R_nir为经过大气校正的近红外波段反射率，M为土壤线的斜率；

步骤1.2，计算影像的植被供水指数，其计算公式为：

VSWI＝NDVI/T_s，

式中，T_s为植被的冠层温度，NDVI是归一化植被指数。VSWI是根据植被覆盖状况的变化来进行反演的，所以不适用于低植被覆盖地区。

步骤1.3，PDI适用于低植被覆盖区，VSWI适用于高植被覆盖区，以NDVI为区分植被覆盖区域高低的阈值，综合上两种指数建立综合干旱指数(CDI，compositive droughtindex)。综合干旱指数的定义公式如下：

式中，CDI_i表示在任意像元点i的综合干旱指数，PDI_i表示在任意像元点i的垂直干旱指数，VSWI_i表示在任意像元点i的植被供水指数，PDI_max和PDI_min分别表示垂直干旱指数最大值和最小值，VSWI_max和VSWI_min分别表示植被供水指数最大值和最小值。

所述步骤2在进行综合干旱指数与微波遥感数据相关性分析时，进一步包括子步骤：

步骤2.1：任取50组CDI像元值，求各组像元的平均值。红外遥感数据分辨率较高，此处取的每组包含若干个CDI像元(具体个数根据所采用红外和微波影像的分辨率比例来确定)，取均值是为了让均值代表此组内的CDI像元值；

步骤2.2：取50个与步骤2.1对应位置的微波遥感土壤水分数据的像元值，与步骤2.1中的50个CDI均值建立拟合关系；

步骤2.3：选取相关性最好的拟合模型函数。

所述步骤3在融合微波/红外影像的结果时，进一步包括子步骤：

步骤3.1，微波遥感数据质量控制。微波土壤水分数据会由于各种原因造成某个像元为空值，则该空值像元使用其周围8个像元值的平均值作为该像元值。按行列顺序逐一计算，计算公式如下：

$x_{\mathrm{\θ}}={\mathrm{\Σ}}_{i=0}^7{x}_i/(8-k),$

式中，x_θ表示该空值像元，x_i表示该像元的邻接像元值，k(k∈[0,7]且k∈Z)为周围8个像元中空值像元的个数。

步骤3.2，影像融合同时对微波数据进行降尺度，并把CDI的空值区补全，得到高分辨率的完整土壤水分产品。

所述步骤3.2在融合CDI和微波土壤水分数据时，进一步包括子步骤：

步骤3.2.1：由步骤2所得，微波遥感数据与综合干旱指数的最佳拟合关系为：θ_M＝f(θ_CDI)，此处假设为线性关系，即

θ_M＝a+b×θ_CDI，

式中，θ_M为红外影像网格尺度的一个像元代表的实际土壤含水量，θ_CDI为综合干旱指数上的一个像元值，a,b为系数；

步骤3.2.2：取红外遥感影像网格上相邻i个像元的土壤含水量的平均值等于对应的微波遥感影像网格上一个像元的土壤含水量的值：

${\mathrm{\θ}}_S=({\mathrm{\θ}}_{M_1}+{\mathrm{\θ}}_{M_2}+\mathrm{.......}+{\mathrm{\θ}}_{M_i})/i=f\left({\mathrm{\θ}}_{M_i}\right),$

式中，θ_S为微波遥感网格尺度上的一个像元的土壤含水量，为红外遥感影像网格尺度上一个像元的土壤含水量；

步骤3.2.3：把步骤3.2.2中的公式代入步骤3.2.1中，有

${\mathrm{\θ}}_S=f\left({\mathrm{\θ}}_{M_i}\right)=f\left(f\right({\mathrm{\θ}}_{{\mathrm{CDI}}_i}\left)\right)=a+b\×f\left({\mathrm{\θ}}_{{\mathrm{CDI}}_i}\right),$

式中，表示CDI指数结果上一个像元的值；

则对于每景影像可列出多个式子：

$\left\{\begin{array}{c}{\mathrm{\θ}}_{S_1}=a+b\×f{\left({\mathrm{\θ}}_{{\mathrm{CDI}}_i}\right)}_1\\ {\mathrm{\θ}}_{S_2}=a+b\×f{\left({\mathrm{\θ}}_{{\mathrm{CDI}}_i}\right)}_2\\ .\\ .\\ .\\ {\mathrm{\θ}}_{S_x}=a+b\×f{\left({\mathrm{\θ}}_{{\mathrm{CDI}}_i}\right)}_x\\ .\\ .\\ .\end{array}\right.,$

上式除了参数a、b都是已知量，通过最小二乘法解算出拟合关系最好的a、b，带入步骤3.2.1，求得每个θ_M的值。

至此，即得到了在红外遥感影像分辨率尺度上每个像元的实际土壤含水量。

本发明具有如下特点和有益效果：

1、本发明建立了综合植被指数。根据地表覆盖状况分成两种情形：稀疏植被覆盖下的土壤和茂密植被覆盖下的土壤。垂直干旱指数PDI比较适用于稀疏植被覆盖下的区域，而植被供水指数VSWI适用于茂密植被覆盖下的地区，因此可以基于研究区域的NDVI值来划分研究区域，从而对复杂地表覆盖加以简化，选择适合不同覆盖类型的反演模型来反演土壤水分，以提高整个区域土壤水分遥感反演精度。

2、本发明将微波遥感数据与红外遥感数据结合，在影像融合的同时对微波遥感数据进行降尺度处理，即弥补了微波遥感数据分辨率低的缺点，同时又弥补了红外遥感数据因云层影响有空缺的情况，并且利用微波遥感数据与红外遥感数据的相关系数，将红外遥感数据反演的CDI指数结果转化为真实的土壤水分值，最终得到完整的土壤水分产品，这种方法能提高大范围区域的陆表土壤水分监测效率，便于土壤水分产品的业务化生产。

附图说明

图1为本发明的方法流程示意图；

图2为降尺度原理示意图；

图3为微波影像质量控制原理示意图。

具体实施方式

下面详细描述本发明的具体实施例。应当注意，这里描述的实施例只用于举例说明，并不用于限制本发明：

1、计算红外遥感影像的综合干旱指数，此例选用MODIS影像，分辨率为1km。

计算影像的垂直干旱指数，式中，R_red为经过大气校正的红光波段反射率；R_nir为经过大气校正的近红外波段反射率；M为土壤线的斜率；

计算影像的植被供水指数，VSWI＝NDVI/T_s，式中，T_s为植被的冠层温度，NDVI是归一化植被指数。

PDI适用于低植被覆盖区，VSWI适用于高植被覆盖区，以NDVI为阈值，综合上两种指数建立综合干旱指数(CDI，compositive drought index)：式中，CDI_i表示在任意像元点i的综合干旱指数，PDI_i表示在任意像元点i的垂直干旱指数，VSWI_i表示在任意像元点i的植被供水指数，PDI_max和PDI_min分别表示垂直干旱指数最大值和最小值，VSWI_max和VSWI_min分别表示植被供水指数最大值和最小值。

2、综合干旱指数与微波遥感数据相关性分析，此例微波影像选用SMOS影像土壤水分数据。

MODIS的空间分辨率为1km，而SMOS的空间分辨率为25km，即一个SMOS像元的网格大小约为25*25个MODIS像元的网格大小，在此，取25*25个MODIS网格像元值的平均值作为一个SMOS相对应位置上的像元值,如图2。

任取30组25*25的CDI(综合干旱指数)像元值，计算其平均值，然后取对应的SMOS尺度上的土壤水分像元值，建立拟合关系，选取相关系数最好的拟合模型函数，此例为线性模型。

3、结合综合干旱指数对SMOS土壤水分数据进行降尺度，得到完整土壤水分产品。

先进行质量控制，SMOS土壤水分数据会由于各种原因造成某个像元为空值，用该空值像元周围8个像元值的平均值为该像元值。按行列顺序逐一计算，某空值像元k为邻接像元中空值像元的个数，k∈[0,7]且k∈Z。如图3。

影像融合同时进行降尺度：

综合植被指数与SMOS土壤水分数据拟合最好的为线性关系，即θ_M＝a+b×θ_CDI①，其中θ_M为MODIS网格尺度的一个像元代表的实际土壤含水量，θ_CDI为综合干旱指数的一个像元值，a,b为系数；

取MODIS网格上相邻i个像元的土壤含水量的平均值为SMOS影像网格上一个像元的土壤含水量的值： 其中θS为SMOS影像网格尺度上的一个像元的土壤含水量，为MODIS影像网格尺度上一个像元的土壤含水量；

把②公式代入①中，有则对于每景影像可列出若干个式子：除了参数a、b都是已知量，通过最小二乘法解算出拟合关系最好的a、b，进而带入公式①，求得每个θ_M的值。

最终得到了在MODIS影像网格尺度上每个像元的实际土壤含水量。即得到了完整的1km分辨率土壤水分产品。

Claims (5)

1.一种结合微波和红外遥感影像的土壤水分反演方法，其特征在于:包括以下步骤：

步骤1：基于红外遥感影像建立以NDVI为阈值的综合PDI指数与VSWI指数的综合干旱指数；

步骤2：进行综合干旱指数与微波遥感土壤水分数据的相关性分析；

步骤3：融合微波和红外影像的结果，得到完整土壤水分产品。

2.如权利要求书1所述的一种结合微波和红外遥感影像的土壤水分反演方法，其特征在于：所述步骤1在基于红外遥感影像建立综合干旱指数，具体包括以下子步骤：

步骤1.1：计算红外影像的垂直干旱指数，其计算公式为：

$PDI=\frac{1}{\sqrt{M^2+1}}(R_{red}+{\mathrm{MR}}_{nir}),$

式中，R_red为经过大气校正的红光波段反射率，R_nir为经过大气校正的近红外波段反射率，M为土壤线的斜率；

步骤1.2:计算红外影像的植被供水指数，其计算公式为：

VSWI＝NDVI/T_s，

式中，T_s为植被的冠层温度，NDVI是归一化植被指数；

步骤1.3:PDI适用于低植被覆盖区，VSWI适用于高植被覆盖区，以NDVI为阈值，综合此两种指数建立综合干旱指数CDI，综合干旱指数的定义公式如下：

${\mathrm{CDI}}_i=\left\{\begin{array}{cc}1-\frac{{\mathrm{PDI}}_i-{\mathrm{PDI}}_{\min }}{{\mathrm{PDI}}_{\max }-{\mathrm{PDI}}_{\min }}& \left({\mathrm{NDVI}}_i\≤0.3\right)\\ \frac{{\mathrm{VSWI}}_i-{\mathrm{VSWI}}_{\min }}{{\mathrm{VSWI}}_{\max }-{\mathrm{VSWI}}_{\min }}& \left({\mathrm{NDVI}}_i>0.3\right)\end{array}\right.,$

式中，CDI_i表示在任意像元点i的综合干旱指数，PDI_i表示在任意像元点i的垂直干旱指数，VSWI_i表示在任意像元点i的植被供水指数，PDI_max和PDI_min分别表示垂直干旱指数最大值和最小值，VSWI_max和VSWI_min分别表示植被供水指数最大值和最小值。

3.如权利要求书2所述的一种结合微波和红外遥感影像的土壤水分反演方法，其特征在于：所述步骤2在进行综合干旱指数与微波遥感数据相关性分析，进一步包括以下子步骤：

步骤2.1：任取50组综合干旱指数CDI像元值；

步骤2.2：取对应的微波遥感土壤水分数据的像元值，建立拟合关系；

步骤2.3：选取相关性最好的拟合模型函数。

4.如权利要求书3所述的一种结合微波和红外遥感影像的土壤水分反演方法，其特征在于：所述步骤3在融合微波和红外影像的结果，进一步包括以下子步骤：

步骤3.1：微波遥感数据质量控制：若微波土壤水分数据中的个像元为空值，则该空值像元使用其周围8个像元值的平均值作为该像元值；计算公式如下：

$x_{\mathrm{\θ}}={\mathrm{\Σ}}_{i=0}^7{x}_i/(8-k),$

式中，x_θ表示该空值像元，x_i表示该像元的邻接像元值，k为周围8个像元中空值的个数,k∈[0,7]且k∈Z；

步骤3.2：影像融合同时进行降尺度，以获得高分辨率的完整土壤水分产品。

5.如权利要求书4所述的一种结合微波和红外遥感影像的土壤水分反演方法，其特征在于：所述步骤3.2在影像融合同时进行降尺度，进一步包括以下子步骤：

步骤3.2.1：由步骤2所得，微波遥感数据与综合干旱指数的最佳拟合关系为：θ_M＝f(θ_CDI)，此处假设为线性关系，即

θ_M＝a+b×θ_CDI，

式中，θ_M为红外影像网格尺度上的一个像元代表的实际土壤含水量，θ_CDI为综合干旱指数上一个像元的值，a,b为系数；

步骤3.2.2：取红外遥感影像网格上相邻i个像元的土壤含水量的平均值为微波遥感影像网格上一个像元的土壤含水量的值：

${\mathrm{\θ}}_S=({\mathrm{\θ}}_{M_1}+{\mathrm{\θ}}_{M_2}+\mathrm{.......}+{\mathrm{\θ}}_{M_i})/i=f\left({\mathrm{\θ}}_{M_i}\right),$

式中，θ_s表示微波遥感网格尺度上的一个像元的土壤含水量，表示红外遥感影像网格尺度上一个像元的土壤含水量；

步骤3.2.3：把步骤3.2.2中的公式代入步骤3.2.1中，有

${\mathrm{\θ}}_S=f\left({\mathrm{\θ}}_{M_i}\right)=f\left(f\right({\mathrm{\θ}}_{{\mathrm{CDI}}_i}\left)\right)=a+b\×f\left({\mathrm{\θ}}_{{\mathrm{CDI}}_i}\right),$

式中，表示CDI指数结果上一个像元的值；

则对于每景影像可列出多个式子：

$\left\{\begin{array}{c}{\mathrm{\θ}}_{S_1}=a+b\×f{\left({\mathrm{\θ}}_{{\mathrm{CDI}}_i}\right)}_1\\ {\mathrm{\θ}}_{S_2}=a+b\×f{\left({\mathrm{\θ}}_{{\mathrm{CDI}}_i}\right)}_2\\ \begin{array}{c}.\\ .\\ .\end{array}\\ {\mathrm{\θ}}_{S_x}=a+b\×f{\left({\mathrm{\θ}}_{{\mathrm{CDI}}_i}\right)}_x\\ \begin{array}{c}.\\ .\\ .\end{array}\end{array}\right.,$

该式中除了参数a、b都是已知量，通过最小二乘法解算出拟合关系最好的a、b，带入步骤3.2.1，求得每个θ_M的值。