本文介绍了如何使用SNR数据进行水位高度反演，涉及数据预处理、计算高度角、方位角和信噪比、PRN号选择以及反演过程。在数据预处理中，通过gfzrnx软件进行数据合并和筛选，并转换为rinex2格式。计算高度角和信噪比时，借助TEQC和RTKLIB解决相关问题。在反演阶段，针对有无潮位对比数据的情况，进行了详细的操作步骤和处理方法。
目录
1 数据预处理
1.1 数据合并与筛选
使用gfzrnx软件合并观测数据
gfzrnx -finp ZJGH23**.20o ZJGH24**.20o -fout ::RX3::
使用gfzrnx软件合并导航电文数据
gfzrnx -finp ZJGH28**.20p -fout ::RX3::
使用gfzrnx软件筛选导航电文数据卫星系统
gfzrnx -finp ZJGH00XXX_R_20202500000_10D_MN.rnx -satsys C -fout ZJGH00XXX_R_20202500000_10D_MN.bds
使用gfzrnx软件筛选观测数据卫星系统
gfzrnx -finp ZJGH00XXX_R_20202370527_06D_15S_MO.rnx -satsys G -fout ZJGH00XXX_R_20202370527_06D_15S_MO.gps
1.2 数据格式转换
使用gfzrnx软件转换为rinex2格式
gfzrnx -finp ZJGH00XXX_R_20202370527_06D_15S_MO.gps -vo 2 -fout ::RX2::
2 计算高度角、方位角、信噪比
2.1 TEQC计算高度角、方位角、信噪比
teqc +qc +plot -nav zjgh2370.20c zjgh2370.20o
2.2 GPS导航电文无法计算高度角问题
GPS导航电文出现错误，TEQC无法计算其高度角和方位角
解决方法：利用RTKLIB_bin 2.4.2 p11将源数据进行Rinex 2.11格式转换
转换之后TEQC可正常计算高度角、方位角、信噪比
2.3 从质量检核文件中提取高度角、方位角、信噪比
%在 matlab下新建一个m文件，将以下代码直接拷贝进去，即可执行。
%需要一个TEQC生成的plot文件作为参数
function out=teqcplot3(files);
%读取TEQC生成的Plot文件，绘制数据图表，支持Copmact、Compact2、Compact3格式
%选取一个TEQC的Plot文件
% 格式说明
% *.sn1 载波 L1的信噪比 Signal to noise ratio (S/N)
% *.sn2 载波 L2的信噪比 Signal to noise ratio (S/N) Carrier L2
% *.iod *.d12 *.d21 电离层延迟观测值变化率(米/秒) Derivative of ionospheric delay observable (m/s)
% *.ion *.i12 *.i21 电离层延迟观测值(米)Ionospheric delay observable (m)
% *.mp1 *.m12 载波L1的多路径误差 Multipath Carrier L1
% *.mp2 *.m21 载波L2的多路径误差 Multipath Carrier L2
% *.azi 卫星方位角（°） Satellite azimuthal data (degrees)
% *.ele 卫星高度角（°） Satellite elevation data (degrees)
if nargin==0
[filen,path]=uigetfile('*.sn1;*.sn2;*sn7;*.iod;*.ion;*.mp1;*.mp2;*.azi;*.ele;*.i12;*i21;*.m12;m21;*.d12;*d21',...
'请选择TEQC报告文件：');
else
[path,filen,ext]=fileparts(files);
path=[path '\'];
filen={
[filen ext]};
end
%读取这个文件
file=char(filen);
%按行读取文件至数组A
[A]=importdata([path file],'\t');
%定义SAT，存放卫星数据
%GPS有32颗卫星，存放序号1-32，
%GLONASS有32颗卫星，存放序号33-64，
%BEIDOU有35颗卫星，存放序号65-99
SAT(1:length(A),1:99)=NaN;
%sats(1:length(A),1:99)=NaN;
%存放采样时间，单位秒
tsec(1:length(A))=NaN;
%读取文件的第一行
filelx=A{
1};
%判断是哪种格式
switch filelx
case 'COMPACT'
%读取数据采样间隔
t_samp=char(A(3));
%读取开始时间
mjl=char(A(4));
%读取数据采样间隔
T_SAMP=str2num(t_samp(max(find(t_samp==' ')):end));
%读取数据采样开始时间
MJL_START=str2num(mjl(max(find(mjl==' ')):end));
%转成时间序列数字，date serial number，从0000年1月1日0时0分0秒开始计算的十进制天数
MJD_START=MJL_START+678941.999999741;
%i为行号
n=1;i=5;
case 'COMPACT2'
%读取数据采用间隔
t_samp=char(A(2));
%读取开始时间
mjl=char(A(3));
%读取数据采样间隔
T_SAMP=str2num(t_samp(max(find(t_samp==' ')):end));
%读取数据采样开始时间
MJL_START=str2num(mjl(max(find(mjl==' ')):end));
%转成时间序列数字，date serial number，从0000年1月1日0时0分0秒开始计算的十进制天数
MJD_START=MJL_START+678941.999999741;
n=1;i=4;
case 'COMPACT3'
%读取开始时间
t_start=char(A(2));
t_start = deblank(t_start);
s = splitstr(t_start,'**',6);
%t_start_time=[char(s{
2}) '年' char(s{
3}) '月' char(s{
4}) '日' char(s{
5}) '时' char(s{
6}) '分' num2str(str2num(char(s{
7})),'%02d') '秒'] ;
%获取采样的开始时间，2013,12，7,03,05,55
t_s_time=[str2num(char(s{
2})), str2num(char(s{
3})),str2num(char(s{
4})),str2num(char(s{
5})),str2num(char(s{
6})),str2num(char(s{
7}))];
n=1;i=3;
otherwise
disp('数据格式存在问题');
return
end
%n=1;i=3;
%sats=str2num(A{
i});
snyggfilen=strrep(filen,'_','-');
%生成进度条
h=waitbar(0,['正在读取数据，请稍等…… ' char(snyggfilen)]);
while i<length(A);
waitbar(i/length(A),h);drawnow
%读取卫星编号数据行 采样时间 卫星数量 卫星编号 COMAPCT3
satbhstr=char(A(i));
s = splitstr(satbhstr,'**',32);
%如果卫星编号数据为0，说明数据缺失，退出本函数
if strtrim(satbhstr)=='0'
disp ('数据为空！');
close(h);
return
end
switch filelx
case 'COMPACT3'
%记录采样时间，单位秒
tsec(n)=str2num(s{
1});
if s{
2}=='-1'
%如果为-1，使用上次的卫星编号行
satbhstr=oldsatbh;
s = splitstr(satbhstr,'**',32);
else
oldsatbh=satbhstr;
end
%获取当前卫星数量
satcount=str2num(s{
2});
case {
'COMPACT','COMPACT2'}
%记录采样时间，单位秒
tsec(n)=(n-1)*T_SAMP;
if s{
1}=='-1'
%如果为-1，使用上次的卫星编号行
satbhstr=oldsatbh;
s = splitstr(satbhstr,'**',32);
else
oldsatbh=satbhstr;
end
%获取当前卫星数量
satcount=str2num(s{
1});
end
%读取卫星数据行
satdatastr=char(A(i+1));
sdata=splitstr(satdatastr,'**',32);
for k=3:satcount+2
switch filelx
case 'COMPACT'
%如果是COMPACT格式，卫星编号前面没有字母，默认是GPS卫星，在前面增加字符G
satbh=['G' char(s{
k-1})];
case 'COMPACT2'
%如果是COMPACT2格式，卫星编号行第2个数据开始是卫星编号
satbh=char(s{
k-1});
case 'COMPACT3'
%如果是COMPACT3格式，卫星编号行第3个数据开始是卫星编号
satbh=char(s{
k});
end
switch satbh(1);
case 'G'
%如果是GPS卫星，获得与编号对应的存储序号，范围1-32
index=str2num(satbh(2:3));
case 'R'
%如果是GLONASS卫星，获得与编号对应的存储序号，范围33-64
index=32+str2num(satbh(2:3));
case 'C'
%如果是BEIDOU卫星，获得与编号对应的存储序号，范围65-99
index=64+str2num(satbh(2:3));
otherwise
%其他情况，获得与编号对应的存储序号，范围1-32
index=str2num(satbh(2:3));
end %switch
%获得对应的卫星数据
sdatastr=char(sdata{
k-2});
%如果卫星数据不是数值型，用0替代
if isempty(str2num(sdatastr))==1
sdatastr='0.000';
end
%将卫星数据存入SAT数组的对应位置
SAT(n,index)=str2num(sdatastr);
end %for
n=n+1;
%根据数据文件结构，一行是卫星编号，下一行就是对应的卫星数据
%程序一次读取2行数据
i=i+2;
end %while
%数据读取完毕，关闭进度条
close(h);
%获取已使用的卫星编号，返回给sat_bh
sat_bh=getsatbh(SAT);
%增加最后一个结束字符
sat_bh_temp=[sat_bh' {'End'}];
sat_bh_end=sat_bh_temp';
%删去所有都为NaN值的列，返回给sat_data
sat_data=delnancol(SAT);
%删去所有NaN值的采样时间列
t_seconds=delnancol(tsec);
%删去所有都为NaN值的行，返回sat_datas
sat_datas=delnanrow(sat_data);
%增加一个nan列
[row,col]=size(sat_datas);
bb(1:row,1)=NaN;
sat_datas=[sat_datas bb];
%将含有NaN值的数据替换为0，返回给s_data
%s_data=repnan20(sat_data);
%将采样时间转成顺序日期serial date number
switch filelx
case {
'COMPACT','COMPACT2'}
t_s_jd_time=MJD_START;
case 'COMPACT3'
t_s_jd_time=datenum(t_s_time);
end
for i=1:length(t_seconds)
%将每次采样时间转换成相应的顺序日期
t_jd_time(i)=t_s_jd_time+t_seconds(i)/60/60/24;
end
%绘制图表++++++++++++++++++++++++++
%根据不同的文件类型，设置坐标轴的范围
[type,maxy,miny]=get_filetype(file);
figure;box on;hold on
[row,col]=size(sat_datas);
%绘制渐变彩色图
pcolor(t_jd_time',[1:col],sat_datas');
set(gcf,'renderer','zbuffer');
shading flat
set(gca,'xticklabel',[t_jd_time]);
set(gca,'xlim',[t_jd_time(1) t_jd_time(end)])
%t_start_times=[t_start_h ':' t_start_m ':' t_start_s];
dateaxis('X',13)
cbar('v',[miny maxy],type);
set(gca,'ylim',[1 col+1])
set(gca,'ytick',[1.5:1:col+0.5])
set(gca,'yticklabel',sat_bh_end);
set(gca,'fontsize',7);
colormap(flipud(jet));
caxis([miny maxy]);
%t_end_time=cal2et(t_s_time,t_seconds(end));
xlabel([datestr(t_jd_time(1)) ' |-------- 采样间隔: ' num2str(t_seconds(2)-t_seconds(1)) ' 秒 --------| ' datestr(t_jd_time(end))])
ylabel('卫星编号')
%timestr=secs2hms(length(sat));
T=title(['TEQC 报告文件: ' strrep(file,'_','-')]);
set(T,'fontsize',8)
%out.(file(end-2:end))=sat;
%out.T_samp=T_SAMP;
out.Start=datestr(t_jd_time(1));
out.Stop=datestr(t_jd_time(end));
%+++++++++++++++++++
最低0.47元/天 解锁文章
扫一扫
2万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
