一种反演积雪面积的范数最小二乘算法

段金亮，张瑞，2，李奎，庞家泰

(1.西南交通大学 地球科学与环境工程学院，成都 611756；2.西南交通大学 高速铁路运营安全空间信息技术国家地方联合工程实验室，成都 611756)

0 引言

积雪是地球表面覆盖的重要组成部分，是地球表面最为活跃的自然元素之一，也是冰冻圈中地理分布最广泛、季节与年际变化最显著的部分[1]。积雪的特征包括积雪面积、积雪分布、雪深、雪粒径等，其中积雪面积是全球能量平衡、气候、水文、生态模型以及积雪定量遥感中的重要输入参数之一[2]。当前，积雪面积的反演主要集中在光学遥感影像上，方法主要有阈值法[3]、监督分类法[4]、光谱混合分析法等。阈值法具有简单快速的特点，但是不同影像的阈值不固定，每景影像都要重新选择阈值，通用性差；监督分类法虽然有一定的鲁棒性，但是需要大量的样本作为支撑条件，这在生产实践中很难满足。同时，阈值法和监督分类法都是从像元尺度来获取积雪面积，由于积雪覆盖范围具有较强的时空变异性特征，以及受影像空间分辨率的限制，积雪通常与土壤、岩石、植被等混合成混合像元[5-6]。如果单独从传统的像元尺度分类算法来考虑，很难解决积雪与其他地物混合的问题，往往难以获得更有效、精度更高的积雪面积产品，所以从亚像元尺度考虑的光谱混合分析法在反演积雪面积上很有优势[7-8]。全约束最小二乘法(fully constrained least squares,FCLS)是遥感影像处理中常用的光谱混合分析模型[9]，其端元光谱库被假定在整景遥感影像中是不变的，然后用其来解混整景遥感影像。但是在实际情况中，端元光谱存在异物同谱和同谱异物的现象，同时端元是沿着影像进行变化的，这种变化导致了端元变化。这种端元变化减低了遥感影像的解混精度，端元变化已经被国内外学者证实是丰度估计的主要误差源[10-11]。

针对遥感影像中端元变化的问题，国内外学者提出了多种光谱混合模型[9-14]，其中比较经典的是多端元光谱混合分析法(multiple endmember spectral mixture analysis，MESMA)。虽然MESMA算法能够很好地处理端元变化问题，但是MESMA算法的运算效率存在一定的局限。同时，由于积雪的时空变化较大，造成了端元变化也随之较大，同时还与其他地物进行混合。基于以上考虑，本文对常规的线性混合模型全约束最小二乘法进行改进，提出范数最小二乘算法(norm least squares,NLS)，通过引入范数来降低积雪端元在可见光波段的相对差异，以减弱端元变化的影响，从而提高反演积雪面积的精度。

1 范数最小二乘算法

在遥感影像的特征提取和分类中，常用的光谱分析模型是线性混合模型(linear mixed model)。其通过提取研究区典型地物的光谱作为端元光谱库，然后应用于线性模型上获取各个地类的丰度值。常用的线性混合模型是全约束最小二乘法[9]，表达如式(1)所示。

(1)

式中：y∈1×L为混合光谱向量；e∈1×L为端元向量；a为丰度值；ε∈1×L为噪声向量；M是端元数目；L为波段的数目。其中，丰度值满足非负、和为1 的限制，如式(2)所示。

(2)

因为光学影像在成像时受到地物物理特性的变化、光照条件和大气环境的影响，存在同物异谱和同谱异物的情况，从而产生端元变化，导致光谱混合分析算法的精度降低。同时，积雪覆盖的MODIS影像的端元变化较大，选择合理的端元存在困难。此外，若应用于解混的端元光谱库的相关性过高，则遥感影像的解混精度也会降低，从而降低反演积雪面积的精度。此时，常规的FCLS算法就有一定局限性，特别是针对积雪这类光谱变化性较大的地物，常规的FCLS算法很难解决积雪的端元变化，导致其算法的反演积雪面积的精度不高。故此处利用范数最小二乘法来减弱端元变化，改写式(1)为式(3)至式(4)所示。

