張進強，蔣夏軍

(南京航空航天大學計算機科學與技術學院，江蘇 南京 210016)

1 引言

長期以來我軍對炮彈的試驗，尤其是針對炮彈外彈道的試驗，仍以實彈射擊為主。但是，實彈射擊試驗本身存在著試驗周期長、成本高、樣本數(shù)量少等特點，已經(jīng)不能完全滿足當前信息化靶場試驗鑒定技術要求。當前對于炮彈外彈道的研究以數(shù)據(jù)處理為主，缺乏以動態(tài)實時直觀的仿真方式對外彈道進行分析和處理。隨著計算機軟硬件的發(fā)展及圖形處理能力的增強，可視化仿真技術漸漸在軍事仿真領域中得到廣泛的應用。

目前，在可視化仿真方面的研究受到越來越多的關注，國內(nèi)外已有的研究工作主要體現(xiàn)在可視化的分類和可視化的實現(xiàn)方面。前者主要從分類的角度來研究各種具體可視化的方法，后者主要從可視化實現(xiàn)的角度、可視化的真實感及機器的負載來討論。從可視化的分類角度展開的研究中，文獻[1]介紹了科學技術可視化、數(shù)據(jù)可視化、信息可視化和知識可視化等分類和研究比較。文獻[2]從文獻分類角度介紹了上述四種方法，并對四種方法進行了分析研究及方法評估。文獻[3]分別從處理對象、數(shù)據(jù)類型及可視數(shù)據(jù)分析技術對可視化進行分類，分析各自的區(qū)別及使用場合。上述僅對可視化方法進行研究比較，但都沒有對相關可視化技術進行深入研究。

在技術實現(xiàn)的角度展開的研究中，文獻[4]使用“開放圖形庫”(OpenGL)結合C++構建三維圖形開發(fā)環(huán)境中，是一個高性能的圖形開發(fā)軟件包，但開發(fā)難度大、周期長。文獻[5]在VS2008開發(fā)環(huán)境下使用WPF和3DS MAX實現(xiàn)虛擬現(xiàn)實技術。WPF作為一套統(tǒng)一圖形化系統(tǒng)，融合在.NET開發(fā)環(huán)境，開發(fā)簡單、快速，但不能滿足多平臺和通用性需求。文獻[6]結合Microsoft Visual Studio .NET2005和VTK開發(fā)了三維可視化系統(tǒng)，提出算法精度與繪制速度相互制約，VTK跨平臺、支持并行處理、執(zhí)行效率高但歸因于二次開發(fā)，無法操控實時數(shù)值仿真的數(shù)據(jù)。

如今，三維視景仿真技術日漸成熟，加之Matlab強大的數(shù)值運算能力和二、三維圖形顯示優(yōu)勢。針對現(xiàn)有工作存在的一些不足，本文提出對模型實現(xiàn)數(shù)值仿真驅(qū)動視景仿真，即通過數(shù)值仿真得到的飛行數(shù)據(jù)實時驅(qū)動二維、三維圖形相結合的視景仿真，以此來全面、直觀地觀察彈丸在空中的實時飛行姿態(tài)、運動軌跡及其落點的散布情況。

本文結構安排如下：第2節(jié)主要對仿真數(shù)學模型進行介紹，其中包括：外彈道數(shù)學建模和隨機風影響的落點散布數(shù)學模型；第3節(jié)主要介紹仿真的實現(xiàn)方法，包括數(shù)值仿真和視景仿真。在數(shù)值仿真部分，介紹了微分方程的各種數(shù)值計算方法，分析各種方法的精度和時耗。在視景仿真部分，創(chuàng)新性地提出圖表相結合的三維視景仿真，并分別介紹了在C#環(huán)境下引用CsGL.dll構建OpenGL三維圖形開發(fā)環(huán)境和Matlab&C#三維網(wǎng)格坐標軸的混編；第4節(jié)是應用案例；第5節(jié)是本文的總結。

2 數(shù)學模型介紹

