基于改進連通區域標記的跳頻信號分選識別 *

郭昭藝，黃祥，孟悅，杜彪，霍丹江

（江蘇方天電力技術有限公司，江蘇 南京 211102）

0 引言

近年來，由于中國低空空域的逐步打開，無人駕駛飛行器也逐步地由軍事領域延伸至民用領域［1］，但在提供方便的同時，卻也造成了“黑飛”“濫飛”等事故時有發生，這將會對電力設施等設備的安全平穩運營造成很大的安全威脅［2-3］。所以，很有必要進行中低空領域反無人機攻擊的關鍵技術研發。

反無人機入侵技術的關鍵是對無人機信號開展偵察識別。為了增加傳輸的穩定性，小型商用和民用無人機通常采用ISM（industrial scientific medical band）頻段的跳頻信號作為遙控指令來進行飛行控制。該頻段中包含的通信信號復雜而擁擠，會形成較多的噪聲和干擾，如何實現復雜電磁環境下跳頻信號的分選識別，進而實現對“黑飛”無人機的檢測，成為反制無人機亟待解決的重難點問題。

當前，一些研究者已經對跳頻信息的測量和分選識別等問題開展了研究。Chung C. D. 等在掌握信號先驗知識條件下對跳頻信號局部參量作出了估算［4-5］，不能達到通信對抗的需求。針對非平穩的跳頻信號，大多數采用時頻分析技術來實現參數盲估計。Barbarossa S. 等采用不同的時頻分析技術獲取相關的峰值曲線來實現參考估算［6-7］，但是該技術需要比較高的信噪比。為此，一些學者給出了一個采用時頻脊線的跳頻時間信息參數估算方案［8-9］，可以在較低峰值信噪比條件下實現參考估算。但是，這些方案均針對接收信道中僅存在一條跳頻信息的特殊情形，而在出現定頻、突發等其他干擾情況時方式均無效。一些學者也從圖像處理的角度研究跳頻信號檢測識別問題［10-13］。文獻［10］對時頻圖像不同尺寸的結構元素完成形態學處理，從而消除了干擾，但該方法必須選定結構元素的尺寸。文獻［11］使用了一種二維模板匹配的方式進行指定跳頻信號的檢測識別，但虛警率較高。文獻［12-13］中將信號時頻矩陣轉換為灰度時頻圖像，經閾值化后采用形態學濾波消除定頻干擾等影響，從而完成了跳頻信號參數估計。以上方式都需要經過人工干預設置合適的門限，僅能估計跳頻信號部分參數，在復雜干擾環境下的檢測識別效果不佳。

本文采用時頻圖像連通區域聚類方法解決復雜電磁環境下混合信號中的跳頻問題。首先，采用基于能量統計的自適應閾值噪聲去除方法去除信號時頻 圖中的背景噪聲；然后，通過連通區域標記將信號區域聚集，接著對各區域進行參數提取、聚類和統計，從而實現了復雜電磁環境下跳頻信號的分選。

1 時頻信號的預處理

1.1 跳頻信號時頻分析

跳頻信號的頻點是隨著時間而不斷變化的，因此可采取時頻分析的方法來描述信號在頻域中的變化特征。本文采用短時傅里葉變換法，該方法不會產生交叉特征干擾，且計算工作量也較小。

對于跳頻信號，接收端和發送端信號的關系可表示為

式中：x(t)為無人機發射的跳頻信號；s(t)為偵察設備接收到的信號；n(t)為噪聲。

對接收信號s(t)進行短時傅里葉變換，可表示為

式中：h(t)為時間長度很短的窗函數，這里取矩形窗。經過短時傅里葉轉換，將接收信息s(t)轉換成時間-頻率二維函數，從而形成二維時頻圖。

1.2 基于局部窗口能量門限統計的去噪

ISM 頻段是一種開放的頻段，接收信號中除了無人機跳頻遙控信號外，還混有大量的噪聲和干擾（定頻干擾和斜變信號干擾）。

