基于蒙特卡洛法的反艦導彈末端機動突防效果研究

范作娥，顧文錦，姜 鵬，馬培蓓

（1.海軍航空工程學院 控制工程系，山東 煙臺 264001；2.91065 部隊，遼寧 葫蘆島 125001）

0 引言

隨著艦空導彈武器系統(tǒng)的不斷發(fā)展，反艦導彈的突防問題變得越來越重要。為了增加對反艦導彈的攔截難度，提高反艦導彈的突防能力，反艦導彈需要進行末端機動。關(guān)于反艦導彈末端機動的控制與導引方法研究，以及各種機動方式的研究，國內(nèi)外已經(jīng)有相關(guān)的公開發(fā)表文獻[1-5]。相應(yīng)地，針對反艦導彈各種機動方式的突防效果的評估，也越來越受到學者們的關(guān)注。文獻[4]以突防概率作為評價標準，研究了突防導彈作正弦擺動突防機動時的突防效果。文獻[5]利用共軛系統(tǒng)和拉氏變換技術(shù)，研究了反艦導彈的螺旋機動和擺式機動的突防效果。這些文獻大多是應(yīng)用伴隨技術(shù)，針對兩維平面內(nèi)導彈與目標之間的相對運動方程，通過求解攔截導彈的脫靶量解析解，研究機動突防持續(xù)時間對穩(wěn)態(tài)脫靶量的影響，從而進一步研究對導彈突防概率有多大貢獻。但是，用伴隨技術(shù)來研究反艦導彈末端機動對突防概率有多大影響的時候，考慮到的實際對抗條件比較少，而且不能比較不同的末端機動方式之間，哪種方式最好。

本文就是從工程應(yīng)用的角度出發(fā)，在考慮盡可能多的實際影響因素的條件下，針對反艦導彈在有艦空導彈攔截的情況下，進行全彈道打靶仿真，采用蒙特卡洛法計算反艦導彈在不同機動方式下（含蛇行機動、擺式機動、螺旋機動）的突防概率，從而比較各種機動方式的優(yōu)劣。當反艦導彈進行末端機動之前，其控制系統(tǒng)的設(shè)計采用的是過載控制方法[6-7]，在末端機動開始之后，采用的是過載控制與位置控制結(jié)合的復合控制方法[8]。

1 反艦導彈末端機動一體化控制

反艦導彈末端機動一體化控制的原理主要是，在導彈進行末端機動的時候，同時引入過載控制與位置控制，分別將它們作為宏觀與微觀的控制手段。在宏觀上，過載控制預先給定了運動的趨勢；在微觀上，由位置反饋形成閉合系統(tǒng)，消除位置上細小的偏差。這樣的復合控制手段在仿真中取得了單獨過載控制與單獨位置控制均無法達到的控制效果，并且魯棒性比較強[8]。

1.1 末端機動期望軌跡描述

末端機動軌跡是導彈質(zhì)心相當于地面坐標系OXYZ的運動軌跡，通常是根據(jù)飛行的距離x 進行控制的，即以變量x為參數(shù)，此時，描述末端機動軌跡的位置參數(shù)方程形式為：

根據(jù)制導理論，在攔截前的恰當時間如果目標實施最大加速度機動，就能夠產(chǎn)生足夠大的脫靶量。通常，目標是不知道攔截導彈系統(tǒng)參數(shù)的，也不知道攔截時間，那么周期性的機動加速度是目標的一種最佳生存策略。本文用到的非平面機動方式（如擺式機動、螺旋機動）同樣是反艦導彈周期性地改變飛行軌跡，機動軌跡在縱向平面和航向平面內(nèi)的投影可分別用如下形式的正弦曲線和余弦曲線來近似表示：

式中：(x1y1z1)為末端機動的初始位置；ξ0為初始相位角，這里取 ξ0=? π /2；ly、lz分別為縱向和航向投影曲線的機動半徑；kyπ、kzπ分別為一個機動周期內(nèi)縱向和航向投影曲線的相位變化值；A稱為節(jié)距，它是指在一個機動周期內(nèi)反艦導彈沿地面坐標系OX軸方向飛行的距離。

定義節(jié)點數(shù)n，它是在末端機動范圍內(nèi)反艦導彈所完成的機動周期的數(shù)值。

