基于SNOMAP算法的积雪面积信息提取及其在ENVI/IDL中的实现

1.       SNOMAP算法简介

SNOMAP算法是Hall等提出的基于Landsat 5 TM数据的积雪识别算法。其物理基础体现在：① 积雪在可见光波段有较高的反射率，在短波红外波段有较强的吸收特征；② 大部分云在可见光波段有较高的反射率，在短波红外波段反射率依然很高。利用该特点，采用NDSI阈值提取方法可以很好的识别积雪。基于Landsat 5 TM数据的NDSI（Normalized Difference Snow Index）计算方法为：

其中， 为积雪在可见光波段的反射率， 为积雪在短波红外波段的反射率。针对Landsat 5 TM数据，对积雪和云反应敏感的第2波段（b₂，绿波段，0.52~0.60 μm）和第5波段（b₅，短波红外波段，1.55~1.75 μm）被用来计算归一化雪指数NDSI。通过在北美的验证，SNOMAP算法设定NDSI阈值为0.4，当像元NDSI值 ≥ 0.4时，该像元被定义为积雪。但积雪和水体在可见光和短波红外波段的反射特征相似，该阈值识别出的积雪中有水体存在。为了进一步识别积雪，利用近红外波段水体强吸收而积雪吸收弱于水体的特点，SNOMAP加入积雪识别的另外一个判别因子，b₄ ≥ 0.11，其中b₄为Landsat 5 TM数据的第4波段（近红外波段，0.76~0.90 μm）的反射率。这样，当满足NDSI ≥ 0.4且b₄ ≥ 0.11时，该像元被识别为积雪。

但考虑到植被覆盖度对积雪识别的影响，SNOAMAP算法中还充分考虑了植被的影响。当0.1≤NDSI≤0.4，但NDVI≥0.38，同时b₄ ≥ 0.11时，该像元同时被识别为积雪。其中基于基于Landsat 5 TM数据的NDVI（Normalized Difference Vegetation Index）计算方法为：

其中， 为积雪在近红外波段的反射率， 为积雪在红光波段的反射率。针对Landsat 5 TM数据，第4波段（b₄，近红外波段，0.76~0.90 μm）和第3波段（b₃，红波段，0.63~0.69μm）被用来计算归一化植被指数NDVI。

2.       SNOMAP算法在ENVI中的实现

SNOMAP算法在ENVI中利用决策树实现，决策树是基于遥感图像数据及其他空间数据，通过专家经验总结、简单的数学统计和归纳的方法，获得分类规则并进行遥感分类。分类的规则易于理解，分类过程也符合人的认知过程，最大的特点是利用多源数据。SNOMAP算法在ENVI中的实现步骤如下：

(1)      准备数据

Landsat 5 TM原始数据是没有经过定标的DN值，而SNOMAP算法需要Landsat 5 TM数据各个波段的反射率，因此运算前需要对原始Landsat 5 TM数据进行几何纠正、大气纠正和传感器定标。

(2)      NDSI、NDVI计算

利用ENVI主菜单﹣> Basic Tools﹣>Band Math，在公式输入栏中输入：

FLOAT (B2-B5)/( B2-B5)以及FLOAT (B2-B5)/( B2-B5)。NDSI和NDVI都是介于0-1之间的数据，需要在计算公示前加FLOAT，以生成浮点型的数据

(3)      决策树建立

利用ENVI主菜单﹣> Classification﹣>Decision Tree﹣> Build New Decision Tree，打开ENVI Decision Tree面板，进行规则定义与决策树的建立。

(4)      决策树执行

在ENVI Decision Tree面板中，选择Options﹣>Execute，执行定义好的决策树。

3.       SNOMAP算法在IDL中的实现

除了ENVI中已经做好的工具菜单，用户在实际的应用系统中也可以利用IDL语言以及ENVI的小部件函数，很方便的实现SNOMAP算法。原始Landsat 5 TM数据进行几何纠正、大气纠正和传感器定标后，SNOMAP算法实现步骤如下：

; 波段读取

Band2  = ENVI_GET_DATA(FID=TM_FID, DIMS=DIMS, POS=[1])

Band3  = ENVI_GET_DATA(FID=TM_FID, DIMS=DIMS, POS=[2])

Band4  = ENVI_GET_DATA(FID=TM_FID, DIMS=DIMS, POS=[3])

Band5  = ENVI_GET_DATA(FID=TM_FID, DIMS=DIMS, POS=[4])

; NDSI计算

NDSI_Add = Band2+Band5

        NDVI   = Band2-Band5

        index  = WHERE(NDSI_Add GT 0,nCount)

        IF(nCount EQ 0) THEN BEGIN

          STR_ERROR = 'Landsat NDSI failed!'

           RETURN

       ENDIF

NDSI[index] = NDSI[index]/NDSI_Add[index]

; NDVI计算

   NDVI_Add = Band4+Band3

       NDVI   = Band4-Band3

      index  = WHERE(NDVI_Add GT 0,nCount)

      IF(nCount EQ 0) THEN BEGIN

         STR_ERROR = 'Landsat NDVI failed!'

         RETURN

      ENDIF

    NDVI[index] = NDVI[index]/NDVI_Add[index]

   ; SNOMAP算法

   SnowMask = BYTARR(NS,NL)

   idx = WHERE(NDSI GT 0.4 and Band4 GT 0.11, nCount

   IF(nCount GT 0) THEN SnowMask[idx] = 1

   idx = WHERE(NDSI LE 0.4 and NDSI GT 0.1 and NDVI GT 0.38 AND Band4 GT 0.11, nCount)

   IF(nCount GT 0) THEN SnowMask[idx] = 1

4.       结果浏览