對接收信號進行短時傅里葉轉換，可以分別獲得如圖1 和圖2 所給出的在無噪聲和有噪聲條件下的信號時頻分布圖。可以看到，噪聲對目標信號的時頻圖影響較大，會影響后續的信號參數估計和信號分選，需要進行降噪處理。

文獻［14］采用基于能量統計的自適應門限去噪方法，該方法用整個時頻分布能量進行統計來估計能量門限，適用于噪聲均勻分布的情況。由于噪聲在不同時間和頻段存在不同分布，因此可采用基于局部窗口能量門限統計的自適應去噪方法。

根據噪聲強度將信號二維時頻能量分布劃分為K個區域，設第i個區域的時頻能量分布矩陣為Ti(x，y)，則局部區域時頻能量分布矩陣的平均能量為

式中：Mi和Ni分別表示局部區域i的頻域點個數和時域點個數。不同區域的去噪門限可表示為

式中：ai為區域i的去噪系數，表示該區域的平均能量的倍數。去噪后的二維時頻分布可表示為

式中：K為局部區域數。在信號二維時頻能量分布中，由于噪聲能量小于目標時頻點能量，且數量遠大于目標信號時頻點，因此可以通過遍歷統計法計算合適的去噪門限，具體步驟如下：

步驟1： 在不同區域中分別對去噪系數ai從0到10 等間隔取值，取值間隔為0.1，然后根據式（4）得到不同區域下的去噪門限thi(k)。

步驟2： 統計不同門限下能量大于閾值點的數目c(k)，可表示如下

步 驟3： 找 出c(k)，k= 1，2，…，Nk下 降 最 大 處的點，即c(k+ 1)較c(k)有最大的下降，則c(k+ 1)對應的thi(k+ 1)為最終的去噪門限thi。

步驟4： 按式（5）計算，得到去噪后的時頻分布，即Fi(x，y)。

2 基于連通區域標記的參數提取

2.1 連通區域標記

去噪后的時頻分布仍然混雜著定頻、斜變等干擾信號。從時頻分布中可以發現，信號分布呈現規則區域的連通特性，每個連通區域的屬性可以很好地表征信號的參數特性，對不同連通區域進行聚類則可以實現信號的分選。

本文采用基于區域生長思想的8 鄰域連通方法，常見的8 鄰域連通如圖3 所示。

圖3 8 鄰 域 連 通 圖Fig. 3 Eight neighborhood connectivity

文獻［15］中使用了二值化的時頻矩陣作為連通算法的輸入，該方式會損失信號的幅度信息，尤其當多目標信號混疊時，并不能充分利用幅度的變化信息。因此，本文中使用了含有一定幅度原始信息的時頻矩陣F作為輸入矩陣，并構建連通標記矩陣B，其大小與F相同，構建隊列queue矩陣并初始化，初始化標記計數label_count=0。

算法實現流程如下：

輸入：去噪后的時頻矩陣F(x，y)，x和y分別為時間軸（橫坐標）和頻率軸（縱坐標）上的坐標。

輸出：標記矩陣B(x，y)。

步驟1： 先依次掃描A，當掃描到非零點p且未被標記時，即label_count=label_count+1 時，在B上將label_count的數值賦于p位置，同時掃描p的鄰域點，如果出現了未被標記的非零點則在B中加以標注，并置于queue中，以作為區域內生長的種子。

步驟2： 若queue不為空時，從queue中取出點p1，掃描p1的8 鄰域點，若存在未被標記的非零點，則在B中進行標記將lobel_count的值賦予該點，并放入queue中。

步驟3： 如果queue為非空矩陣，在queue中提取點p1，并掃描p1的鄰域點，如果出現了未被標記的非零點，則在B中進行標記后將lobel_count的值賦于該點，并存放在隊列queue矩陣中。

步驟4： 當queue為空時，一個連通區域標記完成。

