利用MATLAB编制的GPS单点定位程序
基于M ATLAB 的BDS单点定位程序设计及精度分析
基于M ATLAB 的BDS单点定位程序设计及精度分析作者:王皓来源:《科技创新与生产力》 2018年第11期摘要:本文系统地研究了北斗卫星导航系统(BDS)伪距单点定位的相关理论与方法,利用误差改正后的观测方程进行伪距单点定位,并对电离层误差、对流层误差、多径效应误差、卫星星历误差等进行改正。
利用MATLAB软件编写数据处理程序,得到xls和txt格式的坐标、残差数据、残差图以及卫星的PDOP图、GDOP图、HDOP图和VDOP图。
使用上海佘山国际GPS服务(IGS)观测站的数据进行精度实验,对相关数据进行了计算,得到了有关差值数据的相关统计,并结合中误差以及残差情况分析计算结果。
计算结果表明,该数据处理程序在X方向上的定位精度在15 m以内,在Y方向和Z方向上的定位精度都在10 m以内。
关键词:北斗卫星导航系统;BDS;伪距单点定位;定位误差;误差改正;MATLAB中图分类号:P228.4;O241.1 文献标志码:A DOI:10.3969/j.issn.1674-9146.2018.11.075全球导航卫星系统(Global Navigation Satellite System,GNSS)可以为陆地和海上用户提供全球、全天候、高精度的测距、定速和定时服务,并已成为人类获取位置和时间信息的重要手段。
GNSS在军事、社会与经济发展中都发挥着重要的作用,各大国都争先恐后发展独立自主的导航卫星系统[1]。
中国也紧跟美国、俄罗斯和欧盟的步伐,正在建设中国自己的卫星导航系统——北斗卫星导航系统(BeiDou navigation satellite System,BDS)。
长期以来,GNSS的导航与定位主要依赖于全球定位系统(Global Positioning System,GPS),有研究指出其定位精度可以实现静态厘米级、动态分米级。
BDS的建立将打破GPS在我国的统治地位[2],基于BDS的测量试验有待于更深层次的开展,区域系统的单点定位应用的相应技术指标还没有得到全面、系统性的论证。
MATLAB技术GPS数据处理教程
MATLAB技术GPS数据处理教程导言GPS(全球定位系统)以其高精度的定位功能在各个领域得到广泛应用。
在我们的日常生活中,GPS定位已经成为了司空见惯的技术,无论是导航系统、运输管理还是环境监测等。
本文旨在介绍使用MATLAB技术处理GPS数据的方法和技巧,帮助读者更好地利用现有的GPS数据进行分析和应用。
一、GPS数据处理的基本概念在开始介绍MATLAB技术处理GPS数据之前,我们首先需要了解一些GPS数据处理的基本概念。
GPS接收器通过接收卫星发射的信号来计算接收器与卫星之间的距离,进而确定接收器的位置。
GPS数据通常包括接收器的经纬度坐标、时间戳和位置精度等信息。
二、数据准备与读取要处理GPS数据,首先需要准备好相应的数据文件。
通常,GPS数据会以文本文件(如CSV或TXT)的形式进行存储。
在MATLAB中,可以使用load或readtable函数来读取存储在文本文件中的GPS数据,并将其转换为MATLAB中的矩阵或表格形式,方便进行后续的处理。
三、数据预处理与清洗在进行GPS数据处理之前,常常需要对数据进行预处理和清洗,以去除异常值和噪音,确保数据的准确性。
MATLAB提供了一系列的函数和工具箱,例如数据滤波技术、异常值检测算法等,可以用于数据的预处理和清洗。
此外,还可以利用MATLAB内置的绘图工具对数据进行可视化,快速发现异常数据。
四、数据分析与可视化数据处理的核心在于数据分析和可视化。
MATLAB提供了丰富的函数和工具箱,可以对GPS数据进行多维度的分析和可视化。
例如,可以利用MATLAB的统计工具箱对数据进行统计分析,找出数据中的趋势和规律。
同时,也可以使用MATLAB的绘图函数绘制位置轨迹图、速度图以及加速度图等,直观地展示数据的空间和时序特征。
五、轨迹重建与预测GPS数据处理的另一个重要应用是轨迹重建和预测。
基于历史GPS数据,我们可以通过建立合适的模型,使用MATLAB的预测函数对未来的轨迹进行预测。
matlab单点定位代码
matlab单点定位代码在MATLAB中进行单点定位,通常涉及到的是使用某种形式的信号处理技术,如最小二乘法、卡尔曼滤波等,来估计信号源的位置。
下面是一个简单的示例,演示了如何使用最小二乘法进行单点定位。
假设我们有一个已知位置的接收器数组,以及接收到的信号强度(或其他测量值)。
我们可以使用这些信息来估计信号源的位置。
以下是一个示例代码:matlab复制代码% 假设接收器坐标和测量值rx_coords = [1, 2; 3, 4; 5, 6]; % 接收器坐标(x, y)measured_values = [30; 25; 20]; % 接收到的信号强度或测量值% 计算接收器之间的距离distances = pdist2(rx_coords(:,1), rx_coords(:,1), 'euclidean');% 创建最小二乘问题X = [measured_values; -1*measured_values];Y = distances(:);% 最小二乘求解% 这里假设我们已经知道测量值之间的线性关系为 Y = X*beta% 因此我们使用pinv函数求解最小二乘问题,找到系数betabeta = pinv(X'*X) * X' * Y;% 使用得到的系数来估计信号源位置estimated_source_x = beta(1);estimated_source_y = beta(2);% 输出估计的位置fprintf('Estimated source position: (%f, %f)\n', estimated_source_x,estimated_source_y);这个代码只是一个简单示例,实际上信号传播和接收可能受到许多因素的影响,例如环境中的障碍物、多径效应等。
这些因素可能需要更复杂的模型和算法来处理。
此外,对于更复杂的情况,可能需要使用更高级的定位技术,如多点定位、指纹地图等。
【精品】用MATLAB计算GPS卫星位置-最新文档资料
用M A T L A B计算G P S 卫星位置-最新文档资料用MATLAB计算GPS卫星位置GPS定位的基本原理简单来说就是在WGS-84空间直角坐标系中,确定未知点与GPS卫星的空间几何关系。
因此利用GPS 进行导航和测量时,卫星是作为位置已知的高空观测目标。
那么如何精确快速的解算出卫星在空间运行的轨迹即其轨道是实现未知点快速定位的关键。
1 标准格式RINEX格式简述在进行GPS数据处理时,由于接收机出自于不同厂家,所以厂家设计的数据格式也是五花八门的,但是在实际中,很多时候需要把来自不同型号的接收机的数据放在一块进行处理,这就需要数据格式的统一,为了解决这种矛盾,RINEX(英文全称为:The Receiver Independent Exchange Format)格式则应运而生,该格式存储数据的类型是文本文件,数据记录格式是独立于接收机的出自厂家和具体型号的。
由此可见,其特点是:由于是通用格式,所以可将不同型号接收机收集的数据进行统一处理,并且大多数大型数据处理软件都能够识别处理,此外也适用于多种型号的接收机联合作业,通用性很强。
RINEX标准文件里不是单一的一个文件,而是包括如下几种类型的文件[1]。
(1)观测数据文件(ssssdddf.yyo),记录的是GPS观测值信息,(OBServation data,简写OBS,为接收机记录的伪距、相位观测值;O文件,如XG012191.10O)。
(2)导航电文文件(ssssdddf.yyn),记录的是GPS卫星星历信息(NAVavigation data,简写NAV,记录实时发布的广播星历;N文件,如XG012191.10N)。
(3)气象数据文件(ssssdddf.yym),主要是在测站处所测定的气象数据(METerological data,简写MET,记录气象仪器观测的温、压、湿度状况;M文件,如XG012191.10M)。
(4)GLONASS导航电文文件(ssssdddf.yyg),记录的是地球同步卫星的导航电文。
GPS卫星运动及定位matlab仿真
GPS卫星运动及定位matlab仿真摘要全球定位系统是具有全球性、全能性、全天候优势的导航定位、定时和测速系统,现在在全球很多领域获得了应用。
GPS卫星的定位是一个比较复杂的系统,其包含参数众多,如时间系统、空间坐标系统等。
此次设计是针对卫星运动定位的matlab仿真实现,因要求不高,所以对卫星运动做了理想化处理,摄动力对卫星的影响忽略不计(所以为无摄运动),采用开普勒定律及最小二乘法计算其轨道参数,对其运动规律进行简略分析,并使用matlab编程仿真实现了卫星的运功轨道平面、运动动态、可见卫星的分布及利用可见卫星计算出用户位置。
通过此次设计,对于GPS卫星有了初步的认识,对于静态单点定位、伪距等相关概念有一定了解。
关键字:GPS卫星无摄运动伪距matlab仿真The movement and location of GPS satellite onMA TLABAbstract:Global positioning system is a global, versatility, all-weather advantage of navigation and positioning, timing and speed system, now there has many application in many fields.GPS satellite positioning is a complex system, which includes many parameters, such as time and space coordinates system.This design is based on the matlab simulation of satellite motion and location, because demand is not high, so to do the idealized satellite movement, and ignore the disturbed motion ( so call it non-disturbed motion ).Using the Kepler and least-square method for calculating the parameters of orbital motion, for the characteristics of motion to make a simple analysis, and use the matlab simulation to program achieve the orbital plane of satellite, the dynamic motion, the distribution of visible satellites and using visible satellites to calculate the users‟ home.Through the design have primary understanding for the GPS satellite, and understanding the static single-point, pseudorange and so on.Key words:GPS satellite non-disturbed motion pseudorange matlab simulation目录第一章前言 (4)1.1课题背景 (4)1.2本课题研究的意义和方法 (5)1.3GPS前景 (5)第二章 GPS测量原理 (7)2.1伪距测量的原理 (7)2.1.1 计算卫星位置 (8)2.1.2 用户位置的计算 (8)2.1.3 最小二乘法介绍 (8)2.2载波相位测量原理 (9)第三章 GPS的坐标、时间系统 (13)3.1坐标系统 (13)3.1.1 天球坐标系 (13)3.1.2 地球坐标系 (15)3.2时间系统 (16)3.2.1 世界时系统 (17)3.2.2 原子时系统 (18)3.2.3动力学时系统 (19)3.2.4协调世界时 (19)3.2.5 GPS时间系统 (19)第四章卫星运动基本定律及其求解 (21)4.1开普勒第一定律 (21)4.2开普勒第二定律 (22)4.3开普勒第三定律 (23)4.4卫星的无摄运动参数 (23)4.5真近点角的概念及其求解 (24)4.6卫星瞬时位置的求解 (25)第五章 GPS的MATLAB仿真 (28)5.1卫星可见性的估算 (28)5.2GPS卫星运动的MATLAB仿真 (29)结论 (41)致谢 (43)参考文献 (44)附录.............................................................................................. 错误!未定义书签。
matlab伪距单点定位
matlab伪距单点定位Matlab伪距单点定位伪距单点定位是一种利用卫星信号进行定位的方法,通过测量卫星信号的传播时间差来计算接收器与卫星之间的距离,并利用多个卫星的距离信息进行定位。
Matlab作为一种强大的数学计算工具,可以方便地实现伪距单点定位算法。
伪距单点定位的原理是利用接收器接收到的卫星信号的传播时间差来计算接收器与卫星之间的距离。
接收器通过测量卫星信号的到达时间差来计算伪距,然后利用伪距信息进行定位。
伪距是接收器接收到卫星信号的传播时间与光速之间的乘积,即伪距=传播时间×光速。
在实际应用中,接收器通常能够接收到多个卫星的信号,因此可以利用多个卫星的伪距信息进行定位。
伪距单点定位的核心是通过多个卫星的伪距信息求解接收器的位置坐标。
这个问题可以表示为一个数学模型,通过最小二乘法求解,得到接收器的位置坐标。
在Matlab中实现伪距单点定位算法需要以下几个步骤:1. 数据预处理:首先需要将接收器接收到的卫星信号数据进行预处理,包括数据解码、信号强度计算等。
2. 卫星位置计算:利用卫星星历数据,计算卫星在给定时刻的位置。
3. 伪距计算:通过测量卫星信号的传播时间差,计算接收器与卫星之间的伪距。
4. 伪距单点定位:利用多个卫星的伪距信息,通过最小二乘法求解接收器的位置坐标。
5. 定位结果分析:对定位结果进行分析和评估,包括精度评估、误差分析等。
在实际应用中,伪距单点定位算法还需要考虑多种误差的影响,包括钟差误差、大气延迟误差、多径效应等。
这些误差会对定位结果产生影响,需要进行误差补偿和校正。
Matlab伪距单点定位是一种利用卫星信号进行定位的方法,通过测量卫星信号的传播时间差来计算接收器与卫星之间的距离,并利用多个卫星的距离信息进行定位。
Matlab作为强大的数学计算工具,可以方便地实现伪距单点定位算法。
伪距单点定位的实现主要包括数据预处理、卫星位置计算、伪距计算、伪距单点定位和定位结果分析等步骤。
基于MATLAB的GPS水准拟合方法及应用
摘要基于MATLAB的GPS水准拟合方法及应用摘要GPS高程拟合是GPS研究领域的一个热点。
提高GPS的定位精度以及GPS 高程转换精度,可以得到较高精度的GPS水准高程。
实践证明根据实际情况选择正确的转换方法获得的GPS水准高程的精度可以广泛地应用到工程、变形监测等各个方面中去,这将拓宽GPS的应用领域,真正的实现GPS的三维优越性。
使备受青睐的GPS有更好的应用前景。
本文系统的研究了GPS高程拟合原理,分析了各种拟合模型,论述了影响GPS高程的误差源和改正方法,以及精度评定的标准。
论文首先介绍了高程系统起算边及相互关系,其中分别介绍了大地水准面、正高与正高系统,似大地准面、正常高与正常高系统,参考椭球面、大地高与大地高系统,并用作图形象的表示出正高、正常高与大地高之间的关系。
其次分析和探讨了GPS高程拟合的方法,将其分为三类:多项式曲线拟合,多项式曲面拟合和多面函数拟合,然后分别介绍了这三种方法的基本原理和计算步骤。
最后通过实例分析了这三种拟合方法的优缺点。
关键词:高程系统;GPS高程拟合;多面函数拟合;精度分析IABSTRACTGPS Leveling Fitting Method and ApplicationBased on the MATLABABSTRACTGPS elevation fitting is a hot research field about GPS.Improving GPS positioning accuracy of GPS elevation conversion and precision, and can geting a high accuracy of GPS level elevation.Practice has proved according to the actual situation of choosing the right method for the conversion of the accuracy of GPS level elevation can be widely applied to engineering, deformation monitoring and other aspects, this will broaden the application field of GPS, really realize the GPS 3 d superiority.Make very exciting and popular better application prospect of GPS. This paper studies the elevation of the system GPS fitting principle and the analysis of the various fitting model, discusses the influence of the error sources and GPS elevation correction method, and evaluate accuracy standard.It firstly introduced vertical system date edge and mutual relationship, which respectively introduced the geoid, ortho metric height and ortho metric height system,Like the earth must face, normal high and normal high system,reference ellipsoid, like the earth must face, normal high and normal high system reference ellipsoid, the earth and the earth high high system, and the image of the drawing showed ortho metric height normal high and the earth, the relationship between the high.Secondly analyzed and discussed on GPS elevation fitting method, which is divided into three categories: polynomial curve fitting, polynomial fitting of surface and many-sided function fitting, then respectively introduces the three basic principle and method of the calculation procedure.Finally the example analyzed the three kind of fitting and the advantages and disadvantages of the methods.Keywords: Vertical system; GPS elevation fitting; The function fitting; Precision analysisII目录摘要 (I)ABSTRACT (II)目录...................................................................................................................... I II 第一章绪论. (1)1.1背景及意义 (1)1.1.1本课题的背景 (1)1.1.2本课题的意义 (1)1.2国内外研究现状、水平 (2)1.3本课题研究的主要内容 (3)第二章GPS高程的基本理论 (5)2.1高程系统 (5)2.11.大地水准面、正高与正高系统 (5)2.1.2.似大地准面、正常高与正常高系统 (5)2.1.3参考椭球面、大地高与大地高系统 (6)2.1.4 正高、正常高与大地高之间的关系 (6)2.3本章小结 (7)第三章GPS高程拟合的方法研究 (8)3.1 GPS水准高程拟合的基本原理 (8)3.2多项式曲线拟合法 (8)3.3多项式曲面拟合法 (10)3.4多面函数拟合法 (11)3.4本章小结 (13)第四章案例分析 (15)4.1精度评定 (15)4.1.1适用范围 (15)4.1.2选出合适的高程异常已知点 (15)4.1.3 高程异常已知点的数量 (15)4.1.4 GPS拟合精度分析 (15)III4.2 MATLAB软件特点及应用 (16)4.3狭长线性区域拟合实例 (17)4.4丘陵地区高程异常拟合实例 (20)4.5本章小结 (26)第五章总结与展望 (27)5.1总结 (27)5.2展望 (27)致谢 (29)参考文献 (30)附录 (32)多项式曲线拟合MATLAB程序 (32)多项式曲面拟合MATLAB程序 (32)多面函数拟合MATLAB程序 (33)IV南京工业大学本科生毕业设计(论文)第一章绪论1.1背景及意义1.1.1本课题的背景传统的高程测量采用几何水准测量的方式,但是在地形复杂、高差较大的地区进行水准测量非常困难,大面积水准测量往往要耗费巨大的人力、物力,且效率极低,在大量的手工记录过程中难免出现漏记、错记等情况。
matlab编程定位方法
matlab编程定位方法
在MATLAB中,可以使用多种方法进行定位,这取决于你想要解决的具体
问题。
下面是一些基本的定位方法:
1. 二维图像定位:对于在二维图像中的目标定位,可以使用边缘检测、特征提取、模板匹配等技术。
例如,使用 `edge` 函数进行边缘检测,使用
`regionprops` 函数提取区域属性,使用 `immatch` 函数进行模板匹配等。
2. 三维空间定位:对于三维空间中的目标定位,可能需要结合多个传感器(如GPS、IMU、轮速传感器等)的数据。
这通常涉及到数据融合和卡尔曼滤波等技术。
3. 声音定位:对于声音定位,可以使用声波传播模型和阵列处理技术。
例如,使用 `audioread` 函数读取声音信号,使用 `beamform` 函数进行波束形成等。
4. 机器人定位:在机器人定位中,可以使用 SLAM(Simultaneous Localization and Mapping)技术。
MATLAB 提供了一些工具箱,如Robotics System Toolbox,可以帮助实现这一目标。
5. 雷达/激光雷达定位:对于雷达或激光雷达数据,可以使用点云处理技术。
MATLAB 的Point Cloud Processing Toolbox 可以提供一些基本的工具。
以上只是一些基本的定位方法,具体的实现会根据你的具体需求和数据类型有所不同。
如果你能提供更具体的信息(如数据类型、目标特性、应用场景等),我可以给出更详细的建议。
利用MATLAB编制的GPS单点定位程序
function t=TimetoJD(Y,M,D,h,f,s)if(Y>=80)Y=Y+1900;elseY=Y+2000;endif M<=2Y=Y-1;M=M+12;endJD=fix(365.25*Y)+fix(30.6001*(M+1))+D+h/24+f/1440+s/24/3600+1720981.5;t=JD-2444244.5;function [head,obs]=ReadObsData%读接收机观测数据文件%HeadODat :a structure stores header information if o-file% .ApproXYZ[3]; //approximate coordinate% .interval; //intervals of two adjacent epochs% .SiteName; //point name% .Ant_H; //antenna height% .Ant_E; //east offset of the antenna center% .Ant_N; //north offset of then antenna center% .TimeOB; //first epoch time in modifuied Julian time% .TimeOE; //last epoch time in modifuied Julian time% .SumOType; //number of types of observables% .SumOO[SumOType]; //type of observables 0-L1,1-L2,2-C1,3-P1,4-P2,5-D1,6-D2,%ObsODat :a structure stores observables by one and one epoch% .TimeOEpp[2]; //reciever time of epoch% .SatSum; //number of satellites% .SatCode[SatSum]; //satellites' PRN% .Obs_FreL1[SatSum];% .Obs_FreL2[SatSum];% .Obs_RangeC1[SatSum];% .Obs_RangeP1[SatSum];% .Obs_RangeP2[SatSum];global HeadODat;global ObsODat;[fname,fpath]=uigetfile('*.*','选择一个O文件');O_filename=strcat(fpath,fname);fid1=fopen(O_filename,'rt');if (fid1==-1)msgbox('file invalide','warning','warn');return;end%将文件头数据存入结构体HeadODat中t=0;while(t<100)s=fgets(fid1);t=t+1;L=size(s,2);if L<81s(L+1:81)=' ';endinstrS=s(61:81);%测站点近似坐标if strncmp(instrS,'APPROX POSITION XYZ',19)HeadODat.ApproXYZ=zeros(1,3);HeadODat.ApproXYZ(1,1)=str2num(s(1:14));HeadODat.ApproXYZ(1,2)=str2num(s(15:28));HeadODat.ApproXYZ(1,3)=str2num(s(29:42));%历元间隔;elseif strncmp(instrS,'INTERVAL',8)HeadODat.interval=str2num(s(5:11));%测站点号elseif strncmp(instrS,'MARKER NAME',11)HeadODat.SiteName=s(1:4)%天线中心改正elseif strncmp(instrS,'ANTENNA: DELTA H/E/N',20)HeadODat.Ant_H=str2num(s(1:14));HeadODat.Ant_E=str2num(s(15:28));HeadODat.Ant_N=str2num(s(29:42));%第一个历元时间elseif strncmp(instrS,'TIME OF FIRST OBS',17)year=str2num(s(1:6));month=str2num(s(7:12));day=str2num(s(13:18));hour=str2num(s(19:24));minute=str2num(s(25:30));second=str2num(s(31:42));HeadODat.TimeOB=TimetoJD(year,month,day,hour,minute,second);%最后一个历元时间elseif strncmp(instrS,'TIME OF LAST OBS',16)year=str2num(s(1:6));month=str2num(s(7:12));day=str2num(s(13:18));hour=str2num(s(19:24));minute=str2num(s(25:30));second=str2num(s(31:42));HeadODat.TimeOE=TimetoJD(year,month,day,hour,minute,second);%观测值类型elseif strncmp(instrS,'# / TYPES OF OBSERV',19)HeadODat.SumOType=str2num(s(1:6));HeadODat.SumOO=ones(1,HeadODat.SumOType)*-1;for k=1:HeadODat.SumOTypef=s(k*6+5:k*6+6);if strcmp(f,'L1')HeadODat.SumOO(1,k)=0;elseif strcmp(f,'L2')HeadODat.SumOO(1,k)=1;elseif strcmp(f,'C1')HeadODat.SumOO(1,k)=2;elseif strcmp(f,'P1')HeadODat.SumOO(1,k)=3;elseif strcmp(f,'P2')HeadODat.SumOO(1,k)=4;elseif strcmp(f,'D1')HeadODat.SumOO(1,k)=5;elseif strcmp(f,'D2')HeadODat.SumOO(1,k)=6;endend%头文件结束elseif strncmp(instrS,'END OF HEADER',13)break;elsecontinue;endend%观测数据结构体%观测数据结构t=0;while ~feof(fid1)%每个历元的第一行数据,时间和观测到的卫星号s=fgets(fid1);t=t+1;year=str2num(s(1:3));month=str2num(s(4:6));day=str2num(s(7:9));hour=str2num(s(10:12));minute=str2num(s(13:15));second=str2num(s(16:26));%历元时间ObsODat(t).TimeOEp=[year,month,day,hour,minute,second];ObsODat(t).TimeOEpp=TimetoJD(year,month,day,hour,minute,second);%该历元观测到的卫星数ObsODat(t).SatSum=str2num(s(30:32));%该历元观测到的卫星号ObsODat(t).SatCode=zeros(1,ObsODat(t).SatSum);ObsODat(t).Obs_FreL1=zeros(1,ObsODat(t).SatSum);ObsODat(t).Obs_FreL2=zeros(1,ObsODat(t).SatSum);ObsODat(t).Obs_RangeC1=zeros(1,ObsODat(t).SatSum);ObsODat(t).Obs_RangeP1=zeros(1,ObsODat(t).SatSum);ObsODat(t).Obs_RangeP2=zeros(1,ObsODat(t).SatSum);for k=1:ObsODat(t).SatSumf=s(31+k*3:32+k*3);ObsODat(t).SatCode(1,k)=str2num(f);end%每个历元的观测数据,按卫星号先后顺序分行存for k=1:ObsODat(t).SatSums=fgets(fid1);%判断一个卫星的观测数据是否分两行存储,如果为两行,则再读入一行if HeadODat.SumOType>5sg=fgets(fid1);s=strcat(s,sg);endL=size(s,2);%补充数据长度if L<HeadODat.SumOType*16s(L+1:HeadODat.SumOType*16)=' ';end%对观测数据判断其类型,并存储到相应的数组中for j=1:HeadODat.SumOTypestemp=s(j*16-15:j*16);stemp=deblank(stemp);if isempty(stemp)continue;endstempNum=str2num(stemp);stempLength=size(stempNum,2);if stempLength>1stempNum=stempNum(1,1);endif HeadODat.SumOO(1,j)==0ObsODat(t).Obs_FreL1(1,k)=stempNum;elseif HeadODat.SumOO(1,j)==1ObsODat(t).Obs_FreL2(1,k)=stempNum;elseif HeadODat.SumOO(1,j)==2ObsODat(t).Obs_RangeC1(1,k)=stempNum;elseif HeadODat.SumOO(1,j)==3ObsODat(t).Obs_RangeP1(1,k)=stempNum;elseif HeadODat.SumOO(1,j)==4ObsODat(t).Obs_RangeP2(1,k)=stempNum;elsecontinue;endend%完成一个卫星的所有观测数据存储end%完成一个历元的数据存储end%完成所有历元的数据存储head=HeadODat;obs=ObsODat;returnfunction EphDat=ReadGpsDataglobal EphDatEphDat=struct;[filename,pathname]=uigetfile('*.**N','读取GPS广播星历文件'); fid1=fopen(strcat(pathname,filename),'rt');if(fid1==-1)msgbox('Input File or Path is not correct','warning ','warn');return;endwhile(1)temp=fgetl(fid1);header=findstr(temp,'END OF HEADER');if(~isempty(header))break;endendi=1;while(1)temp=fgetl(fid1);break;endEphDat(i).PRN=str2double(temp(1:2));year=str2double(temp(3:5));EphDat(i).toc=TimetoJD(year,str2double(temp(6:8)),str2double(temp(9:11)),str2double(temp(12: 14)),str2double(temp(15:17)),str2double(temp(18:22)));EphDat(i).a0=str2num(temp(23:41));EphDat(i).a1=str2num(temp(42:60));EphDat(i).a2=str2num(temp(61:79));temp=fgetl(fid1);EphDat(i).idoe=str2num(temp(4:22));EphDat(i).Crs=str2num(temp(23:41));EphDat(i).dn=str2num(temp(42:60));EphDat(i).M0=str2num(temp(61:79));temp=fgetl(fid1);EphDat(i).Cuc=str2num(temp(4:22));EphDat(i).e=str2num(temp(23:41));EphDat(i).Cus=str2num(temp(42:60));EphDat(i).sqa=str2num(temp(61:79));temp=fgetl(fid1);EphDat(i).toe=str2num(temp(4:22));EphDat(i).Cic=str2num(temp(23:41));EphDat(i).omg0=str2num(temp(42:60));EphDat(i).Cis=str2num(temp(61:79));temp=fgetl(fid1);EphDat(i).i0=str2num(temp(4:22));EphDat(i).Crc=str2num(temp(23:41));EphDat(i).w=str2num(temp(42:60));EphDat(i).omg=str2num(temp(61:79));temp=fgetl(fid1);EphDat(i).iDot=str2num(temp(4:22));EphDat(i).cflgl2=str2num(temp(23:41));EphDat(i).weekno=str2num(temp(42:60));EphDat(i).pflgl2=str2num(temp(61:79));temp=fgetl(fid1);EphDat(i).svacc=str2num(temp(4:22));EphDat(i).svhlth=str2num(temp(23:41));EphDat(i).tgd=str2num(temp(42:60));EphDat(i).aodc=str2num(temp(61:79));temp=fgetl(fid1);EphDat(i).ttm=str2num(temp(4:22));if(i~=1) %删除重复数据if (EphDat(i-1).PRN~=EphDat(i).PRN)break;elseif abs(EphDat(i-1).toc-EphDat(i).toc)<0.001i=i - 1;endendendi=i + 1;endfunction DDDWformat longclear allticglobal HeadODat;global ObsODat;global EphDat;%先读N文件,再读O文件EphDat=ReadGpsData;[HeadODat,ObsODat]=ReadObsData;%多个历元加权平均求测站点坐标N = size(EphDat,2);c=299792458;epochNUM=size(ObsODat,2);%观测数据的个数Xk0=HeadODat.ApproXYZ(1,1); %测站点的近似坐标Yk0=HeadODat.ApproXYZ(1,2);Zk0=HeadODat.ApproXYZ(1,3);for k=1:epochNUMtr=ObsODat(k).TimeOEp;%历元接收数据时间Snum=ObsODat(1,k).SatSum; %该历元观测到的卫星数if Snum<4break; %去除无解的历元endCode=ObsODat(k).SatCode;%该历元观测到的卫星号组R=ObsODat(k).Obs_RangeC1 ; %取出C1观测值,列向量%卫星发射时间的迭代计算for j=1:Snumif R(j) == 0continue;endt=TimetoJD(tr(1),tr(2),tr(3),tr(4),tr(5),tr(6));t2 = mod(t,7)*24*3600;%gps second%卫星钟差dd=-1;for i=1:Nif(EphDat(i).PRN==Code(j)&abs(t-EphDat(i).toc)<0.0417*2)%读入时间相近的星历文件a0= EphDat(i).a0;a1=EphDat(i).a1;a2= EphDat(i).a2;toe=EphDat(i).toe;dt=a0+(t2-toe)*a1+(t2-toe)^2*a2;%卫星钟差dd=1;break;endendif dd==-1msgbox('没有相关星历文件');return;endtt = tr;ts = tr(6) - R(j)/c;%用秒进行迭代sign = 1;while(sign>1E-8)tt(6) = ts;[Xs1,Ys1,Zs1]= CalPos(Code(j),tt);ss1 = sqrt((Xk0-Xs1)^2+(Yk0-Ys1)^2+(Zk0-Zs1)*(Zk0-Zs1));%卫星位置加地球自转改正qe=0.000072921151467; %地球自转角速度delt=ss1/c;Rotation=[cos(qe*delt),sin(qe*delt),0;-sin(qe*delt),cos(qe*delt),0;0,0,1];h=Rotation*[Xs1;Ys1;Zs1] ;Xs=h(1);Ys=h(2);Zs=h(3);ss = sqrt((Xk0-Xs)^2+(Yk0-Ys)^2+(Zk0-Zs)*(Zk0-Zs));ts = tr(6)-ss/c;sign = norm(ts-tt(6));endaxk = (Xk0-Xs)/ss;ayk = (Yk0-Ys)/ss;azk = (Zk0-Zs)/ss;EAB=CAL2POL(Xk0,Yk0,Zk0,Xs,Ys,Zs);if j==1A = [axk,ayk,azk,1];L = R(j)-ss++c*dt-2.47/(sin(EAB)+0.0121);elseA = [A;axk,ayk,azk,1]; %构造误差方程L = [L;(R(j)-ss++c*dt)-2.47/(sin(EAB)+0.0121)];%常数项endendif Snum==4dx=inv(A)*L;aaaa(k).Q=inv(A'*A);Px(k) = 1/aaaa(k).Q(1,1);Py(k) = 1/aaaa(k).Q(2,2);Pz(k) = 1/aaaa(k).Q(3,3);xchange(k) = dx(1);ychange(k) = dx(2);zchange(k) = dx(3);elseif Snum > 4dx = inv(A'*A)*A'*L; %构造法方程并求解aaaa(k).Q=inv(A'*A);Px(k) = 1/aaaa(k).Q(1,1);Py(k) = 1/aaaa(k).Q(2,2);Pz(k) = 1/aaaa(k).Q(3,3);xchange(k) = dx(1);ychange(k) = dx(2);zchange(k) = dx(3);endendXp=sum(Px.*(Xk0+xchange))/sum(Px)Yp=sum(Py.*(Yk0+ychange))/sum(Py)Zp=sum(Pz.*(Zk0+zchange))/sum(Pz)figure(1);subplot(3,1,1);plot(xchange,'black--');subplot(3,1,2);plot(ychange,'black--');subplot(3,1,3);plot(zchange,'black--');tocfunction [x,y,z]= CalPos(PRN,time)global EphDatt1=TimetoJD(time(1),time(2),time(3),time(4),time(5),time(6));t2=TimetoJD(time(1),time(2),time(3),0,0,0);if isempty(EphDat)EphDat=ReadGpsData;endsz=size(EphDat,2);gg=0;%判断最近的时间for i=1:szif(EphDat(i).PRN==PRN & abs(EphDat(i).toc-t1)<0.0417*2)G=EphDat(i);gg=1;break;endendt3=t2-G.weekno*7;% 减整周数t=t3*24*3600+time(4)*3600+time(5)*60+time(6);%化为GPS秒u=3.986004418E+14; %地球引力常数we=7.2921151467E-5; %地球自转速度a=G.sqa^2; %地球长半轴n0=sqrt(u/a^3); %计算的平运动tk=t-G.toe; %从参考历元开始的时间if tk>=302400tk=tk-604800;elseif tk<-302400tk=tk + 604800;endn=n0+G.dn; %改正后的运动Mk=G.M0+n*tk;Ek=Mk; %偏近点角弧度for i=1:4Ek=Mk+G.e*sin(Ek);endfk=2*atan((sqrt((1+G.e)/(1-G.e))*tan(Ek/2))); %真近点角if(fk<0)fk=fk+2*pi;endcosf=(cos(Ek)-G.e)/(1-G.e*cos(Ek));sinf=(sqrt(1-G.e^2)*sin(Ek))/(1-G.e*cos(Ek));uk=fk+G.w; %升交角矩及轨道摄动改正项iuk=G.Cuc*cos(2*uk)+G.Cus*sin(2*uk);irk=G.Crc*cos(2*uk)+G.Crs*sin(2*uk);iik=G.Cic*cos(2*uk)+G.Cis*sin(2*uk);uk=uk+iuk ; %改正后的升角距r=a*(1-G.e*cos(Ek))+irk; %改正后的地心相径i=G.i0+iik+G.iDot*tk; %改正后的倾角xx=r*cos(uk); %在轨道面中的位置yy=r*sin(uk);omg=G.omg0+(G.omg-we)*tk-we*G.toe; %改正后的升交点经度x=xx*cos(omg)-yy*cos(i)*sin(omg);y=xx*sin(omg)+yy*cos(i)*cos(omg);z=yy*sin(i);CalPos=[x,y,z];%打印结果% x=num2str(x,'%12.6f');% y=num2str(y,'%12.6f');% z=num2str(z,'%12.6f');% fprintf('卫星号PRN:%2.0f',PRN); fprintf('\n');% fprintf('时间time:%4.0f%4.0f%4.0f%4.0f%4.0f%4.1f',time); fprintf('\n');% fprintf('所求卫星在地心坐标系中的空间坐标为:\n');% fprintf('X = ');fprintf(x,'\n');fprintf('m\n');% fprintf('Y = ');fprintf(y,'\n');fprintf('m\n');% fprintf('Z = ');fprintf(z,'\n');fprintf('m\n');%returnfunction EAB=CAL2POL(X,Y,Z,XB,YB,ZB)format longL=atan(Y/X);R=sqrt(X*X+Y*Y+Z*Z);a=6378137;e=sqrt(0.0066943799013);seta=atan(Z/sqrt(X*X+Y*Y));B=asin(sqrt((a*a-R*R*cos(seta)*cos(seta))/(a*a-R*R*cos(seta)*cos(seta)*e*e)))B1=atan(tan(seta)*(1+a*e*e*sin(B)/Z/sqrt(1-e*e*sin(B)*sin(B))));while abs(B1-B)>0.01*pi/180/3600B=B1;B1=atan(tan(seta)*(1+a*e*e*sin(B)/Z/sqrt(1-e*e*sin(B)^2)));endH=R*cos(seta)/cos(B)-(a/sqrt(1-e*e*sin(B)*sin(B)));L=L*180/pi;B=B*180/pi;DX=XB-X;DY=YB-Y;DZ=ZB-Z;NB=-sin(B)*cos(L)*DX-sin(B)*sin(L)*DY+cos(B)*DZ; EB=-sin(L)*DX+cos(L)*DY;UB=cos(B)*cos(L)*DX+cos(B)*sin(L)*DY+sin(B)*DZ; SAB=sqrt(NB*NB+UB*UB+EB*EB);EAB=asin(UB/SAB);。
运用matlab计算gps卫星的坐标位置
运用matlab计算gps卫星的坐标位置clearclcformat longtp=input('tp=');toc=input('toc=');a0=input('a0=');a1=input('a1=');a2=input('a2=');toe=input('toe=');M0=input('M0=');a=input('长半径a=');deltan=input('卫星平均角速度之差deltan=');e=input('e=');w=input('w=');Cuc=input('Cuc=');Cus=input('Cus=');Cic=input('Cic=');Cis=input('Cis=');Crc=input('Crc=');Crs=input('Crs=');i0=input('i0=');I=input('轨道倾角变化率I=');OM0=input('OM0=');OM=input('升交点赤径变化率OM=');tt=a0+a1*(tp-toc)+a2*(tp-toc);t=tp-tt;tk=t-toe;u=3.986005e14;n0=(sqrt(u))/(a*a*a);n=n0+deltan;Mk=M0+n*tk;{n=MK;ek0=0;ek1=n+e*sin(ek0);ek2=n+e*sin(ek1);EK=ek2;}Dk=1;Ek=0;n1=0;while abs(Ek-Dk)>0.0000000001n1=n1+1;Ek=Dk;Dk=Mk+e*sin(Ek);endEk=Dk;Vk=atan((sqrt(1-e*e)*sin(Ek)/(cos(Ek)-e)); if sin(Ek)>0 & cos(Ek)-e<0Vk=pi-Vk;elseif sin(Ek)<0 & cos(Ek)-e<0Vk=pi+Vk;elseif sin(Ek)<0 & cos(Ek)-e>0Vk=2*pi-Vk;endFaik=Vk+w;SigmaU=Cuc*cos(2*Faik)+Cus*sin(2*Faik); SigmaR=Crc*cos(2*Faik)+Crs*sin(2*Faik);SigmaI=Cic*cos(2*Faik)+Cis*sin(2*Faik);Uk=Faik+SigmaU;Rk=a*(1-e*cos(Ek))+SigmaR;Ik=i0+SigmaI+I*tk;X0=Rk*cos(Uk);Y0=Rk*sin(Uk);we=7.29211567e-5;OMK=OM0+(OM-we)*tk-we*toe;Xk=X0*cos(OMK)-Y0*cos(Ik)*sin(OMK);Yk=X0*sin(OMK)+Y0*cos(Ik)*cos(OMK);Zk=Y0*sin(Ik);disp(['卫星钟差改正dt=',num2str(tt)])disp(['归化时刻tk=',num2str(tk)])disp(['平均运行角速度n=',num2str(n)])disp(['卫星平近点角Mk=',num2str(Mk)])disp(['偏近点角Ek=',num2str(Ek)])disp(['真近点角Vk=',num2str(Vk)])disp(['升交距角Faik=',num2str(Faik)])disp(['摄动改正项SigmaU=',num2str(SigmaU)]) disp(['摄动改正项SigmaR=',num2str(SigmaR)]) disp(['摄动改正项SigmaI=',num2str(SigmaI)]) disp('经过摄动改正项:')disp(['升交距角Uk=',num2str(Uk)])disp(['卫星矢径Rk=',num2str(Rk)])disp(['轨道倾角Ik=',num2str(Ik)])disp('卫星在轨道平面坐标系的坐标')disp(['X0=',num2str(X0)])disp(['Y0=',num2str(Y0)])disp(['观测时刻升交点经度OMK=',num2str(OMK)]) disp('卫星在地心固定坐标系中的直角坐标')disp(['Xk=',num2str(Xk)]) disp(['Yk=',num2str(Yk)]) disp(['Zk=',num2str(Zk)])。
如何使用Matlab进行卫星导航与定位
如何使用Matlab进行卫星导航与定位卫星导航与定位是一门广泛应用于航空、航海、地理信息系统等领域的技术。
在这个时代,人们越来越依赖全球定位系统(GPS)来获取准确的位置信息。
而Matlab作为一种功能强大的数学软件工具,可以帮助我们进行卫星导航与定位的模拟和算法设计。
本文将介绍如何使用Matlab进行卫星导航与定位的相关内容。
在开始之前,我们需要明确一些基础知识。
首先,我们需要了解GPS工作的原理。
GPS系统主要由卫星、地面控制站和接收器三部分组成。
卫星发送信号,接收器接收并计算信号传播的时间差来确定位置。
其次,我们需要了解GPS的信号类型。
GPS信号包括C/A码和P码两种类型,其中C/A码用于民用接收器,P码用于军用接收器。
最后,我们需要了解GPS的误差来源。
GPS的定位误差主要来自多普勒效应、大气延迟、钟差等因素。
使用Matlab进行卫星导航与定位需要用到一些特定的工具和函数。
首先,我们需要导入相应的工具箱。
Matlab提供了GPS工具箱和导航工具箱,可以帮助我们进行相关的计算和模拟。
在导入工具箱之后,我们可以使用相关函数进行卫星导航与定位的模拟和计算。
例如,可以使用gpscoordinates函数来计算卫星的位置信息,使用gps2utc函数来进行时间转换等。
在进行卫星导航与定位的模拟和计算之前,我们需要准备相关的数据。
首先,我们需要获取卫星的星历数据。
星历数据包括卫星的位置、速度和加速度等信息,可以帮助我们计算卫星的轨道和位置。
其次,我们需要获取接收器的观测数据。
观测数据包括接收器接收到的信号的时间和信号的强度等信息,可以帮助我们计算信号传播的时间差和定位误差。
最后,我们需要获取地球的几何形状数据。
地球的几何形状数据包括地球的椭球体参数和地球的地理坐标系等信息,可以帮助我们进行地球坐标和地理位置的转换。
在得到所需的数据之后,我们就可以开始使用Matlab进行卫星导航与定位的模拟和计算了。
多系统伪距单点定位matlab代码
多系统伪距单点定位是一种通过接收来自多个卫星的信号并利用其伪距信息进行单点定位的方法。
这种定位方法适用于GPS、GLONASS、Beidou和Galileo等多个全球卫星定位系统,可以提高定位的准确性和鲁棒性。
在实际应用中,通过Matlab代码实现多系统伪距单点定位是一种常见的方法。
本文将介绍多系统伪距单点定位的原理,以及利用Matlab代码实现这一定位方法的步骤和技巧。
一、多系统伪距单点定位原理多系统伪距单点定位主要依赖于接收来自多个卫星的信号,并通过测量信号的伪距值来进行定位。
接收机接收到来自不同卫星的信号,并测量其与接收机的距离,然后根据测得的伪距值和卫星的位置信息,利用三角定位原理计算出接收机的位置。
在实际应用中,由于不同卫星系统的误差和干扰,需要对各个卫星信号的伪距值进行合理的处理,如加权平均、差分定位等,以提高定位的精度和鲁棒性。
二、Matlab代码实现多系统伪距单点定位步骤1. 读取卫星星历数据和接收机观测数据在Matlab中,可以利用现有的卫星星历数据和接收机观测数据进行多系统伪距单点定位。
通过Matlab的文件读取函数读取卫星星历数据和接收机的观测数据,然后将其存储为Matlab中的数据结构,并进行初步的数据处理。
2. 卫星信号伪距解算利用读取的卫星星历数据和接收机观测数据,可以利用现有的卫星信号伪距解算算法,在Matlab中编写代码对不同卫星信号的伪距值进行计算和处理。
这一步骤需要考虑卫星信号的误差和干扰因素,对伪距值进行合理的加权、滤波和校正,以提高定位的精度和鲁棒性。
3. 接收机定位计算在完成卫星信号伪距解算之后,可以利用Matlab中的三角定位算法对接收机的位置进行计算。
根据处理后的卫星信号伪距值和卫星位置信息,编写Matlab代码实现接收机的位置计算,并利用可视化工具展示定位结果,以便对定位精度和鲁棒性进行评估和分析。
4. 定位精度分析与优化可以利用Matlab对定位结果进行精度分析和优化。
matlab伪距单点定位
matlab伪距单点定位Matlab伪距单点定位伪距单点定位是一种常用的定位方法,通过测量接收器与多个卫星之间的信号传播时间差,从而确定接收器的位置。
在Matlab中,我们可以利用伪距单点定位算法来实现定位功能。
我们需要了解一些基本概念。
伪距是指卫星信号从发射到接收的时间差,乘以光速即可得到距离。
伪距单点定位算法的基本原理是通过多个卫星的伪距观测值,利用数学模型来计算接收器的位置坐标。
在Matlab中,我们可以使用GPS工具箱来实现伪距单点定位。
首先,我们需要获取卫星的位置信息以及接收器的伪距观测值。
这些数据可以通过GPS接收器获取或者从GPS导航数据中提取。
接下来,我们可以使用伪距单点定位算法来计算接收器的位置。
算法的核心思想是利用卫星的位置信息和伪距观测值,通过最小二乘法求解接收器的位置坐标。
具体的计算步骤如下:1. 根据卫星位置信息和接收器的伪距观测值,构建一个误差函数。
误差函数描述了接收器位置与观测值之间的差异。
2. 使用最小二乘法求解误差函数的最小值。
最小二乘法是一种常用的优化方法,通过最小化误差函数来确定最优解。
3. 根据最小二乘法的结果,得到接收器的位置坐标。
这个坐标表示了接收器在三维空间中的位置。
在Matlab中,我们可以使用函数lsqnonlin来实现最小二乘法的求解。
这个函数可以通过迭代的方式,不断优化误差函数,直到找到最优解。
除了伪距单点定位算法,Matlab还提供了其他定位算法的实现,例如差分定位、扩展卡尔曼滤波等。
这些算法可以根据具体的应用需求进行选择。
总结起来,Matlab提供了强大的工具箱和函数,可以方便地实现伪距单点定位算法。
通过合理地利用这些工具,我们可以快速准确地实现定位功能。
同时,我们也需要注意数据的质量和算法的选择,以保证定位结果的准确性和稳定性。
伪距单点定位是一种常用的定位方法,Matlab提供了丰富的工具和函数来实现该算法。
通过合理地利用这些工具,我们可以轻松实现定位功能,并且可以根据具体需求选择适合的定位算法。
用MATLAB计算GPS卫星位置
用MATLAB计算GPS卫星位置摘要:本文主要介绍了GPS测量数据的常用格式RINEX标准文件格式,并利用MA TLAB工具计算出所观测卫星里的五颗卫星(14、20、29、31和32五颗)在283个历元的瞬时位置,即所观测时间段里五颗卫星在WGS-84坐标下的空间运行轨迹。
关键词:RINEX标准文件WGS-84下卫星位置MATLAB工具GPS定位的基本原理简单来说就是在WGS-84空间直角坐标系中,确定未知点与GPS卫星的空间几何关系。
因此利用GPS进行导航和测量时,卫星是作为位置已知的高空观测目标。
那么如何精确快速的解算出卫星在空间运行的轨迹即其轨道是实现未知点快速定位的关键。
1 标准格式RINEX格式简述在进行GPS数据处理时,由于接收机出自于不同厂家,所以厂家设计的数据格式也是五花八门的,但是在实际中,很多时候需要把来自不同型号的接收机的数据放在一块进行处理,这就需要数据格式的统一,为了解决这种矛盾,RINEX(英文全称为:The Receiver Independent Exchange Format)格式则应运而生,该格式存储数据的类型是文本文件,数据记录格式是独立于接收机的出自厂家和具体型号的。
由此可见,其特点是:由于是通用格式,所以可将不同型号接收机收集的数据进行统一处理,并且大多数大型数据处理软件都能够识别处理,此外也适用于多种型号的接收机联合作业,通用性很强。
RINEX标准文件里不是单一的一个文件,而是包括如下几种类型的文件[1]。
2 卫星坐标的计算步骤由于在GPS定位和导航的时候,用户都是把GPS卫星的位置作为已知量来对待,并且GPS定位所用的坐标系是世界大地坐标系WGS-84。
所以就先必须根据GPS接收机观测的相应星历数据,解算出GPS卫星在WGS-84坐标系中的瞬时位置。
为了后面计算方便,先对广播星历中涉及到的计算卫星坐标的一些轨道参数进行说明,如表4所示。
由于每隔两个小时,GPS接收机收到的广播星历才更新一次,所以用户在根据接收机收到的卫星导航电文汇总的广播星历参数推算GPS的瞬时坐标的时候,一定要选取与GPS卫星的瞬时坐标时刻最相近的那组广播星历数据[2],否则误差将会很大。
GPS精密单点定位
基于Matlab 编程的GPS 伪距单点定位摘 要: GPS 伪距单点定位速度快、不存在整周模糊度,因此具有很大的应用价值。
给出GPS 伪距单点定位、相关改正的计算模型以及精度评定方法即精度因子的计算方法,通过实例并结合matlab 编程,实时的得到单历元伪距单点定位的坐标,并给出单历元伪距单点定位坐标精度优于10m 的结论。
关键词:GPS ;matlab ;单点定位;程序实现;伪距算法1 前言全球定位系统,简称GPS ,就是利用GPS 定位卫星,在全球范围内实时进行定位、导航的系统具有自动连续、全球覆盖、全天候、保密性强,实时定位速度快等特点,其应用领域正在不断的得到拓展 。
而伪距单点定位因其速度快捷 、灵活方便且无多值性问题等特点,能够很好地 满足实时测量的要求,被广泛的用于车辆、舰船 、飞机 的导航,地质矿产的野外勘测以及海洋捕鱼等领域。
但是在实际定位中, 由于卫星钟差、地球自转、对流层折射以及电离层折射等因素对伪距单点定位的影响,其本身的定位精 度受到一定的限制,难以达到很好的精度。
所以,对伪距单点定位进行研究很有必要。
理论上,只需要以三颗 以上GPS 卫星到地面观测点之间的距离 ,进行空间距离后方交会,就可以求得观测点的三维坐标 。
考虑到接收机钟差对 定位的影响,在平差计算时将接收机差作为参数与位置参数一起进行联合解算 。
本文给出了伪距单点定位 的计算模型 以及一些改正的计算公式,通过实例并结合matlab 编程对定位精度进行了一定的分析 。
2 GPS 伪距单点定位原理在任意历元,根据伪距观观测方程,对所测得伪距观测值可列出以下误差方程:V B x l δ=- (1)其中 V 为伪距观测值的改正值,(,,,)T x x y z c t δδδδδ=为接收机相位中心坐标近似值和接收机钟差的改正数,B 为误差方程系数:T000000其他矩阵s s s X X Y Y Z Z B ρρρ⎛⎫--- ⎪= ⎪ ⎪⎝⎭ (2) 式中的000(X ,,)Y Z 为本历元接收机的近似坐标,(X ,,)s s s Y Z 为卫星坐标,具体计算可参考[1]。
用MATLAB计算GPS卫星位置-最新文档资料
用MATLAB计算GPS卫星位置GPS定位的基本原理简单来说就是在WGS-84空间直角坐标系中,确定未知点与GPS卫星的空间几何关系。
因此利用GPS进行导航和测量时,卫星是作为位置已知的高空观测目标。
那么如何精确快速的解算出卫星在空间运行的轨迹即其轨道是实现未知点快速定位的关键。
1 标准格式RINEX格式简述在进行GPS数据处理时,由于接收机出自于不同厂家,所以厂家设计的数据格式也是五花八门的,但是在实际中,很多时候需要把来自不同型号的接收机的数据放在一块进行处理,这就需要数据格式的统一,为了解决这种矛盾,RINEX(英文全称为:The Receiver Independent Exchange Format)格式则应运而生,该格式存储数据的类型是文本文件,数据记录格式是独立于接收机的出自厂家和具体型号的。
由此可见,其特点是:由于是通用格式,所以可将不同型号接收机收集的数据进行统一处理,并且大多数大型数据处理软件都能够识别处理,此外也适用于多种型号的接收机联合作业,通用性很强。
RINEX标准文件里不是单一的一个文件,而是包括如下几种类型的文件[1]。
(1)观测数据文件(ssssdddf.yyo),记录的是GPS观测值信息,(OBServation data,简写OBS,为接收机记录的伪距、相位观测值;O文件,如XG012191.10O)。
(2)导航电文文件(ssssdddf.yyn),记录的是GPS卫星星历信息(NAVavigation data,简写NAV,记录实时发布的广播星历;N文件,如XG012191.10N)。
(3)气象数据文件(ssssdddf.yym),主要是在测站处所测定的气象数据(METerological data,简写MET,记录气象仪器观测的温、压、湿度状况;M文件,如XG012191.10M)。
(4)GLONASS导航电文文件(ssssdddf.yyg),记录的是地球同步卫星的导航电文。
基于MATLAB的GPS信号的仿真研究
基于MATLAB的GPS信号的仿真研究一、本文概述随着全球定位系统(GPS)技术的广泛应用,其在导航、定位、授时等领域的重要性日益凸显。
为了更好地理解GPS信号的特性,提高GPS接收机的设计水平和性能,对GPS信号进行仿真研究显得尤为重要。
本文旨在探讨基于MATLAB的GPS信号仿真方法,分析GPS信号的特点,以及如何利用MATLAB这一强大的数值计算环境和图形化编程工具,对GPS信号进行高效、精确的仿真。
文章首先介绍了GPS系统的发展历程、基本原理和信号特性,为后续的信号仿真提供了理论基础。
随后,详细阐述了GPS信号仿真的一般流程,包括信号生成、传播模型、噪声添加等关键环节。
在此基础上,重点介绍了如何利用MATLAB编写GPS信号仿真程序,包括信号生成、传播模型建立、噪声模拟等方面的具体实现方法。
文章还通过实际案例,展示了基于MATLAB的GPS信号仿真在接收机设计、性能评估等方面的应用。
通过仿真实验,可以深入了解GPS信号在不同环境下的传播特性,为接收机算法优化和性能提升提供有力支持。
本文的研究不仅有助于加深对GPS信号特性和仿真方法的理解,也为GPS接收机的研究和开发提供了一种有效的技术手段。
通过MATLAB的仿真研究,可以更加直观地揭示GPS信号的本质规律,为实际应用提供有力指导。
二、GPS信号原理及特性全球定位系统(GPS)是一种基于卫星的无线电导航系统,它利用一组在地球轨道上运行的卫星来提供全球范围内的定位和时间服务。
每个GPS卫星都不断地向地球表面发射射频信号,这些信号被地面上的接收器接收并处理,从而确定接收器的三维位置和速度,以及精确的时间信息。
GPS卫星发射的信号是L波段的射频信号,分为两个频段:L1(142 MHz)和L2(160 MHz)。
每个频段都包含两种类型的信号:C/A码(粗捕获码)和P码(精密码)。
C/A码是对公众开放的,用于民用和商业应用,而P码则用于军事和特定的高精度应用。
matlab伪距单点定位使用牛顿迭代的最小二乘法
matlab伪距单点定位使用牛顿迭代的最小二乘法【中文文章】【标题】Matlab伪距单点定位:使用牛顿迭代的最小二乘法【导言】在全球定位系统(GPS)中,伪距单点定位是一种常见的定位方法。
它通过测量接收器与多颗卫星之间的距离,来确定接收器的位置。
牛顿迭代的最小二乘法是一种数学优化方法,可以帮助提高定位精度。
本文将深入探讨matlab伪距单点定位如何使用牛顿迭代的最小二乘法,从而实现高精度的定位结果。
【正文】1. 介绍GPS定位系统是基于卫星信号的定位系统,它的基本原理是通过计算接收器与卫星之间的距离来确定接收器的位置。
而伪距单点定位是一种常用的GPS定位方法,它通过测量接收器与多颗卫星之间的伪距来计算位置。
在实际应用中,我们经常会遇到接收器与卫星定位存在一定误差的情况,因此需要通过数学方法来提高定位的精度。
2. 最小二乘法原理最小二乘法是一种常用的数学优化方法,可以通过最小化观测值与预测值之间的误差平方和来确定最优解。
在伪距单点定位中,我们可以使用最小二乘法来优化接收器的位置估计。
3. Matlab实现在Matlab中,我们可以使用牛顿迭代的最小二乘法来进行伪距单点定位。
我们需要定义目标函数和其梯度函数。
目标函数即观测值与预测值之间的误差平方和,梯度函数即目标函数对接收器位置的梯度。
通过不断迭代更新接收器位置,我们可以逐渐接近最优解。
4. 算法步骤以下是使用牛顿迭代的最小二乘法进行matlab伪距单点定位的步骤:(1) 初始化接收器位置的初值。
(2) 计算目标函数和其梯度函数的值。
(3) 判断目标函数和梯度函数的收敛条件是否满足,若满足则停止迭代;否则,继续下一步。
(4) 更新接收器位置的值。
(5) 返回第(2)步。
5. 个人观点与理解在进行伪距单点定位时,使用牛顿迭代的最小二乘法可以帮助提高定位的精度和稳定性。
通过不断迭代更新接收器位置,我们可以逐渐接近最优解。
然而,牛顿迭代的最小二乘法也存在一些局限性,例如对初值敏感、迭代次数较多等。
基于MATLAB的GPS信号仿真123
配套的完整源程序代码见百度文库请搜索《基于MATLAB的GPS信号仿真完整源代码123》摘要全球定位定位系统(GPS)是新一代的精密卫星导航定位系统,近年来在民用和军用领域发挥着精确制导的作用。
随着科学技术的发展,GPS导航和定位技术已向高精度、高动态的方向发展。
GPS系统的广泛应用,促使各国政府大力提高发展本国导航定位系统,这就要求我们全面透彻地研究GPS定位系统,为我国的定位导航应用作出贡献。
本文主要研究GPS信号的生成,而在GPS信号的生成过程中,伪码的算法很重要。
m序列是伪码生成的基础,本文首先详细阐述了m序列的生成原理,并进行了相应的自相关性和互相关性分析。
在m序列的基础上,再对C/A码的生成原理进行了详细的介绍。
P码的生成相对来说比较复杂,是伪码的一个重点,本文对其进行了详细的分析。
在MATLAB仿真软件的平台下,本文成功的完成了GPS信号的生成,并且形象的展示了伪码的相关特性,很好的对本文的相关内容进行了仿真。
关键词:全球定位系统;伪码;MATLAB仿真;相关性AbstractAs the new generation of the satellite navigation systems, Global Positioning System (GPS) is more and more important in military and civil fields in recent years. Though the development of science and technology, the technology of navigation and orientation has been progressed to the direction of great precision and dynamic.With the wide application of GPS locating, every country develop our best to make contribution to the positioning and navigation industry of our country.This paper mainly studies the generation of GPS signal, and the algorithm of PN code is very important in GPS signal generation. As m sequence is the basis of PN code generation, this paper firstly explains the m sequence’s generative principle in detail and the characteristics of autocorrelation and cross correlation. On the basis of m series, then this paper introduced C/A code generation principle in detail.As P code is one of the focuses of PN code, and the generation of P code is relatively complex, this paper makes detailed analysis about it.In the platform of MATLAB software simulation, this paper successfully completes the GPS signal generation, and demonstrates the correlation characteristics of PN code visually.It has done a good simulation of the content about this paper.Keywords:Global Positioning System; PN code; MATLAB simulation; correlation characte-ristics目录1 绪论 (1)1.1GPS的应用 (1)1.2本文研究的主要内容 (2)1.3本文研究的目的和意义 (3)2 GPS信号理论及MATLAB软件简介 (4)2.1GPS系统简介 (4)2.2GPS信号结构 (5)2.2.1 载波 (7)2.2.2 伪码 (9)2.2.3 导航电文 (10)2.3伪码 (12)2.3.1 m序列 (12)2.3.2 易捕码C/A码 (16)2.3.3 精码P码 (20)2.4扩频与调制 (24)2.4.1 扩频通信 (24)2.4.2 GPS信号的调制 (25)2.5MATLAB软件简介 (26)2.5.1 MATLAB简介 (26)2.5.2 M文件简介 (27)3 GPS信号仿真程序设计 (33)3.1数据码 (33)3.2C/A码 (34)3.2.1 C/A码的生成及扩频调制 (34)3.2.2 C/A码的相关性分析 (35)3.3P码 (36)3.3.1 P码产生及扩频调制 (36)3.3.2 P码的相关性分析 (39)4 结果分析 (40)4.1C/A码 (40)4.1.1 C/A码的产生及扩频调制 (40)4.1.2 C/A码的相关性分析 (40)4.2P码 (41)4.2.1 P码的产生及扩频调制 (41)4.2.2 P码相关性分析 (43)结论 (45)致谢 (46)参考文献 (47)附录A 英文原文 (48)附录B 汉语翻译 (57)附录C 部分仿真程序代码 (64)1 绪论1.1 GPS的应用全球定位系统(Global Positioning System -GPS)是美国从本世纪70年代开始研制,历时20年,耗资200亿美元,于1994年全面建成,具有在海、陆、空进行全方位实时三维导航与定位能力的新一代卫星导航与定位系统。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
function t=TimetoJD(Y,M,D,h,f,s)if(Y>=80)Y=Y+1900;elseY=Y+2000;endif M<=2Y=Y-1;M=M+12;endJD=fix(365.25*Y)+fix(30.6001*(M+1))+D+h/24+f/1440+s/24/3600+1720981.5;t=JD-2444244.5;function [head,obs]=ReadObsData%读接收机观测数据文件%HeadODat :a structure stores header information if o-file% .ApproXYZ[3]; //approximate coordinate% .interval; //intervals of two adjacent epochs% .SiteName; //point name% .Ant_H; //antenna height% .Ant_E; //east offset of the antenna center% .Ant_N; //north offset of then antenna center% .TimeOB; //first epoch time in modifuied Julian time% .TimeOE; //last epoch time in modifuied Julian time% .SumOType; //number of types of observables% .SumOO[SumOType]; //type of observables 0-L1,1-L2,2-C1,3-P1,4-P2,5-D1,6-D2,%ObsODat :a structure stores observables by one and one epoch% .TimeOEpp[2]; //reciever time of epoch% .SatSum; //number of satellites% .SatCode[SatSum]; //satellites' PRN% .Obs_FreL1[SatSum];% .Obs_FreL2[SatSum];% .Obs_RangeC1[SatSum];% .Obs_RangeP1[SatSum];% .Obs_RangeP2[SatSum];global HeadODat;global ObsODat;[fname,fpath]=uigetfile('*.*','选择一个O文件');O_filename=strcat(fpath,fname);fid1=fopen(O_filename,'rt');if (fid1==-1)msgbox('file invalide','warning','warn');return;end%将文件头数据存入结构体HeadODat中t=0;while(t<100)s=fgets(fid1);t=t+1;L=size(s,2);if L<81s(L+1:81)=' ';endinstrS=s(61:81);%测站点近似坐标if strncmp(instrS,'APPROX POSITION XYZ',19)HeadODat.ApproXYZ=zeros(1,3);HeadODat.ApproXYZ(1,1)=str2num(s(1:14));HeadODat.ApproXYZ(1,2)=str2num(s(15:28));HeadODat.ApproXYZ(1,3)=str2num(s(29:42));%历元间隔;elseif strncmp(instrS,'INTERVAL',8)HeadODat.interval=str2num(s(5:11));%测站点号elseif strncmp(instrS,'MARKER NAME',11)HeadODat.SiteName=s(1:4)%天线中心改正elseif strncmp(instrS,'ANTENNA: DELTA H/E/N',20)HeadODat.Ant_H=str2num(s(1:14));HeadODat.Ant_E=str2num(s(15:28));HeadODat.Ant_N=str2num(s(29:42));%第一个历元时间elseif strncmp(instrS,'TIME OF FIRST OBS',17)year=str2num(s(1:6));month=str2num(s(7:12));day=str2num(s(13:18));hour=str2num(s(19:24));minute=str2num(s(25:30));second=str2num(s(31:42));HeadODat.TimeOB=TimetoJD(year,month,day,hour,minute,second);%最后一个历元时间elseif strncmp(instrS,'TIME OF LAST OBS',16)year=str2num(s(1:6));month=str2num(s(7:12));day=str2num(s(13:18));hour=str2num(s(19:24));minute=str2num(s(25:30));second=str2num(s(31:42));HeadODat.TimeOE=TimetoJD(year,month,day,hour,minute,second);%观测值类型elseif strncmp(instrS,'# / TYPES OF OBSERV',19)HeadODat.SumOType=str2num(s(1:6));HeadODat.SumOO=ones(1,HeadODat.SumOType)*-1;for k=1:HeadODat.SumOTypef=s(k*6+5:k*6+6);if strcmp(f,'L1')HeadODat.SumOO(1,k)=0;elseif strcmp(f,'L2')HeadODat.SumOO(1,k)=1;elseif strcmp(f,'C1')HeadODat.SumOO(1,k)=2;elseif strcmp(f,'P1')HeadODat.SumOO(1,k)=3;elseif strcmp(f,'P2')HeadODat.SumOO(1,k)=4;elseif strcmp(f,'D1')HeadODat.SumOO(1,k)=5;elseif strcmp(f,'D2')HeadODat.SumOO(1,k)=6;endend%头文件结束elseif strncmp(instrS,'END OF HEADER',13)break;elsecontinue;endend%观测数据结构体%观测数据结构t=0;while ~feof(fid1)%每个历元的第一行数据,时间和观测到的卫星号s=fgets(fid1);t=t+1;year=str2num(s(1:3));month=str2num(s(4:6));day=str2num(s(7:9));hour=str2num(s(10:12));minute=str2num(s(13:15));second=str2num(s(16:26));%历元时间ObsODat(t).TimeOEp=[year,month,day,hour,minute,second];ObsODat(t).TimeOEpp=TimetoJD(year,month,day,hour,minute,second);%该历元观测到的卫星数ObsODat(t).SatSum=str2num(s(30:32));%该历元观测到的卫星号ObsODat(t).SatCode=zeros(1,ObsODat(t).SatSum);ObsODat(t).Obs_FreL1=zeros(1,ObsODat(t).SatSum);ObsODat(t).Obs_FreL2=zeros(1,ObsODat(t).SatSum);ObsODat(t).Obs_RangeC1=zeros(1,ObsODat(t).SatSum);ObsODat(t).Obs_RangeP1=zeros(1,ObsODat(t).SatSum);ObsODat(t).Obs_RangeP2=zeros(1,ObsODat(t).SatSum);for k=1:ObsODat(t).SatSumf=s(31+k*3:32+k*3);ObsODat(t).SatCode(1,k)=str2num(f);end%每个历元的观测数据,按卫星号先后顺序分行存for k=1:ObsODat(t).SatSums=fgets(fid1);%判断一个卫星的观测数据是否分两行存储,如果为两行,则再读入一行if HeadODat.SumOType>5sg=fgets(fid1);s=strcat(s,sg);endL=size(s,2);%补充数据长度if L<HeadODat.SumOType*16s(L+1:HeadODat.SumOType*16)=' ';end%对观测数据判断其类型,并存储到相应的数组中for j=1:HeadODat.SumOTypestemp=s(j*16-15:j*16);stemp=deblank(stemp);if isempty(stemp)continue;endstempNum=str2num(stemp);stempLength=size(stempNum,2);if stempLength>1stempNum=stempNum(1,1);endif HeadODat.SumOO(1,j)==0ObsODat(t).Obs_FreL1(1,k)=stempNum;elseif HeadODat.SumOO(1,j)==1ObsODat(t).Obs_FreL2(1,k)=stempNum;elseif HeadODat.SumOO(1,j)==2ObsODat(t).Obs_RangeC1(1,k)=stempNum;elseif HeadODat.SumOO(1,j)==3ObsODat(t).Obs_RangeP1(1,k)=stempNum;elseif HeadODat.SumOO(1,j)==4ObsODat(t).Obs_RangeP2(1,k)=stempNum;elsecontinue;endend%完成一个卫星的所有观测数据存储end%完成一个历元的数据存储end%完成所有历元的数据存储head=HeadODat;obs=ObsODat;returnfunction EphDat=ReadGpsDataglobal EphDatEphDat=struct;[filename,pathname]=uigetfile('*.**N','读取GPS广播星历文件'); fid1=fopen(strcat(pathname,filename),'rt');if(fid1==-1)msgbox('Input File or Path is not correct','warning ','warn');return;endwhile(1)temp=fgetl(fid1);header=findstr(temp,'END OF HEADER');if(~isempty(header))break;endendi=1;while(1)temp=fgetl(fid1);break;endEphDat(i).PRN=str2double(temp(1:2));year=str2double(temp(3:5));EphDat(i).toc=TimetoJD(year,str2double(temp(6:8)),str2double(temp(9:11)),str2double(temp(12: 14)),str2double(temp(15:17)),str2double(temp(18:22)));EphDat(i).a0=str2num(temp(23:41));EphDat(i).a1=str2num(temp(42:60));EphDat(i).a2=str2num(temp(61:79));temp=fgetl(fid1);EphDat(i).idoe=str2num(temp(4:22));EphDat(i).Crs=str2num(temp(23:41));EphDat(i).dn=str2num(temp(42:60));EphDat(i).M0=str2num(temp(61:79));temp=fgetl(fid1);EphDat(i).Cuc=str2num(temp(4:22));EphDat(i).e=str2num(temp(23:41));EphDat(i).Cus=str2num(temp(42:60));EphDat(i).sqa=str2num(temp(61:79));temp=fgetl(fid1);EphDat(i).toe=str2num(temp(4:22));EphDat(i).Cic=str2num(temp(23:41));EphDat(i).omg0=str2num(temp(42:60));EphDat(i).Cis=str2num(temp(61:79));temp=fgetl(fid1);EphDat(i).i0=str2num(temp(4:22));EphDat(i).Crc=str2num(temp(23:41));EphDat(i).w=str2num(temp(42:60));EphDat(i).omg=str2num(temp(61:79));temp=fgetl(fid1);EphDat(i).iDot=str2num(temp(4:22));EphDat(i).cflgl2=str2num(temp(23:41));EphDat(i).weekno=str2num(temp(42:60));EphDat(i).pflgl2=str2num(temp(61:79));temp=fgetl(fid1);EphDat(i).svacc=str2num(temp(4:22));EphDat(i).svhlth=str2num(temp(23:41));EphDat(i).tgd=str2num(temp(42:60));EphDat(i).aodc=str2num(temp(61:79));temp=fgetl(fid1);EphDat(i).ttm=str2num(temp(4:22));if(i~=1) %删除重复数据if (EphDat(i-1).PRN~=EphDat(i).PRN)break;elseif abs(EphDat(i-1).toc-EphDat(i).toc)<0.001i=i - 1;endendendi=i + 1;endfunction DDDWformat longclear allticglobal HeadODat;global ObsODat;global EphDat;%先读N文件,再读O文件EphDat=ReadGpsData;[HeadODat,ObsODat]=ReadObsData;%多个历元加权平均求测站点坐标N = size(EphDat,2);c=299792458;epochNUM=size(ObsODat,2);%观测数据的个数Xk0=HeadODat.ApproXYZ(1,1); %测站点的近似坐标Yk0=HeadODat.ApproXYZ(1,2);Zk0=HeadODat.ApproXYZ(1,3);for k=1:epochNUMtr=ObsODat(k).TimeOEp;%历元接收数据时间Snum=ObsODat(1,k).SatSum; %该历元观测到的卫星数if Snum<4break; %去除无解的历元endCode=ObsODat(k).SatCode;%该历元观测到的卫星号组R=ObsODat(k).Obs_RangeC1 ; %取出C1观测值,列向量%卫星发射时间的迭代计算for j=1:Snumif R(j) == 0continue;endt=TimetoJD(tr(1),tr(2),tr(3),tr(4),tr(5),tr(6));t2 = mod(t,7)*24*3600;%gps second%卫星钟差dd=-1;for i=1:Nif(EphDat(i).PRN==Code(j)&abs(t-EphDat(i).toc)<0.0417*2)%读入时间相近的星历文件a0= EphDat(i).a0;a1=EphDat(i).a1;a2= EphDat(i).a2;toe=EphDat(i).toe;dt=a0+(t2-toe)*a1+(t2-toe)^2*a2;%卫星钟差dd=1;break;endendif dd==-1msgbox('没有相关星历文件');return;endtt = tr;ts = tr(6) - R(j)/c;%用秒进行迭代sign = 1;while(sign>1E-8)tt(6) = ts;[Xs1,Ys1,Zs1]= CalPos(Code(j),tt);ss1 = sqrt((Xk0-Xs1)^2+(Yk0-Ys1)^2+(Zk0-Zs1)*(Zk0-Zs1));%卫星位置加地球自转改正qe=0.000072921151467; %地球自转角速度delt=ss1/c;Rotation=[cos(qe*delt),sin(qe*delt),0;-sin(qe*delt),cos(qe*delt),0;0,0,1];h=Rotation*[Xs1;Ys1;Zs1] ;Xs=h(1);Ys=h(2);Zs=h(3);ss = sqrt((Xk0-Xs)^2+(Yk0-Ys)^2+(Zk0-Zs)*(Zk0-Zs));ts = tr(6)-ss/c;sign = norm(ts-tt(6));endaxk = (Xk0-Xs)/ss;ayk = (Yk0-Ys)/ss;azk = (Zk0-Zs)/ss;EAB=CAL2POL(Xk0,Yk0,Zk0,Xs,Ys,Zs);if j==1A = [axk,ayk,azk,1];L = R(j)-ss++c*dt-2.47/(sin(EAB)+0.0121);elseA = [A;axk,ayk,azk,1]; %构造误差方程L = [L;(R(j)-ss++c*dt)-2.47/(sin(EAB)+0.0121)];%常数项endendif Snum==4dx=inv(A)*L;aaaa(k).Q=inv(A'*A);Px(k) = 1/aaaa(k).Q(1,1);Py(k) = 1/aaaa(k).Q(2,2);Pz(k) = 1/aaaa(k).Q(3,3);xchange(k) = dx(1);ychange(k) = dx(2);zchange(k) = dx(3);elseif Snum > 4dx = inv(A'*A)*A'*L; %构造法方程并求解aaaa(k).Q=inv(A'*A);Px(k) = 1/aaaa(k).Q(1,1);Py(k) = 1/aaaa(k).Q(2,2);Pz(k) = 1/aaaa(k).Q(3,3);xchange(k) = dx(1);ychange(k) = dx(2);zchange(k) = dx(3);endendXp=sum(Px.*(Xk0+xchange))/sum(Px)Yp=sum(Py.*(Yk0+ychange))/sum(Py)Zp=sum(Pz.*(Zk0+zchange))/sum(Pz)figure(1);subplot(3,1,1);plot(xchange,'black--');subplot(3,1,2);plot(ychange,'black--');subplot(3,1,3);plot(zchange,'black--');tocfunction [x,y,z]= CalPos(PRN,time)global EphDatt1=TimetoJD(time(1),time(2),time(3),time(4),time(5),time(6));t2=TimetoJD(time(1),time(2),time(3),0,0,0);if isempty(EphDat)EphDat=ReadGpsData;endsz=size(EphDat,2);gg=0;%判断最近的时间for i=1:szif(EphDat(i).PRN==PRN & abs(EphDat(i).toc-t1)<0.0417*2)G=EphDat(i);gg=1;break;endendt3=t2-G.weekno*7;% 减整周数t=t3*24*3600+time(4)*3600+time(5)*60+time(6);%化为GPS秒u=3.986004418E+14; %地球引力常数we=7.2921151467E-5; %地球自转速度a=G.sqa^2; %地球长半轴n0=sqrt(u/a^3); %计算的平运动tk=t-G.toe; %从参考历元开始的时间if tk>=302400tk=tk-604800;elseif tk<-302400tk=tk + 604800;endn=n0+G.dn; %改正后的运动Mk=G.M0+n*tk;Ek=Mk; %偏近点角弧度for i=1:4Ek=Mk+G.e*sin(Ek);endfk=2*atan((sqrt((1+G.e)/(1-G.e))*tan(Ek/2))); %真近点角if(fk<0)fk=fk+2*pi;endcosf=(cos(Ek)-G.e)/(1-G.e*cos(Ek));sinf=(sqrt(1-G.e^2)*sin(Ek))/(1-G.e*cos(Ek));uk=fk+G.w; %升交角矩及轨道摄动改正项iuk=G.Cuc*cos(2*uk)+G.Cus*sin(2*uk);irk=G.Crc*cos(2*uk)+G.Crs*sin(2*uk);iik=G.Cic*cos(2*uk)+G.Cis*sin(2*uk);uk=uk+iuk ; %改正后的升角距r=a*(1-G.e*cos(Ek))+irk; %改正后的地心相径i=G.i0+iik+G.iDot*tk; %改正后的倾角xx=r*cos(uk); %在轨道面中的位置yy=r*sin(uk);omg=G.omg0+(G.omg-we)*tk-we*G.toe; %改正后的升交点经度x=xx*cos(omg)-yy*cos(i)*sin(omg);y=xx*sin(omg)+yy*cos(i)*cos(omg);z=yy*sin(i);CalPos=[x,y,z];%打印结果% x=num2str(x,'%12.6f');% y=num2str(y,'%12.6f');% z=num2str(z,'%12.6f');% fprintf('卫星号PRN:%2.0f',PRN); fprintf('\n');% fprintf('时间time:%4.0f%4.0f%4.0f%4.0f%4.0f%4.1f',time); fprintf('\n');% fprintf('所求卫星在地心坐标系中的空间坐标为:\n');% fprintf('X = ');fprintf(x,'\n');fprintf('m\n');% fprintf('Y = ');fprintf(y,'\n');fprintf('m\n');% fprintf('Z = ');fprintf(z,'\n');fprintf('m\n');%returnfunction EAB=CAL2POL(X,Y,Z,XB,YB,ZB)format longL=atan(Y/X);R=sqrt(X*X+Y*Y+Z*Z);a=6378137;e=sqrt(0.0066943799013);seta=atan(Z/sqrt(X*X+Y*Y));B=asin(sqrt((a*a-R*R*cos(seta)*cos(seta))/(a*a-R*R*cos(seta)*cos(seta)*e*e)))B1=atan(tan(seta)*(1+a*e*e*sin(B)/Z/sqrt(1-e*e*sin(B)*sin(B))));while abs(B1-B)>0.01*pi/180/3600B=B1;B1=atan(tan(seta)*(1+a*e*e*sin(B)/Z/sqrt(1-e*e*sin(B)^2)));endH=R*cos(seta)/cos(B)-(a/sqrt(1-e*e*sin(B)*sin(B)));L=L*180/pi;B=B*180/pi;DX=XB-X;DY=YB-Y;DZ=ZB-Z;NB=-sin(B)*cos(L)*DX-sin(B)*sin(L)*DY+cos(B)*DZ; EB=-sin(L)*DX+cos(L)*DY;UB=cos(B)*cos(L)*DX+cos(B)*sin(L)*DY+sin(B)*DZ; SAB=sqrt(NB*NB+UB*UB+EB*EB);EAB=asin(UB/SAB);。