单炮多发同时弹着的局部插值算法与并行加速

程李东,姜 毅,牛钰森

(北京理工大学 宇航学院,北京 100081)

20世纪80年代后期,美国人机工程研究所针对野战炮兵射击提出了单炮多发同时弹着法(multiple round simultaneous impact,MRSI),通过调整火炮的发射装药和射角,使连续发射的多发炮弹同时击中目标[1-3],其技术构想如图1所示[2]。

目标受到火炮的首群炮弹打击后会迅速采取规避或保护措施,如车辆机动、人员卧倒或进入掩体等,后续射击对目标的毁伤概率大大降低。统计资料表明,暴露步兵遭受炮击后2 s内有一半就地卧倒,8 s内全部转为卧姿或隐蔽状态,后续炮弹杀伤力将降低50%～80%[4]。采用MRSI射击方法可在不增加火炮数量的条件下将首轮打击的炮弹密集度增大数倍,大幅提升对目标的毁伤概率,进而提高炮兵部队的作战效率和生存能力。

MRSI射击是现代火炮的发展趋势,国外先进的自行加榴炮均具有同时弹着功能[5]。美国在1988～1993年间实现了AS90,M109A6等155 mm自行火炮的四弹同时弹着,同时测试和对比了标准装药模块和液体装药的MRSI性能,1995年研究了先进牵引式火炮系统(ATCA)的MRSI性能,部分火炮的四弹MRSI可覆盖其打击范围的85%[3]。2001年,陆欣等将MRSI的概念引进国内[2],祝军利等给出了国内某型火炮在13～30 km射程范围内的MRSI弹道分布[5],随后李开龙、程恭等研究了MRSI射击方法在大口径舰炮上的应用,提出了弹道解算流程并分析了射击效力[6-7],近年来相关研究逐渐转向增程修正弹等新型可控炮弹[8-10]。

自动化火炮的火控系统担负目标定位与跟踪、气象数据的测定和接收、射击诸元的解算、瞄准与射击控制等任务,其中射击诸元的解算直接关系到打击的速度和精度[11]。在机动作战的条件下,战场情形瞬息万变,这对弹道计算机的解算速度提出了较高的要求,如装备在我国69-II式坦克上的简化轻型坦克火控系统要求弹道计算时间小于1 s[1]。现有的MRIS弹道算法[6-7]为满足计算耗时的要求,均采用稀疏射击序列计算,再结合线性插值的方式确定装药和射角。由于弹道方程具有非线性,稀疏序列插值将造成较大的弹道偏差。

为解决原有算法弹道偏差较大的问题,本文提出一种局部插值算法(local interpolation algorithm,LIA),并基于MPI编程模型实现多核并行计算,在大幅提升射角计算精度的同时满足战场对火控计算时间的要求。最后将研究结果应用于一定条件下5～15 km射程内的MRSI弹道计算,以检验算法的有效性和高效性。

1 数学模型

炮兵实际使用的弹道方程组是高维微分方程组,如美军使用的修正质点弹道方程组[3,12],该弹道方程组考虑了地球自转和弹体转动带来的影响。本文重点研究弹道算法和并行加速方法,因此仅以标准条件下二维质点弹道方程[4,13]为例,但研究结果可以直接应用于修正质点弹道方程、旋转稳定炮弹弹道方程[14]、高空远程弹箭精确弹道模型[15]等复杂弹道方程的解算。

1.1 标准条件下二维质点弹道方程

二维质点弹道方程考虑了空气阻力,同时计及高度变化对大气物理特性的影响：

虚温τ按标准气象条件计算,g=9.81 m/s2,ρ0,.N=1.206 kg/m3,d为榴弹弹径,阻力系数cx0,N(Ma)由阻力定律确定,i为所选阻力定律下的弹形系数,m为弹体质量,p为大气压强。

1.2 真空二维弹道斜射程公式

真空弹道为抛物线弹道,不考虑空气阻力,其射程总是大于同一装药在同一射角下的空气弹道射程。根据外弹道学,真空弹道斜射程公式为[13]

