基于Matlab Simulink的彈道仿真方法

黃　成，于　鵬，趙錫旺（. 9550部隊(duì)，遼寧大連60；.海軍裝備部，北京0008；. 970部隊(duì)，山東威海6409）

?

基于Matlab Simulink的彈道仿真方法

黃成1，于鵬2，趙錫旺3
（1. 91550部隊(duì)，遼寧大連116023；2.海軍裝備部，北京100083；2. 92270部隊(duì)，山東威海264309）

摘要：彈道仿真軟件有助于飛航導(dǎo)彈武器飛行試驗(yàn)結(jié)果的分析與評(píng)定。為此，設(shè)計(jì)了基于Matlab Simulink 的彈道仿真計(jì)算方法，建立了六自由度空間彈道數(shù)學(xué)模型，利用導(dǎo)彈相關(guān)的氣動(dòng)數(shù)據(jù)和控制參數(shù)測(cè)試了算法，實(shí)現(xiàn)了軟件，仿真結(jié)果證明了該模型的準(zhǔn)確性和可信度。

關(guān)鍵詞：彈道仿真；數(shù)學(xué)模型；仿真方法

飛行試驗(yàn)是導(dǎo)彈研制定型過(guò)程中的一個(gè)關(guān)鍵環(huán)節(jié)。從以往的導(dǎo)彈飛行試驗(yàn)來(lái)看，當(dāng)導(dǎo)彈在飛行試驗(yàn)中出現(xiàn)故障或問(wèn)題時(shí)，試驗(yàn)現(xiàn)場(chǎng)常因欠缺彈道仿真軟件而陷入被動(dòng)，這在一定程度上影響了試驗(yàn)進(jìn)展。近年來(lái)，國(guó)內(nèi)關(guān)于導(dǎo)彈的仿真，大多集中于仿真訓(xùn)練和仿真視景的研究[1-5]，而對(duì)于輔助導(dǎo)彈飛行試驗(yàn)結(jié)果分析、評(píng)定的仿真研究相對(duì)較少。為加強(qiáng)對(duì)飛航導(dǎo)彈武器系統(tǒng)飛行試驗(yàn)結(jié)果的分析與評(píng)定，亟需深入開(kāi)展便利、有效的相關(guān)彈道仿真方法研究。

對(duì)飛航導(dǎo)彈的全彈道仿真來(lái)講，仿真研究的全過(guò)程涉及到多個(gè)功能模塊的設(shè)計(jì)難題，如參數(shù)裝定模塊、運(yùn)算模塊、結(jié)果處理模塊、運(yùn)行控制模塊等，各仿真模塊間存在著密切關(guān)聯(lián)。運(yùn)用Matlab Simulink建模工具，采用模塊化程序設(shè)計(jì)，允許用戶(hù)根據(jù)任務(wù)的需要，建立多種子模型，各子模型間采用模塊拼裝方式，允許自由組合，通過(guò)界面交互問(wèn)答確定仿真模型、算法及輸出形式。力求在模型的組合、算法的選擇及氣動(dòng)力數(shù)據(jù)的插值方面提供靈活便利的操作方式[6]。

1　仿真模型的建立

1.1坐標(biāo)系變換關(guān)系

坐標(biāo)系包括：指北地理坐標(biāo)系、彈體坐標(biāo)系、速度坐標(biāo)系、游移方位坐標(biāo)、彈道坐標(biāo)系和風(fēng)速坐標(biāo)系，相互之間通過(guò)歐拉角方程轉(zhuǎn)換。坐標(biāo)系定義和各個(gè)坐標(biāo)系之間轉(zhuǎn)換關(guān)系詳見(jiàn)文獻(xiàn)[7]。

ip,、為彈體繞慣性坐標(biāo)系的轉(zhuǎn)動(dòng)角速度在彈體軸上的3個(gè)投影；M、M、Mz為氣動(dòng)力矩的分量；Mp、Mp、Mzp為發(fā)動(dòng)機(jī)推力偏心產(chǎn)生的力矩分量；ΔM、ΔM、ΔMz為外界干擾力矩的分量；ψ、?、γ為偏航角、俯仰角、滾動(dòng)角；為導(dǎo)彈彈體相對(duì)于平臺(tái)坐標(biāo)系的姿態(tài)速率。

