基于同质区域分割的高光谱图像混合噪声估计

孟 玉

(西北工业大学计算机学院，陕西 西安 710129)

0 引言

高光谱图像处理技术被广泛用于地质勘探、气象监测、环境保护等领域，然而人们在获取HS数据的同时，却不可避免地引入了噪声，如何有效地估计及去除噪声［1］对HS图像处理尤为重要。

空间光谱维去相关法［2］(Spectral and Spatial Decorrelation Method，SSDC)是一种经典的评噪方法，虽然该方法较稳定，但因强硬地将图像分为等大的小块，无法应对具有复杂纹理的HS图像，且评估结果误差较大。局部均匀块标准差法［3］(Local Homogeneous Block Standard Deviations Method，LHBSD)和同质区域划分光谱维去相关法［4-5］(Homogeneous Regions Division and Spectral Decorrelation Method，HRDSDC)则是基于同质区域进行评噪，对复杂纹理的HS图像依然适用，但HR的划分结果并不稳定，存在随机性，影响了噪声估计结果。Martin-Herrero［6］对HR划分方法进行了评估和改进，但改进后的算法仍要高深度递归，计算负担过重。

传统方法［7］只适用于估计信号无关(Signal Independent，SI)噪声，而电子器件敏感性的提升为HS图像引入了光电噪声，也即信号相关［8］(Signal Depend-ent，SD)噪声，导致传统方法难于进行估计。Qian［9-10］等人提出的三维非局部均值滤波方法评噪和去噪的效果显著，而Acito［11］等人也提出了一种自动化的评噪方法，通过矩和基于梯度的多变量迭代方法求解噪声模型的最大似然函数，继而计算SD和SI的方差。这些方法具有较高的稳定性和精确性，但复杂度较高，计算负担较大。

为了解决传统算法的混合噪声估计问题，同时为了提高估计算法的性能，提出一种基于HR分割的混合噪声估计算法。算法的主要思想是在HR分割区域内进行MLR估计，从而得到含较强空间和光谱相关性的纯信号估计，而该信号与原图像间的残差即为混合噪声。得到混合噪声后，可以利用信噪比和比例因子计算得到其内部详细的SD和SI的方差。该方法避免了复杂的预计算和后处理，不仅降低了算法复杂度，且能准确估计混合噪声及其内部参数。同时为了得到较好的HR分割并提升算法性能，提出一种最近邻划分的无监督HR分割算法，该算法结合HS的空间和光谱特性对图像进行扫描并分割，无需递归调用，算法高效稳定且结果要优于其它传统方法。

1 同性质区域分割

同性质区域［5］定义为在空间上相邻，且光谱特性相同或相似的大片连通区域，因此其内部的像元具有较强的空间和光谱相关性，若除去HR内像元间的相关性后仍存在误差，则该误差即由噪声引起。通过在HR内使用MLR去除强相关性，即可得到误差最小的残差，进而对残差进行估计从而得到噪声。理论上来看，HR的分割结果对噪声的估计有较大影响。

此处使用光谱角距离(Spectral Angle Distance，SAD)［12-13］来度量像元间的光谱相似性。若将单个像元的光谱特性描述为一个n维向量，n是像元光谱的波段数，则每一个向量具有确定的长度和方向，分别代表了像元的亮度和光谱特征，由此得到2个像元间的SAD为:

其中t和r分别代表目标像元向量和参考像元向量。光谱角θ的单位为弧度，取值范围为0～π/2，两个像元间的光谱角越小，它们的光谱特征就越相似。由式(1)很容易看出，当t和r有尺度变化时，光谱角θ依然不变，这种尺度不变的特性可以很好地满足刻画高光谱图像像元相似性的要求。

由于HR内的像元具有相似的光谱特征，故可将HR内所有像元向量的加权平均作为该区域的光谱特征。设参考区域R={r1，r2，…，rN}为包含N个像元的同质区域，目标像元t与区域R的SAD为:

同性质区域划分算法步骤如下:

(1)创建一个与原始图像同样尺寸的划分模板，并将左上角的像元标记为1，作为第一个区域。按照从左到右，从上到下的顺序寻找未标记的像元，转步骤(2)。

(2)寻找与目标像元上相邻和左相邻的参考像元(至少可以找到一个)。若只找到一个相邻的参考像元，且找到的参考像元是单个像元，转步骤(3);若参考像元属于某个同质区域，转步骤(4)。若找到两个相邻的参考像元，转步骤(5)。

