洪迎盈，吳亞軍，卜朝暉，楊 珣，劉 爻，王 琪

（1.上海理工大學(xué)生物醫(yī)學(xué)工程研究所，上海 200093；2.中國(guó)科學(xué)院上海天文臺(tái)射電天文技術(shù)實(shí)驗(yàn)室，上海 200030）

0 引言

射電天文通過(guò)接收來(lái)自天體的微弱輻射，利用電磁波頻譜研究天體［1］。射電望遠(yuǎn)鏡是一種寬帶無(wú)線電接收機(jī)，用于探測(cè)和分析來(lái)自天體的微弱無(wú)線電，靈敏度遠(yuǎn)高于地面通信接收機(jī)。由于人類電子設(shè)備的使用頻率劇增，電磁環(huán)境變得十分復(fù)雜［2］，會(huì)造成射電天文收到的信號(hào)為天文信號(hào)、系統(tǒng)總噪聲和RFI 信號(hào)的組合［3］，而RFI 信號(hào)通常會(huì)影響天文信號(hào)檢測(cè)精度，甚至導(dǎo)致信號(hào)失效。因此，研究消除射頻干擾的方法具有是十分重要的現(xiàn)實(shí)意義。

消除射頻干擾信息的原則是在保持原始信息完整度的前提下，盡可能去除RFI信號(hào)［4］?，F(xiàn)階段，國(guó)內(nèi)外學(xué)者已經(jīng)提出多種方法消除RFI信號(hào)，包括手動(dòng)標(biāo)記法［5］、閾值處理法［6］、基于機(jī)器學(xué)習(xí)［7-8］的方法等。其中，手動(dòng)標(biāo)記法的缺點(diǎn)在于面對(duì)海量數(shù)據(jù)時(shí)，標(biāo)記效率極低；閾值處理法通過(guò)計(jì)算信號(hào)數(shù)據(jù)點(diǎn)與邊界點(diǎn)的最小差值是否高于閾值以標(biāo)記RFI 信號(hào)，然而閾值的選擇依賴于RFI 信號(hào)和天文信號(hào)的特征差異，容易將微弱的天文信號(hào)誤判為RFI 信號(hào)；基于機(jī)器學(xué)習(xí)的方法依賴于模型建立的準(zhǔn)確性，當(dāng)模型存在誤差時(shí)將影響RFI 信號(hào)的標(biāo)記，并且建模過(guò)程復(fù)雜，訓(xùn)練時(shí)間較長(zhǎng)。

目前研究表明，利用奇異值分解（Singular Value Decomposition，SVD）可有效消除RFI 信號(hào)［9-11］。該方法利用觀測(cè)信號(hào)構(gòu)建相應(yīng)的觀測(cè)信號(hào)矩陣，然后對(duì)矩陣進(jìn)行奇異值分解，令干擾分量對(duì)應(yīng)的奇異值為零，從而消除干擾信號(hào)。然而，在實(shí)際測(cè)量過(guò)程中，RFI信號(hào)與天文信號(hào)并非始終呈正交狀態(tài)，因此在利用奇異值分解時(shí)，會(huì)在一定程度上消除部分天文信號(hào)。

針對(duì)上述問(wèn)題，本文提出奇異值優(yōu)化分割方法對(duì)射頻干擾進(jìn)行消除。首先，對(duì)觀測(cè)信號(hào)的每個(gè)周期進(jìn)行子帶分解，構(gòu)造相應(yīng)的觀測(cè)信號(hào)矩陣；然后，對(duì)觀測(cè)信號(hào)矩陣進(jìn)行奇異值分解，尋找干擾信號(hào)對(duì)應(yīng)的奇異值，按輸出信噪比最大的原則對(duì)奇異值進(jìn)行優(yōu)化分割，得到最優(yōu)分割比和優(yōu)化后的奇異值分量；最后，利用最優(yōu)分割比對(duì)相應(yīng)的奇異值進(jìn)行優(yōu)化以去除RFI 信號(hào)。該方法解決了傳統(tǒng)方法在消除RFI 信號(hào)時(shí)，對(duì)射電天文信號(hào)造成的影響，并在消色散前進(jìn)行射電干擾消除，有效避免射電干擾信號(hào)隨消色散處理發(fā)生更嚴(yán)重的變化。

1 奇異值優(yōu)化分割算法

脈沖星是快速旋轉(zhuǎn)的中子星，在地球上能探測(cè)到其發(fā)射的周期性高頻率無(wú)線電脈沖［12-13］，本文針對(duì)脈沖星觀測(cè)數(shù)據(jù)提出基于奇異值優(yōu)化分割的射頻干擾消除方法，具體原理流程如圖1所示。

