CN103674244B - 一种基于gpu的快速反卷积二维光纤光谱抽取方法 - Google Patents

Abstract

本发明公开了一种基于GPU的快速反卷积二维光纤光谱抽取方法，该方法能快速准确地从二维光纤光谱图像中抽取一维光谱流量。包括如下步骤：1)通过定标灯光谱图像获得点扩展函数的形状和大小的初步估计，并基于估计和实验确定点扩展函数的高斯模型；2)通过与目标光谱图像相对应的信噪比高的平场光谱图像进行光谱追迹，确定轨迹中心；3)通过点扩展函数的高斯模型和光谱轨迹中心，确定卷积矩阵，采用CUDA并行编程环境下的压缩稀疏行存储格式进存储和运算；4)确定二维光谱成像卷积方程；5)采用共轭梯度最小二乘法反卷积求解二维光谱成像卷积方程，得到一维光谱。对反卷积迭代求解过程中计算最耗时的稀疏矩阵与向量相乘，采用GPU进行加速；6)将获得的一维光谱，通过通用接口，保存为灵活图像传输系统格式文件。

Description

一种基于GPU的快速反卷积二维光纤光谱抽取方法

技术领域

本发明属于计算机图像处理技术领域，特别涉及一种基于GPU的快速反卷积二维光纤光谱抽取方法。

背景技术

光纤光谱，或者称为二维光纤光谱：1966年，高琨发明了可以长距离高速通信的光纤，1969年贝尔实验室的威拉德·博伊尔和乔治·史密斯发明了电荷耦合元件(Charge-Coupled Device，CCD)。从此之后，科学家利用光纤把望远镜焦面上天体目标的像引导到光谱仪，通过光谱仪把天体的像色散到二维CCD上，CCD精确的记录了天体目标每个波长处的流量值。

定标灯光谱，一般通过拍摄特制的定标灯光获得。定标灯一般是充有某种特殊元素的特制光源，其光谱只含有元素受激发所产生的已知波长的谱线，连续谱能量基本为零。

平场光谱，一般通过拍摄特制的平场灯光获得。平场灯一般是光谱能量在整个观测波段上分布比较均匀！光谱形状比较平滑的特制光源。在没有合适的平场灯的情况下，也可以拍摄晨昏天光代替平场光谱。

抽谱方法是针对如何把CCD记录的二维光纤光谱图像转化成一维的能被天文学家使用的一维波长-流量光谱的方法。

目前，有四种光纤光谱抽取方法。

第一种是孔径光纤光谱抽取方法，由de Boer和Snijders于1981年提出。具体做法是：首先确定光谱图像的轨迹，也就是通过每一行光谱轮廓中心的曲线，一般是用低次多项式曲线表示出；再次，在给定的波长处(该波长对应着CCD上的某一行)，光谱的空间方向，轨迹的两侧，取一个固定大小的孔径，然后简单地把孔径内所有CCD像素的计数沿空间方向相加，得到的值便认为就是该行对应的波长处的流量值。如果把CCD每一行的流量值顺序排列在一起，就得到了一维光谱。该方法的缺点是无法解决紧挨着的两条光纤之间交叉污染问题和噪声对光纤的影响问题。

第二种是优化的孔径抽谱方法，由Horne和Robertson于1986年分别提出。具体做法是：根据孔径内每个像素的信噪比，给每个像素上的计数一个权重，信噪比高的像素权重高，否则权重低，然后最优化整体信噪比，来抽取流量的算法。由于引入了权重，计数值小的像素权重较小，计数值大的像素权重大，因此在没有交叉污染的情况下，孔径选的越大越好，并且不会引入很大的噪声。该方法能提高信噪比，但不能有效减少紧挨着的两条光纤之间交叉污染。

第三种是轮廓拟合方法，由Piskunov和Valenti于2002年提出。它来源于优化的孔径方法，但它要比优化的孔径方法简单、快速、便捷、容易使用。该方法认为光谱图像空间轮廓具有某种数学解析表达式的形式，例如高斯函数、Voigt函数、B样条轮廓等等。它能解决紧挨着的两条光纤之间交叉污染问题，但缺点是抽取速度非常慢，并且严重依赖于高信噪比的光纤光谱。

