基于概率矩阵分解的高光谱图像锐化方法与流程

本发明涉及一种基于概率矩阵分解的高光谱图像锐化方法，属于遥感图像处理领域。

背景技术：

光谱图像是指对同一拍摄目标在不同光谱波段上形成的多幅图像集合。依照图像成像的光谱分辨率(光谱分辨率，指仪器所能识别的不同光波的波长间隔；高光谱分辨率，指所能分辨不同光波波长间隔小，对于不同的遥感卫星标准也不一样，高光谱分辨率的标准也不同，此处仅作泛指，下同)，光谱图像分为多光谱图像(可分辨的波长差为100nm)、高光谱图像(可分辨的波长差为10nm)和超光谱图像(可分辨的波长差为1nm)。目前，高光谱图像在遥感领域有着许多应用，许多遥感图像的分类识别任务希望图像的空间分辨率尽可能高(高空间分辨率，简称高分辨，指图像所能识别地面的物理距离非常小，单位区域内成像的像素点非常多，高分辨率对于不同的卫星标准也不同，这里泛指所有场景的高分辨图像要求，下同)，这样遥感卫星对地面目标定位识别越精确。然而，高分辨的高光谱图像获取困难：一方面该图像数据量大，无线传输成本高；另一方面，这种图像光学成像困难，图像的空间分辨率和光谱分辨率存在矛盾和折中，即，高空间分辨的图像，光谱分辨率相对较低，反之亦然。因此，对同一目标地点进行遥感，获得的高光谱图像虽然光谱分辨率很高，但是相对于多光谱图像而言，它的空间分辨率低。为了解决上述困难，人们提出一种技术：在遥感卫星上，对同一场景同时拍摄并传输一幅低分辨的高光谱图像和一幅高分辨的多光谱图像，在通信接收端将两幅图像融合生成一幅高分辨率的高光谱图像。融合生成的图像既具备了输入的多光谱图像的高空间分辨率优势，又具备了输入的高光谱图像的高光谱分辨率优势。这种技术既降低了传送数据量，又绕开光学成像难题。这种技术被称为高光谱图像和多光谱图像的融合，有时也被称为高光谱图像的锐化(pan-sharpening)。

迄今为止，针对高光谱图像的锐化算法层出不穷。然而已有算法大都存在以下缺陷：

1、现有算法锐化精度和计算耗时仍有待提升，大多数算法虽然具备高锐化精度，但是计算耗时成倍上升；

2、大多算法采用优化建模，模型存在许多参数，需要通过实验手动调节，一旦调节不当，锐化精度降低，且调节过程复杂耗时。

技术实现要素：

针对现有技术的缺陷，本发明的目的是提出一种基于概率矩阵分解的高光谱图像锐化方法，该方法从贝叶斯建模出发，通过变分法完成贝叶斯估计过程，推理过程能自动学习大多参数，在消耗较低计算资源下完成较高精度的图像锐化。

本发明的思路是：

输入两幅同拍摄条件同目标地点的不同分辨优势的图像——低分辨的高光谱图像和高分辨的多光谱图像，已知高光谱分辨率到多光谱分辨率的响应矩阵，求两幅图像融合生成的高分辨的高光谱图像。

本发明使用如下的符号：

记融合后的高分辨的高光谱图像为输入的两幅图像——低分辨的高光谱图像记为高分辨的多光谱图像记为输入的高光谱分辨率到多光谱分辨率的响应矩阵，即光谱响应矩阵，记为L和l分别表示两幅图像的光谱通道数(光谱通道数是指，在一个光谱波段范围内，按照一定的光波波长间隔进行成像所获得的总图像数。)，n和N分别表示两幅图像各光谱通道下的像素点数。显然l＜＜L，n＜＜N。

根据光学物理成像过程，可得到以下数学关系：

X＝ZB,Y＝FZ (1)

其中矩阵表示高分辨图像到低分辨图像的响应矩阵，该矩阵一般难以获得。根据第一节的问题描述，本发明针对的问题用数学语言可表达为：已知以及和的数学关系(见(1)式)，求

