基于压缩感知法对法拉第色散函数的重构*

于衍川，孙晓辉

(云南大学物理与天文学院，云南 昆明 650500)

法拉第旋转(Faraday Rotation)是线偏振电磁波通过磁化星际介质时偏振角发生旋转的一种物理现

(1)

其中，ne为热电子密度；B‖为视线方向的磁场分量；dl是视线方向上的线元。这样观测到的偏振量P可以写成[1]：

P(λ2)=Q(λ2)+iU(λ2)

(2)

(3)

其中，F(φ)为法拉第色散函数(Faraday Dispersion Function)，是法拉第深度的函数；Q和U为斯托克斯参量，是观测波长(频率)的函数。从(3)式可以看出，偏振量P和法拉第色散函数F(φ)之间是傅里叶变换对。F(φ)表示法拉第深度为φ处的偏振强度和偏振位置角，反映辐射区域以及辐射传播途径上的磁场结构。为了研究星际介质的磁场，需要重构法拉第色散函数。目前已有多种重构法拉第色散函数的方法。文[1]提出了法拉第旋率综合方法(Faraday Rotation Measure Synthesis, RMS)，通过观测偏振量和法拉第色散函数之间的傅里叶关系重构法拉第色散函数。因为观测波段是有限的，所以获得的法拉第色散函数是真实色散函数与窗口函数的卷积，对复杂的源需要解卷积。与法拉第旋率综合方法类似的还有小波方法[2-3]。此外，还有一些假定辐射模型的方法重构法拉第色散函数，如Q-U拟合法[4]。

本文主要讨论压缩感知方法(或称为压缩采样，Compressive Sensing/Sampling，原理见[5-6])。该方法的优点在于利用极少的采样数据也有可能重构原始的完整信号。压缩感知理论以信号的稀疏性为前提，通过编码测量和重构算法实现信号的精确重构。目前，已经提出了基于压缩感知法的法拉第色散函数重构算法，如文[7]采用最小l1范数法建立了适用于法拉第薄源(Faraday Thin，λ2Δφ≪ 1，这里Δφ为源在法拉第深度域上的延展范围，即F(φ)的宽度)、法拉第厚源(Faraday Thick,λ2Δφ≫ 1)和同时含有法拉第薄源和法拉第厚源的算法，模拟结果显示，该方法可以用于重构法拉第色散函数，并与法拉弟旋率综合方法进行比较，发现用压缩感知方法重构的法拉第色散函数在数值和误差上要优于法拉弟旋率综合方法。文[8]采用匹配追踪算法建立了对同时含有法拉第薄源和法拉第厚源进行重构的方法。但这些方法目前还处于理论模拟阶段，其实用性还有待研究，比如文[7]的算法中设置的频率(波长)范围过大，未将噪声添加到实验数据中，算法中各种参数设置对重构结果的影响也没有进行实验和讨论以及实验样本过少等。本文基于文[7]的方法，从实际观测角度进行数值模拟，对以上问题进行更深入研究，以此探讨其实用性。这里假定法拉第色散函数是稀疏的(即针对源包含一个或多个法拉第薄成分)，用Python语言编写了重构法拉第色散函数的程序。采用的频率范围为1.1～3.1 GHz，对应澳大利亚望远镜致密阵(Australia Telescope Compact Array, ATCT)上L波段的观测频率范围[9]。

1 重构方法

已知测量矩阵A∈RM × N(M≪N)和未知信号x0在该矩阵的线性投影y∈RM，即y=Ax0，要由测量值y重构未知信号x0。在压缩感知理论中，假定未知信号x0为稀疏的，信号x0可以由测量值y通过求解最小l1范数问题精确重构：

x=min‖x0‖l1s.t.Ax0=y.

(4)

若未知信号x0为非稀疏的，可以通过某种变换(如小波变换)进行稀疏性表示，即x0=ωX，X为该信号在ω变换域下的稀疏性表示，则有

x=min‖X‖l1s.t.Ax0=AωX=y.

(5)

将(2)式写成矩阵形式：

Yf=P，

(6)

其中，f为法拉第色散函数，是一个N维矢量；Y为P和f之间傅里叶变换的离散表示，是一个M×N矩阵；P为观测偏振量，是一个M维矢量，这里M为观测频率通道数，每个通道上接收不同波段的偏振量。

由于法拉第色散函数f是复变函数，可以将其实部和虚部分别写成f.real和f.imag，则(6)式可变为

(7)

(8)

(9)

实验中，将均值为0、标准差σ为1(强度单位任意，下同)的高斯噪声加到Q和U上，并设置在法拉第深度φ空间的分辨率δφ为ΔφFWHM的1/4，则总的格点数N(φ空间中坐标点φi的个数，i=1, 2, …,N)为

(10)

其中，int()表示取整。表1列出了一些主要的实验参数。

表1 实验参数Table1 Experimental parameters

2 结果与讨论

2.1 一般情况

