張亞州，張海龍,3,4*，張 萌，王 杰，冶鑫晨，王萬瓊，李 嘉，杜 旭

(1. 中國科學(xué)院新疆天文臺，新疆 烏魯木齊 830011；2. 中國科學(xué)院大學(xué)，北京 100049；3. 中國科學(xué)院射電天文重點實驗室，江蘇 南京 210033；4. 國家天文科學(xué)數(shù)據(jù)中心，北京 100101)

隨著射電天文技術(shù)的進(jìn)步，射電望遠(yuǎn)鏡向大口徑及陣列方向發(fā)展。中國500 m口徑球面射電望遠(yuǎn)鏡(Five-hundred-meter Aperture Spherical radio Telescope, FAST)[1]、荷蘭低頻陣列(the Low-Frequency Array, LOFAR)[2]以及美國新墨西哥州的甚大天線陣(Very Large Array, VLA)[3]等觀測設(shè)備具備高靈敏度和高分辨能力。望遠(yuǎn)鏡在觀測宇宙微弱天文信號的同時，也接收到干擾信號，強(qiáng)度比微弱天文信號強(qiáng)幾個數(shù)量級的人為射電干擾信號，在射電天文數(shù)據(jù)處理過程中必須得到有效抑制。

射頻干擾來源多種多樣，文[4]將射頻干擾分為外部射頻干擾和內(nèi)部射頻干擾，外部射頻干擾特指天文臺址外產(chǎn)生的無線干擾，常見的有人造衛(wèi)星[5](北斗和全球定位系統(tǒng)導(dǎo)航衛(wèi)星等)、飛機(jī)、附近的基站信號以及廣播信號等；內(nèi)部射頻干擾是指在天文臺址內(nèi)使用的電子設(shè)備，如各類數(shù)字終端設(shè)備、計算機(jī)、微波爐、咖啡機(jī)、熒光燈以及其他電子儀器等產(chǎn)生的射頻干擾。

射電天文學(xué)發(fā)展過程中，射頻干擾抑制一直是研究熱點。文[6]提到天文信號通常為寬頻帶且在時間尺度上平滑變化，絕大部分射頻干擾與天文信號相比有明顯的差異。針對射頻干擾抑制，文[7]提出了CUSUM(Cumulative Sum)方法，該方法簡單高效，但依賴于單點采樣，對噪聲比較敏感。為了克服這一局限性，文[8]提出了基于計算相鄰采樣點的組合閾值表示算法SumThreshold，文[9]利用奇異值分解來解決巨型米波射電望遠(yuǎn)鏡(Giant Metrewave Radio Telescope, GMRT)寬帶射頻干擾問題。這些算法本質(zhì)上基于閾值實現(xiàn)，處理幅值遠(yuǎn)大于天文信號的射頻干擾信號時效果較好，但對于特征與天文信號接近的射頻干擾信號，閾值法可能誤刪天文信號，導(dǎo)致信息丟失，影響天文研究的準(zhǔn)確性和有效性。

瞬時射頻干擾是某些設(shè)備或系統(tǒng)正常運(yùn)行時釋放的射頻信號。繼電器、電動機(jī)和開關(guān)電源都會產(chǎn)生瞬時干擾。由于瞬時干擾信號是寬頻、脈沖的，而且持續(xù)時間非常短，難以識別來源。文[10]提出了基于字典的瞬時射頻干擾分類算法，使用隱馬爾科夫模型將瞬態(tài)射頻干擾識別為子序列，并能從觀測數(shù)據(jù)中提取瞬時射頻干擾序列。

1 射頻干擾信號特征提取算法

1.1 信號周期計算

大部分射頻干擾信號在時間尺度上呈現(xiàn)明顯的周期波動。為了分析射頻干擾信號的特征信息，我們需要對信號分窗處理，窗口大小與信號周期信息密切相關(guān)，計算信號周期是分析特征信息的重要前提。

信號周期計算算法有傳統(tǒng)的倒譜法、自相關(guān)函數(shù)法、平均幅度差函數(shù)法和小波變換法等[11]。自相關(guān)函數(shù)法是一種非常有效的低信噪比信號周期估計算法[12]?；谧韵嚓P(guān)函數(shù)法，本文設(shè)計了一種自相關(guān)卷積滑動算法計算射頻干擾信號的周期。算法過程為

(1)對一段長度為N的射頻干擾信號振幅取絕對值；

(2)從射頻干擾信號上截取一段長度為M的信號作為滑動窗口信號；

(3)滑動窗口與射頻干擾信號從起始位置進(jìn)行卷積，每次向后滑動d個采樣點；

(4)重復(fù)步驟(3)，直至滑動窗口移動到射頻干擾信號尾部；

(5)利用卷積后結(jié)果的峰值信息計算信號的周期。

