CN101799561A - 一种基于减灾小卫星的雪灾遥感监测模拟评估方法 - Google Patents

Abstract

本发明公开了一种基于减灾小卫星的雪灾遥感监测模拟评估方法，包括以下步骤：根据积雪光谱特征识别积雪，形成积雪遥感图像数据；建立雪灾案例数据库和雪灾统计数据库，所述雪灾案例数据库包括：成灾时间、受灾地点、积雪日数、归一化植被指数NDVI、受灾人口数和积雪面积，所述雪灾统计数据库包括：总降水量、最高气温、最低气温、日均温度、风速和地表温度；对识别出的积雪做积雪深度反演，并将积雪深度分别存入雪灾案例数据库和雪灾统计数据库；根据雪灾案例数据库和雪灾统计数据库中数据对雪灾进行模拟评价。本发明为雪灾监测提供了新的、适用性强的遥感数据源，对中国雪灾的减灾救灾工作具有重要意义。

Description

一种基于减灾小卫星的雪灾遥感监测模拟评估方法

技术领域

本发明涉及卫星遥感技术领域，特别涉及一种基于减灾小卫星的雪灾遥感监测模拟评估方法。

背景技术

中国是世界上少数几个航天大国之一，但在雪灾遥感方面却存在不少问题。一方面，大量的遥感数据需要靠国外进口，不仅费用高，实时获取困难，而且国外数据在很多方面不能满足中国雪灾管理工作的实际需要。另一方面，中国自主的卫星遥感数据在雪灾管理方面的应用还不够广泛和深入，迫切需要加大研究力度，对于中国目前已发射的卫星及2007年发射的环境与灾害小卫星星座都需要迫切解决这一问题。

小卫星星座上天后，对雪灾具备怎样的监测能力，有待于进行预先研究，同时，提前结合卫星有效载荷的技术指标对雪灾监测和灾情评估的有关关键技术进行针对性的研究，是尽快发挥小卫星星座效益的有力保障。因此需要提前利用遥感数据(如ASTER、MODIS等)模拟小卫星星座光学星数据，获取不同波段的积雪光谱特征信息，根据积雪在不同波段光谱差异，建立积雪遥感提取模型，并建立在GIS支持下的雪灾灾情综合评估模型。总之，要充分模拟小卫星星座的各种载荷，使之满足雪灾管理的需求，是一项难度较大且亟待解决的任务。

发明内容

(一)解决的技术问题

本发明的目的是提供一种基于减灾小卫星的雪灾遥感监测模拟评估方法，以解决中国自主的卫星遥感数据在雪灾监测方面的不足的问题。

(二)技术方案

一种基于减灾小卫星的雪灾遥感监测模拟评估方法，包括以下步骤：

S1：根据积雪光谱特征识别积雪，得到积雪遥感图像；

S2：建立雪灾案例数据库和雪灾统计数据库，所述雪灾案例数据库包括：成灾时间、受灾地点、积雪日数、归一化植被指数NDVI、受灾人口数和积雪面积，所述雪灾统计数据库包括：总降水量、最高气温、最低气温、日均温度、风速和地表温度；

S3：对识别出的积雪遥感图像做积雪深度反演，并将积雪深度分别存入雪灾案例数据库和雪灾统计数据库；

S4：根据雪灾案例数据库和雪灾统计数据库中数据对雪灾进行模拟评价。

其中，步骤S1包括利用中分辨率成像光谱仪MODIS数据模拟小卫星数据识别积雪的步骤，具体包括：

S1011：用MODIS第4和第2波段光谱代替小卫星CCD相机第2和第4波段光谱，用MODIS第6波段光谱代替小卫星红外相机第2波段光谱；

S1012：根据MODIS第4波段反射率和第6波段反射率计算归一化雪被指数NDSI，公式如下：

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

其中，ρ₄和ρ₆为MODIS第4波段的反射率和第6波段的反射率；

S1013：计算MODIS第29波段与第31波段的亮温差T29-T31；

S1014：若T29-T31≥ΔT，并NDSI≥0.4、ρ₄＞0.1和ρ₂＞0.11时执行步骤S1015，其中ρ₂为MODIS第2波段反射率；

S1015：利用陆地/水掩膜数据排除水体范围，去除湖冰，得到陆地积雪和冰川遥感图像数据。

其中，所述步骤S1还包括：

