基于RENIX3.02格式读取卫星观测数据O文件MATLAB源代码

function[APXYZDeltaHENsignal_typeobsdata]=read_obs(filename)%%%function[signal_typeobsdata]=read_obs(filename)%程序描述：读取观测数据文件，rinex格式%---------------------------------------------------------------%--Inputs:%--filename%--Outputs:%--obsdata%--observationmessage%---modifiedbyZhang2017-12-1%---------------------------------------------------------------%-----Openthefile-------%fid=fopen(filename);disp(['读取卫星观测数据文件......:',filename]);%[filename,pathname]=uigetfile('*.**o','读取卫星观测数据文件');%fid=fopen(strcat(pathname,filename),'rt');%if(fid==-1)%error('errortoopenthe'+filename);%end%%-----读取头部-----%hang=zeros(1,80);head_lines=0;while1head_lines=head_lines+1;hang=fgetl(fid);head(head_lines).str=str2mat(hang);if(~isempty(findstr(hang,'APPROXPOSITIONXYZ')))%obtainapproximatemarkerposition(WGS-84)，获得观测点大概位置（WGS-84）APXYZ=[0,0,0];APXYZ(1)=str2num(hang(1:14));APXYZ(2)=str2num(hang(15:28));APXYZ(3)=str2num(hang(29:42));endif(~isempty(findstr(hang,'ANTENNA:DELTAH/E/N')))%obtainAntennaheight、Eccentricitiesofantennacenterrelativetomarkertotheeastandnorth(allunitsinmeters)%天线高度：天线底部相交于标记点的高度，天线中心相对标记点东向和北向距离（单位：米）DeltaHEN=[0,0,0];DeltaHEN(1)=str2num(hang(1:14));DeltaHEN(2)=str2num(hang(15:28));DeltaHEN(3)=str2num(hang(29:42));endif(~isempty(findstr(hang,'WAVELENGTHFACTL1/2')))%获取L1和L2载波的缺省系数，1:?整周；2:半周；0(L2载波):单频接收机L1_Wavalength=str2num(hang(1:6));L2_Wavalength=str2num(hang(7:12));endif(~isempty(findstr(hang,'#/TYPESOFOBSERV')))%文件中观测类型数量%格式例如：21L1P1C1L2P2D1D2S1S2#/TYPESOFOBSERV%L5C5D5S5L7C7D7S7L8#/TYPESOFOBSERV%C8D8S8#/TYPESOFOBSERVtype=str2num(hang(1:6));if(type=9)%观测类型数量少于9个fort=1:typesignal_type(t).types=hang(6*t+5:6*t+6);endelse%观测类型数量多于9个，分多行显示fort=1:9%由于第1行已经读取，先存储第1行signal_type(t).types=hang(6*t+5:6*t+6);endform=1:floor(type/9)%再循环存储其他行head_lines=head_lines+1;hang=fgetl(fid);head(head_lines).str=str2mat(hang);fort=9*m+1:9*m+9if(t=type)signal_type(t).types=hang(6*(t-9*m)+5:6*(t-9*m)+6);endendendendend%文件中观测类型数量读取完毕，相应说明：%TYPESOFOBSERV：在本数据文件中所存储不同观测值类型的数量和观测值类型列表，如果超过9种观测值类型,则使用续行。%说明：%L1,L2：L1和L2上的相位观测值%C1：采用L1上C/A码所测定的伪距；%P1,P2：采用L1、L2上的P码所测定的伪距；%D1,D2：L1和L2上的多普勒频率；%T1,T2：子午卫星的150(T1)和400MHz(T2)信号上的多普勒积分；%S1,S2：接收机所给出的L1、L2相位观测值的原始信号强度或SNR值。%观测值的单位:载波相位为周,伪距为m,多普勒为Hz,子午卫星为周,SNR等则与接收机有关。if(~isempty(findstr(hang,'INTERVAL')))%观测时间间隔，单位：秒interval=str2num(hang(1:6));endif(~isempty(findstr(hang,'ENDOFHEADER')))break;endend%ifhead(1).str(21)~='O'%error('Thisfileisnotanavgationmessagefile.Pleaseseleteacorrectfile');%elseifhead(1).str(41)~='G'%error('ThisisnotaGPSsystem');%end%%--------读取观测数据-------------data_num=0;while1%data_num=data_num+1;%-----Firstlineofobservationdata-----hang=fgetl(fid);if(hang==-1)%endoffilebreak;enda=blanks(80);a(1:length(hang))=hang;hang=a;disp(hang);if(isstrprop(hang(2),'digit')==1)data_num=data_num+1;%--------Epoch/GPStimesystem--------obsdata(data_num).year=2000+str2double(hang(1:3));obsdata(data_num).month=str2num(hang(4:6));obsdata(data_num).day=str2num(hang(7:9));obsdata(data_num).hour=str2num(hang(10:12));obsdata(data_num).min=str2num(hang(13:15));obsdata(data_num).sec=str2num(hang(16:26));%-----------Epochflag星历标志------------------%--0:OK%--1:powerfailurebetweenpreviousandcurrentepoch从前一历元到当前历元观测失败?%--1:Eventflag事件标志obsdata(data_num).Epoch_flag=str2num(hang(27:29));%--Nemberofsatellitesincurrentepoch当前历元的卫星数量------obsdata(data_num).sat_num=str2num(hang(30:32));%--ListofPRNs(sat.numbers)incurrentepoch当前观测到的卫星列表---%例如：171130000.0000000025G05G08G10G13G15G18G20G21G24G28G30R04%R05R06R07R09R14R15R16R22R23E02E07E08%E30ifobsdata(data_num).sat_num=12%当前观测到的卫星如果少于12颗forl=0:obsdata(data_num).sat_num-1obsdata(data_num).sat(l+1)=hang(33+l*3);%卫星系统，如Gobsdata(data_num).prn(l+1)=str2double(hang(34+l*3:35+l*3));%卫星编号，如05endelse%当前观测到的卫星如果多于12颗forl=0:11%先存储第1行12颗卫星编号obsdata(data_num).sat(l+1)=hang(33+l*3);obsdata(data_num).prn(l+1)=str2double(hang(34+l*3:35+l*3));endform=1:(ceil(obsdata(data_num).sat_num/12)-1)%再循环存储其他行hang=fgetl(fid);a=blanks(80);a(1:length(hang))=hang;hang=a;forl=0:11if((l+12*m+1)=obsdata(data_num).sat_num)obsdata(data_num).sat(l+12*m+1)=hang(33+l*3);obsdata(data_num).prn(l+12*m+1)=str2num(hang(34+l*3:35+l*3));endendendend%----------------------Observations读取卫星观测数值，每颗卫星的数值个数由type观测类型数量决定-------------------------%%观测值格式：m(F14.3,I1,I1)，观测值占14位，带3位小数；LLI占1位，范围:?0-7；信号强度占1位，设为?1-9级%每个观测值占16位，每行5个观测值forl=1:obsdata(data_num).sat_numhang=fgetl(fid);a=blanks(80);a(1:length(hang))=hang;hang=a;if(hang==-1)%endoffilebreak;endiftype=5%观测类型数量少于5个fort=1:typeifisspace(hang(16*(t-1)+1:16*(t-1)+14))==[11111111111111]%isspace：测试参数是否是空格obsdata(data_num).PRN(l).record(t)=0;elseobsdata(data_num).PRN(l).record(t)=str2num(hang(16*(t-1)+1:16*(t-1)+14));endifisspace(hang(16*(t-1)+15))==1obsdata(data_num).PRN(l).LLI(t)=0;elseobsdata(data_num).PRN(l).LLI(t)=str2num(hang(16*(t-1)+15));endifisspace(hang(16*(t-1)+16))==1obsdata(data_num).PRN(l).Signal_strenth(t)=0;elseobsdata(data_num).PRN(l).Signal_strenth(t)=str2num(hang(16*(t-1)+16));endendelsefort=1:5%先存储第一行ifisspace(hang(16*(t-1)+1:16*(t-1)+14))==[11111111111111