一种基于无限高斯混合模型的高光图图像解混方法

摘要

本发明公开了一种基于无限高斯混合模型的高光图图像解混方法，假设高光谱图像中的像元满足无限混合模型，这种假定与传统的线性模型相比，尤其在高分辨率的高光谱图像应用中，无限混合模型更能反映出图像像元的复杂性，为了降低计算复杂性，使用合理的降维策略；为了确定高斯组分的个数，利用虚拟维度估计组分个数，进而扩展为高斯组分个数的范围；为了求解无限混合模型，与传统的求解方法不同，本发明采用TTS策略有效的确定了高斯组分的个数，使用Metropolis-within-Gibbs方法确定无限混合模型中的参数和参超数，通过参数和超参数的采样，可以有效地得到混合的像元的组分机器所对应的丰度。

说明书

技术领域

本发明属于图像处理技术领域，涉及基于高光谱图像解混方法，尤其涉及 一种基于无限高斯混合模型的高光图图像解混方法。

背景技术

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

近红外外高光谱因其快速无损等特性被广泛应用于食品、医药、石油化工 等行业。然而由于目前的大多数高光谱图像都是由多个不同的物质(端元)混合合 成，为了更加精准的每个混合成分进行分析，就需要对高光谱图像进行解混分 析，通常需要假设高光谱图像满足线性混合模型(LMM)，该模型中的端元丰度需 要满足非负(ANC)和和为1的限制(ASC)。通常情况下，解混过程包括端元 提取和丰都反演两个步骤。对于端元提取而言，又可以分为端元数据确定和端 元提取两部分。就端元数目确定，第一类方法是基于像素的相关系数矩阵和协 方差矩阵，常见的有主成分分析(PCA)、Harsanyi-Farrand-Chang(HFC)、Akaike 信息准则的等方法，这些方法在低维度图像中工作正常，但是对于高维度的图 像的处理效果不好；第二类方法就是通过子空间的最小化来确定端元。对于端 元提取来讲，主要可以分为监督方法和非监督方法。监督方法假设所有的端元 都是已知的，主要包括定点成分分析(VCA),自动端元提取(AEE),纯像元指 数(PPI)，N-FINDR及迭代误差分析(IEA)，这些方法主要从几何视角进行分 析，但是上述方法必须要求该几何体中需要至少存在一个端元。当算法中没有 纯端元的情况下，最小体积转化(MVT)及其相类似的方法(迭代限制端元(ICES)) 采取包含所有数据的最大的单纯形。这种方法的局限性在于必需存在N-1个端 元(N为端元总数)，但在真实的高混合的数据集中，这种假设不理想。当所有 的端元提取后，通常利用全限制的最小二乘预测(FCLS)或最大似然分析对端 元进行丰度反演。当端元和其对应的丰富不确定，高光谱的解混问题就可以看 成是盲信号分离问题，常见的方法包括包括独立主成分(ICA)和非负矩阵分析 (NMF)。对于ICA来讲，其要求的端元之间相互独立在实际图像中不现实。对 于NMF来讲，求解NMF中的矩阵问题容易陷入最小解问题。

发明内容

针对近红外高光谱图像中的非纯像元不存在，高光谱解混中的ANC、ASC 和端元未知的限制，在高光谱数据模型满足经典的高斯混合模型情况下，本发 明提出一种基于无限高斯混合模型的高光图图像解混方法，通过使用分层贝叶 斯模型对高斯模型中的参数和非参数进行估计，从而可以有效地得到混合的像 元的组分机器所对应的丰度。

为了解决上述技术问题，本发明的技术方案如下：

一种基于无限高斯混合模型的高光图图像解混方法，

11)对高光谱图像进行降维处理，得到处理后的降维数据；

12)利用虚拟维度的方法确定高斯组分个数的大小，并得出高斯组分个数的范 围，针对每个高斯组分个数，利用K-means方法，进行分别聚类，对于每个聚 类的群组，使用PPI方法，提取每个群组中的最纯的像元作为高斯混合模型中 的期望向量；