S1021：利用NOAA/AVHRR的第1波段的反射率和第2波段的反射率计算归一化植被指数NDVI，公式如下：

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

其中ρ₁和ρ₂为NOAA/AVHRR第1波段的反射率和第2波段的反射率；

S1022：计算NOAA/AVHRR第3波段与第4波段的亮温差T3-T4；

S1023：判断条件ρ₂≥a、ρ₅≤b、NDVI≤c和T3-T4≤d，若满足，则得到的遥感图像为积雪图像，其中a和b为反射率阈值，c为NDVI阈值，d为亮温差阈值。

其中，所述步骤S3包括：

S301：选取MODIS的1～7波段的反射率、积雪指数NDSI、积雪表面温度Ts、坡度slope、坡向aspect和高程elevation作为积雪深度反演指标。

S302：计算所述积雪深度反演指标，建立以坡度分段的反演模型，反演积雪深度；

其中，步骤S4包括：

S401：将雪灾案例数据库中数据和雪灾统计数据库中数据分成致灾因子、孕灾环境和抗灾条件三个评价指标，并建立雪灾理论评价指标体系；

S402：去除雪灾理论评价指标体系冗余指标，建立雪灾实际评价指标体系，所述实际评价体系包括孕灾环境和致灾因子，所述孕灾环境包括指标：积雪日数、积雪面积和NDVI，所述致灾因子包括指标：降水总量、最高气温、最低气温、风速和积雪深度；

S403：根据所述雪灾实际评价指标体系的评价指标建立雪灾灾情综合评估模型。

其中，所述步骤S403包括：

S4031：利用z-score算法归一化所述评价指标，公式如下：

$y_i=\frac{x_i-\stackrel{\‾}{x}}{s}$

其中，x_i为归一化前的评价指标，

$\stackrel{\‾}{x}=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i$

为归一化前的评价指标平均值，

$s=\sqrt{\frac{1}{n-1}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(x_i-\stackrel{\‾}{x})}^2}$

为归一化前的评价指标标准差，y_i为归一化后的评价指标；

S4032：根据归一化后的评价指标利用多元回归方法确定评价指标的权重；

S4033：根据评价指标权重，按如下公式计算雪灾评判系数以评价雪灾，

$\mathrm{SNOWD}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\μ}}_i\·{\mathrm{\ω}}_i,(i=\mathrm{1,2,3},...,n)$

其中，SNOWD为被评价区域的雪灾评判系数，μ_i为各个指标的归一化值，ω_i为各个指标的权重，n为评价指标的个数。

(三)有益效果

1、减灾小卫星数据具有较高分辨率，将提高雪灾灾情全天候、全天时、宽覆盖的动态监测能力，为有效地进行灾情评估提供支持与保障。

2、减灾小卫星能够使用建立的雪灾灾情综合评估模型进行雪灾灾情的全面评估，为救灾提供决策依据。

3、减灾小卫星CCD相机和红外相机在波段设置上基本上覆盖了现有雪灾遥感监测模型使用的波段，因此能够在现有模型的基础上进行雪灾遥感监测。

4、减灾小卫星能够满足现有雪灾遥感监测模型参数的获取，实现牧区雪灾遥感监测及其业务运行。

5、HJ-1小卫星星座的建设填补了我国雪灾遥感监测领域信息源供应自主化的空白，是我国牧区雪灾监测与防灾减灾救灾措施实施的重要基础。减灾小卫星星座数据能够满足不同层次、不同方面的区域业务化雪灾监测工作的需要。

附图说明

图1是根据本发明的一种基于减灾小卫星的雪灾遥感监测模拟评估方法流程图；

图2是积雪光谱曲线图；

图3是根据本发明方法中利用MODIS数据模拟小卫星数据进行冰雪分布信息提取算法流程；

图4是根据本发明方法中利用NOAA/AVHRR数据模拟小卫星数据进行冰雪分布信息提取算法流程；

图5是根据本发明方法的牧区雪灾理论评价指标体系结构图；

图6是根据本发明方法的牧区雪灾实际评价指标体系结构图；

图7是根据本发明方法的雪灾评估模型评估精度验证图。

具体实施方式

本发明提出的一种基于减灾小卫星的雪灾遥感监测模拟评估方法，结合附图和实施例说明如下。