圖1 (a)k7ing-3545 kHz噪聲信號時域圖；(b)k7ing-3545 kHz噪聲信號周期圖

1.2 射頻干擾信號特征提取

峰值是射頻干擾信號重要的特征信息，能夠刻畫信號波動的變化細(xì)節(jié)。射頻干擾信號特征提取過程的具體步驟為

(1)選取射頻干擾信號窗口；

(2)對該窗口信號的振幅取絕對值；

(3)獲取該信號窗口的峰值信息；

(4)對離散峰值信息進(jìn)行平滑預(yù)處理；

(5)處理后的峰值信息作為射頻干擾信號的特征信息。

算法第1步選取射頻干擾信號窗口，對一段包含射頻干擾的信號提取特征信息時，我們需要選取一段時間跨度盡可能短且至少包含一個完整射頻干擾周期信息的信號段，根據(jù)1.1節(jié)提出的射頻干擾信號周期計算方法確定信號周期。射頻干擾信號相對于其他信號振幅較大，起始和終止位置明確，算法將信號振幅突增的位置作為窗口的起始點。射頻干擾信號特征提取窗口大小與計算得到的周期相關(guān)，計算機(jī)顯示器產(chǎn)生的干擾信號窗口選擇如圖2，窗口大小為D。第2步對信號的振幅求絕對值，以獲取更多的峰值信息。第3步利用振幅特性，獲取窗口的峰值信息。第4步對獲取的峰值信號進(jìn)行預(yù)處理，將相鄰且振幅大小相差較小的值用均值代替，平滑峰值曲線，減少局部波動，縮小同類別信號差異，提高識別效率。圖3(a)為射頻干擾信號峰值信息折線圖，波動趨勢不平穩(wěn)，不利于后續(xù)特征識別。平滑處理后射頻干擾信號的峰值曲線如圖3(b)，曲線相對平穩(wěn)，能夠縮小同類信號的差異。平滑處理后的峰值信息作為識別射頻干擾信號的特征模板，可以進(jìn)行相似度計算。

圖2 射頻干擾信號窗口選取

圖3 (a)射頻干擾信號的峰值曲線；(b)平滑后射頻干擾信號的峰值曲線

2 射頻干擾信號特征相似度算法

對未知信號進(jìn)行特征提取后，計算與已知特征模板的相似程度，實現(xiàn)對未知射頻干擾信號的識別和分類。算法采用分段打分策略計算兩個序列的相似程度，每個序列根據(jù)相鄰離散點可以分為多個區(qū)段，比較兩個序列對應(yīng)區(qū)段的趨勢是否相同，進(jìn)行記分，趨勢相同加分，否則減分。

在相似度計算前，要確保這兩個序列處于對齊狀態(tài)，即確定兩個序列的區(qū)段是否一一對應(yīng)。待識別信號具體細(xì)節(jié)未知，提取峰值信息后，離散的峰值信息幾乎不可能與特征模板對齊，即相位信息沒有對齊，因此，峰值位置信息等都需要重新調(diào)整，為計算相似度做準(zhǔn)備。

文[13]提出了一種測量時間序列相似性的方法，即動態(tài)時間規(guī)整，該方法能夠比較不同長度的時間序列[14]，已廣泛應(yīng)用于語音信號處理。該算法計算得到的歐幾里得距離越小，表明兩種聲音模式的相似度越高[15]。利用動態(tài)時間規(guī)整算法計算兩個長度不同序列的歐幾里得距離最小時，兩個序列的相位點對應(yīng)情況如圖4。

圖4 使用動態(tài)時間規(guī)整算法后兩個序列的相位點對應(yīng)情況

基于動態(tài)時間規(guī)整算法，本文設(shè)計了射頻干擾信號特征相似度計算算法，具體過程為

(1)通過動態(tài)時間規(guī)整算法計算兩組序列在最短歐幾里得距離條件下離散點的對應(yīng)情況；

(2)計算未知信號離散點的權(quán)重值；

(3)分段比較，如果對應(yīng)段的變化趨勢相同，相似度分?jǐn)?shù)增加，否則相似度分?jǐn)?shù)減??；

(4)累計各段的分?jǐn)?shù)，總分?jǐn)?shù)代表兩個序列的相似程度。

射頻干擾信號特征相似度計算算法第1步利用動態(tài)時間規(guī)整算法求得待比較序列在最小歐幾里得距離下離散點的對應(yīng)情況，為后續(xù)計算做準(zhǔn)備。第2步計算被比較序列的離散點權(quán)值，計算方法為

(1)

其中，D為序列的離散點點數(shù)；S(t)為序列第t個離散點的幅值；W(t)為序列第t個離散點的權(quán)重，離散值越大，權(quán)重越大，在信號特征中的代表性越強(qiáng)。為了使相似度介于-1到1之間，第一個和最后一個的權(quán)重設(shè)置為0。