1.1 信號(hào)子帶分解

由于射電觀測(cè)的天文信號(hào)通常呈現(xiàn)寬帶、平滑的狀態(tài)，而射頻干擾信號(hào)經(jīng)常在時(shí)頻平面上呈現(xiàn)為高強(qiáng)度像素，因此選擇在時(shí)頻域上去除射頻信號(hào)［14-15］。本文利用多相復(fù)數(shù)調(diào)制濾波器組，按脈沖星信號(hào)的每個(gè)周期，將觀測(cè)信號(hào)分解成子帶信號(hào)，構(gòu)建觀測(cè)信號(hào)矩陣。

可將一維射電觀測(cè)信號(hào)x(n)通過(guò)一個(gè)濾波器組分解成M個(gè)子帶信號(hào)x0(n)，x1(n)，…，xM-1(n)，構(gòu)造M×N的信號(hào)矩陣A，如式（1）所示：

Fig.1 Flow of radio frequency interference cancellation method based on singular value optimization segmentation圖1 基于奇異值優(yōu)化分割的射頻干擾消除流程

式中，M為濾波器組通道數(shù)，N為觀測(cè)信號(hào)長(zhǎng)度，n=0，1，…，N-1。

采用多相復(fù)數(shù)調(diào)制濾波器組將觀測(cè)信號(hào)分解為221個(gè)子帶信號(hào)，每個(gè)子帶的帶寬為0.5MHz，整個(gè)濾波器組的最高起始中心頻率為2 300MHz，射電觀測(cè)信號(hào)經(jīng)子帶分解后得到觀測(cè)信號(hào)矩陣。

1.2 射頻干擾消除

通過(guò)濾波器組構(gòu)建觀測(cè)信號(hào)矩陣，基于奇異值優(yōu)化分割方法消除射電觀測(cè)信號(hào)中存在的RFI 信號(hào)。其中，奇異值分解是一個(gè)能適用于任意矩陣的分解方法［16］，M×N觀測(cè)信號(hào)矩陣A經(jīng)過(guò)SVD 分解后可表示為：

式中，Σr為主對(duì)角矩陣，主對(duì)角線的值為奇異值，并按照大小進(jìn)行排列，r為矩陣A的秩。從大到小的奇異值反映了天文信號(hào)、系統(tǒng)總噪聲和RFI 信號(hào)在觀測(cè)信號(hào)矩陣中的能量分布，若奇異值越大，則表示對(duì)應(yīng)成分的能量、占觀測(cè)信號(hào)的占比越大。首先，利用奇異值分布可得到RFI 信號(hào)分量對(duì)應(yīng)的奇異值，然后將其置零，即可得到去除干擾后的觀測(cè)信號(hào)。

由矩陣U、V中的行和列構(gòu)成標(biāo)準(zhǔn)正交基可得?i ∈[1，…，M]，存在，因此U、V中每行、每列至少存在一個(gè)值非零，若將奇異值σi置零，則會(huì)改變A中的大小，信號(hào)矩陣中每個(gè)值的大小會(huì)發(fā)生變化［5］。在實(shí)際射電觀測(cè)中，由于大部分RFI 信號(hào)和天文信號(hào)非正交，若觀測(cè)信號(hào)中RFI 信號(hào)對(duì)應(yīng)奇異值為零，在消除RFI 信號(hào)的時(shí)會(huì)損害天文信號(hào)，導(dǎo)致觀測(cè)信號(hào)信噪比降低。

為此，需要先對(duì)RFI 信號(hào)對(duì)應(yīng)的奇異值進(jìn)行優(yōu)化分割，改變?nèi)≈荡笮。x取輸出信噪比最大時(shí)的最優(yōu)值，算法具體實(shí)現(xiàn)流程如圖2所示。

Fig.2 Flow of SVD optimal segmentation algorithm圖2 奇異值優(yōu)化分割算法流程

首先，利用奇異值差分譜理論［17-18］尋找RFI 信號(hào)對(duì)應(yīng)的奇異值σRFI，并進(jìn)行優(yōu)化分割處理；然后，將σRFI劃分為N等分，按分割比k/N的原則得到σRFI的N+1 個(gè)分割值σk=()σRFI，k=0，1，…，N；最后，用σk替代σRFI恢復(fù)出射電觀測(cè)信號(hào)，取得輸出信噪比最大時(shí)的最優(yōu)分割比()opt。