1.3子模型的建立

1）控制系統(tǒng)模型?？刂葡到y(tǒng)通常采用的方案為PID調(diào)節(jié)規(guī)律：式（4）中：Uδ、Uδ、Uδz為舵機(jī)的輸入電壓值；ψPR、?PR為程序偏航角、程序俯仰角；、PR為巡航高度、程序高度；kij為控制參數(shù)集。

2）指令裝訂模型。指令裝訂模型是指用戶(hù)根據(jù)仿真任務(wù)的需要預(yù)先裝入的飛行方案，飛行方案參數(shù)由用戶(hù)確定。

3）地形匹配模型。根據(jù)仿真任務(wù)要求，需要采用地形匹配技術(shù)，匹配模型根據(jù)實(shí)際匹配計(jì)算參數(shù)修改指令裝訂參數(shù)。

4）舵機(jī)回路模型。舵機(jī)采取何種反饋信號(hào)由具體情況而定。舵信號(hào)經(jīng)舵回路應(yīng)給出控制導(dǎo)彈運(yùn)動(dòng)的舵偏角信號(hào)。

5）質(zhì)量/慣量模型。數(shù)字仿真用到的導(dǎo)彈總體參數(shù)主要包括：導(dǎo)彈質(zhì)量m，繞重心的轉(zhuǎn)動(dòng)慣量J、J、Jz，質(zhì)心位置cm、cm、Zcm。在滿(mǎn)油（藥）、空油（藥）狀態(tài)下，總體參數(shù)見(jiàn)表1。

表1　導(dǎo)彈總體參數(shù)表Tab.1 Parameter table of missile overall

在飛行中要確定重心位置和轉(zhuǎn)動(dòng)慣量的變化規(guī)律，依線(xiàn)性變化可以得出：

質(zhì)量變化由一、二級(jí)裝藥（燃油）秒消耗量決定：式（7）中：t1為助推器脫離時(shí)間；te為飛行終止時(shí)間；m0為起飛質(zhì)量；mt1為拋掉助推器后的質(zhì)量；R為固體裝藥秒消耗量；P為燃油秒消耗量。

6）氣動(dòng)力模型。導(dǎo)彈在飛行中所受的氣動(dòng)力有升力、阻力與側(cè)向力，作用在導(dǎo)彈上的氣動(dòng)力矩為俯仰力矩Mz、偏航力矩M及滾轉(zhuǎn)力矩M。同時(shí)，仿真計(jì)算中要計(jì)算鉸鏈力矩Mh,、Mh,、Mh,z。

在導(dǎo)彈氣動(dòng)力計(jì)算中，首先要明確參考長(zhǎng)度L，參考面積S。

阻力系數(shù)C、側(cè)向力系數(shù)C、升力系數(shù)Cz和轉(zhuǎn)動(dòng)力矩系數(shù)m主要由以下幾個(gè)方面構(gòu)成：導(dǎo)彈飛行中還包括動(dòng)導(dǎo)數(shù)鉸鏈力矩系數(shù)為：

7）發(fā)動(dòng)機(jī)模型。飛航式導(dǎo)彈一般采用渦噴（扇）發(fā)動(dòng)機(jī)，發(fā)動(dòng)機(jī)推力特性及燃油秒消耗量與飛行馬赫數(shù)Ma、高度h、發(fā)動(dòng)機(jī)轉(zhuǎn)速n、大氣溫度T、總壓恢復(fù)系數(shù)ρ0等參數(shù)有關(guān)：

若發(fā)動(dòng)機(jī)推力數(shù)據(jù)以離散數(shù)據(jù)形式給出，可以采用五維插值得到。

8）指數(shù)型大氣模型。采用指數(shù)型大氣密度變化模型，大氣密度為

若ρ0取海平面的大氣密度，則h0=0；H-為大氣標(biāo)量高度。

