CN104408723A - 一种基于非负矩阵临近的拉曼光谱图像解混方法 - Google Patents

Abstract

本发明公开了一种基于非负矩阵临近的拉曼光谱图像解混方法，该方法针对拉曼光谱图像中的非纯像元不存在，拉曼光谱解混中的ANC、ASC和端元未知的限制，在拉曼光谱图像满足LMM模型情况下，通过对非负矩阵临近添加稀疏等多种正则化因子，从而可以有效地得到混合的像元的组分及其所对应的丰度。相对于现有其他方法，本发明具有以下优势：(1)无需端元(组成物质)的个数作为先验知识；(2)采取迭代策略，求取结果与端元数目无关；(3)加入了图关系和1/2范数，提高了算法的抗噪性和稳定性。

Description

一种基于非负矩阵临近的拉曼光谱图像解混方法

技术领域

本发明属于图像处理技术领域，具体涉及一种基于非负矩阵临近的拉曼光谱图像解混方法。

背景技术

拉曼光谱是一种三维图像，包括普通二维平面图像信息和波长信息。在对目标的空间特征成像的同时，对每个空间像元经过色散形成几十个乃至几百个窄波段以进行连续的光谱覆盖。一个拉曼图像为有若干个波长对应的二维图像组成的三维拉曼图像。

拉曼光谱因其对于不同物质都会产生唯一性的光谱特性，而被广泛应用细胞工程(细胞壁结构成分探测)等行业。然而目前的大多数拉曼图像都是由多个不同的物质(端元)混合合成，为了更加精准的每个混合成分进行分析，就需要对拉曼光谱图像进行解混分析，通常需要假设拉曼光谱图像满足线性混合模型(LMM)，该模型中的端元丰度需要满足非负(ANC)以及和为1的限制(ASC)。通常情况下，解混过程包括端元提取和丰都反演两个步骤。对于端元提取来讲，主要可以分为监督方法和非监督方法。监督方法假设所有的端元都是已知的，主要包括定点成分分析，自动端元提取，纯像元指数以及迭代误差分析，这些方法主要从几何视角进行分析，但是上述方法必须要求该几何体中需要至少存在一个端元。当算法中没有纯端元的情况下，最小体积转化及其相类似的方法(迭代限制端元)采取包含所有数据的最大的单纯形。这种方法的局限性在于必需存在N-1个端元(N为端元总数)，但在真实的高混合的数据集中，这种假设不理想。当所有的端元提取后，通常利用全限制的最小二乘预测或最大似然分析对端元进行丰度反演。当端元和其对应的丰富不确定，高光谱的解混问题就可以看成是盲信号分离问题，常见的方法包括包括独立主成分和非负矩阵分析。对于独立主成分方法来讲，其要求的端元之间相互独立在实际图像中不现实。

发明内容

针对现有技术所存在的上述技术问题，本发明提供了一种基于非负矩阵临近的拉曼光谱图像解混方法，通过对非负矩阵临近添加稀疏等多种正则化因子，从而可以有效地得到混合像元的组分及其所对应的丰度。

一种基于非负矩阵临近的拉曼光谱图像解混方法，包括如下步骤：

(1)对拉曼光谱图像进行降噪和去荧光处理；

(2)对于处理后图像中的任一像素，计算该像素与其他像素之间的权重，从而得到权重矩阵W；

(3)根据所述的权重矩阵W，计算出图关系的正则化因子P；

(4)对拉曼光谱图像的解混模型引入图关系正则化因子P和稀疏化因子，得到如下解混优化模型：

其中：f(u,v,Γ)为关于u、v和Г的拉格朗日函数，u为拉曼光谱图像的基本元素向量，v为拉曼光谱图像的残余向量，M为处理后的拉曼光谱图像，|| ||_F为F-范数，Tr()表示矩阵的迹，T表示矩阵转置，Г为拉格朗日算子矩阵，μ和λ均为拉格朗日算子，和φ均为正则化参数；Г_ij为拉格朗日算子矩阵的第i行第j列元素值，()_ij为()内矩阵的第i行第j列元素值，i和j均为自然数且1≤i≤n，1≤j≤n，n为拉曼光谱图像的像素个数；

