基于多进程并行加速的太阳高分辨图像重建方法*

邓 涛，陈 东，代红兵，王新华,3

(1. 云南大学信息学院，云南 昆明 650504；2. 中国科学院云南天文台，云南 昆明 650216；3. 中国科学院大学，北京 100049)

太阳是影响人类活动最大的恒星，尤其是太阳活动对地球环境、气候和天气等的影响[1]。由于受到大气湍流的影响，太阳光线在穿过地球大气层时产生波前误差，光路发生偏转，观测的太阳图像出现不同程度的偏移、抖动、模糊等，导致地基望远镜的观测图像质量下降[2]。为了消除大气湍流对望远镜成像结果的影响，通常采用空间望远镜、自适应光学以及图像重建技术等方式获取太阳高分辨图像。

目前图像重建技术常用的算法主要有相位差法、多帧盲反卷积法、斑点干涉术、K-T算法、斑点掩模法、简单位移叠加法、迭代位移叠加法以及选帧位移叠加法等。这些算法都是通过大量的短曝光图像重建太阳高分辨图像，往往因为数据量大、算法复杂导致无法达到实时重建的需求。

近年，国内外在太阳高分辨图像重建算法并行化研究方面做了很多工作，并获得一定的加速比。文[3]基于双节点集群，提出一种帧选择和斑点掩模法的并行计算方法，采用多线程方法在57 s内重建一幅256 × 256 pixel图像。文[4]利用改写于KISIP(Kiepenheuer Institut für Sonnenphysik)的算法，采用图形处理器一次处理单组分块数据，实现4.2 s内重建225个128 × 128 pixel子块的相位。文[5]从1 m新真空太阳望远镜(New Vacuum Solar Telescope, NVST)的TiO通道选取一组子块图像(100 × 256 × 256 pixel)，采用OpenMP方法实现了一组子块图像的并行化，重建单帧256 × 256 pixel的子块图像，运行时间减少至2.7 s，获得2.5倍加速。文[6]采用图形处理器的统一计算设备架构(Compute Unified Device Architecture, CUDA)对光球TiO通道中的一组子块(100 × 256 × 256 pixel)实现并行化，基于统一计算设备架构方法重建单组子图的运行时间减少到约0.7 s。文[7]基于统一计算设备架构在斑点掩模算法中实现单个子块图形处理器内的并行化，采用并行方法比纯中央处理器运行的串行算法加速比达到7。

综上所述，斑点掩模法重建太阳高分辨图像时的图像并行化研究已经取得一定的成果，但大多数局限于单纯的中央处理器或图形处理器加速，而且图形处理器一次只能处理一组分块数据，没有完全发挥中央处理器/图形处理器的并行化能力。如何将多组分块数据分配到图形处理器同时并行处理，进一步提高中央处理器和图形处理器的并行计算能力和资源效率，本文提出了基于多进程并行加速的太阳高分辨图像重建方法。

1 太阳高分辨图像重建方法

1 m新真空太阳望远镜[8]坐落在抚仙湖畔，主要的观测波段有G-band(430.5 nm)、Hα(656.28 nm)和TiO(705.8 nm)。1 m新真空太阳望远镜太阳高分辨图像重建使用了两个层次：Level1位移叠加法重建色球图像，Level1+斑点掩模法重建光球或色球图像。Level1+计算复杂度比Level1高，但Level1+在视宁度好的时候重建的图像质量更好，信噪比更高。

1 m新真空太阳望远镜Level1+太阳高分辨图像重建流程有图像预处理、图像初对齐、视宁度估计、图像分块处理和子块拼接。其中，图像分块处理主要采用斑点干涉术和斑点掩模法重建太阳高分辨图像的模和相位。图像分块处理包含若干环节，其中相位递推环节数据量大，计算复杂，并且需要考虑多个递推路径的整合，在Level1+的重建流程中最为耗时。因此，子块处理的好坏严重影响重建的效果和时效。

1.1 振幅重建

在满足等晕区条件下，短曝光像是目标和系统的点扩散函数的卷积

i(x,y)=o(x,y)⊗p(x,y)，

(1)

其中，⊗为卷积符号，时域中的卷积对应于频域中的乘积。(1)式在频域中满足

I(u,v)=O(u,v)P(u,v)，

(2)

其中，I(u,v)，O(u,v)和P(u,v)分别为时域中对应项的傅里叶变换。