以炮彈外彈道方程組為基礎對外彈道曲線進行擬合，以蒙特卡洛法為基礎進行炮彈落點散布值的求解。

2.1 外彈道數(shù)學建模

目前常用的外彈道方程組主要有質(zhì)點彈道方程組和剛體彈道方程組。其中，剛體彈道方程組雖然計算精度高，但其計算復雜，運算時間過長，影響系統(tǒng)反應時間，不能很好地滿足實時性要求。因此，在考慮計算精度和計算速度的基礎上，本文采用質(zhì)點彈道方程組。

2.1.1 彈丸質(zhì)心運動模型

由于彈丸在飛行過程中會受到起始擾動、陣風、推力偏心、質(zhì)量分布不均衡、彈重偏差及氣象偏差等各種因素的影響，從而使彈道軌跡及落點的預測變得復雜無序。為了使外彈道方程簡單易懂，作出如下假設：

(1)整個飛行過程中，彈丸的章動角δ恒為零；

(2)彈丸的外形和質(zhì)量分布是軸對稱體；

(3)射擊場所地面是平面；

(4)重力加速度g方向鉛直向下，且大小不變；

(5)忽略由于地球自轉(zhuǎn)而產(chǎn)生的作用在飛行炮彈上的科氏慣性力。

在上述基本條件假設下，根據(jù)外彈道理論，考慮彈道條件(初速、彈重和藥溫)和氣象條件(氣溫、氣壓、縱風和橫風)非標準時的彈丸質(zhì)心運動方程組，由下式表示[7]：

(1)

其中，t為彈丸飛行時間，x、y、z為彈丸離開炮口后在飛行過程中的坐標(x,y,z)的水平分量、垂直分量和橫向分量；c為彈丸系數(shù)；θ為速度V的水平傾角；Vx、Vy、Vz分別為彈丸速度V在水平、垂直和橫向方向上的分量；g為重力加速度；H1(y)為空氣密度函數(shù)；G(v,c1)為空氣阻力函數(shù)。積分的初始條件為：t=0,vx=v0cosθ,vy=v0sinθ,vz=0,x=y=z=0。

其中阻力函數(shù)為[7]：

G(vr,c1)=4.737×10-4vrcxon(Ma),

Ma=vr/c1

(2)

2.1.2 彈道條件

彈道條件非標準主要包括初速、彈重和藥溫等。彈道條件非標準可折合成初速和彈道系數(shù)兩個參量的變化問題，因此可直接應用標準條件下的彈丸質(zhì)心運動微分方程組，只是決定空氣彈道的三個參量Cb、V0、θ0中的Cb和V0分別改為Cb+ΔCb和V0+ΔV0,其中ΔCb和ΔV0為彈道條件非標準時的折合量。

2.1.3 氣象條件

氣溫非標準條件規(guī)定，非標準氣象條件下氣溫隨高度變化服從標準分布，結合我國地理條件，改進氣溫標準模型得到非標準氣象條件下氣溫隨高度變化的規(guī)律滿足如下模型[8]：

(3)

其中，Tvon為標準條件下的地面氣溫，Tv1表示各高度處的不符合標準定律的氣溫,ΔTv為實際溫度與標準定律對應的氣溫偏差；G1=6.328×10-3K·m-1,B1=1.172×10-6K·m-2。

當氣溫、氣壓不符合標準定律時，氣壓函數(shù)的表達式為[8]：

(4)

其中，P0為地面氣壓值；π1(y=9300)、π2(y=12 000)分別表示氣壓函數(shù)在9 300 m和12 000 m時的值；A1=230 K。此時的空氣密度函數(shù)用H1(y)表示，其計算公式[8]為：

(5)

同理，當氣溫不符合標準定律時，對應的聲速也不符合標準定律，其計算公式為：

(6)

所以，氣溫、氣壓非標準時的彈道方程可直接應用標準條件下的彈道方程進行計算，只是以H1(y)代替H(y)、以實際聲速C1(y)代替標準聲速C(y)而已。在式(1)結合式(2)～式(6)的基礎上，針對選定的高炮彈藥，在氣溫、氣壓、射角及橫側風一定時，通過數(shù)值計算方法，即可確定以時間t為自變量的任意點的射擊諸元。