如图1所示，步骤S1中根据积雪光谱特征识别积雪，典型积雪的光谱曲线如图2所示。经分析，积雪在可见光波段(0.38μm-0.72μm)和近红外波段(0.72μm-1.30μm)，均表现出高的反射率，但在近红外波段(1.30μm-3.00μm)，尤其是1.55μm和2.10μm波段附近由于雪吸收太阳辐射反射率很低。同时，鉴于积雪提取主要在于区分云、地表物等信息，经研究发现积雪、云和地表物的光谱有如下特征，如表1所示。在MODIS的多光谱数据中，第二通道(0.858μm)和第四通道(0.555μm)为可见光通道，可用于区别云雪和其它地物；第六通道(1.640μm)近红外通道，可以用于云和雪的区分。在NOAA/AVHRR数据的多光谱数据中，第一通道(0.58μm-0.68μm)和第二通道(0.72μm-1.1μm)为可见光和近红外通道，可以区分云雪和其它地物。根据云的波谱特征，云NDSI值很低，因此NDSI可以作为云、雪的判别指数。但是通过试验发现，水体的NDSI值也非常高，因此极易将水体也误判为积雪像元。水体在近红外波段0.76μm附近的反射率下降幅度很大，而雪仍然具有30％以上的反射率，因此可以利用该波段判别雪与水。对于MODIS而言，即以第二通道作为水体判别的一个因子。另外，积雪在0.55μm附近(MODIS第四通道)的反射率大于10％，可以剔除具有较高NDSI但不是雪的像元。薄云的光谱特性与地表积雪光谱特性非常相似，仅用以上条件会将薄云像元错误判别为积雪。经研究，在NOAA第三通道与第四通道之差影像上，积雪的值明显低于薄云在内的云的信息，因此可用3、4通道之差区分积雪与云(见图5-4-b)；在MODIS第29通道与第31通道的亮温差(T29-T31)影像上，薄云亮温差要普遍高于晴空区域，但是厚云的亮温差比晴空区域还要低，因此，可以利用亮温差区分薄云与其他像元。另外，对于MODIS，第31通道厚云顶部亮温要低于晴空区域亮温，因此将第31通道的亮温T31也作为一个判别因子，以防止将厚云像元误判为雪像元。

表5-1  积雪、云与其它地物波谱特征比较

类型	可见光和近红外波段(0.72μm-1.30μm)	近红外波段(1.30μm-3.00μm)
			积雪	高反射率	低反射率
云	高反射率	高反射率
			一般地物	低反射率(相对于积雪)	高反射率(相对于积雪)

本步骤中利用可通过MODIS数据和NOAA/AVHRR数据分别模拟小卫星数据进行冰雪分布信息的提取。环境与灾害监测预报小卫星星座(以下简称“小卫星”)卫星的重访周期为2天或4天，MODIS、NOAA的卫星重访周期分别为1-2天、1天。MODIS的1、2、4、6、20、31、32通道的波段设置分别与小卫星的CCD相机的3、4、2和红外相机的2、3、4通道的波段设置基本相似(见表2)，相关研究表明，这些波段正是用于积雪遥感监测的最佳波段，加之卫星重访周期基本一致，因此可以用MODIS数据模拟小卫星数据监测雪情。另外，由于中国从2002年5月开始接收EOS卫星的MODIS数据，难以完全满足长期雪灾研究的需要，因此，需要其它数据作为补充。研究发现，NOAA卫星AVHRR传感器的1、2、3、4和5通道的波段设置分别与用于雪灾监测的小卫星CCD相机的2/3、4/红外1通道和红外相机的4通道的波段设置相似(见表3)，因此可用NOAA/AVHRR数据辅助模拟小卫星数据。综合以上因素，本发明拟采用MODIS数据做主要数据源，NOAA/AVHRR数据为辅，模拟小卫星数据进行积雪监测。

表3  小卫星与EOS/MODIS数据传感器波段设置比较表

表3  小卫星与NOAA/AVHRR数据传感器波段设置比较

与小卫星CCD相机第2、4波段对应的分别是MODIS第4、第2波段；与小卫星红外相机第2波段对应的是MODIS第6波段，其算法流程如图3所示，结果为空间分辨率500米的图像数据。用MODIS第4和第2波段光谱代替小卫星CCD相机第2和第4波段光谱，用MODIS第6波段光谱代替小卫星红外相机第2波段光谱；根据MODIS第4波段反射率和第6波段反射率计算归一化雪被指数NDSI，公式如下：

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

