基于Simple Thresholding和CUSUM聯(lián)合算法的L波段太陽射電流量可觀測頻段分析

楊 航，董 亮，何樂生*

(1. 云南大學，云南 昆明 650500；2. 中國科學院云南天文臺，云南 昆明 650216)

太陽活動是空間天氣的源頭，日冕物質(zhì)拋射(Coronal Mass Ejection, CME)[1]事件是影響日地空間環(huán)境最主要的因素，射電觀測反映不同日冕高度的等離子體信息[2]，包括密度、磁場等。等離子體拋射時，其密度隨著遠離太陽而降低，隨著日冕物質(zhì)拋射遠離太陽，其發(fā)射的無線電頻率也隨之降低。為此，太陽物理界一致認為地基太陽射電望遠鏡能觀測到的無線電爆發(fā)頻率越低，日冕物質(zhì)拋射離地球越近。在不同頻率上的太陽射電觀測可以反演日冕高度的動態(tài)信息[3-5]。

同時，長期太陽射電流量監(jiān)測發(fā)現(xiàn)，不同頻段的流量變化與太陽活動周高度相關，是重要的空間天氣觀測指標。但另一方面，射電的觀測頻段與民用無線電使用頻段重合，極易干擾觀測數(shù)據(jù)。在射電望遠鏡建設之前的無線電環(huán)境測試與評估尤為重要，與脈沖星、活動星系核(Active Galactic Nucleus, AGN)光變等需要寬帶觀測不同[6-7]。由于高時間分辨率要求，太陽射電精密流量觀測時對一些偶發(fā)性的干擾并不敏感，這個短時 “壞點” 可以通過后期與太陽光學等數(shù)據(jù)進行比較，排除爆發(fā)可能。

L波段太陽射電爆發(fā)是導航通信潛在的干擾因素，與其他空間天氣事件不同，太陽爆發(fā)的強信號以光速傳播，是作用最快的空間天氣事件，從理論上來說，整個迎日面的導航接收機均有可能受到干擾[8-10]。文[11-13]提出了采用L波段多通道精密太陽射電望遠鏡同步監(jiān)測的辦法，實時監(jiān)測L波段太陽射電流量，在爆發(fā)早期發(fā)出預警信號。

但是L波段存在大量無線電干擾，在建設望遠鏡前，我們需要評估MHz級的觀測帶內(nèi)無線電干擾信號的分布規(guī)律，判斷可觀測太陽時間，對于長時間被干擾占用的信道進行排查，以防發(fā)生誤報等情況。針對即將在云南天文臺建設的L波段精密流量太陽射電望遠鏡(針對中高日冕進行研究和太陽射電爆發(fā)干擾導航通信預警，并作為廊坊太陽射電望遠鏡的備份)，本文介紹了對臺址的無線電測試結果和觀測頻段遴選方法以及遴選介于北斗B1, B2和B3頻點，全球定位系統(tǒng)L1和L2頻點之間7個5 MHz無線電干擾較少的無線電通帶。

1 系統(tǒng)組成

1.1 系統(tǒng)組成

測試系統(tǒng)由HL050天線、低噪聲放大器、穩(wěn)壓電源、FSU26頻譜儀和計算機組成，測試系統(tǒng)框圖如圖1。

系統(tǒng)框圖中，由HL050天線將采集的信號通過低噪聲放大器，傳輸?shù)紽SU26頻譜儀轉(zhuǎn)化成頻譜，并把數(shù)據(jù)傳輸給計算機進行顯示和存儲。其中穩(wěn)壓電源為放大器提供穩(wěn)定電壓，保證數(shù)據(jù)的準確性，計算機通過LabView圖形化編程，實現(xiàn)控制頻段范圍和極化方式，主要分為兩部分：HL050天線頻段控制和頻譜儀控制程序[14]。

圖1 測試系統(tǒng)Fig.1 Test system

1.2 測試設備及參數(shù)

HL050天線應用廣泛，主要參數(shù)見表1，其優(yōu)點是結構簡單、頻帶寬、功率容量大、調(diào)整和使用方便。

頻譜儀采用德國羅德施瓦茨公司的FSU26，工作頻率20 Hz～26.5 GHz，平均噪聲電平為-158 dBm，參數(shù)見表1。