(3)利用式(1)计算SAD，若SAD小于某一阈值μ则将目标像元与参考像元归为同一同质区域，并在划分模板上目标像元的位置标记上参考像元的标号;否则将目标像元划分为新的区域，划分模板上标记为当前最大标号加1，转步骤(6)。

(4)利用式(2)计算SAD，若SAD小于某一阈值μ则将目标像元归入参考的同质区域，并在划分模板上目标像元位置处标记上参考区域的标号;否则将目标像元划分为新的区域，划分模板上标记为当前最大标号加1，转步骤(6)。

(5)像元间或像元与区域间的SAD计算及处理参考步骤(3)和步骤(4)。当计算得到的2个SAD值都小于阈值μ时，将参与计算的两个参考像元或区域和目标像元合并为一个区域，并更改标记。否则将目标像元划分为新的区域，划分模板上标记为当前最大标号加1，转步骤(6)。

(6)继续转步骤(2)进行划分，寻找未标记像元，直到所有像元划分完毕，算法结束。

上述算法中的阈值μ是SAD的分割门限，需手动设置，关于它对算法的影响和讨论将在后续的实验部分给出。

2 混合噪声估计算法

此处假设HS图像仅受混合噪声的影响，且混合噪声是作为零均值的高斯随机过程存在于图像中的。一种普遍使用的混合噪声模型［11］可表示为:

其中σ2(n)为混合噪声的方差，μ(f)为图像均值，σ2(u)和σ2(m)分别为SD和SI噪声的方差。此处的目的就是估计得到除μ(f)之外的其它所有参数。

多元线性回归(MLR)［11］是一种普遍使用的估计方法，它利用了HS图像的高空间/光谱维相关性。假设HS图像共L波段，每波段共有M个像元，若将第 l波段的图像整理成行向量 Xl= ［，...，］，其中为第l波段第m个像元的DN值。再将某个像元m在其它L-1个波段上的DN值写成列向量则某波段所有像元的向量可以组成矩阵Fl=由此得到第l波段的可用信号估计:

其中Wl为第l波段的线性混合系数。

对HS图像的每个波段使用MLR进行估计，即可得到可用信号估计，其均值与图像均值μ(f)相等。再对残差rl=Xl－进行无偏估计，即可得到混合噪声的标准差σ(n)。

为了对SD和SI的标准差进行估计，可在给定图像混合噪声信噪比SNRl=θ的情况下引入SD和SI的比例因子α，有:

其中第l波段SD和SI噪声的信噪比为［11］:

结合式(3)和式(5)、式(6)即可解出:

综上所述，给出混合噪声的估计算法:

(1)首先对HS图像进行分割得到HR，为了消除小区域和杂点的影响，实验中可只保留像元个数大于某阈值的HR参与计算。

(2)对每个波段内的HR分别使用MLR进行估计，并将所有计算结果的统计平均值作为本波段的最终估计，即可得到每个波段的混合噪声参数σ2(n)和图像均值μ(f)。

(3)给出混合噪声的信噪比θ和比例因子α，通过式(7)计算得到每个波段的SD和SI噪声的标准差 σSD和 σSI，算法结束。

3 实验结果与分析

3.1 在仿真数据上实验

从美国地质勘探局(United States Geological Survey)的光谱库中获取5种物质，其光谱特征如图1所示。实验选取其中350 nm到2500 nm(间隔10 nm)共216个波段，将其两两混合(从左到右，从上到下依次为:明矾石+方解石、黑松+方解石、柳树+方解石、混凝土+方解石)，并依据HS线性混合模型产生200×200×216的HS仿真数据，如图2所示。区域划分模板如图3所示。实验在酷睿(R)E6500处理器(2.94 GHz)、2 GB 内存、32位 Windows 7操作系统、Matlab环境下进行。

图1 光谱特征图

图2 实验数据第60波段(940 nm)的仿真图像

图3 区域划分模板

图4 α=3时混合噪声(θ=30 dB)估计对比

图5 α取不同值(θ=30 dB)SD和SI噪声估计对比(实线为实际噪声，点为估计噪声)

表1 混合噪声标准差估计RMSE(θ=30 dB)

表2 SD和SI噪声标准差估计RMSE(θ=30 dB)

实验加入θ为30 dB，α分别为1/3、1和3的3组混合噪声到仿真数据中，然后将 SSDC［2］、LHBSD［3］、HRDSDC［5］和本文的方法进行对比，其中一组估计结果如图4所示。此处使用相对均方误差［11］(Relative Mean Square Error，RMSE)对估计方法进行定量评估:

图5和表2给出本文方法在3组数据间对SD和SI噪声的估计结果和RMSE对比，可以看出估计结果较优，且α仅影响到SD和SI在混合噪声中的比例关系，对估计结果的精度影响不大。

3.2 在AVIRIS数据上实验

选取美国加利福尼亚州萨利纳斯山谷的AVIRIS数据进行实验，其图像分辨率为512×217，空间分辨率为3.7 m，包含224个波段(已去除水扭曲波段)，其第100波段灰度图如图6(a)所示。

图6 (a)AVIRIS数据第100波段灰度图

实验对比了 HRDSDC［5］、Martin-Herrero［6］和本文方法的HR分割结果，并给出本文方法在不同阈值时的HR分割结果。对比图6(c)、(d)、(e)在(b)中的黑框部分可以看出，阈值相同时本文方法的分割结果较优。对比图6(e)、(f)、(g)可看出，在阈值过小时易出现过分割现象，阈值过大时则易造成错分少分的现象，可见阈值的选择对HR分割有较大影响。经过多次实验，最终选择阈值为0.05时可得到较满意的HR分割。

图7 阈值不同时多种方法间混合噪声估计结果对比

图7显示了取不同阈值时各种方法间的混合噪声估计结果对比。可以看出，HRDSDC受阈值影响较大，其改进方法受影响较小，而本文方法在不同阈值下都能获得较一致的估计结果，可见算法相当稳定，对HR分割结果的依赖性较低。

4 结束语

本文提出了一种基于同性质区域分割的高光谱图像混合噪声估计方法。在HR分割方面，解决了HRDSDC具有的随机性问题，同时在算法上避免了大量的递归调用。在评噪方面，解决了传统评噪方法不能估计SD噪声的问题，同时降低了算法复杂度。实验表明，本文提出的方法在稳定性和精确性上都要优于一般传统方法，具有一定的工程应用价值。

［1］张兵，高连如.高光谱图像分类与目标探测［M］.北京:科学出版社，2011:25-46.

［2］Roger R E，Arnold J F.Reliably estimating the noise in AVIRIS hyperspectral images［J］.International Journal of Remote Sensing，1996，17(10):1951-1962.

［3］高连如，张兵，张霞，等.基于局部标准差的遥感图像噪声评估方法研究［J］.遥感学报，2007，11(2):201-208.

［4］Martin-Herrero J.Hybrid object labelling in digital images［J］.Machine Vision and Applications，2007，18(1):1-15.

［5］Gao L R，Zhang B，Zhang X，et al.A new operational method for estimating noise in hyperspectral images［J］.IEEE Geoscience and Remote Sensing Letters，2008，5(1):83-87.

［6］Martin-Herrero J.Comments on“a new operational method for estimating noise in hyperspectral images”［J］.IEEE Geoscience and Remote Sensing Letters，2008，5(4):705-709.

［7］Rakwatin P，Takeuchi W，Yasuoka Y.Stripe noise reduction in MODIS data by combining histogram matching with facet filter［J］.IEEE Transactions on Geoscience and Remote Sensing，2007，45(6):1844-1856.

［8］Uss M L，Vozel B，Lukin V V，et al.Local signal-dependent noise variance estimation from hyperspectral textural images［J］.IEEE Journal of Selected Topics in Signal Processing，2011，5(3):469-486.

［9］Qian Yuntao，Shen Yanhao，Ye Minchao，et al.3-D nonlocal means filter with noise estimation for hyperspectral imagery denoising［C］//Proceedings of the 2012 IEEE International Geoscience and Remote Sensing Symposium.2012:1345-1348.

［10］Qian Yuntao，Ye Minchao.Hyperspectral imagery restoration using nonlocal spectral-spatial structured sparse representation with noise estimation［J］.IEEE Journal of Select-ed Topics in Applied Earth Observations and Remote Sensing，2013，6(2):499-515.

［11］Acito N，Diani M，Corsini G.Signal-dependent noise modeling and model parameter estimation in hyperspectral images［J］.IEEE Transactions on Geoscience and Remote Sensing，2011，49(8):2957-2971.

［12］何中海，何彬彬.基于权重光谱角制图的高光谱矿物填图方法［J］.光谱学与光谱分析，2011，31(8):2200-2204.

［13］Andreou C，Karathanassi V.A novel multiple endmember spectral mixture analysis using spectral angle distance［C］//Proceedings of the 2012 IEEE International Geoscience and Remote Sensing Symposium.2012:4110-4113.