2.2 落點散布數(shù)學建模

隨機干擾(如隨機風和起始擾動)是影響炮彈落點散步的重要因素，將風場作為典型的隨機過程，建立隨機風場的數(shù)學模型并結合外彈道數(shù)學模型，采用經(jīng)典的統(tǒng)計實驗法——蒙特卡羅法進行隨機風作用下炮彈落點散布的模擬代替射擊實驗，以隨機風影響下炮彈落點模擬結果的均值和方差反映實際射擊過程的散布特征，以此分析隨機風對炮彈落點散布的影響。

由于風力和風向千變?nèi)f化，無法事先預知，因而在實際工作中主要考慮橫風、側風和鉛直風，除特殊地區(qū)(山谷)外鉛直風一般較小不予考慮。

仿真思想主要是對于隨機條件，用符合正態(tài)分布的偽隨機數(shù)表示，然后將偽隨機數(shù)代入彈丸的外彈道方程中進行計算，最后通過二、三維坐標圖實現(xiàn)彈道飛行軌跡及落點的可視化。

隨機仿真可得各發(fā)炮彈的縱向和橫向落點坐標，分別記為(x1,x2,…,xn)和(z1,z2,…,zn)，n為仿真次數(shù)。炮彈的縱向、橫向平均落點坐標，散布中心的差值分別為[9]：

通過對仿真結果的統(tǒng)計分析，可以求出落點與目標的射程及側偏偏差的均值，依據(jù)射表查找相應的修正量，從而提高炮彈的射擊精度。

在正確的外彈道模型及隨機風的數(shù)學模型的基礎上，利用蒙特卡羅仿真技術預測彈丸落點的方法，在測得隨機風的實時狀態(tài)下，通過模擬系統(tǒng)得出彈丸落點的方差和數(shù)學期望，分析風對彈丸落點的影響，得出對炮彈射擊精度影響的強弱程度。對火炮而言，最重要的彈道諸元是炮彈的射程及側偏。其中炮彈的射程主要由初速、射角和縱風引起的，而側偏則由橫風引起。

3 仿真實現(xiàn)

3.1 數(shù)值仿真

求微分方程的解有兩種方法，對于一些典型的微分方程(如線性方程、某些特殊的一階非線性方程等)可以求出其解析解，并根據(jù)初值的條件給出表達式，但遇到復雜的微分方程模型時，解析解很難求出，因此只能求方程的近似解或數(shù)值解。目前微分方程常用的數(shù)值計算方法主要有三種:歐拉法、改進的歐拉法及龍格-庫塔法。

經(jīng)測試，將上述三種方法代入質(zhì)點彈道方程組，求得各自的運算時間及精確度，如表1所示。

Table 1 Comparison of numerical calculation methods

因此,彈道方程的求解采用四階龍格-庫塔法，仿真系統(tǒng)中步長h取0.05 s，在保證有較好穩(wěn)定性的前提下整個外彈道仿真過程的累積誤差在0.1 m以內(nèi)。而且四階龍格-庫塔法的計算速度較快，這使得外彈道的仿真過程具有較好的實時性，能較好地滿足軍事打擊高精度和實時性的需求。

3.2 三維空間視景仿真

3.2.1 OpenGL三維場景仿真

若想逼真地模擬實彈飛行場景，視景仿真模塊必須進行實時的數(shù)據(jù)更新，如實時的炮彈的速度、坐標、射角及偏流等，即用數(shù)值仿真來驅(qū)動視景仿真。但同樣也需要建立起相應的逼真三維場景，包括靶場及炮彈的飛行過程還有落地的爆炸場景等。由于圖形的繪制效率和繪制效果是相矛盾的，提高真實感，勢必就會降低仿真場景的繪制速度，從而降低了圖形的實時性。因此，只有在滿足實時性的基礎上，盡可能地獲得高質(zhì)量的繪制效果。