表1 測試設備及參數(shù)Table 1 Test equipment and parameters

計算機利用LabView軟件圖形化編程實現(xiàn)，在LabView軟件圖形化編程中，可以調(diào)節(jié)測試的頻段起始頻率、極化方式、分辨率帶寬、積分時間和執(zhí)行次數(shù)，并將采集的信號保存在預設地址中，便于后續(xù)數(shù)據(jù)的處理分析。

1.3 監(jiān)測過程

測試地點為云南省昆明市中國科學院云南天文臺(25°1′48.538″N，102°48′14.180″E)，我們在東、西、南、北4個方向進行L波段的測試，系統(tǒng)頻率范圍為1～2 GHz，頻譜儀掃頻寬度設置為1 000 MHz。頻譜儀顯示帶寬和分辨率帶寬設置為300 kHz，掃描時間30 s。4個方向都測試了4次且時間超過8 h，共測試16次，數(shù)據(jù)均超過800組。具體測試日期及其組數(shù)如表2。

2 數(shù)據(jù)處理

本節(jié)主要介紹采集數(shù)據(jù)的格式、數(shù)據(jù)閾值遴選和本次最終信道的遴選過程。

表2 測試日期及組數(shù)(2022年)Table 2 Date of test and number of groups (in 2022)

2.1 數(shù)據(jù)介紹

本文采集的數(shù)據(jù)記錄格式為TXT，單位為dBi，縮小數(shù)據(jù)的絕對數(shù)值，方便記錄。由于整個測試過程設置頻譜儀輸出兩條線，所以每一組檢測數(shù)據(jù)分為兩列，如圖2。第1列為所采集數(shù)據(jù)的最大保持值，第2列為掃描時間內(nèi)的平均值，為每次檢測時判斷測試系統(tǒng)的穩(wěn)定效果。數(shù)據(jù)儲存為非線性的數(shù)據(jù)格式，這是一種對數(shù)關系，底數(shù)為10，數(shù)學關系為log10(P2/P1)，其中P2/P1表示功率比。后續(xù)數(shù)據(jù)處理需要對數(shù)據(jù)進行線性轉(zhuǎn)化[15-16]。

圖2 測試系統(tǒng)實測圖Fig.2 Actual measurement diagram of the test system

2.2 閾值遴選方案

針對即將在云南天文臺建設的L波段精密流量太陽射電望遠鏡，遴選介于北斗 B1(1 559.052～1 591.788 MHz)，B2(1 166.220～1 217.370 MHz)和B3(1 250.618～1 286.432 MHz)頻點，全球定位系統(tǒng)L1(1 574.397～1 576.443 MHz)和L2(1 226.577～1 228.623 MHz)頻點之間7個5 MHz干擾較少的無線電通帶。已知全球定位系統(tǒng)中L1頻點在北斗B1頻點范圍內(nèi)，且B2與L2頻點之間相差不超過10 MHz，無法遴選兩處相對潔凈的5 MHz通帶，B2右邊遴選的通帶與L2左邊遴選的通帶重合，一共7個通道(經(jīng)過無線電環(huán)境監(jiān)測帶內(nèi)無干擾信號)，全球定位系統(tǒng)和北斗系統(tǒng)的所有通信頻點都穿插在監(jiān)測頻點之間[17-20]。

本次測試目的是從L波段中遴選潔凈率高于95%的可用信道。遴選方式是通過繪制頻譜瀑布圖，直觀地選出較為干凈的信道，通過遍歷觀測數(shù)據(jù)超過閾值的占比確定頻譜潔凈率，當遴選的信道的潔凈率都大于95%，說明該信道可用。處理流程如下：

(1)選取一個短時 “寧靜” 的觀測頻段作為閾值計算數(shù)據(jù)樣本；

(2)在上述數(shù)據(jù)樣本中通過本文提出的閾值算法計算得到閾值；

(3)以此閾值遍歷其余觀測時間段數(shù)據(jù)，對比閾值，超過閾值時間和低于閾值時間進行對比，求得該時間段的干擾信號占比；