一般情况下，为求解(1)式对应的方程数远少于未知数的个数，数学上将这种问题称为欠定问题。显然，满足(1)式的矩阵有无穷多个，而这些解并不一定都是高分辨的高光谱图像。我们需要更多方程约束才能获得满足现实需求的解。由于代表一副图像，而图像在数据层面上具备一些平常数据不具备的统计学规律。如果我们能挖掘出这种统计学规律，我们就能获得额外的约束，用以求解较理想的融合图像。因此，如何对图像的统计学规律进行建模，成为问题求解的关键。

研究表明，高光谱图像中像素点的光谱矢量具有低维的结构特征。因此，我们做出假设：高分辨的高光谱图像中像素光谱矢量可认为由很少的隐藏光谱特征矢量线性叠加形成。表示为矩阵形式即：

Z＝U'T+R (2)

其中，为隐藏特征矢量构成的矩阵；为对应的线性叠加系数构成的矩阵；矩阵的行数r由使用者自行设定；为残余项，一般趋近于0。在这样的假设下，高分辨的高光谱图像的求解问题转化为对分解矩阵U，T和R的估计问题。

显然，如果准确估计分解矩阵将直接影响融合生成的高分辨的高光谱图像的质量。为了快速准确地估计分解矩阵，我们提出对上述问题进行详细建模，建模过程划分为三个阶段：第一个阶段，先对输入的低分辨的高光谱图像和高分辨的多光谱图像进行图像预处理，并分析处理后图像与待生成的高分辨高光谱图像之间的数学关系；第二阶段，我们借助第一阶段的数学关系分析进行贝叶斯建模；第三个阶段，对贝叶斯模型进行推理，为获得分解矩阵的后验概率分布，引入变分法，完成后验概率分布表达式的快速近似计算，并利用近似计算的后验概率分布表达式获得分解矩阵的均值，从而获得融合生成的高分辨的高光谱图像。

下面对三个阶段分辨进行详细分析：

(1)图像预处理及数学关系分析

由于低分辨的高光谱图像的空间分辨率远低于待求解的高分辨的高光谱图像为了便于数据处理，我们预先对进行插值，得到插值后的图像插值后的图像与的关系满足，其中矩阵表示插值过程等效的运算矩阵。结合(1)式，可知，为简写方便，令G＝BC,则结合(2)式，令W＝TG，N＝RG，可得如下式子：

另一方面，高分辨的多光谱图像Y中包含的成分，为了便于后续的模型建模和学习，我们对多光谱图像进行预处理，消去高分辨多光谱图像Y中的成分：由光学成像关系(1)可知，表示光谱响应矩阵，结合(2)式，令V＝T-W，易得：