步驟5： 直至整個矩陣掃描完畢，則得到標記矩陣B。

連通標記處理后，2 類矩陣示意圖如表1 所示。

表1 時頻矩陣F 和標記矩陣B 示意圖Table 1 Time-frequency matrix F and label matrix B

2.2 時頻參數提取

在對信號時頻分布進行連通標記后，可以對跳頻周期、中心頻點、帶寬等參數進行提取。傳統的時頻參數提取方法是對連通區域構造最小矩形邊界，矩形邊界可用向量(x，y，L，H)表征，x和y表示左下角的橫、縱坐標，L和H表示x軸長度和y軸長度，則信號時頻數據可以提取為起始時間x、結束時間x+L、跳頻周期L、中心頻點y+H/2、帶寬H，如圖4所示。

圖4 信號時頻參數提取示意圖Fig. 4 Extraction of time-frequency parameters of signal

2.3 連通區域優化

信號時頻分布的連通區域形成后，由于還存在定頻、斜變干擾，以及去噪引起的信號連通斷裂的情況，影響后續的時頻參數提取，因此還需要對連通區域進一步優化。

（1） 連通區域片段關聯

在信噪比較低的條件下，采用基于能量統計的自適應閾值噪聲去除方法會使部分信號區域的幅度置0，從而造成同一類信號連通區域的斷裂，需要根據斷裂片段的特征相似性對其進行關聯，然后完成拼接。源于同一類信號的斷裂區域關聯條件包括中心頻率的相似性、帶寬的相似性和時間接續性（即前一段信號連通區域的結束時刻與后一段信號連通區域的開始時刻在一個較小的范圍內）。

設第i個信號段的中心頻率為yi+Hi/2，終止時刻為xi+Li，帶寬為Hi；第j個信號段的中心頻點為yj+Hj/2，起始時刻為xj，帶寬為Hj，滿足前后信號連通區域片段的關聯條件為

關聯門限設置為Δt= 10，ΔH= 1，Δf= 5，將符合關聯條件的信號連通片段之間的時頻幅度設置為前一信號片段終止時刻幅度和后一信號片段起始時刻幅度的均值，可表示為

完成斷裂區域的幅度填充后，則完成了不同連通區域片段的關聯和拼接。

（2） 連通區域的干擾抑制

由于信號時頻分布中存在定頻、斜變信號干擾，容易與目標信號（跳頻信號）發生碰撞，導致連通區域形狀的畸變。干擾和跳頻信號在連通區域標記圖上的分布如圖5 所示，最終產生較大的時頻參數估計誤差。

圖5 干擾和調頻信號連通區域標記分布圖Fig. 5 Label distribution map of interference and frequency-hopping signals connected region

觀察定頻和斜變信號干擾的分布特征可以發現，2 類干擾與目標信號的連通區域標記分布有較大差異。其中，定頻信號干擾在固定頻率范圍內的較長時段內連續分布，即時間軸占比率非常高，而跳頻信號時間軸占比率較低；斜變信號干擾在連通區域形成的最小矩形邊界中占零比非常高，而跳頻信號占零比相對較低。基于以上分析的特征差異，可對連通區域標記分布圖進行干擾抑制。

1） 定頻干擾抑制

對所有連通區域進行參數提取，當提取的跳頻周期滿足Li≥Lmax（Lmax為已知跳頻信號數據庫中的最大跳頻周期）時，則判斷該連通區域為定頻干擾，將其所在區域置0。

2） 斜變信號干擾抑制

首先對所有連通區域計算整體占零比，如果大于一定范圍時，則該連通區域內可能存在斜變信號，可表示為

式中：Δn= 20%。

對滿足式（9）的連通區域計算每一行的占零比，行占零比大于80%的判為斜變信號，將非0 位置的信號置0。

2.4 多目標信號混疊下的參數提取

