反艦導彈輻射源行為分析中的貝葉斯方法*

崔 凱，宋新超，秦長海

（1.江蘇科技大學電子信息學院，江蘇 鎮江 212003；2.中國船舶重工集團第七二三研究所，江蘇 揚州 225101）

0 引言

由于科技的提升使得現在的反艦導彈精確打擊能力不斷加強，這也要求海軍艦隊必須升級對抗反艦導彈的防御系統。世界各國經過一段時間的研究，已經有部分防御方式應用到實際，比如密集陣［1］、反導導彈［2］、誘餌箔條［3］、電子干擾等。

當反艦導彈發起進攻時，其輻射源通常會經歷4 種行為［4］，如下：

1）搜索模式（Search Mode）。導彈發射出去一段時間后，導引頭雷達開機，用于對目標進行搜索。

2）分類識別模式（Classify Mode）。當搜索到有目標時，會對目標進行識別并確認，此過程通常只有10 s 不到。

3）跟蹤模式（Tracking Mode）。一旦制導系統確定了需要攻擊的目標，便會進入該模式跟蹤目標。

4）干擾尋的（HOJ，Home On Jamming）模式。如果導彈制導系統檢測到目標正在釋放電子干擾，其輻射源將會轉為該模式并進一步提高反艦導彈在電子干擾下對目標的精確打擊概率。

當艦隊正在受反艦導彈威脅時，需要采取有效的防御措施。為了能夠達到很好的防御效果，首先要對反艦導彈輻射源的行為以及行為變化進行判斷分析，本文的研究重點就是以貝葉斯模型中常用的蒙特卡洛馬爾可夫鏈算法（Markov Chain Monte Carlo，MCMC）對反艦導彈輻射源的行為以及行為變化進行有效分析預測。雖然反艦導彈輻射源有4 種行為模式，但在電子對抗研究領域，可以將搜索模式與分類識別模式歸為搜索狀態，將跟蹤模式和干擾尋的模式歸為跟蹤狀態，本文也重點研究反艦導彈輻射源的這兩種行為狀態。

1 戰斗中反艦導彈攻防體系分析

以下是反艦導彈發起進攻以及艦隊自我防御的示意圖（圖1）［5］。HVU（High Value Units）指海上高值目標，為反艦導彈重點攻擊對象，戰斗過程分為以下幾步驟：

1）導彈發射出去一段時間后，其輻射源開機并開始搜索海面目標。

2）在距離目標10 n mile，制導系統已經成功定位目標并且導彈輻射源轉為跟蹤模式，反艦導彈開始追蹤目標。

3）離目標還有8 n mile 左右時，HVU 上的設備開始釋放出電子干擾來迷惑反艦導彈的制導系統。

4）距離艦隊6 n mile 左右，此時的反艦導彈正受電子干擾影響，其輻射源可能已經轉為干擾尋的（HOJ）模式并指引導彈繼續追蹤HVU，也可能會追蹤電子干擾產生的虛擬目標。

5）如果導彈追蹤虛擬目標或者誘餌箔條，HVU就會停止釋放電子干擾。如果導彈仍然追蹤HVU，則密集陣待命。

圖1 反艦導彈進攻，艦隊防御示意圖

在對抗反艦導彈的防御過程中，最重要的一步就是釋放電子干擾。釋放干擾的時間不能過早，也不能過遲。如果過早釋放干擾脈沖，會使反艦導彈輻射源提前轉為干擾尋的（HOJ）模式，該模式下制導系統仍然會指引反艦導彈精確打擊目標，航母編隊照樣受到威脅。如果太遲釋放電子干擾，則為時已晚，艦隊照樣受到嚴重打擊。

由于反艦導彈輻射源的行為從搜索狀態轉為跟蹤狀態時會有不到10 s 的時間對目標進行分類識別，所以要在其行為發生改變的那一瞬間釋放電子干擾，才能對反艦導彈制導系統達到最佳干擾效果。因此，需要掌握的信息是，反艦導彈輻射源行為開始發生變化時，它與航母編隊的距離是多少。

2 基于貝葉斯推理的基本概念

假設兩個隨機事件為D 和H，貝葉斯定理的公式［6］如式（1）所示：

如果將H 表示假設，D 表示數據，那么貝葉斯理論的定義就是在引入觀測數據的前提下，假設成立的概率為多少。因此，在作貝葉斯分析時，以下幾個參數特別重要：先驗分布p（H）、似然分布p（D|H）和后驗分布p（H|D）。p（D）為證據，也作邊緣似然，但是對數據作貝葉斯分析時很少用到該項，大多數情況下也可以直接忽略。這樣貝葉斯公式就有另一種表達方式：

綜上所述，在給數據作貝葉斯分析前，先要根據樣本的狀態選擇合適的先驗分布和似然分布［7］。現在的反艦導彈通常以主動雷達制導模式為主，本文就以反艦導彈輻射源的脈沖描述字（PDW）和導彈與HVU 之間的距離為原始樣本，通過貝葉斯模型訓練出PDW 不同的后驗分布，根據不同的后驗均值來判斷出反艦導彈輻射源不同的行為狀態以及行為發生變化的位置。