其中，ρ₄和ρ₆为MODIS第4波段的反射率和第6波段的反射率；计算MODIS第29波段与第31波段的亮温差T29-T31；若T29-T31≥ΔT，并NDSI≥0.4、ρ₄＞0.1和ρ₂＞0.11时执行步骤S1015，其中ρ₂为MODIS第2波段反射率；利用陆地/水掩膜数据排除水体范围，去除湖冰，得到陆地积雪和冰川遥感图像数据。

利用NOAA/AVHRR数据分别模拟小卫星数据进行冰雪分布信息的提取流程如下：利用NOAA/AVHRR的第1波段的反射率和第2波段的反射率计算归一化植被指数NDVI，公式如下：

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

其中ρ₁和ρ₂为NOAA/AVHRR第1波段的反射率和第2波段的反射率；计算NOAA/AVHRR第3波段与第4波段的亮温差T3-T4；判断条件ρ₂≥a、ρ₅≤b、NDVI≤c和T3-T4≤d，若满足，则得到的遥感图像为积雪图像，其中a和b为反射率阈值，c为NDVI阈值，d为亮温差阈值。

步骤S2中建立雪灾案例数据库和雪灾统计数据库，以民政部官方统计资料为准，并以气象部门、互联网、草原管理部门以及其它相关部门的数据信息为有益补充，结合野外实测调查，建立了从1999年底至2006年2月覆盖青海、内蒙、西藏、新疆四省区的中国主要牧区雪灾统计数据库和雪灾案例数据库，包含雪灾案例近五十次，较为全面地记录并反映了这一时段内雪灾的分布及成灾机理。

雪灾案例数据库主要反映雪灾灾情信息，包括成灾时间、受灾地点、受灾面积、平均积雪深度、积雪日数、归一化植被指数NDVI、受灾人口数、牲畜伤亡头数、直接财产损失以及受灾描述等属性数据。统计数据库指雪灾发生期间的气象、地面状况等相关信息，包括站点经纬度、时间、平均雪深、日最高温度、日最低温度、日均温度、最大雪深、积雪日数、风速、地表温度等属性数据。

步骤S3中对识别出的积雪做积雪深度反演，并将积雪深度分别存入雪灾案例数据库和雪灾统计数据库。具体包括：

对遥感获得的积雪图像数据预处理，包括：去除图像边缘重叠现象，输出MODIS头文件的控制点GCP，根据控制点对图像进行几何校正，辐射定标。具体地，使用ENVI4.0的Bow-tie correction功能去除Bow-tie现象，经过Bow-tie去除后的影像与原影像对比发现，经过Bow-tie处理的影像边缘基本不存在重叠，效果较好；运用ENVI4.0的Export GCPs功能，从MODIS头文件中选择并输出控制点文件；根据已输出的控制点文件，使用ENVI4.0的Registration功能，进行几何校正。将校正后的影像与矢量数据对比，效果较好；经过辐射定标将影像灰度(DN)值转换为大气顶层反射率值(对于反射波段而言)或光谱辐亮度值(对于辐射波段而言)。从MODIS头文件中读取相关参数，根据以下公式进行辐射校正，具体如下：

反射波段，大气顶层的反射率，即表观反射率计算如下：

ρ_i＝reflectance_scales_i×(DN-reflectance_offsets_i)式中，ρ_i为表观反射率；i表示波段序号；“refectance_scales_i”、“reflectance_offsets_i”分别为影像增益和偏差，在头文件中读取。

辐射波段，光谱辐亮度值计算如下：。

I_i＝radiance_scales_i×(DN-radiance_offsets_i)

式中，Ii为光谱辐亮度值；i表示波段序号。

有时还需要将定标得到的光谱辐亮度值转换为亮温值，需要通过Planck函数的反函数进行运算，如下：

$T_i=\frac{c_2}{{\mathrm{\λ}}_i\ln [(c_1/{\mathrm{\λ}}_i^5{I}_i)+1]}$

式中，c₁＝hc₂，c₂＝hc/k，I_i是光谱辐亮度值，h为普朗克常数，h＝6.55×10^-27尔格秒，k是波尔兹曼常数，k＝1.380658×10^-23J/K，c是光速，c＝2.992458×10⁸m/s，λi一般取第i波段的中心波长即可，如果精度要求非常高，也可以根据第i波段的波段响应函数取该波段的有效波长。其中，预处理中还包括：气象数据的矢量栅格化、地理数据的投影转化等处理。