本系統(tǒng)建立炮彈三維空間模型，在屏幕上顯示具有三維特效的炮彈及靶場，在此基礎上進行炮彈的空中飛行姿態(tài)的動態(tài)顯示。具體實施步驟如下：首先利用3DMAX生成炮彈的3D模型；其次利用C#開發(fā)環(huán)境調(diào)用CsGL.dll構建OpenGL運行環(huán)境模塊；調(diào)用OpenGL建立炮彈的三維空間模型[10]，根據(jù)數(shù)值仿真計算得到的飛行仿真數(shù)據(jù)構建炮彈的空中飛行過程(飛行姿態(tài)、飛行軌跡)；依照貼圖坐標構建靶場背景。

3.2.2 C#構建OpenGL開發(fā)環(huán)境模塊

用第三方插件CsGL，關鍵是使用其提供的兩個DLL：csgl.dll和csgl.native.dll。在一個C#項目中，用“項目”→“添加引用”→“瀏覽”找到csgl.dll引用。在需要使用OpenGL函數(shù)的cs文件的聲明中,添加：using CsGL.OpenGL。并做以下重構：

protected override voidOnSizeChanged(EventArgse){ }/*窗口大小變化時調(diào)用*/

public override voidglDraw(){ }/*OpenGL初始化，繪制想要畫的圖形*/

protected override voidInitGLContext(){ }/*繪制環(huán)境初始化*/

3.2.3 炮彈模型的建立及繪制

采用3DMAX建立炮彈立體模型及模型的貼圖坐標。模型相對容易建立，為使模型具有真實的視覺效果，可以對模型加上各種獨有的材質(zhì)，如顏色、反色特性等效果。OpenGL載入3D模型文件時，因?qū)ο蟮耐獠啃螤钪饕强咳切蝸砟M的，所以需先定義一系列結構體來存放模型的點、三角形及貼圖的坐標和材質(zhì)等。然后，參照結構體讀入的信息，用OpenGL繪圖命令根據(jù)得到的各種結構體數(shù)組來進行繪制，并通過顏色、材質(zhì)及紋理等深化效果。載入后的模型的大小和位置跟預期的效果不一樣，這時需采用GL.glTranslate()平移函數(shù)、GL.glRotate()旋轉(zhuǎn)函數(shù)及GL.glScale()縮放函數(shù)來調(diào)整模型。

3.2.4 炮彈的飛行仿真

炮彈的飛行仿真在視景系統(tǒng)中模擬導彈的四自由度，包括三維空間位置坐標(X,Y,Z)及射角(θ)的實時變化。于是炮彈的飛行仿真可以通過添加定時器實時地更新坐標和射角并調(diào)用GL.glTranslate()平移函數(shù)、GL.glRotate()旋轉(zhuǎn)函數(shù)及SwapBuffer()雙緩沖函數(shù)來實現(xiàn)炮彈的飛行仿真過程。

3.2.5 靶場背景的繪制

背景環(huán)境主要包括天空和地面，本文主要采取了紋理映射的方法，選擇近似的圖片作為紋理，通過載入圖片、定義紋理及根據(jù)幾何坐標映射來得到三維場景。對于可視化仿真要求不太高的場景，可以直接使用紋理映射，操作方便且速度快。

3.3 Matlab三維空間坐標軸的繪制

Matlab是高級工程計算軟件，程序可移植性高，不僅具有強大的數(shù)據(jù)可視化功能，還涵蓋二維和三維的可視化、圖像處理、動畫及表達式作圖。C#語言簡單、功能強大、類型安全，可實現(xiàn)應用程序的便捷開發(fā)，但作圖操作繁瑣復雜。因此，C#和Matlab混編可以揚長避短。

Figure 1 Visual simulation of trajectory圖1 外彈道視景仿真

本系統(tǒng)主要利用Matlab繪制三維空間網(wǎng)格坐標軸并能動態(tài)顯示彈道的飛行軌跡，達到視景仿真與曲線模擬的同步。具體實施步驟如下：首先通過Matlab的m文件來編譯生成.NET的dll文件[11]；其次C#調(diào)用生成的dll文件，以數(shù)組的方式傳入三維空間坐標(X,Y,Z)及時間t，Matlab利用得到的坐標數(shù)據(jù)進行實時繪圖。