考慮到末端機動是根據(jù)x 進行控制的，假設(shè)反艦導彈在飛行到 x=x2時停止機動，則有：n=(x2? x1)A。式中，ly、lz、ky、kz、A為末端機動參數(shù)，根據(jù)具體的末端機動方式進行選取。若選取 ly=0，則表示航向蛇行機動軌跡，z 隨x 按余弦曲線周期性地變化。由此可見，平面機動方式是非平面機動方式的特例。

1.2 末端機動一體化控制模型

由于末端機動中的控制信號包括法向過載和質(zhì)心兩組控制信號。因此，末端機動的一體化控制模型也應(yīng)該包括這兩組控制信號。即內(nèi)環(huán)的法向過載控制信號 n?y和 n?z，以及外環(huán)的質(zhì)心控制信號 y*和z*。具體的控制信號如下[7]：

其中，

公式(3)～(5)中，ny、θ 及vψ是導彈的飛行力學參數(shù)，具體見文獻[9]。上述兩種控制信號相互配合、協(xié)調(diào)控制，即可控制反艦導彈完成預定的各種末端機動彈道。末端機動的起始位置及機動范圍視具體情況定。

2 反艦導彈與艦空導彈對抗模型

2.1 艦空導彈截擊反艦導彈航向平面模型

如圖1所示，OXYZ為地面坐標系，反艦導彈M點的位置為(x,y,z)，攔截彈I點的位置為(xI,yI,zI)，實際的彈（攔截彈I）目（反艦導彈M）之間的距離大小為：

圖1 航向平面攔截關(guān)系

在航向平面內(nèi)，各角度的正方向的規(guī)定一樣，見圖1。圖中，Rxz為彈目距離，qxz為視線角，v 和vI分別為反艦導彈與艦空導彈的速度，ψM和ψI分別為兩者的彈道偏角，ηM和 ηI分別為兩者的前置角，即各自的速度矢量方向轉(zhuǎn)到視線方向所要旋轉(zhuǎn)的角度。

由圖1建立航向平面內(nèi)的攔截幾何關(guān)系：

式(7)、(8)中：

當 x ＞ xI，即攔截導彈尾追反艦導彈情況時，

2.2 艦空導彈截擊反艦導彈縱向平面模型

圖2中 zI軸的指向與OZ軸指向相反是因為艦空導彈彈道偏角| ψI| ＞ π/2。θM和 θI分別為反艦導彈和艦空導彈的彈道傾角，順時針轉(zhuǎn)動為正，其他都是逆時針轉(zhuǎn)動為正，故 θI＜0、視線角 qxy＜0，并且這兩個角度的正方向相反。前置角 ηI＞ 0，ηM＜0。

圖2 縱向平面攔截關(guān)系

縱向平面內(nèi)的攔截幾何關(guān)系：

式(12)、(13)中：

當 x ＞ xI，即攔截導彈尾追反艦導彈情況時，

3 艦空導彈的制導系統(tǒng)

鑒于艦空導彈的相關(guān)信息我們無法詳細知道，故在艦空導彈攔截反艦導彈的過程中，將艦空導彈視為質(zhì)點，只進行質(zhì)心運動，其運動方程[9]：

考慮艦空導彈的制導控制系統(tǒng)動態(tài)過程，將其簡化為3個1 階慣性延遲環(huán)節(jié)，其簡化動態(tài)過程如圖3所示，其中 T1=0.01 s為導引頭的時間常數(shù)；T2=0.02 s為控制回路的時間常數(shù)；T3=0.5 s為艦空導彈彈體的時間常數(shù)。

艦空導彈（攔截彈I）采用如下的比例導引律：

式中，N1與N2為艦空導彈的比例導引系數(shù)。

圖3 制導控制系統(tǒng)簡化圖

4 艦空導彈與反艦導彈的對抗仿真

在進行反艦導彈與艦空導彈的突防對抗仿真過程中，由于反艦導彈采用的是全彈道仿真模型，艦空導彈采用的是簡化模型，所以，為了使仿真結(jié)果盡可能的接近事實，反艦導彈考慮了多種外界及內(nèi)部干擾的影響，艦空導彈也考慮了一些本身的結(jié)構(gòu)約束，具體的敘述見下文。

4.1 反艦導彈所受干擾

1）海浪影響

