一种固定选图数的实时幸运成像算法*

王锦良，李彬华,2

(1. 昆明理工大学信息工程与自动化学院，云南 昆明 650500；2. 昆明理工大学云南省计算机应用技术重点实验室，云南 昆明 650500)

幸运成像技术是一种用于去除大气湍流效应的简单有效的图像恢复技术[1]。由于该方法需要从大量的短曝光图像中选取少量幸运图像进行配准和叠加，计算量大，通常是在天文观测完成之后的一段时间内用专门的计算机进行图像恢复或重建，难以实时处理，不过实时化仍是其发展方向之一[2]。对于大量图像的幸运重建来说，传统的个人计算机无能为力，因此，近年来具有强大并行处理能力的图形处理器(Graphics Processing Unit, GPU)和现场可编程门阵列作为硬件加速器引入图像处理领域[3-5]。虽然图形处理器具有强大的并行处理能力，但其实现的系统体积较为庞大，处理速度与现场可编程门阵列相比没有任何优势，这一结论是文[6]分别在图形处理器和现场可编程门阵列上实现幸运区域融合(Lucky-region Fusion, LRF)算法的实验中得出的。

2008年，文[7-8]首次将幸运成像技术在现场可编程门阵列上实现，并获得了实时的高分辨率图像。FastCam基于实时阈值进行选图，虽然实时性有保证，但可能导致所选图像出现偏差，特别是由于FastCam系统没有剔除包含宇宙射线的图像，当输入宇宙射线图像后，会选中这帧图像并使实时阈值增大，导致后续的好图不再被选中，这些不足使得系统的处理能力和成像效果受到一定限制，而且他们仅介绍了算法的基本框架和实验结果，并没有介绍各个处理模块的具体设计方法和实现技术。

2018年，文[9]首次对幸运成像算法的现场可编程门阵列实现进行了细致的描述。虽然文[9]中的系统在功能上是完善的，但由于在图像传输链路采用数据速率较低的SD卡接口，同时幸运成像算法与基于中央处理器的算法是一致的，虽然各模块功能和时序划分比较独立，但综合的资源利用率和时钟效率较低，导致该实验系统的处理能力和速度都没有达到实时化的要求，性能与FastCam系统有较大差距。

文[10]介绍了一种新的实时幸运成像算法及实现技术，其中采用的滤波、宇宙射线检测与剔除、选图、配准、叠加等子算法与FastCam和文[9]的算法完全不同，由于采用固定选图比的选图方法，所选出的好图与传统方法稍有差异，导致幸运成像结果也有微小的偏差。

本文以文[9]的幸运成像系统为基础，采用千兆网作为图像输入通道，提出了一种具有实时处理能力的基于现场可编程门阵列的幸运成像算法，该算法创新性地采用一种固定选图数且无排序的图像选择策略和一种以行列坐标为基准的图像配准策略，能够有效节省算法处理时间和硬件资源，算法在好图的配准和叠加等环节与文[10]的方法也有差异。然后在一块中小规模的现场可编程门阵列开发板上实现这一算法，并构建了一个紧凑的实时化的高速幸运成像系统，能够实现大量图像的幸运重建。最后利用云南天文台丽江观测站2.4 m天文望远镜拍摄的大量短曝光图像进行现场可编程门阵列算法设计的可行性验证，并对系统的性能进行了分析讨论。

1 实时幸运成像算法的现场可编程门阵列实现

由于研究工作的延续性，本文提出的幸运成像算法的硬件实现平台和设计开发软件环境与文[9]的系统相同，即以Xilinx公司Spartan-6系列XC6SLX16芯片为核心的开发板和一台中央处理器为Xeon E5620的工作站，设计开发环境是ISE Design Suite 14.7，采用ChipScope pro 14.7和ModelSimSE 10.5进行硬件设计的验证与调试。在研究过程中，首先复现了文[9]系统幸运成像算法的现场可编程门阵列实现过程，然后与本系统幸运成像算法的现场可编程门阵列实现过程作对比。

经典的幸运成像算法包括短曝光图像的预处理、选图、配准、叠加4个流程[3]，用现场可编程门阵列实现幸运成像算法，也应该有这几个流程。本文基于实时处理的要求，以新的选图、配准、叠加思想设计了相应的算法及实现模块。