3.3.1 Matlab生成dll文件

在Matlab的Command Window中輸入deploytool，產(chǎn)生deployment tool配置窗口，選擇Matlab Builder NET選項卡，新建工程并添加m文件，編譯成功即可得到dll文件。dll文件為接口，它實際上不含任何實現(xiàn)，都是調(diào)用另一個.ctf文件來完成具體的功能。

3.3.2 C#調(diào)用dll文件

添加對MWArray.dll和生成dll文件的引用，通過實例化MWNumericArrayMX、MY、MZ和T，將數(shù)值仿真計算所得的空間坐標及對應的時間t以數(shù)組形式傳給MX、MY、MZ和T并將實例數(shù)組傳入Matlab中的繪圖函數(shù)。其中，MWNumericArray是MWArray和C#中數(shù)據(jù)的中間類，C#與Matlab之間傳遞參數(shù)需用MWNumericArray轉(zhuǎn)換。

4 應用案例

依據(jù)上述外彈道、落點散布模型，以某型火炮為例，使用C#語言編制彈道程序，結合Matlab混編實現(xiàn)三維坐標軸的作圖及數(shù)值計算，并用CsGL實現(xiàn)空間視景的仿真，可更直觀地了解整個飛行過程及落點的散布。

4.1 炮彈的數(shù)值仿真

數(shù)值仿真計算輸入數(shù)據(jù)包括：初始速度、空間坐標、初始攻角、溫度偏差及橫風側風等。

初始條件為:V0=900 m/s,x=y=z=0,θ0=45°，tem=3,wx=1,wz=1。經(jīng)計算，炮彈飛行80 s，射程達到20 785.87 m，側偏50.15 m，與所提供的實驗數(shù)據(jù)相符。

4.2 外彈道視景仿真

如圖1是根據(jù)上述初始條件經(jīng)過數(shù)值計算得到的飛行數(shù)據(jù)實時繪制的炮彈射擊場景，圖1的左側部分是炮彈射擊的視景仿真，右側是由Matlab構建的二、三維坐標軸同步曲線模擬。圖2根據(jù)實時的飛行數(shù)據(jù)，通過三維仿真模擬出炮彈的初始狀態(tài)以及15 s、37 s、60 s和結束時的狀態(tài)?？梢愿鶕?jù)炮彈的飛行姿態(tài)，判斷炮彈的飛行穩(wěn)定性,由圖2可看出，炮彈在飛行的過程中，沒有發(fā)生翻倒，因此炮彈的飛行較穩(wěn)定。

實際仿真時幀率是使用Fraps軟件進行監(jiān)測求得，約為32 fps，即大約31.2 ms完成一幀。而數(shù)值仿真部分使用龍格-庫塔法，外彈道方程組運算時間約占16 ms，視景仿真部分引用CsGL.dll在C#下構建OpenGL環(huán)境，三維處理時間約占11 ms。因此，該系統(tǒng)既取得了較高的運算精度，又滿足了實時性需求。

Figure 3 Fall point of shell圖3 炮彈的落點散布

Figure 2 Flight posture of shell圖2 炮彈飛行姿態(tài)

4.3 落點散布仿真

數(shù)值仿真計算輸入的數(shù)據(jù)包括：初速、射角、空間坐標、隨機風的分布參數(shù)及模擬次數(shù)。初始條件為：V0=900 m/s,θ0=45°，x=y=z=0,wx(μ,σ2)=(1,5),wy(μ,σ2)=(1,5),times=30。

如圖3所示，左側為仿真30次的結果，包括橫風、側風以及在此影響下的射程及偏流。右側為模擬30次射擊的二、三維落點散布圖，從落點散布圖可以看出，側風對炮彈落點的影響遠大于橫風；炮彈的射程越遠，側風的影響越大。

5 結束語

