现有的大部分暂态源检测方法都是从像素域中连续帧之间的变化来检测瞬变源，而像素域检测方法需要先得到重建图像。基于目标源在被观测时间段内是静态的假设，传统综合孔径成像方法借助地球自转累积变化的基线获得足够的采样覆盖，来进行亮温图像重建[4]。但对于存在瞬变源的场景，该假设不再成立。例如，银心超大黑洞Sgr A*最内层轨道变化周期在分钟量级[5]，传统方法重建出来的图像中源的变化过程被模糊，丢失了时变信息。另外，单独对每帧可见度数据进行重建，虽然保持了时域上变化信息的独立性，但每帧包含的可见度函数覆盖有限，由此得到的重建图像中存在混叠。
对于从可见度函数采样数据中重构出一段时间的变化亮温图像问题，已有文献提出了解决思路，主要包括逐帧重建方法、时域差分约束重建方法和稀疏重建方法等。逐帧重建方法将各帧图像分开重建，或在减去时间上的平均值后再逐帧重建[6]，或假设变化源是时域光滑的，将时域变化展开为泰勒多项式[7]或傅里叶多项式[8]。这些方法可解决瞬变源变化缓慢的重建问题。但很多情况下，图像的大部分区域不变，只有几个局部区域变化较大。这种情况下，对瞬变源进行光滑近似的方法存在2个问题：①重建结果不能很好地表示变化较大的瞬变源；②不变的静态区域在重建图中出现了缓慢的变化。时域差分约束重建方法通过对每帧图像和前后帧差分分别添加?1和?2约束，进行变化源的重建[9]。但该方法中时间上的前后帧之间并未关联起来，且没有考虑前后帧差分本身的稀疏性，并不适用于变基线情况下的瞬变源重建。Wenger等[10]提出的组稀疏暂态源成像算法属于稀疏重建方法，通过使用单个像素的时域信号满足组稀疏性约束在压缩感知框架下通过?1, ∞范数正则优化重建图像。由于组稀疏性的限制，该算法只能处理背景和瞬变源均在像素域稀疏的场景。
因此，本文提出了一种基于背景和目标变化的直和空间上稀疏性的稀疏基线综合孔径动态场景图像重建算法，解决不同变化和背景下瞬变源的综合孔径图像重建问题。
1 干涉成像原理 综合孔径射电干涉测量本质上是通过对空间中不同点收到的信号进行互相关，测出信号源空间分布的傅里叶分量，使用综合孔径成像算法，通过对干涉处理得到的可视度函数的逆傅里叶变换反演出信号的空间分布，得到亮温分布图。
如图 1所示，测量傅里叶分量的基本单元是二元干涉仪。二元干涉仪的功能是将2个天线接收到的信号在保持其原始相位的情况下进行复相乘，不同基线长度的二元干涉仪便可测量出空间频域不同采样点的幅度值。布置多个天线，每2条天线的输出信号都互相干涉，得到多条干涉基线测量的空间频率信息集合，称为可见度函数，其表达式为
 | (1) |
 |
| 图 1 二元干涉仪原理图 Fig. 1 Principle diagram of binary interferometer |
| 图选项 |
式中：V为可见度函数；u和v分别为两天线单元以波长为单位的距离向量的2个分量；ξ和η分别为观测方向余弦；I(ξ, η)为辐射亮温分布。
通过对可见度函数进行完整测量后，作逆傅里叶变换便可得到亮温图像。干涉测量得到的是空间频率域内的一系列离散点，根据傅里叶变换原理，采样点的最大间距决定了观测图像的分辨率，采样点的最小间距决定了观测的无混叠视场。
实际应用中，经常出现天线阵列对可见度函数采样覆盖不完全的情况，造成反演图像中出现混叠伪像。稀疏理论证明，一个非带限信号只要在某组稀疏基下具有稀疏性，即在该稀疏基下信号的大多数系数为零，就可由另一组稀疏基下的少量测量精确重构出来。对欠采样下综合孔径干涉测量过程进行建模，设x为视场内辐射亮温分布I的矢量形式，v为可见度函数V的矢量形式，M为部分傅里叶变换矩阵，即x为原始信号，v为观测矢量，M为观测矩阵，则测量方程可写为
 | (2) |
