CN107945131B - 射电干涉阵列分布式图像重建方法 - Google Patents

Abstract

本发明公开一种射电干涉阵列分布式图像重建方法，将射电干涉阵列图像重建的逆问题转换为等式约束条件下的凸优化问题，再转换为求解该问题的对偶问题，使用原始对偶算法，每次迭代需要的总计算成本较低，从而降低了该问题的求解复杂性；同时，将数据分割成多份，能够精确快速的分割数据；并使用多计算节点进行计算，方法灵活且具有较高的并行性能，可拓展性强。本发明不需要任何的数据后处理步骤，直接得到重建图像，最终得到的图像的动态范围更大，对延展结构的天体显示更加清晰，相比传统的方法更加高效。

Description

射电干涉阵列分布式图像重建方法

技术领域

本发明涉及射电天文成像技术领域，具体涉及一种射电干涉阵列分布式图像重建方法。

背景技术

射电干涉测量对射电成像起了至关重要的作用，它比单口径射电望远镜有更高的分辨率和灵敏度。然而，射电干涉仪受阵列中天线的数量数非常有限，这极大的限制了最终观测采样测量的数量。因此，由射电望远镜观测的采样测量值重建天体的亮度分布，即图像重建，这相当于解决一个病态的逆问题。传统解决方法大多是采用洁化CLEAN算法及其演变算法，演变算法包括Cotton&Schwab提出CLEAN算法简称CS-CLEAN和多尺度CLEAN算法简称MS-CLEAN，并没有采用现代最先进的图像重建技术。

下一代射电干涉仪，例如低频阵列(Low Frequency Array，LOFAR)，默奇森宽视场阵列(Murchison Widefield Array，MWA)，澳大利亚平方公里阵列探路者(AustralianSquare Kilometre Array Pathfinder，ASKAP)和平方公里阵列(Square KilometreArray，SKA)，必定面临海量数据的处理与成像的极大挑战。这些科学大装置追求更远大、更高的科学目标，包括探测宇宙再电离，绘制大规模结构和探索宇宙磁场。若要实现这些科学目标，用来进行图像重建的方法必须是快速准确和具有较高的可拓展性。然而，上述的传统方法均是依赖本地贪婪迭代过程，当进行大数据处理时，将要消耗大量的内存和运行时间，不适合用于大规模并行计算和分布式计算。

发明内容

本发明所要解决的是现有的射电干涉阵列图像重建方法存在运行效率低、可拓展性差、依赖于复杂的参数设置的问题，提供一种射电干涉阵列分布式图像重建方法，其在天文大数据下能够快速的重建出精确的天文图像，并且具有良好的可拓展性能(或可伸缩性能)。

为解决上述问题，本发明是通过以下技术方案实现的：

射电干涉阵列分布式图像重建方法，包括步骤如下：

步骤1、中央节点将初始化参数广播并发送给计算节点；

步骤2、中央节点根据计算节点的个数n_d，将射电干涉仪的测量数据分割为n_d份；并将这n_d份数据块分别发送给n_d个计算节点；

步骤3、每个计算节点根据所得到的数据块计算当前迭代t下的2个最终对偶变量

和

并将所计算出的2个最终对偶变量

和

同时发送给中央节点；

步骤4、中央节点根据各个计算节点所返回的2个最终对偶变量

和

计算当前迭代t下的最终解

步骤5、判断是否满足收敛条件，即当前迭代t达到最大迭代次数或者当前迭代t下的最终解

使得每个数据块均满足

如果未满足收敛条件，则迭代次数t+1，并返回步骤3；

如果满足收敛条件，则当前迭代t下的最终解

即为重建的图像；

其中，Φ_j表示第j个数据块的测量矩阵，y_j表示第j个数据块的可见度数据，∈_j表示第j个数据块的噪声上限，j∈{1,2,...,n_d}，n_d表示计算节点的个数。

上述步骤2中，对测量数据进行分割的具体步骤如下：

步骤2.1、求出计算节点个数n_d的所有约数，并将这些约数按从小到大的顺序排列为一维数组n_div；

步骤2.2、找出一维数组n_div中最先大于或等于

的数v_no，并令u_no＝n_d/v_no；

步骤2.3、对测量数据的uw投影的秩统计量R_uw进行归一化并量化后，得到uw投影的量化数据s_u：

s_u＝ceil(u_no×R_uw/l_uw)

其中，ceil表示向右取整函数，R_uw表示测量数据的uw投影的秩统计量，l_uw表示测量数据的uw投影的长度；