13)对于高光谱图像中的每个像元，基于无限混合模型，采用两状态策略进行 端元数目采样，然后使用Metropolis-within-Gibb对无限混合模型中的参数和超参 数进行估计，通过多次迭代，得到最终的稳定的参数和超参数的估计；

所述高光谱中的像元满足无限高斯混合模型；

高光谱图像满足如式(a)所属的高斯模型：

$y=\underset{r=1}{\overset{R}{\Sigma }}{E}_r{\alpha }_r---\left(a\right);$

其中E_r为独立的高斯向量，y是高光谱图像中的某个像元，R为组成该像元的组 成个数，α_r为组成部分的比例即丰度，其需要满足如式(b)两种限制：

$\left(\begin{array}{cc}{\alpha }_r\ge 0,{\forall }_r=1,...,R& \left(\mathrm{ANC}\right)\\ \underset{r=1}{\overset{R}{\Sigma }}{\alpha }_r=1& \left(\mathrm{ASC}\right)\end{array}\right)---\left(b\right);$

在无限的高斯混合模型中，设定所有的高斯成分都相同，对于每个高斯成分而 言：

E_r|m_r,σ²～N(m_r,σ²I_L)  (c)；

其中m_r＝[m_r,1,...,m_r,L]^T是第r个高斯分布的均值向量，当所有的端元分布中的方 差为单位矩阵σ²I_L，因此，像元的似然函数可以表述为如式(d)所示：

$f\left(y\right|\theta )\propto \frac{1}{{\left[{\sigma }^{2}g\left(\alpha \right)\right]}^{L/2}}\exp (-\frac{{\left|\right|y-k\left(\alpha \right)\left|\right|}^2}{2{\sigma }^{2}g\left(\alpha \right)})---\left(d\right);$

其中θ＝{α,σ²,R,M_R}，||·||是标准的二阶范数，α＝[α₁,...,α_R]，M_R＝[m₁,...,m_R] 是由聚类算法产生的均值向量；

$g\left(\alpha \right)=\underset{r=1}{\overset{R}{\Sigma }}{m}_r{\alpha }_r,k\left(\alpha \right)=\underset{r=1}{\overset{R}{\Sigma }}{\alpha }_r^2---\left(e\right).$

进一步的，所述的步骤(11)中,使用谱聚类方法中使用的图理论，通过分 解高光谱数据的相似矩阵，计算相似矩阵的特征向量，排序后得到所需的降低 维度的数据集合。

进一步的，所述的步骤(12)中，包括如下步骤：

31)高斯组分个数的范围确定

通过虚拟维度估算出可能的高斯成分个数Rsim，为了尽可能的考虑到所有 的取值范围，如式(f)所示，基于估算的高斯成分个数Rsim,计算得到高斯成分 个数的取值范围Rmin和Rmax.

R_max＝floor(min(2R_sim,N))；

R_min＝ceil(max(R_sim/2,1))；  (f)

32)高斯组分均值向集合的确定

对于Rmin到Rmax中各值，都存在与之相对应的高斯组分均值向量集合， 对于R∈[R_min,...,R_max]，利用K-means将观测数据Y方法形成R个群组，针对 聚类后的每个群组，提取最纯的像元，由此构成数目为R的均值向量MR，因此， 对于所有的R值，可以得到高斯组成均值向量集合

进一步的，所述的步骤(13)中，包括如下步骤：

41)对于每个像元来讲，初始化混合模型中参数和超参数，其中包括高斯 组成个数R(1)，高斯成分的均值向量丰度向量α⁽¹⁾，方差σ²⁽¹⁾；

42)对于每次迭代而言，根据两状态策略对高斯组分的个数及其所对应的 丰度进行调整；

43)调整端元成分中涉及的“前进”状态概率和“后退”状态概率进 行调整；

44)利用Metropolis-within-Gibbs和后验密度对α^(t)、σ^2(t)和ω^(t)进行采样，再 执行步骤42)，直到迭代完毕。

