【数学建模】关于导弹追踪轨迹分析的数学模型

数学建模竞赛承 诺 书我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。我们参赛选择的题号是（从A/B中选择一项填写）： A 我们的队号为： 07 参赛队员：1. 张建彬 2. 陆丽娜 3. 王晓龙 指导教师或指导教师组负责人： 数模组 日期： 2009 年 8 月 11 日评阅编号（由评阅老师评阅前进行编号）：数学建模竞赛编 号 专 用 页评阅编号：评阅记录：评阅人评分备注导弹追踪轨迹分析模型内容摘要 本题分析了空对空导弹的追踪问题以及地对空导弹的追踪问题。对于问题1，本文建立了二维平面上的导弹追逐模型。利用在任何时刻导弹的飞行方向指向敌机的位置得出微分方程，利用导弹飞行轨迹的弧长与敌机飞行路程之比等于两者的速度之比列出积分方程。再利用初值条件，并经过严格的数学公式推导求解出导弹追踪敌机的轨迹方程。发射I型空对空导弹击毁敌机的条件：其中为敌机飞行的速度，为型空对空导弹飞行的速度，击中点为对于问题2，由于导弹是由地面向空中发射，而敌机是在空中飞行，因此本文建立了三维空间上的导弹追逐模型，并把该三维空间上的导弹追逐问题转化为二维平面上的导弹追逐问题，而在转化过程中不改变各点之间的相对位置。运用问题1的解决方法求解得出II型地对空导弹追踪敌机的轨迹方程。发射II型地对空导弹击毁敌机的条件为： 其中为敌机飞行的速度，为型地对空导弹飞行的速度，击中点为对于问题3，我们用Matlab编程语言，采用Euler法迭代并进行改进得出结果。关键字：导弹发射追踪 二维模型 三维分析 Matlab编程 欧拉 迭代1问题的提出某边防导弹基地的雷达发现位于其正东 N 公里处有一家来犯敌机正欲逃往正北方向 M 公里处的安全区。该基地的I型空对空追踪导弹和II型地对空追踪导弹均可针对目标随时自动调节追踪方向，截击敌机。但敌机一旦进入安全区后，由于电子干扰作用，I型、II型导弹均将失去追踪目标，无法将敌机击毁。问题1：如此时恰有一架携有I型空对空追踪导弹的、与敌机处于同一飞行高度的巡航飞机在空中，基地即下令巡航飞机发射I型追踪导弹击毁敌机。试确定导弹追踪敌机的轨迹，并在适当的假定下给出发射该种导弹击毁敌机的条件；问题2：如此时在基地即发射II型地对空追踪导弹去击毁敌机，试确定导弹追踪敌机的轨迹，并在适当的假定下，及发射该种导弹击毁敌机的条件；问题3：若导弹的速度可在发射前根据需要设定，对于不同的N、M取值，编写计算机程序（语言不限），利用计算得到的数据说明怎样的发射速度可确保击毁敌机。2.基本假设(1).假设导弹以及敌机的运动为质点运动。(2).假设导弹是匀速率运动的。(3).假设敌机是匀速运动的。(4).假设导弹与敌机的运动速度跟风速和没有关系，在运动的过程中忽略重力及空气阻力的影响。(5).假设导弹的射程无限远。(6).假设导弹没有发射时间误差。3符号的约定O点:边防导弹基地的雷达所在位置坐标I点:巡航飞机的初始位置坐标P点():型空对空追踪导弹所在的位置坐标型空对空追踪导弹得轨迹方程型地对空追踪导弹得轨迹方程（转化为二维空间后）A点敌机所在的初始位置坐标(模型一)A点敌机所在的初始位置坐标(模型二)Q点(N,): 敌机所在的位置坐标(模型一)Q点:敌机所在的位置坐标(模型二)S点: 击中位置坐标（临界状态）（模型一）S点: 击中位置坐标（临界状态）（模型二）: 敌机飞行速度: 型空对空追踪导弹飞行速度:型地对空追踪导弹飞行速度敌机飞行的高度4.问题1的解答4.1模型的分析本问题的难点是需要假设巡航飞机的初始位置坐标，导弹的飞行速度以及敌机的飞行速度。由于是求解在适当的假定下的发射该种导弹击毁敌机的条件，因此需从临界状态出发，即击中点为点S(N,M)。且需用差分法代替微分法求解。4.2模型的假设(1).击中临界点为S点：(N,M)。(2).导弹与敌机的轨迹始终同一平面内。4.3模型的建立敌机位于(N,)。由于导弹头始终对准敌机，故此时直线 PQ 就是导弹的轨迹曲线弧 IP 在点 P 处的切线。yB:(0,M)IxO(0,0)图一：二维平面追机分析图 其中，下同即 (4.1) 弧的长度计算公式为：该公式的推导过程如下： 又根据题意，弧 IP 的长度为 |AQ| 的 倍，于是： (4.2)根据(4.1),( 4.2),消去，得： (4.3)满足条件： (4.4) 设，降阶得到一阶可分离变量方程： (4.5)其中 ，计算可得： (4.6)将代入得： (4.7)4.4模型求解与分析(1).现从特殊情况出发，以该导弹基地上空O(0,0)为空对空导弹发射点，即该微分方程定解条件为：，则可解得：(4.8)此方程仍满足临界条件：将代入(4.8)式得： 解得：，由于，为大于0的数，因此取：由上式可得以下结论：.若型导弹从基地出发，并在临界点S（N,M）击中敌机，则导弹的速度需为敌机速度的倍。 若型导弹从基地出发，要保证该导弹一定在敌机到达安全区的途中命中敌机，则该导弹的速度必须不小于敌机飞行速度得倍。(2).再推广到其他一般情况，即型导弹的发射地点为任意点I（a,b）,导弹命中敌机仍为临界状态下，即导弹最终命中点为点S（N,M），则再将定解条件代入（4.7）式得： (4.9)则可得： (4.10)(4.11)以上等式还满足：,将其代入(4.11)式可得：（4.12）由（4.12）式可以算出此种情况下的，记为。由此可得出以下结论：. 若型导弹从任意点I（a,b）发射，要使其在临界点S（N,M）击中敌机，则导弹的速度需为敌机速度的倍。. 若型导弹从任意点I（a,b）发射，要保证其必定在途中命中敌机，则导弹的速度必须不小于敌机速度的倍。根据（4.11）式，只要基地在初始状态下已知型导弹的发射地点I（a,b）、型导弹与敌机的速度比，就能计算得出导弹命中敌机所需的纵向距离（Y轴方向），只要满足，就能保证该型导弹一定命中敌机，以及命中所需时间和命中地点均可计算得出。5.问题2的解答5.1模型的分析本问题的难点在于三维空间的导弹追逐问题的求解。本模型通过数量关系，将三维空间合理转化为二维空间，简化求解。并需要假设敌机飞机的高度。之后便可以借助问题1的解答进行求解。5.2模型的假设(1).敌机飞行的高度为H(2).击中临界点依然为敌机到达安全区的位置M5.3模型的建立首先看三维空间：S:(N,M,H)yzQ:(N,H)A:(N,0,H)P:(x,y,z)B:(0,M,0)O:(0,0,0)xO图二：三维空间分析如图二：建立三维坐标系(x,y,z),原点为O:(0,0,0),击中的临界点为S: (N,M,H),敌机的初始点为A: (N,0,H),敌机任意时刻的坐标为Q:(N,Y,H)，导弹的坐标为P:(x,y,z),B点的坐标为(0,M,0)。敌机位于(N,Y ,H),由于导弹头始终对准敌机，导弹头的速度向量平行于敌机与导弹头位置的差向量，即(5.1)其中：，表示两者的比例关系。由(5.1)得， (5.2)经mathmatic求解(5.2)：得： (5.3)由此可得，导弹轨迹点中为一个常数，只要H确定，即敌机的飞行高度确定，就是一个已知的确定的常数。可进一步推得，导弹的轨迹必定在一个平面上，而不是在一个球面上。因为若导弹的轨迹在球面上，则不可能为一个常数。而此平面就为图2中的平面。导弹的轨迹始终在此平面上，而不会离开此平面。于是，三维空间可转化为二维平面。SyABSBOAxOyB:(0,M)xO(0,0)图三：二维平面的分析 如图三，由三维平面转化为二维平面后，点的坐标发生的相应的改变。敌机的初始位置为点，任意时刻t的坐标为点，导弹的坐标为点，击中的临界点为点。敌机位于，由于导弹头始终对准敌机，故此时直线 PQ 就是导弹的轨迹曲线弧 OP 在点 P 处的切线。即 (5.4) 弧的长度计算公式为：(5.5)该公式的推导过程如下： 又根据题意，弧 OP 的长度为 |AQ| 的 倍，于是：(5.6)根据(5.4),(5.6),消去，得： (5.7) 设，降阶得到一阶可分离变量方程： (5.8)其中 ，定解条件为：，则计算可得： (5.9)将代入得：(5.10)5.4模型求解与分析(1).微分方程(5.9)的定解条件为：，则可解得：(5.11)此方程仍满足临界条件：将代入(4.11)式得： 解得：，由于，为大于0的数，因此取：由上式可得以下结论：.若II型导弹在临界点击中敌机，则导弹的速度需为敌机速度的倍。.若要保证II型导弹一定在敌机到达安全区的途中命中敌机，则该导弹的速度必须不小于敌机飞行速度得倍。(2).再推广到其他一般情况，即II型导弹的最终命中点为点，则：将定解条件代入（5.11）式得： 解得：，由于，为大于0的数，因此取：由上可得出以下一般结论：从基地发射的型导弹，若要保证其在途中一定命中敌机，则该导弹的速度与敌机的速度之比需满足：其中：6.问题3的解答6.1模型一的数值求解6.1.1运用欧拉公式求解由于t时刻导弹的位置为，其速度可以由水平分速度和垂直分速度构成：（6.1）导弹的方向指向敌机，所以即 （6.2）可以由此推出一阶微分方程组：（6.3）采用欧拉法进行数值求解得：其中为步长（6.4）计算到为止。我们不妨设敌机的速度，导弹的速度，巡航飞机的初始位置为（0，0），敌机的初始位置为（100，0），取步长用 matlab编程如下：clear;clc;l=0.01;%时间步长k=1;N=100;v1=680;v0=340;t(1)=0;y(1)=0;x(1)=0;%初始值while x=100y(k+1)=y(k)+v1*l/sqrt(1+(N-x(k)/(v0*t(k)-y(k)2);x(k+1)=x(k)+v1*l/sqrt(1+(v0*t(k)-y(k)/(N-x(k)2);t(k+1)= l*k;k=k+1;endplot(x,y,x(1):0.05:x(end),100);title(导弹追击敌机图);grid on;t=t(end)x=x(end)y=y(end)图四：二维导弹追机图（1）由欧拉公式解得=0.1900，=100.0183，=61.3351，即型导弹将在（100.0183，61.3351）处命中敌机，所需时间为0.19。将a=0，b=0，将，代入式（4.12）得到，与数据求解所得结果61.3351较为接近。如果取步长，所得结果为66.4994，与模型值高度接近。命中过程中个时间点导弹位置数据表如下所示：（）表(一)t00.010.020.030.040.050.06x06.813.595520.37627.128733.837840.4837y000.24790.76211.56252.67074.1101t0.070.080.090.10.110.120.13x47.042553.484759.774265.867171.710677.241382.3849y5.90538.082110.667113.686317.16421.120225.568t0.140.150.160.170.180.19x87.056591.161994.60297.27999.1063100.0183y30.509335.930141.795848.046654.596561.3351对给定速度的条件下，对于不同的M，N值，只要改变程序里的M，N的初始值即可。6.1.2改进欧拉公式进行数值求解当时，导弹位于（6.5）其中为步长，为导弹飞行倾角推导过程如下：时，敌机的位置为,导弹在原点指向。时，敌机的位置为，导弹得位置为，而，导弹得飞行倾角为：。时，敌机得位置为，导弹得位置为， 由得表达式可以写出，得表达式，此时导弹得飞行倾角为以此类推.当时，导弹位于（6.6）Matlab程序：clear;clc;l=0.01;%时间步长k=1;N=100;v1=680;v0=340;t(1)=0;y(1)=0;x(1)=0;%初始值while x=100x(k+1)=x(k)+v1*l*(N-x(k)/sqrt(k*v0*l-y(k)2+(N-x(k)2);y(k+1)=y(k)+v1*l*(k*v0*l-y(k)/sqrt(k*v0*l-y(k)2+(N-x(k)2);t(k+1)= l*k;k=k+1;endplot(x,y,x(1):0.05:x(end),100);title(导弹追击敌机图);grid on;t=t(end)x=x(end)y=y(end)图五：二维导弹追机图（2）由改进欧拉公式解得=0.2000,=100.0924,=70.4072，即导弹命中敌机的地点为（100.0924，70.4072），命中所需时间为0.2000将a=0，b=0，将，代入式（4.12）得到，与数据求解所得结果70.4072较为接近。如果取步长，所得结果为66.7260，与模型值高度接近。命中过程中个时间点导弹位置数据表如下所示：（）表（二）t00.010.020.030.040.050.06x06.796113.579220.338627.060933.7340.3263y00.23110.70911.45152.47663.80455.4564t0.070.080.090.10.110.120.13x46.826153.200859.415865.429671.193176.648481.7286y7.45449.821612.58115.754919.363523.42327.9431t0.140.150.160.170.180.190.2x86.358290.454593.930696.701298.688699.8316100.0924y32.923838.351644.195950.405956.90963.612270.40726.2模型二的数值求解 由于t时刻导弹的位置为，其速度可以由方向分速度、方向分速度和方向分速度构成：（6.7）导弹的方向始终指向敌机，所以有：（6.8）其中r 为系数可以由此推出一阶微分方程组：（6.9）采用欧拉法进行数值求解得：（6.10）再不妨设敌机飞行高度H=500Matlab程序：clear;clc;l=0.01;%时间步长k=1;N=100;v0=340;v2=680;H=500;t(1)=0;y(1)=0;x(1)=0;z(1)=0;%初始值while x=100y(k+1)=y(k)+v2*l/sqrt(1+(1+H2/N2)*(N-x(k)/(v0*t(k)-y(k)2);x(k+1)=x(k)+v2*l/sqrt(1+H2/N2+(v0*t(k)-y(k)/(N-x(k)2);z(k+1)=z(k)+(H/N)*(x(k+1)-x(k);t(k+1)=l*k;k=k+1;endplot3(x,y,z,r);title(三维情况下导弹追击敌机图);grid on;t=t(end)x=x(end)y=y(end)z=z(end) 图六：三维导弹追机图用欧拉公式计算得到：=1.0000，=100.0039，=338.9198，=500.0194，即导弹命中敌机的地点位置为（100.0039，338.9198，500.0194），所需时间为1.0000 将，代入式（5.11）得到 数值计算所得结果338.9198很接近。命中过程中个时间点导弹位置数据表如下所示：（）表（三）t00.050.10.150.20.250.3x06.667913.331619.981326.60533.187839.7119y001.12383.72937.642213.055820.0809z033.339766.658299.9066133.0248165.9392198.5595t0.350.40.450.50.550.60.65x46.15552.490158.684264.697270.480975.977881.1201y28.835839.44452.032466.72783.6464102.8922124.5361z230.775262.4505293.4209323.486352.4045379.8888405.6006t0.70.750.80.850.90.951x85.830390.021693.601596.476698.56199.7849100.1043y148.602175.0455203.73234.4069266.703300.1254344.0863z429.1515450.1081468.0073482.3831492.8049498.9244500.5213对给定速度的条件下，对于不同的M，N，H值，只要改变程序里的M，N，H的初始值即可7.模型的优缺点与改进方向7.1优点：(1)使用经典的对微积分方式对其求解的思想，也即先降阶然后再通过转化的方式，得到方程的解，然后采用差分方程，用欧拉法，利用数学实验工具Matlab软件实现求解的工程，比较简单，与得到的方程的解进行对比，验证了模型的正确性。(2)使用所建立模型可以在忽略从发现到发射时间误差的前提下（即如上述模型假设情况下）比较准确的估计导弹击中敌机时刻和以及击中的位置。7.2缺点：(1)未考虑导弹的射程。(2)实际情况下，导弹要受空气阻力，速度是变化的，而且敌机的速度也并不一定是恒定的，导弹和飞机也不可能一直在同一平面上。 7.3模型实践性改进：本文还可以求出导弹射中敌机的概率。(1)当 t无限趋于0的时候，导弹击中敌机的概率很大，趋近于1(2)当Y趋向于M ，即t无限趋向于 的时候，追踪导弹的概率趋向于0(3)函数是时间t的减函数。拟建立线性概率函数：其中 其中 t=0表示的是追踪导弹近乎在敌机没开始逃窜的情况下将其击落。 表示的是在无限趋于敌机安全区时追到敌机将其击毁条件下导弹的飞行路程。当然，还可以对线性概率函数做出合理的修正。另外，本文同样适用于航母与护航舰搜寻人质问题，潜艇决策问题。本题还可还以引申求解两种情况下导弹在安全区内击毁敌机的条件。本模型在实际生活生产与军事上有较好的实际意义。 8.参考文献1 刘慧颖 MTATLAB R2006a基础教程M. 北京：清华大学出版社，20072 徐全智，杨晋浩 数学建模（第二版）M 北京： 高等教育出版社，20083 李天胜，徐军 数学模型应用实例M 合肥： 合肥工业大学出版社，20074 谢进， 李大美 MATLAB与计算方法实验M 武汉：武汉大学出版社 ， 20095 朱道元 数学建模案例精选M. 北京:：科学出版社,2003.6 张天笑，杨奋强.MATLAB7.x基础教程M. 西安：西安电子科技大学出版社，2008. 7 赵静，但琦 数学建模与数学实验M 北京：高等教育出版社 20088 豆丁，导弹追踪飞机问题的初步探索T 2009年8月11日9 豆丁，导弹追踪（07.3T 2009年8月11日10 吴建国 数学建模案例精编M 北京： 中国水利水电出版社 20059.附录 （略）19