CN105890757B - 利用反卷积方法抽取多光纤光谱的一维谱 - Google Patents

Abstract

本发明公开了一种根据二维多光纤光谱图像形成原理，用反卷积方法来抽取可靠的一维光谱的方法。包括光纤光谱中心迹线的获得，每条光纤光谱上所有PSF的获得，目标函数的构造，噪声的抑制方法和目标函数的串行和并行求解算法。利用强的连续光谱获取光谱中心迹线，利用单根发射线光谱来获取所有波长处的PSF；根据噪声和光谱情况构造目标函数；给出了目标函数的串行和并行求解方法。本发明能纠正经典抽谱方法的失真，也解决反卷积方法由于存储量巨大而无法计算的问题，能得到更为真实的一维谱。同时该计算方法很灵活，内存用量少，扩展性强，容易用于并行计算，也是一个计算稀疏矩阵的新算法。

Description

利用反卷积方法抽取多光纤光谱的一维谱

技术领域

本发明涉及多光纤光谱的一维谱抽取技术领域，和大规模稀疏矩阵串行和并行计算领域。

背景技术

在天文领域中，为了同时获取尽可能多的天体目标的光谱，都采用了多光纤光谱技术。也就是，每一条光纤对准一个天体目标，多个天体的光被各自的光纤引导到光谱仪，然后在CCD上色散成二维光谱。采用多光纤光谱技术，一张CCD可以同时记录几百条光谱。这些光谱是以二维谱图像的形式被CCD记录下来的，而用于数据分析的光谱一般是一维谱，即横坐标是波长，纵坐标是流量。从而需要一个算法，把CCD记录的二维谱图像，转化成可以用于分析的一维谱。

目前流行的方法是把垂直于每一条二维光谱色散方向的每一行的流量按照一定加权规则加起来，从而获得一维谱。这就要求加权规则在每一行甚至所有行都是固定不变的，但这个假设是无法成立的。

从另一方面看，CCD记录的二维谱的图像，是进入光学系统之前的一维谱卷积上光学系统的点源扩展函数(PSF)加上各种噪声得到的。若在二维谱图像上，每条光纤光谱上所有位置处的PSF已知，则可以用二维谱图像反卷积掉已知的PSF来获得的一维谱。这种抽谱方法更符合物理过程，从而更为可靠。

发明内容

本发明的目的在于根据二维多光纤光谱图像形成原理，用反卷积方法来抽取更为可靠的一维光谱，从而矫正当前抽谱方法的失真。

为实现上述目的，本发明提出一种用反卷积方法从光纤光谱的二维CCD图像中抽取一维光谱的方法。该方法包括：光纤光谱中心迹线的获得，每条光纤光谱所有波长处PSF的获得，目标函数的构造，噪声的抑制方法、目标函数的串行和并行求解算法。

所述的二维光谱中心迹线的获得方法是，从较强的平场定标灯的二维光谱，利用重心法，求出每一条光纤光谱的每一行的重心位置，然后把每一行重心用低阶多项式拟合，得到每一条光纤光谱的中心轨迹。

所述的PSF获取方法是，对于波长定标灯的每一条光纤光谱图像，首先获取单根发射线的离散轮廓。对于发射线比较弱的情况，可以延长曝光来提高该发射线的信噪比。然后对离散的轮廓进行归一化，就获得了该发射线在它所处CCD位置处的离散的PSF。另外，也可以通过激光梳来获得，只不过这种方法比较昂贵。如果想使用光滑的PSF，可以通过利用B样条曲面插值离散的PSF来获得。我们把通过这种方法获得的PSF称为基本的PSF。

所述的目标函数的构造方法是，根据图像的噪声类型来构造目标函数。

高斯噪声情况：

泊松噪声情况：

其中，N_x表示图像总列数；N_y表示图像总行数；N_f表示图像记录的光纤光谱总数；F(k，m)表示图像上位于第k行m列的CCD像素的计数；PSF_i，j(k，m)表示位于第i条光纤光谱在第j行处的PSF在CCD上第k行第m列像素上的计数；c_i，j表示第i条光纤光谱在第.j行处的流量值，该值是要求解的值。