由于不同跳頻信號的出現時刻（起始時刻）、跳頻周期、中心頻率和帶寬等參數的差異，可能在時頻分布上存在時間和頻率維的混疊，造成連通區域的變形，使得參數提取發生非常大的誤差，其示意如圖6 所示。

圖6 多跳頻信號混疊下的連通區域標記分布圖Fig. 6 Label distribution map of connected regions under aliasing of multi-frequency-hopping signals

雖然2 個跳頻信號在時間和頻率維度上發生了混疊，但由于去噪后的時頻分布矩陣保留了信號幅度信息，且混疊區2 個跳頻信號幅度明顯大于非混疊區的單個信號幅度，因此可以利用幅度的突變性來重構不同跳頻信號的連通區域，從而分離出不同跳頻信號，更準確地提取信號參數。多目標信號混疊下的參數提取步驟如下：

（1） 根據連通區域標記分布圖確定最小矩形邊界，并按照2.2節的方法提取參數向量(x0，y0，L0，H0)。

（2） 逐行計算后一時刻和前一時刻的幅度差，當幅度差有明顯增大時，將后一時刻的點設置為混疊列邊界點，可表示為

式中：σA為連通區域最小矩形邊界內非零幅度值的均方根誤差；k為尺度因子，這里取3。將滿足式（10）的(x+ 1，y)點設置為混疊邊界點（列方向），每一行計算完成后，得到混疊列邊界的所有點。

（3） 逐列計算后一頻率點與前一頻率點的幅度差，當幅度差有明顯增大時，將后一頻率點的點設置為混疊行邊界點，可表示為

將滿足式（11）的(x，y+ 1)點設置為混疊行邊界點（行方向），每一列計算完成后，得到混疊行邊界的所有點。

（4） 完成混疊信號的連通區域重構后，重新進行連通區域標記，分離2 個連通區域，最后將帶有幅度信息的二維時頻圖像進行二值處理，并根據二值圖像信息提取時頻參數，示意如圖7 所示。

圖7 連通區域重構后的參數提取示意圖Fig. 7 Parameter extraction after reconstruction of connected regions

3 仿真結果及分析

3.1 信號分選

由于不同類型無人機采用的遙控信號（跳頻信號）具有不同的跳頻周期和調制帶寬，因此，在對所有連通區域提取信號時頻參數后，通過DBSCAN（density-based spatial clustering of applications with noise）聚類算法實現信號分選。該算法是一種經典的密度聚類算法，通過把類描述為緊密連接的最大節點的最大集合，把具有足夠高密度區域以外的任意形狀劃分為一類［15］，其優勢是不需要預設類別，運用該算法完成信號跳頻周期和調制帶寬參數的分選。

除了跳頻周期和調制帶寬可以作為分選識別的參數，不同跳頻信號的出現時刻也是分選識別的重要參數。尤其是在多個同類型無人機目標條件下，信號具有相同的跳頻周期和帶寬，但信號出現時刻不同，因此在利用跳頻周期和帶寬完成聚類后，再利用信號出現時刻進行聚類，仍采用DBSCAN聚類方法。

3.2 仿真條件

（1） 場景1

跳頻信號1（某型無人機遙控信號）：跳頻周期為1 ms，跳頻頻率集{50，75，100，125，80，125，100，75，50，75，100，125，80，125，100，75，50，75，100，125} MHz，BPSK 調制帶寬為1.8 MHz；

跳頻信號2（其他外來無人機）：跳頻周期為2 ms，跳頻頻率集為{190，70，140，120，155，30，60，115，90，180} MHz，BPSK 調制帶寬為2 MHz。

干擾及噪聲

定頻干擾參數包括：干擾頻率分別為143 MHz和162 MHz，帶寬分別為0.8 MHz 和4 MHz。掃頻干擾：起始頻率為20 MHz，跳頻斜率20 GHz/s，重復周期為250 Hz。噪聲：高斯噪聲，信噪比為6 dB。