(4)利用算法計算得到的閾值對本文遴選的7個5 MHz干擾較少的無線電通帶進行評估。

2.2.1 閾值計算方案

閾值分析法因?qū)崿F(xiàn)簡單、檢測結果精度較高而廣泛應用。代表性閾值分析法有CUSUM法[21]和Simple Thresholding法[22]。在檢測某行(列)觀測數(shù)據(jù)時，CUSUM算法通過估算累積樣本的方差與平均值得到閾值，Simple Thresholding算法使用該行(列)的中位數(shù)作為閾值[23]。

本文首先監(jiān)測取對數(shù)之后的數(shù)據(jù)，取對數(shù)不會改變數(shù)據(jù)的性質(zhì)和相關關系，且壓縮了變量的尺度，數(shù)據(jù)更加平穩(wěn)，也消弱了模型的共線性、異方差性等，所以計算閾值時需要將數(shù)據(jù)線性化，公式為

(1)

樣本方差是表示射頻干擾的統(tǒng)計異常的重要估計值。由于絕對中位差(Median Absolute Deviation, MAD)和中位數(shù)是剔除異常值的最有效方法，因此，本文提出一種基于Simple Thresholding算法與CUSUM算法的改進閾值算法對臺址的L波段無線電測試結果和觀測頻段進行遴選，如果所選取信道閾值水平超過Th，視為射頻干擾。Th計算公式為

(2)

MAD=1.4286×med1≤i≤n{|xi-med(xi)|}，

(3)

(3)式中，n=2m時，med=0.5(xm+xm+1)；n=2m+1時，med=x(m+1)；MAD是Simple Thresholding算法使用的方差估計值，中位數(shù)是原始總體數(shù)組xi的中位數(shù)。然后用此中位數(shù)減去原始樣本的每個元素，創(chuàng)建與原始元素大小相同的新修改樣本。然后計算這個新樣本的中位數(shù)，并用恒定的比例因子1.428 6相乘[24]，使該估計與預期的高斯分布一致。得到線性計算的閾值Th之后，需要將閾值取對數(shù)轉(zhuǎn)換在監(jiān)測的數(shù)據(jù)中計算通帶的潔凈率，通過

Th*=10-3log10(Th)

(4)

得到閾值Th*。整個閾值算法框圖以及算法偽代碼如圖3。

圖3 閾值算法框圖及算法偽代碼。(a)閾值算法框圖；(b)閾值算法偽代碼

2.2.2 信道選取與數(shù)據(jù)比較

首先選取一個短時 “寧靜” 的觀測頻段作為閾值計算數(shù)據(jù)樣本，從表2中選取此次檢測組數(shù)最多的一次作為信道的初選，西向測試在第4次04月27日時共采集了2 219組數(shù)據(jù)，從西向第4次數(shù)據(jù)的頻譜瀑布圖中選取短時 “寧靜” 的觀測頻段，頻譜圖如圖4～圖6； 其次在上述數(shù)據(jù)樣本中通過閾值算法計算得到閾值，以此閾值遍歷其余觀測時段數(shù)據(jù)，對比閾值，超過閾值時間和低于閾值時間進行對比，求得該時段的干擾信號所占比例，對比結果如表3；最后利用算法計算得到的閾值對本文遴選的7個5 MHz無線電干擾較少的無線電通帶進行評估。

圖4 L波段頻譜瀑布圖Fig.4 L-band spectral waterfall

圖4中橫坐標為L波段1～2 GHz，縱坐標為數(shù)據(jù)采樣的組數(shù)，豎軸為采樣的信號強度，因為能較為直觀地從圖中看出本文所需頻段的初步干擾情況，故此次L波段的所有測試均采用此種繪制方法。

圖5 L波段頻譜右視圖Fig.5 L-band spectrum right view

圖5為L波段頻譜瀑布圖的右視圖，圖中能清晰看到所有測試組數(shù)的信號強度，且能直觀看到本文所需的頻段的干擾情況，能為本文尋找短時 “寧靜” 的觀測頻段確定初步篩選工作。圖6為L波段頻譜瀑布圖的俯視圖，圖中標記北斗 B1, B2和B3頻點及全球定位系統(tǒng) L1和L2頻點位置，為后續(xù)通帶的遴選起輔助作用。