选取MODIS的1～7波段的反射率、归一化雪被指数NDSI、积雪表面温度Ts、坡度slope、坡向aspect和高程elevation作为积雪深度反演指标；计算所述积雪深度反演指标，建立以坡度分段的反演模型，反演积雪深度。

计算所述积雪深度反演指标，建立以坡度分段的反演模型，反演积雪深度。将中国牧区气象站点分布图、全国1∶25万DEM、坡度图和坡向图进行叠加，提取每个气象站点所对应的第1～7波段的反射率、NDSI、Ts、坡度、坡向和高程。然后将雪深与以上指标进行相关分析和回归分析，经过分析发现，按照坡度进行分段，每个坡度范围的积雪深度(SD)与各项指标具有较好的相关性。

1、不按坡度分段

利用逐步回归分析法，剔除第7波段的反射率，建立如下积雪深度估算方程：

SD＝-246.68ρ₁+87.57ρ₂+64.108ρ₃+143.84ρ₄-96.365ρ₅+128.167ρ₆-0.124T_s+19.425NDSI-0.00268elevation+0.21slope+0.004742aspect    (1)

线性回归的复相关系数R＝0.534，决定系数R₂＝0.285，调整后的复相关系数Ra＝0.219。

2、按坡度分段

(1)坡度为0

利用逐步回归分析法，剔除第1波段的反射率数据，以及积雪表面温度，建立如下积雪深度估算方程：

SD＝134.979ρ₂+132.224ρ₃-216.066ρ₄-320.702ρ₅+918.546ρ₆-739.980ρ₇+38.337NDSI-0.00142elevation-0.0064aspect

                                                    (2)

线性回归的复相关系数R＝0.642，决定系数R₂＝0.412，调整后的复相关系数Ra＝0.265。

(2)坡度介于0～10间

利用逐步回归分析法，剔除第3、5波段的反射率数据，建立如下积雪深度估算方程：

SD＝-630.411ρ₁+357.620ρ₂+317.224ρ₄-213.281ρ₆+318.230ρ₇-0.208Ts+43.135NDSI-0.00294elevation+0.544slope-0.003575aspect

                                                    (3)

线性回归的复相关系数R＝0.792，决定系数R2＝0.628，调整后的复相关系数Ra＝0.511。

(3)坡度大于10

利用逐步回归分析法，剔除第4、6、7波段的反射率数据，建立如下积雪深度估算方程：

SD＝1211.224ρ₁-1202.375ρ₂-241.546ρ₃+240.161ρ₅-0.262Ts+100.834NDSI-0.000957elevation+1.336slope+0.0349aspect

                                                    (4)

线性回归的复相关系数R＝0.885，决定系数R2＝0.783，调整后的复相关系数Ra＝0.539。

以上复相关系数R和Ra是表示SD的实际观察值与由p个自变量预测的值的相关情况。决定系数R2的大小决定了相关的密切程度，当R2越接近1时，表示相关的方程式参考价值越高，相反，越接近0时，表示参考价值越低。所以这两个值是反映估算方程的准确程度。

总体而言，按照坡度数据进行分段反演积雪厚度，相关性较高，效果较好。因此，本发明基于式(2)、(3)和(4)，建立雪深反演模型。并将积雪深度存储在雪灾案例数据库和雪灾统计数据库。

步骤S4中根据雪灾案例数据库和雪灾统计数据库中数据对雪灾进行模拟评价。具体包括：

将雪灾案例数据库中数据和雪灾统计数据库中数据分成致灾因子、孕灾环境和抗灾条件三个评价指标，并建立雪灾理论评价指标体系，该体系结构包括2个层次，孕灾环境：包括积雪持续日数、积雪面积、载蓄量、冷季草场超载比例、牧草高度；致灾因子：包括平均雪深、日降水量、一次降雪过程的降水总量、日均温度、最高气温、最低气温、风速、低温持续时间；抗灾条件：包括减灾设施、饲料储量等。体系结构图如图5所示。