若存在一组稀疏基B满足：
 | (3) |
测量方程等价为
 | (4) |
使得x在B下的系数s具有稀疏性，通过寻找满足式(4)的最稀疏矢量s，可将x重构出[11-13]。
2 背景和变化的稀疏性 本节构造了一个描述瞬变源静态背景和变化直和空间上的稀疏域变换表示。考虑到典型的包含瞬变的射电图像中大多数像素很少或没有时域变化，可变像素是稀疏的。可将N帧图像序列表示为首帧和N－1个连续帧间差分：
 | (5) |
式中：xi为第i帧图像；di为第i+1帧与第i帧之间的差分。
式(5)可写为
 | (6) |
根据上述分析，首帧图像与帧间差分图像分别在小波域和像素域不同域上稀疏。图 2显示了大小为32×32、背景强度为0.1，包含变化点目标的两帧图所分别对应的像素基系数和小波基系数分布。首帧图像包含静态背景，像素域零值系数频次为0，而小波域零值系数频次很大。帧间差分图像由瞬变源构成，这些瞬变源通常是点源，其像素域零值系数高于小波域。与之对应，表 1列出了非零元素的频次即?0范数。首帧图像小波基系数?0范数||B1Tx10||为338，远小于像素基系数?0范数1 024，说明首帧图像在小波域下更稀疏。帧间差分图像像素基系数?0范数||B2Tdi0||为30，远小于小波基系数?0范数208，说明帧间差分图像在像素域下更稀疏。其中，B1表示小波基变换矩阵，B2表示像素基变换矩阵。
 |
| 图 2 首帧图像和帧间差分图像的像素基系数、小波基系数的分布 Fig. 2 Distribution of coefficients of the first frame and difference in pixel and wavelet domainss |
| 图选项 |
表 1 首帧图像和帧间差分图像的像素基系数、小波基系数的?0范数 Table 1 ?0 norm of coefficients of the first frame and difference in pixel and wavelet domains
| 图像 | 像素基系数?0范数 | 小波基系数?0范数 |
| 首帧图像 | 1024 | 338 |
| 第2帧图像与首帧图像差分图像 | 30 | 208 |
表选项
除首帧可见度函数只包含首帧图像亮温信息外，其他帧可见度函数同时包含首帧和帧间差分的亮温信息。这使得各帧可见度函数的测量方程不独立，无法表示为首帧和帧间差分的单独方程。首帧和帧间差分耦合在多帧测量方程中可写为直和的形式。首帧和帧间差分分别在小波域和像素域上稀疏，则直和空间上的?0范数如下：
 | (7) |
3 动态场景图像重建算法 设N帧观测数据对应的可见度函数采样矩阵分别为S1, S2, …, SN，傅里叶变换矩阵为F，可见度数据分别为v1, v2, …, vN，则观测方程可写为
 | (8) |
将式(8)中的直和运算表示为矩阵和矢量运算，即
 | (9) |
把式(6)代入式(9), 并写成如下矩阵形式：
 | (10) |
观测矩阵M为
 | (11) |
写成直和形式为
 | (12) |
图像重建过程即通过观测方程(12)求解直和矢量x1⊕d1⊕…⊕dN－1的过程。空间频率域采样不足时，v1⊕v2⊕…⊕vN欠采样，M为非满秩矩阵，式(12)为病态方程，x1⊕d1⊕…⊕dN－1的解不唯一。常用求解方法是正则化，通过将病态方程转换为最小平方优化问题，并增加约束条件，得到欠定方程的唯一解。本文引入合理性度量函数f(s)为约束条件，将x求解过程表达为目标函数的最小化：
 | (13) |
