CN104154998A - 一种基于压缩感知的计算多光谱成像图谱的重构方法 - Google Patents

Abstract

本发明公开了一种基于压缩感知的计算多光谱成像图谱的重构方法。本发明的方法基于由望远镜成像模块，数字微反射镜DMD及控制模块，光学汇聚透镜，滤光片轮，光电倍增管PMT，数据采集模块和多光谱图像重构模块组成的系统实现。系统按照预先所设置调制模板的数学形式调制目标场景的空间信息，再经由后续数据计算方法反演得到目标场景的多光谱图像。本发明的优点是：无需任何扫描，目标场景图像重构所需的数据量少，探测灵敏度高，结构简单。

Description

一种基于压缩感知的计算多光谱成像图谱的重构方法

技术领域

本发明涉及计算光谱成像技术以及多光谱图像的重构算法，信号处理等领域。特别涉及一种基于压缩感知的计算多光谱成像图谱的重构方法。

背景技术

多光谱，高光谱成像技术(参见文献1)以物质的光谱分析理论为基础，涉及光学系统设计、成像技术、光电探测、信号处理与信息挖掘、光谱信息传输理论、地物波谱特性研究等领域，是当前重要的空间对地观测技术手段之一。由于高光谱成像技术能同时获取目标场景的二维空间信息和光谱信息，因此在食品安全检测、地物目标的检测识别、土壤中的金属污染的检测、矿物种类分析、洪涝灾害预测、刑事侦查、艺术品诊断、军事应用等领域都有广泛的应用。

传统的高光谱成像仪的成像按光栅分光方式主要有光机扫描式、推帚式两种。光机扫描式成像光谱仪以光机扫描方式工作，扫描镜从刈幅的一端扫至另一端，从而使不同位置的地物目标发出的光进入光学系统成像。推帚式成像光谱仪以固体自扫描方式，使用面阵探测器单元对二维地物目标进行扫描，其空间维像元数与地面给定刈幅的采样元相同，光谱维像元数与给定光谱通道数相符。光机扫描式成像光谱仪由于具有扫描的运动部件，从而容易导致系统的不稳定性。同时光机扫描过程中穿轨空间像元不同时获取将导致图谱后处理难度极大。另外，探测积分时间短将极大地限制空间分辨率和光谱分辨率。推帚式成像仪由于其光学系统结构的原因，成像视场很小，同时定标非常难，光学系统复杂，存在光谱弯曲的情况。而且，由于两种成像方式都是采用光栅分光，单一狭缝的使用将会不可避免的带来空间分辨率和光谱分辨率之间的矛盾。

多光谱成像系统相比较高光谱成像，牺牲了光谱分辨率，但系统结构非常简单，在一些不需要高光谱分辨的场合仍然具有广泛的应用。传统多光谱成像大多采用面阵或者线阵探测器，利用滤光片实现多光谱成像。如果采用线阵或者单元探测器，需要在空间上进行扫描才能得到目标场景的完整像。就目前半导体技术而言，可见谱段的面阵探测器已经非常成熟。但红外谱段，尤其是中远红外谱段的国产化探测器还无法大规模集成。而且目前国际上最先进的红外谱段的面阵探测器仍然存在非均匀性等问题。

计算光谱成像技术(参见文献2)是近几年国内外提出的一门新兴的成像技术，相比较于传统的光谱成像系统,计算光谱成像使用特殊的空间光调制器作为编码孔径对景物目标进行编码,按照预先所设置编码孔径的数学形式调制、捕捉景物空间信息和光谱信息成像。最后基于压缩感知理论(参见文献3、4、5)，再经由后续数据计算方法反演得到最终多光谱图像。能够采用单像素探测器获取数据。将其应用在红外或者中远红外谱段具有重要的研究意义。

基于以上背景，本发明提出一种基于压缩感知的计算多光谱成像图谱的重构方法。

参考文献：

[1]王建宇，舒嵘，刘银年，马艳华.成像光谱技术导论.北京:科学出版社，2011.

[2]Gehm M E,John R,Brady J D,Willett M R,Schulz J T.2008 Opt.Express17 14013.

[3]Donoho D L.Compressed sensing[J].IEEE Transactions on InformationTheory，2006，52(4)：1289-1306.

[4]Candès E，Romberg J，Tao T.Robust uncertainty principles：exact signalreconstruction from highly incomplete frequency information[J].IEEETransactions on Information Theory，2006，52(2)：489-509.