式中:v0为炮口初速,α为瞄准角,ε为炮目高低角。

最大射程:

2 约束条件

2.1 射程重叠区间与装药分级

榴弹炮为达到射程要求,需要对发射装药进行一定的分级。从同时弹着射击的角度来说,分级越多可实现同时弹着的弹道数量越多。但同时又必须尽可能简化战场操作,并考虑补给等问题,因此各级装药的射程重叠量通常为8%～12%,对于30 km射程的火炮,其装药常分为7级[4-5]。

2.2 最小发射时间间隔

2次射击之间需要完成弹药的再装填、关闭舱门、重新瞄准等操作,且射击之前需要保证炮身冷却至容许温度以下[16]。不同型号的火炮完成这些操作所需的时间不同,按常规火炮8 min-1的射速计算[5],2次射击之间的最小时间间隔为Δtmin=8.75 s。这一间隔将作为弹道筛选的条件,即相邻的2条同时弹着弹道的飞行时间差不小于Δtmin。

2.3 最小与最大射角

火炮的最小与最大射角受外弹道设计、火炮射界和炮身纵向倾角的共同限制,如图2所示。

图2中,OA与水平线平行,OB为目标视线,OC与炮台平行;θ为火炮射角,φ为炮身纵向倾角,为射击高低角。

一方面根据弹丸飞行动力学特性和跟随性确定最大射高,进而确定理论最大射角θt,max[5]。另一方面,火炮射界受炮架结构的限制,射界越大炮架质量越大,榴弹炮的最小高低角γmin一般取-3°～-5°,最大高低角γmax可达70°甚至更大[1]。据此,实际火炮的射角范围为

γmin+φ≤θ≤min{θt,max,γmax+φ}

2.4 命中精度

以人员杀伤为例,榴弹对人员的杀伤有2种方式:一是爆炸后飞散的破片造成人体组织穿透、断离或撕裂,二是爆炸冲击波对人体造成冲击损伤[17]。2种杀伤方式的杀伤效果都随爆炸距离的增大而迅速减小,前者是由于受空气阻力影响,破片速度随飞行距离增大而迅速减小,后者是由于冲击压力在空气中迅速衰减。可对人员造成杀伤的最大爆炸距离即榴弹对人员的杀伤半径。

由于弹道数学模型、阵地位置与高度测量、敌方距离和方位测量、气象修正、初速修正以及数值求解过程中均存在误差[18],实际采用MRSI方式射击时一轮炮弹几乎不可能同时到达目标点,最先和最后达到的炮弹的落地时差表征了MRSI弹道的时间精度。为保证杀伤效率,该时差应尽可能小。

3 局部插值算法

确定同时弹着射击方案的思路是通过求解弹道方程组,选择出覆盖给定射程的装药,然后确定出各装药对给定射程的射角,最后根据限制条件得出射击方案[2]。为满足计算耗时的要求,已有文献[6-7]的算法是给定一个较稀疏的离散射击序列:

(v0,j,θj,I),j=1,2,…,M,I=0,1,2,…,Nj

式中:M为装药种数,Nj+1为第j号装药待计算的射角数量。求解该射击序列对应的射程序列,再通过线性插值的方式确定可以命中目标的射击序列:

(v0,j,θj,b),b=0,1,2,…,n

为克服上述算法弹道偏差较大的问题,可采用局部插值算法。对于第j号装药：

①根据真空弹道斜射程公式判断其最大射程是否可以覆盖目标；

②在火炮最大容许射角范围[θmin,θmax]内采用二分法确定真空弹道达到或超过目标射程的射角范围[θj,min,θj,max],θj,min和θj,max作为空气弹道计算的射角下边界和上边界；