9）地球模型。將地球當(dāng)作一個(gè)旋成圓球體，圓球半徑為常數(shù)（RE=6371.221km），即假定在地球表面不同經(jīng)緯度的曲率半徑相同；地球自轉(zhuǎn)角速度；引力加速度；μ=3 .986×1014m3/s2表示引力常數(shù)與地球質(zhì)量之積，g0=9.81 m/s2表示海拔高度為0時(shí)的重力加速度，方向指向地心。

10）陣風(fēng)模型。陣風(fēng)模型采用“1 - cos”模型，陣風(fēng)速度的3個(gè)分量均可以用下面的模型描述，以垂直風(fēng)速為例：

m對(duì)象的各固有頻率。

11）海浪模型。采用典型一維PM海浪譜作為海浪模型[8]。

2　基于Matlab的程序設(shè)計(jì)

隨著計(jì)算機(jī)技術(shù)的快速發(fā)展，利用各種高級(jí)語(yǔ)言基本上都可以開(kāi)發(fā)出導(dǎo)彈飛行仿真平臺(tái)，但較為合適的語(yǔ)言應(yīng)該是Matlab、C或者Fortran等計(jì)算功能強(qiáng)大的語(yǔ)言。

本文采用的是Matlab提供的動(dòng)態(tài)系統(tǒng)仿真工具Simulink，它是眾多仿真軟件中功能最強(qiáng)大、最優(yōu)秀、最容易使用的一種[9]。在Simulink中，對(duì)系統(tǒng)建模將變得簡(jiǎn)單，而且仿真過(guò)程是交互的。因此，可以很隨意地改變仿真參數(shù)，并且立即可以得到修改后的仿真結(jié)果。另外，使用Matlab中的各種分析工具，還可以對(duì)仿真結(jié)果進(jìn)行分析和可視化。導(dǎo)彈空間彈道數(shù)學(xué)模型結(jié)構(gòu)框圖如圖1所示，基于Simulink的雷達(dá)模型如圖2所示。

圖1　導(dǎo)彈空間彈道數(shù)學(xué)模型結(jié)構(gòu)框圖Fig.1 Structure diagram of missile ballistic math model

圖2　基于Simulink的雷達(dá)模型Fig.2 Radar model based on the Simulink

3　模型校驗(yàn)及仿真結(jié)果

3.1模型校驗(yàn)

采用頻譜分析法對(duì)全系統(tǒng)模型直接進(jìn)行模型驗(yàn)證，即將相同初始條件下全系統(tǒng)仿真模型的運(yùn)行結(jié)果與真實(shí)彈道數(shù)據(jù)進(jìn)行比較，從而反應(yīng)出模型的準(zhǔn)確性和可信度[10-11]。對(duì)海浪和大氣紊流等擾動(dòng)模型的校驗(yàn)即對(duì)仿真得到的數(shù)據(jù)進(jìn)行譜估計(jì)，將得到的估計(jì)譜與理論譜進(jìn)行擬合，根據(jù)擬合結(jié)果來(lái)判斷仿真模型與理論模型的近似程度，見(jiàn)圖3、4。

圖3　彈-浪夾角為0時(shí)譜擬合情況Fig.3 Situation of spectral fitting when the angle between missile and waves to 0

圖4　彈-浪夾角為π/2時(shí)譜擬合情況Fig.4 Situation of spectral fitting when the angle between missile and waves to π/2

3.2仿真計(jì)算

利用文中算法，對(duì)某飛航式導(dǎo)彈進(jìn)行了多條彈道的仿真計(jì)算，并結(jié)合歷次的靶場(chǎng)試驗(yàn)進(jìn)行分析和研究，驗(yàn)證了算法；同時(shí)也對(duì)模型進(jìn)行了有效性驗(yàn)證[12]。圖5是相同初始條件下，導(dǎo)彈仿真結(jié)果與實(shí)際飛行彈道對(duì)照曲線(xiàn)。

圖5　仿真結(jié)果與實(shí)際飛行彈道對(duì)照曲線(xiàn)Fig.5 Contrast curve between simulation result and actual flight trajector