步骤2.4、将uw投影的量化数据s_u分别等于1,2,...,u_no的下标记录下来，并根据这些下标将测量数据的投影数据uw和vw，以及可见度数据y均分割成u_no份；

步骤2.5、对测量数据的vw投影的秩统计量R_vw进行归一化并量化后，得到vw投影的量化数据s_v：

s_v＝ceil(v_no×R_vw/l_vw)

其中，ceil表示向右取整函数，R_vw表示测量数据的vw投影的秩统计量，l_vw表示测量数据的vw投影的长度；

步骤2.6、将vw投影的量化数据s_v分别等于1,2,...,v_no的下标记录下来，并根据这些下标将步骤2.4所得的被分割成u_no份的测量数据的投影数据uw和vw，以及可见度数据y均再进一步分割成v_no份，最终得到n_d份数据块。

上述步骤3中，第j个计算节点对当前迭代t下的第一最终对偶变量

的计算过程为：

首先，计算当前迭代t下的第一初始对偶变量

然后，计算当前迭代t下的第一中间对偶变量

最后，计算当前迭代t下的第一最终对偶变量

其中，

表示第j个计算节点在上次迭代t-1下的第一中间对偶变量，Ψ表示小波基，Ψ_i表示小波基Ψ中的第i个基，

表示上次迭代t-1下的最终解，下标_κ表示收敛速度系数，

表示谱范数，λ表示松弛系数，i∈{1,2,...,n_b}，n_b表示小波基的个数，j∈{1,2,...,n_d}，n_d表示计算节点的个数。

上述步骤3中，第j个计算节点对当前迭代t下的第二最终对偶变量

的计算过程为：

首先，计算当前迭代t下的第二初始对偶变量

然后，计算当前迭代t下的第二中间对偶变量

最后，计算当前迭代t下的第二最终对偶变量

其中，

表示第j个计算节点在上次迭代t-1下的第二中间对偶变量，G_j表示第j个计算节点的紧支撑核，M_j表示第j个计算节点的数据选取矩阵，F表示傅里叶变换矩阵，Z表示补零矩阵，

表示上一次迭代t-1下的最终解，y_j表第j个数据块的可见度数据，

表示第j个计算节点求解松弛变量的近似算子，λ表示松弛系数，上标^T表示转置，j∈{1,2,...,n_d}，n_d表示计算节点的个数。

上述步骤4中，中央节点对当前迭代t下的最终解

的计算过程为：

首先，计算当前迭代t下的初始解

然后，计算当前迭代t下的中间解x^(t)：

最后，计算当前迭代t下的最终解

其中，P_C表示指示函数的近似算子，x^(t-1)表示上次迭代t-1下的中间解，τ表示解的收敛常量，Z表示补零矩阵，F表示傅里叶矩阵，

表示第j个计算节点计算节点的第二收敛常量，M_j表示第j个计算节点的数据选取矩阵，σ_j表示第j个计算节点的第一收敛常量，

表示第j个计算节点在当前迭代t下的第二最终对偶变量，

表示第j个计算节点在当前迭代t下的第一最终对偶变量，λ表示松弛系数，n_d表示计算节点的个数，n_b表示小波基的个数，上标T表示转置。

与现有技术相比，本发明具有如下特点：

1、将射电干涉阵列图像重建的逆问题转换为等式约束条件下的凸优化问题，再转换为求解该问题的对偶问题，使用原始对偶算法，每次迭代需要的总计算成本较低，从而降低了该问题的求解复杂性；

2、将数据分割成多份，能够精确快速的分割数据；并使用多计算节点进行计算，方法灵活且具有较高的并行性能，可拓展性强；

3、不需要任何的数据后处理步骤，直接得到重建图像，最终得到的图像的动态范围更大，对延展结构的天体显示更加清晰，相比传统的方法更加高效。

附图说明

图1为射电干涉阵列分布式图像重建方法总体框图。

图2为几种方法图像重建结果对比图。其中(a)本发明；(b)CLEAN；(c)CS_CLEAN；(d)MS_CLEAN。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实例，并参照附图，对本发明进一步详细说明。

本实例使用的原始数据是由Karl G.Jansky甚大阵观测超新星遗迹3C391得到的数据再使用通用天文应用软件CASA进行校准之后的数据，数据格式为Measurement Set，文件名为3c391_ctm_mosaic_spw0.ms。

首先，查看文件3c391_ctm_mosaic_spw0.ms中的数据，选择目标射电源视场，即视场ID编号为0；使用CASA软件将目标射电源视场相关的观测数据从文件3c391_ctm_mosaic_spw0.ms中分离并创建1个子文件3c391_field0.ms。