本发明的有益效果在于：本发明假设高光谱图像中的像元满足无限混合模 型，这种假定与传统的线性模型相比，尤其在高分辨率的高光谱图像应用中， 无限混合模型更能反映出图像像元的复杂性。为了降低计算复杂性，使用合理 的降维策略；为了确定高斯组分的个数，利用虚拟维度估计组分个数，进而扩 展为高斯组分个数的范围；为了求解无限混合模型，与传统的求解方法不同，本 文采用TTS策略有效的确定了高斯组分的个数，使用Metropolis-within-Gibbs方 法确定无限混合模型中的参数和参超数，通过参数和超参数的采样，可以有效 地得到混合的像元的组分机器所对应的丰度。

附图说明

图1为无限混合模型中参数和超参数关系图。

具体实施方式

下面将结合附图和具体实施例对本发明做进一步的说明。

一种基于无限高斯混合模型的高光图图像解混方法，

11)对高光谱图像进行降维处理，得到处理后的降维数据；

12)利用虚拟维度的方法确定高斯组分个数的大小，并得出高斯组分个数的范 围，针对每个高斯组分个数，利用K-means方法，进行分别聚类，对于每个聚 类的群组，使用PPI方法，提取每个群组中的最纯的像元作为高斯混合模型中 的期望向量；

13)对于高光谱图像中的每个像元，基于无限混合模型，采用两状态策略进行 端元数目采样，然后使用Metropolis-within-Gibb对无限混合模型中的参数和超参 数进行估计，通过多次迭代，得到最终的稳定的参数和超参数的估计；

所述高光谱中的像元满足无限高斯混合模型；

高光谱图像满足如式(a)所属的高斯模型：

$y=\underset{r=1}{\overset{R}{\Sigma }}{E}_r{\alpha }_r---\left(a\right);$

其中E_r为独立的高斯向量，y是高光谱图像中的某个像元，R为组成该像元的组 成个数，α_r为组成部分的比例即丰度，其需要满足如式(b)两种限制：

$\left(\begin{array}{cc}{\alpha }_r\ge 0,{\forall }_r=1,...,R& \left(\mathrm{ANC}\right)\\ \underset{r=1}{\overset{R}{\Sigma }}{\alpha }_r=1& \left(\mathrm{ASC}\right)\end{array}\right)---\left(b\right);$

在无限的高斯混合模型中，设定所有的高斯成分都相同，对于每个高斯成分而 言：

E_r|m_r,σ²～N(m_r,σ²I_L)  (c)；

其中m_r＝[m_r,1,...,m_r,L]^T是第r个高斯分布的均值向量，当所有的端元分布中的方 差为单位矩阵σ²I_L，因此，像元的似然函数可以表述为如式(d)所示：

$f\left(y\right|\theta )\propto \frac{1}{{\left[{\sigma }^{2}g\left(\alpha \right)\right]}^{L/2}}\exp (-\frac{{\left|\right|y-k\left(\alpha \right)\left|\right|}^2}{2{\sigma }^{2}g\left(\alpha \right)})---\left(d\right);$

其中θ＝{α,σ²,R,M_R}，||·||是标准的二阶范数，α＝[α₁,...,α_R]，M_R＝[m₁,...,m_R] 是由聚类算法产生的均值向量；

$g\left(\alpha \right)=\underset{r=1}{\overset{R}{\Sigma }}{m}_r{\alpha }_r,k\left(\alpha \right)=\underset{r=1}{\overset{R}{\Sigma }}{\alpha }_r^2---\left(e\right).$

所述的步骤(11)中,使用谱聚类方法中使用的图理论，通过分解高光谱数据的 相似矩阵，计算相似矩阵的特征向量，排序后得到所需的降低维度的数据集合。 所述的步骤(12)中，包括如下步骤：

31)高斯组分个数的范围确定

通过虚拟维度估算出可能的高斯成分个数Rsim，为了尽可能的考虑到所有 的取值范围，如式(f)所示，基于估算的高斯成分个数Rsim,计算得到高斯成分 个数的取值范围Rmin和Rmax.

R_max＝floor(min(2R_sim,N))；

R_min＝ceil(max(R_sim/2,1))；  (f)

