基于二次散射的高光谱遥感图像光谱非线性混合模型

余先川，李建广，2，徐金东，张立保，胡 丹

(1.北京师范大学信息科学与技术学院，北京 100875;2.中国传媒大学信息工程学院，北京 100024)

0 引言

高光谱遥感图像(简称高光谱图像)是光谱分辨率在10 nm数量级范围内的图像[1]。由于其光谱分辨率高，能够提供更为丰富的地球表面信息，因此倍受国内外学者的关注和应用[2]。然而，高光谱图像的空间分辨率普遍较低，混合像元现象严重，这是传统的像元级图像分类和面积量测精度难以达到使用要求的主要原因。为了提高高光谱图像的分类精度，就必须解决混合像元的分解问题，使遥感应用由像元级达到亚像元级。光谱解混的目的就是得到每个混合像元内纯地物(端元)的光谱及每种地物所占的比例(丰度)。

光谱解混模型从本质上分为线性混合模型和非线性混合模型2种，区分二者的标准就是光子有没有在地物间发生多次散射。通常宏观尺度上的混合可认为是线性的，而微观尺度上的混合被认为是非线性的[3]。目前，基于线形混合模型的光谱解混理论方法已得到了广泛的研究与应用，但由于受实际地物间复杂关系以及大气散射的影响，光谱混合都是非线性的，这就使得应用传统的基于线性光谱混合模型的解混结果难以满足精度要求[4－10]。

为了在一定程度上解决线性混合模型解混精度不高的问题，本文根据地物散射特点，提出了一种基于二次散射的非线性混合模型——二次散射模型(secondary scattering model，SSM)，并应用该模型进行了解混实验，对实验结果进行了分析。

1 非线性混合模型

1.1 广义的非线性混合模型

为了克服线性光谱混合模型的上述不足，许多学者利用非线性模型对高光谱图像的光谱混合方式进行了描述。

鉴于线性混合是非线性混合在光学多次散射被忽略情况下的特例[11]，本文定义了广义的非线性光谱混合模型，其表达式为

式中:y=[y1，…，yL]T为 L波段的高光谱图像;f是一未知的非线性函数向量;M=[m1，…，mR]为端元光谱矩阵;a=[a1，…，aR]T为端元的丰度向量(R为端元的个数);n为0均值且方差为σ2的高斯噪声向量。

由于f函数是未知的，因此基于非线性光谱混合模型的光谱解混是个病态问题，需要给出合适且确定的非线性混合模型才能进行解混。为此，本文提出了基于二次散射的非线性混合模型即二次散射模型。

1.2 二次散射模型(SSM)

光子在进入传感器前一般会在地物间发生多次碰撞和散射，为了使计算简化，本文忽略了3次以上的散射，只考虑光子在进入传感器前与地物发生2次撞击的情况，提出二次散射模型，表达式为

1.3 模型推导

假设某一混合像元内只存在A和B 2种地物，且地物的空间位置及比例都未知;设2种地物的光谱向量分别为m1和m2，丰度向量分别为a1和a2，即 M=[m1，m2]，a=[a1，a2]T。

光谱混合方式分线性和二次散射2大类，如图1所示。

图1 地物散射情况Fig.1 Ground object scattering case

从图1可以看出，当光子入射到某种地物上后，将以概率c直接反射到传感器，以概率d(d=1－c)散射到其他地物上。以概率d散射的这部分光子又分别以概率a1和a2碰撞到地物A或B后再反射到传感器中。该混合像元得到的光谱为

式中:前6项分别对应图1(a)—(f)6种情况;第7 项为误差项。则式(3)可进一步化简为

1.4 二次散射系数的确定

光谱混合从形式上分为聚合式、整合式和致密式3种，如图2所示。多次散射大多发生在致密式分布的地物中。

二次散射模型是半确定性模型，需要首先确定二次散射系数d。二次散射系数能直观反映地物分布的复杂程度:地物分布越简单，二次散射系数越低(极限情况d=0，为线性混合);相反，地物分布越复杂，非线性程度越强，则二次散射系数越高。

图2 光谱混合的3种形式Fig.2 3 forms of spectrum mixture

