【作者】种敏琪，秦斌杰

上海交通大学生命科学技术学院生物医学工程系，上海，200240

光谱分离(Spectal Unmixing)技术在卫星遥感成像、生物荧光成像领域是一个重要的技术。国内外在卫星遥感成像光谱分离的研究已经有几十年的历史，但还没有很成熟的算法解决，这个领域公认的难点譬如如何有效地提取端元[1-3]。在遥感成像的高光谱场景中，各个像元的信息大都是由不同组分混合而成的，所谓端元提取就是提取出不同组分纯指纹光谱的一个处理过程。基于这些不同组分物质的纯指纹光谱，通过谱分离算法再把遥感场景分解成不同丰度(abundance fraction)的多种物质。

近年来，随着生命科学技术的发展，生物荧光成像谱分离领域[4]也有极大的发展。生物荧光成像大都利用检测目标荧光蛋白的方法，观察活体小动物体内单细胞水平的生物学演变过程。由于动物组织以及胃内食物中存在弹性蛋白和胶原蛋白等类物质，这些蛋白质的激发波长与发射波长与GFP等目标荧光蛋白类似。因此，在这些谱段进行荧光成像，荧光图像均存在微弱的自发荧光。另外，如果动物体内有多标记目标荧光组分时，则在一定的激发光照射下，在高灵敏CCD采集的多通道图像像元上会有多个荧光组分的混合，基于这些现象，谱分离技术就成了荧光成像的必然要求。

卫星遥感成像时，往往在很多个成像通道进行连续谱段高光谱图像采集，这可以确保有足够的冗余信息准确分辨出每一种组分在高光谱图像中的丰度。而在生物荧光成像时，往往感兴趣的在体目标荧光组分数目一般在1～2个左右，再加上引起小动物自发荧光的组分数目也可以设定为1～2个左右，因此，原则上用于生物荧光成像谱分离的多光谱图像，只需在3～5个谱段采集就行。因此，本研究一个重要的目的就是验证，在遥感成像研究领域常用的端元提取和线性谱分离算法是否适用于生物荧光成像的谱分离。

1 原理及方法

荧光成像谱分离的主要原理是：注射了荧光基团的小动物在多个波长激发光(如图1用波长为488和波长为594的激发光分别激发照射荧光染料Alexa Fluor 488 Dye 和Alexa Fluor 594 Dye[5])的激发照射下，对应不同谱段的发射光图像会在多个成像通道被采集。基于每个成像像元是不同组分荧光基团的线性叠加假设，线性谱分离技术从混合的像元中分离出不同光谱特性荧光组分的贡献比例（所谓丰度contribution fractions），进而根据这些荧光组分的丰度推测目标荧光蛋白在小动物体内目标区域的表达程度。

线性谱分离算法的主要步骤[6]是：首先通过如主分量变化、最小噪声分离[7](MNF, minimum noise fraction，MNF变换等同于噪声调节主分量变换，Noise Adjusted Principle Component transform[8])等变换技术，对多/高光谱图像进行有效降维，使计算量降低；对构成混合像元各特定组分的参考光谱进行估计，这里也称为端元估计(endmember determination)，包括端元数量的估计以及各端元参考光谱的估计；接着就是求逆的过程，即依据多谱段输出图像和各端元参考光谱，估计出各特定组分的丰度。具体来说，线性谱分离满足线性叠加原理：

这里b为像素的亮度，a为端元本身的强度，s为端元丰度，w为噪声，m为端元总数量，满足上式对于图像中的每一个像素均适用，将一个通道采集图像的每个像素信息添加为矩阵的一列。对于多个通道的采集图像来说，上面的公式可以标识为矩阵形式如下：

其中B为采集的多通道图像的像素亮度矩阵，S为待求取的端元丰度矩阵，A为端元参考光谱矩阵（即各组分在各个荧光波段下采集到的信号强度，该数据虽可以参考出厂标注的荧光染料发射光谱，如图1a，b实线显示的即为Alexa Fluor 488 Dye 和Alexa Fluor 594 Dye的发射光谱，但一般需要从采集图像中估计出来才更有现实意义），W为噪声矩阵。通过对多通道图像进行线性方程组的求解，即可获得每个荧光基团端元在不同像素中的丰度，从而分离出不同荧光基团在小动物体内空间位置及表达强度等的更进一步信息。

