MCR-BANDS用于高光谱成像分辨中旋转模糊度的评估
邵常艳, 张欣*, 张卓勇*
首都师范大学化学系, 北京 100048
*通讯联系人 E-mail: xinzhang@cnu.edu.cn; zhangzhuoyong@cnu.edu.cn

作者简介: 邵常艳, 女, 1992年生, 首都师范大学化学系硕士研究生 E-mail: 2160702075@cnu.edu.cn

摘要

自建模曲线分辨用来将双线性光谱数据矩阵分解成具有明确物理或化学意义的曲线, 一方面反映了复杂体系中各个主成分对应的纯光谱, 同时也能够解析得到其对应的相对浓度。 这一方法在高光谱解析中充分发挥了其优势, 成为了高光谱分析中的重要方法之一。 然而, 双线性结构数据多元分辨模型在约束不充分的条件下往往不能获得唯一解, 这一问题是由顺序模糊、 尺度模糊和旋转模糊引起, 其中旋转模糊最难消除。 在高光谱成像中, 如果不能明确浓度分布情况, 将导致在成像领域难以确认目标的准确位置或形貌。 为了充分了解旋转模糊带来的影响, 评估非唯一解情况下可行解的范围, 并进一步为实际应用中需要解析获得的每个主成分的纯光谱或波谱信号, 以及其对应的浓度信息提供客观的评估依据, 在以往的研究中, 研究人员分别使用网格法、 蒙特卡洛法进行抽样, 以计算曲线分辨结果中可行解的范围, 也有科研人员使用了几何多边形内部和外围面积的形式表示结果, 但是这些结果往往面临运算时间过长, 或者无法实现高维可视化而不适用于大于四个主成分的数据等问题, 而且通常这些方法很难将曲线分辨过程中施加的除非负约束之外的其他约束方法加以利用, 导致可行解范围计算不准确。 为了解决以上问题, 采用MCR-BANDS对MCR-ALS(多元曲线分辨-交替最小二乘)分辨获得的结果进行了旋转模糊程度的评估, 并将其应用到遥感高光谱成像的解析中。 首先以美国地质勘探局矿物光谱库中的纯光谱为基础的模拟数据集对MCR-ALS和MCR-BANDS的结果进行了评测, 在模拟数据中能够方便地控制噪声的影响, 控制选取主成分之间的纯光谱差异、 仿照真实环境中浓度渐变特征等因素, 考查了特定条件下MCR结果中旋转模糊的水平。 随后为了证实所用方法的可行性, 进一步采用MCR-ALS分析了机载可见/红外成像光谱仪(AVIRIS)获得的遥感高光谱图像数据, 并首次采用MCR-BANDS对MCR-ALS的分辨结果的旋转模糊进行了分析, 实现了对遥感高光谱数据成像浓度分布的受到旋转模糊影响的可视化表示。 可以发现真实解和MCR-ALS获得的可行解均在MCR-BANDS计算得到的可行解范围之内。 结果表明, MCR-BANDS方法基于最大和最小的信号贡献对旋转模糊的范围进行计算, 能够适用于不同主成分的体系中, 并且完美对接MCR-ALS中使用到的诸如非负、 单峰、 封闭和选择性约束等。 MCR-BANDS的分析结果可以为MCR-ALS的解析结果提供相应的旋转模糊水平估计, 有利于对MCR-ALS结果的解释; 在充分约束条件下, 能够有效减少甚至消除旋转模糊对MCR-ALS分辨结果的影响, 为精确确定遥感高光谱中解析得到的目标物位置提供了客观的范围。

关键词: MCR-BANDS; MCR-ALS; 旋转模糊; 遥感; 高光谱
中图分类号:O657.39 文献标志码:A
Evaluation of Rotation Ambiguity by MCR-BANDS on Hyperspectral Imaging
SHAO Chang-yan, ZHANG Xin*, ZHANG Zhuo-yong*
Department of Chemistry, Capital Normal University, Beijing 100048, China
Abstract

