CN102881003B - 一种消除ccd天文图像中宇宙射线的方法 - Google Patents

Abstract

本发明公开了一种消除CCD天文图像中宇宙射线的方法，包括以下步骤：（1）由原始图像I进行子采样得到子采样图像，对Laplacian算子进行改进，将子采样图像进行放大后与改进后的Laplacian算子进行卷积运算，然后进行去负值处理和恢复原始大小得到Laplacian图像L′；（3）识别Laplacian图像L′中的宇宙射线；（4）去除Laplacian图像L′中的宇宙射线；（5）利用分数阶微分的归整后的Tiansi算子对去除了宇宙射线的图像进行图像边缘增强处理。本发明方法大幅度地提高了宇宙射线的识别率，使星体中受到宇宙射线影响的部位能够最大程度的保留下来，提高了图像的处理质量。

Description

一种消除CCD天文图像中宇宙射线的方法

技术领域

本发明涉及一种CCD天文图像，特别涉及一种消除CCD天文图像中宇宙射线的方法。

背景技术

CCD(Charge-Coupled Device，电荷耦合元件)天文图像在采集过程中会受到各种噪声的影响，其中大能量的宇宙射线噪声严重影响到CCD天文图像的质量。宇宙射线指的是来自于宇宙中的一种具有相当大能量的带电粒子流，主要由质子、氦核、铁核等裸原子核组成，由于宇宙射线的分布没有规律，是随机分布的，有可能位于待测星的星象中，导致星象的位置和光亮度出现偏差，导致无法清晰正确的识别宇宙中的星体，因此有效的识别和消除宇宙射线噪声对于CCD天文图像的提取来说非常重要。

早期最简单的消除宇宙射线的方法为多图像处理算法，指的是多次获取同一区域的图像，因为一般情况下一个像素点或者一个区域的像素点只会在有限次内（一般都在3次以下）被宇宙射线影响，所以可以通过对同一区域获取多张图像，将那些好的图像的像素点代替被宇宙线污染的像素点。但是，在有些情况下获取同一目标的多幅图像时不可能的，甚至随着时间推移，拍摄环境的变化会对天光和目标光谱的位置和强度造成影响。

因此，许多研究人员提出了一种基于单幅图像消除宇宙射线的方法，例如Van Dokkum提出了一种基于Laplacian（拉普拉斯）边缘检测的宇宙射线去除算法，该算法是通过二阶偏微分的Laplacian算子检测星象的边缘以及宇宙射线射线陡峭的边缘，通过星象对称性的规律独立出宇宙射线，在无需对周围环境进行对比的情况下有效地探测并去除宇宙射线噪声。但仅仅利用该算法是不够的，因为简单的Laplacian算子会检测出不连续的边缘；当检测出宇宙射线只是通过简单的中值滤波代替，会对原图像造成一定的信息缺失，星象的有用信息会遭到衰减。随着成像技术的发展进步，获取的CCD天文图像的分辨率也越来越高，所以对于消除宇宙射线的算法要求也越来越高，对于其批量处理大型CCD天文图像来说，单CPU的串行实现效率相对来说也比较低。

发明内容

本发明的目的在于克服现有技术的缺点与不足，提供一种消除CCD天文图像中宇宙射线的方法，该方法大幅度地提高了宇宙射线的识别率，有效地去除了CCD天文图像中的宇宙射线，提高了图像处理的质量。

本发明的目的通过下述技术方案实现：一种消除CCD天文图像中宇宙射线的方法，包括以下步骤：

（1）由原始图像I进行子采样得到子采样图像；

（2）将传统的3×3矩阵模型的Laplacian算子改进为5×5矩阵模型的Laplacian算子，将子采样图像进行放大后与改进后的Laplacian算子进行卷积运算，将卷积运算后得到的Laplacian图像进行去负值处理，然后恢复到原始图像I的大小，得到与原始图像I的分辨率大小相同的Laplacian图像L′；

（3）识别Laplacian图像L′中的宇宙射线，包括以下步骤：

（3-1）先识别宇宙射线和大亮星体的星象：构造关联于原始图像I的噪声模型N，设定第一阈值σ_lim；

（3-2）将Laplacian图像L′的值与图像的像素点对应位置的噪声值作比值，得到像素点的信噪比S，S反映了Laplacian图像L′中每个像素含有噪声量的一个比值；

（3-3）对各像素点的信噪比S做中值滤波处理得到信噪比S′，判断像素点的信噪比S′的值是否大于第一阈值；若信噪比S′大于第一阈值，则该像素点为候选的宇宙射线；

（3-4）根据星象的对称性进一步识别宇宙射线和点光源星体的星象：构造关联于原始图像I的精细结构模型F，设定第二阈值f_lim；

（3-5）将Laplacian图像L′的值与精细结构模型F作比值，得到比值T；

（3-6）判断T值是否大于第二阈值；若T值大于第二阈值f_lim，则该像素点为候选的宇宙射线；

（3-7）将信噪比S′值大于第一阈值且T值大于第二阈值的像素点判定为宇宙射线；

（4）将在步骤（3）中识别出来的宇宙射线通过中值滤波消除，得到图像L^v，设定第三阈值，判断已识别出来的宇宙射线像素数量与原始图像I总像素数量的比值相对于第三阈值X的大小：

若已经识别出来的宇宙射线像素数量与原始图像I总像素数量的比值小于或等于第三阈值X，则步骤（1）中的原始图像I为L^v，重复迭代执行步骤（1）-（4），

若已经识别出来的宇宙射线像素数量与原始图像I总像素数量的比值大于第三阈值X，则执行步骤（5）；

（5）利用分数阶微分的Tiansi算子归整后得到的Tiansi′算子与步骤（5）中的图像L^v进行卷积运算，对L^v图像进行边缘增强处理，得到最终的图像。