表5列出了2017年广东省、江苏省、北京市与上海市自然科学基金资助SCI论文合作发文量排名前10的国家和地区。广东省与江苏省、北京市、上海市自然科学基金资助SCI论文排名前10的合作国家和地区分布有相似之处，如与美国、澳大利亚、英国、加拿大、新加坡、日本、法国等经济、科技发达国家的合作较多。同时，4个省市自然科学基金资助SCI论文数排名前10的合作国家和地区也有各自的特点，如：除美国、澳大利亚、英国、加拿大、新加坡、日本、法国外，广东省与中国台湾地区和爱尔兰合作较多；江苏省与韩国和沙特阿拉伯王国合作较多；北京市与中国台湾地区和新西兰合作较多；上海市与韩国和比利时合作较多。

功率谱统计为

(3)

1.2 相位重建

目标斑点图的三重相关为

(4)

其中，ik(x),im(x)和in(x)为3个强度分布。目标斑点图的三重自相关的傅里叶变换

I(3)(u,v)=I(u)I(v)I(-u,-v)，

(5)

其中，I(u)是i(x)的傅里叶变换；u和v是二维空间频率。

目标斑点图的平均重谱

〈I(3)(u,v)〉=O(3)(u,v)〈P(3)(u,v)〉，

(6)

其中，〈〉为系综平均；〈P(3)(u,v)〉为平均斑点掩模法传递函数。在得到平均重谱后，由低频到高频的相位元逐步递推，恢复目标斑点图的全部相位。

2 多进程并行加速方法

2.1 中央处理器/图形处理器混合计算方法

为了满足整个视场线性空间平移不变性，我们需要把预处理和初对齐后的图像分割成一个个子块，然后相同位置的子块合并成子块组。1 m新真空太阳望远镜Level1+的太阳高分辨图像重建中现有的并行加速方法大多局限于中央处理器/图形处理器混合计算方法，即图形处理器一次只能处理单个子块组，不同子块组之间还是串行计算。基于中央处理器/图形处理器混合计算方法的图像分块处理流程如图1。

虽然以上方法比纯中央处理器图像分块串行处理的时间少，但这种方法存在两个问题：(1)子块组间的处理依次进行(串行)，在图形处理器处理时，中央处理器处在空闲状态，利用率不高；(2)图形处理器一次只处理一个子块组，利用率不高。为此，本文提出了基于多进程图像分块处理并行加速的方法。

2.2 多进程并行加速方法

针对现有方法存在的问题，本文引入多进程，提高中央处理器的利用率。多进程可以使图形处理器处理更多的子块组，同时提高中央处理器和图形处理器的并行化处理能力。待处理的多帧图像分割成很多子块后，子块合并为若干子块组，把所有的子块组加入任务列表，每个进程顺序选择任务列表中的一个子块组，多个进程同时操作图形处理器并行处理多个子块组。基于多进程并行加速方法的图像分块处理流程如图2。

图2 基于多进程并行加速方法的图像分块处理流程Fig.2 Image block processing flow based on multi-process parallel acceleration method

图1 基于中央处理器/图形处理器混合计算方法的图像分块处理流程

多进程并行加速方法的图像分块处理流程为：

(1)创建任务列表。创建并初始化列表，将图像分块后的子块组加入任务列表。

(2)创建进程池。创建并初始化一个进程池，并在进程池中添加适当的子进程数量。

(3)分配进程任务。由主进程依次把任务列表中的任务分配给处于空闲状态的子进程。

(4)传递参数。主要传递子块组处理过程需要的数据和参数。

(5)处理子块组。进程池开启多少子进程，就有多少个子块组同时并行计算。子块对齐、模的重建、平均重谱的计算以及初始相位的计算等都在图形处理器并行完成。图形处理器完成初始相位计算后，把初始相位数据传递给中央处理器的相应子进程开始相位递推，完成相位递推后，子进程将相位数据从中央处理器传递到图形处理器，最后进行模和相位的合成。

(6)再次分配任务。处理完子块组后处于空闲状态的子进程都返回进程池，进程池将空闲子进程信息反馈给主进程，若子任务列表仍有待执行的子块组，主进程将待执行的子块组分配给进程池中空闲的子进程。

(7)关闭进程池。主进程反复检测子块组任务列表，当任务列表中没有待处理子块组时，等待所有子进程执行完毕，进程池将所有空闲子进程信息反馈给主进程，释放所有子进程并关闭进程池。

由上述过程可知，各个子块组相互独立，互不影响，进程池中每个子进程处理相应的子块组，同时并行运行。子块对齐、模的重建、初始相位计算以及模和相位的合成都是由多个子进程在图形处理器并行完成的，但每个子块组的相位递推是由相应的子进程在中央处理器并行计算的。值得注意的是，在创建和初始化进程池时，进程并不是越多越好，创建新进程会耗费系统资源，所以进程池中只能添加合适数量的子进程，数量取决于中央处理器和图形处理器的资源能力。此外，为避免不必要的系统开销，程序并不立即撤回完成任务的子进程，而是在任务列表还有任务时选择进程复用，再次分配任务。