③给定某较小的射角增量Δθ,从θj,min开始向上逐个射角求解弹道。即给定(v0,j,θj,min+bΔθ),采用四阶Runge-Kutta方法[19]求解弹道微分方程组,若存在命中目标(弹道落点与目标的距离不超过杀伤半径)的弹道,则直接保存其弹道参数(v0,j,θj,min+bΔθ,tj,g)；若不存在命中目标的弹道,则计算直到第1条超过目标射程的弹道时停止向上遍历,对最后2条弹道的参数做线性插值,得该装药对目标的低伸弹道；

④从θj,max开始向下逐个求解弹道,直到找到命中目标的弹道或第1条超过目标射程的弹道,结束对该装药的弹道求解,采用与③相同的方法得到该装药的曲射弹道；

⑤对全部装药求得的所有弹道按弹道飞行时间tj,i递增排列,去除与前1条弹道的飞行时间差小于Δtmin的弹道,得到MRSI弹道的射击序列。

记低伸弹道射角为θj,r,曲射弹道射角为θj,q,采用局部插值算法可将第j号装药求解的射角范围限制于θj,use=[θj,min,θj,r]∪[θj,q,θj,max]。在本文第5节的计算条件下,第5号装药对10 km射程的θ5,use只占火炮容许射角范围的4.58%,在计算量相当条件下可以将原有算法的稀疏弹道序列加密20倍左右。第7号装药在45 km附近的弹道偏差小于1 m,弹道计算精度得到了很大提升。第5号装药对10 km射程的边界射角弹道曲线如图3。θ5,use中包含61条弹道,与文献[7]所给计算序列包含的射角数量51相差不大。

图3中弹道曲线Ⅰ和Ⅱ为抛物线低伸和曲射弹道,射角分别为θj,min和θj,max；曲线Ⅲ～Ⅶ为空气弹道,其中曲线Ⅲ和Ⅳ为低伸弹道,射角分别为θj,min和θj,r；曲线Ⅴ和Ⅵ为曲射弹道,射角分别为θj,max和θj,q；曲线Ⅶ为射角在[θj,r,θj,q]范围内的弹道。

采用局部插值算法实现MRSI弹道计算的流程如图4所示。

4 并行计算

在高性能并行计算中,基于消息传递的MPI(message passing interface)编程模型是科研和工程领域中的事实标准[20],本文采用MPI实现MRSI弹道的并行计算。在进行大量空气弹道解算之前,需要先根据抛物弹道斜射程公式判断各装药是否能达到目标射程,并确定可达到目标射程的装药的空气弹道上、下边界,所有装药的空气弹道计算结束后需要对命中目标的弹道排序、筛选以确定MRSI弹道射击序列,以上任务计算量较小,单个CPU便可迅速完成；大量空气弹道的解算则由多个CPU并行完成。因此本文采用MPI模型中的1↔N集合通信模型[21],在所有参与计算的CPU中,存在一个Host(主CPU),负责数据前处理、任务分配和最终的弹道序列生成；其余均为Device(设备CPU),只负责大量空气弹道的并行解算。本文研究所用计算平台的CPU为6核12线程,型号为Intel(R)Core(TM)i7-5820k,主频3.30 GHz；内存32 GB；操作系统为Windows Server 2012R2 Standard。

4.1 任务分配方案

为确定MRSI射击序列,需对局部插值算法确定的射角范围内的密集弹道序列进行解算,为此可有以下3种任务分配方案供选择(设有k个Device)。

① 每个Device负责一种装药的弹道解算,M种装药需要k=M；

② 依次计算各装药的弹道,每个装药的θj,use被均分为k段,每个Device计算其中一段；

③ 依次计算各装药的弹道,各Device计算的弹道按插空分配,设装药j需要计算的射角序列为

[θj,0,θj,1,θj,2,…,θj,Nj]

则第h个Device负责计算的射角序列为

[θj,h,θj,h+k,θj,h+2k,…,θj,h+fk],h+fk≤Nj

式中:f为满足h+fk≤Nj的最大自然数。

