基于壓縮感知法對(duì)法拉第色散函數(shù)的重構(gòu)*

于衍川，孫曉輝

(云南大學(xué)物理與天文學(xué)院，云南 昆明 650500)

法拉第旋轉(zhuǎn)(Faraday Rotation)是線偏振電磁波通過磁化星際介質(zhì)時(shí)偏振角發(fā)生旋轉(zhuǎn)的一種物理現(xiàn)

(1)

其中，ne為熱電子密度；B‖為視線方向的磁場分量；dl是視線方向上的線元。這樣觀測到的偏振量P可以寫成[1]：

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

(2)

(3)

其中，F(xiàn)(φ)為法拉第色散函數(shù)(Faraday Dispersion Function)，是法拉第深度的函數(shù)；Q和U為斯托克斯參量，是觀測波長(頻率)的函數(shù)。從(3)式可以看出，偏振量P和法拉第色散函數(shù)F(φ)之間是傅里葉變換對(duì)。F(φ)表示法拉第深度為φ處的偏振強(qiáng)度和偏振位置角，反映輻射區(qū)域以及輻射傳播途徑上的磁場結(jié)構(gòu)。為了研究星際介質(zhì)的磁場，需要重構(gòu)法拉第色散函數(shù)。目前已有多種重構(gòu)法拉第色散函數(shù)的方法。文[1]提出了法拉第旋率綜合方法(Faraday Rotation Measure Synthesis, RMS)，通過觀測偏振量和法拉第色散函數(shù)之間的傅里葉關(guān)系重構(gòu)法拉第色散函數(shù)。因?yàn)橛^測波段是有限的，所以獲得的法拉第色散函數(shù)是真實(shí)色散函數(shù)與窗口函數(shù)的卷積，對(duì)復(fù)雜的源需要解卷積。與法拉第旋率綜合方法類似的還有小波方法[2-3]。此外，還有一些假定輻射模型的方法重構(gòu)法拉第色散函數(shù)，如Q-U擬合法[4]。

本文主要討論壓縮感知方法(或稱為壓縮采樣，Compressive Sensing/Sampling，原理見[5-6])。該方法的優(yōu)點(diǎn)在于利用極少的采樣數(shù)據(jù)也有可能重構(gòu)原始的完整信號(hào)。壓縮感知理論以信號(hào)的稀疏性為前提，通過編碼測量和重構(gòu)算法實(shí)現(xiàn)信號(hào)的精確重構(gòu)。目前，已經(jīng)提出了基于壓縮感知法的法拉第色散函數(shù)重構(gòu)算法，如文[7]采用最小l1范數(shù)法建立了適用于法拉第薄源(Faraday Thin，λ2Δφ? 1，這里Δφ為源在法拉第深度域上的延展范圍，即F(φ)的寬度)、法拉第厚源(Faraday Thick,λ2Δφ? 1)和同時(shí)含有法拉第薄源和法拉第厚源的算法，模擬結(jié)果顯示，該方法可以用于重構(gòu)法拉第色散函數(shù)，并與法拉弟旋率綜合方法進(jìn)行比較，發(fā)現(xiàn)用壓縮感知方法重構(gòu)的法拉第色散函數(shù)在數(shù)值和誤差上要優(yōu)于法拉弟旋率綜合方法。文[8]采用匹配追蹤算法建立了對(duì)同時(shí)含有法拉第薄源和法拉第厚源進(jìn)行重構(gòu)的方法。但這些方法目前還處于理論模擬階段，其實(shí)用性還有待研究，比如文[7]的算法中設(shè)置的頻率(波長)范圍過大，未將噪聲添加到實(shí)驗(yàn)數(shù)據(jù)中，算法中各種參數(shù)設(shè)置對(duì)重構(gòu)結(jié)果的影響也沒有進(jìn)行實(shí)驗(yàn)和討論以及實(shí)驗(yàn)樣本過少等。本文基于文[7]的方法，從實(shí)際觀測角度進(jìn)行數(shù)值模擬，對(duì)以上問題進(jìn)行更深入研究，以此探討其實(shí)用性。這里假定法拉第色散函數(shù)是稀疏的(即針對(duì)源包含一個(gè)或多個(gè)法拉第薄成分)，用Python語言編寫了重構(gòu)法拉第色散函數(shù)的程序。采用的頻率范圍為1.1～3.1 GHz，對(duì)應(yīng)澳大利亞望遠(yuǎn)鏡致密陣(Australia Telescope Compact Array, ATCT)上L波段的觀測頻率范圍[9]。

1 重構(gòu)方法