Self-modeling curve resolution can resolve the bilinear spectral datasets as the profiles of pure signal and their contributions which can be explained easily in physical and chemical meaning. With the advantage that their results can provide the pure spectra and the corresponding relative concentrations in the analyzed complex systems, MCR methods have been widely applied in the analysis of hyperspectral imaging. However, when the constraints applied are not strong enough, multivariate curve resolution models for bilinear data always suffer the problem of order ambiguity, scale ambiguity and rotation ambiguity which induce non-unique solution. Rotation ambiguity is the most difficult to be removed. To investigate the level of rotation ambiguity and provide the range of feasible solutions, in the published works, the researcher used grid search or Monte Carlo random sampling to display some feasible solutions fulfill the bilinear model under certain constraints. In this way, the concentrations and pure spectra results resolved by MCR can be better explained for their application by providing a range of feasible solutions. Polygons projected by feasible solutions based on geometry were also employed to illustrate the feasible solutions and the extension of rotation ambiguity. These methods normally are time consuming and cannot be used in the system with more than four components. More importantly, they cannot apply different constraints based on the properties of the analyzed samples, except non-negativity. In this work, we applied MCR-BANDS to evaluate the level of rotation ambiguity for resolved by MCR-ALS on the remote sensing hyperspectral imaging. In the first part, the mineral spectra selected from United States Geological Survey Committee were used for simulating a hyperspectral imaging dataset. In the simulated dataset, the noise level can be controlled and the differences of the spectral features between different components were easily identified. The concentrations of different components were simulated the real conditions, with gradual changes in the space. The rotation ambiguity was evaluated in the simulated data by using MCR-BANDS. To better explain the application of MCR-BANDS, this method was used to analyze a remote sensing dataset collected by Airborne Visible Infrared Imaging Spectrometer (AVIRIS), and the affection of rotation ambiguity on different components was visually displayed as concentration distributions in maps. The results show that MCR-BANDS can provide the level of feasible solutions of MCR-ALS by using maximum and minimum signal contribution functions (SCCF). This method can be applied to the system with any number of components, and can use almost all of the constraints which are chosen in MCR-ALS, like non-negativity, unimodality, closure, selectivity/local rank etc. The concentration distribution results from maximum and minimum SCCF are helpful to locate the specific targets in the remote sensing hyperspectral imaging.

Keyword: MCR-BANDS; MCR-ALS; Rotation ambiguity; Remote sensing; Hyperspectral imaging
引 言

自建模曲线分辨(SMCR)在1971年由Lawton和Sylvestre[1]提出, 用来将双线性两组分光谱数据矩阵分解成物理或化学意义明确的曲线, 例如在最小二乘非负约束条件下曲线分辨的方法可将仪器测量获得的数据分解为浓度曲线和光谱吸收曲线。 这一方法在高光谱解析中充分发挥了其优势, 成为高光谱分析中的重要方法之一。 然而, 很多双线性模型在约束不充分的情况下, 得到的解并不是唯一的, 如果没有进一步的约束, 此类方法只能得到纯组分曲线的可行解区域中收敛时对应的某一个可行解。 而与此相对应的特征值分解方法(SVD)通常情况下都可以获得唯一解, 然而此方法在化学数据的分析中, 往往无法给出便于解释其中理化信息的分解曲线。

对于大部分SMCR方法, 满足双线性模型

D=CST+E

式(1)中D(n× m)代表经仪器测量获得的原始光谱数据矩阵, n代表样品个数, m代表波谱数据; C(n× k)为纯组分的浓度分布的矩阵, k代表选择主成分个数; S(k× m)为纯组分的光谱矩阵矩阵, E为误差矩阵。

存在一个可逆矩阵T可以将上式写为D=CTT-1ST+E, 令Cnew=CT, SnewT=T-1ST, 则D=Cnew SnewT+E。 而所有可行解在数学意义上均能同等水平地拟合原始数据。