式中：第1项为数据项，保证由解s测量得到的可见度函数与实际测量结果v1⊕v2⊕…⊕vN一致；第2项为正则约束项，保证重构图像的合理性；λ为权重系数，平衡测量数据的真实性和图像的合理性。
f(s)定义为s的?1范数：
 | (14) |
f(s)越小，则x1⊕d1⊕…⊕dN－1在基B下的系数s越稀疏，图像越合理。
式(13)的求解可转换为迭代求解如下优化问题sk, k=0, 1, …[14]。
 | (15) |
式中:z和αk分别为迭代过程中的中间变量和权重更新比例；uk表示如下：
 | (16) |
给定稀疏基B1和B2，B1将首帧x1转换为稀疏矢量，B2将帧间差分d2转换为稀疏矢量。例如，最小化小波基非零系数个数要求临近像素具有相近强度，适用于扩展源，最小化像素基非零稀疏个数要求最小化非零像素个数，适用于紧致源[15]。
重建算法流程如下：
s=BTMTv1⊕v2⊕…⊕vN(初始化稀疏域系数)
α=1(初始化步长)
repeat
r=v1⊕v2⊕…⊕vN-MBs(更新剩余可见度函数)
λ=BTMTr∞(更新正则项系数)
repeat (最小化式(14))
r=v1⊕v2⊕…⊕vN-MB(s′-s)(更新剩余可见度函数)
s=s′(更新稀疏域系数)
until e≤ε
until r2≤δ
合理性度量项的权重系数
4 实验和结果分析 本节通过仿真实验验证了3种不同场景下算法的有效性，对比不同方法的重构准确性。通过相对均方根误差(Root Mean Square Error, RMSE)定量评估重建图像的质量，并通过多次随机实验的统计结果来验证本文算法的性能。对比方法中，静态重建方法从全部帧可见度数据中恢复一幅图像，?1方法分别在每帧可见度数据上采用?1范数正则优化得到该帧的重建图像，?2方法在?1方法的基础上增加了?2正则项来保持时域光滑性，?1, ∞方法在全部帧可见度数据上采用?1, ∞范数正则优化重建出各帧图像。
实验1仿真在小区域背景下瞬变源，见图 3中第2列。通过在变化的随机分布基线上，仿真连续16帧不同的32像素×32像素的综合孔径图像，获得每帧66个可见度数据(见图 3中第1列)。使用5种不同方法重构得到的图像分别对这些可见度数据进行重建，重建图像见图 3中第3~7列。所有方法均选像素基为稀疏基，λ=0.1，迭代10 000次。静态重建方法从全部可见度函数恢复一张图片忽略了不同时刻可见度函数中的时变信息，所以可恢复小区域背景，但丢失了点源的运动信息。?1方法由于使用的单帧可见度数据不足以恢复图像，导致重图中出现一些伪像。?2方法虽然可重建出光滑的时域变化，但正则项中只约束了时域变化的光滑性而没有约束其稀疏性，使得重建图像中静态背景部分也出现了变化。diff-?1算法的重建图像质量最高。
 |
| 图 3 不同方法对小区域背景下瞬变源的重建结果 Fig. 3 Reconstruction results of a transient source on a small region of background by different methods |
| 图选项 |
实验1中不同方法的相对均方根误差随迭代次数的变化曲线如图 4所示。可以看到，随着迭代的进行，所有方法的相对均方根误差都收敛到一个较低的值，但diff-?1算法具有最小的收敛值，其他方法中静态重建方法具有更小的收敛值。由于?1方法和?2方法只使用单帧可见度函数采样，而静态重建方法和diff-?1算法使用多帧的可见度函数采样，包含了更多信息，重构结果具有更小的相对均方根误差。diff-?1算法除了使用多帧采样数据外，还建立了帧间亮温变化和可见度函数之间的直接测量方程，并对变化添加了稀疏约束，同时利用了多帧可见度函数的信息和差分的稀疏性，所以其具有最低的相对均方根误差。
 |
