基于现场可编程门阵列的幸运成像算法的实现*

赵盼孜，李彬华，毛栊哗，陶 勇

(昆明理工大学信息工程与自动化学院，云南 昆明 650051)

地基大口径光学望远镜对天体成像的分辨率受制于大气湍流[1]。为了消除大气湍流对目标成像质量的影响，常见且有效的方法是采用自适应光学(Adaptive Optics, AO)技术，以及事后处理方式的图像复原技术，幸运成像技术就是其中之一。幸运成像技术利用部分短曝光图像中包含的目标高分辨率的结构信息，在观测后进行图像重建，以消除大气湍流的干扰从而获得高分辨率图像[1-2]。二十一世纪初出现的电子倍增电荷耦合器件(Electron Multiplying Charge-Coupled Devices, EMCCD)成像技术具有高分辨率、高读出速度、低噪声等优点，适合于曝光时间在毫秒及以上量级的微光成像领域[3]。该技术的出现使幸运成像技术获得了长足的进步，并在天文观测领域取得了极大的成功，2013年以前发表的论文已超过200篇[2]。典型的幸运成像系统包括英国剑桥大学的LuckyCam、德国马普研究所的AstraLux和西班牙加那利群岛天文研究所的FastCam，它们在双星或多星系统中暗伴星的发现、天测天力参数计算等方面，科学产出极为丰富[4-5]。

幸运成像技术是一种简单可行的图像复原方法，但由于图像复原是在观测完成后的一段时间内进行，其缺点明显，即天文观测人员对于所拍摄图像的实时信息了解不多，难以及时发现并纠正观测中可能存在的偏差或错误。解决此问题的一种办法是改进算法，增加硬件的处理能力[3]，将幸运成像技术实时或准实时化。

传统的中央处理器(Central Processing Unit, CPU)对图像数据的处理采用串行方式，对于大量图像的幸运重建来说，难以实现实时或准实时化。相较于中央处理器，现场可编程门阵列(Field Programmable Gate Array, FPGA)器件有着先天的并行性和灵活性优势，能提供强大的并行计算能力，可加速数据或信号的处理，特别是图像处理[6]，在天文观测和数据设备中使用越来越多[7]。将现场可编程门阵列用于幸运成像，是将该技术实时或准实时化的一种有效方法。2008年FastCam项目研究过程中，首次将幸运成像技术在现场可编程门阵列上实现，并获得了实时的高分辨率图像[1,8]；文[9]将幸运成像算法在现场可编程门阵列和图形处理器中实现。不过，他们并没有对算法的现场可编程门阵列实现进行细致的描述。现场可编程门阵列与幸运成像算法结合的研究工作国内尚无报道。

本文对传统的基于中央处理器和程序设计语言的幸运成像基本算法进行了简短的分析，结合现场可编程门阵列的特点，提出了一种适合现场可编程门阵列处理的幸运成像方案，用Verilog硬件描述语言(HDL)进行各功能模块的设计，并将幸运成像算法成功地在现场可编程门阵列上实现，最后利用中国科学院云南天文台丽江观测站2.4 m天文望远镜拍摄的大量目标短曝光图像进行硬件算法设计的可行性验证。

1 基于现场可编程门阵列的幸运成像算法处理流程

幸运成像算法的基本思想是按照一定的评价标准，从一系列短曝光图像中选出符合评价标准的图像或区域，之后再对其进行配准叠加处理，以获得不受或少受大气湍流影响的高分辨率图像[2]。本文构建的幸运成像系统在实现过程中按照基本的幸运成像算法的设计要求进行，只是为了适应本系统所用现场可编程门阵列的处理逻辑与资源限制，在不改变幸运成像算法基本处理流程(选图、配准、叠加)的基础上将该算法用Verilog HDL重新进行了设计。

1.1 幸运图像的选取