以Python为例的实现步骤如下：(1)多进程管理Multiprocessing模块提供Pool类，通过import命令导入Multiprocessing模块，使用multiprocessing.Pool函数导入Pool类。(2)创建并初始化任务列表blocklist=[]，将子块组加入任务列表。(3)设置指定的工作子进程数目(num)，使用pool=Pool(num)，按照指定数量子进程创建和初始化进程池，供用户调用。(4)使用pool.apply_async()，pool.map()，pool.apply()和pool.map_async()等方法，调用子块组处理函数并传递子块组任务列表参数，然后提交给进程池。当新子块组任务请求提交到进程池时，如果Pool不满，则创建一个新进程执行该任务请求，如果Pool已满，告知进程池请求等待。(5)通过Cupy模块中cupy.asarray()方法和cp.asnumpy()方法进行中央处理器和图形处理器数据传输，每个子块组计算完成后，由列表收集结果。(6)当任务列表没有待处理子块组时，使用pool.close()方法，进程池不再接受新的任务；当所有子块组计算完成后，工作进程退出。(7)最后使用pool.join()方法，等待进程池中所有的子进程结束，返回主进程。关键过程的代码实现如表1。

表1 关键过程的代码实现Table 1 Code implementation of key processes

3 结果和分析

3.1 实验环境

硬件：Intel(R)Core(TM)i5-8400 CPU@2.80 GHz(2 801 MHz)的处理器(6核)，8 GB随机存取存储器(Random Access Memory, RAM)，NVIDIA GeForce GTX 1050 Ti，4 095 MiB。软件：Microsoft Windows 10，Spyder 4.0.13，CUDA 10.2，Python 3.7.6。

3.2 并行加速结果

以图像分块处理为例，我们对第2节两种方法进行测时比较，实验结果如表2。

表2 两种方法处理时间Table 2 Processing time of two methods

为了验证代码的正确性，本文选取2020年6月6日1 m新真空太阳望远镜Hα波段的观测数据10组，每组100帧(1 028 × 1 024 pixel)，每帧采用重叠方式分割成96 × 96 pixel，将每帧相同位置对应的子块合并成一组，划分为625个子块组。我们还选取2020年6月8日1 m新真空太阳望远镜TiO波段的观测数据10组，由于受显存大小的限制，每组选取50帧(2 160 × 2 560 pixel)，每帧采用重叠方式分割成128 × 128 pixel，每帧相同位置对应的子块合并成一组，划分为1 050个子块组。

实验结果表明，Hα波段和TiO波段数据在采用多进程并行加速方法时，图像分块处理的平均时间相对较少，加速比分别为4.69和4.71。由于在中央处理器/图形处理器混合计算方法中，子块组之间仍然是串行计算，图形处理器一次只能处理单组分块数据；而在多进程并行加速方法中，子块组之间是并行计算，图形处理器能同时处理多组分块数据。基于多进程并行加速方法可提高中央处理器和图形处理器的资源利用效率，显著提高图像分块处理的速度。

在多进程并行加速方法中，不同的进程数有不同的图像分块处理时间，如图3。图3显示了不同的进程数对图像分块处理并行加速的耗时情况，随着进程池中子进程数量的增加，Hα波段的所有子块处理时间降低到157.76 s，TiO波段的所有子块处理时间降低到203.32 s，并行加速的效果明显。但是在使用6个进程以后，受到中央处理器数量、图形处理器显存以及多进程调度开销的影响，图像分块处理时间上下浮动并呈现上升的趋势。在选择合适的子进程后，使用多进程并行加速方法处理时，中央处理器和图形处理器利用率都提高了，中央处理器利用率可以达到100%。

图3 图像分块处理在不同进程数下的耗时情况

4 结 语

针对现有方法子块组间低效的串行处理导致中央处理器和图形处理器利用率不高的问题，本文提出了基于多进程并行加速太阳高分辨图像重建的方法，利用多核和多进程技术，有效提高了中央处理器/图形处理器的利用率和重建速度，研究可以为天文数据并行化处理提供借鉴参考。进一步提高中央处理器/图形处理器的并行化程度，仍然有许多亟待突破的关键问题，其中，中央处理器承担的相位递推压力最大，耗费图像分块处理过程80%的时间。下一步我们考虑基于相位递推的特点进行并行化，同时对相关算法进行优化改进，使中央处理器/图形处理器计算负载达到均衡。此外，相关研究还需要在消息传递接口(Message Passing Interface, MPI)/图形处理器异构环境中进行验证。