崔 凱，宋新超，秦長海

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

0 引言

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

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

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

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

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

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

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

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

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

1）導(dǎo)彈發(fā)射出去一段時(shí)間后，其輻射源開機(jī)并開始搜索海面目標(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 左右時(shí)，HVU 上的設(shè)備開始釋放出電子干擾來迷惑反艦導(dǎo)彈的制導(dǎo)系統(tǒng)。

4）距離艦隊(duì)6 n mile 左右，此時(shí)的反艦導(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)彈進(jìn)攻，艦隊(duì)防御示意圖

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

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

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

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

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

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

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

3.1 MCMC 基本原理

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

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

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

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

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

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

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)生點(diǎn)序列X1，X2，…，Xn。

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

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

3.3 Metropolis-Hastings 算法

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

4 實(shí)驗(yàn)仿真與分析

4.1 實(shí)驗(yàn)環(huán)境介紹

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

4.2 實(shí)驗(yàn)過程分析

假設(shè)反艦導(dǎo)彈正在對HVU 發(fā)起進(jìn)攻，同時(shí)艦隊(duì)上的雷達(dá)以及偵察設(shè)備檢測到導(dǎo)彈輻射源的PDW 值和位置（導(dǎo)彈與HVU 的距離）并記錄在數(shù)據(jù)庫，仿真實(shí)驗(yàn)以載頻（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 方法進(jìn)行分析，首先對載頻（RF）進(jìn)行數(shù)據(jù)分析。

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

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

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

4.3 實(shí)驗(yàn)結(jié)果分析

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

圖5 disRF 的后驗(yàn)分布

圖6 disPA 的后驗(yàn)分布

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

圖7 disPRI 的后驗(yàn)分布

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

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

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

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

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

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

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

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

4.4 實(shí)驗(yàn)誤差分析及總結(jié)

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

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

圖11 disPRI3 的后驗(yàn)分布

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

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

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

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

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

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

5 結(jié)論

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

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