多基雷达测飞行目标及攻击模型.doc

多基雷达测飞行目标及攻击模型
摘要
本题为典型的微分方程模型。
问题一，通过分析得出多基雷达定位方程组，讨论当雷达数目分别为1，2，3时的各种不同情况，得出至少需要三个雷达才能定位飞行物。
问题二，利用控制变量法分别讨论坐标误差和距离误差对定位精度的影响，可知当雷达间距离很近，雷达距目标较远以及雷达间距离比较远，同时和目标的距离又比较近时，误差比较大。
问题三，将问题一中所得三元方程组简化为二元方程组，两次利用最小二乘法得出飞行物的坐标分别为甲（-4.23，1317.26，10997.6），乙（50844.97，-5.55，10883.2），丙（44.17，106865.6，14690.13）。由计算过程可知，尽量使同一组雷达的排列不再同一条直线上，以及最大限度的利用雷达所测数据，都是提高雷达定位精度的方法之一。
问题四，由题意得出Ι型导弹追踪乙机轨迹图，由此得出微分方程，即导弹追踪敌机轨迹方程。当导弹击中乙机时，有,此时，则导弹飞行距离为。因此导弹击中乙机的条件为。由数据得出敌机被击中的时刻为98.44s之后，位置为（50844.97，33464.55，10883）。
问题五，经分析，导弹轨迹位于斜平面内，通过坐标转换，得出ΙΙ型地对空导弹追踪敌机的轨迹，经计算，得出在时乙机被击毁，其被击毁的条件为。
问题六，计算机编程的算法及相关程序见附录。经计算，得出敌机在197.22s之后于（50844.97，67054.78，10883.4）处被击中，并给出了ΙΙ型导弹的轨迹图。
问题七，当敌机不变速时，由题意得出ΙΙΙ型导弹轨迹方程，从而得出丙机在126.53s之后于（28119.3，106865.6，14690.1）处被击中。当敌机做俯冲逃离时，利用问题六中所建斜平面建立微分方程求得敌机和导弹间的距离，再利用坐标转换，可求得敌机在126.21s之后于(28167.96,104234.91,10347.50)处被击中。