高斯噪声下，目标函数的矩阵形式写法：

泊松噪声下，目标函数的矩阵形式写法：

其中，F是N_x×N_y行的列向量，F(k，m)是它的第k×N_x+m个元素；A是N_x×N_y行，N_f×N_y列的矩阵，PSF_i，j(k，m)是它的第k×N_x+m行，第i×N_y+j列元素；C是一个N_f×N_y行的列向量，c_i，j是它第i×N_y+j个元素；W是一个N_x×N_y行N_x×N_y列的对角矩阵，它的第k×N_x+m行，k×N_x+m列元素是F(k，m)，其它元素为0。

所述的噪声抑制方法，是利用Tikhonov正则项来抑制噪声。可以使用所求变量的0、1、2次等导数，及这些导数的组合来抑制噪声。矩阵形式的目标函数如下：

高斯噪声下，目标函数为：

泊松噪声下，目标函数为：

其中，Γ是Tikhonov矩阵，源自导数的离散化；α是权重对角矩阵；Tikhonov项(ΓC)^Tα(ΓC)对噪声有抑制作用。

上面的公式的推导，都假设噪声是独立同分布的。若噪声情况复杂，则刻画噪声的对称正定矩阵W会更复杂。无论何种情况，目标函数可以写为如下统一形式：

所述的目标函数的串行求解方法，是利用矩阵分块和迭代的方式求解目标函数(1)。由于都是稀疏矩阵，矩阵在分块计算过程中，利用了其系数特性进行存储和计算。具体的计算步骤如下：

步骤0：给定ε>0，子矩阵块尺度k，和向量C的长度N。0＜k＜＜N。

步骤1：设初始值C₀：＝0，C中将要计算的变量块的初始位置s：＝0。对于目标函数(1)，设R：＝F，res：＝R^TWR，RES：＝res，f₁：＝0。

步骤2：若s+k>N，则s：＝N转到步骤3；否则，记和分别是矩阵A和Γ的第s+ 1，...，s+k列组成的子矩阵，c_i是C中的第i 个元素，(c_s+1，...，c_s+k)^T：＝(c_s+1，...，c_s+k)^T+X_s，f₁：＝f₁+r₁，RES：＝RES- r₁，转到步骤2。

步骤3：若s-k＜0，s：＝0，转到步骤4；否则，记和分别是矩阵A和Γ的第s-k+ 1，...，s列组成的子矩阵，c_i是C中的第i 个元素，(c_s-k+1，...，c_s)^T：＝(c_s-k+1，...，c_s)^T+X_s，f₁：＝f₁+r₁，RES：＝RES- r₁，转到步骤3。

步骤4：若f₁＜ε，停止；否则f₁：＝0，转到步骤2。

所述的目标函数并行求解方法是，在计算过程中，把要求解的变量向量分成多个变量块，每次并行计算那些没有交叉污染的变量块，直到算法收敛到最小点。具体步骤如下：

步骤0：给定ε>0，子矩阵块尺度k，和向量C的长度N。0＜k＜＜N。

步骤1：设全局变量的初始值C₀：＝0，R：＝F，res：＝R^TWR，RES：＝res，f₁：＝0。把变量C中按照变量编号顺序分为块。除了最后一块，每块有k个变量。记这些块为令是这些矩阵块的集合。

步骤2：在集合中选取变量块，这些变量块在矩阵A中对应的子矩阵的非0行号两两之间都不相同。记这些变量块的集合为Δ。

步骤3：对于每个用一个计算核计算。具体如下：记和分别是矩阵A和Γ中 与相对应列组成的子矩阵块。令 f₁：＝f₁+r_s，RES：＝RES-r_s。

步骤4：若为空集，转到步骤6；否则，转到步骤2。

步骤5：重新把变量C进行划分成新的变量块，使得每一个具有k个元素块与上次块划分的所有具有k个的元素块至少有个变量不同。记这些块为令是这些矩阵块的集合。转到步骤2。