反艦導彈飛行過程中，通過無線電高度表來測量導彈距海面的高度，送入慣導回路，計算得到組合高度，與裝訂高度進行比較，再根據(jù)高度控制方程進行高度控制。由此可見當海面波動時，無線電高度表的測量高度會因測量面的變化而變化，從而對導彈的高度控制產(chǎn)生影響，特別是在較大風浪時，導彈有可能因觸浪而導致反艦導彈掉海，因而必須對海浪進行模擬。在工程上，一般把海浪看作一種具有各態(tài)歷經(jīng)性的平穩(wěn)正態(tài)隨機過程，利用頻譜方法來進行研究。本文采用如式(20)所示的一維PM海浪譜作為海浪模型[10]：

式中：HS為海浪有效波高，單位m（由海情定，本文采用四級海情 HS=2.5 m）；ω為海浪角頻率。

根據(jù)式(20)，采用波浪疊加法對海浪進行仿真，即可得出海浪的瞬時波動 ξ (t)，則導彈相對海平面的瞬時飛行高度為

式中：yw(t)為考慮海浪波動的情況下導彈的飛行高度；y (t)為不考慮海浪波動的情況下導彈的飛行高度。

將導彈相對海平面的瞬時飛行高度代入導彈的導航模塊，從而可模擬海浪對導彈高度控制的影響。

2）建模不確定性影響

由于反艦導彈在飛行過程中，其飛行運動參數(shù)受干擾或建模不確定性等影響，其運動參數(shù)需要進行一定的攝動，故對其采用確定型隨機處理，使運動參數(shù)相對標稱運動參數(shù)進行最大幅度±30%的攝動。

3）高度表測量誤差

反艦導彈的控制系統(tǒng)的組成除了慣性導航系統(tǒng)外，還有必要的彈上測量儀表，例如高度表。反艦導彈在飛行過程中，通過雷達高度表來測量導彈距海面的高度，以此與裝訂高度進行比較。因此，高度表的測量誤差應(yīng)該在仿真對抗中被考慮進去，對測量高度值采用確定型隨機處理，進行±1%的攝動。

4）導彈結(jié)構(gòu)裝備誤差

由于導彈本身的控制系統(tǒng)和導航系統(tǒng)就是一個非常復雜的綜合系統(tǒng)，由其本身的結(jié)構(gòu)裝備引起的干擾和誤差，對其力和力矩肯定產(chǎn)生一定的影響，從而對導彈的控制系統(tǒng)也會產(chǎn)生一定的影響。因此，在對抗仿真過程中，對干擾力和力矩也采用確定型隨機處理，進行±1%的攝動。

5）指揮儀初始方位誤差

由于指揮儀的初始方位也是有儀器進行測量得來的，其初始方位角毫無疑問也會帶有誤差，故對其也采用確定型隨機處理，對其測量值進行±1%的攝動。

6）陀螺漂移影響

反艦導彈的控制系統(tǒng)中，因陀螺儀在測量角度的過程中帶有漂移，故陀螺漂移所引起的ψ值有所變化，在對抗仿真中考慮陀螺漂移速度1 n mile/h。

4.2 艦空導彈所受約束

1）導引頭前置角約束

設(shè)定艦空導彈的導引頭有前置角約束范圍：?5 0° ～ + 50°。當前置角達到此限制范圍之外時，導引頭失效，艦空導彈沿失效點的方向作直線運動。

2）導引頭盲區(qū)約束

設(shè)定艦空導彈的導引頭盲區(qū)為500 m，即導彈進入死區(qū)后，導引頭失效，艦空導彈沿失效點的方向作直線運動。

4.3 仿真條件

當反艦導彈在地面坐標系內(nèi)沿OX軸方向飛行到Xstart的時候，開始以縱向ly和航向lz投影曲線機動半徑分別進行不同方式的末端機動（航向蛇行、縱向蛇行、擺式機動和螺旋機動），與此同時，艦空導彈由初始位置（xI，yI，zI）開始，以相同的初始飛行速度vI進行攔截。反艦導彈開始進行變軌機動時的初始相位角為0ξ，且在一個機動周期內(nèi)，反艦導彈沿OX軸方向飛行的距離為A。

值得注意的是，反艦導彈進行不同方式的變軌機動時，各種變軌方式在航向和縱向投影曲線的最大機動半徑是相同的，即ly=lz=200 m。只有進行擺式機動的時候例外，其縱向投影半徑為ly=50 m。設(shè)Xstart=240 km，vI=850 m/s，ξ0=?π/2，A=7 000 m。

4.4 仿真結(jié)果與分析

在上述的仿真條件下，進行艦空導彈攔截反艦導彈的對抗仿真，當艦空導彈的攔截初始位置不同時，給出如下兩種仿真結(jié)果：