32)高斯组分均值向集合的确定

对于Rmin到Rmax中各值，都存在与之相对应的高斯组分均值向量集合， 对于R∈[R_min,...,R_max]，利用K-means将观测数据Y方法形成R个群组，针对 聚类后的每个群组，提取最纯的像元，由此构成数目为R的均值向量MR，因此， 对于所有的R值，可以得到高斯组成均值向量集合

所述的步骤(13)中，包括如下步骤：

41)对于每个像元来讲，初始化混合模型中参数和超参数，其中包括高斯 组成个数R(1)，高斯成分的均值向量丰度向量α⁽¹⁾，方差σ²⁽¹⁾；

42)对于每次迭代而言，根据两状态策略对高斯组分的个数及其所对应的 丰度进行调整；

43)调整端元成分中涉及的“前进”状态概率和“后退”状态概率进 行调整；

44)利用Metropolis-within-Gibbs和后验密度对α^(t)、σ^2(t)和ω^(t)进行采样，再 执行步骤42)，直到迭代完毕。

实施一：

本实例采用小青菜的农药残留检测为例，对本发明提出的方法的具体实施 方式进行叙述。

仪器准备

实验设备由电子计算机、高光谱仪、卤素灯、矫正黑、白板。高光谱仪使 用美国ASD(Analytical Spectral Device)公司的Handheld Field Spec光谱仪，光谱 采样间隔为1.5nm,采样范围为380nm～1030nm，采用漫反射方式进行样本光谱采 样；采用与光谱仪器配套的14.5卤素灯，进行光谱采集前须使用矫正黑、白板 对高光谱仪进行常规矫正。

材料准备

在实验温室中采集新鲜的小白菜叶片两片(标记为叶A和叶B)，洗涤，烘 干，放置在实验室台；配置两种不同浓度(农A，农B)的农药(嘧霉胺，杭州， 中国)，其中农A配置为1.8253g农药和10ml纯净水配比，农B配置为1.7431g 农药和20ml纯净水配比。首先对四种物品进行高光谱图像采集，接着，使用农 A对叶A进行涂点，使用农B对叶B进行涂点，形成两个有效地的样本(样A 和样B)。

高光谱图像预处理

对所有的高光谱图像，采用高光谱图像处理软件ENVI5.0对大小为的感兴 趣区域的区域高光谱图像进行选取。为了保证实验的准确性，剔除出受仪器影 响和光照影响的前150个波长所对应的光谱图像，只选取151-512总共342个光 波对应的清晰高光谱图像。

本实施例中待处理的高光谱图像为L个波长大小为的高光谱图像，图像大 小为N＝200×200。本实施例基于无限高斯混合模型的高光谱解混方法，其实现 过程如下所示：

1)降维过程

降维方法针对输入的高光谱像元数据Y＝[y₁,...,y_N],设置降维后的 维度为K＝3；首先计算距离的相似度矩阵S和S对应的归一化的拉普拉斯矩阵L。 通过求解矩阵L的前k个特征向量，并由此前k个特征向量组成了降维后的数 据矩阵X＝[x₁,...,x_N]。

2)高斯组分个数范围确定

首先使用虚拟维度确定高斯组分的个数的预估值Rsim，对于降维后的矩阵 X＝[x₁,...,x_N]，首先计算矩阵X的相关系数矩阵G和协方差矩阵V，并使用 和λ＝{λ₁,...,λ_N}分别表示G和V的特征向量。如式(7)所示，对 于像元X中到每个像元l,设定二值状态的假设

$H_0:{\mu }_l={\bar{\lambda }}_l-{\lambda }_l=0$

$H_1:{\mu }_l={\bar{\lambda }}_l-{\lambda }_l>0---\left(1\right)$

如果H1为真，则Rsim＝Rsim+1。根据公式(6)，根据Rsim得到Rmin和Rmax。

3)无限高斯混合模型的参数和超参数的估计

对无限高斯混合模型中模型的估计基于对参数的后验概率采样，为了得到 模型的后验概率，本发明首先定义各个参数的先验概率。