[5]Candès E.Compressive sampling[C].International Congress ofMathematics，2006：1433-1452.

发明内容

本发明的目的是提供一种基于压缩感知的计算多光谱成像图谱的重构方法。在探测器方面，采用单元探测器，解决国产中远红外谱段探测器无法大规模集成以及探测器自身非均匀差等瓶颈问题。在数据获取方面，基于压缩感知理论，采用少量的数据即可重构得到目标的多光谱图像，在采样的过程中就已经压缩了数据，缓解传统多光谱成像中大数据量的采集、传输、存储压力。

本发明的方法使用数字微反射镜DMD作为编码孔径对景物目标进行编码,按照预先所设置调制模板的数学形式调制、捕捉景物空间信息和光谱信息成像。最后基于压缩感知理论，再经由后续数据计算方法反演得到最终目标场景的多光谱图像。

本发明提出的解决思路如下：

如图1所示，该发明的系统包括：望远镜成像模块1，数字微反射镜DMD及控制模块2，光学汇聚透镜3，滤光片轮4，单像素光电倍增管PMT 5，数据采集模块6，多光谱图像重构模块7。其特征在于：望远镜成像模块1采用的望远镜焦距为304.8mm，口径为101.6mm；数字微反射镜DMD及控制模块2采用的DMD其技术指标为：1024×768像素，像素大小为13.69微米；光学汇聚透镜3采用的焦距为40mm；滤光片轮4根据谱段需要组合不同的滤光片；单像素光电倍增管PMT 5工作谱段在300nm-900nm，像元大小8mm，输出最大暗电压20mv；数据采集模块6采用的采集卡其量化位数为16位，采样率250kS/s。

系统各模块之间的工作流程如下：

目标场景通过望远镜成像模块1成像于数字微反射镜DMD上。设定目标场景的空间信息被划分为m×n像素，光谱信息被划分为L个谱段，其图谱信息依次表示为φ(x,y,λ₁),φ(x,y,λ₂)…φ(x,y,λ_L)，这些函数的取值等于目标场景上特定像素点的对应谱段的强度值；

数字微反射镜DMD及控制模块2每次加载一个调制模板，通过调制模板改变DMD上每个小微镜的翻转状态从而达到调制目标场景的目的；

DMD第一次调制，令调制模板为θ₁(m,n)，具体取值为一个事先设定好的m×n阶随机矩阵，矩阵元素的取值为0或者1，所有的元素服从高斯随机分布。保持调制模板不变，调制后的场景经过光学汇聚透镜3汇聚，控制滤光片轮4，使其允许通过的谱段依次为λ₁,λ₂…λ_L。对于每一个谱段的光信号被单像素光电倍增管PMT 5探测后，经数据采集模块6采集，得到的电压信号依次用f₁(x,y,λ₁),f₁(x,y,λ₂)…f₁(x,y,λ_L)表示；

依次类推，DMD第k次调制，记调制模板为θ_k(m,n)，同第一次调制的过程，经数据采集模块6采集，得到的电压信号依次用f_k(x,y,λ₁),f_k(x,y,λ₂)…f_k(x,y,λ_L)表示；

基于压缩感知理论，本发明的方法提出：调制模板的个数即k，同时也就是观测矩阵的行数，取值由目标场景重构所需要的空间分辨率和目标场景的稀疏度决定，k远远小于N，取值范围为N为目标场景的总像素量，即N＝m×n；由于k远远小于N，因此对于每一个谱段，相当于用数据采集模块6采集到的少量k个信号值去重构目标场景的N个像素分辨率的图像。言外之意也就是说在数据采集模块6采集的过程中已经不知不觉的压缩了N-k个数据，从而缓解了传统多光谱成像中大数据量的采集、传输、存储压力；

因此，本发明系统的多光谱图像重构模块7对数据采集模块6采集到的信号进行处理，最后得到目标场景的多光谱图像；处理步骤如下：

1)对目标场景的第一个谱段的图像重构，将数据采集模块6采集到的信号整理写成如下(1)式：