1）當艦空導彈初始攔截位置（xI，yI，zI）為xI=265 km、yI=80 m、zI=1 000 m時，進行不同機動方式的攔截對抗仿真，脫靶量數(shù)據(jù)結(jié)果見表1。

表1 攔截導彈對不同機動方式反艦導彈的脫靶量之一

2）當艦空導彈的初始攔截位置（xI，yI，zI）變?yōu)閤I=244 km、yI=80 m、zI=19.66 km時，進行不同機動方式的攔截對抗仿真，脫靶量數(shù)據(jù)結(jié)果見表2。

表2 攔截導彈對不同機動方式反艦導彈的脫靶量之二

由表1、2 中的仿真數(shù)據(jù)看出，應(yīng)用過載控制方法設(shè)計的反艦導彈，在有艦空導彈攔截的情況下，進行不同的末端機動方式進行突防的時候，攔截導彈的脫靶量要比不進行末端機動的時候高。同樣，反艦導彈進行帶有縱向機動的變軌方式，攔截導彈的脫靶量要明顯高于只進行航向時的機動方式。

但是由于縱向蛇行、擺式機動和螺旋機動這三種方式，其脫靶量差距不是很大，而且沒有明顯的規(guī)律可尋，所以要想比較那種方式比較好，僅僅通過相同條件的一兩次仿真，不足以說明問題，這將是下一節(jié)要解決的問題。

5 基于蒙特卡洛法的突防概率研究

本節(jié)主要是在前面模型和仿真基礎(chǔ)上，在反艦導彈進行不同形式的末端機動時候，艦空導彈由不同的初始位置進行攔截，在計算機上進行多次打靶仿真，采用蒙特卡洛法，針對打靶數(shù)據(jù)計算各種末端機動方式的突防概率，比較末端直線運動、航向蛇行機動、縱向蛇行機動、螺旋機動和擺式機動，那種機動方式的突防概率最高。

5.1 仿真條件

利用蒙特卡洛法來計算反艦導彈的突防概率，比較各種末端機動方式的優(yōu)劣，需要不同的初始仿真條件下，在計算機上進行多次打靶仿真，具體的初始條件的選取，如圖4所示。圖中所示坐標系為地面坐標系中的XOZ 平面，M為反艦導彈的位置，其坐標為（xstart，0），I為攔截彈的位置，其坐標為(xI,zI)，反艦導彈與艦空導彈之間的距離為R1，艦空導彈的初始方位角為ψL。

圖4 初始攔截位置圖

當反艦導彈沿OX軸方向飛行到xstart時，開始進行不同方式的變軌機動，與此同時，艦空導彈由初始攔截位置 (xI,zI)開始，以飛行速度vI對反艦導彈進行攔截，仿真中用到的反艦導彈的相關(guān)參數(shù)同第4節(jié)相同，不同的只是艦空導彈的初始攔截位置及初始方位角不同。關(guān)于艦空導彈初始攔截條件的選擇，遵循如下兩種選擇方法：

1）針對每一種末端機動方式，在艦空導彈的初始方位角ψL相等的情況下，按照不同的反艦導彈與艦空導彈之間的距離 R1，來定義艦空導彈的初始攔截位置。其中，R1分別選擇為10 km、20 km 和30 km，這樣就定義了3個初始攔截條件。

2）針對每一種末端機動方式，在反艦導彈與艦空導彈之間的距離 R1相等的情況下，按照不同的艦空導彈的初始方位角ψL，來定義艦空導彈的初始攔截位置。其中，ψL在(0 ～360°) 范圍內(nèi)，由0°開始，每隔5°選擇一次，定義為一個初始攔截條件，這樣就定義了72個初始條件。

聯(lián)合上述兩種初始攔截條件的選擇方法，可以得出每一種末端機動方式共有216種不同的初始仿真條件，即每一種末端機動方式都要進行216次打靶仿真。由于每種末端機動方式都需要進行216次打靶，所以，要想比較各種末端機動方式的優(yōu)劣（直線飛行、航向蛇行、縱向蛇行、擺式、螺旋機動），就需要總共進行1 080次打靶仿真。

5.2 仿真結(jié)果與分析

在上述初始仿真條件下，進行攔截對抗仿真，總共得到艦空導彈1 080個脫靶量。以殺傷半徑5 m來定義命中區(qū)，即艦空導彈的脫靶量只要小于5 m，就當作艦空導彈成功攔截反艦導彈。針對這些記錄的不同機動方式的打靶數(shù)據(jù)，計算整理反艦導彈的突防概率如表3所示。