去除雪灾理论评价指标体系冗余指标，建立雪灾实际评价指标体系，如图6所示。综合考虑指标易获取性、可操作性、代表性，本着精简、优化的原则，从牧区雪灾灾情理论评价指标中筛选前人重点研究频度最高，且对雪灾灾害形成最为重要的牧区雪灾评价指标，同时借助多元回归后向剔除法，去除冗余或对灾情贡献不大等指标，从理论评价指标中选取实际作为灾情业务化评价指标，所述实际评价体系包括孕灾环境和致灾因子，所述孕灾环境包括指标：积雪日数、积雪面积和NDVI，所述致灾因子包括指标：降水总量、最高气温、最低气温、风速和积雪深度。其中，月累计降水量、日均温度、最高气温、最低气温、风速等指标由气象观测所得，从统计数据库获取。积雪面积和积雪日数由遥感影像监测统计计算所得，根据积雪识别方法获取多日积雪信息，再利用统计法计算积雪持续日数，用GIS分析方法计算积雪面积。

反演积雪深度，根据式(2)、(3)和(4)反演积雪深度；

NDVI，对于同一密度草场，牧草越高和牧草生物量相对越大，因此可用草场生物量代替牧草高度作为雪灾灾情评价指标。经前人研究证明，NDVI根据绿色植物叶绿素对蓝光、红光有一强吸收峰值和对近红外有一强反射值这一特性计算出来的，对植被盖度的敏感性较高，是用来衡量生物量大小的简单、易行的植被指数，计算公式如下：

$\mathrm{NDVI}=\frac{{\mathrm{\ρ}}_{\mathrm{NIR}}-{\mathrm{\ρ}}_R}{{\mathrm{\ρ}}_{\mathrm{NIR}}+{\mathrm{\ρ}}_R}$

式中，ρ_NIR为近红外波段反射率；ρ_R为红光波段反射率。因此，本发明用NDVI替代牧区草高作为雪灾灾情评价指标能较为客观准确地反映了地面植物(牧草)生物量。

根据所述雪灾实际评价指标体系的评价指标建立雪灾灾情综合评估模型。在此之前需要确定雪灾灾情综合评价参数。评价灾害损失程度的参数很多，有直接经济损失、受灾人口数、受灾牲口数、死亡牲口数等，要客观真实反映灾害损失程度，选取评价参数时需综合考虑以上因素。

直接经济损失是受灾人口、受灾牲口、死亡牲口等指标共同作用而产生的，因此直接经济损失是其它灾害损失参数的综合，可以概括作为雪灾灾情程度的评价参数。但介于民政系统统计上报现状，较难全面、准确的统计直接经济损失，往往出现数据不全等问题，因此仅以直接经济损失作为评价参数将对评价结果带来较大误差。为避免直接使用经济损失数据带来的评价误差，本发明选取可综合表示直接经济损失的受灾人口、牲口、死亡牲口共同作为评价依据。具体算法是：将统计上报的直接经济损失数据(Y)作为因变量，其它雪灾损失参数作为自变量(x)，进行多元统计回归分析，结果如下：

表4   雪灾灾情综合评价参数a

      相关系数(a)

表5   雪灾灾情综合评价参数b

R	R2	调整的R2
			.934(a)	.873	.778

由上述结果看出，直接经济损失Y与受灾人口和死亡牲口呈较好正相关关系，与受灾牲口关系较小。因此，可用受灾人口和死亡牲口数拟合直接经济损失量，从而作为雪灾灾情综合评价参数，以此衡量灾害损失程度。根据回归分析结果(表4与表5)，雪灾灾情拟合经济损失量(人民币万元)的计算方程如下：

Y＝a×受灾人口+b×死亡牲口数-c

式中，a为相关系数(212[万元万/人])；b为雪灾灾情综合评价参数(245[万元/万头])；c为误差(416[万元])。

利用z-score算法归一化所述评价指标，公式如下：

$y_i=\frac{x_i-\stackrel{\‾}{x}}{s}$

其中，x_i为归一化前的评价指标，

$\stackrel{\‾}{x}=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i$

为归一化前的评价指标平均值，

$s=\sqrt{\frac{1}{n-1}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(x_i-\stackrel{\‾}{x})}^2}$

为归一化前的评价指标标准差，y_i为归一化后的评价指标。

根据归一化后的评价指标利用多元回归方法确定评价指标的权重。变量回归系数(Beta值)在多元回归中常被用来比较变量间的重要性，可作为衡量各自变量对因变量贡献(权重)的依据。本发明将拟合经济损失量与经过归一化处理的评价指标进行回归，生成指标系数(Beta值)(见表6)，以此作为在利用层次分析法进行每层次评估中，对该层指标进行相对重要性判断的依据。