$\left\{\begin{array}{c}{f}_1(x,y,{\mathrm{\λ}}_1)={\mathrm{\θ}}_1(m,n)\·\mathrm{\φ}(x,y,{\mathrm{\λ}}_1)+e_{11}\\ f_2(x,y,{\mathrm{\λ}}_1)={\mathrm{\θ}}_2(m,n)\·\mathrm{\φ}(x,y,{\mathrm{\λ}}_1)+e_{21}\\ \·\\ \·\\ \·\\ f_k(x,y,{\mathrm{\λ}}_1)={\mathrm{\θ}}_k(m,n)\·\mathrm{\φ}(x,y,{\mathrm{\λ}}_1)+e_{k1}\end{array}\right.---\left(1\right)$

上式中，e₁₁,e₂₁…e_k1为单像素光电倍增管PMT 5的噪声；

将(1)式用矩阵方程表示为如下(2)式：

F＝Θ·Φ+E  (2)

上式中，F是由信号f₁(x,y,λ₁),f₂(x,y,λ₁)…f_k(x,y,λ₁)组成的k×1矩阵；Θ为k×N矩阵，行数k即为调制次数，列数N＝m×n为调制模板θ_k(m,n)的元数个数，Θ的每一行由对应的θ_k(m,n)重新排列而成；Φ为由：φ(x,y,λ₁)组成的N×1矩阵；E为k×1噪声矩阵。显然，(2)式实际上是一个病态方程。直接求解很明显有无穷多个解。但压缩感知理论指出，只要Φ是稀疏的，或者在某种正交变换的表示下具有稀疏性，那么求解(2)式将会有特殊的优化求解方法。稀疏性的意思是指其中包含大量的趋于零的数据，只有少量的非零值；

对于自然场景，一般情况下，可以在一些正交变换的表示下具有稀疏性。例如：傅里叶变换，离散余弦变换等。对于Φ，在离散余弦变换下，将其稀疏表示为如下(3)式：

对于Φ，在离散余弦变换下，将其稀疏表示为如下(3)式：

Φ＝Ψ·α  (3)

上式中，α为Φ的稀疏表示，它是一个N×1矩阵，稀疏度为：S，也就是说α中仅含有S个非0值；Ψ是N×N阶离散余弦变换矩阵；

于是，可以将(2)式重新表示为如下(4)式所示：

F＝Θ·Φ+E＝Θ·Ψ·α+E＝T·α+E  (4)

上式中，T为k×N矩阵，在压缩感知理论中将其称为传感矩阵。(4)式中，只有α为未知数；

图像重构的方法就是求解(4)式中的稀疏系数α。很明显(4)式实际上是一个病态方程。直接求解有无穷多个解，因此将其转化为如下式(5)的优化问题：

$\widehat{\mathrm{\α}}=\arg \min {\left|\right|\mathrm{\α}\left|\right|}_{L_1},\mathrm{st}.F=T\·\mathrm{\α}---\left(5\right)$

上式中，L₁表示1范数，为α的最优近似解；

(5)式的优化求解算法步骤如下：

第一步：初始化一个空矩阵I＝[]，残差矩阵R＝F；

第二歩：将残差R与T中的每一列分别做内积，并找到内积最大的那一列，将本列取出并添加到矩阵I中；

第三歩：更新残差，R＝F-I·(I^T·I)^-1·I^T·F,其中I^T为I的转置矩阵(I^T·I)^-1为(I^T·I)的逆矩阵；

第四步：不断顺序循环第二歩和第三步，如果残差R满足：则退出循环，然后转到第五步。其中为矩阵R中的所有元素做平方然后求和，r表示预先设定的误差门限，一般取r<0.5；

第五步：最终(5)式求得的解为如下(6)式：

$\widehat{\mathrm{\&alpha;}}={(I^T\&CenterDot;I)}^{-1}\&CenterDot;I^T\&CenterDot;F---\left(6\right)$

最终求得的第一个谱段的图像信息表示为如下(7)式：

$\mathrm{\&Phi;}=\mathrm{\&Psi;}\&CenterDot;\widehat{\mathrm{\&alpha;}}---\left(7\right)$

将(7)式中的N×1阶矩阵Φ重新排列成m×n阶矩阵即可得到该谱段目标场景的二维像；

2)对于第二个谱段的图像重构，将数据采集模块6采集到的信号重新整理，也就是将步骤1)中的(1)式写成如下(8)式：