对于高斯组分中的方差σ²，选取逆Gamma作为方差的共轭先验

σ²|ω～inv_gamma(β,ω)  (2)

其中β和ω为形态参数和尺度参数。在本发明中，设置β＝1，并且假定超参数ω 的先验为无信息的Jeffery先验，如式(9)所示：

$f\left(\omega \right)\propto {\left|\det \right[I\left(\omega \right)\left]\right|}^{1/2}\propto \frac{1}{\omega }---\left(3\right)$

其中I(·)是Fisher信息矩阵，det[·]为求解的行列式操作。

对于高斯组分中的丰度_α，由于受限与ANC和ASC，本发明使用对称的狄 利克雷分布作为其共轭先验。

$f(\alpha /R)~\mathrm{Dirichlet}(1/R,...,1/R)=\frac{\Gamma \left(1\right)}{{\left[\Gamma (1/R)\right]}^R}\underset{j=1}{\overset{R}{\Pi }}{\alpha }_j^{1/R-1}---\left(4\right)$

对于高斯组分数目R，其取值范围为Rmin和Rmax，本发明采取等概率的 均匀分布作为R的先验。

$p(R=k)=\frac{1}{R_{\max }-R_{\min }},k\in (R_{\min },R_{\max })---\left(5\right)$

确定了参数的先验分布，使用似然函数与先验函数的到参数的后验模型， 分层贝叶斯中的参数和超参数之间的关系如图1所示，需要确定的参数为 θ＝{α,σ²,R,M_R}，其参数后验如(12)所示：

$\left(\begin{array}{c}\left(\theta \right|y)\propto \int f\left(y\right|\alpha ,{\sigma }^2,R,M)f\left(\alpha \right|R)p\left(M_R\right|R)\\ ...p\left(R\right)f\left({\sigma }^2\right|\omega )f\left(\omega \right)\mathrm{d\omega }\end{array}\right)---\left(6\right)$

通过计算，可以得到后验分布为

$\left(\begin{array}{c}f\left(\theta \right|y)\propto \frac{1}{(R_{\max }-R_{\min })}\times \frac{1}{g{\left(\alpha \right)}^{L/2}{\sigma }^{2R+L}}\\ ...\times \frac{\Gamma \left(1\right)}{\Gamma {(1/R)}^R}\underset{j=1}{\overset{R}{\Pi }}{\alpha }_j^{1/k-1}\times \exp (-\frac{{\left|\right|y-k\left(\alpha \right)\left|\right|}^2}{2{\sigma }^{2}g\left(\alpha \right)})\end{array}\right)---\left(7\right)$

3-1)初始化首次迭代的初始化参数

在初始化过程中，对于高斯组分R，在Rmin和Rmax范围内随机选取高斯 组分个数R(1)，对于高斯组分的均值向量，根据R(1)在高斯组分均值向量集合M 选取MR(1)作为对于丰度向量α⁽¹⁾，按照其先验分布(10)进行初始化。 设置超参数ω⁽¹⁾＝10^-2，并根据式(9)初始化参数σ²⁽¹⁾，设定“前进”策略被执 行的累计次数Z＝0，设定“后退”策略被执行的累计次数Λ＝0；

3-2)二状态策略(TTS策略)

由于无限高斯混合模型中高斯组分个数不确定，并且在每次迭代过程中会 产生变化，因此，本文设计了一种二状态策略(TTS)对 “前进”和“后退”的两种概率策略进行调整。

“前进”策略指的是组分个数R增1，“后退”策略指的是组分个数R减1。 对于前进策略，由于增加了新的组分，需要添加新的丰度，添加新的丰度值 φ'～Beta(1,R^(t))，根据对新的组分进行更新。对 于后退策略，由于减少了某个旧的组分，需要去除与之对应的旧的丰度，随机 选取需要剔除的高斯组分,根据对新的丰度 进行更新，其中rindex为选取将被剔除的高斯组分，F为除rindex之外的其他丰 度的求和$F={\Sigma }_{r\ne \mathrm{rindex}}{\alpha }_r^{\left(t\right)}.$