步骤6：若f₁＜ε，停止；否则f₁：＝0，转到步骤5。

本发明能够利用任意形状的PSF反卷积出任意大小的CCD光纤光谱图像。目标函数的求解方法既可以用并行计算和也可以用串行计算方法进行计算。由于符合二维光纤光谱形成的物理过程，本发明抽取的一维谱相对于传统的抽谱方法，其信噪比和分辨率会大大提高，同时更加真实。本发明可以很自然地处理两根相邻的具有较严重相互污染的光纤光谱图像，这要比传统方法更优越。本发明扩展性很强，可以通过设置矩阵W来刻画不同的噪声，也可以通过修改Tikhonov矩阵α来根据实际情况更加灵活的对噪声进行抑制。

附图说明

图1是本发明利用反卷积方法抽取多光纤光谱的一维谱的流程图。

图2是本发明利用反卷积方法抽取多光纤光谱的一维谱的串行计算示意图。

图3是本发明利用反卷积方法抽取多光纤光谱的一维谱的并行计算示意图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明进一步详细说明。

图1是本发明可利用反卷积方法抽取多光纤光谱的一维谱的流程图。具体包括：光纤光谱中心迹线的获得，每条光纤光谱所有波长处PSF的获得，目标函数的构造，噪声的抑制方法、目标函数的串行和并行求解算法。

首先，我们需要知道每条光纤光谱的中心迹线，从而知道每条光纤光谱在CCD上的位置。为此，需要有一个能发出较强连续谱的平场灯，例如碘钨灯，照射所有光纤的输入端。灯光会通过光学系统在CCD上留下每条光纤的二维光谱图像。利用重心法，求出每一条光纤光谱的每一行的重心位置，然后把每一行重心用低阶多项式拟合，得到每一条光纤光谱的中心轨迹。

其次，需要获取每条光纤光谱所有波长处的PSF。为此，需要有一个能发出多条单根发射线的、且这些单根发射线尽可能均匀分布在整条光纤光谱图像中的灯。例如，汞镉灯。用这种灯照射光纤输入端，光纤会通过光学系统在CCD上留下每条光纤的灯谱图像，然后我们可以沿着光谱中心迹线方向，根据一定的尺寸，依次提取每条光纤光谱上的所有单根发射线。对于比较弱的单根合发射线，可以通过延长曝光时间来提高该发射线的信噪比。然后我们对离散的发射线轮廓进行归一化，就获得了该发射线在它所处CCD位置处的离散的PSF。另外，我们可以通过激光梳来获得，只不过这种方法比较昂贵。如果想使用光滑的PSF，可以通过利用B样条曲面插值离散的PSF来获得把光滑的PSF。我们把通过这种方法获得的PSF称为基本的PSF，这些基本的PSF在计算过程中存储在电脑内存中。对于每一条光纤光谱图像，在没有灯谱发射线的位置，可以通过与它相邻的两个基本PSF线性插值获得。这些非基本的PSF在需要的时候才计算，不存储在电脑内存中。

然后，根据噪声模型设定目标函数中的W，根据光纤光谱的亮度和光纤光谱不同波长处流量相对强度变化，设定目标函数中的α，从而确定了目标函数。

再次，需要求解目标函数。图2是串行求解流程图，具体如下：

步骤0：给定ε＞0，子矩阵块尺度k，和向量C的长度N。0＜k＜＜N。

步骤1：设初始值C₀：＝0，C中将要计算变量块的初始位置s：＝0。对于目标函数(1)，设R：＝F，res：＝R^TWR，RES：＝res，f₁：＝0。

步骤2：若s+k>N，则s：＝N转到步骤3；否则，根据基本PSF计算出A的第s+1，...，s+k 列组成的子矩阵，并记为而是Γ的第s+1，...，s+k组成的子矩阵， c_i是C中的第i 个元素，(c_s+1，...，c_s+k)^T：＝(c_s+1，...，c_s+k)^T+X_s，f₁：＝f₁+r₁，RES：＝RES- r₁，转到步骤2。