$\left\{\begin{array}{c}{f}_1(x,y,{\mathrm{\&lambda;}}_2)={\mathrm{\&theta;}}_1(m,n)\&CenterDot;\mathrm{\&phi;}(x,y,{\mathrm{\&lambda;}}_2)+e_{12}\\ f_2(x,y,{\mathrm{\&lambda;}}_2)={\mathrm{\&theta;}}_2(m,n)\&CenterDot;\mathrm{\&phi;}(x,y,{\mathrm{\&lambda;}}_2)+e_{22}\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ f_k(x,y,{\mathrm{\&lambda;}}_2)={\mathrm{\&theta;}}_k(m,n)\&CenterDot;\mathrm{\&phi;}(x,y,{\mathrm{\&lambda;}}_2)+e_{k2}\end{array}\right.---\left(8\right)$

上式中，e₁₂,e₂₂…e_k2为单像素光电倍增管PMT(5)的噪声；

同理，依次类推，对于第L个谱段的图像重构，将步骤1)中的(1)式写成如下(9)式：

$\left\{\begin{array}{c}{f}_1(x,y,{\mathrm{\&lambda;}}_L)={\mathrm{\&theta;}}_1(m,n)\&CenterDot;\mathrm{\&phi;}(x,y,{\mathrm{\&lambda;}}_L)+e_{1L}\\ f_2(x,y,{\mathrm{\&lambda;}}_L)={\mathrm{\&theta;}}_2(m,n)\&CenterDot;\mathrm{\&phi;}(x,y,{\mathrm{\&lambda;}}_L)+e_{2L}\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ f_k(x,y,{\mathrm{\&lambda;}}_L)={\mathrm{\&theta;}}_k(m,n)\&CenterDot;\mathrm{\&phi;}(x,y,{\mathrm{\&lambda;}}_L)+e_{\mathrm{kL}}\end{array}\right.---\left(9\right)$

上式中，e_1L,e_2L…e_kL为单像素光电倍增管PMT(5)的噪声；

所以，对第2到第L个谱段的图像重构，采用同步骤1)中后续的处理方法，最终即可得到目标场景的多光谱图像。

本发明的优点在于：

(1)本发明系统设计采用DMD调制目标场景的图谱信息，采用单像素光电倍增管PMT获取数据，探测灵敏度高，无需任何扫描，结构简单。

(2)本发明基于压缩感知理论，利用少量的采样数据即可复原完整的多光谱图像。从采样开始就大大压缩了采样数据，极大地改善了系统的数据采集、传输。存储压力。

附图说明

图1是一种基于压缩感知的计算多光谱成像系统及图谱重构方法。其中：1是望远镜成像模块；2是数字微反射镜DMD及控制模块；3是光学汇聚透镜；4是滤光片轮；5是单像素光电倍增管PMT；6是数据采集模块；7是多光谱图像重构模块；

图2是本发明具体实施系统示意图。

图3是基于图2的原理样机，利用本发明的谱图重构方法得到的多光谱图像。其中：A是目标场景；1是目标场景A重构后650nm谱段的图像；2是重构后540nm谱段的图像；3是重构后获得的470nm谱段的图像；4是用前面三个谱段合成的三谱段图像；5是目标场景的全景图。

具体实施方式

下面结合图1给出本发明的一个较好实例，主要作进一步详细说明，而非用来限定本发明的范围。

(1)首先确定系统各模块所用元器件的基本技术参数。具体如下：望远镜成像模块1采用爱蒙特光学(深圳)有限公司的望远镜，选定焦距为304.8mm，口径为101.6mm；数字微反射镜DMD及控制模块2采用美国TI公司生产的DMD，技术指标为1024×768像素，像素大小为13.69μm，控制板采用与之配合的TI-Discovery-4100；光学汇聚透镜3焦距为40mm；滤光片轮4根据需要组合三个谱段的滤光片，谱段分别为650nm，540nm，470nm；单像素光电倍增管PMT 5采用Hamamatsu Photonics K.K公司的H5784-01产品，工作谱段在300nm-900nm，像元大小8mm，输出最大暗电压20mv；数据采集模块6采用美国NI公司的NI-USB-6211采集卡，技术指标为：量化位数16位，采样率250kS/s。

如图2，基于图1的原理系统并利用上述所给的元器件搭建的系统。

(2)望远镜成像模块1将目标场景成像于DMD上。数字微反射镜DMD及控制模块2每次加载一个调制模板发送到DMD，调制模板是一个大小为256×192，并且服从高斯随机分布的0,1矩阵，这些调制模板事先生成好，本次实施生成6000个模板。然而由于步骤(1)所选取的DMD为1024×768像素，为了能让DMD正确识别，实际加载到DMD的模板选取大小为1024×768，超出256×192的元素全部用“0”补全。通过调制，使DMD的微镜处于一定的开关状态，其中“开”用“1”表示，“关”用“0”表示(具体使微镜翻转+12°(开)和-12°(关))。等微镜开关状态稳定以后，即可实现对目标场景的调制，翻转-12°反射的光线被丢弃，翻转+12°反射的光线接着进入后续的光学系统。