表6  指标回归系数

根据表6中的标准系数列表，确定各指标对灾害损失的贡献和相对权重，对各个自变量的重要程度分别进行定量打分，一定程度上解决了传统层次分析模型中使用专家打分方法确定权重的主观经验性和系统不稳定性，提高了评估的精度，相对重要性矩阵详见表7-表9。

表7  第一层各因子相对重要性矩阵表

	孕灾环境
		致灾因子	1

表8  第二层孕灾环境各因子相对重要性矩阵

	积雪面积	NDVI
			积雪日数	0.5	2
积雪面积		4

表9  第二层致灾因子各因子相对重要性矩阵

	最高气温	最低气温	风速	积雪深度
					总降水量	0.5	0.333	1	1
最高气温		0.5	2	2
					最低气温			3	2
风速				1

根据各因子相对重要性矩阵构造出判断矩阵，并通过求解各判断矩阵的特征值，确定该层各元素的权重，并得到各层元素对目标层的组合权重，见表10。

表10  雪灾灾情综合评价指数模型各因子权重

第一级	第二级	权重
			孕灾环境		0.50
	积雪面积	0.286
				积雪日数	0.143
	NDVI	0.071
			致灾因子		0.50
	最低气温	0.186
				最高气温	0.119
	积雪深度	0.069
				总降水量	0.063
	风速	0.063

根据评价指标权重，按如下公式计算雪灾评判系数以评价雪灾，在数理上，雪灾灾情综合评价指数表示各个评价指标的加权。从表达的实际意义上来讲，雪灾灾情综合评价指数则表示雪灾事件中，各评价指标对雪灾灾害损失的贡献度总和。所以，各评价指标加权，可作为雪灾灾情评判系数衡量雪灾程。

$\mathrm{SNOWD}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\μ}}_i\·{\mathrm{\ω}}_i,(i=\mathrm{1,2,3},...,n)$

其中，SNOWD为被评价区域的雪灾评判系数，μ_i为各个指标的归一化值，ω_i为各个指标的权重，n为评价指标的个数。

将全国牧区雪灾灾情等级分为四级，即：轻灾、中灾、重灾、特大雪灾。结合雪灾样本中灾害描述的灾害严重程度，参考内蒙古自治区气象中心畜牧气象灾害标准(DB15/T255-1997)的白灾等级划分标准，依据直接经济损失值和雪灾灾情评判系数，制定等级评定标准，见表11：

表11  中国牧区雪灾灾情等级评定表

SNOWD(评判系数)	灾害等级	雪灾灾情拟合经济损失量
			-1.0～-0.25	轻灾	对应雪灾拟合经济损失量在2000万以下
-0.25～0.0	中灾	对应雪灾拟合经济损失量在2000万-5000万
			0.0～0.5	重灾	对应雪灾拟合经济损失量在5000万-1.8亿
0.5～1.0	特大雪灾	对应雪灾拟合经济损失量大于1.8亿

为了评价上述评估模型的准确程度，进行模型验证计算。

以雪灾案例库中未参与模型建立的、直接经济损失量分散在各区段上的10次雪灾样本灾进行计算，通过与实际雪灾描述对照，发现有两个雪灾样本按雪灾灾情评判系数公式和表11进行评估的结果与实际上报结果不一致，分级出现偏差，如图7所示。

由此可见，以上雪灾灾情综合评估模型准确程度可达到80％左右，可以业务化使用，辅助灾害评估及救灾决策。

以上实施方式仅用于说明本发明，而并非对本发明的限制，有关技术领域的普通技术人员，在不脱离本发明的精神和范围的情况下，还可以做出各种变化和变型，因此所有等同的技术方案也属于本发明的范畴，本发明的专利保护范围应由权利要求限定。

Claims (6)

1.一种基于减灾小卫星的雪灾遥感监测模拟评估方法，其特征在于，包括以下步骤：

S1：根据积雪光谱特征识别积雪，形成积雪遥感图像数据；

S2：建立雪灾案例数据库和雪灾统计数据库，所述雪灾案例数据库包括：成灾时间、受灾地点、积雪日数、归一化植被指数NDVI、受灾人口数和积雪面积，所述雪灾统计数据库包括：总降水量、最高气温、最低气温、日均温度、风速和地表温度；

S3：对识别出的积雪做积雪深度反演，并将积雪深度分别存入雪灾案例数据库和雪灾统计数据库；