3 基于MCMC 的反艦導彈輻射源行為變化的估計

3.1 MCMC 基本原理

MCMC 是所有貝葉斯推理算法中常用模型［8］，其基本思路是：選擇合適的抽樣算法從條件概率分布中進行采樣，結束后得到樣本的后驗概率分布，并將后驗結果作為正態化常數。

此方法分為馬爾可夫鏈和蒙特卡洛兩部分［9］，馬爾可夫鏈是表達數據狀態轉移的隨機過程，而蒙特卡洛方法是對數據的抽樣過程。MCMC 的原理［10］是將隨機過程馬爾可夫鏈引入蒙特卡洛抽樣中進行模擬，在模擬過程中對樣本進行抽樣訓練得到另一個馬爾可夫鏈。本文就將雷達及偵察設備捕獲到的反艦導彈輻射源的PDW 數據作為馬爾可夫鏈狀態樣本，當經過多次模擬后，再基于這些樣本進行隨機抽樣評估。

假設隨機過程{Xn，n∈N}的初始狀態為X0，自然數集N 表示隨機過程的時間參數，I 為隨機過程的狀態參數集，令狀態ai∈I 為隨機過程所有可能狀態。對于任意n∈N 以及a1，a2，a3，…，an+1∈I，條件概率滿足以下公式：

則稱{Xn，n∈N}為馬爾可夫鏈，且當n→∞時，Xn與初始X0狀態無關。在馬爾可夫鏈中，隨機變量的分布收斂于穩態分布π（x），稱為馬爾可夫鏈的收斂，那么在反艦導彈輻射源的PDW 值在后驗分布中，其概率密度函數f（x）關于π（x）的期望為：

3.2 反艦導彈輻射源行為變化估計的全過程

根據上述MCMC 的原理，通過MCMC 對反艦導彈輻射源行為及行為變化進行估計的步驟如下：

1）以反艦導彈輻射源的PDW 數據為訓練樣本X，以反艦導彈與HVU 的距離為行為發生變化參數，參數集I 為搜索狀態和跟蹤狀態兩種行為。

2）選取某一處的PDW 值記為X0，從X0開始，利用馬爾可夫鏈方法對飛行中的反艦導彈PDW 值產生點序列X1，X2，…，Xn。

3）馬爾可夫鏈經過多次迭代后，各處PDW 的邊際分布都是平穩分布，因此，最終馬爾可夫鏈收斂。收斂后的后驗結果就是反艦導彈輻射源PDW的期望值。

4）根據PDW 在不同行為狀態下的兩種不同的期望值，得出反艦導彈輻射源行為發生變化的位置。

3.3 Metropolis-Hastings 算法

蒙特卡洛方法是各類抽樣算法組成的，其中Metropolis-Hastings 算法是一種比較常見的抽樣法，Metropolis-Hastings 算法的步驟如下［11-12］：

4 實驗仿真與分析

4.1 實驗環境介紹

本文的實驗環境是基于python3.6 語言，用于貝葉斯分析的python 庫PyMC3 3.0 和用于可視化的Python 庫Matplotlib 2.0.2。

4.2 實驗過程分析

假設反艦導彈正在對HVU 發起進攻，同時艦隊上的雷達以及偵察設備檢測到導彈輻射源的PDW 值和位置（導彈與HVU 的距離）并記錄在數據庫，仿真實驗以載頻（RF）、脈沖幅度（PA）和脈沖重復周期（PRI）3 類PDW 參數為原始樣本。它們的關系分別如圖2～圖4 所示。

圖2 反艦導彈在不同位置的載頻

圖3 反艦導彈在不同位置的幅度

圖4 反艦導彈在不同位置的PRI

如果僅僅觀察3 張關系圖，或許不太容易直接觀察出反艦導彈在何處其輻射源行為發生變化，接下來就借助MCMC 方法進行分析，首先對載頻（RF）進行數據分析。

對于距離disRF，很難選擇它的先驗。在貝葉斯分析中，如果先驗分布比較難選擇，通常就用一定范圍內的均勻分布表示先驗。因此，dis 的先驗為均勻分布：

似然的選擇：從圖中可以看出數據是離散形式，而泊松分布能模擬計數類型的數據，設Di表示反艦導彈在i 處的載頻信號大小。因此，似然的選擇為：

選好合適的先驗和似然后就可以進行仿真了，本次實驗用Metropolis-Hastings 算法對數據進行采樣，采樣次數設為20 000 次。載頻分析完后，再繼續對PA 和PRI 進行分析，方法和步驟與上述相同。用PA 和PRI 作原始樣本時，分別假設反艦導彈與HVU 的距離為disPA 和disPRI。

4.3 實驗結果分析

由于要掌握反艦導彈輻射源行為發生變化時的具體位置，因此，只要觀察disRF，disPA 以及disPRI 的后驗結果。結果如圖5～圖7 所示。

圖5 disRF 的后驗分布

圖6 disPA 的后驗分布