残余项N的元素几乎为0，假设该矩阵中各元素服从0均值的高斯分布，则FN中光谱响应矩阵改变的N的分布特性，不便于建模。为了降低建模复杂度，需要对残差图像E进行变换，使变换后的噪声也服从各向同性的高斯分布，变换方法推导如下：令变换矩阵为为使ΦFN服从和N相同各向同性的高斯分布，根据高维高斯分布的性质，我们需要使ΦFF'Φ'＝I，I为单位阵。由奇异值分解原理(Singular Value Decomposition Method)可得，(FF')^-1＝QDQ'，其中Q为正交矩阵，满足QQ'＝Q'Q＝I，D为对角阵。显而易见，当ΦFF'Φ'＝I成立。因此，通过矩阵对图像矩阵E进行变换，将E中的各向异性高斯噪声转变为和N同分布的各向同性高斯噪声。变换后的图像记为图像矩阵令结合(3)式，可得：

通过上述的图像预处理与分析，我们得到两幅图像及变换矩阵：基于插值的高分辨的高光谱图像和消去高光谱图像成分后的高分辨多光谱变换图像另外，我们根据输入的光谱响应矩阵获得变换矩阵及变换后的光谱响应矩阵分解矩阵和是待求变量，其行数r需要根据实验手动设定。

通过数学分析，我们将矩阵T分解为W和V＝T-W。结合(2)式可得：

由于R非常小，一般情况下可认为趋近于0，因此R(I-G)≈0，对Z可通过下式进行估计，即：

通过上式易知，高分辨的高光谱融合图像可通过插值图像和分解矩阵乘积进行叠加。插值图像通过预处理即可获得，因此，对高分辨的高光谱融合图像的估计转化为对分解矩阵U和V的估计。本方法采样最小均方误差估计(MMSE估计)。为估计分解矩阵，我们需要在已有的数学关系基础上，建立贝叶斯模型，通过贝叶斯模型推导，获得分解矩阵的后验概率，并根据后验概率获得分解矩阵的均值，完成MMSE的整个估计过程。下面给出贝叶斯建模的过程。

(2)贝叶斯模型的建立

设矩阵N中的元素n_ij服从高斯分布，即由(2)和(4)式，可得和服从如下高斯分布：

为防止分解矩阵过拟合，将均值为零的各向同性高斯分布作为分解矩阵的先验约束，各分解矩阵先验分布如下所示：

令集合Ψ＝{α_n,α_u,α_v,α_w}，该集合的变量起到决定高斯分布的作用，因此数学上将这些变量称为模型的超参数，它们分别表示(9)式高斯分布的精度，为减少参数的手动调节，使模型超参数在贝叶斯推理过程中自动学习，我们引入超先验约束超参数(超参数的先验概率分布约束称为超先验约束)。超先验分布如下所示：

α_x～Γ(α_x|a_x,b_x),α_x∈Ψ (10)

上式中Γ(α_x|a_x,b_x)表示Gamma分布，其概率密度函数为：

超参数中的系数a_x和b_x一般取值非常小，从而使它们不影响后续的贝叶斯推理过程，在本方法中，这些参数取值为a_x＝b_x＝10^-6。

上述的贝叶斯模型可以通过有向概率图表示，见附图。通过概率图，易知，模型联合概率密度分布函数可表示为：

(3)变分贝叶斯方法

回顾建模的目的，我们需要估计分解矩阵U和V。在贝叶斯框架下，估计方法分为最大后验估计法(MAP估计)和最小均方误差估计法(MMSE估计)。前者以后验概率函数的峰值点作为估计结果，后者以后验概率的均值作为估计结果。本方法采用最小均方误差估计法，将分解矩阵的均值作为估计量，从而获得重建的融合图像。下面给出分解矩阵的后验概率推导公式：

通过(11)式可知，分母的积分十分难求，显然，在该公式下难以获得显示表达式。因此，退而求其次，我们引入变分法，近似估计上述后验概率。

变分贝叶斯方法的原理是假设近似估计的后验概率函数可分解为以下形式，其中q(·)表示对应变量的条件后验概率：

为了缩小估计后验概率函数和真实后验概率函数之间的误差，我们需要最优化以下表达式：

上式中，表示和的Kullback–Leibler距离，通过优化函数变量q(Ξ),Ξ∈Ω和q(α),α∈Ψ，以缩短估计的后验概率和真实后验概率的Kullback–Leibler距离。

通过变分分析，容易得到(推导过程见书Matthew J.Beal,Variational algorithms for approximate Bayesian inference,University of London,2003.)：

上式中表示对应表达式均值计算，将集合{Ω,Ψ}中的变量逐个代入到上式，并进行化简，可得到各个函数变量q(Ξ),Ξ∈Ω和q(α),α∈Ψ的计算表达式：

①q(W)的计算：

②q(V)的计算：

③q(U)的计算：

④q(α_n)的计算：

⑤q(α_u)的计算：

⑥q(α_v)的计算：

⑦q(α_w)的计算：

通过设定迭代次数Δ，反复完成①～⑦式的计算，即可获得分解矩阵和的均值，最终通过(6)式，完成高分辨的高光谱融合图像的估计。

本发明的特征在于所述高光谱图像锐化是指在通信接收端把输入的一幅低分辨的高光谱图像和一幅高分辨的多光谱图像融合生成一幅高分辨的高光谱图像，从而折中地解决遥感领域中图像的空间分辨率和光谱分辨率矛盾的一种方法。

本发明提出的基于概率矩阵分解的高光谱图像锐化方法，与已有锐化方法相比，其优点是：

1、本发明对输入高光谱图像和多光谱图像进行预处理。这种预处理将相对低频分量和高频分量分开存于两幅图片，便于后续进行有效的信息提取和学习。

2、本发明能根据输入图像的特点自动学习模型中大多参数的取值，避免对模型参数过多的手动调节。

3、本发明采用变分贝叶斯方法，使模型在很少的迭代次数中达到收敛，降低了计算消耗，实验证实锐化后的图像仍然具备较高准确度。

附图说明

图1是本发明的程序原理框图。

图2是本发明的软件核心结构图。

图3是方法的贝叶斯模型概率图。其中各字母表示建模过程采用的矩阵变量，箭头表示始端随机变量决定末端随机变量。该贝叶斯模型建模原理详见具体实施方法第二节。

图4是实施举例的两幅光谱图像：(a)图代表低分辨的高光谱图像，(b)图代表高分辨的多光谱图像。低分辨的高光谱图像长、宽、光谱通道数为64×64×188，高分辨的多光谱图像长、宽、光谱通道数为：512×512×6，详细介绍见具体实施方式。

图5是为两幅高分辨的高光谱图像：(a)为真实高分辨的高光谱图像，(b)为本发明合成的高分辨的高光谱图像(右)，这两幅图像的长、宽、光谱通道数相同，均为512×512×188。

图6是四副高光谱图像对应的彩色显示图像，(a)图代表真实高分辨的高光谱图像，(b)图代表本方法合成的高分辨的高光谱图像、(c)图代表CNMF方法合成的高分辨的高光谱图，(d)图代表HySure方法合成的高分辨的高光谱图。通过观察易得，(b)图比其他方法合成的图像更接近(a)图的真实图像，即本方法合成的图像精度更高。相关方法的介绍见具体实施方式第二节实施举例。

具体实施方式

一、实施步骤

本发明提出的基于概率矩阵分解的高光谱图像锐化方法，包括以下步骤：

步骤(1)：输入要求参数，实现计算机初始化

用一台机载可见光/红外成像光谱仪对同一目标拍摄一幅低分辨的高光谱图像简称为X，以及一幅高分辨的多光谱图像简称为Y，其中：

L为所述图像X的光谱通道数，l为所述图像Y的光谱通道数，l＜＜L,

n为所述图像X的各光谱波段的成像像素点数，N为所述图像Y的各光谱波段的成像像素点数，n＜＜N；

输入多光谱摄像头对应高光谱图像的频率响应矩阵简称光谱响应矩阵F；

手动设置分解矩阵维数r和算法迭代次数Δ，其中参数r控制分解矩阵的行数；

令所述融合后的高分辨率的高光谱图像为简称为待求图像Z，光学成像原理表明：X＝ZB，Y＝FZ，其中为高分辨图像到低分辨图像的响应矩阵；

步骤(2)：对所述图像X用线性插值运算矩阵C插值，且是一个已知固定常数矩阵。由于X＝ZB，因此，得到的插值图像满足即所述插值图像是待求图像Z经分辨率的降低后插值形成的图像；

步骤(3)：在下述假设和推导下，待求图像Z可近似为：

假设：所述图像Z中的像素光谱矢量是由少量隐藏的光谱特征矢量经过线性叠加形成，这些光谱特征矢量为矩阵U的每一行，像素光谱矢量对应的叠加系数为矩阵T的每一行，则有Z＝U'T+R，其中r表示分解矩阵U的行数，也表示隐藏光谱特征矢量的个数，该参数在步骤(1)中人为预先设定；残余项趋近于0，故设R为各向同性的高斯噪声；

令G＝BC,W＝TG,V＝T-W,N＝RG,则若近似认为所述残余项R满足R≈N，则待求图像Z满足其中

步骤(4)：按以下步骤计算所述图像Y的残差变换图像其中

步骤(4.1)：对所述图像Y进行预处理，消去其中只含在所述插值图像中的高光谱成分，得到残差图像矩阵

光学成像原理表明：Y＝FZ，所以

步骤(4.2)：通过残差图像变换矩阵对所述残差图像矩阵E进行变换，以便把所述残差图像E中各向异性高斯噪声转变为和所述残余项同分布的各向同性高斯噪声，其中：分解矩阵Q和分解矩阵D通过对矩阵(FF')^-1进行奇异值分解获得，Q满足QDQ'＝(FF')^-1，D为对角矩阵，表示对对角矩阵D中元素进行算术平方根运算。通过残差图像变换矩阵Φ，可得到变换后的光谱响应矩阵以及残差变换图像

步骤(5)：根据计算所得的残差变换图像和线性插值图像用变分贝叶斯方法迭代求取待求图像Z的分解矩阵U,V和W，得到U'V，从而估算出待求图像Z，所用变分贝叶斯方法步骤如下：

步骤(5.1)：根据分解矩阵的维数r、所述图像Y的各光谱波段的成像像素点数N和所述图像X的光谱通道数L，见步骤(1)，确定分解矩阵和的大小,在此基础上初始化下列参数：

1_r×N表示大小为r×N的全1矩阵；

I_r表示大小为r×r的单位矩阵；

步骤(5.2)：总迭代次数用Δ表示，第t次迭代的各分解矩阵的取值表示为(·)^(t)，符号的上划线表示符号的概率意义上的均值，t＝1,2,3,...,Δ；

步骤(5.3)：输入第t-1次迭代，即上一次迭代获得的参数按下式计算第t次迭代参数：和

其中I_L表示大小为L×L的单位矩阵；I_Lr表示大小为Lr×Lr的单位矩阵；L为所述图像X的光谱通道数；r为步骤(1)输入的分解矩阵的维数；表示两个矩阵的Kronecker积；vec[·]表示矩阵的矢量化表示形式；为一种运算，定义如下：

其中a_ij表示矩阵A第i行第j列的元素；当r＝1时，该运算值为a_i，j；否则为一个大小为r×r的矩阵；

步骤(5.4)：输入第t-1次迭代(上一次迭代)的获得的参数以及步骤(5.3)中第t次迭代得到的参数根据下列式子，计算获得和

步骤(5.5)：输入第t-1次迭代，即上一次迭代的获得的参数以及步骤(5.3)中第t次迭代得到的参数按下式计算和

步骤(5.6)：输入步骤(5.3)～步骤(5.5)所计算获得的所有参数，按下式计算和

其中，||·||_F表示矩阵的Frobenius范数，tr{·}表示矩阵的迹运算，二者定义见张贤达著《矩阵分析与应用》或者说明书内容详细介绍，N、L和r的定义见步骤(1)；

步骤(5.7)：根据步骤(5.2)～步骤(5.6)的顺序，重复迭代Δ次，得到和再按下式计算输出锐化后的高分辨率的高光谱图像

二、实施举例

1、数据来源和评价方法

本方案的测试图像来源于美国航空航天局AVIRIS(Airborne Visible Infrared Imaging Spectrometer,机载可见光/红外成像光谱仪)免费数据库，示例光谱图像(见附图)拍摄于美国Indian Pine北部的一块农业区，成像波段范围在370nm～2508nm，光谱分辨率达到10nm，光谱通道数为224；由于空气中的水分能吸收波长为480,560,660,830,1650and 2220nm的光，造成第1-2,105-115,150-170和223-224光谱通道成像质量差，因此实际使用图像仅采用除这36个光谱通道以外其他通道的图像，即实际使用图像的光谱通道数为188个。通过模拟美国陆地卫星USGS/NASA Landsat 7的成像系统，可获得高分辨的多光谱图像和低分辨的高光谱图像。采用本发明方法将获得的低分辨的高光谱图像和高分辨的多光谱图像融合，将融合后的高分辨的高光谱图像和真实的高光谱图像进行做差对比，以两幅图像的均方根误差作为评价指标(均方根误差，英文名称：Root Mean Square Error，英文简写RMSE，该误差作为衡量算法性能高低标准，误差越小，算法性能越高)。RMSE的计算式如下所示：

2、实施效果和方法对比

表1本方法和当前其他主流方法的性能对比

采用上述数据核评价方法，此处给出本方法和当前其他高光谱图像锐化算法的对比，对比结果如表1所示，输入的低分辨的高光谱图像和高分辨的多光谱图像见附图，真实的高分辨的高光谱图像和本算法合成的高分辨高光谱图像见附图，本算法生成的图像和其他方法生成的图像对比见附图。

在表1中，其他三种方法分别出处如下：

CNMF方法出自“N.Yokoya,T.Yairi,and A.Iwasaki,“Coupled nonnegative matrix factorization unmixing for hyperspectral and multispectral data fusion,”IEEE Trans.Geosci.and Remote Sens.,vol.50,no.2,pp.528–537,2012.”；

HySure方法出自“M.J.Bioucas-Dias,L.B.Almeida,et al.A convex formulation for hyperspectral image superresolution via subspace-based regularization[J].IEEE Transactions on Geoscience and Remote Sensing,2015,53(6):3373-3388.”。

通过表1和附图可知，本发明的方法在消耗较低的计算时间下，融合生成的高分辨的高光谱图像精度远高于当前其他方法。