假设某混合区由2种纯地物组成，如果某个像元只存在其中一种地物，可认为二次散射系数d=0(线性情况)。假设像元内有2种地物且均等分布，则地物的交界地带最多，可认为此种情况下二次散射系数最大(图3(a));如果2种地物不均等分布，比如地物A占10%，地物B占90%，则认为此种情况地物间存在交界，但交界地带相对较少，二次散射系数较低(图3(b))。

图3 地物分布复杂度Fig.3 Complexity of ground object distribution

地物丰度方差能表征地物的复杂程度，当像元中只存在一种纯地物，则丰度方差最大，对应的地物复杂度最低，即二次散射系数最小(d=0);当像元中地物分布均等时，则地物丰度方差为0，但地物分布最复杂，二次散射系数最大。可见，二次散射系数与地物丰度方差成反比。

如果忽略地物组分的颗粒度、粗糙度及透明度等微观因素，则二次散射系数可定义为

式中:d(i，j)∈[0，dmax]，dmax为最大的二次散射系数(本文设为0.5);var(i，j)为地物丰度方差，(i，j)为图像坐标;varmax为各混合像元中地物丰度方差的最大值。

2 方法实验及结果分析

2.1 应用模拟数据的实验及结果分析

2.1.1 模拟数据的生成

从ENVI USGS Spectral Libraries中找出天然碱、泻利盐和黝帘石3种地物光谱数据[12]。该光谱数据的波长区间为395～2 560 nm，光谱分辨率为10 nm，波段数为210个。分别根据光谱线性混合模型和二次散射模型生成大小为100像元×100像元且包含210个波段的高光谱图像，并增加不同信噪比的高斯白噪声。3种纯地物的光谱曲线如图4所示。

图4 3种纯地物的光谱曲线Fig.4 Spectral curves of 3 pure mineral

根据模拟数据生成的3种纯地物丰度图像如图5所示。

图5 模拟生成的3种纯地物丰度图Fig.5 Simulated abundance image of 3 minerals

将每种纯地物丰度图像变成一个1×10 000的行向量，然后将3个行向量组成一个3×10 000的地物丰度矩阵。3个地物的光谱数据矩阵为210行×3列，每一列对应一种纯地物的光谱向量，每种纯地物光谱向量有210个采样点。

首先，对模拟数据进行光谱线型混合:将光谱矩阵与地物丰度矩阵相乘，得到210×10 000的混合矩阵。该矩阵表示高光谱图像每1行为1个波段，共210个波段。图6为第100和150波段的图像，不同颜色代表不同反射率值。

图6 线性混合模拟生成的高光谱图像Fig.6 Simulated hyperspectral image based on linear mixture

根据式(1)对模拟数据进行光谱非线性混合。图7是二次散射系数d=0.3时模拟得到的第100和150波段的图像。不同颜色代表不同反射率值。

图7 非线性混合模拟生成的高光谱图像Fig.7 Simulated hyperspectral image based on nonlinear mixture

在实验中，对模拟得到的数据加入不同信噪比的噪声，信噪比计算公式为

式中var(s)和var(n)分别代表信号和噪声方差。

图8为非线性光谱混合后再加入50 dB噪声后的第100波段图像。

图8 非线性光谱混合后加入50 dB噪声的第100波段图像Fig.8 100thband simulated hyperspectral image based on nonlinear mixture(RSNR=50 dB)

2.1.2 实验方法及结果

将模拟得到的数据分别用全约束最小二乘法(fully constraind least squares，FCLS)[13]和非线性变换后的最小二乘法进行解混，得到估计的地物丰度数据。利用式(7)计算估计值的均方误差[14]，即

式中:yn表示各纯地物的实际丰度值向量;y^n表示各纯地物的估计丰度值向量;N为采样点个数，在本实验中N=60 000。

表1为在相同信噪比、不同二次散射系数情况下，基于线性和非线性模型解混得到的均方误差。

表1 相同信噪比(RSNR=50)、不同散射系数下的解混均方误差Tab.1 MSE of solved abundance at same RSNR(RSNR=50)but different scatter coefficient

