反艦導(dǎo)彈輻射源行為分析中的貝葉斯方法*

崔 凱，宋新超，秦長海

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

0 引言

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

當(dāng)反艦導(dǎo)彈發(fā)起進攻時，其輻射源通常會經(jīng)歷4 種行為［4］，如下：

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

2）分類識別模式（Classify Mode）。當(dāng)搜索到有目標(biāo)時，會對目標(biāo)進行識別并確認(rèn)，此過程通常只有10 s 不到。

3）跟蹤模式（Tracking Mode）。一旦制導(dǎo)系統(tǒng)確定了需要攻擊的目標(biāo)，便會進入該模式跟蹤目標(biāo)。

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

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

1 戰(zhàn)斗中反艦導(dǎo)彈攻防體系分析

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

1）導(dǎo)彈發(fā)射出去一段時間后，其輻射源開機并開始搜索海面目標(biāo)。

2）在距離目標(biāo)10 n mile，制導(dǎo)系統(tǒng)已經(jīng)成功定位目標(biāo)并且導(dǎo)彈輻射源轉(zhuǎn)為跟蹤模式，反艦導(dǎo)彈開始追蹤目標(biāo)。

3）離目標(biāo)還有8 n mile 左右時，HVU 上的設(shè)備開始釋放出電子干擾來迷惑反艦導(dǎo)彈的制導(dǎo)系統(tǒng)。

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

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

圖1 反艦導(dǎo)彈進攻，艦隊防御示意圖

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

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

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

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

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

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

3 基于MCMC 的反艦導(dǎo)彈輻射源行為變化的估計

3.1 MCMC 基本原理

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

此方法分為馬爾可夫鏈和蒙特卡洛兩部分［9］，馬爾可夫鏈?zhǔn)潜磉_數(shù)據(jù)狀態(tài)轉(zhuǎn)移的隨機過程，而蒙特卡洛方法是對數(shù)據(jù)的抽樣過程。MCMC 的原理［10］是將隨機過程馬爾可夫鏈引入蒙特卡洛抽樣中進行模擬，在模擬過程中對樣本進行抽樣訓(xùn)練得到另一個馬爾可夫鏈。本文就將雷達及偵察設(shè)備捕獲到的反艦導(dǎo)彈輻射源的PDW 數(shù)據(jù)作為馬爾可夫鏈狀態(tài)樣本，當(dāng)經(jīng)過多次模擬后，再基于這些樣本進行隨機抽樣評估。

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

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

3.2 反艦導(dǎo)彈輻射源行為變化估計的全過程

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

1）以反艦導(dǎo)彈輻射源的PDW 數(shù)據(jù)為訓(xùn)練樣本X，以反艦導(dǎo)彈與HVU 的距離為行為發(fā)生變化參數(shù)，參數(shù)集I 為搜索狀態(tài)和跟蹤狀態(tài)兩種行為。

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

3）馬爾可夫鏈經(jīng)過多次迭代后，各處PDW 的邊際分布都是平穩(wěn)分布，因此，最終馬爾可夫鏈?zhǔn)諗?。收斂后的后驗結(jié)果就是反艦導(dǎo)彈輻射源PDW的期望值。

4）根據(jù)PDW 在不同行為狀態(tài)下的兩種不同的期望值，得出反艦導(dǎo)彈輻射源行為發(fā)生變化的位置。

3.3 Metropolis-Hastings 算法

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

4 實驗仿真與分析

4.1 實驗環(huán)境介紹

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

4.2 實驗過程分析

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

圖2 反艦導(dǎo)彈在不同位置的載頻

圖3 反艦導(dǎo)彈在不同位置的幅度

圖4 反艦導(dǎo)彈在不同位置的PRI

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

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

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

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

4.3 實驗結(jié)果分析

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

圖5 disRF 的后驗分布

圖6 disPA 的后驗分布

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

圖7 disPRI 的后驗分布

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

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

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

2）以所有PRI 值與它們均值差的絕對值為訓(xùn)練樣本，步驟與4.2 節(jié)所述相同，這時假設(shè)反艦導(dǎo)彈與HVU 的距離為disPRI2，它的后驗如圖9所示：

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

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

表1 disRF，disPA 和disPRI 的后驗結(jié)果

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

4.4 實驗誤差分析及總結(jié)

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

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

圖11 disPRI3 的后驗分布

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

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

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

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

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

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

5 結(jié)論

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

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