继Lawton和Sylvestre工作后对SMCR方法的优化分成了两个方向: (1)进一步发展可用来评估非唯一解的方法来获得多余两个组分的可行解区域。 (2)利用迭代方法和充分约束来得到唯一解。 这两个方向的研究从1972年以来不断都有文献报道[2, 3, 4, 5, 6]。 Lawton和Sylvestre最早解决了两个主成分体系中确定可行解范围的问题。 但是这个方法很难在多于两个主成分的体系中拓展。 有研究[2, 7]则发展了Lawton Sylvestre方法, 在三个主成分体系中实现了相同的约束, 然而这一体系很难理解, 大众接受程度低, 这就解释了为什么很多化学计量学家倾向于开发一些SMCR的近似方法[8]。 受这种趋势的影响, Rajkó 和Istvá n[4]重新研究了Borgen的方法, 在之后的研究中, 他们在Borgen的方法的基础上作了进一步的发展, 使用几何图形的方式清楚地揭示了可行解范围的特征, 以多边形内部或外部区域范围的形式表现了出来。 但是这一几何方法难以实现单峰约束和选择性约束, 同时也仅能在一定的噪声范围内实现, 这样的条件制约了该方法的实际应用。 Leger和Wentzell使用了蒙特卡洛方法以寻找可行解的范围, 理论上是可以适用于任意主成分数量的体系中的, 在此方法中, 他们使用了随机的方式生成T矩阵以代表其中不同的可行解[9]。 正因为如此, 想要获得精确地可行解范围, 需要很长的计算时间。 近年来, Miron等提出了通过将极化分析应用于多组分混合体系的方法并以此分析了拉曼光谱, 然而此方法也只能使用非负约束[10]。 Hamid Abdollahi等的工作中将三个主成分体系的几何表示可行解的方法扩展到了四个主成分的体系中[11]。 伴随着主成分个数的增加, MCR中的可行解范围也变得更难计算, 而且其可行解范围也会进一步扩大, 以几何多边形的形式表示变得尤为困难。 Sawall等在之前的基础上开发了FAC-PACK工具箱[12], 同时实现了通过某些已知条件进行特定成分约束的目的, 从而实现了根据已知体系中的特定物质信息降低最终多元曲线分辨方法结果中旋转模糊程度的目的。 此方法可以输入特定物质成分的光谱信息, 通过计算可以减少浓度曲线中的旋转模糊, 反之亦然。 但是此方法到目前为止最多也仅适用于四个主成分的分析体系。

Tauler在Gemperline的方法的基础上发展出了Multivariate curve resolution-bands(MCR-BANDS)方法, 基于最大和最小的信号贡献对旋转模糊的范围进行了计算, 能够适用于不同主成分的体系, 并且完美对接Multivariate curve resolution-alternating least squares(MCR-ALS)中使用到的诸如非负、 单峰、 封闭和选择性约束等[13]。 在以往的研究中主要报道了一些简单的模拟数据中的应用, 并与其他方法的性能进行了对比。 本文将此方法应用到了高光谱成像的实际数据中, 为高光谱中解析主成分对应的浓度分布提供了客观评价依据。

1 实验部分
1.1 方法

MCR-BANDS可以用于计算MCR-ALS模型在特定约束下相关联的旋转模糊度的程度[14]

D=CinicSinicT=CinicTminTmin-1SinicT=Cmin, kSmin, kT=CintTmaxTmax-1SinicT=Cmax, kSmax, kT(2)

式(2)中计算得到的CmaxCmin是MCR可行解中能够获得的最大和最小的相对浓度, SmaxSmin则是对应于第k个主成分的最大和最小相对浓度时分辨得到的相应纯组分信号。 MCR-BANDS的工作就是, 优化获得对应于最大和最小浓度信号时相对应的T矩阵。 该方法基于信号分量贡献函数(信号分量贡献函数SCCF)的最大值和最小值, 信号分量贡献函数SCCFk定义为

SCCFk=ckskTCST(3)

式(3)中分式的分子部分对应的是经MCR-ALS分辨得到的第k个主成分对应的纯浓度和纯信号(光谱、 质谱等)可解释的信号, 分母部分对应的是MCR-ALS分辨得到的整个数据中所有主成分的信号总量。 MCR-BANDS通过计算Δ SCCF(最大和最小SCCF值之差)表示旋转模糊的程度, 估算可行解的范围。 每个主成分对应的Δ SCCF值均在0~1之间, 越接近于零说明可行解范围越小, 越接近于1说明可行解范围越大。 当Δ SCCF等于零时, MCR-ALS可以获得唯一解。

1.2 数据