方案①为粗粒度并行,方案②和方案③为细粒度并行。为减少计算耗时,需要尽可能提高并行算法的并行度。由于各装药的炮弹初速差别较大,θj,use占火炮最大容许射角范围的比例相差也会较大,所以按照方案①分配时,各Device的计算量明显不均匀,算法的并行度较低。根据外弹道学,曲射弹道的弹道飞行时间总是大于低伸弹道。按方案②分配计算任务时,由于各Device求解弹道方程的时间步长为统一值,曲射弹道的时间步数量总会大于低伸弹道,算法并行度低。按方案③分配任务时,由于同一装药相邻射角的弹道差别不大,各Device之间可实现最大程度的并行。综合上述分析,任务分配应采用方案③。

图5为采用不同的任务分配方案计算10 km射程的MRSI弹道所消耗的计算时间,数值实验结果与上述分析相吻合。计算耗时为tc,计算中保证Device数量k=M,其余计算条件见第5节。

4.2 并行效率

由图6可知,随着Device数量的增加,并行效率在递增,但递增的速度减缓。这是由于随着Device数量的增加,弹道计算所消耗的时间减少,数据处理消耗的时间占总耗时的比重上升,且用于通信和数据处理的时间消耗会增加。由此可以预测,存在某个Device数Mη使得并行效率取最大值,之后再增加Device的数量,并行效率将降低,实际系统的Mη可由实验确定。从并行效率、硬件复杂程度、可靠性和经济性等角度考虑,弹道计算机的Device数量不应该超过Mη。

5 应用实例

为验证局部插值算法和并行算法的有效性,取如下条件做MRSI弹道仿真计算。

①榴弹弹径d=122 mm,cx0,N(Ma)采用43年阻力定律,弹形系数i=1.035；

②各级装药的炮弹初速如表1所示。

表1 装药初速表

③两次射击之间的最小时间间隔Δtmin=8.75 s；

④火炮最大容许射角范围为5°～85°；

⑤命中精度：榴弹杀伤半径为10 m,MRSI弹道时间精度为0,即计算状态下所有炮弹同时落在距离目标10 m以内；

⑥射程范围5～15 km,计算间隔为1 km,目标相对我方阵地高度为100 m；

⑦Δθ=0.06°,Runge-Kutta法求解弹道方程时,为保证弹道计算精度,时间步长取1 ms。

计算结果：15 km射程计算耗时最大,为14.16 s,与文献[2]相符,所得MRSI弹道数量Q随射程X的分布如图7。文献[3]的弹道数量分布在2～6条,与本文相符。

图7中,10 km射程的MRSI弹道射击序列见表2,j为装药号,θ为射角,t为弹道飞行时间,ts为射击时刻。

表2 10 km射程MRSI弹道射击序列

10 km射程的MRSI弹道曲线如图8,其中编号Ⅰ～Ⅴ对应表2中的序列号1～5.

6 结论

本文分析了已有的MRSI弹道算法和弹道分布规律,在此基础上做了以下研究工作：

①提出了一种局部插值算法,可大幅提高弹道计算精度,同时引入MPI并行计算,保证计算耗时满足战场要求；

②在并行算法设计中,分析了不同任务分配方式的并行度,结合数值实验结果,选择并行度最高的第③方案；

③根据并行效率与CPU核心数量的关系,指出弹道计算机的CPU数量不应该超过Mη；

④应用局部插值算法结合并行计算,在一定条件下计算了5～15 km射程的MRSI弹道,验证了算法的有效性和高效性。

虽然本文的研究对象只是普通自动化榴弹炮,但局部插值算法和并行计算方法可直接推广到基于瞄准点排布的面积杀伤法[7]和一维修正弹[8-10]的MRSI弹道计算中。近年来,通用图形处理器(general purpose graphics processing units,GPGPU)并行加速技术已经在计算传热学[22]、计算流体力学[23]以及信号处理[24]等多个领域取得了应用。相比于CPU,GPGPU更适合于进行大规模重复性计算,未来可以考虑将GPGPU并行加速应用于MRSI弹道求解,以解决更复杂的新型可控炮弹[25]MRSI弹道解算问题。

