基于HMM 的逆雷達輻射源狀態識別推理方法

朱夢韜，張露瑤，李瑞，楊靜

（1.北京理工大學 信息與電子學院, 北京 100081；2.電子信息系統復雜電磁環境效應國家重點實驗室, 河南, 洛陽471003；3.電磁空間認知與智能控制技術實驗室, 北京 100191）

近年來，隨著人工智能等先進技術在雷達對抗場景中的雷達方和干擾方系統中獲得越來越多應用，對抗雙方系統的精準對抗博弈能力不斷提升.通過引入機器學習類算法模型，干擾方可以對多功能雷達（multi-function radar, MFR）的工作狀態以及動態特性實現快速準確識別與估計，用以支撐后續干擾效果評估和自適應干擾決策[1].典型研究如基于觀測雷達狀態序列構建反映雷達各層次行為特征的狀態識別模型，如隱馬爾可夫模型（hidden Markov models，HMM），深度學習模型等，實現對雷達行為意圖的識別與動態特性估計[2-11]；雷達方則往往基于機器學習框架構建各類神經網絡模型實現對復雜干擾信號的高效精準識別，從而有效提高抗干擾性能[12-18].對抗雙方系統都亟須尋求根據所掌握先驗信息及部分可觀測結果獲取更豐富和更準確的對手系統信息的技術手段，以便在對抗博弈過程中取得優勢.因此，針對無法直接觀測得到的對手系統內部細節信息（如干擾方狀態識別的性能）進行逆向推理技術研究開始進入研究者的視線.

針對對手系統的逆向推理研究是一個逐步擴展深入的過程.最初是根據給定的系統結構及決策（控制）策略逆向分析出反映系統動作的最優性能指標值[19].隨后研究者開始關注從對手系統外在動作層面對其內部系統信號處理環節狀態進行推理，有研究假設在對抗場景中的我方狀態轉移概率模型及給定對方策略下的我方觀測到動作概率模型、對方觀測似然概率模型均已知，并據此假設通過已知的我方系統信息及可觀測的對方系統信息逆向分析得到對方對我方的估計結果[20].同樣基于上述概率模型均已知的假設條件下，最近擴展到了更深層次的控制策略（目標函數）[21]的逆向推理，即針對對方系統中效用函數進行逆向分析的研究.李瑞等[22]構建了雷達對抗場景中的逆向推理總體框架，并展現了涵蓋雙方系統各級功能模塊的逆向分析交互過程.

針對自主系統的逆向推理實現方法研究方面，KRISHNAMURTHY 等[20]學者將具有認知能力的無人機或雷達系統作為對手，提出利用逆濾波方法估計對方系統在觀測條件下對我方系統狀態的后驗分布，并根據估計結果設計了有效對抗自主系統的反自主系統.MATTILA 等[23]基于構建的隱馬爾可夫模型結合貝葉斯濾波器輸出的后驗分布序列，采用逆濾波方法逆向推理出對方傳感器精度及所觀測序列.文獻[24 - 26]引入在針對智能體行為（策略）逆向推理中廣泛應用的逆強化學習技術，利用所觀測到的智能體行為軌跡或決策策略逆向分析智能體效用函數.KRISHNAMURTHY 等[27]從干擾方角度結合顯示偏好理論證明了利用逆濾波方法能夠估計雷達方的估計結果和傳感器精度以達到有效逆跟蹤的目的.

綜上所述，逆向推理技術具有在非合作條件下分析發覺對手系統內部更多有效信息的潛力，但已有研究[20]需假設非合作系統內部的各個處理環節已知，從而逆向推理技術轉變為根據觀測對非合作系統中變量節點的后驗計算問題.本文從提升雷達系統抗干擾能力出發，考慮干擾方具備根據觀測得到的雷達信號序列進行雷達工作狀態識別及預測，并基于識別預測結果調整其干擾動作的能力.由于非合作系統內部處理環節[20]實際對我方未知，本文將考慮通過構建隱變量的方式，從可觀測雷達和非合作方交互序列入手，設計了對干擾系統對我方雷達輻射源狀態識別結果的逆向估計方法，從而獲取干擾方狀態識別模塊對我方雷達狀態識別的準確度.該逆向分析結果也可用于輔助對未來時刻干擾方動作的預測等后續任務.雷達對抗仿真實驗結果驗證了所提方法的可行性和準確性.