图1 Alexa Fluor 488 Dye 和Alexa Fluor 594 的吸收光谱和发射光谱[5]Fig.1 The absorption and emission spectra for Alexa Fluor 488 Dye and 594 Dye[5]

在图2中，我们形象地简述一两个波段下谱分离成像的背景。实心点代表着纯像素的端元点，它们在波段i和波段j上是光亮度的最大值或者最小值。而空心圆点代表着普通的混合像素点，它们是经由三个端元进行线性叠加而得到的。由于公式1的约束条件，这些混合点一定落在三个端元所围绕的凸包之内。

图2 两个波段下谱分离成像的原理简图Fig.2 A schematic diagram of spectral unmixing for two-band multispectral imaging

首先，我们需要对图像进行降噪和降维的处理，MNF算法在降维的同时还可以对噪声进行初步的抑制，效果较好。谱分离成像的过程分为两步：第一步是端元的确定，通过所有的点在不同波段下的亮度，求出端元点数，这也是本篇文章研究的重点；第二步通过公式(2)，在已知端元的情况下，我们可以利用最小二乘法[12,13]解出矩阵S。从而获知每个混合像素点中各个端元的组成比例，达到谱分离成像的目的。

现在流行的端元分离算法主要有：convex cone analysis、iterative error analysis(IEA)、 pixel purity index(PPI)、N-FINDR以及orthogonal subspace projection等几种[1-3]。由于我们需要在仅输入每个像素点亮度的情况下分离出端元点，并且希望算法是能在非人为控制下自动完成的。因此，通过进行分析比较，我们发现适合的经典端元分离方法主要有2种，基于像素点特征的PPI[9,10]以及基于矩阵特性的N-FINDR[9]算法。现在对这两种算法的原理及流程进行分别阐述。

1.1 PPI算法[9-10]

PPI算法的原理是在每一个频谱段，总有一个端元会产生最强的亮度，而其它的像素点的亮度必定小于等于最强端元产生的亮度。基于这个原理，我们排除不可能成为端元的点，获取可能的端元数目及其分布。

为了达到这个目的，我们采取向量内积的方法，通过随机生成大量的相互独立的矢量，计算每个像素点在这些矢量上的内积，筛选其中内积最大或者最小的，作为可能的备选端元。然后不断重复这个过程，直到所有的剩余像素点都可以在某一方向向量上有着最大或者最小的内积为止，这些剩余像素点就是端元。如下图3中，画的点在参考矢量上拥有着最大或者最小的投影，因此在本轮筛选中将其作为可能的端元而保留。

图3 PPI算法原理Fig.3 The basic principle of PPI algorithm

我们可以使用Auto-PPI(APPI)算法全自动地完成整个筛选过程，其算法简单描述如下：

1）随机生成多个线性无关向量；

2）对图像中的每一个点，将其多光谱向量与随机生成向量做内积运算；

3）统计内积运算结果，剔除在任一随机向量上内积均不为最大值的像素点；

4）如果剔除了新的点，则跳回1；

5）剩余像素点为所需端元。

PPI流程如图4所示。

1.2 N-FINDR算法[11]

N-FINDR算法是一种在已知端元数量的前提下比较好的算法。它利用了各个波段下端元点之间的亮度差值大于混合点之间亮度差的原理，通过尝试不同点相互之间的组合进行比较，从而确定更加接近端元的点集。而在图像中存在端元的条件下，这些结果点集就是端元点的集合。

根据以上的介绍，算法的实现步骤主要是通过随机选取的像素作为现有可能像素，然后对备选像素的矩阵做行列式，如果在替换了新的像素后行列式的值增加，那么就用新的像素替换原有像素，直到所有的像素点均尝试过为止。此算法描述如下：

图4 PPI算法的流程图Fig.4 The flow chart of PPI algorithm

1）随机选取点作为初始行列式并计算其值；

2）从剩余的点中选取一个替代原有点并计算新的行列式的值；

3）若行列式值变大，则替换这个点，否则，舍弃这个点；

4）跳回1直到所有的点被尝试过；

5）最终行列式取最大值的像素点为所需端元。N-FINDR流程如图5所示。

