CN102034337B - 草原雪灾遥感监测与灾情评估系统及方法 - Google Patents

Abstract

本发明公开了一种草原雪灾遥感监测与灾情评估系统及方法，其中该系统包括：预处理模块，用于读取并处理MODIS L1B数据、AMSR-E L1B数据，得到预处理结果；积雪区域融合模块，用于根据预处理结果进行数据融合处理，得到单日总积雪区域；积雪深度估算模块，用于得到单日积雪深度；积雪覆盖率估算模块，用于得到单日总积雪覆盖率；草群高度估算模块，用于得到单日草群高度；积雪合成模块，用于得到积雪持续日数、平均积雪覆盖率、平均积雪深度、平均草群高度；雪灾等级评价模块，用于对雪灾等级进行评价。本发明能够对草原积雪情况进行全天候动态监测，并能够根据草原雪灾等级做出预警。

Description

草原雪灾遥感监测与灾情评估系统及方法

技术领域

本发明涉及利用MODIS(Moderate Resolution Imaging Spectrometer：中分辨率成像光谱仪)和AMSR-E(Advanced Microwave Scanning Radiometer-EOS：地球观测系统-高级微波扫描辐射计)L1B数据进行草原雪灾遥感监测与灾情评估的技术，特别是涉及利用MODIS和AMSR-E L1B数据融合对草原进行全天候积雪遥感监测与灾情评估的系统及其方法。

背景技术

中国草原面积占国土总面积的40％左右，约4亿hm²，其中90％以上又集中分布在新疆、内蒙古、西藏和青海四省区，这些地区纬度高、海拔高，气候寒冷、干燥，自然灾害频繁，其中尤以雪灾对牲畜越冬危害最重，每次雪灾发生，夺去牲畜轻则几十万头，重则几百万头。长期以来，雪灾对我国草原地区农牧业生产的持续、稳定发展造成了极其严重的危害。鉴于草原雪灾对我国草原畜牧业发展的严重影响，有必要全面、快速地监测雪灾成灾区域、准确评估雪灾对草原地区造成的影响，以便最大限度地减少雪灾造成的损失。

目前基于遥感的方法进行草原雪灾的监测与灾情评估，主要包括积雪识别、积雪深度反演、草群高度反演、雪灾等级评估等方法。

积雪识别可分为光学传感器积雪识别和被动微波传感器积雪识别。对于光学传感器数据积雪识别，一般采用归一化差分积雪指数(Normalized DifferenceSnow Index，NDSI)算法，并结合多阈值的决策树方法。这种方法主要是依据积雪在可见光范围内较高的反射率、在中红外区极低的反射率这种特征，采用比值和差值的积雪指数形式来识别积雪，参照现有技术文献1(Dorothy K.Hall，George A.Riggs and Vincent V.Salomonsont.Development of Methods forMapping Global Snow Cover Using Moderate Resolution ImagingSpectroradiometer Data.Remote Sensing of Environment，1995，54：127-140.)。对于被动微波传感器积雪识别，一般采用多通道亮温差值法、多阈值决策树法等。这种方法的依据是积雪对不同频率的微波辐射能量散射与吸收的不同。由于积雪颗粒对高频能量的散射能力更强，这样造成低频与高频通道的亮温差为正值，参照现有技术文献2(Norman C.Grody，Alan N.Basist.GlobalIdentification of Snowcover Using SSM/I Measurements.IEEE Transactions onGeoscience and Remote Sensing，1996，34(1)：237-249.)。以上两种方法各有各的优点，光学传感器数据空间分辨率较高，同时识别精度也较高，但其受天气(主要是云)的影响比较严重；被动微波传感器数据积雪识别精度比较低，但其具有全天候工作的能力。

对于积雪深度遥感反演，目前采用的方法可分为物理模型法和统计模型法。物理模型法是以理论为基础，研究微波在积雪中的传播规律，进而得到观测亮温与积雪参数之间的关系。而统计模型法通常使用多波段亮温观测值与实测积雪深度值，采用统计模式来求出最佳的反演因子。这种方法主要用于在不了解积雪物理特征的情况下，直接建立辐射亮温与积雪深度之间的经验模型。物理模型法估算积雪深度精度较高，但是其模型复杂，同时模型中需要输入很多已知的积雪参数，因此在实际应用中，更多地使用相对简单的统计模型方法，这里尤以双通道差值法应用最为广泛，参照现有技术文献3(A.T.C Chang，J.LFoster and D.K Hall.Nimbus 7 SMMR derived global snow cover patterns.Annalsof Glaciology，1987，9：39-44.)。

对于草群高度的遥感估算，目前主要使用的方法是基辐射传输模型的单通道算法和基于植被指数模型的统计算法。单通道算法通过观测不同草群高度下的二向性反射分布函数(Bi-directional Reflectance Distribution Function，BRDF)，进而通过使用辐射传输模型，模拟可见光、近红外通道反射率与BRDF、以及草群高度之间的关系，最后建立估算模型，参照非现有技术文献4(KONDA ASAKO，YAMAMOTO HIROKAZU，KAJIWARA KOJI，HONDAYOSHIAKI.A Study on Estimation of Grass Height based on BRDF Model usingSatellite Data.Journal of the Japan Society of Photogrammetry and RemoteSensing，2001，40(6)：15-24.)。使用这种方法时，由于需要实测的BRDF，因此在实际应用中会受到限制。当采用植被指数模型算法时，这相对简单，这时需要针对不同的草地类型、草群生长期建立模型。

对于草原雪灾灾情等级评估，目前主要是基于层次分析法，使用最大积雪深度、最低气温、降雪时间、积雪日数、低于多年月平均气温的延续天数以及降雪过程的降水总量等相关因子，并综合考虑雪情、草情、畜情和气象等多方面的因子，进而构建草原雪灾等级以及灾情评价指标体系。现有技术依据上述指标，提出了牧区雪灾等级国家标准，参照国家标准文献1(国家标准文献1：GB/T 20482-2006，牧区雪灾等级)。

发明内容

本发明所要解决的技术问题在于提供一种草原雪灾遥感监测与灾情评估系统及方法，用于对草原积雪情况进行全天候动态监测，并能够根据草原雪灾等级做出预警。

为了实现上述目的，本发明提供了一种草原雪灾遥感监测与灾情评估系统，其特征在于，包括：

预处理模块，用于读取并处理MODIS L1B数据、AMSR-E L1B数据，得到MODIS反射率、MODIS亮度温度数据、AMSR-E亮度温度数据；

积雪区域融合模块，连接所述预处理模块，用于根据所述MODIS反射率、所述MODIS亮度温度数据进行MODIS积雪像元识别，根据所述AMSR-E亮度温度数据进行AMSR-E积雪像元识别，并对得到的识别结果进行数据融合处理，得到单日积雪区域；

积雪深度估算模块，连接所述预处理模块，用于根据所述AMSR-E亮度温度数据进行AMSR-E积雪深度估算，得到单日积雪深度；

积雪覆盖率估算模块，连接所述预处理模块，用于根据所述MODIS反射率、所述MODIS亮度温度数据进行MODIS积雪覆盖率估算，得到单日积雪覆盖率；

草群高度估算模块，连接所述预处理模块，用于根据所述MODIS反射率、所述MODIS亮度温度数据进行MODIS草群高度估算，得到单日草群高度；