(3)

(4)

式中：x∈1×L为被用于‖·‖p处理的值。

为减弱端元变化对解混结果的影响，这里每个端元光谱除以它们的2范数，这样压缩了数据空间，缩小了端元间相对差异，减弱了端元变化的绝对量，从而减弱了端元的变异性，提高了解混精度，同时保留了它们在高维特征空间的相对位置。

2 实验与分析

2.1 数据来源及预处理

MODIS影像的光谱覆盖度广(36个光谱波段)，空间分辨率适中，且重访周期短，获取数据便捷，适合用于大范围的积雪面积的提取。为此本文对MODIS数据进行反演实验，验证NLS算法反演积雪面积的可行性，同时为验证其反演的积雪面积精度，选择同一时间段的TM影像进行精度评价。TM影像编号为：LT51450362009307KHC00，其成像时间为：2009年11月3日。该影像云量小、地物类型单一，有大范围的积雪覆盖，适合积雪面积反演算法的研究。

首先，将HDF格式的MODIS影像转换成ENVI的标准格式，重新投影MOD09GA产品，使其与TM影像的坐标系统保持一致；其次，按照波长大小对波段进行升序排列重新合并。由于MODIS影像存在大量的坏像元，同时第5波段存在明显的条带，为了减弱它们的影响，采用空间与光谱结合的去噪算法[15]，其参数设置如下：[0.1，0.2，0.2，40]，分别对应λ，u，v，MaxIter[16]。然后，对TM影像进行预处理，去除大气误差获取地物反射率，其处理过程参考文献[17]。

因为无法获取研究区真正的积雪面积，此处采用空间分辨率高的TM影像制作研究区的积雪覆盖度的参考值，通过几何配准让MODIS影像与TM影像处于同一参考坐标系下，然后用TM影像去裁剪MODIS影像得到研究区域如图1所示(波长对应于TM影像的7、4、3 波长)。

图1 研究区的MODIS和TM影像

积雪面积验证数据需要从分类后的TM影像中获得。首先，对TM影像进行辐射定标和大气校正后，结合天地图和TM影像目视判读研究区地物主要包括积雪、裸岩、水体、土壤和阴影5个大类；其次，通过对TM影像进行目视判读，同时参考文献[18]选择分类样本和验证样本；再次，利用ENVI软件中的支持向量机算法(support vector machine algorithm，SVM)进行分类，其参数设置参考文献[19]；然后，使用验证样本对分类结果进行精度验证，其分类结果的总体精度达到88.5%，Kappa系数达到0.86。获得TM影像的分类结果后，通过目视判读改正分类错误的积雪像元，判读时参考天地图的高分辨率影像逐一改正。最后,把经过改正后的积雪类别保存出来。由于TM影像的空间分辨率(30 m)与MODIS影像的空间分辨率(500 m)不一样，为满足后期反演积雪面积算法的精度验证实验需要，对分类后的积雪类别影像进行像元聚合，使积雪类别影像与MODIS影像的空间分辨率保持一致，从而获得研究区的积雪覆盖度的参考值。研究区的积雪覆盖度参考影像如图2所示。

图2 参考的积雪覆盖度

光谱混合分析算法在反演积雪面积前，需要提前构建端元光谱库，这里使用空间分块自动端元提取算法[20]从研究区的MODIS影像上提取积雪、裸岩、水体、土壤和阴影5类端元，然后剔除异常的端元，最后每个类别得到15个端元。由于本文主要反演积雪面积，所以展示积雪光谱曲线如图3所示，图中积雪的光谱曲线变化较为明显，主要是可见光谱波段，最大变化达到0.28，近红外波段变化相对较小。通过对光谱进行范数计算后，如图4所示，引入范数计算，缩小端元变化的绝对量，减少端元间的相对差异，从而在一定程度上减弱积雪的光谱变化，提高反演的精度。

图3 积雪端元光谱库

图4 求解范数后的积雪光谱库

2.2 精度验证

