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

楊 航，董 亮，何樂生*

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

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

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

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

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

1 系統組成

1.1 系統組成

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

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

圖1 測試系統Fig.1 Test system

1.2 測試設備及參數

HL050天線應用廣泛，主要參數見表1，其優點是結構簡單、頻帶寬、功率容量大、調整和使用方便。

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

表1 測試設備及參數Table 1 Test equipment and parameters

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

1.3 監測過程

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

2 數據處理

本節主要介紹采集數據的格式、數據閾值遴選和本次最終信道的遴選過程。

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

2.1 數據介紹

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

圖2 測試系統實測圖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)頻點，全球定位系統L1(1 574.397～1 576.443 MHz)和L2(1 226.577～1 228.623 MHz)頻點之間7個5 MHz干擾較少的無線電通帶。已知全球定位系統中L1頻點在北斗B1頻點范圍內，且B2與L2頻點之間相差不超過10 MHz，無法遴選兩處相對潔凈的5 MHz通帶，B2右邊遴選的通帶與L2左邊遴選的通帶重合，一共7個通道(經過無線電環境監測帶內無干擾信號)，全球定位系統和北斗系統的所有通信頻點都穿插在監測頻點之間[17-20]。

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

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

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

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

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

2.2.1 閾值計算方案

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

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

(1)

樣本方差是表示射頻干擾的統計異常的重要估計值。由于絕對中位差(Median Absolute Deviation, MAD)和中位數是剔除異常值的最有效方法，因此，本文提出一種基于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算法使用的方差估計值，中位數是原始總體數組xi的中位數。然后用此中位數減去原始樣本的每個元素，創建與原始元素大小相同的新修改樣本。然后計算這個新樣本的中位數，并用恒定的比例因子1.428 6相乘[24]，使該估計與預期的高斯分布一致。得到線性計算的閾值Th之后，需要將閾值取對數轉換在監測的數據中計算通帶的潔凈率，通過

Th*=10-3log10(Th)

(4)

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

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

2.2.2 信道選取與數據比較

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

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

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

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

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

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

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

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

本文結合圖5～圖7以及表3初步遴選介于北斗 B1, B2和B3頻點以及全球定位系統 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可知，北斗系統B1、全球定位系統L1頻點遴選通帶為1 551～1 555 MHz和1 596～1 600 MHz，潔凈率為98.329%和98.301%；北斗系統B2頻點遴選通帶為1 161～1 165 MHz和1 221～1 225 MHz，潔凈率為98.315%和98.335%；北斗系統B3頻點遴選通帶為1 246～1 250 MHz和1 291～1 295 MHz，潔凈率為98.224%和97.650%；全球定位系統L2頻點遴選通帶為1 231～1 235 MHz，潔凈率為98.260%，且所遴選通帶的潔凈率都大于95%以上，符合太陽觀測需求。

3 總結與展望

針對即將在云南天文臺建設的L波段精密流量太陽射電望遠鏡，本文提出一種基于Simple Thresholding算法和CUSUM算法的改進閾值算法，并對臺址的L波段無線電測試和觀測頻段遴選出介于北斗 B1, B2和B3頻點，全球定位系統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%，均符合太陽觀測需求，為下一步接收機的設計和信號處理提供了依據。