第3步進(jìn)行序列相似度計算，采用分段打分累計的方式。通過動態(tài)時間規(guī)整算法計算每一區(qū)段的對應(yīng)情況，采用分段比較方式，如圖4。兩個序列的區(qū)段有趨勢相同、趨勢不同和一對多3種情況，兩個區(qū)段的趨勢相同時加分，計算公式為

score=score+M[i] ，

(2)

其中，score是兩個序列相似度的累積分?jǐn)?shù)；M[i]是第i個區(qū)段的權(quán)重，大小為區(qū)段的兩個離散點權(quán)重的平均值，即0.5×(W[i]+W[i+1])。兩個區(qū)段的趨勢不同時減分，計算公式為

score=score-M[i].

(3)

一對多情況不進(jìn)行加減分。計算并累計各段的分?jǐn)?shù)后得到兩個序列相似度值。

3 實驗結(jié)果與分析

本文選擇ARRL官網(wǎng)提供的干擾源作為實驗數(shù)據(jù)，選擇7種射頻干擾源進(jìn)行交叉測試，基本信息如表1。為方便對射頻干擾信號進(jìn)行特征提取及后續(xù)識別測試，利用1.1節(jié)提出的射頻干擾信號周期計算法得到各實驗數(shù)據(jù)周期采樣點數(shù)，以便下一步特征提取及識別測試。

表1 射頻干擾源的基本信息

隨機(jī)選取每個信號的10個周期作為測試數(shù)據(jù)，提取相應(yīng)候選特征模板與其他信號進(jìn)行交叉測試，選擇與相同信號相似度大且與其他信號相似度較低的候選模板為最優(yōu)模板。如圖5(a)選取k7_2作為k7ing-3545 kHz的特征模板，(b)選取ks_9作為ks2am-streetlight的特征模板，(c)選取mo_4作為monitor的特征模板，(d)選取n6_10作為n6rce-sps-carrier的特征模板。圖6(a)選取pl_8作為plc-4的特征模板，(b)選取au_7作為ausoth的特征模板，(c)選取ot_4作為18120oth25的特征模板。

圖5 (a)選擇k7ing-3545 kHz不同模板的相似度結(jié)果；(b)選擇ks2am-streetlight不同模板的相似度結(jié)果；(c)選擇monitor不同模板的相似度結(jié)果；(d)選擇n6rce-sps-carrier不同模板的相似度結(jié)果

圖6 (a)選擇plc-4不同模板的相似度結(jié)果；(b)選擇ausoth不同模板的相似度結(jié)果；(c)選擇18120oth25不同模板的相似度結(jié)果

提取每個信號最優(yōu)特征模板所在組的數(shù)據(jù)，結(jié)果如表2。相同信號的相似度在所在行最大，表明提取的射頻干擾特征模板能夠較正確地識別射頻干擾來源。算法的相似度計算基于特征模板得出，權(quán)值計算依賴于特征模板，即計算兩個序列的相似度時，權(quán)值的大小只與特征模板有關(guān)，因此與對角線互為對稱的數(shù)值不同。同一信號在不同周期內(nèi)存在差別，造成同一信號的相似度小于1。

表2 射頻干擾信號相似度驗證平均值結(jié)果

表2中個別數(shù)據(jù)偏高，如n6rce-sps-carrier行plc-4列的數(shù)值為0.726。通過分析得知，待比較的一組序列振幅波動趨勢比較一致，導(dǎo)致相似度偏高，后續(xù)工作將繼續(xù)優(yōu)化算法，以提高精度。

4 總 結(jié)

本文設(shè)計了滑動卷積周期計算算法，利用卷積后的峰值區(qū)間完成了信號平均周期計算；基于設(shè)計的峰值提取算法實現(xiàn)了信號特征模板提取；基于動態(tài)時間規(guī)整算法和打分策略設(shè)計了信號特征識別算法，實現(xiàn)了未知信號的識別和分類。本文分別對從ARRL官網(wǎng)下載的射頻干擾數(shù)據(jù)進(jìn)行互相關(guān)與自相關(guān)計算，實驗結(jié)果表明，信號生成的模板與原始信號的相似度明顯高于其他信號，說明本文提出的算法可以有效提取特征信息并生成信號特征模板，且利用特征識別算法能對射頻干擾進(jìn)行正確分類。本文提出的方法可以對射頻干擾進(jìn)行細(xì)粒度識別和標(biāo)記，有望為射頻干擾特征識別和標(biāo)記提供一種新的解決方案。

本文中算法實現(xiàn)代碼以及實驗數(shù)據(jù)已在碼云倉庫(https://gitee.com/zyazhou/rfi-feature-recognition.git)開源。