由于大多數(shù)RFI 信號(hào)的強(qiáng)度相較于天文信號(hào)和噪聲更高，在通常情況下較大的奇異值主要反映包含能量較大的RFI 信號(hào)，小的奇異值則主要反映天文信號(hào)和背景噪聲。因此根據(jù)奇異值差分譜理論［17-18］，通過(guò)奇異值分布得到RFI 信號(hào)對(duì)應(yīng)的奇異值。圖3 為天馬射電望遠(yuǎn)鏡觀測(cè)到的PSR B0329+54 脈沖星數(shù)據(jù)進(jìn)行SVD 分解得到奇異值差分譜。

由圖3 可見(jiàn)，峰值出現(xiàn)在第1 個(gè)奇異值上，與其它奇異值形成了較大差異。由上述理論可判斷最大的奇異值即為RFI 信號(hào)對(duì)應(yīng)的奇異值，然后按照輸出信噪比最大的原則對(duì)其進(jìn)行優(yōu)化分割，并保留剩余的奇異值，恢復(fù)處理后的射電觀測(cè)信號(hào)，即為去除RFI信號(hào)后的射電觀測(cè)信號(hào)。

Fig.3 Singular value difference spectrum of actual radio observation signal圖3 實(shí)際射電觀測(cè)信號(hào)奇異值差分譜

1.3 評(píng)價(jià)指標(biāo)

為客觀、定量的評(píng)估本文方法對(duì)干擾抑制的有效性，將射頻干擾抑制程度G_proc 和脈沖星信號(hào)損耗參數(shù)R_proc 作為射頻干擾抑制效果的評(píng)估指標(biāo)。其中，Gproc用來(lái)定量評(píng)估RFI 信號(hào)的抑制程度；Rproc用來(lái)檢測(cè)對(duì)RFI 信號(hào)進(jìn)行抑制后，產(chǎn)生的有效信號(hào)損失。具體計(jì)算公式如下：

式中，SNR為信噪比。

2 仿真實(shí)驗(yàn)

2.1 射電觀測(cè)信號(hào)模型

為驗(yàn)證方法的有效性，根據(jù)實(shí)際脈沖星觀測(cè)數(shù)據(jù)特點(diǎn)建立射電觀測(cè)信號(hào)模型。射電觀測(cè)信號(hào)模型一般表示為：

式中，xsig(n)為脈沖星信號(hào)，xRFI(n)為RFI 信號(hào)，xsig(t)為窄脈沖信號(hào)。系統(tǒng)總噪聲xsys(n)為天空背景噪聲和接收機(jī)噪聲的組合，兩者都為零均值不相關(guān)的高斯分布隨機(jī)信號(hào)。

在射電觀測(cè)中，RFI信號(hào)可能由電視廣播、調(diào)頻無(wú)線電傳輸和手機(jī)導(dǎo)航通訊等電子設(shè)備產(chǎn)生［19］，在二維時(shí)頻圖中的特點(diǎn)為占用一個(gè)或多個(gè)頻率通道，幾乎持續(xù)整個(gè)時(shí)間范圍［20］，具體數(shù)學(xué)表達(dá)式為：

式中，a(n)為高斯包絡(luò)，ω0為載波頻率，φ為其在0～2π之間均勻分布的隨機(jī)相位。

采用合成的射電觀測(cè)信號(hào)分析消除RFI信號(hào)后的輸出信噪比，以確定分割比N。根據(jù)式（7）中的射電觀測(cè)信號(hào)模型，基于相應(yīng)的模型參數(shù)模擬射電觀測(cè)信號(hào)中的脈沖星信號(hào)、系統(tǒng)總噪聲和RFI 信號(hào)，利用表1的射電觀測(cè)信號(hào)仿真參數(shù)，建立觀測(cè)信號(hào)模型并分析消除RFI信號(hào)的效果。

Table 1 Simulation parameter setting表1 仿真參數(shù)設(shè)置

2.2 射頻干擾消除

利用復(fù)數(shù)調(diào)制濾波器組對(duì)仿真信號(hào)進(jìn)行子帶分解，得到221 個(gè)子帶信號(hào)，以此構(gòu)建信號(hào)矩陣，對(duì)RFI 信號(hào)進(jìn)行處理。圖4 為SVD 分解仿真信號(hào)得到的奇異值差分譜，峰值出現(xiàn)在第1個(gè)奇異值上，即為RFI信號(hào)對(duì)應(yīng)的奇異值。

Fig.4 Simulation signal singular value difference spectrum圖4 仿真信號(hào)奇異值差分譜