3.3仿真結(jié)論

通過(guò)多次的仿真計(jì)算，得出下述結(jié)論：

1）采用模塊化的彈道仿真模型建立方法，按面向應(yīng)用及面向?qū)ο笏枷朐O(shè)計(jì)可以實(shí)現(xiàn)軟件；

?

2）頻譜分析法在系統(tǒng)仿真模型校驗(yàn)過(guò)程中得到了成功應(yīng)用，結(jié)合靶場(chǎng)試驗(yàn)數(shù)據(jù)可以進(jìn)行算法驗(yàn)證和模型有效性驗(yàn)證；

3）利用Matlab Simulink進(jìn)行系統(tǒng)仿真，可以滿(mǎn)足試驗(yàn)分析的需求，而且在系統(tǒng)模型建立、程序設(shè)計(jì)、模型驗(yàn)證、仿真試驗(yàn)和數(shù)據(jù)分析處理等方面都提供了極大的便利。

4　結(jié)束語(yǔ)

對(duì)于飛航式導(dǎo)彈全彈道仿真來(lái)講，研究過(guò)程中要優(yōu)化各個(gè)功能模塊，提高仿真效率和仿真軟件的可用性，注重仿真環(huán)境的一體化，這樣可以對(duì)仿真資源進(jìn)行統(tǒng)一的管理。

在本方法的數(shù)學(xué)公式組織過(guò)程中，既參考了經(jīng)典的數(shù)字仿真理論和已有的成功經(jīng)驗(yàn)，也采用了近幾年來(lái)的一些新技術(shù)及新標(biāo)準(zhǔn)。模型設(shè)計(jì)務(wù)求方法詳盡、條理清楚。本仿真方法具有通用性，更新相關(guān)模塊就可對(duì)新型號(hào)飛航式導(dǎo)彈進(jìn)行彈道仿真計(jì)算。

參考文獻(xiàn)：

[1]李小龍.地地戰(zhàn)術(shù)導(dǎo)彈測(cè)試發(fā)控系統(tǒng)訓(xùn)練仿真器研究[D].哈爾濱：哈爾濱工業(yè)大學(xué)，2002. LIIAOLONG. Research of training simulator in tactical missile test control sstem[D]. Harbin：Harbin Institute of Technolog，2002.（in Chinese）

[2]郭會(huì)明，楊明，李雷，等.地空導(dǎo)彈訓(xùn)練仿真系統(tǒng)的設(shè)計(jì)與實(shí)現(xiàn)[J].系統(tǒng)仿真學(xué)報(bào)，2003，26（1）：69-87. GUO HUIMING，ANG MING，LI LEI，et al. Design and implementation of the training simulator for the surface- to- air missile[J]. Journal of Sstem Simulation，2003，26（1）：69-87.（in Chinese）

[3]郭希維，蘇群星，谷宏強(qiáng).導(dǎo)彈訓(xùn)練仿真系統(tǒng)中信號(hào)的采集與處理[J].兵工自動(dòng)化，2004，23（2）：4-6. GUOIWEI，SU QUNIN，GU HONGQIANG. Signal acquisition and processing in training simulation sstem for missile[J]. Ordance IndustrAutomation，2004，23 （2）：4-6.（in Chinese）

[4]向哲，李善高，邱發(fā)廷，等.反艦導(dǎo)彈靶場(chǎng)試驗(yàn)視景仿真技術(shù)[J].海軍航空工程學(xué)院學(xué)報(bào)，2010，25（2）：235-237.IANG ZHE，LI SHANGAO，QIU FATING，et al. Scene simulation technologof certain anti- ship missile in range test[J]. Journal of Naval Aeronautical and Astronautical Universit，2010，25（2）：235-237.（in Chinese）

[5]彭國(guó)圣，蘇洪潮，李續(xù)武.基于Vega的某地空導(dǎo)彈視景技術(shù)研究[J].戰(zhàn)術(shù)導(dǎo)彈技術(shù)，2003（5）：63-67. PENG GUOSHENG，SU HONGCHAO，LIUWU. Research on scene technologof certain air-to-ground missile based on vega[J]. Tactical Missile Technolog，2003 （5）：63-67.（in Chinese）