然后，读取子文件3c391_field0.ms获取观测频率通道、uvw数据、可见度数据的实部和虚部以及可见度数据对应的权重。

接着，根据观测频率通道和偏振划分uvw数据和可见度数据及权重，分离出1个偏振的uvw数据和可见度数据及权重，然后再分离出1个频率通道的UVW数据和可见度数据及权重。将这些分离后的可见度数据转换为复数形式后，与分离后的uvw数据和可见度权重一起按列方式保存到一个文本文件中。

最后，读取上述得到的文本文件uvw数据中的uv数据，计算出uv值的最大范围bmax即

的最大值，计算出uw和vw：vw＝v×π/(bmax×dl)，uw＝u×π/(bmax×dl)，其中dl是1个像素的大小；读取出文本文件中的可见度数据y_I及其权重weights，将y_I与weights相乘得到新的可见度数据y。

本发明所处理的数据是已经进行了校准的数据，因此不需要考虑方向依赖响应。

假设待重建的图像

(N维实数信号)在某个基(如小波基)Ψ中可以用系数向量

稀疏表示(

有且仅有D＜＜N个非零系数)，即：x＝Ψα,||α||₀＝D＜＜N。

本发明通过将由可见度测量值

(M维实数信号)重建图像

的病态逆问题：y＝Φx+n，转换为求解下式的凸优化最小化问题的解重建图像x：

其中，Ψ^T是Ψ的转置；Φ是测量矩阵(或测量算子)，其表达式为Φ＝GFZ；这里

表示内插核，用于去网格；F表示傅里叶变换算子；矩阵

表示补零或零填充，用于解释过采样和缩放图像以补偿不完全的内插；l表示字典Ψ下的稀疏先验

范数l＝||·||₁；函数f表示图像恢复解的真实性和正数性要求，表达式为：

h函数确保数据保真约束残差在由噪声水平∈定义的

范数球内，表达式为：

这里∈是残差的

范数的上限值。

矩阵(算子)

表示n_b个稀疏基的集合。

为了实现分布式(并行)计算本发明将原始数据分割成n_d个数据块，有：

其中，G由紧支撑核

和矩阵

用于选取出数据分块j部分的离散傅里叶平面。

因此，对于每个数据分块有：

y_j＝Φ_jx+n_j

其中，n_j是关于y_j的噪声。

所以本发明通过求解下式的最小化问题，即原始问题：

其中，l＝||·||₁，

∈_j表示第j个数据分块的噪声临界值；γ是额外调优参数；对于稀疏先验，

范数是可分离的稀疏基可以表示为：

其中，

i∈{1,...,n_b}；

是Ψ_i的转置。

本发明使用原始对偶算法寻找上述最小化问题的解，转换为求解原始问题的对偶问题，再利用对偶变量求解出原始变量，上述原始问题的对偶问题定义为：

其中，u_i和v_j是对偶变量；符号*表示取共轭。

一种基于原始对偶算法的射电干涉阵列分布式图像重建方法，具体的迭代步骤如下：

步骤1：中央节点初始化第1次迭代参数：初始图像x⁽¹⁾，第j数据块的对偶变量为

和

收敛常量

其中

表示第j数据块测量矩阵Φ_j的谱范数，τ＝0.49；

变量λ；将初始化的参数广播并发送给计算节点。稀疏基Ψ_i采用的是小波基。

在本实施例中，设迭代次数变量为t并初始化t＝1，初始化图像解x^(t)为512×512的全零矩阵；收敛常量

其中

表示第j数据块测量矩阵Φ_j的谱范数，τ＝0.49；

变量λ；根据uw和vw计算出栅格化矩阵G、傅里叶矩阵F和补零矩阵Z。

步骤2：中央节点根据计算节点个数对数据y进行分块，并将分块后的数据y_j按照j∈{1,2,...,n_d}的顺序发送到n_d个计算节点每个计算节点获得1份数据块。其中数据分割的具体方法如下：

(1)求出n_d的所有约数并按从小到大的顺序排列记为一维数组n_div；

(2)找出n_div中最先大于或等于

的数并记为v_no，令u_no＝n_d/v_no；

(3)计算uw的秩统计量R_uw，秩统计量进行归一化并量化，量化范围为：1～u_no，计算公式为：s_u＝ceil(u_no×R_uw/l_uw)，其中ceil是向右取整函数，l_uw为uw的长度；