关键词：微分方程模型 最小二乘法 图文结合 坐标变换
1 问题重述
.在我防空指挥部的上空发现有一可疑的飞行物，需要对其进行精确定位。常用的定位方法是基于多基雷达的测量方法。每个雷达都可以测量自身的坐标以及它到飞行物距离，其中为雷达的总数。但存在测量误差，给精确定位带来了困难。通过一组雷达位置坐标和飞行物到各雷达的距离测量，我们可以确定目标的空间飞行物的坐标。
现在，我防空指挥部多部雷达发现有一架来路不明的飞机乙（第二组数据），经分析确认是一架敌机后，即命令正处于同一高度进行巡逻的我方飞机甲（第一组数据）发射I型空对空追踪导弹将其击毁（追踪导弹可针对目标随时自动调节追踪方向）。假定雷达发现敌机时，该机正位于我防空指挥部正东N公里高空处，并欲在同一高度上向位于其正北方向M公里处的安全区逃窜（由于电子干扰的作用敌机一旦进入安全区后.导弹将失去追踪目标，无法将其击毁）。（如果没有明确给出数值，可假定敌机速度v = 1 马赫数 ， 导弹u = 2马赫数）。
设距离误差服从正态分布，坐标误差服从正态分布。在这个假定下完成以下工作。
一、至少需要几个雷达才能定位飞行物？
二、在最少雷达的条件下，分析并比较距离误差和坐标误差对定位精度影响。
三、在实际情况中，往往使用更多雷达进行精确定位，试设计一种定位算法。对三组雷达得到的测量数据，计算飞行物的坐标。并给出控制雷达定位精度的建议。（数据见附录）
四、在适当的假设下，确定导弹追踪敌机的轨迹及发射I型空对空导弹击毁敌机的条件。并且利用给出数据，计算出敌机被击中的时刻以及当时敌机被击毁的位置。
五、若当时命令设在防空指挥部的地面导弹基地发射II型地对空追踪导弹截击敌机，假定敌机始终距地面高度为h公里飞行，其他假定同上一问中所述，重新确定此时II型地对空导弹追踪敌机的轨迹及击毁敌机的条件。并且利用给出数据，计算出敌机被击中的时刻以及当时敌机被击毁的位置。
六、若敌机的飞行速度 v 、其位置 N 和追踪导弹速度 u 均为给定的常数；针对情况1中敌机被导弹击中的条件下，给出一个计算机编程的算法及相关程序，以计算出敌机被击中的时刻以及当时敌机的位置。 据此，在 v = 1 马赫数 ， N = 100公里 , u = 2马赫数时，利用上述的程序算出具体敌机被击中的时刻以及当时敌机被击毁的位置。
七、若在发现飞机乙的同时接指挥部通知：又发现一架敌方的飞机丙（第三组数据）以速度800km/h向东水平飞行。该飞机甲立即发射III型空对空导弹，导弹速度3200km/h，方向指向敌机。
（1）若敌机不变速，则何时何地击中？
（2）当导弹离敌机约25km时，敌机发现遭受攻击，立即以30度俯角俯冲，速度1000km/h，方向与III型空对空导弹垂直，何时何地击中。
2 问题分析
问题一，在不考虑各种误差的情况下，由多基雷达定位目标方程组出发，我们讨论了当雷达数目为1，2，3时的不同情况，通过讨论得出了至少需要三个雷达才能定位飞行物的结论。
问题二，由题意列出相关表达式，利用控制变量法，分别讨论距离误差和测量误差对定位精度的影响。
问题三，由三元多基雷达定位方程组，可以得到二元一次超静定方程组，利用最小二乘法求出其近似最优解，最后再反带回多基雷达定位方程组，再次利用最小二乘法得出值，至此，得出飞行物的坐标。由分析过程可以总结出控制雷达定位精度的方法。
问题四，根据题意得出导弹追踪乙机轨迹图，由图及约束条件出发，得出导弹追踪乙机轨迹模型，由导弹发射瞬间的函数关系得出导弹击毁敌机的条件，由问题三所得数据算出参数，并计算出敌机被击中时刻及被击毁位置。
问题五，由题意可知，导弹轨迹位于空间斜平面内，通过坐标转换，得出导弹追踪敌机轨迹方程，经过求解得出敌机被击中时刻及被击中位置。
问题六，由题意得出计算机编程的算法及相关程序，代入数据得出敌机被击中的时刻及当时敌机的位置。
问题七，当敌机不变速时，经坐标转换可以得到ΙΙΙ型导弹轨迹方程，从而计算出敌机被击中时刻及位置。当敌机做俯冲逃离时，利用问题六中所建斜平面建立微分方程，求得飞机和导弹间的距离时。再利用坐标转换，即可求得敌机被击中的时刻及位置。
3 模型假设
1 距离误差服从正态分布，坐标误差服从正态分布；
2 敌机始终在距离地面高度为H公里出飞行；
3 导弹和敌机的速度一经确定即为匀速飞行；
4 不考虑天气等自然因素的影响；

5 假设每次导弹发射时刻均为0。
4 符号说明
： 定位需要雷达数目；
： 雷达到目标的测量距离；
： 雷达的自身测量坐标；
： 甲飞机初始时刻位置；
： 乙飞机初始时刻位置;
： 丙飞机初始时刻位置。
5 模型的建立与求解
5.1 问题一，最少雷达数目求解
在不考虑各种误差的情况下，本问题多基雷达定位目标的方程如下所示;

其中代表雷达数目。显然当时雷达只是定位目标在距离雷达自身半径为的球面上，因此无法定位。为了阐述方便我们定义这样的球面为定位球面。
当时，由于目标的客观存在，两个定位球面必然相交。目标就在该相交线上。当且仅当目标位于两雷达之间时两定位球面才可相交于一点（下图一），否则相交部分为一圆（下图二）。如下图所示：

当时，存在两种情况。三个雷达共线时其定位球面仍然相交于一个圆（下图三），当三个不共线时三个定位球面必然相交于一个点，此点正是目标所在点（下图四）。