在模拟实验中，假定一个法拉第色散函数包含5个法拉第薄源,这些源分别是：F(-164.1)=3 + 6i，F(-82.0)=-4 + 3i，F(0.0)=2 + 3i，F(68.4)=-8 - 5i，F(136.7)=5 + 3i，对应的信噪比分别为：6.7， 5.0， 3.6， 9.4， 5.8。对包含这5个法拉第薄源的法拉第色散函数进行了反解，结果如图1(a)，

(11)

其中，p为偏振强度。从图1可以看出，没有噪声时该模型准确地重构了原始图像，而加入噪声时，重构结果出现了偏差，|dwm|分别为5.80 × 10-5，11.78。

通过压缩感知方法重构的法拉第色散函数随着法拉第薄源偏振强度的大小、相位以及两个法拉第薄源的法拉第深度之间间隔Δφ的不同而不同，甚至不能准确地重构出原始图像。在这里分两种情况进行测试说明：(1)假定上述5个法拉第薄源的法拉第深度之间的间隔Δφ不变，改变其大小和相位，如F(-164.1)=3 - 6i，F(-82.0)=-4 + 3i，F(0.0)=5 + 3i，F(68.4)=8 - 5i,F(136.7)=-5 + 4i，对应的信噪比分别为：6.7，5.0，5.8，9.4，6.4，重构结果如图1(b)；(2)假定大小和相位不变，改变其法拉第深度之间的间隔Δφ，如：F(-150.4)=3 + 6i，F(-68.4)=-4 + 3i，F(0.0)=2 + 3i，F(109.4)=-8 - 5i，F(164.1)=5 + 3i，对应的信噪比分别为6.7，5.0，3.6，9.4，5.8，重构结果如图1(c)。

当假定间隔Δφ不变时，从图1可以看出，前3个源并没有很好地重构出原始图像，其相位和大小存在一定误差，其中，|dwm|分别为6.31，17.03。而假定大小和相位不变时，后面4个源都与原始图像出现了较大偏差，其中，|dwm|分别为8.18，17.11。同时也发现，在不加入噪声时，将这些法拉第薄源的法拉第深度之间的间隔Δφ设置得足够大时，即Δφ> 90 rad·m-2时，无论改变其大小、相位还是间隔Δφ，该模型均能准确地重构原始图像，即使在加入噪声时，产生的误差也很微小。当然，以上测试的这些源的法拉第深度都是设置在该方法的格点N上进行的。

通常，在实际应用中，事先并不知道法拉第色散函数的法拉第深度以及两个源的法拉第深度之间的间隔Δφ，要重构的法拉第色散函数各个成分的法拉第深度是随机的，一般不在这种方法设置的格点N上，并且Δφ也是随机的，所以也无法确定该模型是否能够很好地重构法拉第色散函数。随机选取一个包含5个法拉第薄源的色散函数(这些法拉第薄源的法拉第深度不在格点N上)：F(-114)=-6 - 1i，F(-60)=7 - 4i，F(48)=3 - 5i，F(102)=-3.5 + 3i，F(210)=9 + 3i，对应的信噪比分别为：6.1，8.1，5.8，4.6，9.5。其测试结果如图1(d)， 其中，|dwm|分别为0.40，10.20。

每次模拟都是随机产生的高斯噪声，所以每次重构的结果不同。在实验中，对于法拉第色散函数的重构结果随着源的个数、大小、相位和噪声的不同出现不同程度的误差，同时也与该方法中设置的参数有关，如观测通道M，φ空间的分辨率δφ。

图1 法拉第色散函数的重构结果。每幅图从上往下，上图为原始数据，中图为无噪声重构结果，下图为有噪声重构结果，粗实线表示振幅，细实线表示实部，细虚线表示虚部

Fig.1 The reconstruction result of Faradqy dispersion function. From top to bottom in each panel: the original data, the noiseless reconstruction result and the with-noise reconstruction result. The thick solid line shows the amplitude, the thin solid line as real part, the dashed line as imaginary part

2.2 大样本实验

为了更直观地观察该方法的重构效果，从噪声、两个源的法拉第深度之间的间隔Δφ、φ空间的分辨率δφ(δφ的不同会导致总的格点数N不同)3方面对重构结果的影响分别进行了大样本统计学实验。以重构的源的数目和φ的加权平均的差值描述重构结果的好坏[10]。

首先考察信噪比对重构结果的影响，这里假定法拉第色散函数只包含一个法拉第薄源，固定其法拉第深度为φ=20 rad·m-2，随机模拟了1 000个信噪比范围在(3, 60)之间的源(这里噪声依然设置为1，信噪比相当于法拉第色散函数的偏振强度)，φ空间的分辨率δφ依然为ΔφFWHM的1/4，由于只含有一个法拉第薄源，直接以φ的差值|φrecorved-φoriginal|描述重构效果，并排除信噪比小于2的信号。结果如图2(a)。从图中可以看出，对于含有一个源的情况，可以很准确地重构出源的数目，而φ的差值总是在两个值上，即重构的φrecorved出现在相邻的两个格点上。为了减小重构误差，对重构结果在峰值附近进行最小二乘拟合，结果如图2(b)。经过拟合后，对于大多数源来说，φ的差值缩小为原来的1/2～1/4。由于重构的法拉第深度φ总在相邻的两个格点上，是否也可以通过缩小φ空间分辨率的方式减小φ的重构误差呢？