(4)将s_u分别等于1,2,...,u_no的下标记录下来，uw、vw和y分别按照这些下标分割成u_no份；

(5)将(4)中分割的出的数据按照(3)的计算方法计算出(4)之后的vw的秩统计量，并进行归一化和量化得到s_v；

(6)将s_v分别等于1,2,...,v_no的下标记录下来，(4)分割之后的每份数据按照这些下标分割成v_no份，最终得到n_d份数据。

在本实施例中，将uw、vw和y分割成n_d份，每份数据长度尽量相等，并发送到n_d个计算节点；计算每个数据块的噪声上限

步骤3：n_d个计算节点同时进行计算对偶变量，各个计算节点得到的最终对偶变量

和

同时发送给中央节点。

①第j个计算节点对当前迭代t下的第一最终对偶变量

的计算过程为：

首先，计算当前迭代t下的第一初始对偶变量

然后，计算当前迭代t下的第一中间对偶变量

最后，计算当前迭代t下的第一最终对偶变量

其中，

表示第j个计算节点在上次迭代t-1下的第一中间对偶变量，Ψ表示小波基，Ψ_i表示小波基Ψ中的第i个基，

表示上次迭代t-1下的最终解，下标_κ表示收敛速度系数，

表示谱范数，λ表示松弛系数，i∈{1,2,...,n_b}，n_b表示小波基的个数，j∈{1,2,...,n_d}，n_d表示计算节点的个数。

②第j个计算节点对当前迭代t下的第二最终对偶变量

的计算过程为：

首先，计算当前迭代t下的第二初始对偶变量

然后，计算当前迭代t下的第二中间对偶变量

最后，计算当前迭代t下的第二最终对偶变量

其中，

表示第j个计算节点在上次迭代t-1下的第二中间对偶变量，G_j表示第j个计算节点的紧支撑核，M_j表示第j个计算节点的数据选取矩阵，F表示傅里叶变换矩阵，Z表示补零矩阵，

表示上一次迭代t-1下的最终解，y_j表第j个数据块的可见度数据，

表示第j个计算节点求解松弛变量的近似算子，λ表示松弛系数，上标T表示转置，j∈{1,2,...,n_d}，n_d表示计算节点的个数。

步骤4：中央节点计算最终解

首先，计算当前迭代t下的初始解

然后，计算当前迭代t下的中间解x^(t)：

最后，计算当前迭代t下的最终解

其中，P_C表示指示函数的近似算子，x^(t-1)表示上次迭代t-1下的中间解，τ表示解的收敛常量，Z表示补零矩阵，F表示傅里叶矩阵，

表示第j个计算节点计算节点的第二收敛常量，M_j表示第j个计算节点的数据选取矩阵，σ_j表示第j个计算节点的第一收敛常量，

表示第j个计算节点在当前迭代t下的第二最终对偶变量，

表示第j个计算节点在当前迭代t下的第一最终对偶变量，λ表示松弛系数，n_d表示计算节点的个数，n_b表示小波基的个数，上标T表示转置。

步骤5：重复步骤3和4，直至

或最后一次迭代结束，最后得到的

即为重建的图像。

图2分别为本发明算法与CLEAN、CS-CLEAN、MS-CLEAN算法得到的图像重建结果，迭代次数均为100次。这些图像的动态范围分别是：15266.0000、65.1900、167.9851和131.2030。可以看出本发明算法的得到图像动态范围最高，图像质量最高，CLEAN算法的动态范围最小，图像质量最低，说明了CLEAN算法对延展结构的天体重建效果很差，表明本发明算法图像重建效果最好。

需要说明的是，尽管以上本发明所述的实施例是说明性的，但这并非是对本发明的限制，因此本发明并不局限于上述具体实施方式中。在不脱离本发明原理的情况下，凡是本领域技术人员在本发明的启示下获得的其它实施方式，均视为在本发明的保护之内。

Claims (4)

1.射电干涉阵列分布式图像重建方法，其特征是，包括步骤如下：

步骤1、中央节点将初始化参数广播并发送给计算节点；

步骤2、中央节点根据计算节点的个数n_d，将射电干涉仪的测量数据分割为n_d份；并将这n_d份数据块分别发送给n_d个计算节点；

步骤2.1、求出计算节点个数n_d的所有约数，并将这些约数按从小到大的顺序排列为一维数组n_div；

步骤2.2、找出一维数组n_div中最先大于或等于

的数v_no，并令u_no＝n_d/v_no；

步骤2.3、对测量数据的uw投影的秩统计量R_uw进行归一化并量化后，得到uw投影的量化数据s_u：