积雪合成模块，连接所述积雪区域融合模块、所述积雪深度估算模块、所述积雪覆盖率估算模块、所述草群高度估算模块，用于根据所述单日积雪区域、所述单日积雪覆盖率、所述单日积雪深度、所述单日草群高度，对多日的总积雪区域、总积雪覆盖率、积雪深度、草群高度进行统计合成，得到积雪持续日数、平均积雪覆盖率、平均积雪深度、平均草群高度；

雪灾等级评价模块，连接所述积雪合成模块，用于根据所述积雪持续日数、所述平均积雪覆盖率、所述平均积雪深度、所述平均草群高度，对草原区域的雪灾等级进行评价。

为了实现上述目的，本发明还提供了应用上述系统的草原雪灾遥感监测与灾情评估方法，其特征在于，包括：

数据预处理步骤，读取并处理MODIS L1B数据、AMSR-E L1B数据，得到MODIS反射率、MODIS亮度温度数据、AMSR-E亮度温度数据；

积雪区域融合步骤，根据所述MODIS反射率、所述MODIS亮度温度数据进行MODIS积雪像元识别，根据所述AMSR-E亮度温度数据进行AMSR-E积雪像元识别，并对得到的识别结果进行数据融合处理，得到单日积雪区域；

积雪深度估算步骤，根据所述AMSR-E亮度温度数据进行AMSR-E积雪深度估算，得到单日积雪深度；

积雪覆盖率估算步骤，根据所述MODIS反射率、所述MODIS亮度温度数据进行MODIS积雪覆盖率估算，得到单日积雪覆盖率；

草群高度估算步骤，根据所述MODIS反射率、所述MODIS亮度温度数据进行MODIS草群高度估算，得到单日草群高度；

积雪合成步骤，根据所述单日积雪区域、所述单日积雪覆盖率、所述单日积雪深度、所述单日草群高度，对多日的积雪区域、积雪覆盖率、积雪深度、草群高度进行统计合成，得到积雪持续日数、平均积雪覆盖率、平均积雪深度、平均草群高度；

雪灾等级评价步骤，根据所述积雪持续日数、所述平均积雪覆盖率、所述平均积雪深度、所述平均草群高度，对草原区域的雪灾等级进行评价。

所述的草原雪灾遥感监测与灾情评估方法，其中，所述预处理步骤中，进一步包括：

以如下公式得到所述MODIS反射率：

R＝scale(SI-offset)

式中，R为MODIS反射率，scale，offset为比例系数，可从HDF属性数据中获得，SI为MODIS L1B数据中单个像元的灰度值；

以如下公式得到所述MODIS亮度温度数据：

$\mathrm{Tb}=\frac{\mathrm{hc}}{\mathrm{\λk}}\·\frac{1}{\frac{2{\mathrm{hc}}^2}{L{\mathrm{\λ}}^5}+1}$

式中，L为黑体的辐射亮度(Wm^-2μm^-1sr^-1)；h是Plank常数，其值为6.63×10^-34Js；k是Boltzmann常数，其值为1.38×10^-23JK^-1；c为光速(3×10⁸ms^-1)；λ为波长(μm)；Tb为MODIS亮度温度(K)；

以如下公式得到所述AMSR-E亮度温度数据：

Tb＝scale×SI+offset

式中，Tb为AMSR-E亮度温度(K)，scale，offset为比例系数，可从HDF属性数据中获得，SI为AMSR-E L1B数据中单个像元的灰度值。

所述的草原雪灾遥感监测与灾情评估方法，其中，所述积雪区域融合步骤中，进一步包括：

根据所述MODIS反射率、所述MODIS亮度温度数据进行云识别，获取无云区域MODIS L1B数据；

根据所述无云区域MODIS L1B数据进行MODIS积雪像元识别，得到MODIS识别积雪区域；

根据所述AMSR-E亮度温度数据进行AMSR-E积雪像元识别，得到AMSR-E识别积雪区域；

对所述MODIS识别积雪区域、所述AMSR-E识别积雪区域进行多源数据融合，得到所述单日积雪区域。

所述的草原雪灾遥感监测与灾情评估方法，其中，所述云识别步骤中，进一步包括：

通过判断所述MODIS L1B数据中每个像元是否满足云识别条件获取所述无云区域MODIS L1B数据，云识别条件的表达式为：

Tb₃₅＜220K

Tb₂₇＜226K

Tb₃₁-Tb₂₂＜-12K

ρ₂₆＞0.035

ρ₁＞0.18

式中，ρ₁，ρ₂₆分别为MODIS L1B数据第1、26波段的反射率；Tb₂₂，Tb₂₇，T₃₁，T₃₅分别为MODIS L1B数据第22、27、31、35波段的亮度温度；

所述MODIS积雪像元识别步骤中，进一步包括：

通过判断所述无云区域MODIS L1B数据中每个像元是否满足积雪像元识别条件得到所述MODIS识别积雪区域，该积雪像元识别条件的表达式为：