由表1可以看出，随着二次散射系数的增大，即地物分布复杂情况的加重，基于线性混合模型解混得到的均方误差值增大。可见，地物分布越复杂，线性混合模型的适用性越差;而基于非线性混合模型的解混效果好得多。同时可以发现，随着二次散射系数的提高，非线性解混效果也变得更好。

表2为在相同二次散射系数(d=0.3)、不同信噪比情况下，基于线性和非线性模型解混得到的均方误差。

表2 相同二次散射系数(d=0.3)、不同信噪比下的解混均方误差Fig.2 MSE of solved abundance at same d(d=0.3)but different RSNR

由表2可以看出，不管是基于线性混合模型还是基于非线性混合模型，信噪比RSNR越低，解混效果越差。与线性解混结果相比，非线性的解混效果更优，提高的百分比随着信噪比的增大而增大，表明非线性模型抗噪性更好。

图9为在二次散射系数d=0.3、信噪比RSNR=40情况下，基于线性混合模型解混得到的地物A丰度图像及其误差图像。

图9 基于线性混合模型解混得到的地物A丰度图像(左)及其误差图像(右)(d=0.3，RSNR=40)Fig.9 Solved abundance image(left)of mineral A and it’s error image(right)based on linear mixture(d=0.3，RSNR=40)

图10 为在二次散射系数 d=0.3、信噪比为RSNR=40情况下，基于非线性混合模型解混得到的地物A丰度图像及其误差图像。

图10 基于非线性混合模型解混得到的地物A丰度图像(左)及其误差图像(右)(d=0.3，RSNR=40)Fig.10 Solved abundance image(left)of mineral A and it’s error image(right)based on nonlinear mixture(d=0.3，RSNR=40)

基于图9，10可以发现，基于非线性混合模型的解混误差明显小于基于线性混合模型的解混误差。

2.2 应用AVIRIS数据的实验及结果分析

2.2.1 数据的选取

选取美国内华达州南部Cuprite地区的AVIRIS高光谱数据[15]。图像大小为250像元×191像元，包含224个波段，波长区间为0.389～2.467μm。除去低信噪比和水气吸收波段(波段1～2，104～113，148～167，221～224)，最后可用数据为188个波段数据。图11为第100波段高光谱图像。

图11 Cuprite地区AVIRIS第100波段图像Fig.11 100thband image of AVIRIS of Cuprite district

Cuprite区主要分布6种矿物[16]，分别为Desert-Varnish(沙漠岩漆)、Alunite(明矾石)、Montmorillonite(蒙脱石)、Chalcedony(玉髓)、Kaolinite(高岭石)#1和Kaolinite(高岭石)#2。

2.2.2 实验方法及结果

利用顶点成分分析算法[17]提取6种矿物端元，基于线性混合模型、利用FCLS方法[13]得到的6种矿物丰度图像如图12所示。

图12－1 Cuprite地区6种矿物线性解混结果Fig.12－1 Solved abundance results of 6 minerals based on linear mixture in Cuprite district

图12－2 Cuprite地区6种矿物线性解混结果Fig.12－2 Solved abundance results of 6 minerals based on linear mixture in Cuprite district

利用式(5)确定基于该解混结果的地物复杂度，设定最大二次散射系数dmax=0.5，得到的地物复杂度分布如图13所示。

基于非线性混合模型(式(1))，使用上面得到的二次散射系数对每个混合像元进行解混，得到的解混结果如图14所示。

对比图12，14可以看出，基于二次散射模型的解混效果更好，图像上道路影像特征清晰(图14);而基于线性模型解混的道路影像特征则模糊(图12)。主要原因是道路与周围地物是主要边界地带，而边界地带容易发生二次散射。

图13 Cuprite地区地物分布复杂度图Fig.13 Complexity of ground object distribution in Cuprite district

图14 Cuprite地区6种矿物非线性解混结果Fig.14 Solved abundance results of 6 minerals based on nonlinear mixture in Cuprite district

图15 显示了在不同散射系数下(d=0～0.5)， 蒙脱石的分布情况。

图15 不同二次散射系数下的蒙脱石分布情况Fig.15 Solved abundance of Montmorillonite in different secondary scatter coefficient