1.1 算法的总体设计

1.1.1 算法概述

一个设计合理的基于现场可编程门阵列的实时算法，基本思想与传统的基于中央处理器的算法是一致的，但必须摒弃传统的串行处理的习惯思维，尽可能多地采用并行处理思想进行设计和实现。所以，在设计新的基于现场可编程门阵列的实时幸运成像算法时，仍然采用传统的基于个人计算机或图形处理器的幸运成像算法[1,3]的关键思想——选图、配准和叠加，但数据处理方式与步骤则按并行思想进行。首先将图像数据分成两路，一路进行图像的高斯滤波及数据暂存，一路是宇宙射线检测和选图、配准；然后并行叠加，得到高分辨率图像；最后，采用锐化操作对高分辨率图像进行增强，并送往视频图形阵列(Video Graphics Array, VGA)显示。算法的关键步骤描述如下：

(1)图像滤波，同时并行地查找图像峰值；

(2)保存滤波后图像，同时并行地检测剔除宇宙射线；

(3)按固定的选图数进行选图；

(4)按图像序号和行列坐标进行配准；

(5)并行叠加；

(6)对叠加图像做反转变换进行增强并保存重建图像；

(7)送往视频图形阵列显示。

需要说明的是，由于现场可编程门阵列具有高度的并行处理能力，不仅在上述(1)、(2)两步中是并行处理的，后面的步骤也可以跟图像输入、滤波等过程同时进行，这是基于现场可编程门阵列的处理算法与基于中央处理器的幸运成像算法的最大不同点。

传统的基于中央处理器的幸运成像算法是一种事后图像处理方法，是在所有图像采集并传输完成后，按照固定的选图比选取固定数目的好图。因此，对于传统的幸运成像，固定选图比与固定选图数的实质是一样的。但对于实时处理而言，需要在图像采集和传输的同时进行图像选取，按照固定选图比的方式会丢失少量的好图[10]，而在一定传输帧数下按照固定选图数的方式则不会丢失任何一帧好图。从上述算法流程可见，本算法选图是按照固定的选图数目进行的，虽然这既与传统的采用固定选图比的选图方案[1,3,9]不同，也与采用固定选图比的动态选图方案[10]不同，但当输入图像达到设定值时，所选出的幸运图像(即好图)与传统的采用固定选图比选出的好图结果是一致的。由于采用固定的选图数目，也没有选择阈值，所以有别于FastCam中的实时阈值选图方案，这是本算法的一个创新点。以上算法中各步骤的详细处理过程构成相应的一些子算法，它们与现场可编程门阵列的结构相关，将在后面相关小节中介绍。

1.1.2 现场可编程门阵列系统的总体设计

现场可编程门阵列是一个半定制硬件，基于这种器件的算法设计，与传统的基于中央处理器的算法设计差别巨大，因为现场可编程门阵列算法的实现与硬件的结构、资源的种类和大小等密切相关。对于上面所述的幸运成像算法，主要的步骤可用一个现场可编程门阵列模块电路实现，再结合图像输入输出方式，实现的系统主要包括千兆以太网接收模块，高斯滤波模块，宇宙射线剔除模块，数据输入缓存模块，数据读写控制与存储模块，短曝光图像的选图、配准、叠加模块，图像灰度变换模块，重建图像保存模块以及视频图形阵列驱动显示模块，系统的总体结构如图1，其中千兆网、滤波和射线剔除这3个模块是文[9]的系统中没有的，选图、配准和叠加这3个模块重新进行了设计，灰度变换模块采用反转变换这种简单易实现的锐化操作对高分辨率图像进行增强，其他模块做了少量的优化。本系统DDR3工作频率为312.5 MHz(工作频率为625时，系统也可以正常工作)。

图1 实时幸运成像算法的总体结构框图Fig.1 Overall structure of the real-time lucky imaging system

由于受现场可编程门阵列芯片RAM资源以及硬件逻辑资源的限制，系统目前的目标与FastCam一致，可以使用2 000帧128×128像素的短曝光图像作为原始图像进行选图，在配准过程中截取64×64像素的图片进行配准，并最终以64×64像素显示幸运成像结果。

