61549829GPS.zip

%--------－----函数：读取观测值文件的函数-------------------------------------------------- function [pObsHeader,pObsText,numEpoch]=ReadObsFile(file) format long; %------------------------- % numEpoch 有效历元数 % pObsHeader.nTypeNum %观测值类型个数（如C1，P1，P2，L1，L2等） % pObsHeader.typeObs %typeObs结构体为typeObs.C1、typeObs.P1、typeObs.P2、typeObs.L1、typeObs.L2、typeObs.D1、typeObs.D2 % pObsHeader.name 测站名 % pObsHeader.X0、pObsHeader.Y0、pObsHeader.Z0 %测站近似坐标 % pObsHeader.flag_Rdeltat %历元时标、码伪距、相位是否使用实时确定出的接收机钟差改正，0=否，1=是，默认0 % pObsText(numEpoch+1).epochFlag %历元标志 % pObsText(numEpoch+1).svNumOfEpoch %历元卫星数 % pObsText(numEpoch+1).GPSymd %历元年月日 % pObsText(numEpoch+1).GPSweekno %GPS周 % pObsText(numEpoch+1).GPSweeksec %周秒 % pObsText(numEpoch+1).Hour %小时 % pObsText(numEpoch+1).Min %分 % pObsText(numEpoch+1).Second %秒 % pObsText(numEpoch+1).svPrnOfEpoch(i) %历元第i颗卫星的卫星号 % pObsText(numEpoch+1).Rdeltat %历元接收机钟差改正 % pObsText(numEpoch+1).data(i,j) %历元第i颗卫星第j类观测值的数值 % pObsText(numEpoch+1).LLI(i,j) %失锁标识符，0和空格表示正常，见P249，10.1节Rinex格式 表10-3 % pObsText(numEpoch+1).SN(i,j) %信号强度，即信噪比，1表示最小信号强度，9表示最大信号强度，5表示良好信噪比的阈值，0或空格表示未给出，见P249，10.1节Rinex格式 表10-3 %------------------------- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%打开观测文件%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% fid = fopen(file); %%%fopen(file)表示打开文件名为file的文件 if fid<=0 %%%文件打开成功返回大于0的数 errordlg('文件打开错误','Open File Error') else tline=fgetl(fid); %开始读取文件的第一行 如 tline= "2.10 OBSERVATION DATA G (GPS) RINEX VERSION / TYPE" end if (strfind(tline,'OBSERVATION')>0)&&(strfind(tline,'RINEX')>0) %判断是否为观测值文件 flagReadData=0; %用flagReadData标志读文件，0表示文件为观测值文件，且开始读头文件部分,-1表示文件读完，1表示头文件读完开始读数据部分，2表示打开的不是观测值文件 else errordlg('打开的不是观测值文件','Open File Error') flagReadData=2; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%开始读观测文件rinex的头文件，进入while循环，若flagReadData==0 ,表示程序正在读头文件%%%%%%%%%%%%%%%% while flagReadData==0 tline=fgetl(fid);%读文件的一行，将读的一行作为字符串存到tline变量中，fid指向第二行 if strfind(tline,'TYPES OF OBSERV')>0 %%%表示该行在描述观测值中包含的观测值类型 typeObs.C1 = -1; typeObs.P1 = -1;typeObs.P2 = -1; typeObs.L1 = -1;typeObs.L2 = -1; %所有类型观测值取值为-1，说明没有被观测，当下面检查到存在该类型观测值，将值改为大于0的数 typeObs.D1 = -1;typeObs.D2 = -1; typeObs.S1 = -1;typeObs.S2 = -1; %%%以下代码用以获取有几种观测值，用pObsHeader.nTypeNum来存储观测值类型的个数； 每种观测值出现的顺序 [s,tline] = strtok(tline); % 5 L1 L2 C1 P1 P2 # / TYPES OF OBSERV 分离该行 s=5, tline= L1 L2 C1 P1 P2 nTypeNum = str2num(s); %得到O文件中的观测类型个数 如该例中有5种类型的观测值 L1 L2 C1 P1 P2 for i=1:nTypeNum [s,tline] = strtok(tline); if strcmp(s,'C1')>0 typeObs.C1 = i; elseif strcmp(s,'P1')>0 typeObs.P1 = i; elseif strcmp(s,'P2')>0 typeObs.P2 = i; elseif strcmp(s,'L1')>0 typeObs.L1 = i; elseif strcmp(s,'L2')>0 typeObs.L2 = i; elseif strcmp(s,'D1')>0 typeObs.D1 = i; elseif strcmp(s,'D2')>0 typeObs.D2 = i; elseif strcmp(s,'S1')>0 typeObs.S1 =i; elseif strcmp(s,'S2')>0 typeObs.S2 =i; end end pObsHeader.nTypeNum=nTypeNum; pObsHeader.typeObs=typeObs; end if strfind(tline,'APPROX POSITION XYZ')>0 %表明该行描述的是测站近似坐标 0.0000 0.0000 0.0000 APPROX POSITION XYZ [s,tline] = strtok(tline); pObsHeader.X0 = str2num(s); [s,tline] = strtok(tline); pObsHeader.Y0 = str2num(s); [s,tline] = strtok(tline); pObsHeader.Z0 = str2num(s); end if strfind(tline,'MARKER NAME')>0 %获取测站名 [s,tline] = strtok(tline); pObsHeader.name= s; end if strfind(tline,'RCV CLOCK OFFS APPL')>0 [s,tline] = strtok(tline); pObsHeader.flag_Rdeltat = str2num(s);%历元时标、码伪距、相位是否使用实时确定出的接收机钟差改正，0=否，1=是，默认0 end if strfind(tline,'END OF HEADER')>0 %%%%%%%%%%%%%%%%%%%%%%%表示头文件已读完 flagReadData=1; end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%flagReadData==1，开始读数据文件%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% numEpoch=0;%用来记录历元个数 while flagReadData==1 tline=fgetl(fid); %%%%读取数据文件的第一行 if tline==-1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%表示文件读完 flagReadData=-1; %%%%给读取文件标志赋-1 elseif(length(tline)>=10)&&(strcmp(tline(1:10),' ')~=1)%%%%%%%%%%%%%%%%%%%%%%如果文件没有读完且读的字符串是观测数据 %%%%%%%%%开始读取rinex数据的数据头，即对历元的时间，状态，编号等的描述 %%%%读RINEX数据头为04 1 2 0 0 30.0000000 0 8G 2G25G 1G16G 3G 6G20G14 tline=strrep(tline,'G',' '); [s,tline] = strtok(tline); iYear = str2num(s); s1 = strcat(s,'/'); [s,tline] = strtok(tline); iMonth = str2num(s); s1 = strcat(s1,s,'/'); [s,tline] = strtok(tline); iDay = str2num(s); s1 = strcat(s1,s,'/'); [s,tline] = strtok(tline); iHour = str2num(s); s1 = strcat(s1,s,':'); [s,tline] = strtok(tline); iMin = str2num(s); s1 = strcat(s1,s,':'); [s,tline] = strtok(tline); fSecond = str2num(s); s1 = strcat(s1,s); [s,tline] = strtok(tline); pObsText(numEpoch+1).epochFlag = str2num(s); %历元标志，0表示信号正常，其他表示信号有中断等 [s,tline] = strtok(tline); pObsText(numEpoch+1).svNumOfEpoch = str2num(s); %表示该历元观测到的卫星数 if iYear<50 iYear=iYear+2000; end if iYear>50&&iYear<100 iYear=iYear+1900; end if fSecond==60 fSecond=0; iMin=iMin+1; end if iMin==60 iMin=0; iHour=iHour+1; end if iHour==24 iHour=0; iDay=iDay+1; end [weekno,weeksec,dayofy] =TimeConvertWeekTime(iYear,iMonth,iDay,iHour,iMin,fSecond);%%%%%%%%%%%给出年月日，时分秒，求GPS周，GPS秒等 pObsText(numEpoch+1).GPSweekno = weekno; pObsText(numEpoch+1).GPSweeksec = weeksec; pObsText(numEpoch+1).GPSymd = s1; pObsText(numEpoch+1).Hour = iHour; pObsText(numEpoch+1).Min = iMin; pObsText(numEpoch+1).Second = fSecond; if pObsText(numEpoch+1).svNumOfEpoch>12 %如果该历元的接收到的卫星数目大于12，记录内容要换行 for i=1:12 [s,tline] = strtok(tline); pObsText(numEpoch+1).svPrnOfEpoch(i) = str2num(s); end if isempty( str2num(tline) ) ~= 1 pObsText(numEpoch+1).Rdeltat = str2num(tline); % 接收机钟差改正 else pObsTe