基于matlab的GPS单点定位程序开发

开发初衷

本人是测绘专业的学生，在初学GPS时遇到了很多代码问题，因此希望这篇文章能够帮助广大的测绘和导航专业的学子

代码开始

开发的第一部分是：文件读取,文件读取包括两部分，观测文件和星历文件读取。

首先介绍第一部分文件读取：读取星历文件（代码如下）

读取GPS星历（可以是混合星历）

// An highlighted blockfunction eph=loadxingli(filename)%%读取rinex3的星历文件tic;fprintf('Opening RINEX file: %s\n', filename);finp = fopen(filename, 'rt');
%%%%星历索引值no = 2:3;           year_idx = 5:8;        month_idx = 10:11;    day_idx = 13:14;hour_idx = 16:17;   min_idx = 19:20;       sec_idx = 22:23;      f1 = 5:23;          f2 = 24:42;            f3 = 43:61;           f4 = 62:80;
%%%跳过文件头
while 1line=fgetl(finp);if contains(line,'END OF HEADER')break;end
end
%%预设数组
sat=1:32;nav=cell(1,32);
eph=[];SEP=zeros(size(sat));
block_init =zeros(42,1);
while ~feof(finp)%%每次循环8行line=fgetl(finp);system=line(1);if system=='G'prn = str2double(line(no));idx = find(prn ==sat);SEP(idx) = SEP(idx) + 1;tt = zeros(6,1);tt(1) = round(str2double(line(year_idx)));tt(2) = round(str2double(line(month_idx)));tt(3) = round(str2double(line((day_idx))));tt(4) = str2double(line(hour_idx));tt(5) = str2double(line(min_idx));tt(6) = str2double(line(sec_idx));[GPSWeekNo, GPSSecond, DOY, DOW] = greg2gps(tt');mTime = datenum(tt');block=block_init;block(1:11) = [tt; GPSWeekNo; GPSSecond; DOY; DOW; mTime];%%%第一行block(12:14) = [str2double(line(f2));str2double(line(f3));str2double(line(f4))];%%%第二行line=fgetl(finp);block(15:18) = [str2double(line(f1));str2double(line(f2));...str2double(line(f3));str2double(line(f4))];%%%第三行line=fgetl(finp);block(19:22) =  [str2double(line(f1));str2double(line(f2));...str2double(line(f3));(str2double(line(f4)))^2];%%%% Line 4line=fgetl(finp);block(23:26) = [str2double(line(f1));str2double(line(f2));...str2double(line(f3));str2double(line(f4))];%%%% Line 5line=fgetl(finp);block(27:30) =  [str2double(line(f1));str2double(line(f2));...str2double(line(f3));str2double(line(f4))];%%%% Line 6line=fgetl(finp);block(31:34) =  [str2double(line(f1));str2double(line(f2));...str2double(line(f3));str2double(line(f4))];        %%%% Line 7line=fgetl(finp);block(35:38) =  [str2double(line(f1));str2double(line(f2));...str2double(line(f3));str2double(line(f4))];%%%% Line 8(最后两个可以直接设置为NaN，如果报错)line=fgetl(finp);block(39:42) =  [str2double(line(f1));str2double(line(f2));...str2double(line(f3));str2double(line(f4))];  %%%存到数组中nav{1,prn}(:,SEP(idx))=block;end
end
toc;
eph=nav;
%%%%星历读取完成

关于子程序（greg2gps可以私聊我获得）

datenum函数是2018版matlab自带的

读取观测文件（可以是混合文件）