(5)通过更新拉曼光谱图像M和拉格朗日算子λ，循环对上述解混优化模型最小化求解m次，得到拉曼光谱图像中各物质对应的基本元素向量u，进而合并得到拉曼光谱图像的丰度矩阵U；m为拉曼光谱图像中物质的种类个数。

所述的步骤(2)中，根据以下算式计算像素与像素之间的权重，从而得到权重矩阵W；

$W_{\mathrm{ij}}=e^{-\frac{{\left|\right|m_i-m_j\left|\right|}^2}{\mathrm{\σ}}}$

其中：m_i和m_j分别为处理后拉曼光谱图像中第i个像素和第j个像素的波长向量，σ为热核函数的核宽度，W_ij为处理后拉曼光谱图像中第i个像素与第j个像素之间的权重即权重矩阵W的第i行第j列元素值。

所述的步骤(3)中，根据以下算式计算图关系的正则化因子P：

P＝Tr(u^TLu)

其中：L＝D-W，D为对角矩阵且其第i行第i列元素值W_ij为权重矩阵W的第i行第j列元素值。

所述的步骤(5)中，通过以下迭代方程组对解混优化模型进行最小化求解：

$u_k=u_{k-1}.*(Z_{k-1}{v}_{k-1}+\mathrm{\μ}{W}_{u_{k-1}})./(u_{k-1}{v}_{k-1}^T{v}_{k-1}+\mathrm{\μD}{u}_{k-1}+\frac{\mathrm{\λ}}{4}{u}_{k-1}^{-1/2})$

$v_k=v_{k-1}.*Z_{k-1}^T{u}_{k-1}./u_{k-1}{V}_{k-1}^T{v}_{k-1}$

$\begin{array}{cc}{Z}_k=M-{\mathrm{\Γ}}_k& {\mathrm{\Γ}}_k=\max \{0,{\mathrm{\Γ}}_{k-1}+\frac{1}{k-1}(u_{k-1}{v}_{k-1}^T-M)\}\end{array}$

其中：u_k和u_k-1分别为第k次和第k-1次迭代的基本元素向量，v_k和v_k-1分别为第k次和第k-1次迭代的残余向量，Г_k和Г_k-1分别为第k次和第k-1次迭代的拉格朗日算子矩阵，.*和./分别表示点乘运算和点除运算，Z_k和Z_k-1分别为第k次和第k-1次迭代的中间矩阵，k为迭代次数。

所述的步骤(5)中，根据以下算式对拉曼光谱图像M和拉格朗日算子λ进行更新：

$\begin{array}{cc}{M}_t=M_{t-1}-u_{t-1}{v}_{t-1}^T& {\mathrm{\λ}}_t=\frac{1}{\sqrt{L}}\underset{l=1}{\overset{L}{\mathrm{\Σ}}}\frac{\sqrt{n}-{\left|\right|M_t^l\left|\right|}_1/{\left|\right|M_t^l\left|\right|}_2}{\sqrt{n-1}}\end{array}$

其中：M_t和M_t-1分别为第t次和第t-1次最小化求解过程中的拉曼光谱图像，u_t-1和v_t-1分别为第t-1次最小化求解得到的基本元素向量和残余向量，λ_t为第t次最小化求解过程中的拉格朗日算子，为拉曼光谱图像M_t中所有像素对应第l个波长上的像素值组成的n维向量，t和l均为自然数且1＜t≤m，1＜l≤L，L为拉曼光谱图像的波长维度。

本发明针对拉曼光谱图像中的非纯像元不存在，拉曼光谱解混中的ANC、ASC和端元未知的限制，在拉曼光谱图像满足LMM模型情况下，通过对非负矩阵临近(NMU)添加稀疏等多种正则化因子，从而可以有效地得到混合的像元的组分及其所对应的丰度。本发明假设拉曼光谱中的像元满足线性模型，提出了全新的非监督解混算法。为了满足现象模型中的非负值限定，本发明基于非负算法；为了使相似的像素点在分解中仍具备相似特性，本发明提出了图关系用来定义像素点的相似性；为了去除拉曼光谱的噪声和噪声，本发明使用了小波重构和Savitzky-Golay二阶过滤器；为了更好地表示数据的稀疏特性，本发明使用了L_1/2范数替代传统的L₀、L₁和L₂范数。

故相对于现有其他方法，本发明具有以下优势：1)无需端元(组成物质)的个数作为先验知识；2)采取迭代策略，求取结果与端元数目无关；3)加入了图关系和1/2范数，提高了算法的抗噪性和稳定性。

