导航原理大作业
导航原理期末考试题及答案
导航原理期末考试题及答案一、选择题（每题2分，共20分）1. 卫星导航系统主要利用的是以下哪种信号进行定位？A. 无线电信号B. 声波信号C. 光波信号D. 电磁波信号答案：A2. GPS系统是由多少颗卫星组成的？A. 24颗B. 36颗C. 48颗D. 60颗答案：A3. 以下哪个是卫星导航系统中的误差来源？A. 大气延迟B. 卫星轨道误差C. 接收机时钟误差D. 所有以上答案：D4. 卫星导航接收机的定位精度通常受到哪些因素的影响？A. 卫星几何分布B. 信号传播路径C. 接收机硬件性能D. 所有以上答案：D5. 卫星导航系统在什么情况下可以提供三维定位？A. 接收到至少4颗卫星的信号B. 接收到至少3颗卫星的信号C. 接收到至少2颗卫星的信号D. 接收到至少1颗卫星的信号答案：A6. 什么是差分GPS（DGPS）？A. 一种提高GPS定位精度的技术B. 一种用于军事的GPS系统C. 一种GPS信号干扰技术D. 一种GPS信号加密技术答案：A7. 以下哪个不是卫星导航系统的主要组成部分？A. 卫星B. 地面控制站C. 用户设备D. 通信卫星答案：D8. 卫星导航系统通常使用哪种频率的信号？A. 超高频B. 甚高频C. 特高频D. 微波答案：D9. 卫星导航系统的时间同步是通过什么实现的？A. 原子钟B. 石英钟C. 机械钟D. 电子钟答案：A10. 卫星导航系统中，哪个参数用于计算用户和卫星之间的距离？A. 信号频率B. 信号时间C. 信号强度D. 信号相位答案：D二、简答题（每题10分，共30分）1. 简述卫星导航系统的工作原理。
答案：卫星导航系统通过向地球发射无线电信号，用户设备接收来自至少四颗卫星的信号，根据信号传播的时间差计算出用户与每颗卫星之间的距离。
通过这些距离信息，结合卫星的精确位置，用户设备可以计算出自己的三维位置。
2. 解释什么是时间差定位（Time Difference of Arrival, TDoA）。
北京航空航天大学北航仪器科学与光电工程学院先进惯性导航惯导惯性导航原理大作业1
《惯性导航原理》第一次大作业姓名： XXXXX学号： SY13172XX邮箱： XXXXXXXXXXXXXXXXXXX2013/10/20一、 作业目的本次作业基于指北方位的平台惯导系统的相关原理，通过用matlab 对相关数据的处理，来模拟载体的速度、位置等相关参数，深入理解惯性导航系统的原理。
二、 设计要求1. 初始经度为：116.344762072818度；纬度为：39.981430918136度；高度为：40.8236米；初始姿态角为[0 0 0]；初始速度为0米/秒；飞行高度不变。
2. 排列顺序：第一列：数据包序号；第二至四列：分别为东、北、天向陀螺仪角速率信息wib_INSc （单位：rad/s ）；第五至七列：分别为东、北、天向比力信息f_INSc （单位：m/s^2）。
3. 以经度为横轴，纬度为纵轴（单位均转换为：度）作出系统位置坐标曲线图。
4. 以列表形式给出系统纬度、经度、东向速度、北向速度的终点值。
三、 设计原理惯性导航的基本原理是测量载体的加速度，将加速度积分计算出载体的速度，由速度积分算得载体的相对位置。
对于加速度的测量可以用加速度计测得比力（单位质量相对惯性空间的加速度与引力加速度之差），由于是指北方位平台系统，利用加速度计提供的东、北、天三个方向的比力数据（ff f p zp yp x,,）算得当时平台的瞬间速度，进而计算出每个数据采集点平台的即时经纬度，数据由fw 中保存的比力f （单位：2/s m ）和陀螺仪角速率信息w （单位：rad/s ）提供，最终得到平台在运动过程中任意时刻的速度和位置信息。
四、 设计流程1. 通过比力算加速度（1）基于比力方程：惯导系统的基本方程为eie ie e ie e G r dt drdtr d f -⨯⨯+=+)(222ωωω， 其中edt rd 22 是载体相对地球的加速度；e ie dt dr ω2 是哥氏加速度；)(r ie ie ⨯⨯ωω是地球自转引起的向心加速度；eG 是地球引力加速度。
北斗导航的原理和应用实例
北斗导航的原理和应用实例1. 北斗导航的原理北斗导航系统是中国自主研发的卫星导航系统，具有全球覆盖能力。
该系统由一组卫星、地面站和用户终端组成，通过卫星之间的信号传输和用户终端的接收，实现对地球上任意一点的定位。
北斗导航系统的原理如下：•卫星定位：北斗系统使用了全球定位系统（GPS）的原理，通过在卫星上携带精确的时钟，并以高速向地球发送时钟信号，用户终端接收到多个卫星发出的信号后，通过测量信号的传输时间差，将卫星的位置推算出来，从而实现定位功能。
•位置更新：北斗系统中的卫星定位系统会定期向地面站发送信号，地面站将这些信号转发给用户终端，用户终端通过接收这些信号并测量信号传输时间差，从而实现位置更新功能。
•数据传输：北斗系统不仅可以传输定位信息，还可以传输其他各种类型的数据，例如天气信息、交通信息等。
用户终端通过接收卫星发出的信号，获得所需的数据。
2. 北斗导航的应用实例2.1 航海导航北斗导航系统在航海领域的应用非常广泛。
船只可以通过北斗系统获得准确的定位信息，从而实现航线规划、航行管理和预测位置等功能。
此外，北斗系统还可以提供海洋气象信息和海图更新等服务，大大提高了航海的安全性和准确性。
2.2 土地测量北斗导航系统在土地测量领域的应用也十分重要。
使用北斗系统可以高精度地测量大地坐标和高程等数据，为土地测绘、城市规划和土地管理等提供了重要数据支持。
此外，北斗系统还可以配合其他测量工具实现测距、测角等功能，提高了测量的效率和准确性。
2.3 物流管理北斗导航系统在物流管理中发挥着重要的作用。
物流公司可以通过北斗系统追踪货物的位置和运输过程，实时监控货物的流向和货车的运行状态，提高物流管理的效率和准确性。
同时，北斗系统还可以提供实时的天气信息和路况信息，为物流公司的决策提供支持。
2.4 紧急救援北斗导航系统在紧急救援中起到了重要的作用。
当发生灾害或紧急事件时，救援人员可以利用北斗导航系统快速定位受灾地区和受灾人员，提供准确和及时的救援服务。
哈工大惯性技术(导航原理)大作业讲解
Assignment of Inertial Technology 《惯性技术》作业Autumn 2015Assignment 1: 2-DOF response simulation A 2-DOF gyro is subjected to a sinusoidal torque with amplitude of 4 g.cmand frequency of 10 Hz along its outer ring axis. The angular moment of its rotor is 10000 g.cm/s , and the angular inertias in its equatorial plane are both 4 g.cm/s 2. Please simulate the response of the gyro within 1 second, and present whatever you can discover or confirm from the result.In this simulation, we are going to discuss the response of a 2-DOF gyro to sinusoidal torque input. According to the transfer function of the 2-DOF gyro, the outputs can be expressed as: 12222()()()()yx y x y x y J H s M s M s J J s H s J J s H α=-++ 12222()()()()x y x x y x y J H s M s M s J J s H s J J s H β=+++ The original system transfer function is a 2-input, 2-output coupling system. But the given input only exists one input, we can treat the system as 2 separate FIFO systemAs a consequence, we can establish the block diagram of the system in simulink in Fig 1.1. Substitute the parameters into the system and input, then we have the input signals as follow: 0,4sin(20)y ox M M t π==Then the inverse Laplace transform of the output equals the response of the gyro in timedomain as follows:02222000200222200()sin sin ()()()cos cos ()()ox a ox a x a x a ox ox a ox a a a a a M M t t t J J M M M t t t H H H ωαωωωωωωωωωβωωωωωωωω=-+--=+---Fig 1.1 The block diagram of the system in simulinkAnd the simulation results in time domain within 1 second are shown in the follwing pictures. Fig1.2 is the the output of outer ring ()tα. Fig1.3 is the output of inner ring ()tβ. Fig1.4 is the trajectory of 2-DOF gyro’s response to sinusoidal input. As we can see from the Fig1.3,there are obvious sawtooth wave in the output of the inner ring. It’s a unexpected phenomenon in my original theoretical analysis.Fig1.2 The output of inner ring ()tβFig 1.3 The output of outer ring ()tαI believe the sawtooth wave is caused by the nutation. For the frequency of the nutation isobtained as100002500/3974eHrad s HzJω===≈, which is far higher than the frequencyof the applied sinusoidal torque , namelyaωω.Fig 1.4 Trajectory of 2-DOF gyro’s response to sinusoidal input The trajectory of 2-DOF gyro’s response to sinusoidal input are shown in Fig1.4. As we can see, it’s coupling of X and Y channel scope output. The overall shape is an ellipse, which is not perfect for there are so many sawteeth on the top of it.Note that the major axis of ellipse is in the direction of the forced procession, amplitude of which isoxMHω, whereas the minor axis is in the direction of the torsion spring effects, with amplitude oxaMHω.The nutation components are much smaller than that of the forced vibration, which can be eliminated to get the clear static response.2200222()sin sin()cos(1cos)()ox oxa ax aox ox oxa aa a a aM Mt t tJ HM M Mt t tH H Hαωωωωωωβωωωωωωω≈≈-≈-=--To prove it, we eliminate the effects of the nutation namely the quadratic term in the denominator and get Fig 1.5, which is a perfect ellipse. We can conclude that when input to the 2-DOF gyro is sinusoidal torque, the gyro will do an ellipse conical pendulum as a static response, including procession and the torsion spring effects, together with a high-frequency vibration as the dynamic response.Fig 1.5 Trajectory of the gyro’s response without nutationAssignment 2: Single-axis INS simulation2.1 description of the problemA magnetic levitation train is being tested along a track running north-south. It first accelerates and then cruises at a constant speed. Onboard is a single-axis platform INS, working in the way described by the courseware of Unit 5: Basic problems of INS. The motion information and Earth parameters are shown in table 2-1, and the possible error sources are shown in Table2-2.Fig2.1 The sketch map of the single-axis INS problemYou are asked to simulate the operation of the INS within 10,000 seconds, and investigate, first one by one and then altogether, the impact of these error sources on the performance of the INS.The block diagram in the courseware might be of some help. However, there lurks an inconspicuous error, which you have to correct before you can obtain reasonable results.There are one core relevant formula, to get the specific form of its solution, we should substitute the unknown parameters.(1)()c N ay A K y ga=∆++-Firstly, the input signal is accelerometer of the platform, and the velocity of the platform is the integration of the acceleration./py y ydty Rω=+=-⎰The acceleration along Yp may contains two parts:cos gsin gypf y yααα=-≈-When accelerometer errors are concerned, the output of accelerometer will be:(1)N a yp Na K f A=++When gyro errors concerned:'(1)p g pKωωε=++Only αis unknown:[(1)]tg pttptK dtdtααωεϕϕω=+++-∆∆=⎰⎰And the reference block diagram and simulink block diagram are as follows in Fig2.2, Fig2.3. There is a small fault in the reference block, which is that the sign of the marked add operation should be positive instead of negative.The results of the simulation are shown in Fig 2.4 to Fig 2.13.Fig2.2 The reference block diagram in the courseware(unrectified)Fig2.3 The simulink block diagram for Assignment 22.2 results and analysis of the problemFig.2.4 real acceleration,velocity and displacement output without error sourcesKFig2.5 position bias when only accelerometer scale factor error exists0.0001aFig2.6 position bias output when only gyro scale factor error exists 0.0001g K =Fig2.7 position bias when only acceleromete bias error exists 20.0002/N A m s ∆=Fig2.8 position bias output when only initial velocity error exists 200.01/y m s ∆=Fig2.9 position bias output when only initial position error exists 010y m ∆=Fig2.10 position bias output bias when only gyro bias errorexists 0.000000024240.00681/3/h rad s ε=︒=Fig2.11 position bias output when only initial platform misalignment angle exists0.000012120''345radα==Fig2.12 output considering all the error sourcesFig2.13 position bias output considering all error sourcesAs we can see in the above simulation results, if there is no error we can navigate the train ’s motion correctly, which comes from north to the south as shown in Fig2.4, beginning with an constant acceleration within 60 seconds then cruises at a constant speed, approximately 65 m/s. However, the situation will change a lot when different errors put into the simulation. The initial position error 010y m ∆= effects least as Fig.2.8, for this error doesn ’t enter into the closed loop and it won ’t influence the iterative process. The position bias is constant and can be negligible.In the second case, when the accelerometer scale factor error exists, 0.0001a K =, as shown in Fig2.5, the result are stable and almost accurate, the position bias is a sinusoidal output.So it is with the accelerometer bias error situation, 20.0002/N A m s ∆=, in Fig2.7, the initialvelocity error, 200.01/y m s ∆= in Fig2.9, and the initial platform misalignment angle,05''α=, in Fig2.11. However, the influence degrees of the different factors are not in the samemagnitude. The accelerometer scale factor influences the least with magnitude of 5, then the initial velocity larger magnitude of 8, and the accelerometer bias magnitude of 25. The influence of the initial platform misalignment angle is much more significant with a magnitude of 150. All the navigation bias in the second kind case is sinusoidal, which means they ’re limited and negligible as time passes by.In the third case, such as the gyro scale factor error situation, 0.0001g K =, in Fig2.6, and the gyro bias error,0.01/h ε=︒, results in Fig2.10, effects the most significant, the trajectory of the navigation disvergence accumulated as time goes. The position bias is a combination of sinusoidal signal and ramp signal. They also show that the longitudinal and distance errors resulted from gyro drifts are not convergent in time. It means the errors in the gyroscope do most harm to our navigation. And due to the significant influence of the gyro bias errors and the gyroscope scale factor error, results considering all the error sources disverge, and the navigationposition of the motion will be away from the real motion after a enough long time, as shown in Fig2.10. The gyro bias error is the most significant effect factor of all errors. By the time of 10000s, it has reaches 1600m, and it’s nearly the quantity of the position bias considering all error sources. Through contrasting all the results, We can conclude that the gyro bias error is the main component of the whole position bias.Impression of the Whole simulation experimentThrough contrasting all the results, We can conclude that the gyro bias error is the main component of the whole position bias, and the the gyro bias or the drift error do most harm to our navigation. So it is a must for us to weaken or eliminate it anyway. In spite of all the disadvantages discussed above, the INS still shows us a relatively accurate results of single-axis navigation.Assignment 3: SINS simulation3.1 Task descriptionA missile equipped with SINS is initially at the position of 46o NL and 123 o EL, stationary on a launch pad. Three gyros, GX, GY, GZ, and three accelerometers, AX, AY, AZ, are installed along the axes Xb, Yb, Zb of its body frame respectively.Case 1: stationary testThe body frame of the missile initially coincides with the geographical frame, as shown in the figure, with its pitching axis Xb pointing to the east, rolling axis Yb to the north, and azimuth axis Zb upward. Then the body of the missile is made to rotate in 3 steps:(1) -22o around Xb(2) 78o around Yb(3) –16o around ZbFig 3.1 Introduction to assignment 3After that, the body of the missile stops rotating. You are required to compute the final outputs of the three accelerometers on the missile, using quaternion and ignoring the device errors. It is known that the magnitude of gravity acceleration is g = 9.8m/s2.Case 2: flight tes tInitially, the missile is stationary on the launch pad, 400m above the sea level. Its rolling axis is vertical up,and its pitching axis is to the east. Then the missile is fired up. The outputs of the gyros and accelerometers are both pulse numbers. Each gyro pulse is an angular increment of 0.01arc-sec, and each accelerometer pulse is 1e-7g, with g =9.8m/s2. The gyro output frequency is 200Hz, and the accelerometer’s is 10Hz. The outputs of the gyros and accelerometers within 1315s are stored in a MA TLAB data file named imu.mat, containing matrices gm of 263000× 3 and am of 13150× 3 respectively. The format of the data is as shown in the tables, with 10 rows of each matrix selected. Each row represents the outputs of the type of sensors at each sampling time. The Earth can be seen as an ideal sphere, with radius 6371.00km and spinning rate 7.292× 10-5 rad/s, The errors of the sensors are ignored, so is the effect of height on the magnitude of gravity. The outputs of the gyros are to be integrated every 0.005s. The rotation of the geographical frame is to be updated every 0.1s, so are the velocities and positions of the missile.You are required to:(1) compute the final attitude quaternion, longitude, latitude, height, and east, north, vertical velocities of the missile.(2) compute the total horizontal distance traveled by the missile.(3) draw the latitude-versus-longitude trajectory of the missile, with horizontal longitude axis.(4) draw the curve of the height of the missile, with horizontal time axis.Fig 3.2 simplified navigation algorithm for SINS 3.2Procedure code3.2.1 Sub function code：quaternion multiply code:function [q1]=quml(q1,q2);lm=q1(1);p1=q1(2);p2=q1(3);p3=q1(4);q1=[lm -p1 -p2 -p3;p1 lm -p3 p2;p2 p3 lm -p1;p3 -p2 p1 lm]*q2;endquaternion inversion code:function [qni] =qinv(q)q(1)=q(1);q(2)=-q(2);q(3)=-q(3);q(4)=-q(4);qni=q;end3.2.2case1 DCM algorithm：function ans11cz=[cos(-22/180*pi) sin(-22/180*pi) 0 ;-sin(-22/180*pi) cos(-22/180*pi) 0;0 0 1]; %The third rotation DCMcx=[1 0 0 ;0 cos(78/180*pi) sin(78/180*pi);0 -sin(78/180*pi) cos(78/180*pi)]; %The first rotation DCMcy=[cos(-16/180*pi) 0 -sin(-16/180*pi);0 1 0;sin(-16/180*pi) 0 cos(-16/180*pi)]; %The second rotation DCM A=cz*cy*cx*[0;0;-9.8]end3.2.3case1 quaternion algorithm：function ans12g=[0;0;0;-9.8];q1=[cos(-11/180*pi);sin(-11/180*pi);0;0]; %The first rotation quaternionq2=[cos(39/180*pi);0;sin(39/180*pi);0]; %The second rotation quaternionq3=[cos(-8/180*pi);0;0;-sin(-8/180*pi)]; %The third rotation quaternionr=quml(q1,q2); %call the quaternion multiplication subfunction q=quml(r,q3);P1=[q(1) q(2) q(3) q(4);-q(2) q(1) q(4) -q(3);-q(3) -q(4) q(1) q(2);-q(4) q(3) -q(2) q(1)];P2=[q(1) -q(2) -q(3) -q(4);q(2) q(1) q(4) -q(3);q(3) -q(4) q(1) q(2);q(4) q(3) -q(2) q(1)];P=P1*P2;gn=P*g;gn=gn(2:4)end3.2.4 case2 SINS quaternion algorithm code：function ans2clc;clear;%parameters initializing：T=0.005;K=13150;R=6371000; %radius of earthwE=7.292*10^(-5); %spinning rate of earthQ=zeros(4, 263001) ;%quaternion matrix initializinglongitude=zeros(1,13151);latitude=zeros(1,13151);H=zeros(1,13151); %altitude matrixQ(:,1)=[cos(45/180*pi);-sin(45/180*pi);0;0]; % initial quaternion longitude(1)=123; %initial longitudelatitude(1)=46; % initial latitudeH(1)=400; %initial altitudelength=0;g=9.8;vE = zeros(1,13151); %eastern velocityvN = zeros(1,13151); %northern velocityvH = zeros(1,13151); %upward velocityvE(1)=0;vN(1)=0;vH(1)=0;load imu.mat %data loading%main calculation sectionfor N=1:Kq1=zeros(4,11);q1(:,1)=Q(:,N);for n=1:20 % Attitude iteration wx=0.01/(3600*180)*pi*gm((N-1)*10+n,1); % Angle incrementwy=0.01/(3600*180)*pi*gm((N-1)*10+n,2);wz=0.01/(3600*180)*pi*gm((N-1)*10+n,3);w=[wx,wy,wz]';normw=norm(w); % Norm calculationW=[0,-w(1),-w(2),-w(3);w(1),0,w(3),-w(2);w(2),-w(3),0,w(1);w(3),w(2),-w(1),0];I=eye(4);S=1/2-normw^2/48;C=1-normw^2/8;q1(:,n+1)=(C*I+S*W)*q1(:,n);Q(:,N+1)=q1(:,n+1);endWE=-vN(N)/R; % rotational angular velocity component of a geographic coordinate systemWN=vE(N)/R+wE*cos(latitude(N)/180*pi);WH=vE(N)/R*tan(latitude(N)/180*pi)+wE*sin(latitude(N)/180*pi);attitude=[WE,WN,WH]'*T; %correction of the quaternion by updating the rotation of geographic coordinatenormattitude=norm(attitude);e=attitude/normattitude;QG=[cos(normattitude/2);sin(normattitude/2)*e];Q(:,N+1)=quml(qinv(QG),Q(:,N+1));fx=1e-7*9.8*am(N,1); %specific force measured by accelerometerfy=1e-7*9.8*am(N,2);fz=1e-7*9.8*am(N,3);Fb=[fx fy fz]';F=quml(Q(:,N+1),quml([0;Fb],qinv(Q(:,N+1)))); %The specific force is decomposed into geographic coordinate system.FE(N)=F(2);FN(N)=F(3);FU(N)=F(4);%calculate the velocity of the vehicle:VED(N)=FE(N)+vE(N)*vN(N)/R*tan(latitude(N)/180*pi)-(vE(N)/R+2*wE*cos(latitude(N)/ 180*pi))*vH(N)+2*vN(N)*wE*sin(latitude(N)/180*pi);VND(N)=FN(N)-2*vE(N)*wE*sin(latitude(N)/180*pi)-vE(N)*vE(N)/R*tan(latitude(N)/180 *pi)-vN(N)*vH(N)/R;VUD(N)=FU(N)+2*vE(N)*wE*cos(latitude(N)/180*pi)+(vE(N)^2+vN(N)^2)/R-g;%integration and get the relative velocity of vehicle:vE(N+1)=VED(N)*T+vE(N);vN(N+1)=VND(N)*T+vN(N);vH(N+1)=VUD(N)*T+vH(N);% integration and get the position of vehicle:longitude(N+1)=vE(N)/(R*cos(latitude(N)/180*pi))*T/pi*180+longitude(N);latitude(N+1)=vN(N)/R*T/pi*180+latitude(N);H(N+1)=vH(N)*T+H(N);length=sqrt((vE(N))^2+(vN(N))^2)+length;endfigure(1) %picture the longitude-latitude curve of the motion within 1315 secondstitle(' longitude-latitude ');hold ongrid onplot(longitude,latitude);figure(2) % picture the altitude curve of the motion within 1315 secondstitle('altitude');hold ongrid onplot(0:13150,H);3.3Simulation outputs and results analysisIn case 1, the DCM algorithm and quaternion results are the same as follow:A =7.53165.9788-1.8892And the results suggest the solution of DCM and quaternion is equivalence in this problem.In case 2, the latitude-versus-longitude trajectory of the missile(with horizontal longitude axis) and the altitude curve of the missile are shown as follows in Fig 3.3 and Fig 3.4.Fig3.3 the latitude-versus-longitude trajectory of the missile(with horizontal longitude axis)Fig3.4the altitude curve of the missileImpression of the Assignment 3In this assignment, we simulate the missile navigation situation in a real problem, as the problem we have done before. I did this based on my own original code, but to my surprise, it’s harder than I have expected. For the detailed thoughts for my program I have forgotten. I have to recheck the code line by line, But I am still troubled in correcting the initial parameters for many times. It has taught me a lesson, which is never to be egotistical for your ability to memory and the work you have done or understood.。
哈工大惯性技术(导航原理)大作业
Assignment ofInertial Technology《惯性技术》作业My Chinese NameMy Student ID 15S004001Autumn 2015Assignment 1: 2-DOF response simulationA 2-DOF gyro is subjected to a sinusoidal torque with amplitude of 4 g.cm and frequency of 10 Hz along its outer ring axis. The angular moment of its rotor is 10000 g.cm/s , and the angular inertias in its equatorial plane are both 4 g.cm/s 2. Please simulate the response of the gyro within 1 second, and present whatever you can discover or confirm from the result.In this simulation, we are going to discuss the response of a 2-DOF gyro to sinusoidal torque input. According to the transfer function of the 2-DOF gyro, the outputs can be expressed as:12222()()()()yx y x y x y J Hs M s M s J J s H s J J s H α=-++12222()()()()x yx x y x y J Hs M s M s J J s H s J J s H β=+++ The original system transfer function is a 2-input, 2-output coupling system. But the given input only exists one input, we can treat the system as 2 separate FIFO systemAs a consequence, we can establish the block diagram of the system in simulink in Fig 1.1. Substitute the parameters into the system and input, then we have the input signals as follow: 0,4sin(20)y ox M M t π==Then the inverse Laplace transform of the output equals the response of the gyro in timedomain as follows:0222200020222200()sin sin ()()()cos cos ()()ox a oxa x a x a ox ox a ox a a a a a M M t t t J J M M M t t t H H H ωαωωωωωωωωωβωωωωωωωω=-+--=+---Fig 1.1 The block diagram of the system in simulinkAnd the simulation results in time domain within 1 second are shown in the follwing pictures. Fig1.2 is the the output of outer ring ()t α. Fig1.3 is the output of inner ring ()t β. Fig1.4 is the trajectory of 2-DOF gyro ’s response to sinusoidal input. As we can see from the Fig1.3,there are obvious sawtooth wave in the output of the inner ring. It ’s a unexpected phenomenon in my original theoretical analysis.Fig1.2 The output of inner ring ()t βFig 1.3 The output of outer ring ()t αI believe the sawtooth wave is caused by the nutation. For the frequency of the nutation is obtained as010*******/3974e H rad s Hz J ω===≈, which is far higher than the frequencyof the applied sinusoidal torque , namely0a ωω .Fig 1.4 Trajectory of 2-DOF gyro ’s response to sinusoidal inputThe trajectory of 2-DOF gyro ’s response to sinusoidal input are shown in Fig1.4. As we can see, it ’s coupling of X and Y channel scope output. The overall shape is an ellipse, which is not perfect for there are so many sawteeth on the top of it.Note that the major axis of ellipse is in the direction of the forced procession, amplitude of which is 0ox M H ω, whereas the minor axis is in the direction of the torsion spring effects,with amplitude ox aM H ω.The nutation components are much smaller than that of the forced vibration, which can be eliminated to get the clear static response.22002220()sin sin ()cos (1cos )()ox oxaa x a ox ox oxa a a a a aM M t t t J H M M M t t t H H H αωωωωωωβωωωωωωω≈≈-≈-=--To prove it, we eliminate the effects of the nutation namely the quadratic term in the denominator and get Fig 1.5, which is a perfect ellipse. We can conclude that when input to the 2-DOF gyro is sinusoidal torque, the gyro will do an ellipse conical pendulum as a static response, including procession and the torsion spring effects, together with a high-frequency vibration as the dynamic response.Fig 1.5 Trajectory of the gyro’s response without nutationAssignment 2: Single-axis INS simulation2.1 description of the problemA magnetic levitation train is being tested along a track running north-south. It first accelerates and then cruises at a constant speed. Onboard is a single-axis platform INS, working in the way described by the courseware of Unit 5: Basic problems of INS. The motion informationand Earth parameters are shown in table 2-1, and the possible error sources are shown in Table2-2.Fig2.1 The sketch map of the single-axis INS problemYou are asked to simulate the operation of the INS within 10,000 seconds, and investigate,first one by one and then altogether, the impact of these error sources on the performance of theINS.The block diagram in the courseware might be of some help. However, there lurks aninconspicuous error, which you have to correct before you can obtain reasonable results.There are one core relevant formula, to get the specific form of its solution, we should substitute the unknown parameters.(1)()c N a y A K yga =∆++-Firstly, the input signal is accelerometer of the platform, and the velocity of the platform is the integration of the acceleration.0/p y y ydt yR ω=+=-⎰The acceleration along Yp may contains two parts:cos gsin g yp f y yααα=-≈- When accelerometer errors are concerned, the output of accelerometer will be:(1)N a yp N a K f A =++When gyro errors concerned:'(1)p g p K ωωε=++Onlyαis unknown:000[(1)]tg p t tp t K dt dtααωεϕϕω=+++-∆∆=⎰⎰And the reference block diagram and simulink block diagram are as follows in Fig2.2, Fig2.3. There is a small fault in the reference block, which is that the sign of the marked add operation should be positive instead of negative.The results of the simulation are shown in Fig 2.4 to Fig 2.13.Fig2.2 The reference block diagram in the courseware(unrectified)Fig2.3 The simulink block diagram for Assignment 22.2 results and analysis of the problemFig.2.4 real acceleration,velocity and displacement output without error sourcesKFig2.5 position bias when only accelerometer scale factor error exists0.0001aFig2.6 position bias output when only gyro scale factor error exists 0.0001g K =Fig2.7 position bias when only acceleromete bias error exists 20.0002/N A m s ∆=Fig2.8 position bias output when only initial velocity error exists 200.01/ym s ∆=Fig2.9 position bias output when only initial position error exists 010y m ∆=Fig2.10 position bias output bias when only gyro bias error exists 0.000000024240.00681/3/h rad s ε=︒=Fig2.11 position bias output when only initial platform misalignment angle exists0.000012120''345radα==Fig2.12 output considering all the error sourcesFig2.13 position bias output considering all error sourcesAs we can see in the above simulation results, if there is no error we can navigate the train ’s motion correctly, which comes from north to the south as shown in Fig2.4, beginning with an constant acceleration within 60 seconds then cruises at a constant speed, approximately 65 m/s. However, the situation will change a lot when different errors put into the simulation. The initial position error 010y m ∆= effects least as Fig.2.8, for this error doesn ’t enter into the closed loop and it won ’t influence the iterative process. The position bias is constant and can be negligible.In the second case, when the accelerometer scale factor error exists, 0.0001a K =, as shown in Fig2.5, the result are stable and almost accurate, the position bias is a sinusoidal output. So it is with the accelerometer bias error situation, 20.0002/N A m s ∆=, in Fig2.7, the initialvelocity error, 200.01/ym s ∆= in Fig2.9, and the initial platform misalignment angle, 05''α=, in Fig2.11. However, the influence degrees of the different factors are not in the samemagnitude. The accelerometer scale factor influences the least with magnitude of 5, then the initial velocity larger magnitude of 8, and the accelerometer bias magnitude of 25. The influence of the initial platform misalignment angle is much more significant with a magnitude of 150. All the navigation bias in the second kind case is sinusoidal, which means they ’re limited and negligible as time passes by.In the third case, such as the gyro scale factor error situation, 0.0001g K =, in Fig2.6, and the gyro bias error,0.01/h ε=︒, results in Fig2.10, effects the most significant, the trajectory of the navigation disvergence accumulated as time goes. The position bias is a combination of sinusoidal signal and ramp signal. They also show that the longitudinal and distance errors resulted from gyro drifts are not convergent in time. It means the errors in the gyroscope do most harm to our navigation. And due to the significant influence of the gyro bias errors and the gyroscope scale factor error, results considering all the error sources disverge, and the navigationposition of the motion will be away from the real motion after a enough long time, as shown in Fig2.10. The gyro bias error is the most significant effect factor of all errors. By the time of 10000s, it has reaches 1600m, and it’s nearly the quantity of the position bias considering all error sources. Through contrasting all the results, We can conclude that the gyro bias error is the main component of the whole position bias.Impression of the Whole simulation experimentThrough contrasting all the results, We can conclude that the gyro bias error is the main component of the whole position bias, and the the gyro bias or the drift error do most harm to our navigation. So it is a must for us to weaken or eliminate it anyway. In spite of all the disadvantages discussed above, the INS still shows us a relatively accurate results of single-axis navigation.Assignment 3: SINS simulation3.1 Task descriptionA missile equipped with SINS is initially at the position of 46o NL and 123 o EL, stationary on a launch pad. Three gyros, GX, GY, GZ, and three accelerometers, AX, AY, AZ, are installed along the axes Xb, Yb, Zb of its body frame respectively.Case 1: stationary testThe body frame of the missile initially coincides with the geographical frame, as shown in the figure, with its pitching axis Xb pointing to the east, rolling axis Yb to the north, and azimuth axis Zb upward. Then the body of the missile is made to rotate in 3 steps:(1) -22o around Xb(2) 78o around Yb(3) –16o around ZbFig 3.1 Introduction to assignment 3After that, the body of the missile stops rotating. You are required to compute the final outputs of the three accelerometers on the missile, using quaternion and ignoring the device errors. It is known that the magnitude of gravity acceleration is g = 9.8m/s2.Case 2: flight tes tInitially, the missile is stationary on the launch pad, 400m above the sea level. Its rolling axis is vertical up,and its pitching axis is to the east. Then the missile is fired up. The outputs of the gyros and accelerometers are both pulse numbers. Each gyro pulse is an angular increment of 0.01arc-sec, and each accelerometer pulse is 1e-7g, with g =9.8m/s2. The gyro output frequency is 200Hz, and the accelerometer’s is 10Hz. The outputs of the gyros and accelerometers within 1315s are stored in a MA TLAB data file named imu.mat, containing matrices gm of 263000× 3 and am of 13150× 3 respectively. The format of the data is as shown in the tables, with 10 rows of each matrix selected. Each row represents the outputs of the type of sensors at each sampling time. The Earth can be seen as an ideal sphere, with radius 6371.00km and spinning rate 7.292× 10-5 rad/s, The errors of the sensors are ignored, so is the effect of height on the magnitude of gravity. The outputs of the gyros are to be integrated every 0.005s. The rotation of the geographical frame is to be updated every 0.1s, so are the velocities and positions of the missile.You are required to:(1) compute the final attitude quaternion, longitude, latitude, height, and east, north, vertical velocities of the missile.(2) compute the total horizontal distance traveled by the missile.(3) draw the latitude-versus-longitude trajectory of the missile, with horizontal longitude axis.(4) draw the curve of the height of the missile, with horizontal time axis.Fig 3.2 simplified navigation algorithm for SINS 3.2Procedure code3.2.1 Sub function code：quaternion multiply code:function [q1]=quml(q1,q2);lm=q1(1);p1=q1(2);p2=q1(3);p3=q1(4);q1=[lm -p1 -p2 -p3;p1 lm -p3 p2;p2 p3 lm -p1;p3 -p2 p1 lm]*q2;endquaternion inversion code:function [qni] =qinv(q)q(1)=q(1);q(2)=-q(2);q(3)=-q(3);q(4)=-q(4);qni=q;end3.2.2case1 DCM algorithm：function ans11cz=[cos(-22/180*pi) sin(-22/180*pi) 0 ;-sin(-22/180*pi) cos(-22/180*pi) 0;0 0 1]; %The third rotation DCMcx=[1 0 0 ;0 cos(78/180*pi) sin(78/180*pi);0 -sin(78/180*pi) cos(78/180*pi)]; %The first rotation DCMcy=[cos(-16/180*pi) 0 -sin(-16/180*pi);0 1 0;sin(-16/180*pi) 0 cos(-16/180*pi)]; %The second rotation DCM A=cz*cy*cx*[0;0;-9.8]end3.2.3case1 quaternion algorithm：function ans12g=[0;0;0;-9.8];q1=[cos(-11/180*pi);sin(-11/180*pi);0;0]; %The first rotation quaternionq2=[cos(39/180*pi);0;sin(39/180*pi);0]; %The second rotation quaternionq3=[cos(-8/180*pi);0;0;-sin(-8/180*pi)]; %The third rotation quaternionr=quml(q1,q2); %call the quaternion multiplication subfunction q=quml(r,q3);P1=[q(1) q(2) q(3) q(4);-q(2) q(1) q(4) -q(3);-q(3) -q(4) q(1) q(2);-q(4) q(3) -q(2) q(1)];P2=[q(1) -q(2) -q(3) -q(4);q(2) q(1) q(4) -q(3);q(3) -q(4) q(1) q(2);q(4) q(3) -q(2) q(1)];P=P1*P2;gn=P*g;gn=gn(2:4)end3.2.4 case2 SINS quaternion algorithm code：function ans2clc;clear;%parameters initializing：T=0.005;K=13150;R=6371000; %radius of earthwE=7.292*10^(-5); %spinning rate of earthQ=zeros(4, 263001) ;%quaternion matrix initializinglongitude=zeros(1,13151);latitude=zeros(1,13151);H=zeros(1,13151); %altitude matrixQ(:,1)=[cos(45/180*pi);-sin(45/180*pi);0;0]; % initial quaternion longitude(1)=123; %initial longitudelatitude(1)=46; % initial latitudeH(1)=400; %initial altitudelength=0;g=9.8;vE = zeros(1,13151); %eastern velocityvN = zeros(1,13151); %northern velocityvH = zeros(1,13151); %upward velocityvE(1)=0;vN(1)=0;vH(1)=0;load imu.mat %data loading%main calculation sectionfor N=1:Kq1=zeros(4,11);q1(:,1)=Q(:,N);for n=1:20 % Attitude iteration wx=0.01/(3600*180)*pi*gm((N-1)*10+n,1); % Angle incrementwy=0.01/(3600*180)*pi*gm((N-1)*10+n,2);wz=0.01/(3600*180)*pi*gm((N-1)*10+n,3);w=[wx,wy,wz]';normw=norm(w); % Norm calculationW=[0,-w(1),-w(2),-w(3);w(1),0,w(3),-w(2);w(2),-w(3),0,w(1);w(3),w(2),-w(1),0];I=eye(4);S=1/2-normw^2/48;C=1-normw^2/8;q1(:,n+1)=(C*I+S*W)*q1(:,n);Q(:,N+1)=q1(:,n+1);endWE=-vN(N)/R; % rotational angular velocity component of a geographic coordinate systemWN=vE(N)/R+wE*cos(latitude(N)/180*pi);WH=vE(N)/R*tan(latitude(N)/180*pi)+wE*sin(latitude(N)/180*pi);attitude=[WE,WN,WH]'*T; %correction of the quaternion by updating the rotation of geographic coordinatenormattitude=norm(attitude);e=attitude/normattitude;QG=[cos(normattitude/2);sin(normattitude/2)*e];Q(:,N+1)=quml(qinv(QG),Q(:,N+1));fx=1e-7*9.8*am(N,1); %specific force measured by accelerometerfy=1e-7*9.8*am(N,2);fz=1e-7*9.8*am(N,3);Fb=[fx fy fz]';F=quml(Q(:,N+1),quml([0;Fb],qinv(Q(:,N+1)))); %The specific force is decomposed into geographic coordinate system.FE(N)=F(2);FN(N)=F(3);FU(N)=F(4);%calculate the velocity of the vehicle:VED(N)=FE(N)+vE(N)*vN(N)/R*tan(latitude(N)/180*pi)-(vE(N)/R+2*wE*cos(latitude(N)/ 180*pi))*vH(N)+2*vN(N)*wE*sin(latitude(N)/180*pi);VND(N)=FN(N)-2*vE(N)*wE*sin(latitude(N)/180*pi)-vE(N)*vE(N)/R*tan(latitude(N)/180 *pi)-vN(N)*vH(N)/R;VUD(N)=FU(N)+2*vE(N)*wE*cos(latitude(N)/180*pi)+(vE(N)^2+vN(N)^2)/R-g;%integration and get the relative velocity of vehicle:vE(N+1)=VED(N)*T+vE(N);vN(N+1)=VND(N)*T+vN(N);vH(N+1)=VUD(N)*T+vH(N);% integration and get the position of vehicle:longitude(N+1)=vE(N)/(R*cos(latitude(N)/180*pi))*T/pi*180+longitude(N);latitude(N+1)=vN(N)/R*T/pi*180+latitude(N);H(N+1)=vH(N)*T+H(N);length=sqrt((vE(N))^2+(vN(N))^2)+length;endfigure(1) %picture the longitude-latitude curve of the motion within 1315 secondstitle(' longitude-latitude ');hold ongrid onplot(longitude,latitude);figure(2) % picture the altitude curve of the motion within 1315 secondstitle('altitude');hold ongrid onplot(0:13150,H);3.3Simulation outputs and results analysisIn case 1, the DCM algorithm and quaternion results are the same as follow:A =7.53165.9788-1.8892And the results suggest the solution of DCM and quaternion is equivalence in this problem.In case 2, the latitude-versus-longitude trajectory of the missile(with horizontal longitude axis) and the altitude curve of the missile are shown as follows in Fig 3.3 and Fig 3.4.Fig3.3 the latitude-versus-longitude trajectory of the missile(with horizontal longitude axis)Fig3.4the altitude curve of the missileImpression of the Assignment 3In this assignment, we simulate the missile navigation situation in a real problem, as the problem we have done before. I did this based on my own original code, but to my surprise, it’s harder than I have expected. For the detailed thoughts for my program I have forgotten. I have to recheck the code line by line, But I am still troubled in correcting the initial parameters for many times. It has taught me a lesson, which is never to be egotistical for your ability to memory and thework you have done or understood.。
哈工大导航原理作业
Assignments of Inertial Navigation 《惯性导航》作业1Autumn 2018Assignment 1: Aircraft Displaying on the EarthAttached is a group of MA TLAB programs for displaying an aircraft near the surface of the Earth, as shown in figure 1.1 which is drawn by running the program main.m.If you input the longitude, latitude, heading, pitching, and rolling angles (all in degree), and click the updating button, then the aircraft are expected to be displayed at a proper position on the Earth, and with proper attitude. The coordinates of each vertex of the aircraft 3D model occupies a row in the matrix vtx.Before you click the updating button for the first time, the aircraft is hidden at the center of the Earth. Please rewrite the program redraw.m by adding the necessary rotation and shifting operations to the vertex coordinates of the aircraft. The acquired position and attitude information has been stored in the variables lon, lat, head, pitch and roll. With each updating, the coordinates of the aircraft in vtx in the MA TLAB drawing frame (same as the Earth frame) need to be re-calculated, instead of being kept the same as their values in the original 3D model (vtx0). That is, in the program redraw.m, you need to replace the command line vtx=vtx0 with your own codes for re-calculating vtx.And please explain the rationale behind your rewriting.作业报告一．理论基础a. 旋转变换矩阵在坐标系中，可以通过旋转角度知道一个点在新旧坐标系的坐标 在这样的旋转变化下存在如下的坐标变换矩阵：[x ′y ′z ′]=[1000cosαsinα0−sinαcosα][x y z] 同样的，绕y 轴和z 轴旋转的坐标变换矩阵为：[x ′y ′z ′]=[cosβ0−sinβ010sinβ0cosβ][xy z ] [x ′y ′z ′]=[cosγsinγ0−sinγcosγ0001][x y z ] 其中α、β、γ分别代表绕轴旋转的角度。
现代导航技术大作业(GPS)
GPS导航综述报告学号：1010200219，姓名：赵玲摘要：本文针对GPS导航技术，由GPS的发展历程进一步介绍了GPS系统的组成部分，基本工作原理和定位方法。
最后，在此基础上指出GPS的特点，并根据其特点列举了GPS的应用领域。
关键词：GPS导航、定位、测距一.引言导航的定义是“使运载体或人员从一个地方到另一个地方的科学”。
在日常生活中，我们每一个人都会进行某种形式的导航。
驱车去上班或步行去商店需要我们使用基本的导航技能。
对于我们大多数人来说这些技能需要我们的眼睛、常识和地标。
然而在一些情况下，需要更精确的知道我们的位置、预期的航向或达到期望目的地所需的时间，此时便需要不同于地标的导航装置。
这些导航装置也许是一个简单的时钟，以确定经过已知距离的速度；或者是汽车的里程表，以随时保持行驶的距离。
其他一些导航装置要发射电子信号，因而更复杂一些。
这些导航装置称为无线电导航装置。
全球卫星定位系统（Global Positioning System 简称GPS）是随着现代航天及无线电通讯科学技术的发展建立起来的一个高精度、全天候和全球性的无线电导航定位、定时的多功能系统。
它利用位于距地球2万多公里高的由24颗人造卫星组成的卫星网（即所谓“天网”），向地球不断发射定位及时间信号。
地球上的任何一个GPS 接收机，只要接收到四颗以上的卫星发出的信号，经过计算处理后，就可报出GPS 接收机（目标）的位置（经度、纬度、高度）、时间和运动状态（速度、航向）。
数据会适时地通过无线通讯网链传送至主控制基地中心，而后面具有强大地理信息处理、查询功能的电子地图上进行运动轨迹的显示，并能对准确位置、速度、运动方向、车辆状态等用户感兴趣的参数进行监控和查询，以确保车辆的安全，方便调度管理，提高远营效率。
该系统适用于公安、银行、部队、机场对车辆的监控和调度管理。
二.发展历程1958年底美国海军武器实验室就着手建立为美国军用舰艇导航服务的卫星系统，即“海军导航卫星系统”（Navy Navigation Satellite System－NNSS）。
导航原理测试题及答案
导航原理测试题及答案一、选择题（每题2分，共10分）1. 导航系统通常不包括以下哪项功能？A. 定位B. 通信C. 导航D. 娱乐答案：D2. GPS导航系统的全称是什么？A. Global Positioning SystemB. Global Positioning ServiceC. Global Positioning SatelliteD. Global Positioning Software答案：A3. 北斗导航系统是由哪个国家开发的？A. 美国B. 俄罗斯C. 中国D. 欧盟答案：C4. 以下哪个不是惯性导航系统的优点？A. 不依赖外部信号B. 精度高C. 价格低廉D. 可靠性高答案：C5. 以下哪个是全球卫星导航系统？A. GPSB. GLONASSC. BeiDouD. All of the above答案：D二、填空题（每题2分，共10分）1. 惯性导航系统（INS）是一种不依赖于外部信号的______系统。
答案：自主导航2. GPS系统由三部分组成：空间部分、地面控制部分和______。
答案：用户设备部分3. 北斗卫星导航系统由______、______和地面控制部分构成。
答案：空间段；用户段4. 卫星导航系统通常由______、______和用户终端三部分组成。
答案：卫星段；地面段5. 惯性导航系统误差随时间累积，因此需要______进行校准。
答案：外部参考三、简答题（每题5分，共20分）1. 简述GPS系统的主要组成。
答案：GPS系统主要由空间段、地面控制段和用户设备段三部分组成。
2. 请解释什么是差分GPS（DGPS）？答案：差分GPS是一种通过一个已知位置的基准站发送误差修正信息来提高GPS定位精度的技术。
3. 北斗卫星导航系统相比GPS有哪些优势？答案：北斗卫星导航系统具有短报文通信功能，同时在亚太地区提供更高精度的定位服务。
4. 惯性导航系统在哪些领域有广泛应用？答案：惯性导航系统广泛应用于航空、航天、航海、陆地交通等领域。
导航的工作原理
导航的工作原理
导航是现代科技的产物，它的工作原理是基于卫星定位系统（GPS）和地图数据的结合。
通过GPS技术，导航设备可以准确地确
定车辆、行人或其他移动物体的位置，并且结合地图数据，为用户
提供最佳的行车路线或步行路线。
下面我们将详细介绍导航的工作
原理。
首先，导航设备通过接收来自卫星的信号，确定自身的位置。
卫星定位系统由一系列绕地球轨道运行的卫星组成，它们向地面发
送信号，接收器接收这些信号并计算出自身的位置。
这种技术可以
在全球范围内实现定位，因此导航设备可以在任何地方使用。
其次，导航设备将确定的位置与地图数据进行结合，为用户提
供导航服务。
地图数据包括道路信息、建筑物位置、地形等，它们
被存储在导航设备的内部或通过互联网实时获取。
导航设备根据用
户输入的目的地，结合地图数据计算出最佳的行车路线或步行路线，并在屏幕上显示出来。
此外，导航设备还会根据实时交通信息进行路线优化。
通过接
收交通信息，导航设备可以及时更新路线，避开拥堵的道路，为用
户节省时间和油费。
这种实时更新的功能使得导航设备更加智能、实用。
总的来说，导航的工作原理是基于卫星定位系统和地图数据的结合，通过确定位置、计算路线、实时更新交通信息等功能，为用户提供准确、便捷的导航服务。
随着科技的不断发展，导航设备的功能将会不断增强，为人们的出行带来更多便利和安全保障。
哈尔滨工业大学导航原理大作业报告
Assignment of Principles of Navigation 《导航原理》作业(惯性导航部分2016 秋)My Chinese NameMy Class No.My Student No.Task descriptionIn an fictitious mission, a spaceship is to be lifted from a launching site located at 19o37' NL and 110o57' EL, into a circular orbit 400 kilometers high along the equator. The spaceship is equipped with a strapdown INS whose three gyros, G X , G Y , G Z , and three accelerometers, A X , A Y , A Z , are installed respectively along the axes X b , Y b , Z b of the body frame.Case 1: Stationary testDuring a pre-launching ground test, the body frame of the spaceship initially coincides with the local geographical frame, with its pitching axis X b pointing to the east, rolling axis Y b to the north, and heading axis Z b upwards. Then the body of the spaceship is made to rotate in 3 steps: (1) 80oaround X b(2) 90oaround Y b (3) 170oaround Z bAfter that, the body of the spaceship stops rotating. You are required to compute the final outputs of the three accelerometers in the spaceship, using quaternion and ignoring the device errors. It is assumed that the magnitude of gravity acceleration at the ground level is g 0 = 9.79m/s 2.Case 2: The launching processThe spaceship is installed on the top of an vertically erected rocket. Its initial heading, pitching and rolling angles with respect to the local geographical frame are -90, 90 and 0 degrees respectively. The default rotation sequence is heading → pitching → rolling. The top of the rocket is initially 100m above the sea level. Then the rocket is fired up.The outputs of the gyros and accelerometers in the spaceship are both pulse numbers. Each gyro pulse is an angular increment of 0.01 arcsec, and each accelerometer pulse is 1e -7g 0, with g 0 = 9.79m/s 2. The gyro output fre- quency is 100Hz, and the accelerometer’s is 5Hz.The outputs of the gyros and accelerometers within 1460s are stored in a MATLAB data file named mission.mat, con-taining matrices GGM of 146000×3 from gyros and AAM of 7300×3 from accelerometers respectively. The format of the data in the two matrices is as shown in the tables, with 10 rows of each matrix selected. Each row represents the outputs of the type of sensors at a sampling time.The Earth can be seen as an ideal sphere, with radius R = 6371.00km and spinning rate is 7.292×10-5rad/s, The errors of the gyros and ac- celerometers can be ignored, but the effect ofheight on the magnitude of gravity has to be taken into account. The gravity acceleration at the sea level of the near equator region can be chosen as g 0=9.79 m/s 2, and the magnitude ofgravity acceleration at a height of h can be computed asBesides, the influence of height on the angular rates of the geographical frame and the changing rates oflatitude and longitude should also be considered.The outputs of the gyros can be integrated every 0.01s. The rotation of the geographical frame is to be up- dated every 0.2s, so are the velocities and positions of the spaceship. You are required to:(1) compute the final attitude quaternion, longitude, latitude, height, and east, north, vertical velocities of the spaceship.(2) draw the latitude-versus-longitude trajectory of the spaceship, with horizontal longitude axis. (3) draw the curve of the height of the spaceship, with horizontal time axis. (4) draw the curves of the attitude angles of the spaceship, with horizontal time axis.g h =g 0R 2(R +h )2一 任务一1 利用四元数计算飞船上加速度计的最终输出初始时刻三个加速度计的输出A=[0,0,9.79]。
导航原理测试题及答案
导航原理测试题及答案一、单项选择题（每题2分，共10题）1. 导航系统的基本功能是什么？A. 定位B. 导航C. 定位与导航D. 定位、导航与授时答案：D2. GPS系统由哪三部分组成？A. 空间部分、地面控制部分、用户设备部分B. 空间部分、地面控制部分、卫星通信部分C. 空间部分、地面控制部分、用户定位部分D. 空间部分、地面控制部分、地面接收部分答案：A3. 北斗卫星导航系统属于哪种类型的导航系统？A. 区域性导航系统B. 全球性导航系统C. 区域性与全球性结合的导航系统D. 区域性与局部性结合的导航系统答案：C4. 以下哪个不是导航系统中常用的定位方法？A. 伪距定位B. 载波相位定位C. 多普勒定位D. 惯性定位5. GPS定位的基本原理是什么？A. 通过测量卫星信号的传播时间来确定位置B. 通过测量卫星信号的频率来确定位置C. 通过测量卫星信号的强度来确定位置D. 通过测量卫星信号的相位来确定位置答案：A6. 卫星导航系统能够提供的定位精度通常是多少？A. 米级B. 分米级C. 厘米级D. 毫米级答案：A7. 卫星导航系统中，卫星的轨道高度通常是多少？A. 几百公里B. 几千公里C. 几万公里D. 几十万公里答案：B8. 卫星导航系统中，卫星的轨道倾角通常是多少？A. 0度B. 30度C. 60度D. 90度答案：D9. 卫星导航系统中，卫星的轨道周期通常是多少？B. 一天C. 几个月D. 一年答案：B10. 卫星导航系统中，卫星的轨道面通常是多少？A. 一个B. 两个C. 三个D. 四个答案：C二、多项选择题（每题3分，共5题）1. 以下哪些因素会影响GPS定位精度？A. 卫星的几何分布B. 电离层延迟C. 大气层延迟D. 地面多路径效应答案：ABCD2. 北斗卫星导航系统的主要服务包括哪些？A. 定位服务B. 导航服务C. 授时服务D. 短报文通信服务答案：ABCD3. 卫星导航系统中，哪些是用户设备需要解决的问题？A. 卫星信号的接收B. 卫星信号的解码C. 卫星信号的定位D. 卫星信号的增强答案：ABC4. 卫星导航系统中，哪些是地面控制部分的主要功能？A. 卫星轨道的监测B. 卫星时钟的校准C. 卫星信号的广播D. 用户定位的计算答案：ABC5. 卫星导航系统中，哪些是空间部分的主要组成？A. 导航卫星B. 地面控制站C. 用户接收机D. 地面监测站答案：AD三、判断题（每题1分，共5题）1. GPS系统可以提供全球范围内的定位服务。
哈工大-导航原理作业
Assignments of Inertial Navigation 《惯性导航》作业Autumn 2017Assignment 1: Coordinate transformation for 3D animationfigure 1.1 Interface for rotating animation control of missileAttached is a group of MATLAB programs for 3D animation control of a missile. Initially, the body frame of the missile coincides with the local geographical frame, with its pitching axis X pointing to the east, rolling axis Y to the north, and heading axis Z to the sky, as shown in figure 1.1 which is the controlling in-terface produced by running the program main.m.The three coordinates of each vertex of the 3D model occupies a row in the matrix VTX, The color property and surface information of the model are defined by matrices VTXcolor and faceM respectively. These matrices are loaded from the data file missiledata.mat.Input an angle (in degree) in the "Rotation Angle" box, and click one of the rotation buttons ("Heading Ro-tation", "Pitching Rotation", or "Rolling Rotation"), then it is expected that the missile will rotate from its current attitude for the input angle around the chosen axis of its current body.The animation is to be achieved by the program redraw.m which redraws the missile every once it rotates for one degree, using patch command. During each step of the animation, the current rotation angle of the missile relative to its pre-animation attitude has been generated and stored in one of the variables head, pitch and roll, with only one of them being nonzero. Before each redrawing, you need to re-calculate the coordinates of the missile in VTX resolved in the geographical frame, instead of keeping them the same as their values at the very beginning (VTX0). That is, in the program redraw.m, you need to replace the command line VTX=VTX0with your own codes for re-calculating VTX. Elsewhere, additional codes also might be required to make it work.Please rewrite the program redraw.m, or others if necessary, so that successive rotating animations can be achieved, and explain the rationale behind your rewriting.一、 任务分析本作业需要实现的功能为实现火箭的连续旋转。
导航原理考试题及答案
导航原理考试题及答案一、选择题（每题2分，共20分）1. 卫星导航系统主要利用哪种信号来确定位置？A. 无线电信号B. 声波信号C. 光波信号D. 电磁波信号答案：A2. GPS系统的全称是什么？A. Global Positioning SystemB. General Positioning SystemC. Global Positioning ServiceD. Global Positioning Setup答案：A3. 以下哪个不是卫星导航系统的组成部分？A. 卫星B. 地面监控站C. 接收机D. 天线答案：D4. 卫星导航系统通常至少需要多少颗卫星来确定一个位置？A. 1颗B. 2颗C. 3颗D. 4颗答案：C5. 卫星导航的误差来源不包括以下哪项？A. 大气延迟B. 多路径效应C. 卫星时钟误差D. 地面建筑物答案：D二、填空题（每空1分，共10分）1. 卫星导航系统通过测量卫星与接收机之间的_________来确定位置。
答案：距离2. GPS系统由_________、卫星和用户接收机三部分组成。
答案：地面监控站3. 卫星导航系统的时间同步是通过_________来实现的。
答案：原子钟4. 卫星导航系统中，_________是影响定位精度的主要因素之一。
答案：卫星分布5. 卫星导航系统在室内或城市峡谷中可能受到_________的影响，导致定位不准确。
答案：多路径效应三、简答题（每题10分，共20分）1. 简述卫星导航系统的工作原理。
答案：卫星导航系统通过接收来自至少三颗卫星的信号，利用信号传播的时间和卫星的位置信息，计算出接收机的三维坐标。
接收机接收到的信号包括卫星的轨道参数和发射时间戳，通过计算信号传播的时间差，结合卫星的精确位置，可以解算出接收机的位置。
2. 描述卫星导航系统在军事领域的应用。
答案：卫星导航系统在军事领域有广泛的应用，包括但不限于精确制导、战场定位、部队调度、无人机导航等。
哈工大导航原理大作业
《导航原理》作业(惯性导航部分)一、题目要求A fighter equipped with SINS is initially at the position of ︒35 NL and︒122 EL,stationary on a motionless carrier. Threegyros X G ,Y G ,Z G ,and three accelerometers, X A ,Y A ,Z A are installed along the axes b X ,b Y ,b Z of the body frame respectively. Case 1:stationary onboard testThe body frame of the fighter initially coincides with the geographical frame, as shown in the figure, with its pitching axis b X pointing to the east,rolling axis b Y to the north, and azimuth axis b Z upward. Then the body of the fighter is made to rotate step by step relative to the geographical frame.(1) ︒10around b X (2) ︒30around b Y (3) ︒50-around b ZAfter that, the body of the fighter stops rotating. You are required to compute the final output of the three accelerometers on the fighter, using both DCM and quaternion respectively,and ignoring the device errors. It is known that the magnitude of gravity acceleration is 2/8.9g s m =. Case 2:flight navigationInitially, the fighter is stationary on the motionless carrier with its board 25m above the sea level. Its pitching and rolling axes are both in the local horizon, and its rolling axis is ︒45on the north by east, parallel with the runway onboard. Then the fighter accelerate along therunway and take off from the carrier.The output of the gyros and accelerometers are both pulse numbers,Each gyro pulse is an angular increment of sec arc 1.0-,and each accelerometer pulse is g 6e 1-,with 2/8.9g s m =.The gyro output frequency is 10 Hz,and the accelerometer ’s is 1Hz. The output of gyros and accelerometers within 5400s are stored in MATLAB data files named and , containing matrices gm of 35400⨯ and am of 35400⨯ respectively. The format of data as shown in the tables, with 10 rows of each matrix selected. Each row represents the out of the type of sensors at each sample time.The Earth can be seen as an ideal sphere, with radius and spinning raterad/s 10292.75-⨯, Theerrors of sensors are ignored, so is theeffect of height on the magnitude of gravity. The output of the gyros are to integrated every . The rotation of geographical frame is to be updated every 1s, so are the velocities and position of the figure. You are required to:(1)Compute the final attitude quaternion, longitude, latitude, height, and east, north, vertical velocities of the fighter.(2)Compute the total horizontal distance traveled by the fighter. (3)Draw the latitude-versus-longitude trajectory of the fighter, with horizontal longitude axis.(4)Draw the curve of the height of fighter, with horizontal time axis.二、Case1解答方向余弦阵法（1） 绕Xb 轴转过ψ︒=10ϕ⎪⎪⎪⎭⎫ ⎝⎛︒︒-︒︒=⎪⎪⎪⎭⎫ ⎝⎛-=10cos 10sin 010sin 10cos 0001cos sin 0sin cos 0001ϕϕϕϕϕC（2） 绕Yb 轴转过︒=30θ⎪⎪⎪⎭⎫ ⎝⎛︒︒︒-︒=⎪⎪⎪⎭⎫ ⎝⎛-=30cos 030sin 01030sin 030cos cos 0sin 010sin 0cos θθθθθC（3） 绕Zb 轴转过︒-=50ψ⎪⎪⎪⎭⎫ ⎝⎛︒-︒-︒-︒-=⎪⎪⎪⎭⎫ ⎝⎛=1000）50（cos ）50（sin -0）50（sin ）50（cos 1000cos sin -0sin cos ψψψψψC 所以变换后的坐标X E Y C C C N Z ϕθψζ⎛⎫⎛⎫ ⎪⎪= ⎪ ⎪ ⎪ ⎪⎝⎭⎝⎭由于初始时刻有009.8E N ζ⎛⎫⎛⎫ ⎪ ⎪= ⎪ ⎪ ⎪ ⎪⎝⎭⎝⎭所以计算得三个加速度计的输出分别是 ,/4504.42x s m A -=,/6027.22y s m A -=s m A /3581.8z =2计算程序见附录一四元数法（1）绕Xb 轴转过ψ︒=10ϕ i 2sin2cos q ϕϕϕ+=（2）绕Yb 轴转过︒=30θj 2sin2cosq θθθ+=（3）绕Zb 轴转过︒-=50ψk 2sin2cosq ψψψ+=则合成四元数)2sin 2(cos )2sin 2(cos )2sin 2(cosk j i q q q q ⋅+⋅+⋅+==ψψθθϕϕψθϕ合成四元数的逆1123q p i p j p k λ-=---由公式 q N E qZ YX ⎪⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎭⎫⎝⎛-ξ1计算得三个加速度计的输出分别是 ,/4504.42x s m A -=,/6027.22y s m A -=s m A /3581.8z =2由两种计算方法的计算结果可以看出，方向余弦阵法和四元数法的计算结果是一致的。
导航原理大作业
导航原理作业（惯性导航部分）作业内容需包括：1. 题目内容（即本页）2. 程序设计说明及代码3. 计算或仿真结果及分析4. 仿真中遇到的问题或体会一枚导弹采用捷联惯性导航系统，三个速率陀螺仪Gx, Gy, Gz 和三个加速度计 Ax, Ay, Az 的敏感轴分别沿着着弹体坐标系的Xb, Yb, Zb 轴。
初始时刻该导弹处在北纬45.75度，东经126.63度。
第一种情形：正对导弹进行地面静态测试（导弹质心相对地面静止）。
初始时刻弹体坐标系和地理坐标系重合，如图所示，弹体的Xb 轴指东，Yb 轴指北，Zb 轴指天。
此后弹体坐标系 Xb-Yb-Zb 相对地理坐标系的转动如下：首先，弹体绕 Zb （方位轴）转过 -10 度； 接着，弹体绕 Xb （俯仰轴）转过 15 度； 然后，弹体绕 Yb （滚动轴）转过 20 度； 最后弹体相对地面停止旋转。
请分别用方向余弦矩阵和四元数两种方法计算：弹体经过三次旋转并停止之后，弹体上三个加速度计Ax, Ay, Az 的输出。
取重力加速度的大小g = 9.8m/s 2。
第二种情形：导弹正在飞行中。
初始时刻弹体坐标系仍和地理坐标系重合；且导弹初始高度200m ，初始北向速度 1800 m/s ，初始东向速度和垂直速度都为零。
陀螺仪和加速度计的输出都为脉冲数形式，陀螺输出的每个脉冲代表 0.00001弧度的角增量。
加速度计输出的每个脉冲代表1μg ，1g = 9.8m/s 2。
陀螺仪和加速度计输出的采样频率都为10Hz ，在200秒内三个陀螺仪和三个加速度计的输出存在了数据文件gaout.mat 中，内含一矩阵变量ga ，有2000行，6列。
每一行中的数据代表每个采样时刻三个陀螺Gx, Gy, Gz 和三个加速度计Ax, Ay, Az将地球视为理想的球体，半径 6371.00公里，且不考虑仪表误差，也不考虑弹体高度对重力加速度的影响。
选取弹体的姿态计算周期为0.1秒，速度和位置的计算周期为1秒。
导航系统大作业
导航系统1．简述捷联惯性系统中地理系到机体系的姿态阵bg C 其含义及其功能。
答:含义:导航坐标系g g g O x y z -到机体坐标系b b b O x y z -的一组欧拉角为,,θγψ，导航坐标系经过3次转动到机体坐标系。
g g g x y z 依次沿g O z -、'b O x -、''b O y -旋转角度-ψ、θ、γ后到b b b x y z 。
姿态矩阵中包含了机体的姿态角方位角ψ、俯仰角θ和横滚角γ。
功能：机体陀螺仪输出的角速度信息经过补偿后,积分得到机体坐标系与导航坐标系的姿态信息和姿态转移矩阵。
捷联惯导系统中,加速度计与载体固连,利用姿态阵完成加速度计输出信息从机体坐标到导航坐标的转换。
转换后的加速度计信息经过积分可得到机体在导航坐标系下的速度和位置。
2．画出并用式表达速度三角形（地速、控速、风速)及航迹角、航向角与偏流角之间的关系.答：风速：空气相对于地面的运动速度;空速:飞机相对于空气运动的速度；地速：飞机相对于地面的运动速度。
=+v v v 风地空航向角：机头在水平面投影与真北方向的夹角ϕ;偏流角:空速矢量和地速矢量之间的夹角,用δ表示；航迹角:飞机速度矢量在水平面投影与真北方向的夹角。
航向角ϕ加上偏流角δ等于地速v 地的方位角α。
v 地v 空v 风3．简述惯性导航系统、卫星导航系统、多普勒导航、塔康、VOR/DME 、天文导航其各自的基本工作原理、特点及误差特性。
答:一、惯性导航系统（1）工作原理以牛顿力学定律为基础，以陀螺仪和加速度计为敏感器件进行导航参数解算。
系统根据陀螺仪的输出建立导航坐标系，根据加速度计输出解算出运载体的速度和位置，从而实现姿态和航向解算. （2）特点惯性导航系统不需要任何外来信息，也不会向外辐射任何信息，仅依靠惯性器件就能全天候，全球性的自主三维定位和三维定向,同时具备自主性、隐蔽性和信息的完备性。
（3）误差特性误差随时间积累，短时间导航精度较高。
组合导航大作业
基于组合导航系统的数值积分粒子滤波算法摘要本研究工作中，我们研究数值积分粒子滤波器（CPF）算法来计算GPS / INS组合导航系统的估值。
GPS/ INS组合导航系统的误差模型是非线性的。
CPF算法是建立在数值积分卡尔曼滤波器（CKF）和粒子滤波器（PF）之上的，并且具有两者的优点。
因此CPF可以为高维非线性滤波器问题提供一个系统的解决方案。
CPF通过模拟的方式来展现。
模拟结果表明，当与次优技术例如数值积分卡尔曼滤波器（CKF）比较时，这种方法具有优越的性能，因为在这种情况下CKF有大的初始偏差。
模拟的结果表明，该改进的CPF的表现超越了传统的非线性滤器。
该研究为工程设计和改进提供了支持。
关键字：数值积分规则数值积分卡尔曼滤波（CKF）数值积分粒子滤波器（CPF）组合导航1、引言滤波器在实际系统起着非常重要的作用。
在现实世界中，几乎所有的系统都具有非线性特性。
只有深刻领会非线性系统的本质，合理建立系统模型和非线性网络滤波方法才能有效地帮助分析和解决各种在工程实践中遇到的问题。
卡尔曼滤波器（KF）是用于估算线性系统[1]的最优滤波器。
该GPS / INS导航系统通常采用卡尔曼滤波器来估计所述系统的状态[2，3]。
卡尔曼滤波器的简单计算，递归结构和数学严谨的推导，使其适用和吸引大量实际应用的使用。
然而，许多现实世界的系统都是非线性的。
扩展卡尔曼滤波器（EKF）的开发来帮助这些非线性系统。
但是扩展卡尔曼滤波器是一种次优的非线性滤波器，由于当是线性系统[4，5]时截断了高阶项。
EKF的计算时间类似于卡尔曼滤波器[6]的。
数值积分卡尔曼滤波器（CKF）是没有非线性模型[7]的直链化非线性滤波连接的方法。
在CKF算法是公正和最小方差，这比GPS / INS 动态系统[8] EKF方法更好。
在CKF中被建议使用贬低线性偏置的非线性测量方程[9]。
物理系统往往受到意想不到的偏差或失误[10]。
其结果是，有一个方法，来维持一个精确的和可靠的解决方案[11]是重要的。
北航惯导第一次大作业
《惯性导航原理》第一次大作业一、 原理分析惯导系统为指北方位的平台系统，则利用比力方程以及陀螺提供的东、北、天三个比力数据，即可计算得到在每个数据采集点的平台即时速度，再通过经纬度的计算公式，就可以得到每个数据采集点平台的即时经纬度，以每个数据采集点为下一个采集点的起点，即可对速度和经纬度进行累计计算，从而得到平台在运动过程中任意时刻的速度和位置情况。
运动过程中任意时刻的速度和位置情况。
1.模型公式的推导载体相对地球运动时，载体相对地球运动时，加速度计测得的比力表达式，加速度计测得的比力表达式，加速度计测得的比力表达式，称为比力方程，称为比力方程，称为比力方程，方程如方程如下：下：g V V f epep ieep-´++=)2(vv （1）在指北方案中，平台模拟地理坐标系，将上式中平台坐标系用地理坐标系代入得：入得：t tt ett iettgV f V+´+-=)2(v v （2）系统中测量的是比力分量，将上式写成分量形式系统中测量的是比力分量，将上式写成分量形式=-+ （3） 又因为地球的自转角速率为：又因为地球的自转角速率为：(4)地理坐标系相对于地球坐标系的角速率为：地理坐标系相对于地球坐标系的角速率为：= (5)将（将（44）（5）两个式子带入（）两个式子带入（33）式，即可得到如下方程组：）式，即可得到如下方程组：（6）2.速度计算作业要求只考虑水平通道，作业要求只考虑水平通道，因此只需要计算正东、因此只需要计算正东、因此只需要计算正东、正北两个方向的速度即可。
正北两个方向的速度即可。
正北两个方向的速度即可。
理理论上计算得到t x V 、t y V 后，再积分一次可得到速度值，即后，再积分一次可得到速度值，即ïîïíì+=+=òòt t y t y t ytt x t x tx V dt V V V dt V V 000但在本次计算过程中，三个方向的速度均是从零开始在各时间节点上的累加，并不是t的函数，因此速度计算可以由以下方程组实现：（7）此方程组表示了从第i 个采集点到第（个采集点到第（i+1i+1i+1）个采集点的速度递推公式。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
导航原理实验报告——基于Matlab的地心坐标系到大地坐标系的算法原理探测制导与控制11014201110127雷继松一、算法描述地心坐标系（简称“ECEF ”）到大地坐标系（简称“GC ”）的转换原理主要采用牛顿迭代法，采用牛顿迭代法的优点是高效、高精度以及在距离地心小范围内不受“奇点”（空间-时间上具有无限曲率的一点）及“不收敛”两大因素的影响。
1、所用到的参数R----地球半长轴r----地球半短轴f----地球椭率L----大地纬度λ----大地经度H----相对参考椭球体的高度2、算法给定坐标值Xe,Ye,Ze,椭球参数R ，r,f,迭代精度δ，最大迭代次数Kmax;（1）开始；（2）计算x0，z0和系数A ，B ，C ；).(2,2,)()()(;,;sin sin sin ;sin ;cos ));1/()cot((cos ;cos cos cos 220034022011r R C Rx B rz A A t C B t C B At t f Z z Y X x L h r L h z Z M Y M X f L a L h R L h x M e e e e e e −===−−+++==+=+=+===−=+=+=θλλθθ（3）若z0=0，设L=0，h=x0-R,执行第十步，否则执行第4步；（4）为迭代计算初值t0,初始化k=0;2000000])1[(1)sgn()1(z xf z z x f t −++−−=（5）令k=k+1;（6）根据迭代原理计算t （k ）；)()(34)(,...2,1,0),(/)(23''1C B t C B At t f k t f t f t t k k k k −+++==−=+（7）如果|t(k)-t(k-1)|δ≤,执行步骤9，否则执行步骤8；（8）若k<kmax,执行步骤5，否则算法不收敛，执行步骤12；（9）计算纬度L ；]21)1cot[(2t t f a L −−=（10）计算经度λ；),(2atan e e X Y =λ（11）计算高度h ；2220]21)1[(1)12)(sgn(t t f t t r z L h −−++−=（12）结束算法。
18000个测试点：对于高度h ，从-6300km 到30000km 取20个测试点；对于纬度L ，从-90度到90度取45个测试点；对于经度λ，从-180度到180度取30个测试点；精度δ这里取10e-9，这里计算产生的误差成为欧几里德误差，即给定值和计算值之间的误差；实验所要求的图3得到的坐标转换下迭代数的图形，图4则是相应的欧几里德误差图。
对于所要求实验图5和6，精度δ取2×10e-8。
由于此算法大多数应用都是在地球表面，所以高度h 取-50km 到1000km ，也是取20个点，其他不变，则可以得到图5和图6。
最后一个图是评估算法在整个空间的一个效果，仿真选择接近地心的的测试点进行。
高度h 选择-6378km 到-6300km ，其他不变，和前面一样。
二、代码实现及图形（很多说明直接在代码中描述）本程序是由三个M 文件组成，图形绘制主要由第三个M 文件，由于编写代码时候很多代码相同，所以采用Ctrl+R 和Ctrl+T 的方法注释掉和恢复来进行，当然下面说明时并没有这样，已经是分开来论述，这里要注意一下——第一个M文件：%achieve Newton-Raphson(NR)method;from ECEF to Geodetic Coordinates function[L,nanbuda,h,k]=from_ecef_to_geodetic(xe,ye,ze)%这里编程出现一个问题，为了防止xe,ye,ze的重复定义，这里分别用两个函数来实现；x0=sqrt(xe^2+ye^2);z0=ze;r=6356.7523142;R=6378.137;f=(R-r)/R;iter_accuracy=10^(-9);%画前两个图时的迭代精度%iter_accuracy=2*10^(-8);%画第三个图到最后一个图的迭代精度kmax=20;A=r*z0;B=2*R*x0;C=2*(R^2-r^2);if(z0==0)L=0;h=x0-R;k=0;elsek=1;t0=-(1-f)*x0/z0+sign(z0)*sqrt(1+((1-f)*x0/z0)^2);t1=t0-(A*t0^4+(B+C)*t0^3+(B-C)*t0-A)/(4*A*t0^3+3*(B+C)*t0^2+(B-C)) ;%这里注意matlab定义数组下标是从1开始的，C语言是从0开始的，注意！while(abs(t1-t0)>iter_accuracy)t0=t1;t1=t0-(A*t0^4+(B+C)*t0^3+(B-C)*t0-A)/(4*A*t0^3+3*(B+C)*t0^2+(B-C)) ;k=k+1;if k>=kmax,break;endendL=acot((1-f)*(1-t1^2)/2/t1);h=sign(L)*(z0-r*2*t1/(1+t1^2))*sqrt(1+((1-f)*(1-t1^2)/2/t1)^2); endnanbuda=atan2(ye,xe);第二个M文件function[xe,ye,ze]=transformation(L,nanbuda,h)r=6356.7523142;R=6378.137;f=(R-r)/R;sita=acot(cot(L)/(1-f));M=R*cos(sita)+h*cos(L);xe=M*cos(nanbuda);ye=M*sin(nanbuda);ze=r*sin(sita)+h*sin(L);第三个M文件：h=linspace(-6300,-30000,20);%这里使用linspace，表达在-6300到-30000之间取出20个点；L=linspace(-pi/2,pi/2,45);%这里考虑到sin,cos,tan等的运算，所以把角度转换成弧度形式；nanbuda=linspace(-pi,pi,30);%由于采用矩阵或者向量的形式处理多个点的话，考虑矩阵维数一定要符合基本运算，太复杂，%所以这里决心采用多次循环迭代的方法；%由于我们绘制的图形是空间三维图形，所以我们建立点的时候就考虑三维变量好了；%h=linspace(-50,1000,20);%L=linspace(-pi/2,pi/2,45);%这里考虑到sin,cos,tan等的运算，所以把角度转换成弧度形式；%nanbuda=linspace(-pi,pi,30);%h=linspace(-6378,-6300,20);%L=linspace(-pi/2,pi/2,45);%nanbuda=linspace(-pi,pi,30);R=6378137;r=6356752.3142;k1=zeros(20,45,30);%定义数组来存放上面h、L、nanbuda的值，方便对点进行迭代运算xe1=zeros(20,45,30);ye1=zeros(20,45,30);ze1=zeros(20,45,30);xe2=zeros(20,45,30);%下面三个牛顿迭代后的值的存储数组ye2=zeros(20,45,30);ze2=zeros(20,45,30);eucliddistance=zeros(20,45,30);%存放欧几里德误差的数组for hi=1:20for Li=1:45for nanbudai=1:30%利用三个for循环来实现迭代，相当于是x,y,z三个维度；[xe1(hi,Li,nanbudai),ye1(hi,Li,nanbudai),ze1(hi,Li,nanbudai)]=tra nsformation(L(Li),nanbuda(nanbudai),h(hi));[L2(hi,Li,nanbudai),nanbuda2(hi,Li,nanbudai),h2(hi,Li,nanbudai),k1(hi,Li,nanbudai)]=from_ecef_to_geodetic(xe1(hi,Li,nanbudai),...ye1(hi,Li,nanbudai),ze1(hi,Li,nanbudai));[xe2(hi,Li,nanbudai),ye2(hi,Li,nanbudai),ze2(hi,Li,nanbudai)]=tra nsformation(L2(hi,Li,nanbudai),nanbuda2(hi,Li,nanbudai),h2(hi,Li, nanbudai));eucliddistance(hi,Li,nanbudai)=...%计算欧几里德误差（相当于是给定值和计算值之间的差值）sqrt((xe1(hi,Li,nanbudai)-xe2(hi,Li,nanbudai))^2+... %这里我们前面说，看作各个点的运算，所以(ye1(hi,Li,nanbudai)-ye2(hi,Li,nanbudai))^2+...%误差也就是给定点和计算值点两点之间的距离(ze1(hi,Li,nanbudai)-ze2(hi,Li,nanbudai))^2);%所以采用算数平方根来计算endendendfigurez=k1(:,:,1);[x,y]=meshgrid(h,L);%[X,Y]=meshgrid(x,y)将向量x和y定义的区域转换成矩阵X和Y，%这两个矩阵可以用来表示mesh和surf的三维空间点以及两个变量的赋值。
%其中矩阵X的行向量是向量x的简单复制，而矩阵Y的列向量是向量y的简单复制。
surf(x,y,z')%绘制空间三维图形[x,y]=meshgrid(h,L);z=eucliddistance(:,:,1);figure;surf(x,y,z')-34迭代数1-34-202468欧几里德误差1紧接着绘制第三个图和第四个图时，把第一个M文件的精度改为iter_accuracy=2*10^(-8);第三个M文件的前三行代码改为（即h、L、nanbuda的参数），其他不变：h=linspace(-50,1000,20);L=linspace(-pi/2,pi/2,45);%这里考虑到sin,cos,tan等的运算，所以把角度转换成弧度形式；nanbuda=linspace(-pi,pi,30);绘制图形如下：1000迭代数21000欧几里德误差2最后一个图：绘图的时候要在前两个基础上改动的除了前三行的定义以外，还要继续改变画图的程序语句：其余一样：h=linspace(-6378,-6300,20);L=linspace(-pi/2,pi/2,45);nanbuda=linspace(-pi,pi,30); ...........................................z=k1(:,:,1);[x,y]=meshgrid(h,L);figure;plot3(x,y,z','+')%由于surf绘制绘制曲面图，这里绘制离散的点图我们使用plot3Grid绘制结束后，和实验理论值发生很大的改变，其中理论的算法编程可能出现逻辑问题；-6300不收敛点图三、实验结论和心得此算法和以前的算法相比较，精确性和效率性都有显著提高。