%%读取rinex3的观测文件，并完成定位
%%%观测文件的结构体，不启用
%RNX=struct('obs_body','approx');fprintf('open obs file; %s\n',filename);fid=fopen(filename,'rt');lineidx=0;P1=6:17;L1=21:33;S1=44:49;fprintf('open obs file; %s\n',filename);
fid=fopen(filename,'rt');
lineidx=0;
while 1line=fgetl(fid);lineidx=lineidx+1;if contains(line,'APPROX POSITION XYZ')app_xyz=sscanf(line,'%.4f %.4f %.4f');%计算高度角，随机模型(强调不为零)endif contains(line,'END OF HEADER')break;endend
%%%读取文件主体
while ~feof(fid)line=fgetl(fid);first_str=line(1);switch first_strcase '>'lineidx=lineidx+1;help=0;help_satpos=[];sat_num=str2double(line(34:35));%系数矩阵xishu=zeros(sat_num,4);Z  =zeros(sat_num,1);P  =eye(sat_num);xishu_idx=0;time(1)=str2double(line(3:6));time(2)=str2double(line(8:9));time(3)=str2double(line(11:12));time(4)=str2double(line(14:15));time(5)=str2double(line(17:18));time(6)=str2double(line(20:29));JD = juliandate(datetime(time(1),time(2),time(3)));gpsweek = floor((JD - 2444244.5)/7);dow = JD - (2444244.5 + gpsweek*7);tow = dow*86400 + time(4)*3600 + time(5)*60 + time(6);doy = JD - juliandate(time(:,1),1,1) + 1;t_time=[gpsweek, tow, doy, dow];case 'G'lineidx=lineidx+1;sat_id=str2double(line(2:3));idx=find(sat==sat_id);sep(idx)=sep(idx)+1;obs1=str2double(line(P1));obs2=str2double(line(L1));obs3=str2double(line(S1));xishu_idx=xishu_idx+1;obs=[t_time(1) t_time(2) obs1 obs2 obs3];%%%%%%% 这一部分是选取合适的星历，这是是本程序中需要重点改进的部分eph=nav{1,sat_id}';eph_time=eph(:,8);help=0for i=1:length(eph_time)help=help+1;if eph_time(i)>towbreak;endend  %该方法强调星历必须顺序存储,且保证不能超过两小时eph_idx=help;% eph_idx= eph_idx_ini+fix(abs(obs(2)-eph(help,8))/7200);%两小时星历%%%%%%%%（能够完成自动识别，而不是手动改正）%%%%卫星传播时间矫正time_s=tow-obs1/299792458;            af2 = eph(eph_idx,14);weekno = eph(eph_idx,33);af0 = eph(eph_idx,12);toc    = eph(eph_idx,8);af1 = eph(eph_idx,13);timeeph= weekno*7*86400 + toc;dt  =  check_t(time_s - timeeph);corr = (af2 * dt + af1) * dt + af0;obs(2)=time_s-corr;%%%%%时间矫正完成xyz =getSatPosGPS(obs(1:2),eph(eph_idx,:));help_satpos=[help_satpos;xyz];%%%%构建法方程，可以用最小二乘也可以使用卡尔曼滤波%%卡尔曼滤波读者自己开发，最小二乘需要迭代dist=norm(xyz,app_xyz);xishu(xishu_idx,:)=[-(xyz(1)-app_xyz(1))/dist,... -(xyz(2)-app_xyz(2))/dist,-(xyz(3)-app_xyz(3))/dist,-1];Z(xishu_idx)=obs1-dist;%未考虑电离层，对流层等影响%%%计算卫星高度角temp_enu=(xyz2enu(xyz,app_xyz))';Elevation = atan2(abs(temp_enu),norm(temp_enu(1,1:2)));P(xishu_idx)=0.3*0.3/sind(Elevation);     %%%%if xishu_idx==sat_num%%%历元数据遍历完成delta=(xishu'*P*xishu)\(xishu'*P*Z)X=app_xyz+delta;while 1for i=1:xishu_idxxyz=help_satpos(i,:);dist=norm(xyz,X);xishu(i,:)=[-(xyz(1)-X(1))/dist,... -(xyz(2)-X(2))/dist,-(xyz(3)-X(3))/dist,-1];Z(i)=obs1-dist;enddelta=(xishu'*P*xishu)\(xishu'*P*Z)X=app_xyz+delta;        if norm(X-app_xyz)<0.1%迭代终止条件break;   %iteration end for this epochendapp_xyz=X;endendendend

关于子程序check_t、xyz2enu可以私聊我获得， juliandate和datetime是18版matlab自带的子函数

计算卫星位置（可以参考GPS数据处理）