由于千兆以太网接收模块的代码是由现场可编程门阵列开发板厂商提供的代码修改得出的，其他类似的模块可参阅文[9]，高斯滤波模块与宇宙射线剔除模块可参阅文[10]。以下主要介绍新的思想方法重新设计选图、配准和叠加这3个算法及实现模块。

1.2 选图算法及实现模块设计

传统的选图思想是基于每帧图像的瞬时斯特列尔比进行比较排序，当图像较多时，选图算法计算量巨大，消耗的资源很多，不适合在中小规模的现场可编程门阵列系统上实现。为此，本文采用一种固定选图数目且无需排序的选图思路，符合现场可编程门阵列的逻辑设计，可节省大量的时钟和硬件资源，算法示意图如图2。

图2 选图算法示意图Fig.2 Schematic diagram of the image selection algorithm

图2中m代表固定的选图数目，也就是从n帧天文图像中选出m个峰值进行配准叠加，本系统实现的是从2 000帧天文图像中按1%的选图比[3,7]选出20个峰值进行配准叠加，故系统的m=20，n=2 000。系统每收到一个经过宇宙射线剔除模块后的峰值都暂存到RAM地址m+1上，然后依次读取RAM地址1到m+1，相互比较并记录最小值所在的RAM地址k，一轮比较结束后，把暂存在地址m+1上的峰值写入最小值所在存储地址单元k中，如此往复，直到n帧图像峰值读取完毕，RAM中地址1到m的数据就是从n个图像峰值中的前m个最大值依次读出，选图结束。

本文选图思想与文[10]不同，本文采用固定选图数目的方法，即m是不变的，而文[10]采用固定选图比的方法，即m是变化的。相应地，现场可编程门阵列实现的方法和Verilog HDL代码也大不同。文[10]的优点是可以随图像的输入动态选图，但代码设计复杂，代码量大，特别是选图结果与传统的算法有小的偏差[10]。本文算法的优点是在确定的输入图像数时选图结果与传统的算法完全一致，即不会丢失任何一帧好图，代码设计简单。

1.3 配准算法及实现模块设计

本模块设计的目的是对经过选图模块选出的128×128像素的图像以峰值所在像素为图像中心，截取所在行列以及前31行和后32行、前31列和后32列，形成一帧64×64像素的图像进行叠加处理。设计中的主要问题是如何找到以峰值所在像素为中心的64×64像素图像的第1个像素在DDR3存储器中的地址，用于从DDR3存储器中读取相应的图像像素的灰度值，送入叠加模块处理。本文根据DDR3存储器数据存储原理以及图像配准要求，采用峰值所在图像序号和峰值所在像素位置(即横纵坐标)进行图像配准，得到了符合现场可编程门阵列逻辑实现的配准公式：

addr_r=(pic_cnt-1)×65 536+[(axis_x-32)×128+axis_y-32]×4 ，

(1)

其中，addr_r为以峰值所在像素为中心的64×64像素图像的第1个像素在DDR3存储器中的地址；pic_cnt为信息寄存器中寄存的图像序号；axis_x和axis_y分别为信息寄存器中寄存的峰值所在像素的横坐标和纵坐标。(1)式中乘4是因为一个像素在DDR3存储器中占4个地址单元，65 536(128×128×4)代表1帧图像在DDR3存储器中占用的地址单元的个数。需要注意的是，如果一帧图像中存在两个或多个相同的峰值时，系统取第1个峰值为中心进行配准。

另外，由于本文中好图在DDR3中的存取方法与文[10]的方法不同，配准公式(1)与文[10]的配准公式有明显的区别。

1.4 叠加算法及实现模块设计

本模块实现的功能是将配准之后的图像进行叠加，将叠加后的图像送入灰度变换模块，为了实现图像数据的无缝缓冲和处理任务，并节省缓冲区空间，系统采用一个16位真双口RAM和一个16位简单双口RAM进行 “乒乓” 操作，实现图像的叠加。“乒乓” 操作主要依靠切换两个RAM之间的读写操作实现，即当RAM1进行写操作时，RAM2进行读操作，当一帧图像叠加完毕后，切换状态，即RAM1进行读操作，RAM2进行写操作，直到m帧图像叠加完毕，此时RAM2中存储的就是全部叠加后的图像。