第四种是反卷积方法，它是一种经典的图像处理方法，在图像复原领域得到充分的应用。但直到2010年，才由Bolton和Schlegel应用于模拟二维光纤光谱的抽取。这个方法根据光谱图像的成像原理来抽取光谱，Bolton和Schlegel通过模拟二维光纤光谱数据验证了通过反卷积方法进行抽谱的有效性，但由于真实二维光纤光谱相比模拟二维光纤光谱数据量巨大，将反卷积方法应用于真实二维光纤光谱的抽取会遇到存储和计算的瓶颈，而难以实现。因此，本发明提出一种基于GPU的快速反卷积二维光纤光谱抽取方法。

GPU的英文全称是Graphic Processing Unit，中文翻译为“图形处理器”。GPU通用计算技术在浮点运算、并行计算等方面表现的很好，可以提供数十倍乃至于上百倍于CPU的性能。GPU通用计算方面的标准目前有OPEN CL、CUDA、ATI STREAM。CUDA是一种由NVIDIA推出的通用并行计算架构，该架构使GPU能够解决复杂的计算问题。它包含了CUDA指令集架构(ISA)以及GPU内部的并行计算引擎。开发人员现在可以使用C语言来为CUDA^TM架构编写程序，所编写出的程序可以在支持CUDA^TM的处理器上以超高性能运行。

发明内容

本发明的目的是解决真实二维光纤光谱大数据量的难以计算以及反卷积方法光谱抽取速度慢等问题。

为此，本发明公开了一种基于GPU的快速反卷积二维光纤光谱抽取方法，该方法能快速准确地从二维光纤光谱图像中抽取真实一维光谱，包括如下步骤：

1)通过定标灯光谱图像获得点扩展函数的形状和大小的初步估计，并基于估计和实验确定点扩展函数的高斯模型；

2)通过与目标光谱图像相对应的信噪比高的平场光谱图像进行光谱追迹，确定轨迹中心；

3)通过点扩展函数的高斯模型和光谱轨迹中心，确定点扩展函数扩展卷积的矩阵表示形式，采用CUDA并行编程环境下的压缩稀疏行存储(Compressed Sparse Row，CSR)格式进存储；

4)确定二维光谱成像卷积方程；

S＝AC+ε

方程中S是向量形式表示的二维光谱图像，A是卷积矩阵，C是待抽取的一维光谱，ε是光谱图像上的高斯白噪声。

5)采用共轭梯度最小二乘法(conjugate gradient least-squared method，CGLS)反卷积求解二维光谱成像卷积方程，得到一维光谱，对反卷积迭代求解过程中计算最耗时的卷积矩阵与向量相乘，采用GPU进行加速；

6)将获得的一维光谱，通过通用接口，保存为灵活图像传输系统格式(Flexible Image Transport System，FITS)文件。

附图说明

图1为基于GPU的快速反卷积二维光纤光谱抽取方法的流程图。

图2为定标灯光谱图像示例图。

图3为平场光谱图像示例图。

图4为目标光谱图像示例图。

图5为反卷积二维光纤光谱抽取与孔径抽取方法部分结果比较。

图6为本发明应用于定标灯光谱抽取结果实例图A。

图7为本发明应用于定标灯光谱抽取结果实例图B。

图8为本发明应用于平场光谱抽取结果实例图A。

图9为本发明应用于平场光谱抽取结果实例图B。

图10为本发明应用于目标光谱抽取结果实例图A。

图11为本发明应用于目标光谱抽取结果实例图B。

具体实施方式

下面结合附图对本发明作进一步详细的说明。

图1为基于GPU的快速反卷积二维光纤光谱抽取方法的流程图。

下面通过图1将本发明基于GPU的快速反卷积二维光纤光谱抽取方法进行详尽的描述。

步骤1：通过定标灯光谱图像获得点扩展函数H的形状和大小的初步估计，并基于估计和实验确定点扩展函数的高斯模型如下：