综上，至少需要三个不共线的雷达才能定位目标。
5.2 问题二的建立与求解
5.2.1 误差模型的建立
在有3个雷达测量定位情况下分析并比较距离误差和坐标误差对精度的影响，记3个雷达的测量值为，距离为，。定位目标飞行物的坐标为。
则有

从而可解得

则目标飞行物的坐标（x，y，z）的误差可表示如下：

下面分析各因素的影响。
（1）分析坐标误差对定位精度的影响
此时可作如下假设：
（a）所测得的距离是精准的，即
（b）所测得的是精准的，即且有.
（可类似讨论对精度的影响，下文将不对作进一步讨论分析）
不妨令，由上式解出值为：

考虑到的平等性，可令，从而由的变化对的影响来讨论坐标误差对定位精度的影响，则由

由题意：

则

同理可讨论对精度的影响，此处省略不写。
（2）分析距离误差变化对精度的影响
此时可作如下假设：
所测得的是精准的，即且有，.
由的平等性，可令，同样不妨令，由上式解出值

则

由题意则
,

同理可讨论对精度的影响
5.2.2误差分析
根据a式和b式可知：
1.当雷达间距离很近时，误差比较大；
2.当雷达距目标太远时，误差也比较大；
3.当雷达间距离比较远，同时和目标的距离又比较相近时，误差比较小。、
5.3问题三，飞机空间位置的求解
根据所得方程组（1），用方程组（1）的第一个方程分别与后面方程相减。得到以下二元一次方程组（2）：

令

则以上二元一次方程组退化为：

对于这样的超静定方程没有精确解，我们利用最小二乘法求解其近似最优解。令

方程组可以表示为，其近似最优解为；但是这样只能求得，无法求得.在求解时我们将求得的代入方程组（1）的每一个方程中，分别得到不同的，然后求其平均值。
至此我们可以求得飞机甲，乙，丙的空间位置。如下表所示：
表一 飞机空间坐标结果
飞机编号
空间坐标
x
y
z
甲
-4.23
1317.26
10997.6
乙
50844.97
-5.55
10883.2
丙
44.17
106865.6
14690.13
由上述建模过程可以得知，尽量使同一组雷达的排列不再同一条直线上，以及最大限度的利用雷达所测数据，例如用最小二乘法，都可以提高雷达定位的精度。
5.4 问题四，I型导弹轨迹模型的建立和求解
5.4.1模型的建立
导弹追击乙机轨迹图如下：

图五
由于导弹时刻对准飞机，则有：

由二者的轨迹长度关系，我们得到：

由以上两个方程联立求解，得导弹轨迹的模型如下所示：

5.4.2 模型的求解
由以上两个式子消去得到：

令，等式两边同时对求导：

令，则以上方程降为如下一阶微分方程，分离变量得：

同时积分得到：

（1）
将代入（1）式，进一步求得其解析解如下：
（2）
导弹发射瞬间有，将初始值代入（1）式，则有：

将导弹发射时刻位置坐标代入（2）式得到：

导弹击中乙飞机时，有,此时，导弹飞行距离为。则能够击中乙飞机的条件为。
将问题三所得数据代入方程得到。此时空间坐标为（50844.97，33464.55，10883），假设导弹发射时刻时间，此时。
5.5 问题五，II导弹轨迹模型的建立与求解
图中给出II型导弹击毁乙机的大致轨迹：

图六 II型飞机轨迹示意图
5.5.1 模型的建立
由图我们可以看出导弹轨迹处于平面内，所以我们建立在该平面内求解轨迹方程的模型后再将坐标转换回原坐标系，即可得到空间内导弹轨迹方程。具体做法如下：
我们建立如图新坐标系，模型同问题四相似，如下所示：

5.5.2 模型的求解
以上模型的求解同问题四中求解方法相同，结果如下所示：
(3)
其中

由坐标转换公式我们得到以下坐标转换式：

代入（3）式得到原坐标系中的方程：

新平面坐标在原空间坐标系中的方程为：