在一系列短曝光序列图像中选出 “幸运图像” 是最终获得高分辨图像的关键。然而为获得 “幸运图像” 必须选取一个合适的像质评价函数。对于空间点源目标图像来说，像质评价函数通常采用图像的斯特列尔比(Strehl Ratio, SR)作为评价标准。由于斯特列尔比定义为存在像差时图像的峰值光强与不存在像差时图像的峰值光强之比，因此，计算时需要理想的无像差峰值光强。但这个理想的峰值光强在现实中不容易获得，故通常采用其它的替代方案。一个简单有效的方法是使用瞬时的斯特列尔比作为图像的像质评价函数[2,10]，即只对图像上的峰值光强像素点的灰度值进行斯特列尔比值估算。

对于现场可编程门阵列这类运算资源有限的硬件，计算斯特列尔比只能采用简单的处理方式。具体来说，在系统实现过程中，直接使用有像差时图像的峰值光强像素点的灰度值。由于在配准过程中对每帧图像的峰值光强像素点灰度值的位置有要求，故在系统选图过程中，以图像的峰值光强像素点灰度值及其位置共同对图像进行像质评价。

1.2 图像的配准与叠加

配准是幸运成像处理中最重要的一环，如果配准不当会直接导致叠加后的图像模糊不清，分辨率下降，从而无法将暗弱目标图像显示出来。由于系统使用恒星的图像，所以最常用的配准方法有两种：(1)以整幅图像中的最大灰度值为中心截取所需的成像区域；(2)以图像的质心为配准基点。由于受现场可编程门阵列逻辑及资源的限制，系统采用前一种配准方法。

图像配准完成后，需要对配准图像进行叠加。在系统中由于所选图像数量偏小，所以采用直接叠加方法。叠加后所得图像进行灰度变换(即锐化)可得到易于观察的高分辨率图像。

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

本文所述的幸运成像算法实现的硬件平台是以Xilinx公司Spartan-6系列的XC6SLX16芯片为核心开发板。开发环境为ISE Design Suite 14.7，逻辑设计所用HDL是Verilog，采用ChipScope pro 14.7和ModelSim 10.1d进行硬件设计的验证与调试。

2.1 算法总体设计

由于项目设计的主要目的是在现场可编程门阵列上实现幸运成像算法，因此为读取图片像素数据，首先将天文短曝光图像存储在开发板外挂的安全数据卡(Secure Digital Card, SD卡)中，以便给算法模块读取图像数据。此外，为满足算法对速度及片内存储的要求，同时由于所用芯片逻辑资源和存储空间的限制，故在数据处理方式上采用数据流，即逐像素处理的方式，但在每个模块内以及各个模块间均采用并行的方式对数据进行处理。与此同时在片内数据的存储上为使处理速度和存储器之间达到最佳状态，采用片内的双端口随机存储器(Random Access Memory, RAM)模块作为数据缓存器。然而，由于所用现场可编程门阵列芯片随机存储器资源以及逻辑资源的约束，只能对1 000张128 × 128像素大小的短曝光图像作选图处理，并且在配准过程中只能最大截取64 × 64像素大小的图片进行配准处理，并最终以64 × 64像素大小显示幸运成像结果。

对于整个算法的硬件构架，为了算法实现过程中修改及调试的方便，采用模块化的设计方式。主要模块包括数据的读写、保存模块，短曝光图像的选图、配准、叠加模块以及最终的重建高分辨率图像显示模块，并用Verilog HDL设计实现。基于现场可编程门阵列实现的幸运成像算法的结构框图如图1。

图1 基于现场可编程门阵列的幸运成像算法的结构框图
Fig.1 Structure block diagram of the lucky imaging algorithm based on FPGA