利用模拟数据对MCR-BANDS在高光谱的应用进行了试验。 该数据是模拟2007年USGS(美国地质勘探局, United States Geological Survey, 简称USGS)数据库中光谱图, 此数据库中包含了矿藏、 植被、 土壤等不同地表物质的光谱。 数据集模拟了21× 21个像素的三组分混合体系的成像数据, 共计441条光谱, 每条光谱有224个变量, 经三个主成分浓度(C441× 3)与光谱( S3×224T)相乘获得模拟数据D441× 224。 并在此数据中加入了相当于最大信号响应的0.5%且均值为0的白噪声(符合Gaussian分布)。

另一个数据为真实高光谱成像实验数据, 采集于1997年的内华达州(Nevada)黄铜矿的遥感成像数据, 通过机载可见/红外成像光谱仪(AVIRIS)收集测定(NASA, Jet Propulsion Laboratory, http://aviris.jpl.nasa.gov/data/free data.html), 此数据为开放高光谱数据, 可免费作为科学研究中不同方法的相互对比。 该数据在光谱维度具有在400~2 500 nm波段的224个变量。 数据集的大小为614× 512像素。 为了更好地分析其中的理化信息, 删除一些具有水信号反射过高的光谱变量。

2 结果与讨论

模拟数据中共包含三个主成分, 每种主成分的光谱各不相同, 在整个波段范围内有些波段的吸收又有些相似之处。 每种主成分对应的浓度分布也具有各自的形状。 有些区域中同时出现了多个主成分, 而有些区域中仅有一种成分。 如图1所示, 模拟数据的真实浓度分布为1(b)中所示的三个不同图形。 红色代表浓度高的区域, 蓝色表示浓度低的区域。 首先通过主成分分析能够确定体系中共存在三个主成分。 通过MCR-ALS仅在施加非负约束的情况下, 模型收敛时可以获得一种可行解, 并将其按照像素折叠后即可获得相应区域的相对浓度分布, 即图1(c)。

图1 模拟数据中通过MCR-ALS和MCR-BANDS获得的浓度分布对照Fig.1 Comparison of concentration distribution obtained by MCR-BANDS with those obtained by MCR-ALS

MCR-BANDS对MCR-ALS模型获取的每个分辨的纯组分与其相应光谱的相对信号贡献的最大值和最小值进行估计, 获得的浓度分布分别对应图1(a)和图1(d), 在对应不同SCCF值得浓度分布中, 区域形状有明显的差异, 但是基本位置都在真实浓度分布附近, 由此可见即便是当前的结果不是唯一解, 依然能够为浓度分布提供可靠的依据。 MCR-ALS在采用非负约束条件下的结果与真实值已经非常接近。 MCR-ALS的优化结果和真实的浓度分布则都位于MCR-BANDS的最大最小值区间中。

在MCR-ALS的分辨结果中, 既包含了每个主成分的浓度分布, 还有每个主成分对应的特征光谱(红色曲线)。 MCR-BANDS同样可以获得相应的光谱结果(蓝色曲线)。 图2中分别给出了三个主成分的真实光谱、 MCR-ALS收敛时获得的一种可行解对应的光谱, 以及MCR-BANDS解出的SCCF最大和最小值情况下分别对应的光谱曲线。 实际上可行解的光谱曲线的集合将是一种三维的条带状曲面, 所以在部分的波段直接通过二维平面观察, 有些真实光谱和通过MCR-ALS获得的光谱超出了MCR-BANDS结果的范围, 如果通过三维曲面表示, 它们将仍然位于MCR-BANDS结果中间。

图2 模拟数据中通过MCR-ALS(红色)和MCR-BANDS(蓝色, 分别对应最大和最小SCCF值)获得的三个主成分光谱比较
(a): 第一主成分光谱; (b): 第二主成分光谱; (c): 第三主成分光谱
Fig.2 Comparison of pure spectra obtained by MCR-BANDS (Blue lines, corresponding to maximum and minimum SCCF) with those obtained by MCR-ALS (Red line)
(a): The spectra of the first component; (b): The spectra of the second component; (c): The spectra of the third component

另一个结果就是对内华达地区黄铜矿的遥感高光谱成像的解析。 黄铜区有少量植被覆盖, 地表主要为暴露的岩石, 包括以关键标志性矿物的出现为特征的变化带。 为了更好地分析处理数据, 删除一些具有严重水吸收导致的高反射区域中的光谱波段。

MCR-ALS结合非负约束对该地区高光谱成像的研究曾有报道[15], 在此高光谱成像数据中该地区有五种不同的矿物质浓度及其对应的纯光谱曲线得以分辨出来。

在此基础上仅采用非负约束进一步采用MCR-BANDS方法获得了双线性模型, 在旋转模糊存在的情况下, 最高和最低的浓度分布曲线及其分别对应的光谱曲线, 浓度曲线经折叠重构, 得到了五个主成分分别对应的最大和最小浓度分布图(如图3)。 五个主成分对应的Δ SCCF值分别为0.13, 0.21, 0.48, 0.31和0.05。 此数值越接近于0, 则可行解范围越小。 并且从结果可以看出, 对于同一个体系的不同主成分, 其可行解范围也有很大差别。

图3 内华达州黄铜矿的高光谱成像数据MCR-ALS和MCR-BANDS 分析获得的五个主成分分别对应的相对浓度分布
(a): SCF最大值对应; (b): MCR-ALS结果对应; (c): SCCF最小值对应
Fig.3 The concentration distribution results of the five components obtained by MCR-ALS and MCR-BANDS on the remote sensing dataset of AVIRIS data from Nevada
(a): MCR-BANDS (MAX); (b): MCR-ALS; (c): MCR-BANDS (MIN)

3 结 论

MCR-BANDS能够判断双线性模型遇到的旋转模糊的程度。 结果表明, MCR-ALS可用于解析高光谱成像区域中主成分的纯光谱及其对应的浓度分布图像。 MCR-BANDS能够适应于不同主成分数目复杂混合物体系, 可以方便地施加对应于MCR-ALS分辨过程中使用过的约束方法, 以用于准确估计评测特定约束条件下双线性模型分辨结果中可能存在的旋转模糊的程度。 在高光谱解析中, 旋转模糊的最大最小值能够以浓度分布的图像形式表示, 为提高高光谱中解析得到的主成分对应的浓度分布提供了客观评价依据。

The authors have declared that no competing interests exist.

参考文献
[1] Lawton W H, Sylvestre E A. Technometrics, 1971, 13(3): 617. [本文引用:1]
[2] Tauler R, Barceló D. TrAC Trends in Analytical Chemistry, 1993, 12(8): 319. [本文引用:2]
[3] de Juan A, Tauler R. Analytica Chimica Acta, 2003, 500(1): 195. [本文引用:1]
[4] Rajkó R, István K. Journal of Chemometrics, 2005, 19(8): 448. [本文引用:2]
[5] de Juan A, Tauler R. Journal of Chromatography A, 2007, 1158(1): 184. [本文引用:1]
[6] Hamilton J C, Gemperline P J. Journal of Chemometrics, 1990, 4(1): 1. [本文引用:1]
[7] Borgen O S, Kowalski B R. Analytica Chimica Acta, 1985, 174(Supplement C): 1. [本文引用:1]
[8] Craig Hamilton J, Gemperline P. Mixture Analysis using Factor Analysis. Ⅱ: Self-Modeling Curve Resolution, 1990. [本文引用:1]
[9] Leger M N, Wentzell P D. Chemometrics & Intelligent Laboratory Systems, 2002, 62(2): 171. [本文引用:1]
[10] Miron S, Dossot M, Carteret C, et al. Chemometrics & Intelligent Laboratory Systems, 2011, 105(1): 7. [本文引用:1]
[11] Golshan A, Maeder M, Abdollahi H. Analytica Chimica Acta, 2013, 796(10): 20. [本文引用:1]
[12] Sawall M, Neymeyr K. Journal of Chemometrics, 2015, 27(5): 106. [本文引用:1]
[13] Tauler R. Journal of Chemometrics, 2001, 15(8): 627. [本文引用:1]
[14] Rajko R. Anal Chim Acta, 2009, 645(1-2): 18. [本文引用:1]
[15] Zhang X, Tauler R. Analytica Chimica Acta, 2013, 762(11): 25. [本文引用:1]