由實(shí)驗(yàn)可知，射電觀測(cè)信號(hào)中RFI信號(hào)的抑制程度Gproc與最大奇異值的分割比存在一定的關(guān)系。因此，為了確定最優(yōu)分割比，本文選擇先確定N值的大小，再通過(guò)改變k值，使分割比發(fā)生變化，從而得到不同的輸出信噪比進(jìn)行比較。圖5展示了不同N值對(duì)仿真信號(hào)輸出信噪比的影響。

Fig.5 Influence of k value on output signal-to-noise ratio of simulation signal under different N value圖5 不同N值時(shí)k值對(duì)仿真信號(hào)輸出信噪比的影響

由此可知，不同N值的輸出信噪比差異的相對(duì)比值為：

式中，(SNR)200為N取200 時(shí)，仿真信號(hào)的最大輸出信噪比；(SNR)N為N分別取10、20、50、100、200 時(shí)，仿真信號(hào)的最大輸出信噪比。

由式（9）可分別算得(SNR)N與(SNR)200輸出信噪比差異的相對(duì)比值。如圖6 所示，當(dāng)N＜20 時(shí)，輸出信噪比差異的相對(duì)比值明顯大于N≥20 時(shí)的相對(duì)比值，為了減少計(jì)算量，增加算法的通用性，本文選擇將N=20，然后根據(jù)輸出信噪比隨k值變化曲線（見(jiàn)圖5（b）），可見(jiàn)當(dāng)k=5時(shí)輸出信噪比最大。因此，本文最終設(shè)定最優(yōu)分割比為=520=0.25。

表2 為利用射頻干擾抑制程度Gproc和脈沖星信號(hào)損耗參數(shù)Rproc評(píng)價(jià)本文算法的結(jié)果，利用最優(yōu)分割比消除RFI信號(hào)的結(jié)果如圖7 所示。可見(jiàn)奇異值優(yōu)化分割方法對(duì)RFI信號(hào)的抑制效果明顯，奇異值優(yōu)化分割方法處理后的信號(hào)的輸出信噪比增益明顯，對(duì)天文信號(hào)的影響較小。

Fig.6 Influence of N on the relative ratio of simulated signal SNR difference圖6 N值對(duì)仿真信號(hào)信噪比差異的相對(duì)比值的影響

Table 2 Evaluation parameters of RFI mitigation algorithm表2 干擾消除算法評(píng)價(jià)指標(biāo)

Fig.7 Simulation signal RFI elimination effect圖7 仿真信號(hào)射頻干擾消除效果

3 實(shí)例驗(yàn)證與分析

3.1 實(shí)驗(yàn)數(shù)據(jù)

實(shí)驗(yàn)數(shù)據(jù)來(lái)自上海天馬射電望遠(yuǎn)鏡對(duì)PSR B0329+54脈沖星在2020 年7 月5 日的觀測(cè)結(jié)果［21-23］。觀測(cè)中心頻率2 300MHz，觀測(cè)波段為2 180-2 300MHz，脈沖星周期為0.715s。

3.2 實(shí)驗(yàn)結(jié)果

為驗(yàn)證奇異值優(yōu)化分割方法的真實(shí)性和可行性，利用該方法對(duì)實(shí)際觀測(cè)數(shù)據(jù)進(jìn)行處理，再經(jīng)過(guò)子帶分解構(gòu)建信號(hào)矩陣。如圖8（a）所示，RFI 信號(hào)分別出現(xiàn)在80-83 通道，124-126 通道，181 通道附近，信號(hào)具有強(qiáng)度大、持續(xù)時(shí)間長(zhǎng)的特點(diǎn)。信號(hào)矩陣經(jīng)過(guò)奇異值優(yōu)化分割去干擾后的結(jié)果如圖8（b）所示，射頻干擾已基本被消除，并且射頻干擾處的天文信號(hào)幾乎不受影響。

Fig.8 RFI mitigation for PSR B0329+54圖8 PSR B0329+54脈沖星觀測(cè)信號(hào)射頻干擾消除效果

圖9 為利用奇異值優(yōu)化分割方法消除實(shí)際觀測(cè)信號(hào)中的RFI 信號(hào)的時(shí)域結(jié)果。其中，圖9（a）、（b）均能看到PSR B0329+54 脈沖星的尖峰結(jié)構(gòu)，但圖9（b）中尖峰峰值相較于圖9（a）明顯更高，說(shuō)明RFI信號(hào)已經(jīng)得到抑制。