优选的，所述步骤（1）原始图像I是通过CPU对其采集到的图像进行分割得到的，CPU将分割得到的多个原始图像I传送给CUDA中的GPU，由GPU对各分割得到的原始图像I进行所述步骤（1）到步骤（5）的并行处理；将经过所述步骤（5）处理得到的最终图像再由GPU传送给CPU，由CPU对各接收到的最终图像进行拼接得到完整的图像。

优选的，所述步骤（2）中子采样图像I′经过f_z放大后的图像I^fz为：

I^fz＝f_zI′；

I′为步骤（1）中原始图像I子采样后得到的子采样图像，f_z为子采样因子；

改进后的Laplacian算子的5x5矩阵模型为：

${\▿}^{2}f=\left(\begin{array}{ccccc}0& 1/8& 0& 1/8& 0\\ 1/8& 1/2& 1& 1/2& 1/8\\ 0& 1& 7& 1& 0\\ 1/8& 1/2& 1& 1/2& 1/8\\ 0& 1/8& 0& 1/8& 0\end{array}\right);$

所述步骤（2）中子采样图像进行放大后与改进后的Laplacian算子进行卷积运算得到的Laplacian图像L^fz为：

$L^{\mathrm{fz}}={\▿}^{2}f*I^{\mathrm{fz}};$

所述卷积运算后得到的图像L^fz进行去负值处理后得到的非负值的Laplacian图像L^fz′表达式为：