(3)DMD每加载一个调制模板，即相当于调制一次。每调制一次的同时DMD上处于“开”状态的微镜翻转+12°，使发射光线进入后续的光学系统通过滤光片轮4，然后组合不同谱段的滤光片实现对目标场景的多光谱成像。每一次调制的同时单像素光电倍增管PMT 5探测不同谱段的数据，并且数据采集模块6同步采集所有谱段的数据。

(4)上述实施步骤(1)(2)(3)重复6000次，即调制6000次。最终完成数据采集。最后多光谱图像重构模块7将采集到的所有数据，重新整理组合，并利用本发明提出的重构方法得到目标场景的多光谱图像。为了简要说明，目标场景被划分为256×192像素大小，谱段划分为650nm，540nm，470nm三个谱段。对每一个谱段重构，需要求解256×192个未知数，而采样数据只有6000个，相当于采样的过程中压缩了256×192-6000＝43152个数据。

如图3，基于图2的原理样机，利用本发明的谱图重构方法所得的多光谱图像。

Claims (1)

1.一种基于压缩感知的计算多光谱成像图谱的重构方法，它基于包括望远镜成像模块(1)，数字微反射镜DMD及控制模块(2)，光学汇聚透镜(3)，滤光片轮(4)，单像素光电倍增管PMT(5)，数据采集模块(6)和多光谱图像重构模块(7)的计算多光谱成像系统实现；其特征在于方法如下：

由控制模块(2)加载到数字微反射镜DMD上的调制模板依次为：

θ₁(m,n),θ₂(m,n)…θ_k(m,n)

其中：k为调制次数，θ_k(m,n)为高斯随机分布的m×n阶矩阵；

在每个调制模板调制过程中，由数据采集模块(6)采集，得到的k组数字信号依次为：

f_k(x,y,λ₁),f_k(x,y,λ₂)…f_k(x,y,λ_L)

其中：x,y为目标场景的二维空间信息坐标；λ_L为滤光片轮(4)允许通过的波段；

对上述数据具体处理步骤如下：

1)对目标场景的第一个谱段的图像重构，将数据采集模块(6)采集到的信号整理写成如下(1)式：