为了验证模块设计的正确性，系统采用固定20帧选图(即2 000帧×1%)在现场可编程门阵列上做预处理、选图、配准及叠加操作，并用ChipScope软件对最终结果图中峰值周围部分像素进行抓取，结果如图3(a)，在MATLAB上进行同样的处理，结果如图3(b)。

图3 叠加模块结果。(a)现场可编程门阵列处理的结果；(b)MATLAB处理的结果Fig.3 Superposition module processing results. (a) Results of FPGA processing; (b) Results of MATLAB processing

对比图3(a)和图3(b)可知，经现场可编程门阵列实现的结果与在MATLAB上实现的结果完全一样，说明叠加模块设计正确，同时也验证了本文提出的幸运成像算法的各主要模块及其互联接口的设计是正确的。

2 实验结果及分析

为了验证本文提出的基于现场可编程门阵列的幸运成像算法的性能以及算法实现的正确性，首先用同样的算法在基于个人计算机的MATLAB平台上对天文目标的短曝光图像进行处理，其次在基于现场可编程门阵列开发板的硬件平台上进行同样的处理，最后将相关结果做对比分析。

实验使用的天文目标短曝光系列图像是2016年10月20日在云南天文台丽江观测站对双星HDS 70的实测图像，共10 000帧，这组图像的观测条件和参数见文[3]。在图像帧数以及图像大小的选取上，由于现场可编程门阵列的限制，多次随机从10 000帧图像中抽取2 000帧，并裁剪成128×128像素的天文目标短曝光图像进行处理。

2.1 幸运成像算法在CPU+MATLAB平台上的处理实验

本实验的硬件平台是：Dell Precision T5500图像工作站，其中内存16 GB，中央处理器为Xeon E5620，显卡为NVIDIA GTX1080，运行的软件环境是Windows 7操作系统、MATLAB R2017a。

在MATLAB中对随机抽取的2 000帧128×128像素的图像进行幸运成像算法处理，在处理时只截取64×64像素的目标区域。根据文[3,7]的研究结果，选图比为1%时幸运成像效果最好，也就是从2 000帧图像中选20帧做叠加，所得结果如图4。对图4(a)的高分辨率图像进一步分析，采用二维修正矩算法[11]对目标(主星和伴星)进行计算，得到以像素为单位的主星坐标(33.015，32.147)和伴星坐标(37.906，26.068)，由此计算出两星之间的距离约为7.802像素，再根据望远镜成像系统的比例尺0.043″/pixel[3]，最终得到双星间距为7.802×0.043≈0.335″，这个结果与文[3]的结果相近，实验结果说明用MATLAB实现的幸运成像算法是正确的。

图4 2 000帧短曝光图像中选20帧时所得的幸运成像结果。(a)二维灰度图；(b)三维灰度分布图Fig.4 Results of 20 frames of lucky images selected from 2 000 frames of short exposure.(a) 2D gray image; (b) 3D gray distribution map

2.2 幸运成像算法在现场可编程门阵列上的实验验证

本实验的硬件平台是Xilinx公司Spartan-6系列XC6SLX16芯片为核心的开发板。系统在现场可编程门阵列上对幸运成像算法设计进行验证，实时输入2 000帧原始图像，采用固定20帧好图(相当于总帧数1%的选图比例)。由于验证所用图片总帧数以及选图比例过低，导致最终得到的高分辨率图像像素灰度过低，所以无论在MATLAB上还是现场可编程门阵列上进行算法实现所得的高分辨率图像都比较模糊，故系统在现场可编程门阵列和MATLAB上均对算法所得的高分辨率图像做了相同的灰度变换——锐化，将最终得到的高分辨率图像中感兴趣的目标或区域突出出来，分别得到了如图5(a)，(b)的锐化图像。

图5(a)是通过以太网回传最终高分辨率图像锐化结果(数据)给个人计算机并调用MATLAB的函数显示的结果图，图5(b)是在MATLAB上进行幸运成像算法处理并做相同锐化的结果图。通过对比可以看出两帧图像完全相同，从而验证了系统在现场可编程门阵列上实现的幸运成像算法正确。