1 雷達對抗場景中的逆狀態識別模型

本文結合文獻[22]中構建的雷達對抗場景下的逆向推理框架，以及文獻[2 - 11]中干擾方基于觀測雷達信號序列對雷達工作狀態進行識別的研究成果，從雷達方視角對逆信號處理模塊中逆輻射源狀態識別問題進行建模與求解.基于MFR 前后工作狀態間高度相關性使得干擾方對雷達工作狀態識別的結果序列前后狀態間也具有潛在相關性的時序特性，雷達方可將干擾方用來控制干擾發射信號的干擾方對雷達方真實狀態的識別結果變量視為服從馬爾可夫過程的隱變量，將雷達方對干擾方通過工作狀態識別后采取的干擾樣式輸出動作的識別結果作為雷達方可觀測變量.由此雷達可構建干擾系統內部的“雷達狀態識別結果（隱藏層）-雷達觀測干擾動作（觀測層）”的HMM 模型.具體的，雷達方對干擾方輻射源狀態識別環節逆向推理任務圖如圖1 所示.

圖1 雷達方對干擾方輻射源狀態識別逆向推理任務圖Fig.1 Flow chart of inverse inference for radar to recognize the jammer emitter work mode

上述干擾方和雷達方之間的交互過程的概率模型如下：

式中：k=1,2,···,N表示離散化時間；p(·)表示條件概率密度函數或條件概率質量函數； ～表示服從的分布.其中k時刻MFR 信號對應的真實工作狀態為qk，在干擾方視角通常假定qk的轉移具有一階馬爾可夫性.干擾方基于觀測信號yk完成對k時刻MFR 的真實工作狀態的識別，為識別結果.由干擾方的狀態識別模型決定，為服從某種概率分布的算子,該概率密度分布對應了狀態識別模型的混淆矩陣性質.例如在現有典型干擾方研究中[2-11]，干擾方狀態識別模型輸入雷達脈沖序列，輸出識別得到的雷達工作狀態概率分布，此時.為描述簡便，本文建模時不考慮干擾方機動對我方雷達狀態切換的影響,該假設對本文提出的方法沒有影響.

本文設計的逆狀態識別過程不需要P,E,G的先驗知識，通過對離散的和構建離散隱變量狀態空間HMM 模型，利用隱變量和狀態轉移矩陣、觀測矩陣分別表征交互過程中的中間變量和條件概率分布，憑借觀測數據學習參數P,E,G,具體如下：

設定n為時間序列長度，雷達通過對干擾方實施的真實干擾樣式動作序列U={u1,u2,···,un}進行干擾樣式識別后得到的動作觀測序列O={o1,o2,···,on}估計干擾方對雷達真實工作狀態序列Q={q1,q2,···,qn}的識別結果，得到逆狀態識別序列

針對上述逆狀態識別問題，構建雷達方基于HMM的逆輻射源狀態識別模型，該HMM 記為λ={A,B,π}，其中：

（2）雷達方所觀測到的干擾樣式數目為M，干擾樣式集合記為V={V1,V2,···,VM}.

（3）干擾方狀態識別模型中的狀態轉移概率矩陣為A=[aij]，其中aij=,i,j=1,2,···,N，用于描述干擾方對MFR 工作狀態識別模型中估計得到的各雷達狀態間的轉移概率.

（4）觀測特征概率分布矩陣為B=[bij],bij=P(ok==Si),i=1,2,···,N,j=1,2,···,M，用于在干擾方狀態識別結果條件下雷達觀測到干擾樣式ok的概率.