附图说明

图1为本发明拉曼光谱图像解混方法的流程示意图。

具体实施方式

为了更为具体地描述本发明，下面结合附图及具体实施方式对本发明的技术方案进行详细说明。

实施例1：

本发明基于秩1-NMU算法，提出了图关系稀疏化的非负矩阵临近算法。

如下式所示，基本的秩1-NMU算法可以表述为优化问题：

$f(u,v)=\underset{u\≥0,v\≥0}{\mathrm{\Σ}}\min {\left|\right|M-{\mathrm{uv}}^T\left|\right|}_F^2,\mathrm{such\; that}{\mathrm{uv}}^T\≤M$

其中：为F范数。矩阵为图像数据矩阵，矩阵u和矩阵v分别为长度为m和长度为n的行向量。

如图1所示，本发明基于非负矩阵临近的拉曼光谱图像解混方法，包括如下步骤：

(1)对拉曼图像进行降噪和去荧光处理，得到处理后的图像数据。

(2)对于图像数据中的像素点，根据如下式计算像素点之间的权重：

$W_{\mathrm{ij}}=e^{-\frac{{\left|\right|m_i-m_j\left|\right|}^2}{\mathrm{\σ}}}$

其中：m_i和m_j为图像数据中的任意两个像素点，σ为热核函数的核宽度。

(3)如下式所示，计算图关系正则化因子：

$\begin{array}{c}\frac{1}{2}\underset{i,j=1}{\overset{m}{\mathrm{\Σ}}}{\left|\right|u_i-u_j\left|\right|}^2{W}_{\mathrm{ij}}=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{u}_i{U}_i^T{D}_{\mathrm{ii}}-\underset{i,j=1}{\overset{m}{\mathrm{\Σ}}}{u}_i{u}_i^T{W}_{\mathrm{ij}}\\ =\mathrm{Tr}\left(u^T\mathrm{Du}\right)-\mathrm{Tr}\left(u^T\mathrm{Du}\right)\\ =\mathrm{Tr}\left(u^T\mathrm{Lu}\right)\end{array}$

其中，Tr(·)为矩阵的迹，D是对角矩阵D_ii＝Σ_jW_ij and L＝D-W。

(4)基于图关系和多种正则化，拉曼光谱解混问题可以转换为优化如下式所示的解混模型：

其中：μ，λ为拉格朗日算子。对于上式来说，采用单步更新策略对u和v进行估计，通过多次迭代，得到一次秩后的分解结果。接着，使用M←max(0,-(M-u_kv_k ^T))使用对M进行更新，继续重复一次秩后的分解结果。直到秩达到目标大小，停止迭代。

对上述解混模型进行最小化求解，具体包括如下步骤：

4.1拉格朗日算子的确定

在解混模型中，需要确定μ，Γ，λ三个拉格朗日算子。μ用于控制光滑度，其值可以在实验中确定。Γ可以根据Γ^(k+1)←max(0,Γ^(k)+a_k(uv^T-M))进行更新，其中k为当前的迭代次数。对于λ，可以根据下式进行更新。

$\mathrm{\λ}=\frac{1}{\sqrt{m}}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\frac{\sqrt{n}-{\left|\right|M_i\left|\right|}_1/{\left|\right|M_i\left|\right|}_2}{\sqrt{n-1}}$

4.2u和v的更新策略

在解混模型中，需要考虑在迭代过程中u和v的更新策略，在本实施方式中，我们通过KKT条件对模型进行求解，得到如下式所示的单步更新策略。

$u=u.*(\mathrm{Zv}+\mathrm{\μWu})./({\mathrm{uv}}^{T}v+\mathrm{\μDu}+\frac{\mathrm{\λ}}{4}{u}^{-1/2})$

v＝v.*Z^Tu./uv^T v

4.3解混算法后的后续处理步骤

在算法停止后，本实施方式将图像数据M分解为U和V。其中矩阵U为丰度矩阵，V由于误差扩散原因，无法被当作端元矩阵。为了得到丰度矩阵和端元矩阵。

对于丰度矩阵，由于必须满足LMM模型中的归一化限制，因此对于基本成分矩阵U，使用如下式对其归一化。为了避免零值的出现，我们在归一化过程中添加了极小值ω。