从图15可以发现，基于二次散射模型解混的道路影像特征更加清晰连续，而基于线性混合模型(d=0)的道路影像特征则比较模糊，并出现道路影像局部中断现象。

3 结论

本文对高光谱遥感图像的光谱解混问题进行了研究，得出以下结论:

1)提出了广义的非线性混合模型。

2)针对线性混合模型解混精度不高的问题，提出了一种基于二次散射的非线性混合模型，即二次散射模型(SSM)。

3)基于地物丰度方差，定义了地物混合复杂度。

4)通过计算解混后的地物丰度与实际地物丰度的均方误差发现，与基于线性混合模型解混结果相比，基于二次散射模型的解混精度能提高近一个数量级。

本论文丰富了高光谱遥感图像的光谱解混理论，为光谱解混的进一步应用提供了理论基础。今后将针对不同区域地形地貌特点，对二次散射模型的二次散射系数做进一步研究，以期得到更精确的解混结果。

[1]Schweizer S M，Moura J M F.Hyperspectral imagery:Clutter adaptation in anomaly detection[J].IEEE Tran on Information Theory，2000，46(5):1855－1871.

[2]童庆禧，张 兵，郑兰芬.高光谱遥感——原理、技术与应用[M].北京:高等教育出版社，2006:166－169.Tong Q X，Zhang B，Zheng L F.Hyperspectral remote sensing[M].Beijing:Higher Education Press，2006:166－169.

[3]Keshava N，Mustard J.Spectral unmixing[J].IEEE Signal Processing Mag，2002，19(1):44－57.

[4]Craig M.Minimum volume transforms for remotely sensed data[J].IEEE Trans Geosci and Remote Sensing，1994，32(3):542－552.

[5]Nascimento J M P，Bioucas－ Dias J M.Does independent component analysis play a role in unmixing hyperspectral data?[J].IEEE Trans Geosci and Remote Sensing，2005，43(1):175－187.

[6]Jia S，Qian Y.Constrained nonnegative matrix factorization for hyperspectral unmixing[J].IEEE Trans Geosci and Remote Sensing，2009，47(1):161－173.

[7]Zortea M，Plaza A.Spatial preprocessing for endmember extraction[J].IEEE Trans Geosci and Remote Sensing，2009，47(8):2679－2693.

[8]Ghosh G，Kumar S，Saha S K.Hyperspectral satellite data in mapping salt－affected soils using linear spectral unmixing analysis[J].Journal of the Indian Society of Remote Sensing，2012，40(1):129－136.

[9]Liu K H，Wong E，Du E Y，et al.Kernel－ based linear spectral mixture analysis [J].IEEE Trans Geosci and Remote Sensing，2012，9(1):129－133.

[10]Deng Y B，Fan F L，Chen R R.Extraction and analysis of impervious surfaces based on a spectral unmixing method using Pearl River Delta of China landsat TM/ETM+imagery from 1998 to 2008[J].Sensors，2012，12(2):1846－1862.

[11]Iordache M D，Bioucas－Dias J M.Sparse unmixing of hyperspectral data[J].IEEE Trans Geosci and Remote Sensing，2011，49(6):2014－2039.

[12]RSI(Research Systems Inc.).ENVI user’s guide version 4.0[S].Boulder，CO 80301 USA，2003.

[13]Heinz D C，Chang C I.Fully constrained least squares linear spectral mixture analysis method for material quantification in hyperspectral imagery[J].IEEE Tran GRS，2001，39(3):529－ 545.

[14]Plaza J，Plaza A，Martinez P，et al.Nonlinear mixture models for analyzing laboratory simulated－ forest hyperspectral data[J].Proc SPIE Image and Signal Processing for Remote Sensing IX，2004，52(38):480－487.

[15]AVIRIS data of Cuprite distuict，America.http://aviris.jpl.nasa.gov/html/aviris.freedata.html.

[16]Vane G，Green R，Chrien T G，et al.The airborne visible/infrared imaging spectrometer(AVIRIS)[J].Remote Sensing of the Environment，1993，44(2/3):127－ 143.

[17]Nascimento J，Dias J B.Vertex component analysis:A fast algorithm to unmix hyperspectral data[J].IEEE Trans Geosci and Remote Sens，2005，43(4):175－187.