$H_l(x,y)=\frac{\mathrm{\α}}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_x{\mathrm{\σ}}_y}{e}^{-(\frac{{(x-x_0)}^2}{2{\mathrm{\σ}}_x^2}+\frac{{(y-y_0)}^2}{2{\mathrm{\σ}}_y^2})}---\left(1\right)$

其中，l是波长位置，α为流量控制系数，x，y分别是二维光纤光谱图像的横坐标和纵坐标，σ_x和σ_y分别是水平方向和垂直方向的标准差，本发明中分别设为0.6和0.4；x₀和y₀分别是轨迹中心的横坐标和纵坐标。图2为定标灯光谱图像示例图。

步骤2：通过与目标光谱图像相对应的信噪比高的平场光谱图像进行光谱追迹，确定轨迹中心(x₀，y₀)。其过程如下：提取与目标光谱相对应的平场光谱图像每条光谱的轮廓，将每个轮廓的峰值处定为光纤光谱在该处的轨迹中心。图3为平场光谱图像示例图。图4为目标光谱图像示例图。

步骤3：通过点扩展函数的高斯近似模型和光谱轨迹中心，确定卷积矩阵，采用CUDA并行编程环境下的压缩稀疏行存储(Compressed Sparse Row，CSR)格式进存储。具体过程为：将点扩展函数的高斯模型，用矩阵H表述如下：

$H=\left(\begin{array}{ccccc}{h}_{(-m,m)}& ...& h_{(-m,0)}& ...& h_{(-m,-m)}\\ .& & .& & .\\ .& & .& & .\\ .& & .& & .\\ h_{(0,m)}& ...& h_{\left(\mathrm{0,0}\right)}& ...& h_{(0,-m)}\\ .& & .& & .\\ .& & .& & .\\ .& & .& & .\\ h_{(m,m)}& ...& h_{(m,0)}& ...& h_{(m,-m)}\end{array}\right)---\left(2\right)$

其中，H的大小为(2m+1)，m为一个正整数，本发明中设定为4，h_(0，0)是点扩展函数H的矩阵的中心。

本发明采用了边界条件为0的假定，卷积矩阵A由点扩展函数H矩阵扩展。矩阵A具有Toeplitz矩阵形式，如下：

其中，A是块Toeplitz矩阵，并且T也是Toeplitz矩阵。以T₍₀₎作为例子，并去除符号表示中的相同下标0，T₍₀₎的形式如下：

对矩阵A采用CSR格式进存储的方式为：分别将列指针和非零值显式的存储在索引数组和数据数组中，行指针数组的存储也采取压缩稀疏行存储格式；对于一个M×N稀疏矩阵，行指针数组长度为M+I，并将第i行的偏移量存储在行指针数组中的第i个位置，行指针数组的最后一个值将另外对应于第M+1行，存储矩阵中非零值的数据。

步骤4：确定二维光纤光谱卷积成像方程，其推导过程如下：

设S为向量形式表示的二维光纤光谱图像，当该二维光纤光谱图像中只有1条光纤源成像结果时，设该二维光纤光谱图像的所有像素的个数为N，光谱图像的行数(行的像素个数)为M，S中的第t个像素表示为s_i，点扩展函数H在波长为λ处的第i个像素值表示为A_iλ，待抽取的一维光谱在波长为λ处的流量表示为C_λ，则S中所有像素值计算如下：

$S=\{s_1,s_2,...,s_t,...s_N\},s_t=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\underset{\mathrm{\λ}=1}{\overset{M}{\mathrm{\Σ}}}{A}_{\mathrm{i\λ}}{C}_{\mathrm{\λ}}+{\mathrm{\ϵ}}_i---\left(5\right)$

公式(2)用矩阵形式表达如下：

S_(N×1)＝A_(N×M)C_(M×N)+ε_(N×1)，      (6)

当该二维光纤光谱图像是由D条光纤源成像得到的结果时，设点扩展函数H在波长为λ处的第i个像素值表示为A′_iλ，抽取的一维光谱在波长为λ处的流量表示为C′_λ，则S中所有像素值计算如下：

$S=\{s_1,s_2,...,s_t,...s_N\},s_t=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\underset{\mathrm{\λ}=1}{\overset{M\×D}{\mathrm{\Σ}}}{A_{\mathrm{i\λ}}^{\′}{C}_{\mathrm{\λ}}^{\′}+\mathrm{\ϵ}}_i---\left(7\right)$

