1、飞行器空间的坐标修正摘要飞行器的导航精度问题一直是航空航天领域研究的重要课题，惯性导航系统是一种不依赖于任何外部信息的自主式导航系统，在航空航天领域起着越来越重要的作用。由于其系统结构误差、惯性测量部件误差、标度系数误差等因素的影响，惯性导航系统的积累误差随着时间的推移而逐渐增大，这一问题严重影响到航空航天技术的发展。对于问题一，先对坐标数据进行分析，找出误差分布，利用 Matlab 进行拟合分析，发现 X、Y 坐标的变化是线性的，完全符合x=-65*t+1360.3; y=-70*t+1309.5;随机误差波动较小，而 Z 坐标的波动较为明显，对此，分别采用均值迭代法和极小似然迭代法两种方法
2、进行修正，使 Z 坐标的观测值更逼近真实值，有效消除了随机波动误差。并对两种方法的模型进行比较，得出在消除随机误差过程中极小似然迭代法优于均值迭代法的结论。对于问题二，我们用问题一修正过后的数据去修正飞行器飞行速度所得的坐标，对于修正后的误差用改进后的极小似然迭代法拟合出差值多项式，避免了惯性导航系统中误差随着时间的推移而增大的效应。对于问题三，只要给出某飞行器的具体数据，都可以用问题一及问题二的方法来减小误差，由于时间问题，下面只给出对于给出数据的误差修正的过程，而没有用具体数据来检验模型的效果。关键词：均值迭代法，线性拟合，极小似然迭代法，插值多项式一、问题重述1.1 问题背景飞行器的导航
3、精度问题一直是航空航天领域研究的重要课题，惯性导航系统是一种不依赖于任何外部信息的自主式导航系统，在航空航天领域起着越来越重要的作用。由于其系统结构误差、惯性测量部件误差、标度系数误差等因素的影响，惯性导航系统的积累误差随着时间的推移而逐渐增大，这一问题严重影响到航空航天技术的发展。目前关于定位精度的研究成果主要是从物理技术(例如红外测距)方面来提高定位的精度，近年来，围绕定位坐标精度问题的相关研究也渐渐展开。因此进一步研究飞行器空间坐标修正方法有重要的理论意义和应用价值。本题的目标是利用数学的方法对飞行器的误差进行修正，并利用结果进行飞行器的仿真。 1.2 问题提出 某一观测站测得飞行器空间
4、位置(假设观测站为坐标原点)X(x、y、z)，飞行器的飞行速度 V(x 轴、y 轴、z 轴),飞行器与观测站之间的偏向角 ，俯仰角 以及观测数据的时间间隔 t。所给的各项数据均含有一定的误差，其中观测站的坐标 (0,0,0)不含误差，飞行器的坐标（观测值）可能含有较大误差。请根据所给数据进行如下工作： 问题一：飞行器坐标的数据为观测值，由于电子仪器的精度和噪声干扰等，含有一定的误差波动，建立数学模型对飞行器坐标观测值的随机波动误差进行修正。问题二：由于观测数据的仪器误差，飞行器坐标在长时间的飞行中，坐标数据的观测值由于误差的累积发生漂移，建立数学模型，对飞行器的坐标的这种误差进行修正。(提示：
5、在短时间内，可以视为飞行器坐标含有一定的常量误差，或者飞行器的这种误差是线性变化的)。 问题三：结合具体的飞行器给出误差修正方案。二、模型假设假设一：不考虑地球自转等因素；假设二：每一时间段内飞行器匀速运动，如 0-0.15s 内速度不变；假设三：在短时间内，可以视为飞行器坐标含有一定的常量误差，或者飞行器的这种误差是线性变化的；假设四：观测站为坐标原点（0，0，0）不含误差三：符号说明与解释Xi 观测值 x 方向上坐标alpha 飞行器与观测站之间的偏向角 .Yi 观测值 y 方向上坐标theta 飞行器与观测站之间的俯仰角 Zi 观测值 z 方向上坐标Vx X 方向上速度Vy Y 方向上速
6、度 VZi Z 方向上速度T 时间间隔四、问题分析问题一:数据的观测值受到电子仪器的精度和噪声干扰等，电子仪器的精度产生的误差无法避免，模型一通过均值迭代法尽量减小噪声干扰带来的误差；模型二根据实测数据首先预测误差然后通过迭代方法逐步减小误差，使得实测数据之间达到最大程度的一致，此模型是基于实测数据而创造的。根据 x，y 坐标与时间 t 的数据，通过 Matlab 做函数拟合：x=-65*t+1360.3;x方 向 坐 标 X1数 据 图 像020040060080010001200140016000 0.3 0.6 0.9 1.2 1.5 1.8 2.1 2.4 2.7 3 3.3 3.6
7、3.9 4.2 4.5 4.8 5.1 5.4 5.7 6时 间 t（ s）x方向坐标x方 向 坐 标 数 据 图 像y=-70*t+1309.5;由于拟合函数值并与实测数据完全相同，可知 x，y 方向坐标与时间都是完全线性的，可以认为没有误差。高度的数据存在很大波动，受噪音干扰的影响较大，所以下面准备建立均值迭代法和极小似然迭代法两个模型来消除噪音干扰带来的影响。y方 向 坐 标 Y1数 据 图 像02004006008001000120014000 0.45 0.9 1.35 1.8 2.25 2.7 3.15 3.6 4.05 4.5 4.95 5.4 5.85时 间 t（ s）y方向坐
8、标y方 向 坐 标 数 据 图 像H1高 度 数 据 图 像509.7509.8509.9510510.1510.2510.3510.4510.5510.60 0.45 0.9 1.35 1.8 2.25 2.7 3.15 3.6 4.05 4.5 4.95 5.4 5.85时 间 t（ s）z方向坐标高 度 数 据 图 像问题二：首先从原始数据中发掘新数据，即根据源数据中飞行器高度与垂直速度的关系估测出数据中存在的固有误差，由此固有误差，修正高度和速度；在由修正后的高度和速度值估测新的误差，如此循环往复直至 41 个随机误差尽量达到最小并且随机误差的均方差亦为最小，已达到既修正原不准确数据又
9、保证数据的一致性的效果。这就是接下来模型中将要介绍的极小似然法。问题三:问题要求给出具体飞行器，这里将导弹作为例子来给出误差方案，只要根据问题一问题二的方法来得出结论。五、模型的建立与求解5.1 针对问题一：5.1.1 模型建立与求解模型一:均值迭代法采用均值迭代原理，由上面分析，x、y 方向上对空间坐标影响不大，不需要修正，所以关键是对 h 修正，以下是对 h 采用迭代法后的图像。509.7509.8509.9510510.1510.2510.3510.4510.5510.60 0.3 0.6 0.9 1.2 1.5 1.8 2.1 2.4 2.7 3 3.3 3.6 3.9 4.2 4.5
10、 4.8 5.1 5.4 5.7 6观 测 值迭 代 一 次迭 代 两 次迭 代 三 次由上图可知，h 的数据修正后得到了很好的误差修正，这些数据的方差几近于零，具体的迭代数据见附录五。模型二: 极小似然迭代法由以上分析,小 x,y 坐标基本不存在误差,于是此处仅讨论飞行器高度误差;而所需用到的数据仅有 z 轴(垂直方向)坐标和速度,寻求两组数据的吻合性并尝试使两组数据达到最大程度的一致性,即克服由于观测仪器误差带来的数据不一致性.极小似然的基本思路如下:1.通过高度数据估测出一组垂直方向速度:数据组 1；211,.nVZZ2.通过取垂直方向相邻速度的平均值得出另一垂直方向速度数据组 2；12
11、1:()/,.()/2nV由于两组数据都表示 0-0.15，0.15-0.3，.各时间内的平均速度，因此可以认为两组数据的差作为飞行器垂直方向上的速度误差。因此，第三步如下：3.数据组 1 减数据组 2 作为速度误差：delta=V-eltaV120.5;Zltd速 度 误 差 ： ； 此 处 ， ， 均 为 向 量进 而 求 出 高 度 误 差 ：4.根据上一步得到的 Vdelta 计算飞行器高度值 Z（1） ：由图 3，高度值是一个上下波动的不稳定的值且在均值上下浮动，因而可以用高度的平均值近似等于第一个高度值，其余高度值可加上 Zdelta 来修正：(1)()2(1)(1)sumz0/l
12、ength();. (1);nZdaelt5.修正各时刻速度使之尽可能与高度数据一致：0 时刻的速度值取法与高度值类似,最后时刻的速度可通过加上速度误差来近似 (1)(1)()23(1)()(1)2=sumv/length0; /.3;=+Vdeltan-nnZ6.求向量 Vdelta,Zdelta 各分量的标准差，即采用最小二乘法的思想求 41个高度误差和速度误差的标准差，相应的 Matlab 程序见附录一，如果标准差小于给定数值精度 epsonV 和 epsonZ，则表示执行成功，迭代终止；7.重复以上 6 个步骤，得到(2)3()(2)(),.,;,.;nnZVdeltaelt直到满足条
13、件：标准差小于给定数值精度 epsonV 和 epsonZ，迭代终止；模型二的整个计算过程由 Matlab 完成，编写的相应程序请见附录二。5.1.2 问题求解结果x,y 坐标与时间 t 的关系代入原始数据：高度值 z0,速度值 v0,通过 Matlab 运行程序：z,v,Zdelta,Vdelta,count,Dv,Dz,flag=fun1(z0,v0,1,1,150)fun1 的输入参数意义分别为：原高度数据，原速度数据，迭代精度 1，迭代精度 2，最大迭代次数；输出参数分别为：修正后高度数据，修正后速度数据，高度误差，速度误差，迭代次数，高度误差的标准差，速度误差的标准差，迭代是否成功(
14、1 为成功，0 为失败)；所得具体结果见附录四。相对而言，极小似然迭代法减小误差的效果更好些，误差的方差也很小。5.2 针对问题二：问题一给定了多组实际测量数据，我们根据这些数据采用均值迭代法进行了处理。同时给出了我们自己另外创造的模型：极小似然迭代法，即努力使得41 个随机误差尽量小，均方差也尽量小以达到使各组数据努力保持一致。此处，我们将在问题一所建立模型的基础上，对极小似然迭代法进行改进与扩充：利用 Matlab 函数 interp1 拟合附录三中两种误差，可以得到差值多项式。相应的 Matlab 程序如下：见附录三结果如图：同样：垂直速度误差拟合曲线如图：对于飞行器的 x,y 方向，由
15、问题一的讨论可认为其基本不存在误差。根据以上误差数据可以结合具体高度和速度数据矫正出较为真实的数据。5.3 针对问题三：针对导弹的坐标修正可以采用前面所述的误差修正模型来求解，利用地面观测站测量的数据，得到惯性导航系在任意时刻t的位置（Xi，Yi，Zi），观测站每间隔T时间进行一次数据测量，可测的导弹的空间位置（Xi，Yi，Zi）和瞬时速度（Vix，Viy，Viz）。首先为消除随机误差，对空间位置坐标和瞬时速度的数据进行极小似然迭代分别得到（Xi，Yi，Zi）、 （Vix，Viy，Viz）.根据问题二得到的误差与t之间的差值函数（即问题二中用到的interp1）可估算出误差的近似值,delta
16、 和delta。然后采用与处理高度误差和垂直方向速度误差相同的方法（问题二中的极小似然法模型）估测Xdelta和Xdelta，Ydelta和Ydelta六、模型评价对于问题一，我们采用了均值迭代法和极小似然迭代法处理对于噪声干扰对空间坐标的影响，我们的程序可以设定数据的方差在一定范围内，使得误差更小，但还是有缺陷，我们查阅资料，可能通过联邦卡尔曼滤波能更好的消除噪声干扰对空间坐标产生的影响。对于问题二，我们通过上面修正后的数据与通过速度坐标修正后的数据利用改进后的极小似然迭代法拟合出差值多项式，使得随机误差和误差方差都尽量小，这种方法很巧妙的将飞行器的空间坐标误差进行修正，而且效果非常的好，很
17、好的避开了因常量误差的积累导致飞行器坐标数据发生漂移的现象。对于问题三，因前面我们通过改进后的极小似然迭代法后的数据有很好的修正误差的效果，所以对于具体的飞行器，我们依然采用的是改进后的极小似然迭代法来修正误差！虽然对于噪声干扰对于空间的坐标得到了修正，但对于电子仪器的精度干扰不能很好的修正。七、模型改进本个模型对误差的修正并没有达到真正最优的效果，因为数据的组数较少，加之我们对一些程序设计没有很好的掌握，所以还是略有缺陷，经查阅资料，早消除噪声干扰产生的影响上，联邦卡尔曼滤波法更加有效，如果能在问题一采用的卡尔曼滤波和均值迭代处理中通过程序设计消除噪声影响，加上后面的改进后的极小似然迭代法进
18、行修正效果会更好些。八、参考文献1、徐明华 张燕新 李志林 编 数值分析附录一len=length(Vdelta);averageV=sum(abs(Vdelta)/len;averageZ=sum(abs(Vdelta)/len;Dv=sum(abs(Vdelta(1:len)-averageV).2);Dv=Dv0.5;Dz=sum(abs(Vdelta(1:len)-averageZ).2);Dz=Dz0.5;if DvepsonV break;end附录二Function z,v,Zdelta,Vdelta,count,Dv,Dz,flag=fun1(z0,v0,epsonV,epso
19、nZ,maxCount)z=z0;v=v0;n=length(z0);av1=zeros(1,n-1);Zdelta=zeros(1,n-1);av2=zeros(1,n-1);Vdelta=zeros(1,n-1);Dz=0;Dv=0;flag=0;for count=1:maxCountav1(1:n-1)=(z0(2:n)-z0(1:n-1)/0.15;av2(1:n-1)=(v0(1:n-1)+v0(2:n)/2;Vdelta=av1-av2;Zdelta=Vdelta*0.15;len=length(Vdelta);averageV=sum(abs(Vdelta)/len;avera
20、geZ=sum(abs(Vdelta)/len;Dv=sum(abs(Vdelta(1:len)-averageV).2);Dv=Dv0.5;Dz=sum(abs(Vdelta(1:len)-averageZ).2);Dz=Dz0.5;if DvepsonV break;endz(1)=sum(z0)/length(z0);v(1)=sum(v0)/length(v0);z(2)=z0(1)+Zdelta(1);for i=2:n-1z(i+1)=z0(i)+Zdelta(i);v(i)=(z(i+1)-z(i-1)/0.3;endv(n)=v(n-1)+Vdelta(n-1);z0=z;v0
21、=v;endreturn;附录三t=0.15:0.15:6;z=0.0858 -0.0848 0.0831 -0.0809 0.0783 -0.0751 0.0716 -0.0678 0.0639 0.0598 0.0557 -0.0517 0.0479 -0.0443 0.041 -0.038 0.0354 -0.03320.0313 -0.0299 0.0287 -0.0279 0.0273 -0.0269 0.0266 -0.0264 0.0261-0.0258 0.0252 -0.0245 0.0236 -0.0223 0.0208 -0.019 0.0168 -0.01450.011
22、8 -0.009 0.0061 -0.003;z=abs(z);Z1=interp1(t,z,t,spline);z1=polyval(p,t);plot(t,z,k+,t,z1,r);附录四z(m) vz(m/s) Zdelta(m) Vdelta(m/s)z*(m) v*(m/s)高 度 垂 直 方 向 (高 度 )速 度 高 度 误 差 垂 直 速 度 误 差修 正 后 高 度 修 正 后 垂 直 速 度510.3520237 -0.812188133 0.0858 0.5718 510.5635 -0.0141510.2121547 -0.692034667 -0.0848 -0.56
23、51 510.6475 -0.0095510.2149607 -0.80888028 0.0831 0.5543 510.5606 -0.0183510.1983958 -0.61262966 -0.0809 -0.5397 510.642 -0.0055510.1080095 -0.611442845 0.0783 0.5217 510.559 -0.0219510.3356322 -0.747974657 -0.0751 -0.5007 510.6354 -0.0023510.0836268 -0.43440305 0.0716 0.4774 510.5583 -0.0246510.442
24、1405 -0.603809599 -0.0678 -0.4522 510.628 -0.0002510.2034774 -0.665666953 0.0639 0.4257 510.5582 -0.0262510.1592621 -0.789401121 -0.0598 -0.3986 510.6202 0.0007510.0681463 -0.476828359 0.0557 0.3714 510.5585 -0.0266510.2475029 -0.628712183 -0.0517 -0.3448 510.6122 0.0004510.2802798 -0.718854549 0.04
25、79 0.3192 510.5586 -0.0259510.4076986 -0.544532286 -0.0443 -0.2952 510.6044 -0.0011510.4327193 -0.527395078 0.041 0.2731 510.5583 -0.0243510.2397617 -0.572231206 -0.038 -0.2533 510.5971 -0.0034510.4504262 -0.701085244 0.0354 0.2359 510.5572 -0.0222510.2929935 -0.63122458 -0.0332 -0.2211 510.5905 -0.
26、0062510.3129799 -0.561235904 0.0313 0.2088 510.5554 -0.0198510.4911516 -0.487937208 -0.0299 -0.1991 510.5845 -0.0089510.2900452 -0.570562304 0.0287 0.1916 510.5527 -0.0175510.2421483 -0.521424799 -0.0279 -0.186 510.5793 -0.0111510.3149417 -0.494330988 0.0273 0.1821 510.5494 -0.0158510.0247663 -0.741
27、681622 -0.0269 -0.1794 510.5745 -0.0128510.1027471 -0.723808815 0.0266 0.1775 510.5455 -0.0152510.317599 -0.504467461 -0.0264 -0.1759 510.57 -0.0139510.249558 -0.70365197 0.0261 0.1741 510.5414 -0.0158510.2451786 -0.785142757 -0.0258 -0.1717 510.5652 -0.0141510.1042307 -0.820964576 0.0252 0.1682 510
28、.5371 -0.0171510.1029878 -0.691124582 -0.0245 -0.1635 510.5601 -0.0133510.0710206 -0.559037038 0.0236 0.157 510.5332 -0.0189510.0260389 -0.639731612 -0.0223 -0.1488 510.5544 -0.0116510.0317023 -0.671942967 0.0208 0.1385 510.5297 -0.0208510.4294694 -0.509338316 -0.019 -0.1263 510.5482 -0.0095510.1992
29、947 -0.602289998 0.0168 0.1122 510.5268 -0.0226510.150653 -0.71394819 -0.0145 -0.0964 510.5414 -0.0074510.3240992 -0.496014175 0.0118 0.0789 510.5246 -0.0239510.42705 -0.662702202 -0.009 -0.0602 510.5342 -0.0054510.0885538 -0.596728181 0.0061 0.0405 510.523 -0.0242510.0590776 -0.573856336 -0.003 -0.
30、02 510.527 -0.0036510.4997458 -0.422227046 510.5219 -0.0238附录五迭 代 一 次 迭 代 两 次 迭 代 三 次510.2375 510.1709 510.1978510.3267 510.2601 510.2332510.3295 510.2629 510.236510.3129 510.2463 510.2732510.2226 510.2892 510.2623510.2211 510.2877 510.2608510.1982 510.2648 510.2379510.3276 510.261 510.2341510.318 5
31、10.2514 510.2245510.2738 510.2072 510.2341510.1827 510.2493 510.2224510.133 510.1996 510.2265510.1657 510.2323 510.2592510.2931 510.2265 510.2534510.3182 510.2516 510.2247510.3543 510.2877 510.2608510.3359 510.2693 510.2424510.1784 510.245 510.2719510.1984 510.265 510.2381510.3766 510.31 510.2831510
32、.1755 510.2421 510.269510.1276 510.1942 510.2211510.2004 510.267 510.2401510.1393 510.2059 510.2328510.2173 510.2839 510.257510.203 510.2696 510.2427510.135 510.2016 510.2285510.1306 510.1972 510.2241510.2188 510.2854 510.2585510.2175 510.2841 510.2572510.1856 510.2522 510.2253510.1406 510.2072 510.2341510.1463 510.2129 510.2398510.3149 510.2483 510.2752510.3138 510.2472 510.2741510.2652 510.1986 510.2255510.2095 510.2761 510.2492510.3125 510.2459 510.2728510.2031 510.2697 510.2428510.1736 510.2402 510.2671510.3852 510.3186 510.2917