图2 从左到右为 (a) 无拟合；(b) 拟合后。每幅图从上往下，上为重构数目，下为φ的差值

Fig.2 From left to right: (a) no fitting; (b) fitting. From top to bottom in each panel:the number of reconstructions, the difference ofφ

下面考查两个源的法拉第深度之间的间隔Δφ和φ空间的分辨率δφ对重构结果的影响。假定一个法拉第色散函数包含2个法拉第薄源，固定其中一个源的偏振强度p=13.235，偏振位置角χ=2.749，法拉第深度φ=10.256 rad·m-2。第2个源是随机模拟了1 000个偏振强度在(3, 30)、偏振位置角在(0, 2π)、间隔Δφ在(5, 200)之间的源。同时，对不同分辨率δφ的情况(即δφ取ΔφFWHM的1/2，1/3，1/4，1/5，1/6，1/7，1/8，1/9)分别进行以上实验，这里排除信噪比小于1的信号。结果如图3。

为了比较不同分辨率情况下的重构效果，对不同分辨率下准确重构出源的数目的精确度进行了比较，如图4。从图4可以看出，当δφ取ΔφFWHM的1/3和1/4时，精确度最高，达到83%左右。在以上对噪声影响的实验时，曾提出是否可以通过缩小分辨率的方式减小重构误差的问题，而在这个实验中可以看出，分辨率δφ并不是越小越好。

图3 不同分辨率下间隔Δφ对重构结果的影响。第1排从左到右为(a)δφ取ΔφFWHM的1/2；(b)δφ取ΔφFWHM的1/3；第2排从左到右为 (c)δφ取ΔφFWHM的1/4；(d)δφ取ΔφFWHM的1/5；第3排从左到右为 (e)δφ取ΔφFWHM的1/6；(f)δφ取ΔφFWHM的1/7；第4排从左到右为(g)δφ取ΔφFWHM的1/8；(h)δφ取ΔφFWHM的1/9。每幅图从上往下，上为重构数目，下为φ的加权平均的差值dwm

Fig.3 Influence of the interval Δφat different resolutions. From left to right in first row: (a)δφtakes 1/2 of ΔφFWHM; (b)δφtakes 1/3 of ΔφFWHM, in second row: (c)δφtakes 1/4 of ΔφFWHM; (d)δφtakes 1/5 of ΔφFWHM, in third row: (e)δφtakes 1/6 of ΔφFWHM; (f)δφtakes 1/7 of ΔφFWHM, in fourth row: (g)δφtakes 1/8 of ΔφFWHM; (h)δφtakes 1/9 of ΔφFWHM. From top to bottom in each panel: the number of reconstructions, the difference (dwm) of the weighted average ofφ

图4 不同分辨率下的重构精度比较。x轴表示δφ/ΔφFWHM的比值，y轴表示重构数目的精确度

Fig.4 Comparison of reconstruction accuracy at different resolutions.X-axis represents the ratio ofδφ/ΔφFWHM,andY-axis represents the accuracy of the number of reconstructions

图5 以δφ取ΔφFWHM的1/4时实验为例，随机选取的重构图像
Fig.5 Examples of reconstructed images, whenδφtakes 1/4 of ΔφFWHM

图6dwm随着间隔Δφ的变化
Fig.6 This picture showsdwmvaries with Δφ

3 结 论

在多个频率上测量弥漫射电辐射的偏振是研究银河系磁场的一种有效的方法。目前的法拉第旋率综合法是利用偏振强度与法拉第色散函数之间是傅里叶变换对，把观测到的偏振强度做傅里叶逆变换得到法拉第色散函数[1]。因受观测波段的限制，得出的解需要进行解卷积，所以需要一种更强有力的方法解决这一问题。

本文基于压缩感知的法拉第色散函数重构方法应用到实际观测频率范围中，并对可能影响该方法重构效果的因素分别进行了实验。实验结果表明，该方法的重构效果随着法拉第薄源的个数、源偏振强度的大小、相位以及两个法拉第薄源的法拉第深度之间间隔Δφ的不同而出现较大差异，具有较大的随机性。大样本统计学实验发现，当φ空间的分辨率δφ取ΔφFWHM的1/3和1/4时，重构结果精确度最高，对其重构结果各峰值附近再次进行最小二乘拟合后得到的法拉第深度φ值更接近原始值。通过压缩感知法重构后的再次拟合，可以得到更为精确的法拉第深度φ值，因此该方法可以作为现有方法的一种补充，同时该方法对于重构得出的偏振强度普遍偏低，重构效果具有较大的随机性，所以，还需要不断地改进，提高其稳定性和精确度，以便更好地应用于实际观测数据中。