對信號作短時傅里葉變換，使用窗長為4 096的Hamming 窗，1 024 為步長，畫出時頻圖。信號如圖8，9 所 示，圖8 為 跳 頻 信 號1 的 時 頻 圖，圖9 為 跳頻信號2 的時頻圖。

圖8 跳頻信號1 時頻圖Fig. 8 Time-frequency diagram of frequency-hopping signal 1

圖9 跳頻信號2 時頻圖Fig. 9 Time-frequency diagram of frequency-hopping signal 2

圖10 為加入定頻、掃頻干擾得到的時頻分布，圖11 為混入噪聲的時頻分布。使用基于局部窗口能量統計的去噪方法，得到去噪后的時頻分布如圖12 所示。

圖10 混合信號時頻圖Fig. 10 Time-frequency diagram of mixed signal

圖11 含噪信號時頻圖Fig. 11 Time-frequency diagram of noisy signal

圖12 去噪后的時頻圖Fig. 12 Time-frequency diagram after denoising

對去噪后的時頻矩陣進行連通區域標記，形成初步的連通區域標記分布圖，然后進行連通區域優化，完成連通區域片段關聯和拼接、干擾抑制、多目標信號混疊下連通區域優化，得到如圖13 所示的修正后的二值化連通區域標記圖。

圖13 修正后的連通區域標記圖Fig. 13 Label map of connected regions after correction

針對每個連通區域提取其跳頻周期、調制帶寬和出現時刻，參數歸一化后運用DBSCAN 算法進行聚類。鄰域內樣本分布的密集程度參數(ε，MinPts)設置為(0.02，1)，表示以聚類中心半徑為0.01 的鄰域內點的個數不少于1 個的聚類，根據聚類結果將信號分為2 類，如圖14 所示。2 類信號的參數估計如表2 所示。

表2 兩組連通區域參數統計值Table 2 Parameter statistics of two groups of connected regions

圖14 帶類別信息的時頻矩陣連通區域標記圖Fig. 14 Label map for connected regions of timefrequency matrix with class information

此時與原始信號參數設置進行對比，可以采用決策樹的思想對每個參數進行判斷，從而判斷分選出跳頻1 信號和跳頻2 信號。

（2） 場景2

跳頻信號1：跳頻頻率集、跳頻周期和帶寬同場景1；

跳頻信號2：跳頻周期和帶寬同場景1，跳頻頻 率 集 為{190，74，140，124，155，30，60，124，74，180} MHz。

在跳頻信號中加入定頻、掃頻干擾和噪聲（參數同場景1）得到的時頻分布如圖15 所示，分別采用基于局部能量統計方法去噪，采用本文提出的改進連通區域標記方法完成連通區域的優化和參數提取，最后采用DBSCAN 算法進行聚類，得到如圖16，17 所示的連通區域標記圖和類別信息。2 類信號的參數估計如表3 所示（這里僅列舉平均跳頻周期和平均帶寬參數的統計值）。

表3 2 組連通區域參數統計值Table 3 Parameter statistics of two groups of connected regions

圖15 跳頻+干擾+噪聲混合信號時頻圖Fig. 15 Time-frequency diagram of mixed signal containing frequency-hopping， interference，and noisy signals

圖16 修正后的連通區域標記圖Fig. 16 Label map of connected regions after correction

圖17 帶類別信息的時頻矩陣連通區域標記圖Fig. 17 Label map for connected regions of timefrequency matrixes with class information

（3） 場景3

跳頻信號1（某型無人機遙控信號）：跳頻頻率集、跳頻周期和帶寬同場景1；

跳頻信號2（外來無人機遙控信號1）：跳頻頻率集、跳頻周期和帶寬同場景1；

跳頻信號3（外來無人機遙控信號2）：跳頻周期和帶寬同跳頻信號2，跳頻頻率集為{165，130，75，110，155，45，80，135，60，120} MHz。