已知測量矩陣A∈RM × N(M?N)和未知信號(hào)x0在該矩陣的線性投影y∈RM，即y=Ax0，要由測量值y重構(gòu)未知信號(hào)x0。在壓縮感知理論中，假定未知信號(hào)x0為稀疏的，信號(hào)x0可以由測量值y通過求解最小l1范數(shù)問題精確重構(gòu)：

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

(4)

若未知信號(hào)x0為非稀疏的，可以通過某種變換(如小波變換)進(jìn)行稀疏性表示，即x0=ωX，X為該信號(hào)在ω變換域下的稀疏性表示，則有

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

(5)

將(2)式寫成矩陣形式：

Yf=P，

(6)

其中，f為法拉第色散函數(shù)，是一個(gè)N維矢量；Y為P和f之間傅里葉變換的離散表示，是一個(gè)M×N矩陣；P為觀測偏振量，是一個(gè)M維矢量，這里M為觀測頻率通道數(shù)，每個(gè)通道上接收不同波段的偏振量。

由于法拉第色散函數(shù)f是復(fù)變函數(shù)，可以將其實(shí)部和虛部分別寫成f.real和f.imag，則(6)式可變?yōu)?/p>

(7)

(8)

(9)

實(shí)驗(yàn)中，將均值為0、標(biāo)準(zhǔn)差σ為1(強(qiáng)度單位任意，下同)的高斯噪聲加到Q和U上，并設(shè)置在法拉第深度φ空間的分辨率δφ為ΔφFWHM的1/4，則總的格點(diǎn)數(shù)N(φ空間中坐標(biāo)點(diǎn)φi的個(gè)數(shù)，i=1, 2, …,N)為

(10)

其中，int()表示取整。表1列出了一些主要的實(shí)驗(yàn)參數(shù)。

表1 實(shí)驗(yàn)參數(shù)Table1 Experimental parameters

2 結(jié)果與討論

2.1 一般情況

在模擬實(shí)驗(yàn)中，假定一個(gè)法拉第色散函數(shù)包含5個(gè)法拉第薄源,這些源分別是：F(-164.1)=3 + 6i，F(xiàn)(-82.0)=-4 + 3i，F(xiàn)(0.0)=2 + 3i，F(xiàn)(68.4)=-8 - 5i，F(xiàn)(136.7)=5 + 3i，對(duì)應(yīng)的信噪比分別為：6.7， 5.0， 3.6， 9.4， 5.8。對(duì)包含這5個(gè)法拉第薄源的法拉第色散函數(shù)進(jìn)行了反解，結(jié)果如圖1(a)，

(11)

其中，p為偏振強(qiáng)度。從圖1可以看出，沒有噪聲時(shí)該模型準(zhǔn)確地重構(gòu)了原始圖像，而加入噪聲時(shí)，重構(gòu)結(jié)果出現(xiàn)了偏差，|dwm|分別為5.80 × 10-5，11.78。

通過壓縮感知方法重構(gòu)的法拉第色散函數(shù)隨著法拉第薄源偏振強(qiáng)度的大小、相位以及兩個(gè)法拉第薄源的法拉第深度之間間隔Δφ的不同而不同，甚至不能準(zhǔn)確地重構(gòu)出原始圖像。在這里分兩種情況進(jìn)行測試說明：(1)假定上述5個(gè)法拉第薄源的法拉第深度之間的間隔Δφ不變，改變其大小和相位，如F(-164.1)=3 - 6i，F(xiàn)(-82.0)=-4 + 3i，F(xiàn)(0.0)=5 + 3i，F(xiàn)(68.4)=8 - 5i,F(136.7)=-5 + 4i，對(duì)應(yīng)的信噪比分別為：6.7，5.0，5.8，9.4，6.4，重構(gòu)結(jié)果如圖1(b)；(2)假定大小和相位不變，改變其法拉第深度之間的間隔Δφ，如：F(-150.4)=3 + 6i，F(xiàn)(-68.4)=-4 + 3i，F(xiàn)(0.0)=2 + 3i，F(xiàn)(109.4)=-8 - 5i，F(xiàn)(164.1)=5 + 3i，對(duì)應(yīng)的信噪比分別為6.7，5.0，3.6，9.4，5.8，重構(gòu)結(jié)果如圖1(c)。