表3 不同機動方式的突防概率

對比仿真數(shù)據(jù)表1～3，看出表2、3的仿真結(jié)果完全對應(yīng)，但表1中，縱向蛇行機動的脫靶量為18.455 m，而螺旋機動的脫靶量為15.978 m，這與表3中的仿真結(jié)果不相符。造成這種結(jié)果的主要原因是，反艦導彈在進行縱向蛇行機動時，如果反艦導彈飛行在正弦曲線的峰值點，這時橫向過載最大，就會造成大脫靶量，但它只是蛇行飛行狀態(tài)下的一種個例，其他情況都是以小橫向過載來飛行；而進行螺旋機動的時候，反艦導彈時時刻刻都是以最大過載飛行，因此從統(tǒng)計數(shù)據(jù)中來看，縱向蛇行機動的突防概率仍然小于螺旋機動的突防概率。

由表3還可以看出，同樣是進行平面機動，縱向蛇行機動的突防概率為71.76%，要明顯高于航向蛇行機動的突防概率為41.2%。造成這種結(jié)果的原因是，在對抗仿真的過程中，艦空導彈在攔截進行縱向蛇行機動的反艦導彈的時候，容易掉海，這就造成了縱向機動的突防概率明顯高于航向機動。

由仿真結(jié)果和分析得出，在有艦空導彈攔截的情況下，關(guān)于反艦導彈突防概率有3個結(jié)論：① 反艦導彈進行末端機動的突防概率，要比末端直線飛行的時候大；② 反艦導彈進行縱向平面行機動的突防概率，要比航向平面機動的時候大；③ 反艦導彈以不同的末端機動方式突防的時候，按照突防概率的大小，機動方式由優(yōu)至劣的順序是：螺旋機動—擺式機動—縱向蛇行機動—航向蛇行機動—直線飛行。

6 結(jié)束語

本文從工程應(yīng)用的角度出發(fā)，在考慮盡可能多的實際影響因素及約束條件的影響下，建立了反艦導彈在有艦空導彈攔截情況下的對抗模型，針對上述模型，在艦空導彈不同攔截初始條件下，進行了上千次全彈道打靶仿真，應(yīng)用蒙特卡洛法，以打靶試驗數(shù)據(jù)為依據(jù)，比較直觀全面的比較了各種末端機動方式的突防效果，這在工程上具有較大實用參考價值。但是，本文在建模過程中，艦空導彈采用的是簡化模型，因此，如何將艦空導彈模型擴展為全彈道仿真模型，并在此基礎(chǔ)上考慮艦空彈的實際影響因素，使得仿真數(shù)據(jù)更接近工程應(yīng)用，將是進一步研究的課題。

[1]IMADO F,MIWA S.Missile guidance algorithm against high-g barrel roll maneuvers[J].Journal of Guidance,Control,and Dynamics,1994,17(1):123-128.

[2]顧文錦,畢蘭金,武志東.基于伴隨技術(shù)的超音速反艦導彈末端機動突防的脫靶量分析[J].海軍航空工程學院學報,2009,24(4):400-404.

[3]周狄,鄒昕光,孫德波.導彈機動突防滑模制導律[J].宇航學報,2006,27(2):213-216.

[4]姜玉憲,崔靜.導彈擺動式突防策略的有效性[J].北京航空航天大學學報,2002,28(2):133-136.

[5]顧文錦,趙紅超,王鳳蓮.反艦導彈末端機動的突防效果研究[J].宇航學報,2005,26(6):758-761.

[6]趙紅超,顧文錦,馬登武.飛航導彈基于輸出重定義的變結(jié)構(gòu)控制[J].兵工學報,2005,26(6):783-786.

[7]趙紅超,顧文錦,王瑞奇.反艦導彈的自適應(yīng)全局滑模變結(jié)構(gòu)控制[J].控制工程,2005,12(4):320-322.

[8]顧文錦,趙紅超,王鳳蓮.導彈末端機動的一體化控制模型[J].宇航學報,2004,25(6):677-680.

[9]錢杏芳,林瑞雄,趙亞男.導彈飛行力學[M].北京:北京理工大學出版社,2000:53-57.

[10]賴志昌,金鴻章,李國斌,等.隨機海浪作用下的船舶橫搖減搖預報方法[J].哈爾濱工程大學學報,2001,22(3):13-16.