S4：根据雪灾案例数据库和雪灾统计数据库中数据对雪灾进行模拟评价。

2.如权利要求1所述的基于减灾小卫星的雪灾遥感监测模拟评估方法，其特征在于，步骤S1包括利用中分辨率成像光谱仪MODIS数据模拟小卫星数据识别积雪的步骤，具体包括：

S1011：用MODIS第4和第2波段光谱代替小卫星CCD相机第2和第4波段光谱，用MODIS第6波段光谱代替小卫星红外相机第2波段光谱；

S1012：根据MODIS第4波段反射率和第6波段反射率计算归一化雪被指数NDSI，公式如下：

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

其中，ρ₄和ρ₆为MODIS第4波段的反射率和第6波段的反射率；

S1013：计算MODIS第29波段与第31波段的亮温差T29-T31；

S1014：若T29-T31≥ΔT，并NDSI≥0.4、ρ₄＞0.1和ρ₂＞0.11时执行步骤S1015，其中ρ₂为MODIS第2波段反射率；

S1015：利用陆地/水掩膜数据排除水体范围，去除湖冰，得到陆地积雪和冰川遥感图像数据。

3.如权利要求1所述的基于减灾小卫星的雪灾遥感监测模拟评估方法，其特征在于，所述步骤S1还包括：

S1021：利用NOAA/AVHRR的第1波段的反射率和第2波段的反射率计算归一化植被指数NDVI，公式如下：

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

其中ρ₁和ρ₂为NOAA/AVHRR第1波段的反射率和第2波段的反射率；

S1022：计算NOAA/AVHRR第3波段与第4波段的亮温差T3-T4；

S1023：判断条件ρ₂≥a、ρ₅≤b、NDVI≤c和T3-T4≤d，若满足，则得到的遥感图像为积雪图像，其中a和b为反射率阈值，c为NDVI阈值，d为亮温差阈值。

4.如权利要求1所述的基于减灾小卫星的雪灾遥感监测模拟评估方法，其特征在于，所述步骤S3包括：

S301：对遥感获得的积雪图像数据预处理，包括：去除图像边缘重叠现象，输出MODIS头文件的控制点GCP，根据控制点对图像进行几何校正，辐射定标；

S302：选取MODIS的1～7波段的反射率、归一化雪被指数NDSI、积雪表面温度Ts、坡度slope、坡向aspect和高程elevation作为积雪深度反演指标；

S303：计算所述积雪深度反演指标，建立以坡度分段的反演模型，反演积雪深度。

5.如权利要求1所述的基于减灾小卫星的雪灾遥感监测模拟评估方法，其特征在于，步骤S4包括：

S401：将雪灾案例数据库中数据和雪灾统计数据库中数据分成致灾因子、孕灾环境和抗灾条件三个评价指标，并建立雪灾理论评价指标体系；

S402：去除雪灾理论评价指标体系冗余指标，建立雪灾实际评价指标体系，所述实际评价体系包括孕灾环境和致灾因子，所述孕灾环境包括指标：积雪日数、积雪面积和NDVI，所述致灾因子包括指标：降水总量、最高气温、最低气温、风速和积雪深度；

S403：根据所述雪灾实际评价指标体系的评价指标建立雪灾灾情综合评估模型。

6.如权利要求5所述的基于减灾小卫星的雪灾遥感监测模拟评估方法，其特征在于，所述步骤S403包括：

S4031：利用z-score算法归一化所述评价指标，公式如下：

$y_i=\frac{x_i-\stackrel{\‾}{x}}{s}$

其中，x_i为归一化前的评价指标，

$\stackrel{\‾}{x}=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i$

为归一化前的评价指标平均值，

$s=\sqrt{\frac{1}{n-1}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(x_i-\stackrel{\‾}{x})}^2}$

为归一化前的评价指标标准差，y_i为归一化后的评价指标；

S4032：根据归一化后的评价指标利用多元回归方法确定评价指标的权重；

S4033：根据评价指标权重，按如下公式计算雪灾评判系数以评价雪灾，

$\mathrm{SNOWD}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\μ}}_i\·{\mathrm{\ω}}_i,(i=\mathrm{1,2},...,n)$

其中，SNOWD为被评价区域的雪灾评判系数，μ_i为各个指标的归一化值，ω_i为各个指标的权重，n为评价指标的个数。

Images