图5 N-FINDR算法流程图Fig.5 The flow chart of N-FINDR algorithm

1.3 两种算法的混合

PPI算法对于随机向量有着一定的要求，否则将会产生一些并不真实的端元，尤其对于并不满足严格线性叠加条件的实验数据，某些混合点的亮度值过大，PPI算法会将它们视为端元并计入到计算结果中。N-FINDR算法则涉及了复杂的行列式计算，并且需要预先明确端元数量。而在生物荧光成像这种波段数较少的条件下，它们的缺点则更为明显。为此，我们对两种算法进行了合并，可以得到一个更好的实验结果。合并中，程序先根据适当的阈值进行多次PPI算法的循环，再对剩余点进行N-FINDR算法，来获得一个更加可靠的结果。

2 实验结果及讨论

2.1 C++随机模拟数据

我们利用C++中algorithm库的随机函数生成了由四个虚拟端元组成的1000个像素点矩阵，生成了一个类似于图2的四维矩阵(图2为2维)。在这个矩阵中，包含了4个端元点以及他们在公式(1)约束下生成的大量混合像素点（类似图2），通过对这些随机数据求解，来验证算法的有效性。

这4个虚拟端元为：(100,0,10,0)、(1,100,1,1)、(0,0,100,50)、(10,10,10,100)括号内数字的含义表明端元点在不同波段的亮度值（公式2中的矩阵A），每个混合像素点的亮度正是由这些端元亮度值的线性叠加（通过随机生成的丰度矩阵S）而形成的。通过公式

进行归一化，其中b'为归一化后的值，为像素点本身值（既随机生成像素点的亮度），max为最亮像素点值（最大值为100），min为最暗像素点值（最小值为0）。之后通过PPI，N-FINDR和混合算法3种方法分别计算出端元。计算的结果表明，对于此类完全符合线性叠加原理，并且每个端元之间有着明显的差异的理想情况下，PPI和N-FINDR算法都有着良好的效果。能够准确地分离出4个端元。

2.2 实验成像及处理

我们使用的成像系统是QuickView3000，荧光染料为Alexa Fluor 488 Dye 和Alexa Fluor 594 Dye，其参考光谱见图1。我们分别把两种染料以5种不同的混合比例（100%，0），（75%，25%），（50%，50%），（25%，75%），（0%，100%）置于在96孔板中。采集两幅图像的带通波长范围分别为470～530 nm和565～630 nm，成像结果见图6所示。

图6 两通道多光谱图像采集Fig.6 The two-band images in multispectral imaging

对于如上采集得到的两幅多图像，我们分别通过MATLAB从成像系统进行导入和归一化处理。再通过C++实现PPI和N-FINDR算法，分别对其进行端元分离。结果见表1。

表1 PPI算法和N-FINDR算法的实现结果Tab.1 The results of PPI and N-FINDR algorithms

图7 端元像素点的分布Fig.7 The distribution of endmember pixels

在图7中可以看到，端元1端元3所在像素落在红色圆圈的试管内，它们是真实的端元。端元2端元4所在像素落在蓝色的试管内，它们是假的端元。

我们发现，因为PPI算法由于是自动确定端元，所以计算出来的端元数量是可能变化的。从表1中可以看到，端元2和端元4并非真正的端元，但是它们依旧出现在了PPI的计算结果中，其具体的出现与否取决于PPI算法随机化初始条件。在某些情况下，这些端元并不会被算法本身筛选掉。我们虽然可以通过多次随机化以及增加筛选次数（试验中发现，将筛选轮数由3增加到30可以有效的减少端元2和端元4在结果中出现的频率），但是这样的代价妨碍了算法实现的效率。

而对于N-FINDR算法，虽然其准确率较高，但是非常依赖于事前对于端元数的估算。我们可以看到，对于不同的假定端元数，直接决定了端元2是否存在于最后的计算结果之中。

我们可以看出两种算法在端元分离中的各自优缺点：N-FINDR有着计算速度快，能够准确地计算出最优m个端元的优点，但同时需要预先对m有一个准确的确定；PPI算法则计算速度较慢，并且会挑出一些并不是端元的可能像素点，但是同时却也有着不会因为m的误差而筛除端元的可能。