Fig.9 RFI mitigation for PSR B0329+54圖9 PSR B0329+54脈沖星觀測(cè)信號(hào)射頻干擾消除效果

接下來(lái)，對(duì)射電觀測(cè)信號(hào)進(jìn)行奇異值優(yōu)化分割，當(dāng)N=20 時(shí)，隨著k值的變化，射電觀測(cè)信號(hào)輸出信噪比的變化情況如圖10所示。

Fig.10 Influence of k on the observed signal output SNR in RFI mitigation圖10 k值對(duì)觀測(cè)信號(hào)RFI消除后輸出信噪比的影響

由圖10 可見(jiàn)，當(dāng)k值增大時(shí)，輸出信噪比呈現(xiàn)先增大后減小的趨勢(shì)。當(dāng)k=5 時(shí)，輸出信噪比達(dá)到峰值，射頻干擾抑制效果最好，此時(shí)分割比為5/20。圖11 展示了不同N值對(duì)實(shí)際射電觀測(cè)信號(hào)輸出信噪比的影響。

由圖12 可見(jiàn)，根據(jù)式（9）可計(jì)算實(shí)際射電觀測(cè)信號(hào)中(SNR)N與(SNR)200時(shí)輸出信噪比差異的相對(duì)比值，當(dāng)N＜20 時(shí)輸出信噪比差異的相對(duì)比值明顯大于當(dāng)N≥20時(shí)的相對(duì)比值，因此本文將N值取為20，由輸出信噪比隨k值變化曲線（見(jiàn)圖11（b）），可見(jiàn)當(dāng)k=5 時(shí)輸出信噪比最大。由此可得，最優(yōu)分割比為=520=0.25，該結(jié)果與對(duì)仿真信號(hào)求取的最優(yōu)分割比一致。

圖13 為同一射電源在不同時(shí)間段的多組射電觀測(cè)信號(hào)，利用奇異值優(yōu)化分割在消除RFI 信號(hào)時(shí)發(fā)現(xiàn)，分割比=520 時(shí)信號(hào)輸出信噪比出現(xiàn)極大值。因此，設(shè)定奇異值分割比為=5/20具有通用性。

Fig.11 Influence of k value on output signal-to-noise ratio of observation signal under different N value圖11 不同N值時(shí)k值對(duì)觀察信號(hào)輸出信噪比的影響

Fig.12 Influence of N on the relative ratio of observed signal SNR difference圖12 N值對(duì)觀測(cè)信號(hào)信噪比差異的相對(duì)比值的影響

Fig.13 Influence of k on the output SNR in observed signals at different time圖13 不同時(shí)間段的射電觀測(cè)信號(hào)k值對(duì)輸出信噪比的影響

為評(píng)定本方法對(duì)實(shí)測(cè)信號(hào)射頻干擾的抑制效果，計(jì)算不同時(shí)間段射電觀測(cè)信號(hào)的處理后輸出信噪比增益大小Gproc，具體數(shù)據(jù)見(jiàn)表3。由此可見(jiàn)，奇異值優(yōu)化分割方法對(duì)不同天文觀測(cè)信號(hào)處理后的輸出信噪比增益顯著。

Table 3 Gproc after RFI mitigation in observed signals at different time表3 不同時(shí)間段的射電觀測(cè)信號(hào)射頻干擾消除后增益Gproc值

4 結(jié)語(yǔ)

本文通過(guò)奇異值優(yōu)化分割方法，消除射電觀測(cè)信號(hào)的射頻干擾。相較于傳統(tǒng)方法，該方法具有以下優(yōu)勢(shì)：①避免閾值法造成的誤判問(wèn)題；②避免機(jī)器學(xué)習(xí)方法中復(fù)雜的建模過(guò)程和模型誤差造成的影響；③減小奇異值分解在消除RFI信號(hào)的同時(shí)對(duì)天文信號(hào)造成的影響。

奇異值優(yōu)化分割方法按輸出信噪比最大的原則優(yōu)化RFI 對(duì)應(yīng)的奇異值，得到最優(yōu)分割比和優(yōu)化后的奇異值分量。該方法可用于消色散處理和信號(hào)積分之前，對(duì)于同一射電源的射電觀測(cè)信號(hào)，獲取的最優(yōu)分割比具有通用性，并且去除干擾后的輸出信號(hào)信噪比增益明顯，可達(dá)到1.2-1.4。后續(xù)將利用該方法針對(duì)不同射電源的射電觀測(cè)信號(hào)進(jìn)行去干擾處理，找到更具針對(duì)性和通用性的最優(yōu)分割比，以進(jìn)一步提升使用效率。