步骤3：若s-k＜0，s：＝0，转到步骤4；否则，根据基本PSF计算出A的第s-k+1，...，s 列组成的子矩阵，并记为而是Γ的第s-k+1，...，s组成的子矩阵，c_i是C中的第i个元素(c_s-k+1，...，c_s)^T：＝(c_s-k+1，...，c_s)^T+X_s， f₁：＝f₁+r₁，RES：＝RES-r₁，转到步骤3。

步骤4：若f₁＜ε，停止；否则f₁：＝0，转到步骤2。

最后，目标函数也可以用并行计算求解。图3是并行求解流程图，具体如下：

步骤0：给定ε＞0，子矩阵块尺度k，和向量C的长度N。0＜k＜＜N。

步骤1：设全局变量的初始值C₀：＝0，R：＝F，res：＝R^TWR，Γ_d：＝0，RES：＝res，f₁：＝0。把变量C中按照变量编号顺序分为块。除了最后一块，每块有k个变量。记这些块为令是这些矩阵块的集合。

步骤2：在集合中选取变量块，这些变量块在矩阵A中对应的子矩阵的非0行号两两之间都不相同。记这些变量块的集合为Δ。

步骤3：对于每个用一个计算核计算。具体如下：记和分别是矩阵A和Γ 中与相对应列组成的子矩阵块。令 f₁：＝f₁+r_s，RES：＝RES-r_s。

步骤4：若为空集，转到步骤6；否则，转到步骤2。

步骤5：重新把变量C进行划分成新的变量块，使得每一个具有k个元素块与上次块划分的所有具有k个的元素块至少有个变量不同。记这些块为令是这些矩阵块的集合。转到步骤3。

步骤6：若f₁＜ε，停止；否则f₁：＝0，转到步骤5。

在本发明实施中，由于k值较小，所以在计算步骤中所有的矩阵块运算都可以进行稀疏存储和计算，从而能加快计算速度。以上所述的具体实施例，对本发明的目的、技术方案进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (6)

1.一种用反卷积方法从光纤光谱的二维CCD图像中抽取一维光谱的方法，其特点是：根据已知的点源扩散函数PSF，通过计算目标函数的最小值，来反卷积整个光纤光谱图像，从而获得一维光谱，该方法包括：光纤光谱中心迹线的获得，每条光纤光谱所有波长处PSF的获得，目标函数的构造，噪声的抑制方法、目标函数的串行和并行求解方法；图像的噪声分为高斯噪声或者泊松噪声的情况来构造目标函数：

高斯噪声情况：

泊松噪声情况：

其中，N_x表示图像总列数；N_y表示图像总行数；N_f表示图像记录的光纤光谱总数；F(k，m)表示图像上位于第k行m列的CCD像素的计数；PSF_i，j(k，m)表示第i条光纤光谱在第.j行处的PSF在CCD上第k行第m列像素上的计数；c_i，j表示第i条光纤光谱在第j行处的流量值，该值是要求解的值；

高斯噪声下，目标函数的矩阵形式写法：

泊松噪声下，目标函数的矩阵形式写法：

其中，F是N_x×N_y行的列，F(k，m)是它的第k×N_x+m个元素；A是N_x×N_y行，N_f×N_y列的矩阵，PSF_i，j(k，m)是它的第k×N_x+m行，第i×N_y+j列元素；C是一个N_f×N_y行的列向量，c_i，j是它第i×N_y+j个元素；W是一个N_x×N_y行N_x×N_y列的对角矩阵，它的第k×N_x+m行，k×N_x+m列元素是F(k，m)，其它元素为0。

2.根据权利要求1所述的方法；PSF的获取方法如下，对于波长定标灯的每一条光纤光谱图像，首先获取单根发射线的离散轮廓，对于发射线较弱的情况，延长曝光来提高该发射线的信噪比，然后对离散的轮廓进行归一化，就获得了该发射线在它所处CCD位置处的离散的PSF，或者通过激光梳来获得，再通过B样条曲面插值上述PSF来获得光滑的PSF，我们把通过这种方法获得的PSF称为基本的PSF；