对FCLS、NLS和MESMA算法反演积雪面积的时间效率进行对比分析，在电脑处理器为Intel(R)Core(TM)i5-4210M CPU@2.60 GHz、运行内存为8.00 GB和64位windows7的操作系统下，对7 000行8 000列大小的数据进行反演计算，最后得到3种模型的运行时间，如表1所示。

表1 3种模型的评价指标

为验证FCLS、NLS和MESMA算法反演的积雪面积的精度，使用TM分类影像作为参考值，同时使用绝对差异值、均方根误差(RMSE)、相关系数(R)、总的积雪面积(TSCA)指标来进行精度评价，其FCLS、NLS和MESMA算法的精度评价指标对比结果如表1所示。针对FCLS、NLS和MESMA算法反演出的结果，对积雪覆盖度进行求和与归一化，其中NLS、MESMA和FCLS反演的积雪覆盖度如图5所示，其与参考值的绝对差异如图6所示，其绝对差异值分段统计结果如表2所示。

图5 NLS、MESMA和FCLS算法反演的积雪覆盖度

图6 NLS、MESMA和FCLS算法的绝对差异图

通过将图5(a)和图5(b)与参考图对比，发现图5(a)十分接近参考影像，同时图6(a)和图6(b)的对比，进一步说明图5(a)的结果十分接近参考值，再结合表1的3种评价指标来看，NLS的评价指标参数都要优于FCLS。而在表2中，FCLS的绝对差异值分布在0.0～0.1区间的只占有44.41%，且在0.2～0.4区间都超过了10%，而NLS的绝对差异值主要分布在0.0～0.1区间，其次分布在0.1～0.2区间，且0.0～0.1区间的值高达63.52%。综上，可以得出NLS反演积雪面积明显优于FCLS，说明引入范数减弱了端元的变化，从而提高了算法的精度。

表2 3种算法的绝对差异值

通过分析表1和表2中MESMA和NLS的评价参数，可以发现其RMSE、R、TSCA的值相对一致，同时其与真值的绝对差异值在各个区间基本相似，可以证明NLS的反演积雪面积效果跟MESMA反演积雪面积的效果非常吻合。通过分析评价指标的值，可以得出，85.71%的绝对差异值都集中在0.0～0.3范围，且在0.0～0.1的范围高达62.82%，最重要的是TSCA更是高达0.95以上，从而可以说明MESMA和NLS反演的积雪面积与真值十分接近，且可信度高，可以用MESMA和NLS来反演MODIS影像的积雪面积。虽然MESMA和NLS都可以用反演积雪面积，但是通过表2的时间指标可以看出来，NLS的时间效率是MESMA的几十倍，所以在这点上NLS要优于MESMA。

综上，NLS反演积雪面积的效果要远优于FCLS，且与MESMA的结果保持一致，都能够很好地来反演积雪面积，且与真值十分接近，说明减弱端元变化能够显著提高反演积雪面积的精度。同时验证了MESMA算法通过直接将同类物质的所有光谱参与光谱分解来减弱端元变化，但是增加算法的时间复杂度，消耗大量的时间。而NLS通过引入范数的方法同样成功消除了端元变化的影响，且反演效果很理想，时间效率比MESMA高几十倍。

3 结束语

本文基于MODIS影像反演积雪面积提出了一种范数最小二乘法，来解决端元变化和积雪光谱的异质性问题，从而提高反演积雪面积的精度。实验结果表明，NLS在反演MODIS影像上的积雪面积方面的可行性强且应用前景广阔。本文得出以下结论。

1)NLS反演积雪面积的精度远远高于FCLS，引入范数减低端元变化和光谱异质性的影响，提高了反演精度。

2)NLS反演的积雪面积和MESMA保持高度一致，说明其解决端元变化的效果跟MESMA一样，且时间效率比MESMA高，其更适用于反演MODIS影像的积雪面积。

NLS算法对于高效反演MODIS影像的积雪覆盖度有较好的效果，但仍有改进空间。此算法没有考虑空间领域信息、波段噪声权重对反演积雪面积的影响，同时其运行效率还有待进一步提高[21-22]。