利用反卷积方法抽取多光纤光谱的一维谱的制作方法
【专利摘要】本发明公开了一种根据二维多光纤光谱图像形成原理，用反卷积方法来抽取可靠的一维光谱的方法。包括光纤光谱中心迹线的获得，每条光纤光谱上所有PSF的获得，目标函数的构造，噪声的抑制方法和目标函数的串行和并行求解算法。利用强的连续光谱获取光谱中心迹线，利用单根发射线光谱来获取所有波长处的PSF；根据噪声和光谱情况构造目标函数；给出了目标函数的串行和并行求解方法。本发明能纠正经典抽谱方法的失真，也解决反卷积方法由于存储量巨大而无法计算的问题，能得到更为真实的一维谱。同时该计算方法很灵活，内存用量少，扩展性强，容易用于并行计算，也是一个计算稀疏矩阵的新算法。
【专利说明】
利用反卷积方法抽取多光纤光谱的一维谱
技术领域
[0001] 本发明涉及多光纤光谱的一维谱抽取技术领域，和大规模稀疏矩阵串行和并行计 算领域。
【背景技术】
[0002] 在天文领域中，为了同时获取尽可能多的天体目标的光谱，都采用了多光纤光谱 技术。也就是，每一条光纤对准一个天体目标，多个天体的光被各自的光纤引导到光谱仪， 然后在CCD上色散成二维光谱。采用多光纤光谱技术，一张 CCD可以同时记录几百条光谱。 这些光谱是以二维谱图像的形式被CCD记录下来的，而用于数据分析的光谱一般是一维 谱，即横坐标是波长，纵坐标是流量。从而需要一个算法，把CCD记录的二维谱图像，转化成 可以用于分析的一维谱。
[0003] 目前流行的方法是把垂直于每一条二维光谱色散方向的每一行的流量按照一定 加权规则加起来，从而获得一维谱。这就要求加权规则在每一行甚至所有行都是固定不变 的，但这个假设是无法成立的。
[0004] 从另一方面看，C⑶记录的二维谱的图像，是进入光学系统之前的一维谱卷积上光 学系统的点源扩展函数（PSF)加上各种噪声得到的。若在二维谱图像上，每条光纤光谱上 所有位置处的PSF已知，则可以用二维谱图像反卷积掉已知的PSF来获得的一维谱。这种 抽谱方法更符合物理过程，从而更为可靠。
【发明内容】
[0005] 本发明的目的在于根据二维多光纤光谱图像形成原理，用反卷积方法来抽取更为 可靠的一维光谱，从而矫正当前抽谱方法的失真。
[0006] 为实现上述目的，本发明提出一种用反卷积方法从光纤光谱的二维CCD图像中抽 取一维光谱的方法。该方法包括：光纤光谱中心迹线的获得，每条光纤光谱所有波长处PSF 的获得，目标函数的构造，噪声的抑制方法、目标函数的串行和并行求解算法。
[0007] 所述的二维光谱中心迹线的获得方法是，从较强的平场定标灯的二维光谱，利用 重心法，求出每一条光纤光谱的每一行的重心位置，然后把每一行重心用低阶多项式拟合， 得到每一条光纤光谱的中心轨迹。
[0008] 所述的PSF获取方法是，对于波长定标灯的每一条光纤光谱图像，首先获取单根 发射线的离散轮廓。对于发射线比较弱的情况，可以延长曝光来提高该发射线的信噪比。然 后对离散的轮廓进行归一化，就获得了该发射线在它所处CCD位置处的离散的PSF。另外， 也可以通过激光梳来获得，只不过这种方法比较昂贵。如果想使用光滑的PSF，可以通过利 用B样条曲面插值离散的PSF来获得。我们把通过这种方法获得的PSF称为基本的PSF。
[0009] 所述的目标函数的构造方法是，根据图像的噪声类型来构造目标函数。
[0010] 高斯噪声情况：
[0011]
[0012]
[0013]
[0014] 其中，Nx表示图像总列数；Ny表示图像总行数；N f表示图像记录的光纤光谱总数； F(k，m)表示图像上位于第k行m列的CCD像素的计数；PSF1, Jk，m)表示位于第i条光纤 光谱在第j行处的PSF在CCD上第k行第m列像素上的计数；Cl, j表示第i条光纤光谱在 第.j行处的流量值，该值是要求解的值。
[0015] 高斯噪声下，目标函数的矩阵形式写法：
[0016]
C
[0017] 泊松噪声下，目标函数的矩阵形式写法：
[0018]
[0019] 其中，F是NxXNy行的列向量，F(k，m)是它的第kXN x+m个元素；A是ΝχΧΝ#， Nf X Ny列的矩阵，PSF L j (k，m)是它的第k X Nx+m行，第i X Ny+j列元素；C是一个Nf X Ny行的 列向量，Cy是它第iXNy+j个元素；W是一个NxXN yRNxXNy列的对角矩阵，它的第kXNx+m 行，kXNx+m列元素是F(k，m)，其它元素为0。
[0020] 所述的噪声抑制方法，是利用Tikhonov正则项来抑制噪声。可以使用所求变量的 0、1、2次等导数，及这些导数的组合来抑制噪声。矩阵形式的目标函数如下：
[0021] 高斯噪声下，目标函数为：
[0022]
[0023]
[0024]
[0025] 其中，Γ是Tikhonov矩阵，源自导数的离散化；α是权重对角矩阵；Tikhonov项 (rC) Ta (rc)对噪声有抑制作用。
[0026] 上面的公式的推导，都假设噪声是独立同分布的。若噪声情况复杂，则刻画噪声的 对称正定矩阵W会更复杂。无论何种情况，目标函数可以写为如下统一形式：
[0027]
⑴
[0028] 所述的目标函数的串行求解方法，是利用矩阵分块和迭代的方式求解目标函数 (1)。由于都是稀疏矩阵，矩阵在分块计算过程中，利用了其系数特性进行存储和计算。具 体的计算步骤如下：
[0029] 步骤0 :给定ε >〇,子矩阵块尺度k，和向量C的长度N。0 < k << N。
[0030] 步骤1 :设初始值C。: = 0, C中将要计算的变量块的初始位置s : = 0。对于目标 函数（1)，设 R : = F，res : = RTWR，RES := res，fy = 0〇
RES-ri，转到步骤3。
[0033] 步骤4 :若f\< ε，停止；否则f 1: = 〇,转到步骤2。
[0034] 所述的目标函数并行求解方法是，在计算过程中，把要求解的变量向量分成多个 变量块，每次并行计算那些没有交叉污染的变量块，直到算法收敛到最小点。具体步骤如 下：
[0035] 步骤0 :给定ε >〇,子矩阵块尺度k，和向量C的长度N。0 < k << N。
[0036] 步骤 1 :设全局变量的初始值 C。: = 0, R : = F，res : = RTWR，RES := res，f1: = 0〇 把变量c中按照变量编号顺序分为pv/^块。除了最后一块，每块有k个变量。记这些块 为= U，...，「7V/^，令p :=丨巧丨是这些矩阵块的集合。
[0037] 步骤2 :在集合P中选取变量块，这些变量块在矩阵A中对应的子矩阵的非0行号 两两之间都不相同。记这些变量块的集合为A。
[0038] 步骤3 :对于每个^ e Δ，用一个计算核计算。具体如下：记:^和方分别是矩阵A和 Γ中与G相对应列组成的子矩阵块。令
，
[0039] 步骤4少:=W Δ，若P ,为空集，转到步骤6 ;否则，转到步骤2。
[0040] 步骤5:重新把变量C进行划分成新的pV/&l变量块，使得每一个具有k个元 素块与上次块划分的所有具有k个的元素块至少有个变量不同。记这些块为 = ，令0:=·^}是这些矩阵块的集合。转到步骤2。
[0041] 步骤6 :若f\< ε，停止；否则f 1:= 〇,转到步骤5。
[0042] 本发明能够利用任意形状的PSF反卷积出任意大小的CXD光纤光谱图像。目标函 数的求解方法既可以用并行计算和也可以用串行计算方法进行计算。由于符合二维光纤光 谱形成的物理过程，本发明抽取的一维谱相对于传统的抽谱方法，其信噪比和分辨率会大 大提高，同时更加真实。本发明可以很自然地处理两根相邻的具有较严重相互污染的光纤 光谱图像，这要比传统方法更优越。本发明扩展性很强，可以通过设置矩阵W来刻画不同的 噪声，也可以通过修改Tikhonov矩阵α来根据实际情况更加灵活的对噪声进行抑制。
【附图说明】
[0043] 图1是本发明利用反卷积方法抽取多光纤光谱的一维谱的流程图。
[0044] 图2是本发明利用反卷积方法抽取多光纤光谱的一维谱的串行计算示意图。
[0045] 图3是本发明利用反卷积方法抽取多光纤光谱的一维谱的并行计算示意图。
【具体实施方式】
[0046] 为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照 附图，对本发明进一步详细说明。
[0047] 图1是本发明可利用反卷积方法抽取多光纤光谱的一维谱的流程图。具体包括： 光纤光谱中心迹线的获得，每条光纤光谱所有波长处PSF的获得，目标函数的构造，噪声的 抑制方法、目标函数的串行和并行求解算法。
[0048] 首先，我们需要知道每条光纤光谱的中心迹线，从而知道每条光纤光谱在CCD上 的位置。为此，需要有一个能发出较强连续谱的平场灯，例如碘钨灯，照射所有光纤的输入 端。灯光会通过光学系统在CCD上留下每条光纤的二维光谱图像。利用重心法，求出每一 条光纤光谱的每一行的重心位置，然后把每一行重心用低阶多项式拟合，得到每一条光纤 光谱的中心轨迹。
[0049] 其次，需要获取每条光纤光谱所有波长处的PSF。为此，需要有一个能发出多条单 根发射线的、且这些单根发射线尽可能均匀分布在整条光纤光谱图像中的灯。例如，汞镉 灯。用这种灯照射光纤输入端，光纤会通过光学系统在C⑶上留下每条光纤的灯谱图像，然 后我们可以沿着光谱中心迹线方向，根据一定的尺寸，依次提取每条光纤光谱上的所有单 根发射线。对于比较弱的单根合发射线，可以通过延长曝光时间来提高该发射线的信噪比。 然后我们对离散的发射线轮廓进行归一化，就获得了该发射线在它所处CCD位置处的离散 的PSF。另外，我们可以通过激光梳来获得，只不过这种方法比较昂贵。如果想使用光滑的 PSF，可以通过利用B样条曲面插值离散的PSF来获得把光滑的PSF。我们把通过这种方法 获得的PSF称为基本的PSF，这些基本的PSF在计算过程中存储在电脑内存中。对于每一条 光纤光谱图像，在没有灯谱发射线的位置，可以通过与它相邻的两个基本PSF线性插值获 得。这些非基本的PSF在需要的时候才计算，不存储在电脑内存中。
[0050] 然后，根据噪声模型设定目标函数中的W，根据光纤光谱的亮度和光纤光谱不同波 长处流量相对强度变化，设定目标函数中的α，从而确定了目标函数。
[0051] 再次，需要求解目标函数。图2是串行求解流程图，具体如下：
[0052] 步骤0 :给定ε > 〇,子矩阵块尺度k，和向量C的长度Ν。0 < k << Ν。
[0053] 步骤1 :设初始值C。: = 0, C中将要计算变量块的初始位置s : = 0。对于目标函 数（1)，设 R : = F，res : = RTWR，RES := res，A: = 0〇
[0054] 步骤2 :若s+k>N，则s : = N转到步骤3 ;否则，根据基本PSF计算出A的第 s+1，. . .，s+k列组成的子矩阵，并记为2，而方是Γ的第s+1，. . .，s+k组成的子矩阵， Zi :=(ATWA + BTaByl(A WR-BT〇TC0) , r, :=2RrWMs-X^WA^s-2C0^aBX s^XjBaBXx, Λ := i? -?，η := + , (^是 C 中的第 i 个元素，（cs+1，…，cs+k)T:= (c S+1，…，（^)4^:=5 + 1^/21 + 1，f1:= f Jr1, RES := RES-ri，转到步骤2。
[0055] 步骤3 :若s-k < 0，s : = 0,转到步骤4 ;否则，根据基本PSF计算出A的第 s-k+1，. . .，s列组成的子矩阵，并记为:i，而；§是Γ的第s-k+l，. . .，s组成的子矩阵， Xi := (A'WA +BTaBy1 (AWR-Bt〇TC0), r,-21^W^sWAXs-ICcrT raBXs -X/'B aBXs, -v4Xs，（^是 C 中的第 i 个元素 （c A!，· · ·，cs)T: = (c s k+1，…，cs)T+Xs， s:= 「A：/2~|- I，A: = f Jr1，RES : = RES-Ir1，转到步骤 3。
[0056] 步骤4 :若f\< ε，停止；否则f 1: = 〇,转到步骤2。
[0057] 最后，目标函数也可以用并行计算求解。图3是并行求解流程图，具体如下：
[0058] 步骤0 :给定ε > 〇,子矩阵块尺度k，和向量C的长度N。0 < k << N。
[0059] 步骤 1 :设全局变量的初始值 C。:= 0，R : = F，res : = RTWR，rd: = 〇，RES : = res， f1:=〇。把变量C中按照变量编号顺序分为「W/Αφ夬。除了最后一±夬，每块有k个变量。 记这些块为5,/ = 1,2,...,「7V/^，令$^丨是这些矩阵块的集合。
[0060] 步骤2 :在集合W中选取变量块，这些变量块在矩阵A中对应的子矩阵的非0行号 两两之间都不相同。记这些变量块的集合为A。
[0061] 步骤3 :对于每个^ e Δ，用一个计算核计算。具体如下：记]和互分别是矩阵A和 Γ中与G相对应列组成的子矩阵块。令 X5 ：= (I1WA +B1 aBy1 (ArWR-Br〇TC0) , rs :=2^WAXs -Xj~AW~AXS-IC^YaBXs -Xj~B aBXs, R-R-AXs, €:=€ + 尤，匕：=f^+bRES:= RES-rs。
[0062] 步骤4 #:= W Δ，若p为空集，转到步骤6 ;否则，转到步骤2。
[0063] 步骤5 :重新把变量C进行划分成新的「Λ7/Γ1 + 1变量块，使得每一个具有k个 元素块与上次块划分的所有具有k个的元素块至少有个变量不同。记这些块为 ^^_ = 1,2,.··,「Λ7Α:]，令是这些矩阵块的集合。转到步骤3。
[0064] 步骤6 :若f\< ε，停止；否则f 1: = 〇,转到步骤5。
[0065] 在本发明实施中，由于k值较小，所以在计算步骤中所有的矩阵块运算都可以进 行稀疏存储和计算，从而能加快计算速度。以上所述的具体实施例，对本发明的目的、技术 方案进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不 用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应 包含在本发明的保护范围之内。
【主权项】
1. 一种用反卷积方法从光纤光谱的二维CCD图像中抽取一维光谱的方法。其特点是： 根据已知的点源扩散函数（简记：PSF)，通过计算目标函数的最小值，来反卷积整个光纤光 谱图像，从而获得一维光谱。该方法包括：光纤光谱中心迹线的获得，每条光纤光谱所有波 长处PSF的获得，目标函数的构造，噪声的抑制方法，目标函数的串行和并行求解方法。2. 根据权利要求1所述的光纤光谱中心迹线的获得方法，其特征在于，从较强的平场 定标灯的二维光谱，利用重心法，求出每一条光纤光谱的每一行的重心位置，然后把每一行 重心用低阶多项式拟合，得到每一条光纤光谱的中心轨迹。3. 根据权利要求1所述的PSF获取方法，其特征在于，对于波长定标灯的每一条光纤光 谱图像，首先获取单根发射线的离散轮廓。对于发射线比较弱的情况，可以延长曝光来提高 该发射线的信噪比。然后对离散的轮廓进行归一化，就获得了该发射线在它所处CCD位置 处的离散的PSF。另外，也可以通过激光梳来获得，只不过这种方法比较昂贵。如果想使用 光滑的PSF，可以通过利用B样条曲面插值离散的PSF来获得。我们把通过这种方法获得的 PSF称为基本的PSF。 对于每一条光纤光谱图像，在没有灯谱发射线的位置，我们可以通过与它相邻的两个 基本PSF线性插值获得。 这些PSF都位于每一条光纤光谱的中心迹线上。4. 根据权利要求1所述的目标函数的构造方法，其特征在于，把图像的噪声分为高斯 噪声和泊松噪声两种情况来构造目标函数。 高斯噪声情况：其中，Nx表示图像总列数；Ny表示图像总行数；Nf表示图像记录的光纤光谱总数；F(k， m)表示图像上位于第k行m列的CCD像素的计数；PSFuG^m)表示位于第i条光纤光谱在 第j行处的PSF在CCD上第k行第m列像素上的计数； Cl, j表示第i条光纤光谱在第j行 处的流量值，该值是要求解的值。 高斯噪声下，目标函数的矩阵形式写法： - AC)T (F - AC) C 泊松噪声下，目标函数的矩阵形式写法： mjn(F - ACf W(F -AC) C ' 其中，F是NxXNy行的列向量，F(k，m)是它的第kXNx+m个元素；A是N xXNy行，NfXNy 列的矩阵，？3匕#，111)是它的第1^凡+111行，第1\%+」列元素；(：是一个队\%行的列向 量，h i是它第iXNy+j个元素；W是一个NxXNyR NxXNy列的对角矩阵，它的第kXNx+m行， kXNx+m列元素是F(k，m)，其它元素为0。5. 根据权利要求1所述的噪声抑制方法，其特征在于，利用Tikhonov正则项来抑制噪 声。可以使用所求变量的〇、1、2次等导数，及这些导数的组合来抑制噪声。矩阵形式的目 标函数如下： 高斯噪声下，目标函数为：其中，Γ是Tikhonov矩阵，源自导数的离散化；α是权重对角矩阵；Tikhonov项（rc) Ta (rc)对噪声有抑制作用。 上面的公式的推导，都假设噪声是独立同分布的。若噪声情况复杂，则刻画噪声的对称 正定矩阵W会更复杂。无论何种情况，目标函数可以写为如下统一形式：(1) 其中，Γ是Tikhonov矩阵，源自导数的离散化；α是权重对角矩阵；Tikhonov项（rc) Ta (rc)对噪声有抑制作用。6. 根据权利要求1所述的目标函数的串行和并行求解方法，其特征在于，利用矩阵分 块和迭代的方式求解目标函数（1)，由于都是稀疏矩阵，矩阵在分块计算过程中，利用了其 系数特性进行存储和计算。7. 根据权利要求6所述的目标函数的串行求解方法，具体的计算步骤如下： 步骤〇:给定ε >〇,子矩阵块尺度k，和向量C的长度N。0<k<<N。 步骤1 :设初始值C。: = 0, C中将要计算的变量块的初始位置S : = 0。对于目标函数 (1)，设 R : = F，res : = RTWR，RES : = res，= 0〇 步骤2 :若s+k > N，则s : = N转到步骤3 ;否则，记:i和云分别是矩阵A和Γ的第 s+1，…，s+k 列组成的子矩阵，尤:=(7#] + 矿6^)_1(7 呢-FarC0)，T?:=i? -IXi, 6 :=2Λ7诙IXi ，(^是 C 中的第 i 个元素， (cs+1，…，cs+k)T:= (c S+1，…，CsJkXvr=Hp/]"! + !，f1:= f Jr1, RES := RES-IT1，转到 步骤2。 步骤3 :若s-k < 0，s : = 0,转到步骤4 ;否则，记j和5分别是矩阵A和Γ的第 s-k+1，· · ·，s 列组成的子矩阵，Xs :=(7r^ +矿咖-7arC0)，i?:=i?-]iXs， η := - (2(^Τ?Χ5 +X/安《万尤），' (^是 C 中的第 i 个元素， (cs k+1，…，cs)T: = (c s k+1，…，cs) T+Xs，S := S -「A:/2"| -1，A: = f Jr1，RES : = RES-IT1，转到 步骤3。 步骤4 :若f\< ε，停止；否则f 1: = 〇,转到步骤2。8.根据权利要求6所述的目标函数的并行求解方法，其特征在于，在计算过程中，把要 求解的变量向量分成多个变量块，每次并行计算那些没有交叉污染的变量块，直到算法收 敛到最小点。具体计算步骤如下： 步骤〇:给定ε >〇,子矩阵块尺度k，和向量C的长度N。0<k<<N。 步骤 1 :设全局变量的初始值 C。: = 0, R : = F，res : = RTWR，RES : = res，f1: = 0。把 变量C中按照变量编号顺序分为块。除了最后一块，每块有k个变量。记这些块为 ，· = 1,2,…,「tv / q，令p :=后}是这些矩阵块的集合。 步骤2 :在集合炉中选取变量块，这些变量块在矩阵A中对应的子矩阵的非0行号两两 之间都不相同。记这些变量块的集合为A。 步骤3 :对于每个Cs e Δ，用一个计算核计算。具体如下：记 j和云分别是矩阵A和Γ中与$相对应列组成的子矩阵块。令 Xi ：= (a W^A+ 1ToB)-1 (AWR-BtaTC0) ， R-=R-AXs ， Q :=Q+XS ， rs :=2RrWAXs-X^~AW~AXS-2C0 rTTaBXs-XsTBTaBX s, f1：= f !+r,, RES ：= RES-rs〇 步骤4少:=<?ΛΔ，若p为空集，转到步骤6 ;否则，转到步骤2。 步骤5 :重新把变量C进行划分成新的变量块，使得每一个具有k个元素 块与上次块划分的所有具有k个的元素块至少有「&/2]-1个变量不同。记这些块为 ^，/ = l,2,...,「iV/q，令是这些矩阵块的集合。转到步骤2。 步骤6 :若f\< ε，停止；否则f 1: = 〇,转到步骤5。
【文档编号】G01J3/28GK105890757SQ201410765824
【公开日】2016年8月24日
【申请日】2014年12月15日
【发明人】李广伟, 张昊彤, 董义乔, 袁海龙, 雷亚娟, 白仲瑞, 杨卉沁
【申请人】中国科学院国家天文台