公式(7)用矩阵形式表达如下：

S_(N×1)＝A′_(N×M×D)C′_(M×D×1)+ε_(N×1)      (8)

在公式(5)到公式(8)中，S是向量表示的二维光谱图像，A和A′是卷积矩阵，C和C′是待抽取的一维光谱，ε是光谱图像上的高斯白噪声，角标括号内符号表示矩阵的大小，符号N是所有像素的个数，M是光谱图像的行数(行的像素个数)，D是光谱图像中的被抽取光纤光谱源的个数。

步骤5：采用共轭梯度最小二乘法(conjugate gradient least-squared method，CGLS)反卷积求解二维光纤光谱卷积成像方程，抽取一维光谱C。对求解过程每次迭代中计算最耗时的矩阵与向量相乘，采用GPU进行加速，具体的加速实现为：

在CUDA并行编程环境下实现y＝Ax的矩阵与向量相乘的并行计算，需要对计算任务进行适当的划分，使之适合CUDA编程模型的特点，充分利用GPU提供的资源，GPU上的任务分配方式如下：

每个线程(thread)负责计算计算输出向量的一个元素，即每个线程负责计算矩阵每行非零元素与向量乘积并求和，这种方式是串行算法到GPU上的直接映射。利用线程之间的并行执行实现并行计算，一般不满足对存储器的合并访问；

每个线程束(warp)负责计算输出向量的一个元素，即每个warp中的线程负责计算矩阵每行非零元素与向量的乘积，然后对warp的中间结果进行规约求和。这种方式满足对存储器的部分合并访问，利用同一个warp内的所有线程总是同步执行的特点，规约操作之前不需要进行同步操作来保证结果的正确性；

每个线程块(block)负责计算输出向量的一个元素，即矩阵每行非零元素与向量的乘积由一个block负责计算，然后对同一个block得到的中间结果进行规约求和。在矩阵与向量乘法计算中，每行元素之间的计算互不相关，不需要进行通信，而在CUDA编程模型中，block之间并行执行，只能通过代价很高的原子操作来实现同步，通过全局存储器进行通信，将一个block映射到矩阵每行元素的计算上，这种方式满足对存储器的部分合并访问，规约之前必须对同一个block中的线程进行同步操作。

步骤6：将获得的一维光谱，通过通用接口，保存为灵活图像传输系统格式(Flexible Image TransportSystem，FITS)文件。

采用GPU加速的共轭梯度最小二乘法并行求解，相比采用CPU的共轭梯度最小二乘法串行求解，最优的加速比为3.52。

图5为本发明与孔径抽取方法部分结果比较。图中从左至右分别是灯谱示例图，本发明抽谱效果图，孔径方法抽谱效果图，本发明抽取残差图和孔径方法抽取残差图。图6至图11为应用本发明抽取光谱结果实例图。

总之，本发明的实施例公布的是其较佳的实施方式，但并不限于此。本领域的普通技术人员极易根据上述实施例，领会本发明的精神，并做出不同的引申和变化，但只要不脱离本发明的精神，都在本发明的保护范围之内。

Claims (4)

1.一种基于GPU的快速反卷积二维光纤光谱抽取方法，其特征在于，包括以下步骤：

1)通过定标灯光谱图像获得点扩展函数的形状和大小的初步估计，并基于估计和实验确定点扩展函数的高斯模型；

2)通过与目标光谱图像相对应的信噪比高的平场光谱图像进行光谱追迹，确定轨迹中心；其过程如下：提取与目标光谱相对应的平场光谱图像每条光谱的轮廓，将每个轮廓的峰值处定为光纤光谱在该处的轨迹中心；

3)通过点扩展函数的高斯模型和光谱轨迹中心，确定点扩展函数卷积的矩阵表示形式，采用CUDA并行编程环境下的压缩稀疏行存储格式进行存储和运算；

4)确定二维光谱成像卷积方程；