$\begin{array}{cc}{u}_j=\frac{u_j}{{\max }_i\left(u_j^i\right)}\&ForAll;j& u^i\&LeftArrow;\frac{u^i}{{\left|\right|u^i\left|\right|}_1+\mathrm{\&omega;}},\mathrm{\&omega;}<<1\end{array}$

其中：为像素点i中基本成分j的反射值。

已知图像数据矩阵M和丰度矩阵U，我们通过求解如下式的非负最小二乘问题得到端元矩阵V。

$\underset{V\&GreaterEqual;0}{\min {\left|\right|M-U{V}^T\left|\right|}_F^2}$

4.4算法的停止策略

停止策略是影响本实施方式执行性能的重要因素。在本实施方式中，我们预设最小残差作为算法停止的临界值，除此之外，我们设置最大迭代次数作为算法停止的临界值，两种临界值满足任意一个，本算法停止。

实施例2：

以下我们采用茶叶植物细胞壁的结构成分为例，对本发明提出的方法的具体实施方式进行叙述。

仪器准备：

实验设备由电子计算机、共焦拉曼光谱显微镜、尼康的浸入式纤维性，可见光谱波段的激光，光谱精度为～4cm^-1，激光电压为10mw。

材料准备：

在实验温室中采集新鲜的西湖龙井叶片1片，洗涤，放置在实验室台；对龙井页面制成切片，提取光谱范围579.335cm^-1～3062.07cm^-1。

拉曼图像预处理：

对拉曼光谱像，裁剪大小为的感兴趣区域(25×25)的细胞区域。为了保证实验的准确性，剔除出受透明器皿影响的波段(1751.2cm-1～3062.07cm-1，1450.7～cm-11471.4cm-1and 1517.8cm-1～1576.8cm-1)，总共选取462个光波对应的清晰拉曼光谱图像。

本实施例中待处理的拉曼图像为m个波长大小为的高光谱图像，图像大小为N＝25×25。本实施例基于非负矩阵临近的拉曼光谱解混方法，其实现过程如下所示：

1)噪音和自荧光性的去除

为了去除噪音和自荧光性，我们使用了小波变化中的‘db1’小波基将波长信号分解到6层，仅仅保留6层的高频信号，并使用该信号对原信号进行重构，有效的去除了基底噪音对图像的影响，接着，我们使用小波重构和Savitzky-Golay二阶过滤器，宽度为7个像素宽度对信号进行过滤，可以有效的去除了自荧光性对图像信号的干扰。

2)拉格朗日算子和初始化参数的选择

在本实例中，需要确定μ，Γ，λ三个拉格朗日算子。μ用于控制光滑度，在实例中，在μ∈(0,0.4)中以0.05为步长进行多次实验，最终得出当μ＝0.1,算法可以取到最优值。Γ可以根据Γ^(k+1)←max(0,Γ^(k)+a_k(uv^T-M))进行更新，其中k为当前的迭代次数。对于λ的更新与实施例1中相同。在初始化过程中，我们通过对数据矩阵M使用奇异值分解得到对于第一次迭代的u和v。

3)算法的停止策略

停止策略是影响本实例执行性能的重要因素。在本实施方式中，我们预设最小残差为0.001，除此之外，我们设置最大迭代次数作为算法停止的临界值，在本实验中，通过实验分析得知，最大迭代次数为100时，算法已经收敛。

利用本实施方式对采集到茶叶细胞壁结构成分进行确定，能够区分被细胞壁结构中各个成分的分布和各个成分的权重。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员，在不脱离本发明构思的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明保护范围内。

Claims (5)

1.一种基于非负矩阵临近的拉曼光谱图像解混方法，包括如下步骤：

(1)对拉曼光谱图像进行降噪和去荧光处理；

(2)对于处理后图像中的任一像素，计算该像素与其他像素之间的权重，从而得到权重矩阵W；

(3)根据所述的权重矩阵W，计算出图关系的正则化因子P；

(4)对拉曼光谱图像的解混模型引入图关系正则化因子P和稀疏化因子，得到如下解混优化模型：