图5 基于现场可编程门阵列和MATLAB实现的幸运成像算法处理所得的结果图像。(a)现场可编程门阵列处理的结果；(b)MATLAB处理的结果Fig.5 Resultant images Obtained by lucky imaging algorithm based on the FPGA and the MATLAB.(a) Result of FPGA processing; (b) Result of MATLAB processing

2.3 对比实验及结果分析

为了体现基于固定选图数的幸运成像算法的优越性，本文设计了两个对比实验。第1个对比实验：分别用基于固定选图数的算法、文[9]的基于排序的算法、文[10]的基于固定选图比的算法以及在MATLAB上实现的相同算法对1 000帧原始图像按1%的选图比例进行幸运成像实验，就算法运行时间、图像数据读取时间和总运行时间进行对比，实验结果如表1。

第2个对比实验：对同帧数同大小的原始图像在基于固定选图数的幸运成像算法、FastCam项目[8]中基于实时阈值的幸运成像算法以及文[10]中基于固定选图比的幸运成像算法的资源消耗进行对比，结果如表2。

表1和表2中，本系统算法的现场可编程门阵列实现和第1代系统算法的现场可编程门阵列实现的硬件平台均为2.2节所述的现场可编程门阵列开发板，CPU+MATLAB实现的硬件平台为2.1节所述的平台。

表1 不同计算平台下1 000帧原始图像的幸运成像算法运行时间(单位： s)Table 1 Running times (in seconds)

表2 不同计算平台下2 048帧128×128像素大小的原始图像幸运成像算法资源消耗Table 2 FPGA on-chip resource consumption of lucky imaging for 2 048 frames of cropped images with 128×128 pixels

从表1可以看出，本文的幸运成像算法用现场可编程门阵列实现后的处理速度比文[9]的系统快27倍，同比提升了96.33%，比CPU+MATLAB实现的算法处理速度快约152倍。虽然与固定选图比算法处理速度几无差异，但本文算法与传统算法结果完全一致，而固定选图比算法与传统算法结果稍有差异。

此外，从表2可以看出，本文的幸运成像算法用现场可编程门阵列实现的片资源消耗仅为FastCam系统中片资源消耗的34.91%，块RAM资源消耗仅为14.28%，相比固定选图比算法实现的资源消耗也略有优势。

通过以上对比分析可知，本文的幸运成像算法的现场可编程门阵列实现在处理速度和资源消耗上都具有一定的优势，图像帧率可达197 frame/s(约5 ms/frame)。天文幸运成像观测通常使用电子倍增电荷耦合器件(Electron Multiplying Charge-Coupled Devices, EMCCD)的相机，这是一种高速低噪声的科学级相机[12]。实际观测时，相机的曝光时间通常不小于20 ms，例如，FastCam系统的曝光时间是30 ms[7]，而丽江2.4 m望远镜观测时的曝光时间是20 ms[3]，所以，本文的基于现场可编程门阵列的幸运成像算法及实现技术可以用于构建实时的幸运成像系统。

3 结 论

本文在仔细分析已有实时幸运成像算法和系统优缺点的基础上，充分利用现场可编程门阵列的并行性与灵活性优势，提出了一个符合现场可编程门阵列开发特点的实时幸运成像算法。该算法在选图方面摒弃传统的排序选图思想，采用固定选图数且无需排序的方式；在图像配准方面，采用好图的二维行列坐标进行配准，公式简洁明了。文中介绍了实时幸运成像算法的流程和各关键步骤对应模块的逻辑设计方法和实现技术，以及模块的测试结果与验证过程，以此为基础，在一块中小规模低成本的现场可编程门阵列开发板上进行了系统的实现，构建了一个具有实时处理能力的幸运成像系统，新的幸运成像算法在中小规模的现场可编程门阵列上可以处理2 000帧以上的图像，达到或者略超FastCam的处理能力，但成像结果与传统算法一致，资源消耗更少。由于处理图像(128×128像素)的平均帧率可高达197 frame/s，所以，本文的幸运成像算法及现场可编程门阵列实现技术可以用于构建实时幸运成像系统。

致谢：中国科学院云南天文台张西亮、季凯帆、金振宇为本研究工作提供了原始天文图像，特此致谢。