5)采用共轭梯度最小二乘法反卷积求解二维光谱成像卷积方程，得到一维光谱，对求解二维光谱成像卷积方程的过程中计算最耗时的点扩展函数矩阵与向量相乘，采用GPU进行加速；

6)将获得的一维光谱，通过通用接口，保存为灵活图像传输系统格式文件。

2.根据权利要求1所述的基于GPU的快速反卷积二维光纤光谱抽取方法，其特征在于，步骤1)中所述确定点扩展函数H的高斯模型为

$H_l(x,y)=\frac{\mathrm{\α}}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_x{\mathrm{\σ}}_y}{e}^{-(\frac{{(x-x_0)}^2}{{2\mathrm{\σ}}_x^2}+\frac{{(y-y_0)}^2}{{2\mathrm{\σ}}_y^2})}---\left(1\right)$

其中，l是波长位置，α为流量控制系数，x，y分别是两维光纤光谱图像的横坐标和纵坐标变量，σ_x和σ_y分别是水平方向和垂直方向的标准差，x₀和y₀分别是轨迹中心的横坐标和纵坐标。

3.根据权利要求1所述的基于GPU的快速反卷积二维光纤光谱抽取方法，其特征在于，步骤4)中所述确定二维光谱成像卷积方程的过程如下：

设S为向量形式表示的二维光纤光谱图像，当该二维光纤光谱图像中只是1条光纤源成像结果时，设该二维光纤光谱图像的所有像素的个数为N，光谱图像的行数为M，S中的第t个像素表示为s_t，点扩展函数H在波长为λ处的第i个像素值表示为A_iλ，待抽取的一维光谱在波长为λ处的流量表示为C_λ，则S中所有像素值计算如下：

$S=\{s_1,s_2,...,s_t,...s_N\},s_t\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\underset{\mathrm{\λ}=1}{\overset{M}{\mathrm{\Σ}}}{A}_{\mathrm{i\λ}}{C}_{\mathrm{\λ}}+{\mathrm{\ϵ}}_t---\left(2\right)$

公式(2)用矩阵形式表达如下：

S_(N×1)＝A_(N×M)C_(M×N)+ε_(N×1)，  (3)

当该二维光纤光谱图像是由D条光纤源成像得到的结果时，设点扩展函数H在波长为λ处的第i个像素值表示为A′_iλ，抽取的一维光谱在波长为λ处的流量表示为C′_λ，则S中所有像素值计算如下：

$S=\{s_1,s_2,...s_t,...s_N\},s_t=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\underset{\mathrm{\λ}=1}{\overset{M\×D}{\mathrm{\Σ}}}{A}_{\mathrm{i\λ}}^{\′}{C}_{\mathrm{\λ}}^{\′}+{\mathrm{\ϵ}}_t---\left(4\right)$

公式(4)用矩阵形式表达如下：

S_(N×1)＝A′_(N×M×D)C′_(M×D×1)+ε_(N×1)  (5)

在公式(2)到公式(5)中，S是向量形式表示的二维光谱图像，A和A′是卷积矩阵，C和C′是待抽取的一维光谱，ε是光谱图像上的高斯白噪声，角标括号内符号表示矩阵的大小，符号N是所有像素的个数，M是光谱图像的行数，D是光谱图像中的被抽取光纤光谱源的个数；公式(5)则为所确定的用矩阵形式表示的二维光谱成像卷积方程。

4.根据权利要求1所述的基于GPU的快速反卷积二维光纤光谱抽取方法，其特征在于，步骤5)中所述采用共轭梯度最小二乘法反卷积求解二维光谱成像卷积方程过程中，对矩阵与向量相乘采用GPU进行加速的方法为：

在CUDA并行编程环境下实现矩阵与向量相乘的并行计算，需要对计算任务进行适当的划分，使之适合CUDA编程模型的特点，能充分利用GPU提供的资源，GPU上的任务分配方式如下：

每个线程负责计算输出向量的一个元素，即每个线程负责计算矩阵每行非零元素与向量乘积并求和，这种方式是串行算法到GPU上的直接映射；通过全局存储器进行通信，将一个线程块映射到矩阵每行元素的计算上，然后对同一个线程块得到的中间结果进行规约求和。