在跳頻信號中加入定頻、掃頻干擾和噪聲（參數同場景1）得到的時頻分布如圖18 所示。

圖18 跳頻+干擾+噪聲混合信號時頻圖Fig.18 Time-frequency diagram of mixed signal containing frequency-hopping，interference， and noisy signals

采用本文提出的改進連通區域標記方法完成連通區域的優化和參數提取，最后采用DBSCAN 算法進行聚類，得到如圖19，20 所示的連通區域標記圖和類別信息。

圖19 修正后的連通區域標記圖Fig. 19 Label map of connected regions after correction

圖20 帶類別信息的時頻矩陣連通區域標記圖Fig. 20 Label map for connected regions of timefrequency matrixes with class information

3.3 不同方法的性能比較

2 種方法比較：基于連通區域標記的信號分選法［16］和本文提出的改進連通區域標記分選法。

基于連通區域標記的信號分選法：基于整體能量統計的自適應門限去噪+常規連通區域標記+信號拼接和拆解（連通區域修正）。

本文提出的改進連通區域標記分選法：基于局部能量統計的自適應門限去噪+常規連通區域標記+改進的信號拼接和拆解+信號混疊條件下的連通區域修正。

比較的性能指標包括參數（平均跳頻周期、平均帶寬）估計誤差和目標正確識別概率。

采用1 000 次蒙特卡羅仿真統計平均跳頻周期、帶寬參數估計誤差，并在不同信噪比條件下統計目標無人機遙控信號（跳頻信號2）的正確識別概率。2 類方法的參數估計誤差如表4 所示。

表4 不同場景下2 種方法的參數估計誤差比較Table 4 Parameter estimation errors of two methods in different scenarios

從表4 結果可以看出，3 個場景下，本文提出的改進連通標記方法提取的參數相對誤差明顯小于常規連通標記法。其中，場景1 下的性能提升是由于改進的連通標記法進行了連通區域片段關聯和干擾抑制，使得整個時頻分布的連通圖更完整和干凈。場景2 和3 下的性能提升是由于采用了基于時頻幅度差異的連通區域重構，將混疊在一起的跳頻信號進行了連通邊界區分，從而使參數估計更加精確。

圖21 為多目標跳頻信號混疊（場景2 和3）環境下，采用1 000 次蒙特卡羅仿真，不同信噪比的類別2 跳頻信號（外來無人機遙控信號1）的識別概率。可以看到，信噪比對跳頻信號識別影響明顯，識別概率隨著信噪比增加而提高。常規連通標記方法下，為達到高識別概率（0.9），信噪比分別需達到5 dB 和7dB 以上，而采用了改進連通區域標記方法后，信噪比只需3 dB 和4.2 dB 以上。且由于參數估計精度更高，在不同信噪比條件下，改進連通區域標記方法的目標識別概率均高于常規方法，能夠更好地識別出外來無人機的遙控信號。

圖21 類別2 跳頻信號的識別概率Fig. 21 Recognition probability of Class 2 frequency-hopping signal

4 結束語

本文針對干擾環境下無人機跳頻遙控信號的分選識別問題開展了研究。首先考慮到噪聲環境的不均勻性，設計了基于局部能量門限統計的去噪方法。然后對去噪后時頻矩陣進行連通區域標記，在常規連通區域標記方法的基礎上，對連通標記方法進行了改進。根據連通區域片段的相關性設計了片段關聯方法，完成信號的拼接。根據定頻、掃頻干擾的時頻分布特征設計了干擾抑制方法，并在多目標跳頻信號混疊情況下提出了基于幅度差異的連通區域重構。最后根據優化的連通區域標記實現對跳頻信號的參數估計和分選識別。仿真實驗表明，本文方法在無多目標混疊和有多目標混疊2 種場景下，信號時頻參數估計精度明顯高于常規連通標記方法，對目標跳頻信號的正確識別率也明顯提高，能夠較好地檢測識別外來無人機。