| 图 4 小区域背景下，迭代次数增加时不同方法对瞬变源重建结果相对均方根误差的变化 Fig. 4 RMSE of reconstruction results of transient source on a small region of background by different methods as number of iteration increasing |
| 图选项 |
在实验1中将瞬变源强度设为随时间变化，通过在每个时刻选择一个0~1之间的随机值进行仿真，得到不同方法的相对均方根误差随迭代次数的变化曲线，如图 5所示。与图 4对比可以看到，瞬变源运动的同时强度变化，对diff-?1算法的重建效果影响不是很大。
 |
| 图 5 小区域背景下，迭代次数增加时不同方法对运动且强度变化的瞬变源重建结果相对均方根误差的变化 Fig. 5 RMSE of reconstruction results of transient source which is moving and varying in intensity on a small region of background by different methods as number of iteration increasing |
| 图选项 |
实验2使用与实验1相同的观测基线，综合孔径图中包含静态、线性光滑变化、不规则变化3种源。图像具有暗背景，所有瞬变源随机分布在不同位置，其中30个静态源强度为1，15个线性光滑变化源强度从1变化到16，15个不规则变化源强度在0 ~2之间随机变化。不同方法的重建结果如图 6所示。不同方法使用与实验1相同的参数。静态重建方法重建图中瞬变源的变化信息由于平均而丢失，?1方法由于基线覆盖不足，重建图像中丢失了大部分瞬变源。?2方法重建图像中只恢复了均匀瞬变源，静态和不规则瞬变源完全丢失。diff-?1算法的重建图像几乎无法与原图没有差异。
 |
| 图 6 不同方法对暗背景下瞬变源的重建结果 Fig. 6 Reconstruction results of transient source with dark background by different methods |
| 图选项 |
对不同的每帧可见度函数采样数，重复20次实验2，使用相对均方根误差的中位数作为重建质量的评价指标，对比了不同方法在不同每帧可见度函数采样数下的重建性能。如图 7所示，随着每帧可见度函数采样数的增加，不同方法重建结果的相对均方根误差先增大，然后逐渐减小并收敛。静态重建方法下降速度最快，但每帧可见度函数采样数达到65后重建质量最差，超过20后，每帧可见度函数采样数的增加对图像重建质量提升效果不大。?1方法的重建质量随着每帧可见度函数采样数的增加持续下降，这是因为?1只使用了单帧可见度数据，每帧可见度函数采样数越多重建效果越好。?2方法介于静态重建方法和?1方法之间，在每帧可见度函数采样数较少时相对均方根误差下降较快，并在每帧可见度函数采样数继续增大时持续下降，这是因为?2方法包含时域上的依赖，可重建出部分变化信息，且随着每帧可见度函数采样数的增加，重建的变化信息更准确。随着每帧可见度函数采样数的增加，?1, ∞方法的下降趋势与?2方法相似，但?1, ∞方法具有更小的相对均方根误差，这是因为?1, ∞方法中的组稀疏约束比?2方法的时域光滑约束能更准确地描述变化信息。diff-?1算法具有最小的相对均方根误差，这是因为diff-?1算法显式要求瞬变源的时变分分布具有稀疏性，更符合瞬变源的变化特征。
 |
| 图 7 暗背景下，每帧可见度函数采样数增加时不同方法对瞬变源重建结果相对均方根误差的变化 Fig. 7 RMSE of reconstruction results of transient source with dark background by different methods as number of sampling of visibility increasing |
| 图选项 |
实验3的实验条件与实验2基本相同，唯一的变化是图像加上了强度为0.1的静态背景。由于像素域不再稀疏，不同方法使用小波基作为稀疏基。不同方法的重建结果如图 8所示。可以看到，除了diff-?1算法可准确重建出目标信息，其他方法重建图像均失真严重。这是因为图像中背景虽是小波域稀疏的，但其中的瞬变源却是像素域稀疏的，diff-?1算法通过在首帧和差分上分别使用小波基和像素基作为稀疏基，能准确地描述这种特征，从而得到无失真的重建图像。?1方法使用的数据不足，导致其重建图中静态背景失真和瞬变源缺失。?2方法对变化约束不强，导致重建图像中背景也出现了变化。小波基下，?1, ∞方法的稀疏约束不能准确描述瞬变源变化，重建图像中出现部分瞬变源周围像素强度增加的情况，另一方面组稀疏约束导致重建图像中静态背景部分区域缺失。
 |