[1] 谈乐斌,张相炎,管红根,等. 火炮概论[M]. 北京:北京理工大学出版社,2005.

TAN Lebin,ZHANG Xiangyan,GUAN Honggen,et al. Introduction to artillery[M]. Beijing:Beijing Institute of Technology Press,2005.

[2] 陆欣,周彦煌. 单炮多发同时弹着的数学模型和数值分析[J]. 火炮发射与控制学报,2001(3):5-6,27.

LU Xin,ZHOU Yanhuang. Mathematic model and numerical simulation of single gun multiple round time-on-target technology[J]. Gun Launch & Control Journal,2001(3):5-6,27. (in Chinese)

[3] KOGLER T M. Single gun,multiple round,time-on-target capability for advanced towed cannon artillery[R]. Maryland:Army Research Lab Aberdeen Proving Ground,1995.

[4] 谢黎焱,廖瑞,王雪琴. 信息化条件下单炮多发同时弹着研究[J]. 指挥控制与仿真,2009,31(4):30-32,36.

XIE Liyan,LIAO Rui,WANG Xueqin. Study on single cannon multiple round time-on-target under informationizated condition[J]. Command Control & Simulation,2009,31(4):30-32,36. (in Chinese)

[5] 祝军利,高效生,侯晓峰. 自行加榴炮多发同时弹着技术[J]. 火炮发射与控制学报,2005(1):6-8.

ZHU Junli,GAO Xiaosheng,HOU Xiaofeng. Technique of MRSI of self-propelled cannon/howitzer[J]. Gun Launch & Control Journal,2005(1):6-8. (in Chinese)

[6] 李开龙,石章松,龚驰,等. 大口径舰炮多发同时弹着射击效力分析[J]. 舰船电子工程,2011,31(9):31-33,42.

LI Kailong,SHI Zhangsong,GONG Chi,et al. Analysis of time-on-target fire effectiveness of large caliber naval gun[J]. Ship Electronic Engineering,2011,31(9):31-33,42. (in Chinese)

[7] 程恭,石章松,王航宇. 基于瞄准点排布模型的多发同时弹着射击诸元解算[J]. 舰船电子工程,2011,31(11):29-32,42.

CHENG gong,SHI Zhangsong,WANG Hangyu. Firing data solutions for multiple rounds simultaneous impact based on aiming point distribution model[J]. Ship Electronic Engineering,2011,31(11):31-33,42. (in Chinese)

[8] 黄义,汪德虎,汪汇川,等. 舰炮发射一维修正弹多发同时弹着研究[J]. 弹箭与制导学报,2012,32(5):127-129.

HUANG Yi,WANG Dehu,WANG Huichuan,et al. The research on MRSI of one-dimensional trajectory correction projectile fired by shipborne gun[J]. Journal of Projectiles,Rockets,Missiles and Guidance,2012,32(5):127-129.(in Chinese)

[9] 刘剑威,王海川. 增程修正弹单炮多发同时弹着火控技术研究[J]. 指挥控制与仿真,2012,43(1):70-72,77.

LIU Jianwei,WANG Haichuan. Research on method of multiple rounds simultaneous impact of extended range and trajectory correction projectile[J]. Command Control & Simulation,2012,43(1):70-72,77. (in Chinese)

[10] 石章松,傅冰,胡献君,等. 基于增程修正弹的同时弹着火控机理[J]. 海军工程大学学报,2013,25(3):7-12.

SHI Zhangsong,FU Bing,HU Xianjun,et al. Time-on-target fire control mechanism based on extended-range trajectory correction projectile[J]. Journal of Naval University of Engineering,2013,25(3):7-12. (in Chinese)

[11] 孟灿,毛征,孟博,等. 弹丸飞行时间与射击诸元关系分析[J]. 兵工自动化,2016,35(11):24-27.

MENG Can,MAO Zheng,MENG Bo,et al. Analysis on relation of projectile flight time and firing data[J]. Ordnance Industry Automation,2016,35(11):24-27. (in Chinese)