（5）初始狀態概率分布為π=[πi] ,πi=,i=1,2,···,N表示初始時刻干擾方狀態識別結果的概率分布.

此外，記真實MFR 工作狀態集合為D={D1,D2,···,DN}，雷達采取的抗干擾措施序列為Z={z1,z2,···,zn}，其中z∈W，W={W1,W2,···,WL}為雷達抗干擾措施集合，L為抗干擾措施總數.

2 HMM 模型下的逆狀態識別推理算法

本節基于上述逆輻射源狀態識別HMM 模型，設計對應的模型參數初始化和模型訓練方法，以及基于該模型的逆狀態識別任務求解方法.首先運用自適應粒子群優化（adaptive particle swarm optimization,APSO）算法獲取HMM 模型的初始參數λ0={A0,B0,π0}.然后通過多觀測序列的Baum-Welch 算法訓練得到最優參數λ*={A*,B*,π*}，使得HMM 模型在該給定參數 λ*的條件下，出現觀測干擾樣式序列O的概率P(O|λ*)最大.最后利用對數Viterbi 算法解決雷達方在已知模型最優參數 λ*和觀測干擾序列O的情況下，得到可解釋該觀測干擾序列O的干擾方識別結果的逆狀態識別序列

2.1 基于自適應粒子群算法的HMM 參數初始化

HMM 模型的初始參數λ0={A0,B0,π0}會直接影響模型的訓練結果.其中A0和 π0對模型訓練結果影響較小，可采用均勻分布或平均分布取值[28].對B0，本文采用自適應粒子群優化（APSO）算法進行參數初始化.

為解決帶有等式約束的HMM 初始參數優化問題，本文采用罰函數法[29]，增加懲罰項后適應度函數表示為

式中：函數f(x)即由前向-后向算法所得的P(O|λ) ，σ為懲罰因子；P(x)為整體懲罰項.

為計算P(x)將等式約束條件轉換為不等式約束，有

其中， ε為等式約束容忍度值.

整體懲罰項P(x)為

式中：Lk表示每個約束懲罰項的權重；ek(x)=max(0,|gk(x)|-ε).

2.2 基于多觀測序列HMM 參數優化的逆狀態識別

2.2.1 多觀測序列的HMM 參數優化

考慮截獲到多條獨立同分布的觀測干擾樣式序列，記為O′={O(1),O(2),···,O(K)}，其中K為觀測序列總數.本節基于多觀測序列Baum-Welch 算法和初始化后的HMM 參數λ0={A0,B0,π0}進行模型參數的優化訓練，具體如下：

記第k個觀測序列為O(k)=，其長度為nk.模型訓練的目標函數為

其中，將已知參數λ={A,B,π}時第k個干擾樣式觀測序列出現的概率P(O(k)|λ)記為Pk.將不同觀測序列頻率累加可得修正后的重估公式[30]為

2.2.2 基于對數Viterbi 算法的隱狀態估計

雷達方對干擾方輻射源狀態識別環節的逆向推理是狀態估計問題，因此在訓練觀測序列O′={O(1),O(2),···,O(K)}上進行參數估計獲得最佳模型 λ*后，針對測試觀測序列需要估計出對應的干擾方輻射源狀態識別結果序列.本文采用對數Viterbi 算法避免由于輸入的觀測序列較長，計算所得概率過小導致的數值下溢問題[31].

針對單觀測序列O={o1,o2,···,on}，定義變量χt(i)為式(12)，表示在已知HMM 參數λ={A,B,π}的條件下，到t時刻為止觀測序列為O={o1,o2,···,ot}情況下，出現概率最大的系統隱狀態序列，且=Si.

χt+1(i)可通過式(13)迭代得到.

3 仿真驗證

3.1 仿真場景設置

本文所用仿真場景如圖2 所示，雷達方為地基多功能相控陣雷達，干擾系統包括遠距支援和自衛干擾兩部分，可在對多功能雷達狀態進行識別后進行干擾決策，執行相應干擾樣式.