幸运成像的现场可编程门阵列的工作流程如下：系统上电后，安全数据卡中存储的天文图像数据会通过串行外设接口(Serial Peripheral Interface, SPI)以6.25 Mb/s的速度不断地向安全数据卡顶层模块传送，经安全数据卡顶层模块处理过的数据经先入先出(First-In First-Out, FIFO)缓冲器向第3代双倍数据率同步动态随机存取存储器(Double-Data-Rate Three Synchronous Dynamic Random Access Memory, DDR3)写数据模块中写入(DDR3速度666 Mb/s)，然后通过DDR3写数据模块的处理，向现场可编程门阵列芯片外部的DDR3存储器中写入图像数据。在写数据模块向DDR3写入像素数据的同时，选图模块不断地对每帧图片的像素值进行最大灰度值求解，然后对求出的所有图片的最大灰度值排序，当排序结束，说明选图结束，产生选图结束信号。

其后，向配准模块发送要截取的原始图片的序号以及该图片最大灰度值的位置参数，以便配准模块向DDR3读数据模块发送所选图片的首像素地址。与此同时，选图结束信号也会发送给DDR3读数据模块，当选图结束信号和所选图像首像素地址同时作用于DDR3读数据模块时，该模块将数据从DDR3存储器中读出满足条件的图片像素进行叠加处理。从图1看出，配准、DDR3读出与叠加这3个模块是串行的，但用HDL设计及现场可编程门阵列实现时，它们是并行的，它们的运行从时间上看是重叠的。

当所选的全部图片叠加完成后，将结果进行分段线性灰度变换的处理，以增强目标区域的可视性。然后将最终高分辨率图像的像素数据保存在片内随机存储器中。当保存完毕后，通过视频图形阵列(Video Graphics Array, VGA)驱动模块驱动，将灰度图在视频图形阵列显示器上显示。

2.2 选图模块的设计

在选图模块的实现上，系统主要采用最大灰度值求解模块以及最大灰度值排序模块共同搭建而成，这两个子模块的搭建主要由比较器、计数器以及片内随机存储器完成。其中，最大灰度值查找模块的状态转移图如图2。

图2 最大灰度值查找的状态转移图
Fig.2 State transfer diagram for looking up the maximum gray value

在图2中，myvaild=1表示安全数据卡中的数据可以取出；cnt表示已读取的一帧图像的像素数目。其中，cnt=0表示一帧图像的全部像素已读完；eve_t=1表示单帧图像的最大灰度值开始取出；pic_num表示图像的帧数；wea_max=1表示可以开始向随机存储器中写数据。从2.1节中可知，在本设计中所用的图像大小为128 × 128，帧数为1 000，所以cnt最大为16 383，pic_num最大为999。

从图2可以看出，该模块的主要功能是对从安全数据卡读出的图片像素值逐个进行比较，找出每帧图像的最大灰度值并将其保存在片内随机存储器中供排序使用。值得注意的是，在每求出一帧图像最大灰度值后，在条件触发时先将其缓存，然后等相应条件到来时将其进行判断，看是否满足配准对最大灰度值位置的要求。满足则在随机存储器相应的地址空间中保存该像素值，不满足则在随机存储器相应的地址空间中存入与该最大灰度值相同比特位数的零值。其次，在向随机存储器中保存每帧图像最大灰度值时，将该最大灰度值所在的图片序号及位置参数一并保存在随机存储器相应地址空间，以便进行并行、串行混合排序。

2.3 配准模块的设计

在配准算法的实现上，采用以最大灰度值为图片中心、截取64 × 64像素大小进行叠加处理，其设计主要通过选图模块中发送的最大灰度值所在图片序号以及其位置参数，再根据(1)式计算出要截取的图片首地址，然后发送给DDR3读数据模块，用从DDR3存储器中读出相应的图片像素值，供叠加模块处理，该模块并不占用片内随机存储器资源，只是少量地使用了逻辑片(Slices)资源。

c3_p0_cmd_addr_r=((pic_cnt-1)×65 536)+

(1)