| 图 8 不同方法对大区域背景下瞬变源的重建结果 Fig. 8 Reconstruction results of transient source on a big region of background by different methods |
| 图选项 |
对实验3进行与实验2相同的重复实验，得到不同方法在不同每帧可见度函数采样数下的重建性能。从图 9中可以看到，所有方法重建结果的相对均方根误差都随着每帧可见度函数采样数的增加而减少。?1方法相对均方根误差最大，主要原因是单帧可见度数据不足够。静态重建方法、?2方法和?1, ∞方法在大于每帧60个可见度函数采样后相对均方根误差接近，在小于60的范围内静态重建方法下降最快。diff-?1算法具有最小的相对均方根误差。
 |
| 图 9 大区域背景下，每帧可见度函数采样数增加时不同方法对瞬变源重建结果相对均方根误差的变化 Fig. 9 RMSE of reconstruction results of transient source on a big region of background by different methods as number of sampling of visibility increasing |
| 图选项 |
在与实验2相同的场景下，每帧取固定的60个UV采样，在强度变化率分别为0, 0.5，…, 5情况下进行重复实验，得到diff-?1算法在不同强度变化率下的误差曲线。从图 10中可以看到，diff-?1算法重建结果的相对均方根误差都随着强度变化率的增大而增加。
 |
| 图 10 大区域背景下，强度变化率增加时diff-?1算法对瞬变源重建结果相对均方根误差的变化 Fig. 10 RMSE of reconstruction results of transient source on a big region of background by diff-?1 algorithm as intensity varying rate increasing |
| 图选项 |
在与实验2相同的场景下，每帧取固定的60个UV采样，在信噪比(Signal to Noise Rate, SNR)分别为1, 2, …, 11 dB情况下进行重复实验，得到diff-?1算法在不同信噪比下的相对均方根误差曲线。从图 11中可以看到，diff-?1算法重建结果的相对均方根误差都随着信噪比的增大而减小。
 |
| 图 11 大区域背景下，信噪比增加时diff-?1算法对瞬变源重建结果相对均方根误差的变化 Fig. 11 RMSE of reconstruction results of transient source on a big region of background by diff-?1 algorithm as SNR increasing |
| 图选项 |
图 12给出了射电望远镜VLA实测数据上不同方法对瞬变源的成像结果。使用已定标过的可见度数据[16]，场景中第6帧存在一个FRB，图像像素大小为128像素×128像素，单帧UV覆盖率约7.67%。原数据UV覆盖固定不随时间变化，各帧具有相同的UV覆盖，通过把该覆盖率为23%的UV覆盖均匀随机划分为10个部分，得到每帧覆盖率约2.3%的变化UV覆盖序列，并取对应的可见度数据序列，从而获得了基线变化的实测数据序列。使用不同方法分别对这些可见度数据进行重建，所有方法均选像素基为稀疏基，λ=0.1，迭代10 000次。图 12中第1列为UV覆盖，第2列为使用未划分前的完整UV测量所得到的重建图像，可以作为参考图像。第3~7列分别为静态重建方法、?1方法、?2方法、?1, ∞方法和diff-?1算法的重建结果。可以看到，静态重建方法从全部可见度函数恢复一张图片忽略了不同时刻可见度函数中的时变信息，完全丢失了瞬变源的信息。?1方法、?2方法、?1, ∞方法的重建图像中虽然得到了瞬变源但强度偏低，diff-?1算法的重建图像质量最高。
 |