在上述的3 組后驗分布里，每組上圖為核密度估計，是平滑后的頻率直方圖，最高點表示頻率最大。下圖是馬爾可夫采樣鏈，如果發現它在某一值附近振蕩，說明收斂狀態較好。初步觀察，disRF 和disPRI 的后驗分布沒有收斂。后驗分布收斂較好的是disPA，仔細觀察其后驗分布，它的核密度估計在18.5 左邊頻率最大，它的馬爾可夫采樣鏈在18 的上方某處振蕩，因此，可以初步判斷反艦導彈離HVU 為18 km～18.5 km 之間時，其輻射源行為發生變化。

圖7 disPRI 的后驗分布

disRF 和disPRI 無法收斂的原因是導彈輻射源行為發生變化時，RF 和PRI 的均值始終保持不變。當繼續觀察最初的PDW 數據時，可以發現RF 的變化范圍不大，但是PRI 的抖動范圍是有變化的，所以是否可以通過PRI 抖動的范圍不同來判斷導彈輻射源行為發生變化的位置，接下來就繼續通過實驗仿真來驗證，步驟如下：

1）首先計算出所有不同位置的PRI 與導彈飛行過程中所有PRI 均值的差的絕對值，結果如圖8所示：

圖8 各處PRI 與PRI 均值的差的絕對值

2）以所有PRI 值與它們均值差的絕對值為訓練樣本，步驟與4.2 節所述相同，這時假設反艦導彈與HVU 的距離為disPRI2，它的后驗如圖9所示：

disPRI2 的后驗結果中，核密度估計在18.5 右邊頻率最大，而馬爾可夫采樣鏈中，也在18～19 之間的某處振蕩，因此，可以判斷反艦導彈離HVU 的距離18.5 km 左右，其輻射源行為發生變化。

除了觀察后驗分布的核密度估計和馬爾可夫采樣鏈外，在PyMC3 中還能使用summary 函數得到具體的后驗分布均值（mean）、標準差（sd）和蒙特卡洛誤差（mc_error），這些值也可以用來分析后驗結果。結果如表1 所示：

表1 disRF，disPA 和disPRI 的后驗結果

在MCMC 中，希望sd 和mc_error 越小越好。如果mc_error 小于sd 的5 %，則可以說明收斂精度高。表中disRF 和disPRI 的這兩個值都偏大，disPA和disPRI2 的這兩個值比較小，且也滿足收斂精度要求。根據表中disPA 和disPRI2 的后驗mean 結果，仍然可以判斷出當反艦導彈距離航母編隊18.5 km 左右時，其輻射源行為發生改變。

4.4 實驗誤差分析及總結

上述仿真中，用于訓練的PDW 的每類參數樣本都很多。如果訓練樣本變少，后驗結果又如何？圖8 中，現在每0.4 km 取一次值，使得原始樣本只剩下70 個左右，設disPRI3 表示此時的反艦導彈與HVU的距離。樣本原始分布圖和仿真后的后驗分布分別如圖10 和圖11 所示：

圖10 減少樣本后的各處PRI 與PRI 均值的差的絕對值

圖11 disPRI3 的后驗分布

disPRI3 的后驗分布中，核密度估計有3 個峰值，馬爾可夫采樣鏈的振蕩范圍也偏大，很顯然與disPRI2 的后驗相比disPRI3 的后驗分布要模糊得多。

通過上述的實驗結果，可以得出以下結論：

1）由于仿真實驗中的反艦導彈輻射源的行為發生變化時，載頻的變化范圍不大，因此，很難通過單一頻域來判斷其輻射源的行為變化。

2）其輻射源的行為發生變化時，幅度的變化是很明顯的，而且幅度的均值也會變化，因此，通過幅度就可以分析出結果。

3）仿真實驗中的反艦導彈輻射源行為發生變化時，PRI 的均值始終不變，因此，僅僅通過PRI 也無法推斷反艦導彈輻射源行為變化，但是根據PRI的抖動范圍不同繼續分析就得到預想的后驗結果。因此，必要時要對原始樣本進行適當的預處理來簡化計算過程。

4）通常反艦導彈輻射源行為發生變化時，輻射源的時域、頻域或能量域會發生變化，因此，在分析其輻射源行為及行為變化時，要從輻射源的時域、頻域和能量域等多種角度綜合分析。

5 結論

本文以MCMC 方法，通過對反艦導彈輻射源的PDW 中的3 類參數進行研究處理，并判斷其行為以及行為在何處開始變化并完成仿真，仿真過程和結果表明，用該方法判斷反艦導彈輻射源的行為以及行為變化理論上是可行的，也充分證明貝葉斯模型在推理應用中的高效性。

但是貝葉斯模型的缺點也很明顯，可以歸納為以下兩點：1）用貝葉斯作推理時，要根據訓練樣本的狀態選擇合適的先驗分布，有時先驗分布不是特別容易判斷，如果先驗分布選擇不當就會導致后驗結果無法收斂；2）另外用于模型訓練的原始樣本也要足夠多，如果原始樣本太少，后驗的收斂效果也不會太好。