圖2 “空對地”干擾場景示意圖Fig.2 Diagram of "air-to-ground" jamming scenario

設置仿真場景中，干擾方有3 種干擾措施分別為密集假目標、噪聲壓制和拖引干擾；雷達方的4 種工作狀態為搜索、確認、跟蹤和識別；雷達方可采取的4 種抗干擾措施為無抗干擾措施、頻率捷變、旁瓣匿影和記憶波門跟蹤.對應地，逆狀態識別HMM模型中的干擾樣式個數M=3，隱狀態個數N=4，抗干擾措施集數L=4.表1 為專家經驗指導下的雷達雙方干擾與抗干擾措施設置.

表1 雷達對抗雙方的干擾與抗干擾措施設置Tab.1 Jamming and anti-jamming measures setting for both sides of radar countermeasure

3.2 仿真數據生成

本文依照第1 節中的交互模型基本思想，設計采用如圖3 所示流程產生仿真所需真實雷達狀態序列Q、干擾方對MFR 工作狀態識別結果序列、干擾方真實干擾樣式序列U及雷達方觀測干擾樣式序列O，用于模型訓練與測試.

圖3 仿真序列生成關系示意圖Fig.3 Diagram of simulation sequence generation

為了簡便，本文假設觀測干擾樣式序列O與干擾方真實干擾樣式序列U相同.

3.2.1 真實觀測特征概率分布矩陣生成方法

本節介紹真實觀測特征概率分布矩陣Br=[bri j]N×M的仿真方法，即brij=P(uk=Vj|qk=Di,zk=Wl).該矩陣的設置目的是為了仿真驗證本文的方法，該矩陣的不同設置實現方法對本文的逆狀態推理方法無影響，本文考慮如下仿真方法.采用如表2 所示的0～9 區間標度對不同干擾樣式的有效程度進行定量化干擾效果評估，標度值越大表示干擾樣式越有效，干擾方就越可能選取該干擾樣式.

表2 干擾有效性標度定義Tab.2 Jamming effectiveness scale definition

干擾方首先對不同雷達狀態、抗干擾措施下，各干擾樣式進行有效性標定，記雷達真實狀態為Di，抗干擾措施為Wl時，干擾樣式Vj的有效性標度為dVj|Wl,Di，選擇該干擾樣式Vj的概率為

為簡化問題，本文假設干擾方認為雷達在當前狀態下按照表1 等概率選取對應抗干擾措施，如當雷達處于搜索狀態時，雷達抗干擾措施為無措施、頻率捷變或旁瓣匿影的概率均為1/3.雷達狀態為Di時采取抗干擾措施Wl的概率為pWl|Di，由全概率公式可得Br=[bri j]N×M為

1）搜索狀態干擾樣式概率.

當雷達處于搜索狀態時，干擾有效性如表3 所示.

表3 雷達搜索狀態下干擾有效性表Tab.3 Jamming effectiveness under radar search mode

根據式(15)、(16)可算得當干擾方認為雷達在搜索狀態的條件下，選取噪聲壓制、密集假目標和波門拖引干擾的概率br11～br13分別為0.611、0.289 和0.

2）確認狀態干擾樣式概率.

雷達處于確認狀態時，使用大能量波束在目標的小區域內照射，可使目標回波信號強度增大.壓制類干擾效果比搜索狀態時所有降低，而密集假目標干擾效果受影響較小.綜上，確認狀態的干擾有效性如表4 所示.

表4 雷達確認狀態下干擾有效性表Tab.4 Jamming effectiveness in radar confirmation mode

計算可得干擾方狀態識別結果為確認狀態的條件下，干擾樣式選取噪聲壓制、密集假目標干擾和波門拖引干擾的概率br21～br23分別為0.546、0.454 和0.

3）跟蹤/識別狀態干擾樣式概率.

雷達跟蹤和識別狀態下的干擾有效性類似，均如表5 所示.