當(dāng)假定間隔Δφ不變時(shí)，從圖1可以看出，前3個(gè)源并沒有很好地重構(gòu)出原始圖像，其相位和大小存在一定誤差，其中，|dwm|分別為6.31，17.03。而假定大小和相位不變時(shí)，后面4個(gè)源都與原始圖像出現(xiàn)了較大偏差，其中，|dwm|分別為8.18，17.11。同時(shí)也發(fā)現(xiàn)，在不加入噪聲時(shí)，將這些法拉第薄源的法拉第深度之間的間隔Δφ設(shè)置得足夠大時(shí)，即Δφ> 90 rad·m-2時(shí)，無論改變其大小、相位還是間隔Δφ，該模型均能準(zhǔn)確地重構(gòu)原始圖像，即使在加入噪聲時(shí)，產(chǎn)生的誤差也很微小。當(dāng)然，以上測試的這些源的法拉第深度都是設(shè)置在該方法的格點(diǎn)N上進(jìn)行的。

通常，在實(shí)際應(yīng)用中，事先并不知道法拉第色散函數(shù)的法拉第深度以及兩個(gè)源的法拉第深度之間的間隔Δφ，要重構(gòu)的法拉第色散函數(shù)各個(gè)成分的法拉第深度是隨機(jī)的，一般不在這種方法設(shè)置的格點(diǎn)N上，并且Δφ也是隨機(jī)的，所以也無法確定該模型是否能夠很好地重構(gòu)法拉第色散函數(shù)。隨機(jī)選取一個(gè)包含5個(gè)法拉第薄源的色散函數(shù)(這些法拉第薄源的法拉第深度不在格點(diǎn)N上)：F(-114)=-6 - 1i，F(xiàn)(-60)=7 - 4i，F(xiàn)(48)=3 - 5i，F(xiàn)(102)=-3.5 + 3i，F(xiàn)(210)=9 + 3i，對(duì)應(yīng)的信噪比分別為：6.1，8.1，5.8，4.6，9.5。其測試結(jié)果如圖1(d)， 其中，|dwm|分別為0.40，10.20。

每次模擬都是隨機(jī)產(chǎn)生的高斯噪聲，所以每次重構(gòu)的結(jié)果不同。在實(shí)驗(yàn)中，對(duì)于法拉第色散函數(shù)的重構(gòu)結(jié)果隨著源的個(gè)數(shù)、大小、相位和噪聲的不同出現(xiàn)不同程度的誤差，同時(shí)也與該方法中設(shè)置的參數(shù)有關(guān)，如觀測通道M，φ空間的分辨率δφ。

圖1 法拉第色散函數(shù)的重構(gòu)結(jié)果。每幅圖從上往下，上圖為原始數(shù)據(jù)，中圖為無噪聲重構(gòu)結(jié)果，下圖為有噪聲重構(gòu)結(jié)果，粗實(shí)線表示振幅，細(xì)實(shí)線表示實(shí)部，細(xì)虛線表示虛部

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 大樣本實(shí)驗(yàn)

為了更直觀地觀察該方法的重構(gòu)效果，從噪聲、兩個(gè)源的法拉第深度之間的間隔Δφ、φ空間的分辨率δφ(δφ的不同會(huì)導(dǎo)致總的格點(diǎn)數(shù)N不同)3方面對(duì)重構(gòu)結(jié)果的影響分別進(jìn)行了大樣本統(tǒng)計(jì)學(xué)實(shí)驗(yàn)。以重構(gòu)的源的數(shù)目和φ的加權(quán)平均的差值描述重構(gòu)結(jié)果的好壞[10]。

首先考察信噪比對(duì)重構(gòu)結(jié)果的影響，這里假定法拉第色散函數(shù)只包含一個(gè)法拉第薄源，固定其法拉第深度為φ=20 rad·m-2，隨機(jī)模擬了1 000個(gè)信噪比范圍在(3, 60)之間的源(這里噪聲依然設(shè)置為1，信噪比相當(dāng)于法拉第色散函數(shù)的偏振強(qiáng)度)，φ空間的分辨率δφ依然為ΔφFWHM的1/4，由于只含有一個(gè)法拉第薄源，直接以φ的差值|φrecorved-φoriginal|描述重構(gòu)效果，并排除信噪比小于2的信號(hào)。結(jié)果如圖2(a)。從圖中可以看出，對(duì)于含有一個(gè)源的情況，可以很準(zhǔn)確地重構(gòu)出源的數(shù)目，而φ的差值總是在兩個(gè)值上，即重構(gòu)的φrecorved出現(xiàn)在相鄰的兩個(gè)格點(diǎn)上。為了減小重構(gòu)誤差，對(duì)重構(gòu)結(jié)果在峰值附近進(jìn)行最小二乘擬合，結(jié)果如圖2(b)。經(jīng)過擬合后，對(duì)于大多數(shù)源來說，φ的差值縮小為原來的1/2～1/4。由于重構(gòu)的法拉第深度φ總在相鄰的兩個(gè)格點(diǎn)上，是否也可以通過縮小φ空間分辨率的方式減小φ的重構(gòu)誤差呢？

