基于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]%isspace：測試參數是否是空格obsdata(data_num).PRN(l).record(t)=0;elseobsdata(data_num).PRN(l).record(t)=str2double(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)=str2double(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)=str2double(hang(16*(t-1)+16));endendform=1:floor(type/5)%再循環存儲其他行hang=fgetl(fid);a=blanks(80);a(1:length(hang))=hang;hang=a;forn=