[6]張平. Matlab基礎(chǔ)與應(yīng)用簡(jiǎn)明教程[M].北京：航空航天大學(xué)出版社，2001：156-158. ZHANG PING. Matlab foundation and application of introductortutorial[M]. Beijing：Beihang UniversitPress，2001：156-158.（in Chinese）

[7]李宗吉，程善政，劉洋.蒙特卡洛模擬法計(jì)算航空自導(dǎo)深彈命中概率[J].彈箭與制導(dǎo)學(xué)報(bào)，2012，32（2）：22-24. LI ZONGJI，CHENG SHANZHEN，LIUANG. The calculation of aerial homing-depth charge hitting probabilitbmonte-carlo method[J]. Journal of Projectiles，Rockets，Missiles and Guidance，2012，32（2）：22-24.（in Chinese）

[8]顧文錦，葉顯武.導(dǎo)彈仿真中海浪的數(shù)字仿真及其動(dòng)畫(huà)生成[J].系統(tǒng)仿真學(xué)報(bào)，1998，10（1）：33-38. GU WENJIN，EIANWU. Digital simulation and animation generation of sea-wave in the simulation of thesubmarine-launched missile[J]. Journal of Sstem Simulation，1998，10（1）：33-38.（in Chinese）

[9]劉云南.基于Matlab和Simulink的導(dǎo)彈彈道仿真與應(yīng)用[J].彈箭與制導(dǎo)學(xué)報(bào)，2011，31（3）：183-185. LIUUNNAN. Model building and application of missile simulation based on matlab and simulink[J]. Journal of Projectiles，Rockets，Missiles and Guidance，2011，31 （3）：183-185.（in Chinese）

[10]周廣東.基于小波變換的非平穩(wěn)脈動(dòng)風(fēng)時(shí)變功率譜估計(jì)方法研究[J].工程力學(xué)，2013，30（3）：89-97. ZHOU GUANGDONG. Estimation method of evolutionarpower spectrum for non-stationarfluctuating wind using wavelet transforms[J]. Engineering Mechanics， 2013，30（3）：89-97.（in Chinese）

[11]廖瑛.系統(tǒng)建模與仿真的校核、驗(yàn)證與確認(rèn)（VV&A）技術(shù)[M].長(zhǎng)沙：國(guó)防科技大學(xué)出版社，2006：110-116. LIAOING. Technologabout verification validation and accreditation of the modeling and simulation[M]. Changsha：National Universitof Defense TechnologPress，2006：110-116.（in Chinese）

[12]吳海霞.蒙特卡羅方法在實(shí)際問(wèn)題中的應(yīng)用[J].太原師范學(xué)院學(xué)報(bào)，2009，8（1）：76-79. WU HAIIA. Monte carlo method and applications in solving special problems[J]. Journal of Taiuan Normal Universit，2009，8（1）：76-79.（in Chinese）

HUANG Cheng1,U Peng2, ZHAOi-wang3
（1. The 91550thUnit of PLA, Dalian Liaoning 116023, China; 2. Naval Equipment Department, Beijing 100083, China; 3. The 92270thUnit of PLA, Weihai Shandong 264309, China）

Abstrraacctt:: The trajectorsimulation software can help to analze and evaluate the flight eperimental results of the winged missile weapon sstem. The trajectorsimulation method was designed based on Matlab Simulink, the mathematical model of missile trajectorwith sidegrees of freedom was established, and the software was programmed and tested using the aerodnamic data and control parameters of the missile. The simulation results is proved the model's accuracand reliabil?it.

作者簡(jiǎn)介：黃成（1978-），男，工程師，大學(xué)。

收稿日期：2014-08-16；

DOI:10.7682/j.issn.1673-1522.2015.02.016

文章編號(hào)：1673-1522（2015）02-0169-05

文獻(xiàn)標(biāo)志碼：A

中圖分類(lèi)號(hào)：TP391.9

修回日期：2014-12-28