function pos = getSatPosGPS(GPStime,eph)TOE      = eph(23);a        = eph(22);ecc      = eph(20);Delta_n  = eph(17);M0       = eph(18);Omega    = eph(25);DOmega   = eph(30);omega    = eph(29);i0       = eph(27);Di       = eph(31);CRS      = eph(16); CRC      = eph(28); CIS      = eph(26); CIC      = eph(24); CUS      = eph(21); CUC      = eph(19);%%%常量const.GM  = 3.986005e14;const.wE  = 7.2921151467e-5;const.a   = 6378137;const.e   = 0.0818191908426215;const.f   = 1/298.257223563;%tol_Kepler = 0.001;tk = GPStime(:,2) - TOE;tk(tk < 0) = tk(tk < 0) + 604800;n0 = sqrt(const.GM*a^-3);n = n0 + Delta_n;Mk = mod(M0 + n*tk, 2*pi);Ek = Mk;for d = 1:length(Mk)Ek(d) = kepler(ecc, Mk(d), tol_Kepler);endEk = mod(Ek, 2*pi);if max(abs(Ek - ecc.*sin(Ek) - Mk)) > tol_Kepler *(pi/648000)%打出来变成斜体error('Kepler equation not computed correct values !!!');endvk = mod(atan2(sqrt(1 - ecc^2)*sin(Ek), cos(Ek) - ecc), 2*pi);r0k = a.*(1 - ecc*cos(Ek));Phik = mod(vk + omega, 2*pi)delta_uk = CUC.*cos(2*Phik) + CUS.*sin(2*Phik);delta_rk = CRC.*cos(2*Phik) + CRS.*sin(2*Phik);delta_ik = CIC.*cos(2*Phik) + CIS.*sin(2*Phik);uk = mod(Phik + delta_uk, 2*pi);rk = r0k + delta_rk;ik = i0 + Di*tk + delta_ik;lamk = Omega + (DOmega - const.wE)*tk - const.wE*TOE;lamk = mod(lamk, 2*pi);pos = [rk.*(cos(lamk).*cos(uk) - sin(lamk).*sin(uk).*cos(ik)),...rk.*(sin(lamk).*cos(uk) + cos(lamk).*sin(uk).*cos(ik)),...rk.*sin(uk).*sin(ik)];

子程序kelper可以私聊我获得

基于matlab的GPS单点定位程序开发（初学者）相关推荐

《基于Windows 7特性的程序开发系列》视频分享
前一阵录制了<基于Windows 7特性的程序开发系列>视频课程,主要针对WinForm.WPF 开发具有Windows 7 特性的程序.现已发布到MSDN Webcast 欢迎大家拍砖. ...
基于VTK的MFC应用程序开发(3)
基于VTK的MFC应用程序开发(3) 分类: VTK应用示例 2013-05-17 13:37 3307人阅读评论(23) 收藏举报目录(?)[+] 之前介绍了基于VTK的单文档应用程序开发,并 ...
基于VTK的MFC应用程序开发(2)
基于VTK的MFC应用程序开发(2) 分类: VTK应用示例 2013-03-29 13:03 6647人阅读评论(18) 收藏举报 MFCVTK图像重采样目录(?)[+] 现在基于VTK的MF ...
基于VTK的MFC应用程序开发(1)
基于VTK的MFC应用程序开发(1) 分类: VTK应用示例 2013-03-24 22:35 4195人阅读评论(28) 收藏举报 CMakeVTKMFC 目录(?)[+] 提到MFC一般都不陌 ...
车牌识别与计算机编程,基于MATLAB的车牌识别程序详解.ppt
基于MATLAB的车牌识别程序详解自定义一个字符函数,用来从车牌区域中提取出7个字符,其中利用切割函数来进行切割. 程序:function [word,result]=getword(d) word ...
小程序 | 基于WAMP的新闻网小程序开发（体验全栈式开发微信小程序）
之前学习微信小程序开发,主要是基于JS.WXML.WXSS的前端开发,对于后端技术不精的我也是使用了微信开发者工具中的云开发功能,但是今天突发奇想,特别想体验一下全栈式开发微信小程序,学习了一下基于W ...
基于matlab的gps信号仿真123,MATLABGPS信号仿真完整源代码.doc
配套毕业设计论文见百度文库请搜索 <基于MATLAB的GPS信号仿真123> 附录仿真程序代码数据码的产生 function datacode=data(x) y=rand(1,x) ...
基于matlab的指纹识别程序
提示:文章写完后,目录可以自动生成,如何生成可参考右边的帮助文档基于matlab的指纹识别程序前言一.程序思路是什么? 二.预处理步骤 1.指纹图像的灰度化处理 2.指纹图像的归一化与分割处理 ...
matlab相关性分析频谱_基于Matlab的相关频谱分析程序教程
基于Matlab的相关频谱分析程序教程 Matlab 信号处理工具箱谱估计专题频谱分析 Spectral estimation(谱估计)的目标是基于一个有限的数据集合描述一个信号的功率(在频率上的 ...

基于matlab的GPS单点定位程序开发（初学者）

开发初衷

代码开始

读取GPS星历（可以是混合星历）

关于子程序（greg2gps可以私聊我获得）

读取观测文件（可以是混合文件）

关于子程序check_t、xyz2enu可以私聊我获得， juliandate和datetime是18版matlab自带的子函数

计算卫星位置（可以参考GPS数据处理）

子程序kelper可以私聊我获得

基于matlab的GPS单点定位程序开发（初学者）相关推荐

最新文章

热门文章