将两式联立即可得到II型导弹轨迹方程：
乙机被击中的时刻为：，则只有满足条件时，才能击毁乙飞机。
5.6 问题六的求解
计算机编程的算法及相关程序【见附录】。
当在马赫数，公里，马赫数时，运用Matlab可解得敌机被击中的位置为。
击中时刻
II型导弹实际轨迹图如下图七：

图七 II型导弹轨迹图
5.7 问题七的求解
5.7.1敌机不变速情况下被击中时间和地点的求解。
由问题三所得数据，我们能够判断出甲丙飞机的空间相对位置，从而很容易作出下图六：

图六 III型导弹击中丙飞机轨迹示意图
由图我们可以看出导弹轨迹处于平面内，所以求解本问题和问题四相同。
在新坐标下的轨迹方程为：
（4）
其中 （5）

将起始时刻数据代入（4）式得到：

由坐标变换公式得到以下坐标关系式：

将其代入（4）式得到在原坐标系下的方程为：

新平面坐标系在原空间坐标系的方程为：

将两式联立即可得到IIΙ型导弹轨迹方程：

将初始数据代入方程得到，丙飞机被击毁时其空间坐标为：
击中时刻为：
III型导弹实际轨迹图如下图八：

图八 III型导弹实际轨迹图
5.7.2 在距离敌机时敌机做俯冲逃离情况下被击中时间地点的求解
该问题我们分为丙机俯冲前后两个部分来求解。
（1）俯冲前我们在上一问题中建立的斜平面内建立如下微分方程：

利用程序【见附录】求解该方程数值解，并且搜索得到飞机和导弹间的距离为时的坐标以及时间。
利用上一问题已经得到的坐标转换公式即可得到在原坐标系中的作标。
（2）

图九 III型导弹追击敌机示意图
如图我们以为坐标原点，导弹追击轨迹和飞机逃离轨迹构成的平面建立坐标系。
设时刻丙飞机坐标为则坐标轴方向在空间内单位向量：

设轴与空间坐标内平面相较于点。则坐标轴方向在空间内单位向量：

由几何关系我们得到以下方程求解。

在该平面内导弹轨迹以及击中时间地点的求解和问题六的求解相同。设击中时刻丙飞机坐标。
求得丙飞机击中时在坐标：。击中时刻。
由坐标转换公式我们得到该坐标系与原坐标系坐标转换式：

代入数据得到击中坐标
时间
6 模型的推广与改进
在本题中，我们只根据题目要求给出了多基雷达测飞行目标及攻击模型，由于时间限制，并未对雷达的排列方式进行深入的讨论。若再加上雷达排列方式的模型限制，问题的求解可能会更加精确。
7 模型的评价
7.1 优点
1 问题一中确定定位飞行物的最少雷达数目时，我们通过讨论的方式论证了1个，2个及在同一条直线上的3个雷达是不满足要求的。
2 在本题中我们利用最小二乘法处理数据，使得所给数据得到了最大限度的利用，提高了求解的精确度。
3 在求解过程中，我们通过图文结合的方式给出了解答，通俗易懂；
4 在计算过程中，我们将空间曲线转化为平面曲线求解，使得计算过程更加简便易行。
7.2 缺点
在处理数据时，我们采用了最小二乘法，使数据得到了最大限度的利用，但并没有考虑数据中的坏值。
参考文献
1 姜启源等.数学模型（第三版）[M].北京：高等教育出版社，2003
2 杨建军.地空导弹武器系统概论.北京：国防工业出版社，2006
3 吴祈宗.运筹学[M].北京：机械工业出版社，2002

附录：
程序：
程序一：
%I型导弹击中乙机轨迹程序
u=680;v=340;
m=v/u;
x0=-4.3;y0=1317.1;
xe0=50844.97;ye0=-5.55;
p0=(ye0-y0)/(xe0-x0);
C1=(p0+sqrt(1+p0^2))/(xe0-x0)^(-m);
fprintf('C1=%f\n',C1)
C2=y0-(xe0-x0)^(1+m)/(2*C1*(m+1))+C1*(xe0-x0)^(1-m)/(2*(1-m));
fprintf('C2=%f\n',C2)
break
ezplot('(50844.97-x)^1.5/659-219.7086*(50844.97-x)^0.5+33464.545961',[-4.3,51844.97])
hold on
y=-5.5:33464.545961
x=50844.97;
plot(x,y,'r')
title('I型导弹击中乙机轨迹')