对于每一条光纤光谱图像，在没有灯谱发射线的位置，我们通过与它相邻的两个基本PSF线性插值获得；

这些PSF都位于同一条光纤光谱的中心迹线上。

3.根据权利要求1所述的方法，其特征在于，利用Tikhonov正则项来抑制噪声，使用所求变量的0、1、2次导数，及这些导数的组合来抑制噪声，矩阵形式的目标函数如下：

高斯噪声下，目标函数为：

泊松噪声下，目标函数为：

其中，Γ是Tikhonov矩阵，源自导数的离散化；α是权重对角矩阵；Tikhonov项(ΓC)^Tα(ΓC)对噪声有抑制作用；

上面的公式的推导，都假设噪声是独立同分布的，若噪声情况复杂，则刻画噪声的对称正定矩阵W会更复杂，无论何种情况，目标函数写为如下统一形式：

其中，Γ是Tikhonov矩阵，源自导数的离散化；α是权重对角矩阵；Tikhonov项(ΓC)^Tα(ΓC)对噪声有抑制作用。

4.根据权利要求3所述的方法，其特征在于，利用矩阵分块和迭代的方式求解目标函数(1)，由于都是稀疏矩阵，矩阵在分块计算过程中，利用了其系数特性进行存储和计算。

5.根据权利要求4所述的方法，具体的计算步骤如下：

步骤0：给定ε＞0，子矩阵块尺度k，和向量C的长度N，0＜k＜＜N；

步骤1：设初始值C₀：＝0，C中将要计算的变量块的初始位置s：＝0，对于目标函数(1)，设R：＝F，res：＝R^TWR，RES：＝res，f₁：＝0；

步骤2：若s+k＞N，则s：＝N转到步骤3；否则，记和分别是矩阵A和Γ的第s+1，...，s+k列组成的子矩阵，

c_i是C中的第i个元素，(c_s+1，...，c_s+k)^T：＝(c_s+1，...，c_s+k)^T+X_s，f₁：＝f₁+r₁，RES：＝RES-r₁，转到步骤2；

步骤3：若s-k＜0，s：＝0，转到步骤4；否则，记和分别是矩阵A和Γ的第s-k+1，...，s列组成的子矩阵，

c_i是C中的第i个元素，(c_s-k+1，...，c_s)^T：＝(c_s-k+1，...，c_s)^T+X_s，f₁：＝f₁+r₁，RES：＝RES-r₁，转到步骤3；

步骤4：若f₁＜ε，停止；否则f₁：＝0，转到步骤2。

6.根据权利要求4所述的方法，其特征在于，在计算过程中，把要求解的变量向量分成多个变量块，每次并行计算那些没有交叉污染的变量块，直到算法收敛到最小点，具体计算步骤如下：

步骤0：给定ε＞0，子矩阵块尺度k，和向量C的长度N，0＜k＜＜N；

步骤1：设全局变量的初始值C₀：＝0，R：＝F，res：＝R^TWR，RES：＝res，f₁：＝0，把变量C中按照变量编号顺序分为块，除了最后一块，每块有k个变量，记这些块为令是这些矩阵块的集合；

步骤2：在集合中选取变量块，这些变量块在矩阵A中对应的子矩阵的非0行号两两之间都不相同，记这些变量块的集合为Δ；

步骤3：对于每个用一个计算核计算，具体如下：记和分别是矩阵A和Γ中与相对应列组成的子矩阵块，令

f₁：＝f₁+r_s，RES：＝RES-r_s；

步骤4：若为空集，转到步骤6；否则，转到步骤2；

步骤5：重新把变量C进行划分成新的变量块，使得每一个具有k个元素块与上次块划分的所有具有k个的元素块至少有个变量不同，记这些块为令是这些矩阵块的集合，转到步骤2；

步骤6：若f₁＜ε，停止；否则f₁：＝0，转到步骤5。