其中，c3_p0_cmd_addr_r表示要读取的DDR3的地址；pic_cnt表示要读取的最大灰度值所在的图片序号；max_cnt表示要读取的最大灰度值在图片中的位置参数；INT(max_cnt/27)表示max_cnt除以27之后只取整数部分，在现场可编程门阵列中用移位寄存器实现。

2.4 叠加模块的设计

对于图像的叠加主要是通过片内随机存储器资源的使用实现的，在叠加模块设计中占用的片内随机存储器个数为12，消耗约37.5%的存储器资源。该模块的结构框图如图3。

图3 叠加模块结构框图
Fig.3 Block diagram of stacking module

3 实验结果及分析

为了验证幸运成像算法在现场可编程门阵列上实现的正确性，拟用MATLAB对相同的幸运成像算法进行编程并将其在计算机上的处理结果与现场可编程门阵列处理结果进行对比。对于现场可编程门阵列上运行幸运成像算法的处理时间问题，只能与计算机的运行时间做比较。所以本文设计了两个实验：(1)与MATLAB处理结果的对比实验；(2)系统自身处理结果与数据的实验。

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

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

实验硬件平台是Dell Precision T5500图像工作站、16 GB内存、NVIDIA GTX1050Ti显卡；运行的软件环境是Windows 7操作系统、MATLAB R2014a。

将随机抽取的1 000帧128 × 128像素大小的图像进行分组，用MATLAB分别对每组图像进行同样的算法处理，在实际处理中只截取64 × 64像素大小的目标区域。根据文[10]的研究结果，选图比取为1%，得到如图4的高分辨率图像和三维灰度分布图。这一结果与文[10]的结果一致。

图4 1%选图比下所得的高分辨率图像。 (a) 二维灰度图；(b) 三维灰度分布图
Fig.4 High-resolution image obtained at 1% selection rate. (a) 2D gray image; (b) 3D gray distribution map

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

硬件平台是以Xilinx公司Spartan-6系列的XC6SLX16芯片为核心的开发板。基于现场可编程门阵列实现的幸运成像算法设计的验证，采用1%的选图比，即从1 000帧选10帧，执行与MATLAB相同幸运成像的算法处理，然后用ChipScope捕捉最终结果图中最大灰度值周围部分像素数据值，得到了如图5(a)图的经现场可编程门阵列实现的幸运成像算法处理所得的高分辨率图像部分像素数据值，并与如图5(b)图的经MATLAB实现的幸运成像算法处理所得的高分辨率图像的相同部分像素值做对比。

图5 分别用现场可编程门阵列与MATLAB实现的幸运成像算法处理所得图像像素数值

(a) 现场可编程门阵列处理的结果; (b) MATLAB处理的结果
Fig.5 Image pixel values processed by FPGA and MATLAB respectively(a) Results of FPGA processing; (b) Results of MATLAB processing

通过图5可知，经现场可编程门阵列实现的幸运成像算法处理所得图像的最终像素数据与在MATLAB上做同样算法处理所得图像像素数据完全相同。

从处理效果上看，最终得到的高分辨率图像像素值较小，直接通过显示器显示图像时，无论在MATLAB上还是现场可编程门阵列做算法处理所得的高分辨率图像都比较模糊。所以，在现场可编程门阵列和MATLAB上均对算法处理所得的高分辨率图像做了相同的灰度变换处理——锐化或图像增强，将最终得到的高分辨率图像中感兴趣的目标或区域突显出来，分别得到了如图6(a)、图6(b)的最终高分辨率图像。

图6 基于现场可编程门阵列和MATLAB实现的幸运成像算法处理所得的高分辨率图像(a) 现场可编程门阵列处理的结果; (b) MATLAB处理的结果
Fig.6 High-resolution images Obtained by lucky imaging algorithm based on the FPGA and the MATLAB(a) Results of FPGA processing; (b) Results of MATLAB processing