| 图 12 不同方法对VLA的FRB测量数据的重建结果 Fig. 12 Reconstruction results of measurement data of FRB on VLA by different methods |
| 图选项 |
5 结论 针对时变稀疏基线下瞬变源重建图像中目标模糊和动态信息丢失问题，本文提出了一种基于背景和目标变化的直和空间上稀疏性的diff-?1算法，用于稀疏基线综合孔径瞬变源图像重建。实验结果表明，对于稀疏基线下的观测数据：
1) diff-?1算法可得到准确的重建图像和动态信息。
2) 对于小区域背景下的运动和瞬变源，diff-?1算法与已有方法效果相当。
3) 对于大区域背景下的运动和瞬变源，diff-?1算法优于已有方法。
参考文献
| [1] | LORIMER D R, KARASTERGIOU A, MCLAUGHLIN M A, et al. On the detectability of extragalactic fast radio transients[J]. Monthly Notices of the Royal Astronomical Society, 2013, 436(1): L5-L9. DOI:10.1093/mnrasl/slt098 |
| [2] | VAN HAARLEM M P, WISE M W, GUNST A W, et al. LOFAR:The low-frequency array[J]. Astronomy & Astrophysics, 2013, 556: A2. |
| [3] | ZARKA P, GIRARD J, TAGGER M, et al.FAR: The LOFAR super station project in Nan?ay[C]//SF2A-2012: Proceedings of the Annual Meeting of the French Society of Astronomy and Astrophysics, 2012: 687-694. |
| [4] | THOMPSON A R, MORAN J M. Interferometry and synthesis in radio astronomy[M]. Berlin: Springer, 2017. |
| [5] | JOHNSON M D, LOEB A, SHIOKAWA H, et al. Measuring the direction and angular velocity of a black hole accretion disk via lagged interferometric covariance[J]. The Astrophysical Journal, 2015, 813(2): 132. |
| [6] | LAW C J, BOWER G C, BURKE-SPOLAOR S, et al. Realfast:Real-time, commensal fast transient surveys with the very large array[J]. The Astrophysical Journal(Supplement Series), 2018, 236(1): 8. DOI:10.3847/1538-4365/aab77b |
| [7] | RAU U.Radio interferometric imaging of spatial structure that varies with time and frequency[C]//Image Reconstruction from Incomplete Data VII.Bellingham: SPIE, 2012: 1-12. |
| [8] | STEWART I M, FENECH D M, MUXLOW T W B. A multiple-beam CLEAN for imaging intra-day variable radio sources[J]. Astronomy & Astrophysics, 2011, 535: A81. |
| [9] | SWINBANK J D, STALEY T D, MOLENAAR G J, et al. The LOFAR transients pipeline[J]. Astronomy and Computing, 2015, 11: 25-48. DOI:10.1016/j.ascom.2015.03.002 |
| [10] | WENGER S, RAU U, MAGNOR M. A group sparsity imaging algorithm for transient radio sources[J]. Astronomy and Computing, 2013, 1: 40-45. DOI:10.1016/j.ascom.2013.02.002 |
| [11] | CANDèS E J, TAO T. Decoding by linear programming[J]. IEEE Transactions on Information Theory, 2005, 51(12): 4203-4215. DOI:10.1109/TIT.2005.858979 |
| [12] | CANDōS E J.Compressive sampling[C]//Proceedings of the International Congress of Mathematicians Madrid, 2006: 1433-1452. |
| [13] | DONOHO D L. Compressed sensing[J]. IEEE Transactions on Information Theory, 2006, 52(4): 1289-1306. DOI:10.1109/TIT.2006.871582 |
| [14] | WRIGHT S J, NOWAK R D, FIGUEIREDO M A. Sparse reconstruction by separable approximation[J]. IEEE Transactions on Signal Processing, 2009, 57(7): 2479-2493. DOI:10.1109/TSP.2009.2016892 |
| [15] | LUSTIG M, DONOHO D L, SANTOS J M, et al. Compressed sensing MRI[J]. IEEE Signal Processing Magazine, 2008, 25(2): 72-82. DOI:10.1109/MSP.2007.914728 |
| [16] | CHATTERJEE S, LAW C J, WHARTON R S, et al. A direct localization of a fast radio burst and its host[J]. Nature, 2017, 541(7635): 58-61. DOI:10.1038/nature20797 |