$\left\{\begin{array}{c}{f}_1(x,y,{\mathrm{\&lambda;}}_1)={\mathrm{\&theta;}}_1(m,n)\&CenterDot;\mathrm{\&phi;}(x,y,{\mathrm{\&lambda;}}_1)+e_{11}\\ f_2(x,y,{\mathrm{\&lambda;}}_1)={\mathrm{\&theta;}}_2(m,n)\&CenterDot;\mathrm{\&phi;}(x,y,{\mathrm{\&lambda;}}_1)+e_{21}\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ f_k(x,y,{\mathrm{\&lambda;}}_1)={\mathrm{\&theta;}}_k(m,n)\&CenterDot;\mathrm{\&phi;}(x,y,{\mathrm{\&lambda;}}_1)+e_{k1}\end{array}\right.---\left(1\right)$

上式中，e₁₁,e₂₁…e_k1为单像素光电倍增管PMT(5)的噪声；φ(x,y,λ₁)为待重构的目标场景的第一个波段的图像，像素大小为：m×n；

将(1)式用矩阵方程表示为如下(2)式：

F＝Θ·Φ+E  (2)

上式中，F是由信号f₁(x,y,λ₁),f₂(x,y,λ₁)…f_k(x,y,λ₁)组成的k×1矩阵；Θ为k×N矩阵，行数k即为调制次数，列数N＝m×n为调制模板θ_k(m,n)的元数个数，Θ的每一行由对应的θ_k(m,n)重新排列而成；Φ是由φ(x,y,λ₁)组成的N×1矩阵；E是由e₁₁,e₂₁…e_k1组成的k×1阶噪声矩阵；

对于Φ，在离散余弦变换下，将其稀疏表示为如下(3)式：

Φ＝Ψ·α  (3)

上式中，α为Φ的稀疏表示，它是一个N×1矩阵；Ψ是N×N阶离散余弦变换矩阵；

于是，可以将(2)式重新表示为如下(4)式所示：

F＝Θ·Φ+E＝Θ·Ψ·α+E＝T·α+E  (4)

上式中，T为k×N矩阵，(4)式中，只有α为未知数；

图像重构的方法就是求解(4)式中的稀疏系数α。将其转化为如下式(5)的优化问题：

$\widehat{\mathrm{\&alpha;}}=\arg \min {\left|\right|\mathrm{\&alpha;}\left|\right|}_{L_1},\mathrm{st}.F=T\&CenterDot;\mathrm{\&alpha;}---\left(5\right)$

上式中，L₁表示1范数，为α的最优近似解；

(5)式的优化求解算法步骤如下：

第一步：初始化一个空矩阵I＝[]，残差矩阵R＝F；

第二歩：将残差R与T中的每一列分别做内积，并找到内积最大的那一列，将本列取出并添加到矩阵I中；

第三歩：更新残差，R＝F-I·(I^T·I)^-1·I^T·F,其中I^T为I的转置矩阵(I^T·I)^-1为(I^T·I)的逆矩阵；

第四步：不断顺序循环第二歩和第三步，如果残差R满足：则退出循环，然后转到第五步。其中为矩阵R中的所有元素做平方然后求和，r为预先设定的误差门限，一般取r<0.5；

第五步：最终(5)式求得的解为如下(6)式：

$\widehat{\mathrm{\&alpha;}}={(I^T\&CenterDot;I)}^{-1}\&CenterDot;I^T\&CenterDot;F---\left(6\right)$

最终求得的第一个谱段的图像信息表示为如下(7)式：

$\mathrm{\&Phi;}=\mathrm{\&Psi;}\&CenterDot;\widehat{\mathrm{\&alpha;}}---\left(7\right)$

将(7)式中的N×1阶矩阵Φ重新排列成m×n阶矩阵即可得到该谱段目标场景的二维像；

2)对于第二个谱段的图像重构，将数据采集模块(6)采集到的信号重新整理，也就是将步骤1)中的(1)式写成如下(8)式：

$\left\{\begin{array}{c}{f}_1(x,y,{\mathrm{\&lambda;}}_2)={\mathrm{\&theta;}}_1(m,n)\&CenterDot;\mathrm{\&phi;}(x,y,{\mathrm{\&lambda;}}_2)+e_{12}\\ f_2(x,y,{\mathrm{\&lambda;}}_2)={\mathrm{\&theta;}}_2(m,n)\&CenterDot;\mathrm{\&phi;}(x,y,{\mathrm{\&lambda;}}_2)+e_{22}\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ f_k(x,y,{\mathrm{\&lambda;}}_2)={\mathrm{\&theta;}}_k(m,n)\&CenterDot;\mathrm{\&phi;}(x,y,{\mathrm{\&lambda;}}_2)+e_{k2}\end{array}\right.---\left(8\right)$

上式中，e₁₂,e₂₂…e_k2为单像素光电倍增管PMT(5)的噪声；φ(x,y,λ₂)为待重构的目标场景的第二个波段的图像，像素大小为：m×n；

同理，依次类推，对于第L个谱段的图像重构，将步骤1)中的(1)式写成如下(9)式：

$\left\{\begin{array}{c}{f}_1(x,y,{\mathrm{\&lambda;}}_L)={\mathrm{\&theta;}}_1(m,n)\&CenterDot;\mathrm{\&phi;}(x,y,{\mathrm{\&lambda;}}_L)+e_{1L}\\ f_2(x,y,{\mathrm{\&lambda;}}_L)={\mathrm{\&theta;}}_2(m,n)\&CenterDot;\mathrm{\&phi;}(x,y,{\mathrm{\&lambda;}}_L)+e_{2L}\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ f_k(x,y,{\mathrm{\&lambda;}}_L)={\mathrm{\&theta;}}_k(m,n)\&CenterDot;\mathrm{\&phi;}(x,y,{\mathrm{\&lambda;}}_L)+e_{\mathrm{kL}}\end{array}\right.---\left(9\right)$

上式中，e_1L,e_2L…e_kL为单像素光电倍增管PMT(5)的噪声；φ(x,y,λ_L)为待重构的目标场景的第L个波段的图像，像素大小为：m×n；

对第2到第L个谱段的图像重构，采用上述步骤1)的处理方法，最终得到目标场景的多光谱图像。