s_u＝ceil(u_no×R_uw/l_uw)

其中，ceil表示向右取整函数，R_uw表示测量数据的uw投影的秩统计量，l_uw表示测量数据的uw投影的长度；

步骤2.4、将uw投影的量化数据s_u分别等于1,2,...,u_no的下标记录下来，并根据这些下标将测量数据的投影数据uw和vw，以及可见度数据y均分割成u_no份；

步骤2.5、对测量数据的vw投影的秩统计量R_vw进行归一化并量化后，得到vw投影的量化数据s_v：

s_v＝ceil(v_no×R_vw/l_vw)

其中，ceil表示向右取整函数，R_vw表示测量数据的vw投影的秩统计量，l_vw表示测量数据的vw投影的长度；

步骤2.6、将vw投影的量化数据s_v分别等于1,2,...,v_no的下标记录下来，并根据这些下标将步骤2.4所得的被分割成u_no份的测量数据的投影数据uw和vw，以及可见度数据y均再进一步分割成v_no份，最终得到n_d份数据块；

步骤3、每个计算节点根据所得到的数据块计算当前迭代t下的2个最终对偶变量

和

并将所计算出的2个最终对偶变量

和

同时发送给中央节点；

步骤4、中央节点根据各个计算节点所返回的2个最终对偶变量

和

计算当前迭代t下的最终解

步骤5、判断是否满足收敛条件，即当前迭代t达到最大迭代次数或者当前迭代t下的最终解

使得每个数据块均满足

如果未满足收敛条件，则迭代次数t+1，并返回步骤3；

如果满足收敛条件，则当前迭代t下的最终解

即为重建的图像；

其中，Φ_j表示第j个数据块的测量矩阵，y_j表示第j个数据块的可见度数据，∈_j表示第j个数据块的噪声上限，j∈{1,2,...,n_d}，n_d表示计算节点的个数。

2.根据权利要求1所述射电干涉阵列分布式图像重建方法，其特征是，步骤3中，第j个计算节点对当前迭代t下的第一最终对偶变量

的计算过程为：

首先，计算当前迭代t下的第一初始对偶变量

然后，计算当前迭代t下的第一中间对偶变量

最后，计算当前迭代t下的第一最终对偶变量

其中，

表示第j个计算节点在上次迭代t-1下的第一中间对偶变量，Ψ表示小波基，Ψ_i表示小波基Ψ中的第i个基，

表示上次迭代t-1下的最终解，下标κ表示收敛速度系数，

表示谱范数，λ表示松弛系数，i∈{1,2,...,n_b}，n_b表示小波基的个数，j∈{1,2,...,n_d}，n_d表示计算节点的个数。

3.根据权利要求1所述射电干涉阵列分布式图像重建方法，其特征是，步骤3中，第j个计算节点对当前迭代t下的第二最终对偶变量

的计算过程为：

首先，计算当前迭代t下的第二初始对偶变量

然后，计算当前迭代t下的第二中间对偶变量

最后，计算当前迭代t下的第二最终对偶变量

其中，

表示第j个计算节点在上次迭代t-1下的第二中间对偶变量，G_j表示第j个计算节点的紧支撑核，M_j表示第j个计算节点的数据选取矩阵，F表示傅里叶变换矩阵，Z表示补零矩阵，

表示上一次迭代t-1下的最终解，y_j表第j个数据块的可见度数据，

表示第j个计算节点求解松弛变量的近似算子，λ表示松弛系数，上标T表示转置，j∈{1,2,...,n_d}，n_d表示计算节点的个数。

4.根据权利要求1～3中任意一项所述射电干涉阵列分布式图像重建方法，其特征是，步骤4中，中央节点对当前迭代t下的最终解

的计算过程为：

首先，计算当前迭代t下的初始解

然后，计算当前迭代t下的中间解x^(t)：

最后，计算当前迭代t下的最终解

其中，P_C表示指示函数的近似算子，x^(t-1)表示上次迭代t-1下的中间解，τ表示解的收敛常量，Z表示补零矩阵，F表示傅里叶矩阵，

表示第j个计算节点计算节点的第二收敛常量，M_j表示第j个计算节点的数据选取矩阵，σ_j表示第j个计算节点的第一收敛常量，

表示第j个计算节点在当前迭代t下的第二最终对偶变量，

表示第j个计算节点在当前迭代t下的第一最终对偶变量，λ表示松弛系数，n_d表示计算节点的个数，n_b表示小波基的个数，上标T表示转置。

Images