多系统精密星历下载与分析
- 概述
精密星歷是由若干衛星跟蹤站的觀測數據,經事后處理算得的供衛星精密定位等使用的衛星軌道信息。
IGS精密星歷采用sp3格式,其存儲方式為ASCII文本文件,內容包括表頭信息以及文件體,文件體中每隔15 min給出1個衛星的位置,有時還給出衛星的速度。它的特點就是提供衛星精確的軌道位置。采樣率為15分鐘,實際解算中可以進行精密鐘差的估計或內插,以提高其可使用的歷元數。
IGS多系統精密星歷下載地址:ftp://cddis.gsfc.nasa.gov/pub/gps/products/mgex/
文件根據GPS周分為很多文件夾,計算GPS周可以通過以下小程序:
https://download.csdn.net/download/gou_hailong/11739153
?
- 命名規則
用的sp3格式的命名規則為:tttwwwwd.sp3
其中:ttt表示精密星歷的管理機構名稱,包括{wum(武漢大學)、tum(慕尼黑大學)}wwww表示GPS周;d表示星期,0表示星期日,1~6表示星期一至星期六。文件名如:igs20121.sp3,其中igs為計算單位名,2012為GPS周,1為星期一。
- 文件格式
文件可以分為文件頭與文件正文兩部分,其中文件頭格式如下:
第一行:文件版本 年 月 日 時 分 秒? 歷元數量 數據類型? 協調系統 軌道類型? 管理機構
第二行:GPS周 第一個歷元周內秒 歷元間隔 約化儒略歷整數 約化儒略歷小數
第三至十二行:衛星數量 衛星序號
第十三至二十二行:衛星序號對應的衛星精度
第二十三至二十四行:未使用
。。。。。。
文件正文格式如下:
‘P’表示‘衛星位置’,‘V’表示‘衛星速度’,‘G’表示‘GPS衛星’,‘R’表示'GLONASS衛星',‘E’表示‘伽利略衛星’,‘C’表示‘北斗衛星’
第一列:位置速度選項?系統類型 衛星序號
第二列:x坐標
第三列:y坐標
第四列:z坐標
第五列:鐘差
文件讀取程序如下:
%讀取精密星歷文件夾中的數據 global SP3_data folders = dir('C:\Users\zhupf\Desktop\IGS數據\SP3\DATA\mix\解壓后'); folders = {folders.name}; folders = setdiff(folders,{'.','..'})'; num = length(folders); t=1; for j = 1:numtline = folders{j};fid=fopen(['C:\Users\zhupf\Desktop\IGS數據\SP3\DATA\mix\解壓后\' tline ],'rt'); Head_data=struct; SP3_data =struct; if(fid==-1)msgbox('輸入的文件或者路徑不正確,無法正確打開sp3文件','警告信息');return; end % 讀取文件頭,并將其存入Head_data結構數組中 Head=[]; for i=1:31 %制度文件頭的前16行line=fgetl(fid);Head=[Head;line];%將sp3文件的文件頭以字符的形式顯示出來,切記此處要保證sp3文件中的頭文件的每行字符數是相同的,否則讀取將失敗 end line=fgetl(fid); %提取文件頭第一行的數據 Head_data.Version=Head(1,1:2); Head_data.P_or_V=Head(1,3); Head_data.start_year=str2num(Head(1,4:7)); Head_data.start_month=str2num(Head(1,9:10)); Head_data.start_day=str2num(Head(1,12:13)); Head_data.start_hour=str2num(Head(1,15:16)); Head_data.start_minute=str2num(Head(1,18:19)); Head_data.start_second=str2num(Head(1,21:31)); Head_data.EpochNum=str2num(Head(1,33:39)); Head_data.data_type=Head(1,41:45); Head_data.Coordinate_system=Head(1,47:51); Head_data.Orbit_type=Head(1,53:55); Head_data.Institution=Head(1,57:60); %第二行 Head_data.GPS_week=str2num(Head(2,4:7)); Head_data.FirstEpoch_weeksecond=str2num(Head(2,9:23)); Head_data.Epoch_interval=str2num(Head(2,25:38)); Head_data.MJD_zhengshu=str2num(Head(2,40:44)); Head_data.MJD_xiaoshu=str2num(Head(2,46:60)); %第三 至 七行-------提取觀測衛星的PRN編號 Head_data.SatNum=str2num(Head(3,5:6)); PRN=strcat(Head(3,10:60),Head(4,10:60),Head(5,10:60),Head(6,10:60),Head(7,10:60)); L=length(PRN)/3; SP3_PRN=[]; for i=1:LSp3_PRN=str2num(PRN((3*i-1):3*i));if Sp3_PRN==0continue;endSP3_PRN=[SP3_PRN Sp3_PRN]; end Head_data.SP3_PRN=SP3_PRN; %第十三 至 十七行-----提取相應衛星的精度 SatAccu=strcat(Head(13,10:60),Head(14,10:60),Head(15,10:60),Head(16,10:60),Head(17,10:60)); L=length(SP3_PRN); SP3_SatAccu=[]; for i=1:LSp3_SatAccu=str2num(SatAccu((3*i-1):3*i));if Sp3_SatAccu==0continue;endSP3_SatAccu=[SP3_SatAccu Sp3_SatAccu]; end Head_data.SP3_SatAccu=SP3_SatAccu; %第二十五,二十六行-----讀取標準差的浮點基數 Head_data.pos_jishu=str2num(Head(25,4:13)); Head_data.Satclock_jishu=str2num(Head(25,15:26)); %--------------------End of the reading of the Fileheader------------------%--------------------讀取衛星坐標及衛星鐘鐘差-------------------------------- %讀取的衛星的三維坐標的單位為km,衛星鐘差的單位為uswhile 1 line=fgetl(fid);result=findstr(line,'EOF');if(isempty(result))SP3_data(t).year=str2num(line(4:7));SP3_data(t).month=str2num(line(9:10));SP3_data(t).day=str2num(line(12:13));SP3_data(t).hour=str2num(line(15:16));SP3_data(t).minute=str2num(line(18:19));SP3_data(t).second=str2num(line(21:22));SP3_data(t).mm=SP3_data(t).hour*60+SP3_data(t).minute+ SP3_data(t).second/60;%將該歷元的時分秒化為分鐘,方便后面和觀測歷元比較for j=1:Head_data.SatNumline=fgetl(fid);SP3_data(t).PRN(j)=str2num(line(3:4));%將每個歷元的衛星編號保存在一個數組中SP3_data(t).x(j)=str2num(line(5:18)); %提取衛星的三維坐標,單位為kmSP3_data(t).y(j)=str2num(line(19:32));SP3_data(t).z(j)=str2num(line(33:46));SP3_data(t).Clock_corr(j)=str2num(line(47:60));%提取衛星鐘差,單位為us%判斷是否有壞的或空缺的位置值if(SP3_data(t).x(j)==0||SP3_data(t).y(j)==0||SP3_data(t).z(j)==0)fprintf('注意:星歷文件中的第%g個歷元的%g號衛星的位置值是壞的或空缺的!\n',t,SP3_data(t).PRN(j))end%判斷是否為壞的或空缺的鐘差值if(SP3_data(t).Clock_corr(j)==999999.999999)fprintf('注意:星歷文件中的第%g個歷元的%g號衛星的鐘差值是壞的或空缺的!\n',t,SP3_data(t).PRN(j))end%判斷X的標準差是否為未知,若是將其置為0if(length(line)>=73 && isempty(str2num(line(62:63))))SP3_data(t).X_mi(j)=0;elseif(length(line)>=73 && ~isempty(str2num(line(62:63))))SP3_data(t).X_mi(j)=str2num(line(62:63));end%判斷Y的標準差是否為未知,若是將其置為0if(length(line)>=73 && isempty(str2num(line(65:66))))SP3_data(t).Y_mi(j)=0;elseif(length(line)>=73 && ~isempty(str2num(line(65:66))))SP3_data(t).Y_mi(j)=str2num(line(65:66));end%判斷Z的標準差是否為未知,若是將其置為0if(length(line)>=73 && isempty(str2num(line(68:69))));SP3_data(t).Z_mi(j)=0;elseif(length(line)>=73 && ~isempty(str2num(line(68:69))))SP3_data(t).Z_mi(j)=str2num(line(68:69));end%判斷衛星鐘差的標準差是否為未知,若是將其置為0if(length(line)>=73 && isempty(str2num(line(71:73))))SP3_data(t).Clock_mi(j)=0;elseif(length(line)>=73 && ~isempty(str2num(line(71:73))))SP3_data(t).Clock_mi(j)=str2num(line(71:73));end%查看某個衛星坐標和星站距離變化if PRN(3*j-2:3*j) == 'G01' %選擇需要查看的衛星序號x(t) = SP3_data(t).x(j);y(t) = SP3_data(t).y(j);z(t) = SP3_data(t).z(j);if t== 96*7QD = [-2665073.9893,4556123.9737,3776821.7794]; %測站坐標d = sqrt((x*1000-QD(1)).^2+(y*1000-QD(2)).^2+(z*1000-QD(3)).^2);figure(1);plot(x,'r','LineWidth',2.5);ylabel('x坐標');xlabel('時間');title('GPS 1號衛星x坐標');set(gca,'Xtick',[0 96 192 288 384 480 576]);set(gca,'Xticklabel',{'星期日','星期一','星期二','星期三','星期四','星期五','星期六'});figure(2);plot(y,'r','LineWidth',2.5);ylabel('y坐標');xlabel('時間');title('GPS 1號衛星y坐標');set(gca,'Xtick',[0 96 192 288 384 480 576]);set(gca,'Xticklabel',{'星期日','星期一','星期二','星期三','星期四','星期五','星期六'});figure(3);plot(z,'r','LineWidth',2.5);ylabel('z坐標');xlabel('時間');title('GPS 1號衛星z坐標');set(gca,'Xtick',[0 96 192 288 384 480 576]);set(gca,'Xticklabel',{'星期日','星期一','星期二','星期三','星期四','星期五','星期六'});figure(4);plot(d,'r','LineWidth',2.5);ylabel('星站距離');xlabel('時間');title('GPS 1號衛星與測站距離'); % legend('星站距離');set(gca,'Xtick',[0 96 192 288 384 480 576]);set(gca,'Xticklabel',{'星期日','星期一','星期二','星期三','星期四','星期五','星期六'});endendendelseif(~isempty(result))break;%判斷讀到文件的末尾,到達文件尾則退出endt=t+1;%記錄文件體的行數 end fclose(fid); end程序運行結果:
?
?
總結
以上是生活随笔為你收集整理的多系统精密星历下载与分析的全部內容,希望文章能夠幫你解決所遇到的問題。
- 上一篇: 前端学习(2110):组件化得开发和实现
- 下一篇: 前端学习(2155):htmlwebpa