其中：f(u,v,Γ)为关于u、v和Г的拉格朗日函数，u为拉曼光谱图像的基本元素向量，v为拉曼光谱图像的残余向量，M为处理后的拉曼光谱图像，||||_F为F-范数，Tr()表示矩阵的迹，T表示矩阵转置，Г为拉格朗日算子矩阵，μ和λ均为拉格朗日算子，和φ均为正则化参数；Г_ij为拉格朗日算子矩阵的第i行第j列元素值，()_ij为()内矩阵的第i行第j列元素值，i和j均为自然数且1≤i≤n，1≤j≤n，n为拉曼光谱图像的像素个数；

(5)通过更新拉曼光谱图像M和拉格朗日算子λ，循环对上述解混优化模型最小化求解m次，得到拉曼光谱图像中各物质对应的基本元素向量u，进而合并得到拉曼光谱图像的丰度矩阵U；m为拉曼光谱图像中物质的种类个数。

2.根据权利要求1所述的拉曼光谱图像解混方法，其特征在于：所述的步骤(2)中，根据以下算式计算像素与像素之间的权重，从而得到权重矩阵W；

$W_{\mathrm{ij}}=e^{-\frac{{\left|\right|m_i-m_j\left|\right|}^2}{\mathrm{\&sigma;}}}$

其中：m_i和m_j分别为处理后拉曼光谱图像中第i个像素和第j个像素的波长向量，σ为热核函数的核宽度，W_ij为处理后拉曼光谱图像中第i个像素与第j个像素之间的权重即权重矩阵W的第i行第j列元素值。

3.根据权利要求1所述的拉曼光谱图像解混方法，其特征在于：所述的步骤(3)中，根据以下算式计算图关系的正则化因子P：

P＝Tr(u^TLu)

其中：L＝D-W，D为对角矩阵且其第i行第i列元素值为权重矩阵W的第i行第j列元素值。

4.根据权利要求1所述的拉曼光谱图像解混方法，其特征在于：所述的步骤(5)中，通过以下迭代方程组对解混优化模型进行最小化求解：

$u_k=u_{k-1}.*(Z_{k-1}{v}_{k-1}+\mathrm{\&mu;W}{u}_{k-1})./(u_{k-1}{v}_{k-1}^T{v}_{k-1}+\mathrm{\&mu;D}{u}_{k-1}+\frac{\mathrm{\&lambda;}}{4}{u}_{k-1}^{-1/2})$

$v_k=v_{k-1}.*Z_{k-1}^T{u}_{k-1}./u_{k-1}{v}_{k-1}^T{v}_{k-1}$

$\begin{array}{cc}{Z}_k=M-{\mathrm{\&Gamma;}}_k& {\mathrm{\&Gamma;}}_k=\max \{0,{\mathrm{\&Gamma;}}_{k-1}+\frac{1}{k-1}(u_{k-1}{v}_{k-1}^T-M)\}\end{array}$

其中：u_k和u_k-1分别为第k次和第k-1次迭代的基本元素向量，v_k和v_k-1分别为第k次和第k-1次迭代的残余向量，Г_k和Г_k-1分别为第k次和第k-1次迭代的拉格朗日算子矩阵，.*和./分别表示点乘运算和点除运算，Z_k和Z_k-1分别为第k次和第k-1次迭代的中间矩阵，k为迭代次数；D为对角矩阵且其第i行第i列元素值W_ij为权重矩阵W的第i行第j列元素值。

5.根据权利要求1所述的拉曼光谱图像解混方法，其特征在于：所述的步骤(5)中，根据以下算式对拉曼光谱图像M和拉格朗日算子λ进行更新：

$\begin{array}{cc}{M}_t=M_{t-1}-u_{t-1}{v}_{t-1}^T& {\mathrm{\&lambda;}}_t=\frac{1}{\sqrt{L}}\underset{l=1}{\overset{L}{\mathrm{\&Sigma;}}}\frac{\sqrt{n}-{\left|\right|M_t^l\left|\right|}_1/{\left|\right|M_t^l\left|\right|}_2}{\sqrt{n-1}}\end{array}$

其中：M_t和M_t-1分别为第t次和第t-1次最小化求解过程中的拉曼光谱图像，u_t-1和v_t-1分别为第t-1次最小化求解得到的基本元素向量和残余向量，λ_t为第t次最小化求解过程中的拉格朗日算子，为拉曼光谱图像M_t中所有像素对应第l个波长上的像素值组成的n维向量，t和l均为自然数且1＜t≤m，1＜l≤L，L为拉曼光谱图像的波长维度。