表5 雷達跟蹤/識別狀態下干擾有效性表Tab.5 Jamming effectiveness in radar tracking / recognition mode

計算可得干擾方狀態識別結果為跟蹤/識別狀態的條件下，干擾樣式選取噪聲壓制、密集假目標干擾和波門拖引干擾的概率分別為0.38、0.224 和0.396.

至此，用于生成干擾方真實干擾樣式序列U和雷達方觀測干擾樣式序列O的真實觀測特征概率分布矩陣如下：

3.2.2 真實雷達狀態轉移概率矩陣生成方法

本節介紹真實雷達狀態轉移概率矩陣Ar=[ari j]N×N的仿真方法.同3.2.1 節，該矩陣的仿真實現僅為一個實例，不影響本文的方法.根據雷達狀態工作原理，有如圖4 所示的雷達狀態轉換圖.與用于生成干擾樣式序列的真實觀測特征概率分布矩陣類似，雷達各工作狀態間的轉移概率受干擾方干擾和雷達方抗干擾措施間的相互作用影響.

圖4 雷達工作狀態轉移示意圖Fig.4 Diagram of radar work mode transition

真實雷達狀態轉移概率矩陣記為Ar=[ari j]N×N,ari j=P(qk+1=Dj|qk=Di)，則有

與干擾樣式序列生成中抗干擾方式設置相似，則上式可寫為

其中,P(qk+1=Dj|zk=Wl,ok=Vm,qk=Di)表示觀測干擾方干擾樣式和雷達方抗干擾措施條件下的雷達狀態轉移概率.下面以跟蹤/識別狀態下的狀態轉移概率為例進行求解.

跟蹤和識別狀態可轉移至搜索、跟蹤或識別3種狀態，與搜索和確認狀態相同，轉為搜索狀態的概率與干擾有效性標度呈正相關，如式(20)所示.當干擾有效性標度大于5 時，從跟蹤/識別狀態轉移至跟蹤狀態的概率應大于轉移至識別狀態的概率，當有效性標度小于5 時，則相反，具體如下:

其中， λ為可調參數，當干擾有效性標度dVj|Wl,Di＞5時，λ=0.6，否則λ=0.4.表6 所示為跟蹤/識別狀態轉移概率.

表6 雷達跟蹤/識別狀態轉移概率表Tab.6 Transition probability in radar tracking / recognition mode

根據雷達方在每種狀態下等概率選取表1 中設置的抗干擾措施，并將表6 所示的狀態轉移概率帶入式(19)可得真實雷達狀態轉移概率矩陣為

3.2.3 干擾方的雷達狀態識別結果序列生成方法

干擾方狀態識別結果序列Q?的生成是在真實雷達狀態序列Q的基礎上假定干擾方狀態識別算法存在一定概率識別錯誤.即若序列長度為nk，干擾方狀態識別準確率為 γ，則干擾方狀態估計結果序列和真實雷達狀態序列在同一時刻相同的時刻數目為nk·γ，其余時刻的雷達狀態被干擾方識別錯誤為其他狀態.

3.3 仿真驗證分析

1）HMM 參數初始化算法分析.

本節仿真通過對比逆狀態識別準確率 ζ，觀察不同參數初始化算法對模型估計準確率的影響，其中ζ為

以前向-后向算法輸出的P(O|λ)作為適應度函數，算法迭代次數上限為500，在初始化過程中PSO 和APSO 優化算法的適應度值變化情況如圖5 所示.

圖5 HMM 參數初始化適應度變化曲線圖Fig.5 Fitness change curve of HMM parameter initialization

由圖5 可知，在基于PSO 和APSO 算法的HMM參數初始化過程中，適應度值隨著算法迭代次數增加而增長并逐漸趨于平穩，其達到平穩的時間可表示算法收斂速度.故APSO 算法初始化方法的收斂速度快于PSO 算法，且APSO 算法的適應度值在收斂后更接近于最優值1.綜上所述，基于APSO 算法的HMM 參數初始化方法的全局尋優能力更好.隨后利用隨機參數初始化、PSO 和APSO 算法初始優化后的參數分別作為HMM 模型訓練的起始參數，按本文所提方法完成逆狀態識別.