图6(a)是经现场可编程门阵列实现的幸运成像算法处理所得的高分辨率图像以1 024 × 768的分辨率在扩展图形阵列(Extended Graphics Array, XGA)显示器上显示并通过手机拍摄得到的，在此采用手机拍摄获得，是因为系统实现所采用的现场可编程门阵列是一个独立的外部系统，它是脱离个人计算机单独工作的，所以基于现场可编程门阵列实现的幸运成像算法处理所得的高分辨率图像要通过单独的扩展图形阵列显示器进行显示，因此它与个人计算机上MATLAB所得高分辨率结果图获取方式不同，它不能像MATLAB一样直接进行图像截屏，而只能采用外部的图像采集设备进行图像获取，所以经现场可编程门阵列实现的幸运成像算法处理所得的高分辨率图像在本文中采用手机拍摄得到。而图6(b)经MATLAB处理后显示的最终高分辨率图像是以1 280 × 768的宽扩展图形阵列(Wide Extended Graphics Array, WXGA)分辨率显示的。由于两图的显示器分辨率不同导致纵横比不同，导致相同图像在显示上稍有差异。此外，图6(a)是手机拍摄获得的，由于手机拍摄角度及分辨率的不同导致图像有些许变化。不过，可看出这两幅图像是基本相同的，说明经现场可编程门阵列实现的幸运成像算法处理后的结果与经MATLAB实现的幸运成像算法处理后的结果一致，从而说明基于现场可编程门阵列实现的幸运成像算法设计正确。

综上所述，无论是从最终高分辨率图像中的像素数据值对比来看，还是从灰度变换后的直观显示最终高分辨率图像对比来看，都说明了基于现场可编程门阵列实现的幸运成像算法设计与实现的正确性。

3.3 现场可编程门阵列实现幸运成像算法的优势

现场可编程门阵列的并行性和灵活性是其最大的优势，因此相较于传统的中央处理器，现场可编程门阵列对同样算法的处理速度肯定快。但对于研究的幸运成像算法来说，具体快多少，必须通过实验说明。为此，分别在现场可编程门阵列和CPU + MATLAB上对相同幸运算法的运算速度进行了测试。实验过程中，统计1 000帧图像的处理时间，选图比采用1%。现场可编程门阵列及CPU + MATLAB平台上均采用2.1节中介绍的算法实现。由此得到了不同平台下相同算法的平均运行时间：对于现场可编程门阵列是2.45 s，对于CPU + MATLAB是46.93 s；运行所得结果(高分辨率图像)分别如图6(a)、图6(b)。显然，相同的算法在不同的平台处理下，得到相同的高分辨率图像，现场可编程门阵列上算法处理速度比CPU + MATLAB平台上算法处理速度快约19倍。如果采用更先进的现场可编程门阵列，速度还将更快。

不过受制于安全数据卡的速度，总体运行时间并没有明显的优势。目前正在设计基于USB3和千兆以太网的高速数据接口，以加速整个幸运成像处理的速度。另外，采用处理能力更强、逻辑资源更多的现场可编程门阵列，可以更多地引入并行计算、减少串行处理过程，充分发挥现场可编程门阵列速度的优势。

4 结 论

本文首先简述了幸运算法的基本思想及处理流程，根据所用现场可编程门阵列的硬件资源，用Verilog HDL设计了幸运成像算法的选图、配准、叠加等关键模块，详细介绍了这些模块的实现方案，并在现场可编程门阵列上实现了这一幸运成像系统；然后，通过将现场可编程门阵列处理的结果与个人计算机上用MATLAB处理得到的结果进行对比，证明了此系统设计和实现过程的正确性；最后，通过对现场可编程门阵列与CPU + MATLAB平台下实现相同算法所需的处理时间进行对比分析，说明了现场可编程门阵列在处理速度上快19倍的优势。虽然本设计方案还有一些可以改进之处，但就幸运成像算法在现场可编程门阵列上的具体实现，为构建实时化的高速幸运成像系统探索了一种可行的技术方法。

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