圖2 從左到右為 (a) 無擬合；(b) 擬合后。每幅圖從上往下，上為重構(gòu)數(shù)目，下為φ的差值

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φ

下面考查兩個(gè)源的法拉第深度之間的間隔Δφ和φ空間的分辨率δφ對(duì)重構(gòu)結(jié)果的影響。假定一個(gè)法拉第色散函數(shù)包含2個(gè)法拉第薄源，固定其中一個(gè)源的偏振強(qiáng)度p=13.235，偏振位置角χ=2.749，法拉第深度φ=10.256 rad·m-2。第2個(gè)源是隨機(jī)模擬了1 000個(gè)偏振強(qiáng)度在(3, 30)、偏振位置角在(0, 2π)、間隔Δφ在(5, 200)之間的源。同時(shí)，對(duì)不同分辨率δφ的情況(即δφ取ΔφFWHM的1/2，1/3，1/4，1/5，1/6，1/7，1/8，1/9)分別進(jìn)行以上實(shí)驗(yàn)，這里排除信噪比小于1的信號(hào)。結(jié)果如圖3。

為了比較不同分辨率情況下的重構(gòu)效果，對(duì)不同分辨率下準(zhǔn)確重構(gòu)出源的數(shù)目的精確度進(jìn)行了比較，如圖4。從圖4可以看出，當(dāng)δφ取ΔφFWHM的1/3和1/4時(shí)，精確度最高，達(dá)到83%左右。在以上對(duì)噪聲影響的實(shí)驗(yàn)時(shí)，曾提出是否可以通過縮小分辨率的方式減小重構(gòu)誤差的問題，而在這個(gè)實(shí)驗(yàn)中可以看出，分辨率δφ并不是越小越好。

圖3 不同分辨率下間隔Δφ對(duì)重構(gòu)結(jié)果的影響。第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。每幅圖從上往下，上為重構(gòu)數(shù)目，下為φ的加權(quán)平均的差值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 不同分辨率下的重構(gòu)精度比較。x軸表示δφ/ΔφFWHM的比值，y軸表示重構(gòu)數(shù)目的精確度

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時(shí)實(shí)驗(yàn)為例，隨機(jī)選取的重構(gòu)圖像
Fig.5 Examples of reconstructed images, whenδφtakes 1/4 of ΔφFWHM

圖6dwm隨著間隔Δφ的變化
Fig.6 This picture showsdwmvaries with Δφ

3 結(jié) 論

在多個(gè)頻率上測量彌漫射電輻射的偏振是研究銀河系磁場的一種有效的方法。目前的法拉第旋率綜合法是利用偏振強(qiáng)度與法拉第色散函數(shù)之間是傅里葉變換對(duì)，把觀測到的偏振強(qiáng)度做傅里葉逆變換得到法拉第色散函數(shù)[1]。因受觀測波段的限制，得出的解需要進(jìn)行解卷積，所以需要一種更強(qiáng)有力的方法解決這一問題。

本文基于壓縮感知的法拉第色散函數(shù)重構(gòu)方法應(yīng)用到實(shí)際觀測頻率范圍中，并對(duì)可能影響該方法重構(gòu)效果的因素分別進(jìn)行了實(shí)驗(yàn)。實(shí)驗(yàn)結(jié)果表明，該方法的重構(gòu)效果隨著法拉第薄源的個(gè)數(shù)、源偏振強(qiáng)度的大小、相位以及兩個(gè)法拉第薄源的法拉第深度之間間隔Δφ的不同而出現(xiàn)較大差異，具有較大的隨機(jī)性。大樣本統(tǒng)計(jì)學(xué)實(shí)驗(yàn)發(fā)現(xiàn)，當(dāng)φ空間的分辨率δφ取ΔφFWHM的1/3和1/4時(shí)，重構(gòu)結(jié)果精確度最高，對(duì)其重構(gòu)結(jié)果各峰值附近再次進(jìn)行最小二乘擬合后得到的法拉第深度φ值更接近原始值。通過壓縮感知法重構(gòu)后的再次擬合，可以得到更為精確的法拉第深度φ值，因此該方法可以作為現(xiàn)有方法的一種補(bǔ)充，同時(shí)該方法對(duì)于重構(gòu)得出的偏振強(qiáng)度普遍偏低，重構(gòu)效果具有較大的隨機(jī)性，所以，還需要不斷地改進(jìn)，提高其穩(wěn)定性和精確度，以便更好地應(yīng)用于實(shí)際觀測數(shù)據(jù)中。