对于二状态策略，设定执行“前进”策略时的概率为执行“后退”策 略时的概率为不执行任何策略的概率为其中执行二 状态策略的具体过程如下所示：随机选取η₁～Uniform(0,1)，如果执 行“前进”操作，Z＝Z+1，如果执行“后退”操作，Λ＝Λ-1, 如果不执行任何操作，其中t为当前的迭代次数。

3-3)调整“前进”策略和“后退”策略

对于“前进”策略和“后退”策略的调整使求解高斯组分的个数在一定程度 上避免陷入局部的最优解。本发明记录了在t次迭代中执行“前进”的个数Z和 执行“后退”的个数Λ。

$f_{R^{\left(t\right)}}=\lambda f_{R^{(t-1)}}+(1-\lambda )\left(\frac{Z}{t}\right)---\left(8\right)$

$b_{R^{\left(t\right)}}=\lambda b_{R^{(t-1)}}+(1-\lambda )\left(\frac{\Lambda }{t}\right)---\left(9\right)$

${\mu }_{R^{\left(t\right)}}=\lambda {\mu }_{R^{(t-1)}}+(1-\lambda )\left(\frac{t-Z-\Lambda }{t}\right)---\left(10\right)$

3-4)参数和超参数采样

对于高斯组分丰度向量a，本发明采用Metropolis-within-Gibbs产生候选样 本，根据后验概率(13)得到a的条件概率分布为

f(α|y,R,σ²,M_R)∝f(y|α,σ²,R,M_R)f(α|R)  (11)

通过计算，(17)可以化简为

$\left(\begin{array}{c}f\left(\alpha \right|y,R,{\sigma }^2,M_R)\propto \frac{1}{g{\left(\alpha \right)}^{L/2}{\sigma }^L}\times \\ ...\frac{\Gamma \left(1\right)}{\Gamma {(1/R)}^R}\underset{j=1}{\overset{R}{\Pi }}{\alpha }_j^{1/k-1}\times \exp (-\frac{{\left|\right|y-k\left(\alpha \right)\left|\right|}^2}{2{\sigma }^{2}g\left(\alpha \right)})\end{array}\right)---\left(12\right)$

对于新端元向量中的元素使用如式(19)的概率作为是否接受或者拒绝产 生新的元素。

$\rho =\min \{\frac{f\left({\xi }_r\right|y,R,{\sigma }^2,M_R,{\alpha }_{\backslash r}^{\left(t\right)})}{f\left({\alpha }_r^{\left(t\right)}\right|y,R,{\sigma }^2,M_R,{\alpha }_{\backslash r}^{\left(t\right)})},1\}---\left(13\right)$

其中R为高斯组分个数，α^(t)是当前的丰度向量，σ^2(t)是方差，是除去的 元素，ξ_r是依据高斯分布产生的候选样本，因此得新的丰度值可以 根据(20)设置为

则最终的候选样本可以根据狄利克雷分布确定新的样本元素 ${\alpha }^{\left(t\right)}~\mathrm{Dirichlet}({\alpha }_1^{{\left(t\right)}^{\prime }},...,{\alpha }_R^{{\left(t\right)}^{\prime }})$

对于新的样本σ²，条件概率分布如

f(σ²|y,R,α,M_R,ω)∝f(y|α,σ²,R,M_R)f(σ²|ω)  (15)

因此，通过计算，(12)可以得到

${\sigma }^2|y,R,\alpha ,M_R,\omega ~\mathrm{invgamma}(\frac{L}{2}+1,\frac{{\left|\right|y-k\left(\alpha \right)\left|\right|}^2}{2g\left(\alpha \right)}+\omega )---\left(16\right)$

对于新的样本ω，新的超参数ω可以通过式(23)得到

$\omega |{\sigma }^2~\mathrm{gamma}(1,\frac{1}{{\sigma }^2})---\left(17\right)$

利用本发明中叙述的方法对采集到受侵染样本样本进行解混，能够区分被 农药侵染的叶片部位和被侵染的各个部分的权重。

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