圖6 L波段頻譜俯視圖Fig.6 L-band spectrum vertical view

通過閾值算法對數(shù)據(jù)閾值的確定，實現(xiàn)對整個L波段監(jiān)測數(shù)據(jù)遴選的多個5 MHz無線電干擾較少的無線電通帶中，最終確定介于北斗系統(tǒng)及全球定位系統(tǒng)頻點附近的7個無線電通帶。信道選取流程如圖7。

本次測試主要是為了遴選介于北斗 B1, B2和B3頻點，全球定位系統(tǒng) L1和L2頻點之間7個5 MHz無線電干擾較少的無線電通帶，首先從圖4～圖6中L波段頻譜圖中初選出5 MHz干擾較少的無線電通帶1 161～1 165 MHz頻段作為初選閾值確定的可用信道，再根據(jù)上述算法確定閾值后，通過閾值遍歷對比程序計算考察信道的潔凈率，得到各組信道監(jiān)測數(shù)據(jù)的潔凈率如表3。從表3可以看出，在16次L波段測試中所遴選的5 MHz干擾較少的無線電通帶1 161～1 165 MHz頻段信道潔凈率都在95%以上，符合太陽觀測需求，所以本文選用此閾值為所遴選的7個5 MHz干擾較少的無線電通帶的閾值。

本文結合圖5～圖7以及表3初步遴選介于北斗 B1, B2和B3頻點以及全球定位系統(tǒng) L1和L2頻點之間7個5 MHz干擾較少的無線電通帶如圖8，從圖8看到本文所選的通帶具體的干擾情況：圖中B2a(1 161～1 165 MHz)和B2b/L2a(1 221～1 225 MHz)為B2頻點附近通帶；B2b/L2a(1 221～1 225 MHz)和L2b(1 231～1 235 MHz)為L2頻點附近通帶；B3a(1 246～1 250 MHz)和B3b(1 291～1 295 MHz)為B3頻點附近通帶；L1a(1 551～1 555 MHz)和L1b(1 596～1 600 MHz)為B1/L1頻段附近通帶。利用上述所求閾值確定遴選的7個5 MHz干擾較少的無線電通帶的潔凈率如表4。

圖7 遴選信道流程框圖Fig.7 Block diagram of the selection channel flow

表3 1 161～1 165 MHz頻段潔凈率Table 3 Cleanliness rates in the 1 161-1 165 MHz band

圖8 L波段頻譜遴選通帶Fig.8 L-band spectrum selection channel

表4 遴選通帶潔凈率Table 4 Selection belt cleanliness rate

由表4可知，北斗系統(tǒng)B1、全球定位系統(tǒng)L1頻點遴選通帶為1 551～1 555 MHz和1 596～1 600 MHz，潔凈率為98.329%和98.301%；北斗系統(tǒng)B2頻點遴選通帶為1 161～1 165 MHz和1 221～1 225 MHz，潔凈率為98.315%和98.335%；北斗系統(tǒng)B3頻點遴選通帶為1 246～1 250 MHz和1 291～1 295 MHz，潔凈率為98.224%和97.650%；全球定位系統(tǒng)L2頻點遴選通帶為1 231～1 235 MHz，潔凈率為98.260%，且所遴選通帶的潔凈率都大于95%以上，符合太陽觀測需求。

3 總結與展望

針對即將在云南天文臺建設的L波段精密流量太陽射電望遠鏡，本文提出一種基于Simple Thresholding算法和CUSUM算法的改進閾值算法，并對臺址的L波段無線電測試和觀測頻段遴選出介于北斗 B1, B2和B3頻點，全球定位系統(tǒng)L1和L2頻點之間7個5 MHz干擾較少的無線電通帶，分別為1 551～1 555 MHz, 1 596～1 600 MHz, 1 161～1 165 MHz, 1 221～1 225 MHz, 1 246～1 250 MHz, 1 291～1 295 MHz和1 231～1 235 MHz，其潔凈率分別為98.329%, 98.301%, 98.315%, 98.335%, 98.224%, 97.650%和98.260%，均符合太陽觀測需求，為下一步接收機的設計和信號處理提供了依據(jù)。