因此，将两种算法的优点相结合，先通过一个较少循环的PPI算法剔除大量的像素点，然后再通过N-FINDR算法对于剩余的像素点进行计算处理（由于像素点减少，此时的计算复杂度大幅下降），可以更快速准确地分离出我们所需要的端元点。

3 结论与展望

在上述实验中，我们分别分析了PPI算法和N-FINDR算法在生物荧光成像实验中的效果，对比了两个算法各自的优缺点。通过生物荧光染料的实验验证，发现将两个算法的优点合并使用，可以得到一个更优的处理结果，这一点为实验结果所证实。考虑到更复杂的试验，由于不同深度的荧光源，在生物体内传播时会发生散射、折射和反射等现象，导致像素亮度的叠加也可能并不是绝对线性叠加的。另外，图像中也可能仅存在混合像元点，而不存在纯的端元像素点，而上述两种算法的实现都要求图像中必须包含纯的端元点才能将其分离。因此，我们计划在后续的试验中尝试并改进其它的谱分离算法来适应生物荧光成像的特点，例如盲分离[14]算法可以在图像中不含有纯端元点以及无需预先估计端元参考光谱的情况下，直接进行谱分离来得出各个端元的丰度。

[1]Maciel Zortea and Antonio Plaza, Spatial Preprocessing for Endmember Extraction, IEEE Transactions on Geoscience and Remote Sensing, 2009, 47(8): 2679-2693.

[2]李素, 李文正, 周建军等, 遥感影像混合像元分解中的端元选择方法综述[J]2007, 23(5): 35-42.

[3]Antonio Plaza and Chein-I Chang.Impact of Initialization on Design of Endmember Extraction Algorithms[J].IEEE Transactions on Geoscience and Remote Sensing, 2006, 44(11):3397-3407.

[4]Heng Xu, and Brad W.Rice.In-vivo fluorescence imaging with a multivariate curve resolution spectral unmixing technique [J].Journal of Biomedical Optics, 2009, 14(6): 064011-1-9.

[5]http://www.invitrogen.com/site/us/en/home/support/Research-Tools/Fluorescence-SpectraViewer.html.

[6]Nirmal Keshava and John F.Mustard, Spectral Unmixing[J],IEEE Signal Processing Magazine, 2002, 44-57.

[7]J.W.Boardman, and F.A.Kruse.Automated spectral analysis:a geological example using AVIRIS data, north Grapevine Mountains, Nevada[A].Proceedings, ERIM Tenth Thematic Conference on Geologic Remote Sensing[C]Environmental Research Institute of Michigan, 1994, 1:407-418

[8]R.E.Roger.A fast way to compute the noise-adjusted principal components transform matrix[J]IEEE Transactions on Geoscience and Remote Sensing, 1994, 32(6):1194-1196.

[9]Mingkai Hsueh and Chein-I Chang.Field programmable gate arrays (FPGA) for pixel purity index using blocks of skewers for endmember extraction in hyperspectral imagery[J].International Journal of High Performance Computing Applications, 2008, 22(4): 408-423.

[10]Farzeen Chaudhry, Chao-Cheng Wu, Weimin Liu.et al.Pixel purity index-based algorithms for endmember extraction from hyperspectral imagery Recent Advances in Hyperspectral Signal and Image Processing [M], 2006

[11]Chao-Cheng Wu, Shihyu Chu and Chein-I Chang.Sequential N-FINDR Algorithms.Proc.of SPIE, 2008, 7086: 70860C.1-70860.12.

[12]Y.H.Hu, H.B.Lee, and F.L.Scarpace, Optimal Linear Spectral Unmixing[J].IEEE Transactions on Geoscience and Remote Sensing, 1999, 37(1):639-644.

[13]D.C.Heinz, and Chein-I Chang.Fully constrained Least sqaures linear spectral mixture analysis method for material quantification in hyperspectral imagry[J].IEEE Transactions on Geoscience and Remote Sensing, 2001, 39(3): 529-545.

[14]N.Dobigeon, S.Moussaoui, M.Coulon, et al Joint bayesian endmember extraction and linear unmixing for hyperspectral imagery[J], IEEE Transactions on Signal Processing, 57(11):2009, 4355-4368.