程序二：
%II型导弹击中敌机轨迹方程
N=100000;
ye0=-5.55;
ze0=10883.3832;
xe0=50844.97;%原坐标里的敌机初始x坐标值
xe=sqrt(N^2+ze0^2);%转换后起始坐标值
fprintf('xe=%f\n',xe)
m=0.5;
c1=(N^2+ze0^2)^(1/4);
fprintf('c1=%f\n',c1)
c2=c1*xe^0.5-xe^1.5/(3*c1);
fprintf('(c2+ye0)=%f\n',c2+ye0)
a=cos(atan(ze0/N));
fprintf('a=%f\n',a);
x=-49156:100:50850;
for i=1:1001
y(i)=(xe-(x(i)+N-xe0)/a)^1.5/(3*c1)-c1*sqrt(xe-(x(i)+N-xe0)/a)+c2+ye0;
%轨迹方程
z(i)=ze0*x(i)/N+ze0-ze0*xe0/N;
end
plot3(x,y,z,'r',50844.97*ones(1,1001),y,10883.3832*ones(1,1001),'g'),grid
t=(c2+ye0)/340

程序三：
%III型导弹击中敌机的轨迹方程
vb=800; %'b'表示丙
vd=3200; %'d'表示导弹
m=vb/vd
xd0=44.23; %导弹起始坐标
yd0=1317.26;
zd0=10997.18;
xb0=44.17; %丙飞机的起始坐标
yb0=106865.6;
zb0=14690.13;
a=cos(atan((zb0-zd0)/(yb0-yd0)));
y1b=sqrt((zb0-zd0)^2+(yb0-yd0)^2);%转换后起始坐标值
c1=y1b^m;
c2=c1*y1b^(1-m)/(2*(1-m))-y1b^(1+m)/(2*c1*(1+m));
fprintf('c2=%f\n',c2)
y=1317:100:106870;
for i=1:1056
x(i)=(y1b-(y(i)-yd0)/a)^(1+m)/(2*c1*(1+m))-c1*(y1b-(y(i)-yd0)/a)^(1-m)/(2*(1-m))+c2+xd0; %轨迹方程
z(i)=(zb0-zd0)*y(i)/(yb0-yd0)+zd0-(zb0-zd0)*yd0/(yb0-yd0);
end
plot3(x,y,z,'r',x,106865.6*ones(1,1056),14690.13*ones(1,1056),'g'),grid

程序四：
%用来求解导弹距离丙机还有25km时的程序
%ff3主程序文件
function Dy=ff3(t,y)
Dy=zeros(2,1);
Dy(1)=3200000/3600*(800000/3600*t-y(1))/sqrt((800000/3600*t-y(1))^2+(105612.93-y(2))^2);
Dy(2)=3200000/3600*(105612.93-y(2))/sqrt((800000/3600*t-y(1))^2+(105612.93-y(2))^2);

[t,y]=ode45('ff3',[0 125],[0 0]);
for i=1:93
r(i)=sqrt((800000/3600*t(i,1)-y(i,1))^2+(105612.93-y(i,2))^2);
if r(i)<25000
xb=800000/3600*t(i,1);
b=t(i,1);
c=y(i,1);
d=y(i,2);
break
end
end
xd0=44.23; %导弹起始坐标
yd0=1317.26;
zd0=10997.18;
xb0=44.17; %丙飞机的起始坐标
yb0=106865.6;
zb0=14690.13;
A=(zb0-zd0)/(yb0-yd0);
a=cos(atan(A));
x=c+xd0;%原坐标系中导弹距离飞机25km时的坐标
y=d*a+yd0;
z=A*y+zd0-A*yd0;
fprintf('b=%f\nxb=%f\nyb0=%f\nzb0=%f\nx=%f\ny=%f\nz=%f\n',b,xb,yb0,zb0,x,y,z);
%输出导弹距离丙飞机25KM的时刻b,以及在此时导弹和飞机在原坐标系中的坐标
22