訓練所需樣本數據根據3.2 節所設置矩陣生成.其中生成的干擾方狀態識別結果準確率均勻分布在區間70%～100%.仿真產生訓練數據序列為50 組，每組序列長度為50.測試樣本與訓練樣本的生成方法相同，序列長度設置為100.對同一測試樣本基于不同參數初始化方法進行逆狀態識別的準確率如表7 所示.

表7 基于不同參數初始化方法的逆狀態識別準確率Tab.7 Inverse mode recognition accuracy of different parameter initialization methods

由表7 可知，對于同一測試樣本，采取隨機參數、PSO 和APSO 算法初始化方法實現的逆狀態識別準確率依次提升.且從平均水平角度觀察，APSO 算法參數初始化方法也優于另外兩種算法.

2）HMM 模型參數優化先驗知識分析.

采用較優的APSO 算法進行參數初始化后，為分析HMM 模型參數優化訓練對先驗知識的需求，利用Baum-Welch 算法模型訓練開展HMM 參數優化實驗.訓練樣本集與1）中所用的樣本集相同.參數優化完成后，將1）中生成的測試樣本輸入到訓練完成的逆狀態識別HMM 模型中進行隱狀態估計，逆狀態識別準確率隨訓練樣本數變化曲線如圖6 所示，圖中準確率為對多組測試序列進行逆狀態識別后的平均準確率.

圖6 不同訓練樣本數下逆狀態識別準確率變化曲線圖Fig.6 Inverse mode recognition accuracy change curve under different training sample amounts

由圖6 可知，訓練樣本數在0～15 階段，隨著訓練樣本數的增加，逆狀態識別準確率迅速提升，當訓練樣本數到達15 時，逆狀態識別準確率可達90%以上.當訓練樣本數達到約25 時，準確率可達到約94%，隨著訓練樣本數進一步增加，逆狀態識別準確率趨于平穩.仿真結果表明在先驗知識較少的情況下，HMM 模型參數優化訓練能夠較準確地對干擾方狀態識別結果進行估計，且參與HMM 參數優化訓練的樣本數量越多，雷達方對干擾方多功能雷達工作狀態識別的逆向推理準確率越高.

3）HMM 模型逆狀態識別穩定性分析.

本驗證基于APSO 算法進行HMM 參數初始化，并利用30 組樣本訓練后的模型參數作為最終的逆狀態識別HMM 模型，將干擾方對MFR 狀態識別準確率為70%～100%下生成的測試樣本分別輸入HMM，可得如圖7 所示的實驗結果.

圖7 不同干擾方狀態識別準確率下逆狀態識別準確率變化曲線圖Fig.7 Inverse mode recognition accuracy change curve under different jammer mode recognition accuracy rates

由圖7 可知，在訓練樣本數為30 的條件下，當干擾方狀態識別準確率變化時，雷達方對干擾方狀態識別結果估計的準確率仍穩定在約94%，說明針對MFR 輻射源工作狀態識別的逆向推理方法不會受干擾方狀態識別準確率的影響，有較好的穩定性.

4 結 論

本文針對對抗場景中提升雷達系統抗干擾能力的實際需求，研究了一種基于隱馬爾可夫模型對干擾方的雷達狀態識別結果實現逆向推理的方法.雷達方可以利用推理得到的干擾系統內部信號處理節點信息，優化調整雷達工作狀態對應的發射波形方案，也可以基于所得HMM 參數，實現對干擾方后續干擾動作的預測，提升系統主動抗干擾的能力.設定雷達對抗場景中的逆狀態識別仿真實驗，驗證了本文任務驅動下模型的可行性及相應算法性能的穩定性.