[12] LIESKE R F,REITER M L. Equations of motion for a modified point mass trajectory[R]. Maryland:Army Ballistic Research Lab Aberdeen Proving Ground,1966.

[13] 韩子鹏. 弹箭外弹道学[M]. 北京:北京理工大学出版社,2008:72-77.

HAN Zipeng. Exterior ballistics of rockets and missiles[M]. Beijing:Beijing Institute of Technology Press,2008:72-77. (in Chinese)

[14] 钱明伟,王良明,郭锡福. 火炮武器高原射击时的弹道特性研究[J]. 弹道学报,2009,21(4):21-25.

QIAN Mingwei,WANG Liangming,GUO Xifu. Ballistic analysis for artillery systems firing on plateau[J]. Journal of Ballistics,2009,21(4):21-25. (in Chinese)

[15] 李臣明,舒敬荣. 高空远程弹箭精确弹道模型[J]. 弹道学报,2008,20(4):8-11.

LI Chenming,SHU Jingrong. Precise exterior ballistic model for long-range missile versus high altitude[J]. Journal of Ballistics,2008,20(4):8-11. (in Chinese)

[16] 朱文芳,王育维,郭映华,等 某火炮身管温度仿真计算及其影响因素分析[J]. 火炮发射与控制学报,2016,37(4):58-62.

ZHU Wenfang,WANG Yuwei,GUO Yinghua,et al. Temperature simulation calculation and analysis of influential factors of a certain gun barrel[J]. Journal of Launch & Control,2016,37(4):58-62. (in Chinese)

[17] 王志军,尹建平. 弹药学[M]. 北京:北京理工大学出版社,2005:106-116.

WANG Zhijun,YIN Jianping. Ammunition[M]. Beijing:Beijing Institute of Technology Press,2005:106-116. (in Chinese)

[18] 武瑞文,王兆胜. 自行火炮武器系统射击精度研究[J]. 兵工学报,2004,25(4):407-409.

WU Ruiwen,WANG Zhaosheng. A study on the firing accuracy of self-propelled artillery systems[J]. Acta Armamentarii,2004,25(4):407-409. (in Chinese)

[19] BURDEN R L,FAIRES J D. Numerical analysis[M]. Boston,Massachusetts:Brooks/Cole,2011:282-291.

[20] 龚春叶,包为民,汤国建,等. 航天领域高性能并行计算研究进展[J]. 计算机工程与科学,2014,36(9):1629-1636.

GONG Chunye,BAO Weimin,TANG Jianguo,et al. Recent progress in high-performance parallel computing of the aerospace aera[J]. Computer Engineering & Science,2014,36(9):1629-1636. (in Chinese)

[21] 张武生,薛巍,李建江,等. MPI并行程序设计实例教程[M]. 北京:清华大学出版社,2009:152-157.

ZHANG Wusheng,XUE Wei,LI Jianjiang,et al. Examples of MPI parallel programming[M]. Beijing:Tsinghua University Press,2009:152-157. (in Chinese)

[22] OUYANG A,TANG Z,ZHOU X,et al. Parallel hybrid PSO with CUDA for lD heat conduction equation[J]. Computers & Fluids,2015,110:198-210.

[23] PICKERING B P,JACKSON C W,SCOGLAND T R W,et al. Directive-based GPU programming for computational fluid dynamics[J]. Computers & Fluids,2015,114:242-253.

[24] LOBEIRAS J,AMOR M,DOALLO R. BPLG:a tuned butterfly processing library for GPU architectures[J]. International Journal of Parallel Programming,2015,43(6):1078-1102.

[25] 克里斯托弗·福斯,库宗波. 智能弹药:野战炮兵精确制导弹药[J]. 国外坦克,2016(10):20-26,31.

FUSS C,KU Zongbo. Intelligent munition:precision-guided munitions of field artillery[J]. Foreign Tanks,2016(10):20-26,31. (in Chinese)