本文實現(xiàn)了基于數(shù)值仿真實時驅(qū)動的彈道可視化仿真平臺，借助此平臺對炮彈的飛行軌跡、實時運動姿態(tài)及落點散布進行了仿真研究。研究表明，彈道可視化仿真系統(tǒng)可對炮彈運動狀態(tài)進行實時的計算及在三維場景中實時顯示其運動狀態(tài)，較完整地反映了炮彈的運動規(guī)律，操作簡便，畫面流暢；落點散布仿真根據(jù)添加的隨機風條件，得出側風對炮彈落點的影響遠大于橫風，炮彈的射程越遠，側風的影響越大，為實彈射擊的修正提供參考。

通過添加內(nèi)彈道模塊及毀傷分析模塊，此仿真系統(tǒng)可以擴展為全彈道的仿真系統(tǒng)。此仿真方法也可用于火箭、導彈等武器的仿真。

[1] Zhang Zhuo,Xuan Lei,Hao Shu-yong.Study and comparison on visualization technology[J].Modern Electronics Technique,2010,33(17):133-138.(in Chinese)

[2] Liu Bo,Xu Xue-wen.Comparison of the visualization taxonomy[J].Journal of Intelligence,2008,27(2):28-30.(in Chinese)

[3] Zhu Yao-hua,Hao Wen-ning,Chen Gang.The summary of visualization theory[J].Computer Knowledge and Technology,2012(6):1.(in Chinese)

[4] Cai Hua,Lu Jing-chao,Xiao Gang.Real-time 3D animated flight simulation based on OpenGL[J].Computer Simulation,2007,24(10):28-30.(in Chinese)

[5] Yang Chao,Xu Xiao-long,Yu Kun,et al. 3D visual control system of underwater vehicle[J].Computer Simulation,2012(4):185-188.(in Chinese)

[6] Schroeder W J. The visualization toolkit users’ guide[M]. New York:Kitware Inc, 2001.

[7] Guo Xi-fu,Zhao Zi-hua.Theory and application of fire control ballistic model[M].Beijing:National Defence Industry Press，1997.(in Chinese)

[8] Zhao Xin-sheng,Shu Jing-rong.The theory and application of ballistic calculation[M].Beijing:Weapon Industry Press，2006.(in Chinese)

[9] Wang Xin,Zhang Feng-ge,Yao Jun,et al.Analyze effect of random wind on shell distribution of rocket using Monte-Carlo method[J].Journal of Projectiles, Rockets, Missiles and Guidance,2009,29(1):198-201.(in Chinese)

[10] The studio of peace dove.OpenGL advanced programming and visualization system[M].Beijing:China Water Power Press,2005.(in Chinese)

[11] Phan J.MATLAB+C#+Book[M].Colorado:LePhan Publishing, 2004.

附中文參考文獻：

[1] 張卓，宣蕾，郝樹勇.可視化技術研究與比較[J].現(xiàn)代電子技術，2010,33(17):133-138.

[2] 劉波，徐學文.可視化分類方法對比研究[J].情報雜志，2008,27(2):28-30.

[3] 朱耀華，郝文寧，陳剛.可視化技術簡述[J].電腦知識與技術，2012(6):1.

[4] 蔡華，盧京潮，肖剛.基于OpenGL的飛行實時三維動畫仿真[J].計算機仿真，2007,24(10):28-30.

[5] 楊超，徐筱龍，余琨，等.水下機器人三維可視化控制系統(tǒng)[J].計算機仿真，2012(4):185-188.

[7] 郭錫福，趙子華.火控彈道模型理論及應用[M].北京：國防工業(yè)出版社，1997.

[8] 趙新生，舒敬榮.彈道解算理論與應用[M].北京：兵器工業(yè)出版社，2006.

[9] 王欣，張鳳閣，姚俊，等.蒙特卡羅法分析隨機風對火箭彈落點散布的影響[J].彈箭與制導學報，2009,29(1):198-201.

[10] 和平鴿工作室.OpenGL高級編程與可視化系統(tǒng)開發(fā)[M].北京：中國水利水電出版社，2005.