$L^{{\mathrm{fz}}^{\&prime;}}=\left\{\begin{array}{ccc}{L}^{\mathrm{fz}}& \mathrm{if}& L^{\mathrm{fz}}\&GreaterEqual;0\\ 0& \mathrm{if}& L^{\mathrm{fz}}<0\end{array}\right.;$

将非负值的Laplacian图像L^fz′的恢复到原始图像I的大小，其中原始图像I的分辨率大小为mxn，放大倍数为f_z，图像I^fz、Laplacian图像L^fz以及Laplacian图像L^fz′的分辨率大小均为f_zm×f_zn；得到分辨率大小为mxn的Lplacian图像L′的像素点L′_i，j表达式为：

$L_{i,j}^{\&prime;}=\left\{\begin{array}{cc}\frac{1}{{f_z}^2}\left(\begin{array}{c}(L_{f_{z}i,f_{z}j}^{{\mathrm{fz}}^{\&prime;}}+L_{f_{z}i-1,f_{z}j}^{{\mathrm{fz}}^{\&prime;}}+...+L_{f_{z}i-b,f_{z}j}^{{\mathrm{fz}}^{\&prime;}})+(L_{f_{z}i,f_{z}j-1}^{{\mathrm{fz}}^{\&prime;}}+L_{f_{z}i-1,f_{z}j-1}^{{\mathrm{fz}}^{\&prime;}}+...+L_{f_{z}i-b,f_{z}j-1}^{{\mathrm{fz}}^{\&prime;}})+\\ ...+(L_{f_{z}i,f_{z}j-b}^{{\mathrm{fz}}^{\&prime;}}+L_{f_{z}i-1,f_{z}j-b}^{{\mathrm{fz}}^{\&prime;}}+...+L_{f_{z}i-b,f_{z}j-b}^{{\mathrm{fz}}^{\&prime;}})\end{array}\right)& f_z\&GreaterEqual;3\\ \frac{1}{{f_z}^2}(L_{f_{z}i,f_{z}j}^{{\mathrm{fz}}^{\&prime;}}+L_{f_{z}i-1,f_{z}j}^{{\mathrm{fz}}^{\&prime;}}+L_{f_{z}i,f_{z}j-1}^{{\mathrm{fz}}^{\&prime;}}+L_{f_{z}i-1,f_{z}j-1}^{{\mathrm{fz}}^{\&prime;}})& f_z=2\end{array}\right.;$

其中i＝1,2,...,m，j=1,2,...,n，b＝f_z-1；

$L_{f_{z}i,f_{z}j}^{{\mathrm{fz}}^{\&prime;}},L_{f_{z}i-1,f_{z}j}^{{\mathrm{fz}}^{\&prime;}},L_{f_{z}i-b,f_{z}j}^{{\mathrm{fz}}^{\&prime;}},L_{f_{z}i,f_{z}j-1}^{{\mathrm{fz}}^{\&prime;}},L_{f_{z}i-1,f_{z}j-1}^{{\mathrm{fz}}^{\&prime;}},L_{f_{z}i-b,f_{z}j-1}^{{\mathrm{fz}}^{\&prime;}},L_{f_{z}i,f_{z}j-b}^{{\mathrm{fz}}^{\&prime;}},L_{f_{z}i-1,f_{z}j-b}^{{\mathrm{fz}}^{\&prime;}},L_{f_{z}i-b,f_{z}j-b}^{{\mathrm{fz}}^{\&prime;}}$

均为Laplacian图像L^fz′的像素点。

更进一步的，所述放大倍数f_z=2，图像I^fz、Laplacian图像L^fz以及Laplacian图像L^fz′的分辨率大小均为2m×2n；

得到分辨率大小为mxn的Lplacian图像L′的像素点L′_i，j表达式为：

$L_{i,j}^{\&prime;}=\frac{1}{4}(L_{2i-\mathrm{1,2}j-1}^{{\mathrm{fz}}^{\&prime;}}+L_{2i-\mathrm{1,2}j}^{{\mathrm{fz}}^{\&prime;}}+L_{2i,2j-1}^{{\mathrm{fz}}^{\&prime;}}+L_{2i,2j}^{{\mathrm{fz}}^{\&prime;}});i=\mathrm{1,2},...,m,j=\mathrm{1,2},...,n;$

其中均为Laplacian图像L^fz′的像素点。优选的，所述步骤（3）中的关联于原始图像I噪声模型N为：

$N=g^{-1}\sqrt{g(M_5*I)+{\mathrm{\&sigma;}}_{\mathrm{rn}}^2};$

其中g是增益因子；σ_rn是读入电子噪声，M₅*I表示对原始图像I进行5×5的中值滤波处理；

所述步骤（3）中信噪比S为：

$S=\frac{{\left(L^{\&prime;}\right)}^{+}}{f_{z}N};$

其中f_z为子采样因子，像素点的S比值越大，该像素点所包含的噪声越大；

所述步骤（3）中S经过中值滤波处理后得到的信噪比S′为：

S′=S-(S*M₅)；

S*M₅表示S通过5×5的中值滤波模板进行中值滤波处理。

更进一步的，所示噪声模型N中的增益因子g为7，读入电子噪声σ_rn为5。

优选的，所述步骤（3）中的关联于原始图像I的精细结构模型F为：

F=(M₃*I)-[(M₃*I)*M₇]；

其中M₃与M₇分别为3×3的中值滤波模板和7×7的中值滤波模板，原始图像I与M₃进行中值滤波操作M₃*I，得到原始图像中的中、低频信息，与M₇进行中值滤波操作(M₃*I)*M₇后得到图像中的低频信息；

所述步骤（3）中的Laplacian图像L′的值与精细结构模型F作比值得到的T值为：

$T=\frac{{\left(L^{\&prime;}\right)}^{+}}{F}.$

优选的，所述步骤（3）中的第一阈值σ_lim为0.5，第二阈值f_lim为1.5；所述步骤（4）中的第三阈值X为0.1%。

优选的，所述步骤(5)中的分数阶微分的Tiansi算子归整后得到Tiansi′算子与消除宇宙射线后的图像L^v进行卷积运算得到的增强后的图像L^T为：

L^T=Tiansi′*L^v；

其中Tiansi算子的表达式为：

Tiansi=8-8e+8×(e²-e)/2=8-12e+4e²=4(e-2)(e-1)；

e代表分数阶的阶次，其取值范围为(0,1)；所述Tiansi′算子为Tiansi算子的每一项均除以4(e-2)(e-1)。

更进一步的，所述分数阶的阶次e的大小根据图像增强需要的大小做调整。

本发明相对于现有技术具有如下的优点及效果：

（1）Van Dokkum提出的基于Laplacian的宇宙射线去除算法的最后是通过中值滤波来处理检测出的宇宙射线和噪声，中值滤波虽然能快速的消除宇宙射线和噪声，但会造成图像的信息丢失，星象的有用细节信息会受到影响。本发明通过引入分数阶微分归整后的Tiansi算子对去除了宇宙射线的图像进行图像边缘增强处理，使星象中受到宇宙射线影响的部位能够最大程度的保留下来，提高了图像的信息量，从而可得到更高质量的图像。另外分数阶微分除了能够一定程度地提升高频信号外，同时也能够很好地保留低频信号部分，在提取图像边缘方面更加精确，能更好的地保留图像的纹理信息，进一步提高图像处理的质量。

（2）本发明通过构造噪声模型和精细结构模型来区别宇宙射线和一个星体中的大亮星体和点光源星体，大幅度地提高了宇宙射线的识别率和算法收敛速度。另外改进后的Laplacian算子矩阵模型除了能够对边缘信息提取更加准确，同时也避免了一些假边缘噪声的出现，避免了Laplacian算子由于可能检查到不连续的边缘，而将星体误判为宇宙射线的情况。

（3）本发明可以采用由NVIDIA公司发布的基于GPU的并行运算平台CUDA进行了图像处理，CUDA是一个基于C语言的基础架构，具备有大量高性能计算指令和良好的编程接口，可以极大地提高图像处理的效率。另外CUDA将大型的CCD天文图像分割成GPU能承受的大型的子图像，使得各算法在CUDA中并行运行，极大提高了本发明图像处理的速度。

附图说明

图1是本发明在CUDA环境下进行图像处理的框图。

图2是本发明消除CCD天文图像中宇宙射线的方法流程图。

图3是本发明在GPU和CPU中图像像素和处理时间的关系图。

图4是各阶微分的调幅调相特性曲线图。

具体实施方式

下面结合实施例及附图对本发明作进一步详细的描述，但本发明的实施方式不限于此。

实施例

如图1和图2所示，本实施例的一种消除CCD天文图像中宇宙射线的方法，包括以下步骤：

（1）由运行在CPU上的主程序采集得到CCD天文图像，然后将原始图像分割成N个原始图像I，CPU将这N个原始图像I传送给CUDA中的GPU，由GPU对这N个原始图像I进行以下并行处理。

（2）将原始图像I进行子采样得到图像I′，对图像I′进行f_z倍放大得到图像I^fz，即I^fz＝f_zI′，其中本实施例的f_z为子采样因子，子采样因子f_z可以任意取值而不会影响处理结果，但算法的运行时间随着子采样因子的增大而呈现指数上升。因此，考虑到算法运行效率问题，本实施例的子采样因子f_z取值为2，其中原始图像I的分辨率为mxn。

（3）将传统只考虑4个方向到8个方向的3×3矩阵模型的Laplacian算子改进为16个方向的5×5矩阵模型的Laplacian算子，将放大后的子图像I^fz与改进后的Laplacian算子进行卷积运算得到Laplacian图像L^fz，再对Laplacian图像L^fz进行去负值处理，然后恢复到原始图像I的大小，最终得到与原始图像I的分辨率大小相同的Laplacian图像L′。

本实施例改进后的Laplacian算子的计算公式：

${\&dtri;}^{2}f=7f(x,y)-[f(x,y-1)+f(x+1,y)+f(x-1,y)+f(x,y+1)]$

$-\frac{1}{2}[f(x-1,y-1)+f(x+1,y-1)+f(x-1,y+1)+f(x+1,y+1)]$

$-\frac{1}{8}[f(x+1,y-2)+f(x-1,y-2)+f(x+2,y+1)+f(x-2,y+1)$

$+(x+2,y-1)+f(x-2,y-1)+f(x+1,y+2)+f(x-1,y+2)]$

改进后的Laplacian算子的5x5矩阵模型：

${\&dtri;}^{2}f=\left(\begin{array}{ccccc}0& 1/8& 0& 1/8& 0\\ 1/8& 1/2& 1& 1/2& 1/8\\ 0& 1& 7& 1& 0\\ 1/8& 1/2& 1& 1/2& 1/8\\ 0& 1/8& 0& 1/8& 0\end{array}\right);$

放大后的子图像I^fz与改进后的Laplacian算子进行卷积运算得到Laplacian图像L^fz为：

$L^{\mathrm{fz}}={\&dtri;}^{2}f*I^{\mathrm{fz}};$

使用改进后的Laplacian算子矩阵模型除了能够对边缘信息提取更加准确，同时也避免了一些假边缘噪声的出现。

经过Laplacian算子卷积运算后的Laplacian图像L^fz中，星象或者宇宙射线边缘的外部会呈负值，而对于宇宙射线和星体边缘和内部的像素点会呈正值。然而由于后续的处理中，负值部分会被误认为是图像的边缘成份或者宇宙射线，即图像边缘相邻的外部像素误认为是宇宙射线，从而影响图像的有用信息。因此为了防止宇宙射线对相邻外部像素的影响，需要对Laplacian图像L^fz中的负值取值为0，从而得到非负值的Laplacian图像L^fz，对Laplacian图像L^fz进行去负值处理后得到的非负值的Laplacian图像L^fz′为：

$L^{{\mathrm{fz}}^{\&prime;}}=\left\{\begin{array}{ccc}{L}^{\mathrm{fz}}& \mathrm{if}& L^{\mathrm{fz}}\&GreaterEqual;0\\ 0& \mathrm{if}& L^{\mathrm{fz}}<0\end{array}\right.;$

将非负值的Laplacian图像L^fz′的恢复到原始图像I的大小，其中原始图像I的分辨率大小为mxn，放大倍数为f_z，所以图像I^fz、Laplacian图像L^fz以及Laplacian图像L^fz′的分辨率大小均为f_zm×f_zn；得到分辨率大小为mxn的Lplacian图像L′各像素点L′_i，j的表达式为：

$L_{i,j}^{\&prime;}=\left\{\begin{array}{cc}\frac{1}{{f_z}^2}\left(\begin{array}{c}(L_{f_{z}i,f_{z}j}^{{\mathrm{fz}}^{\&prime;}}+L_{f_{z}i-1,f_{z}j}^{{\mathrm{fz}}^{\&prime;}}+...+L_{f_{z}i-b,f_{z}j}^{{\mathrm{fz}}^{\&prime;}})+(L_{f_{z}i,f_{z}j-1}^{{\mathrm{fz}}^{\&prime;}}+L_{f_{z}i-1,f_{z}j-1}^{{\mathrm{fz}}^{\&prime;}}+...+L_{f_{z}i-b,f_{z}j-1}^{{\mathrm{fz}}^{\&prime;}})+\\ ...+(L_{f_{z}i,f_{z}j-b}^{{\mathrm{fz}}^{\&prime;}}+L_{f_{z}i-1,f_{z}j-b}^{{\mathrm{fz}}^{\&prime;}}+...+L_{f_{z}i-b,f_{z}j-b}^{{\mathrm{fz}}^{\&prime;}})\end{array}\right)& f_z\&GreaterEqual;3\\ \frac{1}{{f_z}^2}(L_{f_{z}i,f_{z}j}^{{\mathrm{fz}}^{\&prime;}}+L_{f_{z}i-1,f_{z}j}^{{\mathrm{fz}}^{\&prime;}}+L_{f_{z}i,f_{z}j-1}^{{\mathrm{fz}}^{\&prime;}}+L_{f_{z}i-1,f_{z}j-1}^{{\mathrm{fz}}^{\&prime;}})& f_z=2\end{array}\right.;$

其中i＝1,2,...,m，j=1,2,...,n，b＝f_z-1；

$L_{f_{z}i,f_{z}j}^{{\mathrm{fz}}^{\&prime;}},L_{f_{z}i-1,f_{z}j}^{{\mathrm{fz}}^{\&prime;}},L_{f_{z}i-b,f_{z}j}^{{\mathrm{fz}}^{\&prime;}},L_{f_{z}i,f_{z}j-1}^{{\mathrm{fz}}^{\&prime;}},L_{f_{z}i-1,f_{z}j-1}^{{\mathrm{fz}}^{\&prime;}},L_{f_{z}i-b,f_{z}j-1}^{{\mathrm{fz}}^{\&prime;}},L_{f_{z}i,f_{z}j-b}^{{\mathrm{fz}}^{\&prime;}},L_{f_{z}i-1,f_{z}j-b}^{{\mathrm{fz}}^{\&prime;}},L_{f_{z}i-b,f_{z}j-b}^{{\mathrm{fz}}^{\&prime;}}$

均为Laplacian图像L^fz′的像素点。

由于本实施例放大倍数f_z=2，所以图像I^fz、Laplacian图像L^fz以及Laplacian图像L^fz′的分辨率大小均为2m×2n；

得到分辨率大小为mxn的Lplacian图像L′各像素点L′_i,j的表达式为：

$L_{i,j}^{\&prime;}=\frac{1}{4}(L_{2i-\mathrm{1,2}j-1}^{{\mathrm{fz}}^{\&prime;}}+L_{2i-\mathrm{1,2}j}^{{\mathrm{fz}}^{\&prime;}}+L_{2i,2j-1}^{{\mathrm{fz}}^{\&prime;}}+L_{2i,2j}^{{\mathrm{fz}}^{\&prime;}});i=\mathrm{1,2},...,m,j=\mathrm{1,2},...,n;$

其中均为Laplacian图像L^fz′的像素点。

（4）识别Laplacian图像L′中的宇宙射线：

（4-1）首先识别宇宙射线和大亮星体（高亮度成像大的星体）的星象，构造一个关联于原始图像I的噪声模型N，设定第一阈值σ_lim，利用Laplacian图像L′的值与图像的像素点对应位置的噪声值作比值，得到各像素点的信噪比S，再通过设定一个阈值σ_lim进行比较，识别出宇宙射线与大亮星体。

其中噪声模型N为：

$N=g^{-1}\sqrt{g(M_5*I)+{\mathrm{\&sigma;}}_{\mathrm{rn}}^2};$

其中，g是增益因子，其单位为电子ADU^-1，本实施例g=7；σ_rn是读入电子噪声，本实施例σ_rn=5；M₅*I是对原始图像I进行5×5的中值滤波处理。

根据噪声模型N和得到的Laplacian图像L′，即可得到信噪比S为：

$S=\frac{{\left(L^{\&prime;}\right)}^{+}}{f_{z}N};$

其中，f_z为子采样因子，公式中S反映了Laplacian图像L′中每个像素点含有噪声量的一个比值，像素点的噪声越大，其比值也越大，信噪比S值大于第一阈值σ_lim的像素点推测为候选的宇宙射线。

由于经Laplacian算子卷积后的得到的Laplacian图像提高了一个像素，特别是高频像素的灰度，使得Laplacian图像的像素与其相邻域的像素的差别增大(对比于原始图像)。但是Laplacian图像的所示的边缘及宇宙射线(中高频像素)并没有包含边缘的特征信息。由于泊松噪声所成的星象，其边缘均匀平滑，其边缘信噪比值S一般较小，一般不会超过σ_lim，将其误判为宇宙射线的概率比较小。但是对于被扩展的较大的大亮星体，特别是当子采样函数无法很好采样时，其误判为宇宙射线的概率比较大。解决此问题需要通过较大范围的中值滤波，去除S中被扩展的较大的大亮星体，如下公式所示：

S′=S-(S*M₅)；

通过5×5的中值滤波模板，对S进行中值滤波方法S*M₅，得到S中的扩展的较大的大亮星体，由于5×5的中值滤波足够大，使得宇宙射线和噪声不包含在其中。最后得到的S'有效的去除的被扩展的大亮星体，保留了宇宙射线和没被扩展的亮星体，（没被扩展的亮星体指的是亮度高成像小的星体，也称为点光源星体），使得信噪比为S′的图像中对应的像素中宇宙射线的比值增加。该公式中S′反映了Laplacian图像L′中每个像素点含有噪声量的一个比值，像素点噪声越大，S′值也就越大；因此判断像素点的信噪比S′的值是否大于第一阈值，以得到更为精确的宇宙射线的候选像素点；若信噪比S′大于第一阈值，则该像素点为候选的宇宙射线；

（4-2）利用星象的对称性来进一步区别宇宙射线与点光源星体的星象：构造一个关联于原始图像I的精细结构模型F，设定第二阈值f_lim。

F=(M₃*I)-[(M₃*I)*M₇]；

其中，M₃与M₇分别为3×3的中值滤波模板和7×7的中值滤波模板。原始图像I与M₃进行中值滤波操作M₃*I后，得到的是图像中的中、低频信息；原始图像I与M₇进行中值滤波操作(M₃*I)*M₇，得到是图像中的低频信息。两者相减后，得到的精细结构图像F剩下的只有中频信息。

构造关联于Laplacian图像L′比值T，具体为：

$T=\frac{{\left(L^{\&prime;}\right)}^{+}}{F};$

在3×3的中值滤波下，点光源星体虽然被周围相邻像素的中值所替代，但由于其具有对称性的性质，其相邻像素的中值依然保持高灰度值，因此点光源星体所在区域的F值会比较大。宇宙射线尺寸一般小于3×3，在进行中值滤波时，其灰度值会被周围相邻的像素削弱，所以宇宙射线所在区域的F值会比较小。根据上述T的表达式可以得出点光源星体像素点的T值比较小，而宇宙射线像素点的T值比较大，根据设定的第二阈值f_lim区别点光源星体与宇宙射线，把T值大于f_lim的像素判定为候选的宇宙射线。

（4-3）当像素点满足S′>σ_lim且T>f_lim时，可以判断该像素点为宇宙射线，其中σ_lim=0.5和f_lim=1.5。

（5）将在步骤（4）中识别出来的宇宙射线通过Van Dokkum提出的Laplacian边缘检测的宇宙射线去除算法的中值滤波进行消除处理，得到图像L^v，采用的中值滤波的模版为5×5；然后判断已经识别出来的宇宙射线像素数量与图像总像素数量的比值是否小于或等于第三阈值X值，其中本实施例的第三阈值X为0.1%；

若已经识别出来的宇宙射线像素数量与图像总像素数量的比值小于或等于第三阈值X，则步骤（2）中的原始图像I为L^v，重复迭代执行步骤（2）-（5）；

若已经识别出来的宇宙射线像素数量与图像总像素数量的比值大于第三阈值X，则执行步骤（6）；

（6）利用分数阶微分的Tiansi算子（掩摸算子）归整（归一化）后的Tiansi′算子对去除宇宙射线的图像L^v进行边缘增强处理。分数阶微分除了能够一定程度地提升高频信号外同时也能够很好地保留低频信号部分，分数阶微分在提取的图像边缘方面更加精确，能很好地保留图像的纹理信息。

其中基于分数阶微分的Tiansi算子如表1所示：

表1

(e2-e)/2	0	(e2-e)/2	0	(e2-e)/2
					0	-e	-e	-e	0
(e2-e)/2	-e	8	-e	(e2-e)/2
					0	-e	-e	-e	0
(e2-e)/2	0	(e2-e)/2	0	(e2-e)/2

e代表分数阶的阶数，其取值范围为(0,1)。对分数阶微分Tiansi算子作进一步归整处理。由于Tiansi算子与图像进行卷积时，会使得每个像素都得到4(e-2)(e-1)的加权值，如以下公式所示:

Tiansi=8-8e+8×(e²-e)/2=8-12e+4e²=4(e-2)(e-1)；

因此对Tiansi算子的每一项均除以4(e-2)(e-1)，得到归整后的Tiansi算子Tiansi′算子，如表格2所示。其中e代表分数阶的阶次，其大小根据图像增强需要的大小做调整。

表2

e/8(e-2)	0	e/8(e-2)	0	e/8(e  2)
					0	-e/4(e-2)(e-1)	-e/4(e-2)(e-1)	-e/4(e-2)(e-1)	0
e/8(e-2)	-e/4(e-2)(e-1)	2/(e-2)(e-1)	-e/4(e-2)(e-1)	e/8(e-2)
					0	-e/4(e-2)(e-1)	-e/4(e-2)(e-1)	-e/4(e-2)(e-1)	0
e/8(e-2)	0	e/8(e-2)	0	e/8(e-2)

将分数阶微分的归整后得到的Tiansi′算子与消除宇宙射线的图像L^v进行卷积运算，得到最终图像L^T：

L^T＝Tiansi′*L^v；

图像L^v经过与Tiansi′算子卷积运算后，大部分的低频信息部分（平滑的背景和无纹理的图像内部）的像素基本不变或者变化很小，而高、中频信息(图形图像的边缘和纹理细节部分)部分的像素值会得到提升，输出的像素会比原来的像素值大。因此，经过与Tiansi′算子卷积后的图像L^T的高频部分(图像边缘)和中频部分(图像的纹理)特征得到突显，线条更加清晰；而图像低频部分(平滑的背景和无纹理的图像内部)会基本保持不表，从而完成了对图像的增强处理，减少由于算法中的中值滤波处理或算法误判的宇宙射线产生的图像信息丢失。

（7）GPU将经过步骤（6）处理过的N子图像L^T传送给CPU，CPU将接受到的N个子图像进行拼接处理得到完整的图像。

通过本实施例步骤（3）中得出的最终Lplacian图像L′为：

$L^{\&prime;}\&ap;\left\{\begin{array}{cc}{f}_z(I-b)& I-b\&GreaterEqual;0\\ 0& I-b<0\end{array}\right.;$

其中b代表图像的平滑背景(即为图像中的低频部分)，I-b代表图像中的中高频信号部分，即星象边缘、泊松噪声及宇宙射线；f_z为子采样因子。

经过子采样的图像I与经过该进后的Laplacian算子卷积后得到的Laplacian图像L′为f_z(I-b)，是中高频像素的图像，包括星象边缘、泊松噪声及宇宙射线。其中当I-b出现小于0的情况时，该像素设为0值，从而可以得到非负值的Laplacian图像L′。

本实施例步骤（5）中得出的消除宇宙射线后的图像L^v表示为：

$L^v\&ap;\left\{\begin{array}{cc}I*M_5& \mathrm{if}{S}^{\&prime;}>{\mathrm{\&sigma;}}_{\lim }\mathrm{andT}>f_{\lim }\\ I& \mathrm{else}\end{array}\right.;$

其中M₅为5×5的中值滤波模板。即当像素点满足S′>σ_lim且T>f_lim时，该像素点为宇宙射线，通过中值滤波消除后得到图像L^v，当像素点不满足S′>σ_lim且T>f_lim时，该像素点不是宇宙射线，将不会被消除，所以得到的图像L^v基本上和原始图像I一致。

为提高本实施例的图像的处理速度，本实施例采用了由NVIDIA公司发布的基于GPU的并行运算平台CUDA(Compute Unified Device Architecture)来对本实施例的图像进行处理。CUDA是一个基于C语言的基础架构，具备有大量高性能计算指令和良好的编程接口，可以极大地提高本实施例图像处理的效率。将大型的CCD图像分割成GPU能承受的大小的子图像，然后由CUDA中的GPU(GraphicProcessing Unit，图像处理单元)对图像进行并行处理。如图3所示随着图像像素的提高，仅采用CPU处理图像所需要的时间越来越长，而采用GPU处理图像的时间几乎都没有变化。

在本实施例中如果直接利用改进后的Laplacian算子与子采样图像I′进行卷积，在对高亮像素(宇宙射线或者高亮星象)周边的像素进行卷积运算时，该像素会出现负值，然后会在后续的处理中被置为0，导致该像素被衰减。因此，在对图像进行卷积运算之前，必须对图像进行子采样放大，将放大倍数称为子采样因子f_z。然而子采样因子f_z可以任意取值，不会影响实验结果，但算法的运行时间随着子采样因子的增大而呈现指数上升。考虑到算法运行效率问题，本实施例的子采样因子f_z取值为2。

如图4所示为各阶微分的调幅调相特性曲线，从图中可以得出：各阶微分运算都能够提升高频和中频信号(|v|>1)，削弱低频信号(|v|<1)。对于高、中频信号的提升来说，高、中频频信号的提升幅度与微分运算的阶数成正比，即随着微分运算阶数的增大，高、中频信号的提升幅度成非线性的增长。对于低频信号来说，低频信号的削弱幅度与微分运算阶数成正比，即随着微分运算阶数的增大，低频信号的削弱幅度呈现非线性增长。二阶微分运算对高频信号的提升幅度远远大于一阶微分运算。同时，对于低频信号的削弱程度也远远大于一阶微分。然而对于e<1，即分数阶微分来说，对于高频信号也有足够大地提升，同时与整数阶一样，能够一定程度地提升中频信号，而对低频信号削弱程度却远远低于整数阶，即能有效地保留低频信号部分。总之，对于整数阶微分运算来说，能极大提升高频信号的同时，却无法很好地保留低频信号。而对于分数阶微分来说，能够一定程度地提升高频信号的同时也能很好地保留低频信号部分。

上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。

Claims (9)

1.一种消除CCD天文图像中宇宙射线的方法，其特征在于，包括以下步骤：

(1)由原始图像I进行子采样得到子采样图像；

(2)将传统的3×3矩阵模型的Laplacian算子改进为5×5矩阵模型的Laplacian算子，将子采样图像进行放大后与改进后的Laplacian算子进行卷积运算，将卷积运算后得到的Laplacian图像进行去负值处理，然后恢复到原始图像I的大小，得到与原始图像I的分辨率大小相同的Laplacian图像L′；

(3)识别Laplacian图像L′中的宇宙射线，包括以下步骤：

(3-1)先识别宇宙射线和大亮星体的星象：构造关联于原始图像I的噪声模型N，设定第一阈值σ_lim；

(3-2)将Laplacian图像L′的值与图像的像素点对应位置的噪声值作比值，得到像素点的信噪比S，S反映了Laplacian图像L′中每个像素含有噪声量的一个比值；

(3-3)对各像素点的信噪比S做中值滤波处理得到信噪比S′，判断像素点的信噪比S′的值是否大于第一阈值；若信噪比S′大于第一阈值，则该像素点为候选的宇宙射线；

(3-4)根据星象的对称性进一步识别宇宙射线和点光源星体的星象：构造关联于原始图像I的精细结构模型F，设定第二阈值f_lim；

(3-5)将Laplacian图像L′的值与精细结构模型F作比值，得到比值T；

(3-6)判断T值是否大于第二阈值；若T值大于第二阈值f_lim，则该像素点为候选的宇宙射线；

(3-7)将信噪比S′值大于第一阈值且T值大于第二阈值的像素点判定为宇宙射线；

(4)将在步骤(3)中识别出来的宇宙射线通过中值滤波消除，得到图像L^v，设定第三阈值，判断已识别出来的宇宙射线像素数量与原始图像I总像素数量的比值相对于第三阈值X的大小：

若已经识别出来的宇宙射线像素数量与原始图像I总像素数量的比值小于或等于第三阈值X，则步骤(1)中的原始图像I为L^v，重复迭代执行步骤(1)-(4)，

若已经识别出来的宇宙射线像素数量与原始图像I总像素数量的比值大于第三阈值X，则执行步骤(5)；

(5)利用分数阶微分的Tiansi算子归整后得到的Tiansi′算子与步骤(5)中的图像L^v进行卷积运算，对L^v图像进行边缘增强处理，得到最终的图像；

所述步骤(2)中子采样图像I′经过f_z放大后的图像I^fz为：

I^fz＝f_zI′；

I′为步骤(1)中原始图像I子采样后得到的子采样图像，f_z为子采样因子；

改进后的Laplacian算子的5x5矩阵模型为：

${\&dtri;}^{2}f=\left(\begin{array}{ccccc}0& 1/8& 0& 1/8& 0\\ 1/8& 1/2& 1& 1/2& 1/8\\ 0& 1& 7& 1& 0\\ 1/8& 1/2& 1& 1/2& 1/8\\ 0& 1/8& 0& 1/8& 0\end{array}\right);$

所述步骤(2)中子采样图像进行放大后与改进后的Laplacian算子进行卷积运算得到的Laplacian图像L^fz为：

$L^{\mathrm{fz}}={\&dtri;}^{2}f*I^{\mathrm{fz}};$

所述卷积运算后得到的图像L^fz进行去负值处理后得到的非负值的Laplacian图像L^fz′表达式为：

$L^{{\mathrm{fz}}^{\&prime;}}=\left\{\begin{array}{ccc}{L}^{\mathrm{fz}}& \mathrm{if}& L^{\mathrm{fz}}\&GreaterEqual;0\\ 0& \mathrm{if}& L^{\mathrm{fz}}<0\end{array}\right.$

将非负值的Laplacian图像L^fz′的恢复到原始图像I的大小，其中原始图像I的分辨率大小为m×n，放大倍数为f_z，图像I^fz、Laplacian图像L^fz以及Laplacian图像L^fz′的分辨率大小均为f_zm×f_zn；得到分辨率大小为m×n的Lplacian图像L′的像素点L′_i,j表达式为：

$L_{i,j}^{\&prime;}=\left\{\begin{array}{cc}\frac{1}{{f_z}^2}\left(\begin{array}{c}(L_{f_{z}i,f_{z}j}^{{\mathrm{fz}}^{\&prime;}}+L_{f_{z}i-1,f_{z}j}^{f^{z^{\&prime;}}}+...+L_{f_{z}i-b,f_{z}j}^{{\mathrm{fz}}^{\&prime;}})+(L_{f_{z}i,f_{z}j-1}^{{\mathrm{fz}}^{\&prime;}}+L_{f_{z}i-1,f_{z}j-1}^{f^{z^{\&prime;}}}+...+L_{f_{z}i-b,f_{z}j-1}^{{\mathrm{fz}}^{\&prime;}})+\\ ...+(L_{f_{z}i,f_{z}j-b}^{{\mathrm{fz}}^{\&prime;}}+L_{f_{z}i-1,f_{z}j-b}^{f^{z^{\&prime;}}}+...+L_{f_{z}i-b,f_{z}j-b}^{{\mathrm{fz}}^{\&prime;}})\end{array}\right)& f_z\&GreaterEqual;3\\ \frac{1}{{f_z}^2}(L_{f_{z}i,f_{z}j}^{{\mathrm{fz}}^{\&prime;}}+L_{f_{z}i-1,f_{z}j}^{f^{z^{\&prime;}}}+...+L_{f_{z}i,f_{z}j-1}^{{\mathrm{fz}}^{\&prime;}}+L_{f_{z}i-1,f_{z}j-1}^{{\mathrm{fz}}^{\&prime;}})& f_z=2\end{array}\right\};$

其中i＝1,2,...,m，j＝1,2,...,n，b＝f_z-1；

$L_{f_{z}i,f_{z}j}^{{\mathrm{fz}}^{\&prime;}},L_{f_{z}i-1,f_{z}j}^{{\mathrm{fz}}^{\&prime;}},L_{f_{z}i-b,f_{z}j}^{{\mathrm{fz}}^{\&prime;}},L_{f_{z}i,f_{z}j-1}^{{\mathrm{fz}}^{\&prime;}},L_{f_{z}i-1,f_{z}j-1}^{{\mathrm{fz}}^{\&prime;}},L_{f_{z}i-b,f_{z}j-1}^{{\mathrm{fz}}^{\&prime;}},L_{f_{z}i,f_{z}j-b}^{{\mathrm{fz}}^{\&prime;}},L_{f_{z}i-1,f_{z}j-b}^{{\mathrm{fz}}^{\&prime;}},L_{f_{z}i-b,f_{z}j-b}^{{\mathrm{fz}}^{\&prime;}}$ 均为Laplacian图像L^fz′的像素点。

2.根据权利要求1所述的消除CCD天文图像中宇宙射线的方法，其特征在于，所述步骤(1)原始图像I是通过CPU对其采集到的图像进行分割得到的，CPU将分割得到的多个原始图像I传送给CUDA中的GPU，由GPU对各分割得到的原始图像I进行所述步骤(1)到步骤(5)的并行处理；将经过所述步骤(5)处理得到的最终图像再由GPU传送给CPU，由CPU对各接收到的最终图像进行拼接得到完整的图像。

3.根据权利要求1所述的消除CCD天文图像中宇宙射线的方法，其特征在于，所述放大倍数f_z＝2，图像I^fz、Laplacian图像L^fz以及Laplacian图像L^fz′的分辨率大小均为2m×2n；

得到分辨率大小为m×n的Lplacian图像L′的像素点L′_i,j表达式为：

$L_{i,j}^{\&prime;}=\frac{1}{4}(L_{2i-\mathrm{1,2}j-1}^{{\mathrm{fz}}^{\&prime;}}+L_{2i-\mathrm{1,2}j}^{{\mathrm{fz}}^{\&prime;}}+L_{2i,2j-1}^{{\mathrm{fz}}^{\&prime;}}+L_{2i,2j}^{{\mathrm{fz}}^{\&prime;}});i=\mathrm{1,2},...,m,j=\mathrm{1,2},...,n;$

其中均为Laplacian图像L^fz′的像素点。

4.根据权利要求1所述的消除CCD天文图像中宇宙射线的方法，其特征在于，所述步骤(3)中的关联于原始图像I噪声模型N为：

$N=g^{-1}\sqrt{g(M_5*I)+{\mathrm{\&sigma;}}_{\mathrm{rn}}^2};$

其中g是增益因子；σ_rn是读入电子噪声，M₅*I表示对原始图像I进行5×5的中值滤波处理；

所述步骤(3)中信噪比S为：

$S=\frac{{\left(L^{\&prime;}\right)}^{+}}{f_{z}N};$

其中f_z为子采样因子，像素点的S比值越大，该像素点所包含的噪声越大；

所述步骤(3)中S经过中值滤波处理后得到的信噪比S′为：

S′＝S-(S*M₅)；

S*M₅表示S通过5×5的中值滤波模板进行中值滤波处理。

5.根据权利要求4所述的消除CCD天文图像中宇宙射线的方法，其特征在于，所示噪声模型N中的增益因子g为7，读入电子噪声σ_rn为5。

6.根据权利要求1所述的消除CCD天文图像中宇宙射线的方法，其特征在于，所述步骤(3)中的关联于原始图像I的精细结构模型F为：

F＝(M₃*I)-[(M₃*I)*M₇]；

其中M₃与M₇分别为3×3的中值滤波模板和7×7的中值滤波模板，原始图像I与M₃进行中值滤波操作M₃*I，得到原始图像中的中、低频信息，与M₇进行中值滤波操作(M₃*I)*M₇后得到图像中的低频信息；

所述步骤(3)中的Laplacian图像L′的值与精细结构模型F作比值得到的T值为：

$T=\frac{{\left(L^{\&prime;}\right)}^{+}}{F}.$

7.根据权利要求1所述的消除CCD天文图像中宇宙射线的方法，其特征在于，所述步骤(3)中的第一阈值σ_lim为0.5，第二阈值f_lim为1.5；所述步骤(4)中的第三阈值X为0.1％。

8.根据权利要求1所述的消除CCD天文图像中宇宙射线的方法，其特征在于，所述步骤(5)中的分数阶微分的Tiansi算子归整后得到Tiansi′算子与消除宇宙射线后的图像L^v进行卷积运算得到的增强后的图像L^T为：

L^T＝Tiansi′*L^v；

其中Tiansi算子的表达式为：

Tiansi＝8-8e+8×(e²-e)/2＝8-12e+4e²＝4(e-2)(e-1)；

e代表分数阶的阶次，其取值范围为(0,1)；所述Tiansi′算子为Tiansi算子的每一项均除以4(e-2)(e-1)。

9.根据权利要求8所述的消除CCD天文图像中宇宙射线的方法，其特征在于，所述分数阶的阶次e的大小根据图像增强需要的大小做调整。