$\left\{\begin{array}{c}\mathrm{NDSI}\≥0.4\\ {\mathrm{\ρ}}_2>0.11\\ {\mathrm{\ρ}}_4\≥0.1\end{array}\right.$

式中，ρ₂、ρ₄分别为MODIS L1B数据第2、4波段的反射率，NDSI为每个像元的归一化差分积雪指数；

以MODIS L1B数据的第4、6波段来构建NDSI，表达式为：

$\mathrm{NDSI}=\frac{{\mathrm{\ρ}}_4-{\mathrm{\ρ}}_6}{{\mathrm{\ρ}}_4+{\mathrm{\ρ}}_6}$

式中，ρ₄、ρ₆分别为MODIS L1B数据第4、6波段的反射率。

所述的草原雪灾遥感监测与灾情评估方法，其中，所述积雪深度估算步骤中，进一步包括：

根据所述AMSR-E亮度温度数据进行AMSR-E积雪像元识别，得到AMSR-E识别积雪区域；

根据所述AMSR-E识别积雪区域进行AMSR-E积雪深度计算，得到所述单日积雪深度。

所述的草原雪灾遥感监测与灾情评估方法，其中，所述AMSR-E积雪深度计算步骤中，进一步包括：

通过构建积雪深度反演模型得到所述单日积雪深度，该积雪深度反演模型的表达式为：

SD＝a×(Tb_18.7H-Tb_36.5H)+b

式中，SD为单日积雪深度，a、b系数的取值分别为1.59cm/K和0，Tb_18.7H和Tb_36.5H分别表示AMSR-E L1B数据的18.7GHz通道的水平极化亮度温度、36.5GHz通道的水平极化亮度温度。

所述的草原雪灾遥感监测与灾情评估方法，其中，所述积雪覆盖率估算步骤中，进一步包括：

根据所述MODIS反射率、所述MODIS亮度温度数据进行云识别，获取无云区域MODIS L1B数据；

根据所述无云区域MODIS L1B数据进行MODIS积雪覆盖率计算，得到所述单日积雪覆盖率。

所述的草原雪灾遥感监测与灾情评估方法，其中，所述MODIS积雪覆盖率计算步骤中，进一步包括：

通过构建以NDSI为变量的多项式拟合模型得到所述单日积雪覆盖率，该多项式拟合模型的表达式为：

SF＝a×NDSI+b或

SF＝c×NDSI²+d×NDSI+e

式中，SF为单日积雪覆盖率，a、b系数的取值分别为1.21和0.06，c、d、e系数取值分别为0.26、0.37和0.18。

所述的草原雪灾遥感监测与灾情评估方法，其中，所述草群高度估算步骤中，进一步包括：

根据所述MODIS反射率、所述MODIS亮度温度数据进行云识别，获取无云区域MODIS L1B数据；

根据所述无云区域MODIS L1B数据进行MODIS草群高度计算，得到所述单日草群高度。

所述的草原雪灾遥感监测与灾情评估方法，其中，所述MODIS草群高度计算步骤中，进一步包括：

通过NDVI统计模型法得到所述单日草群高度，NDVI统计模型的表达式为如下任一种：

GH＝a×NDVI+b

式中，a、b系数的一般取值为0.7371、-4.6131，GH为单日草群高度；

GH＝a×NDVI²+b×NDVI+c

式中，a、b、c系数的一般取值为0.0062、0.1966、4.9485；

GH＝a×NDVI³+b×NDVI²+c×NDVI+d

式中，a、b、c、d系数的一般取值为0.00008、-0.0046、0.6183、0.4100；

GH＝a×e^b×NDVI

式中，a、b系数的一般取值为5.4666、0.0318；

GH＝a×NDVI^b

式中，a、b系数的一般取值为0.3402、1.1371；

以MODIS L1B数据的第1、2通道来构建NDVI，表达式为：

$\mathrm{NDVI}=\frac{{\mathrm{\ρ}}_2-{\mathrm{\ρ}}_1}{{\mathrm{\ρ}}_2+{\mathrm{\ρ}}_1}$

式中，ρ₁、ρ₂分别为MODIS L1B数据的第1、2通道的反射率。

与现有技术相比，本发明的有益技术效果在于：

本发明提供了一种基于MODIS L1B和AMSR-E L1B数据融合进行全天候积雪监测与灾情评估的方法，该方法能够快速准确地提取积雪覆盖面积，积雪深度、草群高度，且在计算过程中，根据MODIS的云识别结果，将无云区的MODIS监测与有云区的AMSR-E监测结果进行融合，避免了云层对MODISL1B数据的影响，同时也避免了仅使用AMSR-E L1B数据时精度较差的缺陷。

以下结合附图和具体实施例对本发明进行详细描述，但不作为对本发明的限定。

附图说明

图1为本发明的草原雪灾遥感监测与灾情评估系统结构图；

图1a为本发明的积雪区域融合模块的结构图；

图1b为本发明的积雪深度估算模块的结构图；

图1c为本发明的积雪覆盖率估算模块的结构图；

图1d为本发明的草群高度估算模块的结构图；

图2为本发明的草原雪灾遥感监测与灾情评估方法流程图；

图3为本发明的AMSR-E积雪识别方法流程图；

图4为本发明的获取单日总积雪区域的方法流程图；

图5为本发明的获取单日积雪深度的方法流程图；

图6为本发明的获取单日总积雪覆盖率的方法流程图；

图7为本发明的获取单日草群高度的方法流程图。

具体实施方式

下面结合附图和具体实施方式对本发明的技术方案作进一步更详细的描述。

如图1所示，为本发明的草原雪灾遥感监测与灾情评估系统结构图。并结合图1a-图1d所示。该系统100具体包括：预处理模块10、积雪区域融合模块21、积雪深度估算模块22、积雪覆盖率估算模块23、草群高度估算模块24、积雪合成模块30、雪灾等级评价模块40。

预处理模块10，用于从原始L1B级HDF属性数据中，读取相关波段数据，并对该相关波段数据进行辐射标定计算和空间参考变换。该相关波段数据包括MODIS L1B数据、AMSR-E L1B数据。

对于MODIS L1B数据，进行预处理，计算反射率或光谱辐射亮度值，如下：

R＝scale(SI-offset)                (1)

上式(1)中，R为像元反射率或光谱辐射亮度值，scale，offset为比例系数，可从HDF属性数据中获得，SI为像元的灰度值。此外，基于Plank方程，可由光谱辐射亮度值计算目标亮度温度，如下：

$\mathrm{Tb}=\frac{\mathrm{hc}}{\mathrm{\λk}}\·\frac{1}{\frac{2{\mathrm{hc}}^2}{L{\mathrm{\λ}}^5}+1}---\left(2\right)$

上式(2)中，L为黑体的辐射亮度(Wm^-2μm^-1sr^-1)；h是Plank常数，其值为6.63×10^-34Js；k是Boltzmann常数，其值为1.38×10^-23JK^-1；c为光速(3×10⁸ms^-1)；λ为波长(μm)；Tb为亮度温度(K)。

对于AMSR-E L1B数据，进行预处理，计算每个像元的亮度温度，如下：

Tb＝scale×SI+offset               (3)

上式(3)中，Tb为波段亮度温度(K)，scale，offset为比例系数，可从HDF属性数据中获得，SI为像元的灰度值。

积雪区域融合模块21，连接预处理模块10，用于对无云区域的MODIS积雪识别结果与有云区域的AMSR-E积雪识别结果进行多源数据融合，实现完整的全天候的积雪监测，得到单日积雪区域，该单日积雪区域可用于积雪持续日数的统计。该积雪区域融合模块21通过如下方法进行融合：

首先，数据重采样，将AMSR-E 10km空间分辨率的积雪范围图重采样为与MODIS相同的1km空间分辨率；

其次，以像元为单位进行运算，此时参与运算的数据包括MODIS云识别结果、MODIS积雪像元识别结果、AMSR-E积雪像元识别结果，运算方法为：MODIS积雪像元识别结果×(MODIS云识别结果＝＝0)+AMSR-E积雪像元识别结果×(MODIS云识别结果＝＝1)。

以上三种参与运算的数据均为二值数据，也就是像元取值只有两种：0和1。对于MODIS和AMSR-E积雪识别结果，像元值为1表示该像元为积雪像元，像元值为0时表示该像元为未积雪像元；同理，MODIS云识别结果中，像元值为1表示该像元为云，像元值为0则表示该像元未被云覆盖。

积雪深度估算模块22，连接预处理模块10，用于根据预处理模块10预处理后的结果进行积雪深度估算，得到单日积雪深度。

积雪覆盖率估算模块23，连接预处理模块10，用于根据预处理模块10预处理后的结果进行积雪覆盖率估算，得到单日积雪覆盖率。

草群高度估算模块24，连接预处理模块10，用于根据预处理模块10预处理后的结果进行草群高度，得到单日草群高度。

积雪合成模块30，连接积雪区域融合模块21、积雪深度估算模块22、积雪覆盖率估算模块23、草群高度估算模块24，用于根据单日的积雪和草群高度监测结果，对多日的积雪覆盖和草群高度进行统计合成，得到积雪持续日数、平均积雪深度、平均积雪覆盖率和平均草群高度。

其中，积雪持续日数为连续积雪的时间，其合成算法包括：使用当日的遥感数据识别积雪像元；提取前一日积雪持续日数图中的积雪日数，如果某像元当日积雪，则新积雪持续日数图中该像元处的积雪日数等于前一日积雪日数加一，否则该像元处的积雪日数为零(重新开始统计)。

上述平均积雪深度的合成算法包括：根据合成周期内的单日积雪深度图，以单个像元为运算单位，取每个像元的单日积雪深度值进行平均计算，得到多日平均积雪深度；

上述平均积雪覆盖率的合成算法包括：根据合成周期内的单日积雪覆盖率图，以单个像元为运算单位，取每个像元的单日积雪覆盖率值进行平均计算，得到多日平均积雪覆盖率；

上述平均草群高度的合成算法包括：根据合成周期内的单日草群高度图，以单个像元为运算单位，取每个像元的单日草群高度值进行平均计算，得到多日平均草群高度；

雪灾等级评价模块40，连接积雪合成模块30，用于根据积雪合成模块30得到的平均积雪深度、平均草群高度计算牧区雪灾等级评价中的因子积雪厚度/草群高度，根据积雪合成模块30得到的积雪持续日数计算牧区雪灾等级评价中的因子积雪持续时间(日)，根据积雪合成模块30得到的平均积雪覆盖率计算牧区雪灾等级评价中的因子积雪面积/土地面积，并以牧区雪灾等级国家标准(表3)为基础，对草原区域的雪灾等级进行评价。

进一步地，积雪区域融合模块21又包括：

MODIS云识别模块211，用于通过判断遥感影像中每个像元是否满足云识别条件，获取无云区域MODIS L1B数据、有云区域MODIS L1B数据。其判断表达式为：

Tb₃₅＜220K

Tb₂₇＜226K

Tb₃₁-Tb₂₂＜-12K                (4)

ρ₂₆＞0.035

ρ₁＞0.18

上式(4)中，ρ₁，ρ₂₆分别为MODIS L1B数据的第1、26波段的反射率；Tb₂₂，Tb₂₇，Tb₃₁，Tb₃₅分别为MODIS L1B数据的第22、27、31、35波段的亮度温度。

MODIS积雪像元识别模块212，连接MODIS云识别模块211，用于根据无云区域MODIS L1B数据进行MODIS积雪像元识别，获取MODIS积雪识别结果：无云区域的MODIS积雪识别结果、有云区域的MODIS积雪识别结果。

MODIS积雪像元识别是利用MODIS L1B数据的第4、6、7波段来构建归一化差分积雪指数NDSI的步骤，其表达式为：

$\mathrm{NDSI}=\frac{{\mathrm{\ρ}}_4-{\mathrm{\ρ}}_6}{{\mathrm{\ρ}}_4+{\mathrm{\ρ}}_6}---\left(5\right)$

上式(5)中，ρ₄、ρ₆分别为MODIS L1B数据第4、6波段的反射率。由于Aqua MODIS的第6通道的数据有问题，因此在计算Aqua MODIS NDSI时，采用第7通道反射率ρ₇来代替ρ₆进行计算。此外，还包括用于判断每个像元的积雪指数NDSI，绿、近红外波段反射率，是否满足积雪像元的识别条件的步骤：

$\left\{\begin{array}{c}\mathrm{NDSI}\≥0.4\\ {\mathrm{\ρ}}_2>0.11\\ {\mathrm{\ρ}}_4\≥0.1\end{array}\right.---\left(6\right)$

上式(6)中，ρ₂、ρ₄分别为MODIS L1B数据第2、4波段的反射率。

AMSR-E积雪像元识别模块213，用于通过判断AMSR-E L1B数据中每个像元是否满足积雪识别条件，得到AMSR-E积雪识别结果：AMSR-E识别积雪区域、AMSR-E识别未积雪区域。

这里，针对AMSR-E L1B数据，定义散射指数(scatter index)具有如下的形式；

scat＝max(Tb_18.7V-Tb_36.5V-3，Tb_23.8V-Tb_89V-3，Tb_36.5V-Tb_89V-1)        (7)

上式(7)中，scat为散射指数，max为取最大值函数，Tb为亮度温度，Tb右下标中的数字为通道频率，字母为极化方式(V表示垂直极化、H表示水平极化)，如Tb_18.7V表示18.7GHz通道垂直极化亮度温度。该步骤中，首先根据散射指数确定散射体，散射体是根据散射指数判断出来的像元，然后再根据各条件式依次识别出降雨、寒漠、冻土等地表，最后识别出积雪区。

数据融合模块214，连接MODIS积雪像元识别模块212、AMSR-E积雪像元识别模块213，用于对无云区域的MODIS积雪识别结果与有云区域的AMSR-E积雪识别结果进行多源数据融合，实现完整的全天候的积雪监测，得到单日积雪区域，并可用于积雪持续日数的统计。

进一步地，积雪深度估算模块22又包括：

AMSR-E积雪像元识别模块213，用于根据预处理后得到的AMSR-E亮度温度数据进行AMSR-E积雪像元识别，得到AMSR-E积雪识别结果：AMSR-E识别积雪区域、AMSR-E识别未积雪区域；

积雪深度计算模块215，连接AMSR-E积雪像元识别模块213，用于根据AMSR-E识别积雪区域进行AMSR-E积雪深度计算，得到单日积雪深度。其是通过构建积雪深度反演模型来估算的，采用AMSR-E 18.7GHz和36.5GHz通道水平极化亮度温度差(散射指数)进行积雪深度估算，并在此之上构建了积雪深度反演的线性模型，其表达式为：

SD＝a×(Tb_18.7H-Tb_36.5H)+b              (8)

上式(8)中，a、b两个系数的取值分别为1.59cm/K和0，Tb_18.7H和Tb_36.5H分别表示18.7GHz水平极化亮度温度、36.5GHz水平极化亮度温度。

此外，还可以根据草原区域的地域特征，采用不同的模型系数来估算积雪深度。

进一步地，积雪覆盖率估算模块23又包括：

MODIS云识别模块211，用于根据预处理后得到的MODIS反射率与亮度温度数据进行MODIS云识别，得到无云区域MODIS L1B数据、有云区域MODIS L1B数据；

积雪覆盖率计算模块216，连接MODIS云识别模块211，用于根据无云区域MODIS L1B数据进行MODIS积雪覆盖率计算，得到单日积雪覆盖率。其是采用NDSI多项式拟合模型来计算单个像元的积雪覆盖率，即通过构建以NDSI为变量的积雪覆盖率估算多项式模型，用于计算如MODIS 1km(或500米)分辨率数据中单个像元的积雪覆盖率，其表达式为：

SF＝a×NDSI+b             (9)

SF＝c×NDSI²+d×NDSI+e    (10)

上式(9)、(10)中，a、b系数取值分别为1.21和0.06，c、d、e系数取值分别为0.26、0.37和0.18。

针对大面积比较连续(较强的降雪)的积雪区域，在进行积雪覆盖率计算时可采用二次多项式模型；而对于较为稀疏的积雪覆盖(不连续斑点状的积雪覆盖情形)，则采用一次多项式模型。

具体计算方法包括：

A.首先使用多项式模型，以单个像元为运算单位，计算积雪覆盖率；

B.以单个像元为运算单位，将MODIS云识别结果与上步中的积雪覆盖率计算结果进行运算，运算方法：(MODIS云识别结果＝＝0)×积雪覆盖率。由于MODIS云识别结果中单个像元取值只能为0或1(0表示该像元没有云覆盖、1表示该像元被云覆盖)，经过该步骤运算，则仅有未被云覆盖的像元进行积雪覆盖率的运算(有云的区域不能得到正确的计算结果，因此要去除掉)。

进一步地，草群高度估算模块24又包括：

MODIS云识别模块211，用于根据预处理后得到的MODIS反射率与亮度温度数据进行MODIS云识别，得到无云区域MODIS L1B数据、有云区域MODIS L1B数据；

草群高度计算模块217，连接MODIS云识别模块211，用于根据无云区域MODIS L1B数据进行MODIS草群高度计算，得到单日草群高度。其是采用NDVI统计模型法来计算单个像元的草群高度，即通过构建以NDVI为变量的草群高度估算多项式模型、指数模型来计算如MODIS 1km(或500米、250米)分辨率数据中单个像元的草群高度。其利用MODIS L1B数据的第1、2通道来构建归一化差分植被指数NDVI的步骤，其表达式为：

$\mathrm{NDVI}=\frac{{\mathrm{\ρ}}_2-{\mathrm{\ρ}}_1}{{\mathrm{\ρ}}_2+{\mathrm{\ρ}}_1}---\left(11\right)$

上式(11)中，ρ₁、ρ₂分别为MODIS第1、2通道的反射率。

此外，还包括采用归一化差分植被指数NDVI统计模型法来计算单个像元的草群高度(单位：cm)的步骤，NDVI统计模型的表达式为：

线性(Linear)模型：

GH＝a×NDVI+b                      (12)

上式(12)中，a、b系数的一般取值为0.7371、-4.6131。

二次多项式(Quadratic)模型：

GH＝a×NDVI²+b×NDVI+c             (13)

上式(13)中，a、b、c系数的一般取值为0.0062、0.1966、4.9485。

三次多项式(Cubic)模型：

GH＝a×NDVI³+b×NDVI²+c×NDVI+d    (14)

上式(14)中，a、b、c、d系数的一般取值为0.00008、-0.0046、0.6183、0.4100。

指数(Exponential)模型：

GH＝a×e^b×NDVI                    (15)

上式(15)中，a、b系数的一般取值为5.4666、0.0318。

幂(Power)模型：

GH＝a×NDVI^b                       (16)

上式(16)中，a、b系数的一般取值为0.3402、1.1371。

具体计算方法：

A.首先使用NDVI统计模型，以单个像元为运算单位，计算草群高度；

B.以单个像元为运算单位，将MODIS云识别结果与上步中的草群高度计算结果进行运算，运算方法：(MODIS云识别结果＝＝0)×草群高度。由于MODIS云识别结果中单个像元取值只能为0或1(0表示该像元没有云覆盖、1表示该像元被云覆盖)，经过该步骤运算，则仅有未被云覆盖的像元进行草群的运算。

本发明使用的数据是搭载于Terra卫星上的MODIS L1B数据，以及Aqua卫星上的MODIS L1B、AMSR-E L1B数据。使用了MODIS的第1、2、4、6、7波段(空间分辨率为500m)，以及第22、26、27、31、35波段(空间分辨率1km)L1B数据，同时还使用了AMSR-E的18.7、23.8、36.5、89GHz(重采样空间分辨率10km)的L1B数据。原始的MODIS、AMSR-E数据级别为level 1B(L1B)，未经过任何处理。在经过几何纠正和辐射纠正之后，分别生成了表观反射率数据和亮度温度数据。研究中使用到的MODIS和AMSR-E波段属性参数见表1和表2。

表1

表1中的上标“1”是指研究中所有MODIS波段数据重采样至1000米分辨率；上标“2”此处为噪声等效温差。

表2

表2中的上标“1”是指研究中所有AMSR-E通道数据重采样至10km分辨率。

如图2所示，为本发明的草原雪灾遥感监测与灾情评估方法流程图，该方法具体包括如下步骤：

步骤201，数据读取，预处理的步骤，具体是：

从原始L1B级HDF属性数据中，读取相关波段数据，并进行辐射标定计算和空间参考变换。

对于MODIS L1B数据，则计算反射率和光谱辐射亮度值，如下：

R＝scale(SI-offset)           (1)

上式(1)中，R为像元反射率或光谱辐射亮度值，scale，offset为比例系数，可从HDF属性数据中获得，SI为像元的灰度值。此外，基于Plank方程，可由光谱辐射亮度值计算目标亮度温度，如下：

$\mathrm{Tb}=\frac{\mathrm{hc}}{\mathrm{\λk}}\·\frac{1}{\frac{2{\mathrm{hc}}^2}{L{\mathrm{\λ}}^5}+1}---\left(2\right)$

上式(2)中，L为黑体的辐射亮度(Wm^-2μm^-1sr^-1)；h是Plank常数，其值为6.63×10^-34Js；k是Boltzmann常数，其值为1.38×10^-23JK^-1；c为光速(3×10⁸ms^-1)；λ为波长(μm)；Tb为亮度温度(K)。

对于AMSR-E L1B数据，计算每个像元的亮度温度，如下：

Tb＝scale×SI+offset             (3)

上式(3)中，Tb为波段亮度温度(K)，scale，offset为比例系数，可从HDF属性数据中获得，SI为像元的灰度值。

L1B是经过仪器标定，但是没有经过大气校正的数据产品；该数据产品中包含有地理坐标产品，但是“科学数据”和“地理数据”还没有连接；因此该类型的数据需要进行辐射标定计算和地理参考计算(几何校正)。HDF为数据格式；L1B HDF为MODIS和AMSR-E数据共享时的一种数据类型，为卫星下行传输后经解码、预处理以后的初级产品，包含了众多的原始数据信息，因此具体应用时常以该类型数据为基础进行运算。

由于原始的L1B HDF中的科学数据为16bit类型无量纲数据，因此通过步骤201进行运算，计算得到后续步骤中需使用的反射率(R)和亮度温度(Tb)。

步骤202，MODIS云识别决策，用于判断遥感影像中每个像元是否满足云识别条件，其表达式为：

Tb₃₅＜220K

Tb₂₇＜226K

Tb₃₁-Tb₂₂＜-12K             (4)

ρ₂₆＞0.035

ρ₁＞0.18

上式(4)中，ρ₁，ρ₂₆分别为MODIS L1B数据第1、26波段的反射率；Tb₂₂，Tb₂₇，Tb₃₁，Tb₃₅分别为MODIS L1B数据第22、27、31、35波段的亮度温度。

通过步骤202获取无云区域MODIS L1B数据、有云区域MODIS L1B数据。

像元亦称像素或像元点，即影像单元(picture element)，是组成数字化影像的最小单元。遥感影像中的某个像素通常使用波段(通道)、行、列来表示，如pixel(k，i，j)表示第k波段、第i行、第j列的像素值，通常也叫做像元值。

步骤203，根据无云区域MODIS L1B数据进行MODIS积雪像元识别，获取MODIS积雪识别结果，其先是利用MODIS L1B数据的第4、6波段来构建归一化差分积雪指数NDSI的步骤，其表达式为：

$\mathrm{NDSI}=\frac{{\mathrm{\ρ}}_4-{\mathrm{\ρ}}_6}{{\mathrm{\ρ}}_4+{\mathrm{\ρ}}_6}---\left(5\right)$

上式(5)中，ρ₄、ρ₆分别为MODIS L1B数据第4、6波段的反射率。由于Aqua MODIS的第6通道的数据有问题，因此在计算Aqua MODIS NDSI时，采用第7通道反射率ρ₇来代替ρ₆进行计算。此外，还包括用于判断每个像元的归一化差分积雪指数NDSI，绿、近红外波段反射率，是否满足积雪像元的识别条件的步骤：

$\left\{\begin{array}{c}\mathrm{NDSI}\≥0.4\\ {\mathrm{\ρ}}_2>0.11\\ {\mathrm{\ρ}}_4\≥0.1\end{array}\right.---\left(6\right)$

上式(6)中，ρ₂、ρ₄分别为MODIS L1B数据第2、4波段的反射率。

步骤204，AMSR-E积雪识别决策，用于判断AMSR-E L1B数据中每个像元是否满足积雪识别条件。这里，针对AMSR-E数据，定义散射指数(scatterindex)具有如下的形式；

scat＝max(Tb_18.7V-Tb_36.5V-3，Tb_23.8V-Tb_89V-3，Tb_36.5V-Tb_89V-1)      (7)

上式(7)中，scat为散射指数，max为取最大值函数，Tb为亮度温度，Tb右下标中的数字为通道频率，字母为极化方式(V表示垂直极化、H表示水平极化)，如Tb_18.7V表示18.7GHz通道垂直极化亮度温度。该步骤中，首先根据散射指数确定散射体，散射体是根据散射指数判断出来的像元，然后再根据各条件式依次识别出降雨、寒漠、冻土等地表，最后识别出积雪区。其具体实现见图3所示。

该步骤针对AMSR-E L1B数据，在无云和有云区域都能识别出积雪像元。而步骤203中的判断步骤针对MODIS L1B数据，且仅能识别无云区域的积雪像元。

由于AMSR-E数据的空间分辨率较低，因此其积雪识别结果仅作为MODIS数据积雪识别结果的一个补充，即实现有云区域的积雪识别。

步骤205，MODIS-AMSR-E数据融合识别积雪，用于对无云区域的MODIS积雪识别结果与有云区域的AMSR-E积雪识别结果进行多源数据融合，实现完整的全天候的积雪监测，得到单日积雪区域，并可用于积雪持续日数的统计。其具体实现见图4所示。

步骤206，AMSR-E积雪深度(Snow Depth，SD)估算，其是通过构建积雪深度反演模型来估算的，采用AMSR-E L1B数据的18.7GHz和36.5GHz通道水平极化亮度温度差(散射指数)进行积雪深度估算，并在此之上构建了积雪深度反演的线性模型，其表达式为：

SD＝a×(Tb_18.7H-Tb_36.5H)+b            (8)

上式(8)中，a、b两个系数的取值分别为1.59cm/K和0，Tb_18.7H和Tb_36.5H分别表示18.7GHz水平极化亮度温度、36.5GHz水平极化亮度温度。

此外，还可以根据草原区域的地域特征，采用不同的模型系数来估算积雪深度。

具体计算方法：

A.首先使用线性模型，以单个像元为运算单位，计算积雪深度；

B.以单个像元为运算单位，将AMSR-E积雪识别结果与上步中的积雪深度计算结果进行运算，运算方法：(AMSR-E积雪识别结果＝＝1)×积雪深度。由于AMSR-E积雪识别结果中单个像元取值只能为0或1(0表示该像元未积雪、1表示该像元积雪)，经过该步骤运算，则仅有积雪像元进行积雪深度的运算，因此能降低线性模型对某些特定区域的计算误差。

步骤207，MODIS积雪覆盖率(Snow Fractral，SF)估算，其是采用NDSI多项式拟合模型来计算单个像元的积雪覆盖率，即通过构建以NDSI为变量的积雪覆盖率估算多项式模型，用于计算如MODIS 1km(或500米)分辨率数据中单个像元的积雪覆盖率，其表达式为：

SF＝a×NDSI+b             (9)

SF＝c×NDSI²+d×NDSI+e    (10)

上式(9)、(10)中，a、b系数取值分别为1.21和0.06，c、d、e系数取值分别为0.26、0.37和0.18。

步骤208，MODIS草群高度(Grass Height，GH)估算，其是采用NDVI统计模型法来计算单个像元的草群高度，即通过构建以NDVI为变量的草群高度估算多项式模型、指数模型来计算如MODIS 1km(或500米、250米)分辨率数据中单个像元的草群高度。其利用MODIS L1B数据的第1、2通道来构建归一化差分植被指数NDVI的步骤，其表达式为：

$\mathrm{NDVI}=\frac{{\mathrm{\ρ}}_2-{\mathrm{\ρ}}_1}{{\mathrm{\ρ}}_2+{\mathrm{\ρ}}_1}---\left(11\right)$

上式(11)中，ρ₁、ρ₂分别为MODIS第1、2通道的反射率。

此外，还包括采用归一化差分植被指数NDVI统计模型法来计算单个像元的草群高度(单位：cm)的步骤，NDVI统计模型的表达式为如下任意一种模型：

线性(Linear)模型：

GH＝a×NDVI+b                      (12)

上式(12)中，a、b一般取值为0.7371、-4.6131。

二次多项式(Quadratic)模型：

GH＝a×NDVI²+b×NDVI+c             (13)

上式(13)中，a、b、c一般取值为0.0062、0.1966、4.9485。

三次多项式(Cubic)模型：

GH＝a×NDVI³+b×NDVI²+c×NDVI+d    (14)

上式(14)中，a、b、c、d一般取值为0.00008、-0.0046、0.6183、0.4100。

指数(Exponential)模型：

GH＝a×e^b×NDVI                    (15)

上式(15)中，a、b一般取值为5.4666、0.0318。

幂(Power)模型：

GH＝a×NDVI^b             (16)

上式(16)中，a、b一般取值为0.3402、1.1371。

步骤209，积雪合成的步骤，基于步骤205、206、207、208中单日的积雪和草群高度监测结果，对多日的积雪覆盖和草群高度进行统计合成，包括积雪持续日数、平均积雪深度、平均积雪覆盖率和平均草群高度。其中，积雪持续日数为连续积雪的时间，其合成算法包括：使用当日的遥感数据识别积雪像元；提取前一日积雪持续日数图中的积雪日数，如果某像元当日积雪，则新积雪持续日数图中该像元处的积雪日数等于前一日积雪日数加一，否则该像元处的积雪日数为零(重新开始统计)。

上述平均积雪深度合成算法包括：根据合成周期内的单日积雪深度图，以单个像元为运算单位，取每个像元的单日积雪深度值进行平均计算，得到多日平均积雪深度；

上述平均积雪覆盖率合成算法包括：根据合成周期内的单日积雪覆盖率图，以单个像元为运算单位，取每个像元的单日积雪覆盖率值进行平均计算，得到多日平均积雪覆盖率；

上述平均草群高度合成算法包括：根据合成周期内的单日草群高度图，以单个像元为运算单位，取每个像元的单日高群高度值进行平均计算，得到多日平均高群高度；

步骤210，雪灾等级计算，其是基于步骤209中的平均积雪深度、平均草群高度计算牧区雪灾等级评价中的因子积雪厚度/草群高度，基于步骤209中的积雪持续日数计算牧区雪灾等级评价中的因子积雪持续时间(日)，基于步骤209中的平均积雪覆盖率计算牧区雪灾等级评价中的因子积雪面积/土地面积，并以牧区雪灾等级国家标准(表3)为基础，对草原区域的雪灾等级进行评价。

表3：

上表中的因子积雪持续时间(日)即步骤209中的积雪持续日数；

上表中的因子积雪面积/土地面积即步骤209中的平均积雪覆盖率；

上表中的因子积雪厚度/草群高度即步骤209中的平均积雪深度/平均草群高度。

上述步骤中的系数值引自公开发表的论文，一般情况下，上述步骤中的各公式类型不可为其他值，但如有特殊需要，可针对特定的时间和区域对系数稍作修改，所以本发明中还提供自定义修改系数的功能。

在图2中，步骤203与步骤204之间无严格的先后顺序关系，步骤205、步骤206、步骤207、步骤208无严格的先后顺序关系。

如图3所示，是本发明的AMSR-E积雪识别方法流程图，该识别方法是通过判断AMSR-E数据中每个像元是否满足积雪识别条件来实现的，其具体实现步骤如下：

步骤301，根据散射指数确定散射体；

步骤302，根据如下积雪识别条件判断是否为降雨，若满足下列条件，则判断为降雨，否则，继续步骤303；

Tb_23.8V＞260K or Tb_23.8V≥254K and scat≤3K or Tb_23.8V≥168+0.49xTb_89V

步骤303，根据如下积雪识别条件判断是否为寒漠，若满足下列条件，则判断为寒漠，否则，继续步骤304；

Tb_18.7V-Tb_18.7H≥18K

Tb_18.7V-Tb_36.5V≤12K and Tb_36.5V-Tb_89V≤13K

步骤304，根据如下积雪识别条件判断是否为冻土，若满足下列条件，则判断为冻土，否则，继续步骤305；

Tb_18.7V-Tb_36.5V≤5K and Tb_23.8V-Tb_89V≤8K

Tb_18.7V-Tb_18.7H≥8K

步骤305，为积雪。

如图4所示，是本发明的获取单日积雪区域的方法流程图，其是对无云区MODIS积雪识别结果与有云区AMSR-E积雪识别结果进行融合而得到的，具体实现步骤如下：

步骤401，读取MODIS L1B数据、AMSR-E L1B数据，并对其进行数据预处理，得到预处理结果：MODIS反射率与亮度温度数据、AMSR-E亮度温度数据；

步骤402，根据MODIS反射率与亮度温度数据进行MODIS云识别，得到无云区域MODIS L1B数据、有云区域MODIS L1B数据，转入步骤404；

步骤403，根据AMSR-E亮度温度数据进行AMSR-E积雪像元识别，得到AMSR-E积雪识别结果：AMSR-E识别积雪区域、AMSR-E识别未积雪区域，转入步骤405；

步骤404，根据无云区域MODIS L1B数据进行MODIS积雪像元识别，得到MODIS积雪识别结果：MODIS识别积雪区域、MODIS识别未积雪区域，转入步骤405；

步骤405，对MODIS识别积雪区域、AMSR-E识别积雪区域进行融合，得到单日积雪区域。

如图5所示，是本发明的获取单日积雪深度的方法流程图。该流程具体包括如下步骤：

步骤501，读取AMSR-E L1B数据，并对其进行数据预处理，得到AMSR-E亮度温度数据；

步骤502，根据AMSR-E亮度温度数据进行AMSR-E积雪像元识别，得到AMSR-E积雪识别结果：AMSR-E识别积雪区域、AMSR-E识别未积雪区域；

步骤503，根据AMSR-E识别积雪区域进行AMSR-E积雪深度计算，得到单日积雪深度。

如图6所示，是本发明的获取单日积雪覆盖率的方法流程图。该流程具体包括如下步骤：

步骤601，读取MODIS L1B数据，并对其进行数据预处理，得到MODIS反射率与亮度温度数据；

步骤602，根据MODIS反射率与亮度温度数据进行MODIS云识别，得到无云区域MODIS L1B数据、有云区域MODIS L1B数据；

步骤603，根据无云区域MODIS L1B数据进行MODIS积雪覆盖率计算，得到单日积雪覆盖率。

如图7所示，是本发明的获取单日草群高度的方法流程图。该流程具体包括如下步骤：

步骤701，读取MODIS L1B数据，并对其进行数据预处理，得到MODIS反射率与亮度温度数据；

步骤702，根据MODIS反射率与亮度温度数据进行MODIS云识别，得到无云区域MODIS L1B数据、有云区域MODIS L1B数据；

步骤703，根据无云区域MODIS L1B数据进行MODIS草群高度计算，得到单日草群高度。

本发明所提供的草原雪灾遥感监测与灾情评估系统及方法，主要是通过光学传感器MODIS，用于云识别、积雪像元识别(无云区)、积雪覆盖率计算、草群高度计算；通过被动微波传感器AMSR-E，用于积雪像元识别(有云区)、积雪深度计算。通过二者合成的积雪像元识别结果，进行积雪持续日数的计算。可对草原积雪情况进行全天候动态监测，并可对草原雪灾发生的时间和地点进行预警。

当然，本发明还可有其他多种实施例，在不背离本发明精神及其实质的情况下，熟悉本领域的技术人员当可根据本发明作出各种相应的改变和变形，但这些相应的改变和变形都应属于本发明所附的权利要求的保护范围。

Claims (11)

1.一种草原雪灾遥感监测与灾情评估系统，其特征在于，包括：

预处理模块，用于读取并处理MODIS L1B数据、AMSR-E L1B数据，得到MODIS反射率、MODIS亮度温度数据、AMSR-E亮度温度数据；

积雪区域融合模块，连接所述预处理模块，用于根据所述MODIS反射率、所述MODIS亮度温度数据进行MODIS积雪像元识别，根据所述AMSR-E亮度温度数据进行AMSR-E积雪像元识别，并对得到的识别结果进行数据融合处理，得到单日积雪区域；

积雪深度估算模块，连接所述预处理模块，用于根据所述AMSR-E亮度温度数据进行AMSR-E积雪深度估算，得到单日积雪深度；

积雪覆盖率估算模块，连接所述预处理模块，用于根据所述MODIS反射率、所述MODIS亮度温度数据进行MODIS积雪覆盖率估算，得到单日积雪覆盖率；

草群高度估算模块，连接所述预处理模块，用于根据所述MODIS反射率、所述MODIS亮度温度数据进行MODIS草群高度估算，得到单日草群高度；

积雪合成模块，连接所述积雪区域融合模块、所述积雪深度估算模块、所述积雪覆盖率估算模块、所述草群高度估算模块，用于根据所述单日积雪区域、所述单日积雪覆盖率、所述单日积雪深度、所述单日草群高度，对多日的积雪区域、积雪覆盖率、积雪深度、草群高度进行统计合成，得到积雪持续日数、平均积雪覆盖率、平均积雪深度、平均草群高度；

雪灾等级评价模块，连接所述积雪合成模块，用于根据所述积雪持续日数、所述平均积雪覆盖率、所述平均积雪深度、所述平均草群高度，对草原区域的雪灾等级进行评价。

2.根据权利要求1所述的草原雪灾遥感监测与灾情评估系统，其特征在于，所述预处理模块以如下公式得到所述MODIS反射率：

R=scale(SI-offset)

式中，R为MODIS反射率，scale，offset为比例系数，该比例系数从HDF属性数据中获得，HDF为数据类型，SI为像元的灰度值；

所述预处理模块以如下公式得到所述MODIS亮度温度数据：

$\mathrm{Tb}1=\frac{\mathrm{hc}}{\mathrm{\λk}}\·\frac{1}{\frac{{2\mathrm{hc}}^2}{{\mathrm{L\λ}}^5}}$

式中，L为黑体的辐射亮度（Wm^-2μm^-1sr^-1）；h是Plank常数，其值为6.63×10^-34Js；k是Boltzmann常数，其值为1.38×10^-23JK^-1；c为光速（3×10⁸ms^-1）；λ为波长（μm）；Tb1为MODIS亮度温度（K）；

所述预处理模块以如下公式得到所述AMSR-E亮度温度数据：

Tb2=scale×SI+offset

式中，Tb2为AMSR-E亮度温度（K），scale，offset为比例系数，该比例系数从HDF属性数据中获得，SI为像元的灰度值。

3.根据权利要求2所述的草原雪灾遥感监测与灾情评估系统，其特征在于，所述积雪区域融合模块进一步包括：

MODIS云识别模块，用于根据所述MODIS反射率、所述MODIS亮度温度数据进行云识别，区分无云、有云区域MODIS LIB数据；

MODIS积雪像元识别模块，连接所述MODIS云识别模块，用于根据所述无云区域MODIS L1B数据进行MODIS积雪像元识别，得到MODIS识别积雪区域；

AMSR-E积雪像元识别模块，用于根据所述AMSR-E亮度温度数据进行AMSR-E积雪像元识别，得到AMSR-E识别积雪区域；

数据融合模块，连接所述MODIS积雪像元识别模块、所述AMSR-E积雪像元识别模块，用于对所述MODIS识别积雪区域、所述AMSR-E识别积雪区域进行多源数据融合，得到所述单日积雪区域。

4.根据权利要求3所述的草原雪灾遥感监测与灾情评估系统，其特征在于，所述MODIS云识别模块通过判断所述MODIS L1B数据中每个像元是否满足云识别条件获取所述无云区域MODIS L1B数据，云识别条件的表达式为：

Tb₃₅＜220K

Tb₂₇＜226K

Tb₃₁-Tb₂₂＜-12K

ρ₂₆＞0.035

ρ₁＞0.18

式中，ρ₁,ρ₂₆分别为MODIS L1B数据第1、26波段的反射率；Tb₂₂,Tb₂₇,Tb₃₁,Tb₃₅分别为MODIS L1B数据第22、27、31、35波段的亮度温度；

所述MODIS积雪像元识别模块通过判断所述无云区域MODIS L1B数据中每个像元是否满足积雪像元识别条件得到所述MODIS识别积雪区域，该积雪像元识别条件的表达式为：

$\left\{\begin{array}{c}\mathrm{NDSI}\≥0.4\\ {\mathrm{\ρ}}_2>0.11\\ {\mathrm{\ρ}}_4\≥0.1\end{array}\right.$

式中，ρ₂、ρ₄分别为MODIS L1B数据第2、4波段的反射率，NDSI为每个像元的归一化差分积雪指数；

以MODIS L1B数据的第4、6波段来构建NDSI，表达式为：

$\mathrm{NDSI}=\frac{{\mathrm{\ρ}}_4-{\mathrm{\ρ}}_6}{{\mathrm{\ρ}}_4+{\mathrm{\ρ}}_6}$

式中，ρ₄、ρ₆分别为MODIS L1B数据第4、6波段的反射率。

5.根据权利要求1或2所述的草原雪灾遥感监测与灾情评估系统，其特征在于，所述积雪深度估算模块进一步包括：

AMSR-E积雪像元识别模块，用于根据所述AMSR-E亮度温度数据进行AMSR-E积雪像元识别，得到AMSR-E识别积雪区域；

积雪深度计算模块，连接所述AMSR-E积雪像元识别模块，用于根据所述AMSR-E识别积雪区域进行AMSR-E积雪深度计算，得到所述单日积雪深度。

6.根据权利要求5所述的草原雪灾遥感监测与灾情评估系统，其特征在于，所述积雪深度计算模块通过构建积雪深度反演模型得到所述单日积雪深度，该积雪深度反演模型的表达式为：

SD=a×(Tb_18.7H-Tb_36.5H)+b

式中，SD为单日积雪深度，a、b系数的取值分别为1.59cm/K和0，Tb_18.7H和Tb_36.5H分别表示AMSR-E L1B数据的18.7GHz通道的水平极化亮度温度、36.5GHz通道的水平极化亮度温度。

7.根据权利要求4所述的草原雪灾遥感监测与灾情评估系统，其特征在于，所述积雪覆盖率估算模块进一步包括：

MODIS云识别模块，用于根据所述MODIS反射率、所述MODIS亮度温度数据进行云识别，区分无云、有云区域MODIS LIB数据；

积雪覆盖率计算模块，连接所述MODIS云识别模块，用于根据无云区域MODIS L1B数据进行MODIS积雪覆盖率计算，得到所述单日积雪覆盖率。

8.根据权利要求7所述的草原雪灾遥感监测与灾情评估系统，其特征在于，所述积雪覆盖率计算模块通过构建以NDSI为变量的多项式拟合模型得到所述单日积雪覆盖率，该多项式拟合模型的表达式为：

SF＝a×NDSI+b或

SF＝c×NDSI²+d×NDSI+e

式中，SF为单日积雪覆盖率，a、b系数的取值分别为1.21和0.06，c、d、e系数取值分别为0.26、0.37和0.18。

9.根据权利要求2、3、4或6所述的草原雪灾遥感监测与灾情评估系统，其特征在于，所述草群高度估算模块进一步包括：

MODIS云识别模块，用于根据所述MODIS反射率、所述MODIS亮度温度数据进行云识别，获取无云区域MODIS L1B数据；

草群高度计算模块，连接所述MODIS云识别模块，用于根据所述无云区域MODIS L1B数据进行MODIS草群高度计算，得到所述单日草群高度。

10.根据权利要求9所述的草原雪灾遥感监测与灾情评估系统，其特征在于，所述草群高度计算模块通过NDVI统计模型法得到所述单日草群高度，NDVI统计模型的表达式为如下任一种：

GH＝a×NDVI+b

式中，a、b系数的一般取值为0.7371、-4.6131，GH为单日草群高度；

GH＝a×NDVI²+b×NDVI+c

式中，a、b、c系数的一般取值为0.0062、0.1966、4.9485；

GH＝a×NDVI³+b×NDVI²+c×NDVI+d

式中，a、b、c、d系数的一般取值为0.00008、-0.0046、0.6183、0.4100；

GH＝a×e^b×NDVI

式中，a、b系数的一般取值为5.4666、0.0318；

GH＝a×NDVI^b

式中，a、b系数的一般取值为0.3402、1.1371；

以MODIS L1B数据的第1、2通道来构建NDVI，表达式为：

$\mathrm{NDVI}=\frac{{\mathrm{\ρ}}_2-{\mathrm{\ρ}}_1}{{\mathrm{\ρ}}_2+{\mathrm{\ρ}}_1}$

式中，ρ₁、ρ₂分别为MODIS L1B数据的第1、2通道的反射率。

11.一种草原雪灾遥感监测与灾情评估方法，应用于上述权利要求1～10中任一项所述的系统，其特征在于，包括：

数据预处理步骤，读取并处理MODIS L1B数据、AMSR-E L1B数据，得到MODIS反射率、MODIS亮度温度数据、AMSR-E亮度温度数据；

积雪区域融合步骤，根据所述MODIS反射率、所述MODIS亮度温度数据进行MODIS积雪像元识别，根据所述AMSR-E亮度温度数据进行AMSR-E积雪像元识别，并对得到的识别结果进行数据融合处理，得到单日积雪区域；

积雪深度估算步骤，根据所述AMSR-E亮度温度数据进行AMSR-E积雪深度估算，得到单日积雪深度；

积雪覆盖率估算步骤，根据所述MODIS反射率、所述MODIS亮度温度数据进行MODIS积雪覆盖率估算，得到单日积雪覆盖率；

草群高度估算步骤，根据所述MODIS反射率、所述MODIS亮度温度数据进行MODIS草群高度估算，得到单日草群高度；

积雪合成步骤，根据所述单日积雪区域、所述单日积雪覆盖率、所述单日积雪深度、所述单日草群高度，对多日的积雪区域、积雪覆盖率、积雪深度、草群高度进行统计合成，得到积雪持续日数、平均积雪覆盖率、平均积雪深度、平均草群高度；

雪灾等级评价步骤，根据所述积雪持续日数、所述平均积雪覆盖率、所述平均积雪深度、所述平均草群高度，对草原区域的雪灾等级进行评价。

Images