PPP
1.读取文件
1.1RINEX观测文件(.**o)
观测文件这部分在SPP里面已经进行详细的介绍
1.2精密星历文件(.sp3)
调用readsp3()函数
function nav=readsp3(nav,fname)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% only surpport to read ".sp3" or ".eph" file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global glc
%首先是查找文件名,是否是空的,如果没有发现文件名的格式为.sp3、eph等,提示错误,文件格式不支持
if ~isempty(strfind(fname,'.'))
if isempty(strfind(fname,'sp3'))&&isempty(strfind(fname,'SP3'))&&...
isempty(strfind(fname,'eph'))&&isempty(strfind(fname,'EPH'))
error('The file %s is not surpported!!!',fname);
end
end
idx=strfind(fname,glc.sep); fname0=fname(idx(end)+1:end);%这段代码的作用是从字符串 fname 中截取最后一个分隔符 glc.sep 后面的部分,并将结果保存到变量 fname0 中。这可能用于文件名的处理,提取文件名中的文件扩展名或最后一个目录的名称。
fprintf('Info:reading sp3 file %s...',fname0);%读取后缀为.sp3文件的名字,最终在命令行窗口输出Info:reading obs file wum20464.sp3...over
fid=fopen(fname);%打开sp3文件
% read sp3 header
[headinfo,fid]=decode_sp3h(fid);%读取精密星历的头文件
% read sp3 body
nav=decode_sp3b(nav,headinfo,fid);%读取精密星历的主体文件
fprintf('over\n');
return
在读取精密星历文件的时候分为两部分,先通过调用decode_sp3h()函数来读取精密星历的头文件,之后通过调用decode_sp3b()函数来读取精密星历的主体文件。
调用decode_sp3h()函数
function [headinfo,fid]=decode_sp3h(fid)%一个函数 decode_sp3h(fid),用于解析 SP3 文件的头部信息。
%函数的输入参数是文件句柄 fid,表示要解析的 SP3 文件。函数通过逐行读取文件,提取文件头部信息,并将其保存在一个 headinfo 结构体中。
%函数首先定义了一些变量 type、time、ns、sats、tsys 和 bfact,用来存储头部信息的各个字段。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%read file header information
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global glc gls
type=' '; time=gls.gtime; ns=0;
sats=zeros(1,glc.MAXSAT); tsys=''; bfact=zeros(1,2);
n=0; idxns=0; idxc=0; idxf=0; idx_break=0;
while ~feof(fid)%然后,函数通过循环读取文件的每一行,并根据行的不同内容,执行不同的操作来解析头部信息。
line=fgetl(fid);n=n+1;
if n==1%当 n 等于 1 时,解析第一行,获取文件类型 type 和时间信息 time。
type=line(3);
time=str2time(line(4:31));
elseif strcmp(line(1),'+')%当行以 + 开头时,表示卫星编号信息,解析卫星编号,并存储在 sats 数组中。
if n==3,ns=str2double(line(3:6));end
if idxns<ns
j=0;
while j<17
sys=code2sys(line(10+3*j));%在每次迭代中,首先使用函数 code2sys(line(10+3*j)) 将字符串中的字符转换为卫星系统代码(如 GPS、GLO、GAL 等),存储在变量 sys 中。
prn=str2double(line(11+3*j:12+3*j));%使用 str2double(line(11+3*j:12+3*j)) 解析出卫星的 PRN 号码,存储在变量 prn 中。
if sys==glc.SYS_QZS,prn=prn+192;end%如果卫星系统为 QZSS (glc.SYS_QZS),则将 PRN 号码加上 192,因为 QZSS 的 PRN 号码是 GPS PRN 号码加上 192。
if idxns<=glc.MAXSAT
sats(idxns+1)=satno(sys,prn);
idxns=idxns+1;
end%最后,如果索引 idxns 小于等于 glc.MAXSAT,则将卫星的编号存储在数组 sats 的对应位置上,并将索引 idxns 加 1。
j=j+1;
end
end
elseif strcmp(line(1:2),'%c')&&idxc==0%当行以 %c 开头时,表示时钟系统信息,解析并存储在 tsys 变量中。
tsys=line(10:12);
idxc=idxc+1;
elseif strcmp(line(1:2),'%f')&&idxf==0 %当行以 %f 开头时,表示钟差比例因子信息,解析并存储在 bfact 数组中。
bfact(1)=str2double(line(4:13));
bfact(2)=str2double(line(15:26));
idxf=idxf+1;
elseif strcmp(line(1:2),'/*')%当行以 /* 开头时,表示头部信息解析结束,跳出循环。
idx_break=idx_break+1;
if idx_break==4
break;
end
end
end
%最后,将解析得到的头部信息存储在 headinfo 结构体中,并返回该结构体作为函数的输出。
headinfo.type=type;
headinfo.time=time;
headinfo.ns=ns;
headinfo.sats=sats;
headinfo.tsys=tsys;
headinfo.bfact=bfact;
return
调用decode_sp3b()函数
function nav=decode_sp3b(nav,headinfo,fid)
global glc gls
peph=gls.peph;%精密星历
NMAX=10000; nav.peph=repmat(gls.peph,NMAX,1);
if headinfo.type=='P'
ns=headinfo.ns;
else
ns=2*headinfo.ns;
end
%这段代码的作用是初始化导航结构体 nav 中的精密星历部分,并根据头部信息结构体 headinfo 中的数据确定导航结构体中精密星历的数量 ns。
%按照时间步进遍历)是一种处理数据或执行操作的方法。它通常用于按照时间顺序逐个处理一系列数据点或执行一系列操作。在时间序列数据处理中,"traversal by epoch"指的是按照时间的顺序,从最早的数据点开始,逐个处理每个数据点,并按照一定的步进(如固定时间间隔)依次处理下一个数据点,直到遍历完所有的数据点。
while ~feof(fid)
line=fgetl(fid);
if strcmp(line,'EOF'),continue;end
if line(1)~='*'
fprintf('sp3 invalid epoch\n');
continue;
end
time=str2time(line(4:31));
if strcmp(headinfo.tsys,'UTC'),time=utc2gpst(time);end
peph.time=time;
%这段代码的作用是解析 SP3 文件中的每个时刻的数据,并将时间信息存储在精密星历结构体 peph 中的 time 字段中。该代码还进行了一些错误处理,例如检查是否到达文件末尾、验证时间信息的有效性等。
%初始化精密星历的结构体
for i=1:glc.MAXSAT
for j=1:4
peph.pos(i,j)=0;%位置
peph.std(i,j)=0;%标准差
peph.vel(i,j)=0;%速度
peph.vst(i,j)=0;%速度标准差
end
for j=1:3
peph.cov(i,j)=0;%协方差
peph.vco(i,j)=0;%速度协方差
end
end
%这段代码的作用是将精密星历结构体 peph 中每个卫星的位置、标准差、速度、速度标准差、协方差矩阵和速度协方差矩阵等字段初始化为 0,为后续的数据填充和存储做好准备。
for i=1:ns
line=fgetl(fid);
if size(line,2)<4||(line(1)~='P'&&line(1)~='V'),continue;end
if line(2)==' '
sys=glc.SYS_GPS;
else
sys=code2sys(line(2));
end
prn=str2double(line(3:4));
if sys==glc.SYS_QZS,prn=prn+192;end
sat=satno(sys,prn);
if ~sat,continue;end
%此代码段的作用是解析 SP3 文件中每个卫星的信息,并通过解析卫星系统类型和 PRN 号将其转换为卫星号。同时,此代码段还进行了一些错误处理,例如检查行的长度是否符合要求、验证卫星系统类型和 PRN 号的有效性等。
for j=1:4%对于前三个 j 值,执行以下赋值操作:m=2:将变量 m 的值设置为 2。k1=1000:将变量 k1 的值设置为 1000。%base=headinfo.bfact(1):将变量 base 的值设置为 headinfo.bfact(1) 的值,headinfo.bfact 是一个结构体或数组,表示其中的第一个元素。%base=headinfo.bfact(1):将变量 base 的值设置为 headinfo.bfact(1) 的值,headinfo.bfact 是一个结构体或数组,表示其中的第一个元素。%k3=0.1:将变量 k3 的值设置为 0.1。k4=10^-7:将变量 k4 的值设置为 0.0000001。
if j<=3
m=2;k1=1000;base=headinfo.bfact(1);k2=10^-3;k3=0.1;k4=10^-7;
else
m=3;k1=10^-6;base=headinfo.bfact(2);k2=10^-12;k3=10^-10;k4=10^-16;
end%m=3:将变量 m 的值设置为 3。k1=10^-6:将变量 k1 的值设置为 0.000001。base=headinfo.bfact(2):将变量 base 的值设置为 headinfo.bfact(2) 的值,headinfo.bfact 是一个结构体或数组,表示其中的第二个元素。k2=10^-12:将变量 k2 的值设置为 0.000000000001。k3=10^-10:将变量 k3 的值设置为 0.0000000001。
val=str2double(line(4+14*(j-1)+1:4+14*(j-1)+14));
if size(line,2)>61
std=str2double(line(61+3*(j-1)+1:61+3*(j-1)+m));
else
std=0;
end%段代码的作用是解析文本行 line 中特定位置的字符串,并将其转换为数字,分别存储在变量 val 和 std 中。具体的位置计算是根据变量 j 和之前的条件判断得出的。
if line(1)=='P'
if val~=0&&abs(val-999999.999999)>=10^-6
peph.pos(sat,j)=val*k1; v=1;%这两行代码的作用是对变量 peph.pos 和 v 进行赋值操作。根据给定的索引 sat 和 j,将 val 与 k1 相乘的结果赋值给 peph.pos。同时,将变量 v 的值设置为 1。
end
if base>0&&std>0
peph.std(sat,j)=base^std*k2;
end%这段代码的作用是根据不同的条件,对 peph 数据结构进行赋值操作。如果 line 的第一个字符为 'P',并且满足一定条件,将 val 乘以特定的系数并赋值给 peph.pos,将 base 的 std 次方乘以另一个系数并赋值给 peph.std。同时,在满足特定条件时,设置变量 v 的值为 1。
elseif v
if val~=0&&abs(val-999999.999999)>=10^-6
peph.vel(sat,j)=val*k3;
end
if base>0&&std>0
peph.vst(sat,j)=base^std*k4;
end%这段代码是在之前的代码的if line(1)=='P'判断条件不成立的情况下执行的,它执行了以下操作:首先,通过判断变量 v 是否为真(非零)来执行接下来的代码块。
%判断 val 不等于 0,并且 val 的绝对值与 999999.999999 的差值大于等于 10 的负 6 次方(即绝对值大于等于 0.000001)。如果条件满足,执行以下操作:
%将变量 val 乘以 k3,并将结果赋值给 peph.vel(sat, j)。peph.vel 是一个数组或矩阵,并且 sat 和 j 用作索引。
%如果 base 大于 0 且 std 大于 0,则执行以下操作:
%将 base 的 std 次方乘以 k4,并将结果赋值给 peph.vst(sat, j)。peph.vst 是一个数组或矩阵,并且 sat 和 j 用作索引。
%简而言之,这段代码的作用是在 v 为真的情况下,根据一些条件对 peph 数据结构的 vel 和 vst 进行赋值操作。如果条件满足,将 val 乘以特定的系数,并赋值给 peph.vel 或 peph.vst。
end
end
end
if v
if nav.np>size(nav.peph,1)
nav.peph(nav.np+1:nav.np+NMAX,1)=repmat(gls.peph,NMAX,1);
end%nav.np:这是一个变量,表示当前位置;nav.peph:这是一个数组或矩阵,存储了位置信息。size(nav.peph,1):这是一个函数,返回 nav.peph 的行数。nav.np < size(nav.peph,1):这是一个条件判断,检查 nav.np 是否小于 nav.peph 的行数。
nav.peph(nav.np+1)=peph;
nav.np=nav.np+1;
end%这段代码的作用是在 v 为真的情况下,根据一些条件将变量 peph 添加到 nav.peph 中的特定位置。如果条件满足,则根据需要扩展 nav.peph 的大小,并将 peph 添加到 nav.peph 中的下一个空闲位置,同时更新 nav.np 的值。
end
fclose(fid);
if nav.np<size(nav.peph,1)
nav.peph(nav.np+1:end,:)=[];
end%这段代码的目的是将 nav.peph 中超出当前索引范围的部分数据移除。如果 nav.np 较小,导致 nav.peph 中存在多余的行数,则通过删除超出范围的行来保持数据的正确性。
return
1.3精密钟差文件(.clk)
调用readclk()函数
function nav=readclk(nav,opt,fname)
global glc gls
NMAX=10000; nav.pclk=repmat(gls.pclk,NMAX,1);
idx=strfind(fname,glc.sep); fname0=fname(idx(end)+1:end);%这段代码主要是对变量进行赋值操作。它将 gls.pclk 在行方向上复制了 NMAX 次,并将结果赋值给 nav.pclk。然后,通过 strfind 函数找到 fname 中最后一个斜杠的位置,并将斜杠后面的部分赋值给 fname0。
fprintf('Info:reading clk file %s...',fname0);%读取后缀为.clk文件的名字,最终在命令行窗口输出Info:reading obs file wum20464.clk...over
fid=fopen(fname);%代开精密钟差文件
% 读取精密钟差的头文件
while ~feof(fid)
line=fgetl(fid);
label=line(61:end);
if size(line,2)<=60
continue;
elseif ~isempty(strfind(label,'RINEX VERSION / TYPE'))
ver=str2double(line(1:9)); %#示例数据版本为3.0,已经经过对比源文件进行查看
type=line(21);%根据示例数据类型为'C'
elseif ~isempty(strfind(label,'END OF HEADER'))
break;
end
end
if type~='C'%type为C的话,那么表示下载的是正确的clk文件的数据,如果type不为C的话,提示错误“该文件不是clock file”
error('The file %s is not clock file!!!',fname0);
end
mask=set_sysmask(opt.navsys);%这个函数可以用于根据输入参数的不同设置系统屏蔽状态,进而控制系统的行为。
% 读取精密钟差的主体文件
while ~feof(fid)
buff(1:glc.MAXRNXLEN)=' ';
line=fgets(fid);
buff(1:size(line,2))=line;
if size(buff,2)<8,continue; end
if ~strcmp(buff(1:2),'AS'),continue;end
%这段代码的作用是按行读取文件内容,并对读取到的每行进行一些条件判断。它首先将 buff 数组初始化为空格字符,然后使用 fgets 函数逐行读取文件内容,并将其复制到 buff 数组中。接着,它进行一些条件判断,例如判断行的长度是否足够、是否以特定的字符开头等,从而决定是否继续处理该行。
satid=buff(4:8);
sat=satid2no(satid);
[sys,~]=satsys(sat);
if sys==0,continue;end
if mask(sys)==0,continue;end
%这段代码提取了行内容中的卫星标识符,并根据卫星标识符获取卫星编号和导航系统类型。然后,根据导航系统类型进行一些条件判断,例如判断导航系统是否被禁用,从而决定是否继续处理该行。
time=str2time(buff(8:34));%time=str2time(buff(8:34));:将 buff 数组中第 8 到第 34 个元素(即行内容的第 8 到第 34 个字符)作为字符串参数,调用函数 str2time 进行解析和转换,并将转换后的时间结果赋值给变量 time。该行代码的目的是将字符子串解析为时间值。
clk=str2double(buff(41:59));%clk=str2double(buff(41:59));:将 buff 数组中第 41 到第 59 个元素(即行内容的第 41 到第 59 个字符)作为字符串参数,调用函数 str2double 进行解析和转换,并将转换后的数值结果赋值给变量 clk。该行代码的目的是将字符子串解析为双精度浮点数。
std=str2double(buff(61:79));%std=str2double(buff(61:79));:将 buff 数组中第 61 到第 79 个元素(即行内容的第 61 到第 79 个字符)作为字符串参数,调用函数 str2double 进行解析和转换,并将转换后的数值结果赋值给变量 std。该行代码的目的是将字符子串解析为双精度浮点数。
if isnan(clk),clk=0;end
if isnan(std),std=0;end
%这段代码的作用是在变量 clk 和 std 的值为 NaN 时,将它们置为 0,以避免计算或使用这些值时出现错误。NaN 是一种特殊的浮点数值,表示不是一个数字。
if nav.nc<=0%这段代码是用于更新卫星钟差信息的逻辑判断和操作。下面是每行代码的含义:
nav.nc=nav.nc+1;%if nav.nc<=0:如果变量 nav.nc 小于或等于 0,即当前没有卫星钟差信息,执行下面的操作。
%nav.nc=nav.nc+1;:将变量 nav.nc 增加 1,表示有一个卫星钟差信息。
nav.pclk(nav.nc).time=time;%nav.pclk(nav.nc).time=time;:将当前时间 time 赋值给 nav.pclk 数组中的 nav.nc 索引的 time 字段。
elseif abs(timediff(time,nav.pclk(nav.nc).time))>1e-9
if nav.nc+1>size(nav.pclk,1)
nav.pclk(nav.nc+1:nav.nc+NMAX)=repmat(gls.pclk,NMAX,1);
end
%%elseif abs(timediff(time,nav.pclk(nav.nc).time))>1e-9:如果当前时间 time 与 nav.pclk 数组中的 nav.nc 索引的 time 字段之间的时间差值大于 1e-9 秒(即差值超过 1 纳秒),执行下面的操作。
%if nav.nc+1>size(nav.pclk,1):如果 nav.nc+1 大于 nav.pclk 数组的行数(即超过数组的大小),执行下面的操作。
%nav.pclk(nav.nc+1:nav.nc+NMAX)=repmat(gls.pclk,NMAX,1);:将 gls.pclk 数组重复 NMAX 次,并赋值给 nav.pclk 数组中从 nav.nc+1 到 nav.nc+NMAX 索引的位置,即扩展 nav.pclk 数组的大小以容纳更多的卫星钟差信息。
%nav.nc=nav.nc+1;:将 nav.nc 增加 1,表示有一个新的卫星钟差信息。
%nav.pclk(nav.nc).time=time;:将当前时间 time 赋值给 nav.pclk 数组中的 nav.nc 索引的 time 字段。
nav.nc=nav.nc+1;
nav.pclk(nav.nc).time=time;%该代码段的作用是更新卫星钟差信息,在 nav.pclk 数组中加入一个新的时间值,并将 nav.nc 的值增加 1。这样就可以在 nav.pclk 数组中持续追踪和存储卫星钟差信息。
end
nav.pclk(nav.nc).clk(sat,1)=clk;
nav.pclk(nav.nc).std(sat,1)=std;%这段代码的作用是将输入的时钟值和标准差赋值给卫星钟差信息数组 nav.pclk 中对应位置的元素。通过 nav.nc 索引获取当前的卫星钟差信息,然后根据 sat 索引指定具体的卫星,将对应的时钟值 clk 和标准差 std 赋值给相应的字段。这样就更新了卫星钟差信息的时钟值和标准差。
end
if nav.nc<size(nav.pclk,1)
nav.pclk(nav.nc+1:end,:)=[];
end
%该代码段的作用是检查 nav.pclk 数组的大小是否超过 nav.nc 的值(即卫星钟差信息的数量),如果超出,则删除多余的行,以确保数组的大小与 nav.nc 的值相匹配。通过这样的操作,只保留最新的卫星钟差信息,并删除旧的无效信息。
fclose(fid);
fprintf('over\n');
return
在读取精密星历文件的时候主要处理两部分,一部分是精密星历头文件,另外一部分是主体文件。
1.4天线校正文件(.atx)
调用readatx()函数
function nav=readatx(nav,fname)
global glc gls
freq=0; pcv0=gls.pcv; pcvs=gls.pcvs;
NMAX=10000; pcvs.pcv=repmat(gls.pcv,NMAX,1);
idx=strfind(fname,glc.sep); fname0=fname(idx(end)+1:end);
fprintf('Info:reading atx file %s...',fname0);
fid=fopen(fname);%打开天线文件
%这段代码通过设置变量、重复数组、提取字符串以及打印消息等操作,为读取 ATX 文件做好了准备。最后,使用 fopen 函数打开了文件,并返回文件句柄 fid,可以用于后续的文件读取操作。
while ~feof(fid)
line=fgets(fid);
if size(line,2)<=60 || ~isempty(strfind(line(61:end),'COMMENT')),continue;end%找到相关的标签是“COMMENT”。
if ~isempty(strfind(line(61:end),'START OF ANTENNA'))
pcv=pcv0;%如果第61-最后一列中发现有标签“START OF ANTENNA”,将pcv0的值赋给pcv。
end
if ~isempty(strfind(line(61:end),'END OF ANTENNA'))
if pcvs.n+1>size(pcvs.pcv,1)
pcvs.pcv(pcvs.n+1:pcvs.n+NMAX)=repmat(gls.pcv,NMAX,1);
end
pcvs.pcv(pcvs.n+1)=pcv;
pcvs.n=pcvs.n+1;
continue;
end
%代码段会检查读取到的 ATX 文件中的行,如果某行包含字符串 “END OF ANTENNA”,则在特定条件下进行数据处理和更新。这些操作包括扩展数组大小、赋值操作和变量的更新。最后,使用 continue 语句跳过当前循环,继续执行下一轮循环。
%这段代码用于检测并处理读取到的 ATX 文件中的行,如果某行包含字符串 “END OF ANTENNA”,则执行一系列的操作。下面是代码的含义:
if ~isempty(strfind(line(61:end),'TYPE / SERIAL NO'))%查看标签“TYPE / SERIAL NO”
pcv.type=line(1:20);%pcv.type 是一个变量,表示天线相位中心变量的类型
pcv.code=line(21:40);%pcv.code 是一个变量,表示天线的编码(code)
prn=fix(str2num(pcv.code(2:3)));
%这行代码对字符串 pcv.code 进行处理,提取其中的第2个和第3个字符,并将其转换为数值类型,然后使用 fix 函数对结果进行向下取整。最后将得到的结果赋值给变量 prn。
%fix() 是一个MATLAB(或类似语言)中的函数,用于将输入的数值向下取整到最接近的较小整数。
if ~prn,continue;end
if strcmp(pcv.code(4:11),' ')%pcv.code(4:11):这个索引操作提取了 pcv.code 字符串的第4个到第11个字符,即一个8个字符长的子字符串。' ':这个字符串是由8个空格组成的,用于与 pcv.code(4:11) 的子字符串进行比较。
pcv.sat=satid2no(pcv.code);%如果 pcv.code(4:11) 的子字符串与8个空格相等,则调用 satid2no() 函数,并将 pcv.code 作为参数传入。satid2no() 函数可能用于根据 pcv.code 获取卫星的编号或其他相关信息,并将其赋值给 pcv.sat 变量
end
elseif ~isempty(strfind(line(61:end),'VALID FROM'))%有效起始日期(值得是)
time=str2time(line(1:43));%将含有“有效起始日期”标签的第一列到43列存储到time中。
pcv.ts=time; %将time值赋给pcv.ts
continue;
elseif ~isempty(strfind(line(61:end),'VALID UNTIL'))%有效截止日期
time=str2time(line(1:43));%将含有“有效截止日期”标签的第一列到43列存储到time中。
pcv.te=time; %将time值赋给pcv.te
continue;
elseif ~isempty(strfind(line(61:end),'DAZI'))%“DAZI” 代表着差分方位角(Differential Azimuth)。方位角是指某一点相对于参考点的方向角度。在差分GPS中,使用DAZI来表示接收机与参考站之间的方位角差异。
pcv.dazi=str2num(line(3:8)); %#ok将第3-8列的字符以str2num的形式赋值给pcv.dazi
continue;
elseif ~isempty(strfind(line(61:end),'ZEN1 / ZEN2 / DZEN'))%“ZEN1” 和 “ZEN2” 代表测站天顶角(Zenith Angle),即接收机天线指向接收到信号的角度。这些参数用于估计接收机的卫星可见性和观测品质,对于实时运动定位和静态测量都很重要。
%“DZEN” 指的是测站天顶角的变化。在RTKLIB的解算过程中,由于卫星位置的变化或者天线运动引起的变化,接收机的测站天顶角也会发生变化。“DZEN” 用于表示测站天顶角的变化量。
pcv.zen1=str2num(line(3:8)); %#ok第3-8列的值以str2num的形式赋值给pcv.zen1
pcv.zen2=str2num(line(9:14)); %#ok第9-14列的值以str2num的形式赋值给pcv.zen2
pcv.dzen=str2num(line(15:20)); %#ok第15-20列的值以str2num的形式赋值给pcv.zen3
continue;
elseif ~isempty(strfind(line(61:end),'START OF FREQUENCY'))%表示的是开始的频率
%它检查文本行中是否包含"START OF FREQUENCY"的标识,如果是的话,就表示接下来要处理的是频率信息。
[f,count]=sscanf(line(5:6),'%d');
if count<1,continue;end
[csys,count]=sscanf(line(4),'%c');%“csys” 是一个表示坐标系统的变量
if count<1,continue;end
if csys=='G'
%如果坐标系统是 ‘G’,则将频率值直接赋给 freq。
freq=f;
elseif csys=='R'
%如果坐标系统是 ‘R’,则将频率值加上glc.NFREQ(频率数量)后赋给 freq。
freq=f+glc.NFREQ;
elseif csys=='E'%主要是根据频率的特性进行计算
%如果坐标系统是 ‘E’,根据频率值的不同情况进行赋值:
if f==1,freq=f+2*glc.NFREQ;
%如果频率值为 1,将其加上 2 乘以 glc.NFREQ 后赋给 freq。
elseif f==5,freq=2+2*glc.NFREQ;
%如果频率值为 5,将其加上 2 乘以 glc.NFREQ 后再加 2 后赋给 freq。
elseif f==6,freq=3+2*glc.NFREQ;
%如果频率值为 6,将其加上 2 乘以 glc.NFREQ 后再加 3 后赋给 freq。
else ,freq=0;
%其他情况下,将 freq 赋为 0。
end
elseif csys=='C'
%如果坐标系统是 ‘C’,则将频率值加上 3 乘以 glc.NFREQ 后赋给 freq。
freq=f+3*glc.NFREQ;
elseif csys=='J'
%如果坐标系统是 ‘J’,根据频率值的不同情况进行赋值:
if f<5,freq=f+4*glc.NFREQ;
%如果频率值小于 5,将其加上 4 乘以 glc.NFREQ 后赋给 freq。
elseif f==5,freq=3+4*glc.NFREQ;
%如果频率值为 5,将其加上 4 乘以 glc.NFREQ 后再加 3 后赋给 freq。
else ,freq=0;
%其他情况下,将 freq 赋为 0。
end
else
freq=0;
%如果坐标系统不属于上述情况,则将 freq 赋为 0。
end%最终,根据不同的坐标系统和频率值,通过计算得到对应的频率索引值 freq。
elseif ~isempty(strfind(line(61:end),'END OF FREQUENCY'))%表示的是结束的频率
freq=0; %此时将频率为0
elseif ~isempty(strfind(line(61:end),'NORTH / EAST / UP'))%北、东与垂直方向
%if freq<1||freq>glc.NFREQ,continue;end
[neu,count]=decodef(line,3);%调用decodef()函数
if count<3,continue;end
if freq<1,continue;end
if pcv.sat~=0%如果卫星不为0的时候
pcv.off(freq,:)=neu;%将neu值赋给pcv.off()
else
pcv.off(freq,:)=[neu(2),neu(1),neu(3)];%其他的话将neu(2)、neu(1)与neu(3)的值赋给pcv.off()
end
elseif ~isempty(strfind(line,'NOAZI'))%在句子里面查找NOAZI,找到的话执行下面的操作
if freq<1,continue;end
dd=(pcv.zen2-pcv.zen1)/pcv.dzen+1;%其中,pcv.zen2 和 pcv.zen1 是探测器的 Zenith 角度,pcv.dzen 表示 Zenith 角度的间隔。
if dd~=myRound(dd)||dd<=1
%fprintf('error');
continue;
end
if pcv.dazi==0
[val,count]=decodef(line(9:end),fix(dd));%调用decodef()函数
%decodef() 是一个用于解码字符串并返回浮点数数组的函数。
if count<=0
%fprintf('error');
continue
elseif count~=fix(dd)
%fprintf('error');
continue
end
pcv.var(freq,1:count)=val;
else
id=fix(360/pcv.dazi)+1;
%计算 id 的值,该值为 360/pcv.dazi 的取整值加 1。
for i=1:id
line=fgets(fid);
[val,j]=decodef(line(9:end),fix(dd));
if j<=0,continue
elseif j~=fix(dd),continue%fix()函数用于取整数部分
end
pcv.var(freq,(i-1)*fix(dd)+1:(i-1)*fix(dd)+j)=val;
end%复制数组:可以使用这个语句将一个数组的值复制到另一个数组的特定位置。
end%pcv.var(1, 6:9) = val; % 将数组 val 的值复制到 pcv.var 的第一个行中的索引 6 到 9 的元素
end%更新部分数据:可以使用这个语句将特定位置的元素更新为新的值。
%pcv.var(2, 3:6) = val; % 将数组 val 的值更新到 pcv.var 的第二行中的索引 3 到 6 的元素
end
fclose(fid);
fprintf('over\n');
if pcvs.n<size(pcvs.pcv,1)%如果 pcvs.n 的值小于 pcvs.pcv 数组的行数(即数组的大小),则执行条件判断内的语句。
pcvs.pcv(pcvs.n+1:end,:)=[];%在条件判断内部,pcvs.pcv(pcvs.n+1:end,:)=[]; 用于删除 pcvs.pcv 数组中从 pcvs.n+1 行到末尾的所有行。
end
nav.ant_para=pcvs;%pcvs的结果赋值给nav.ant_para里面
return
1.5差分码偏差文件(.DCB)
调用readdcb()函数
function nav=readdcb(nav,obs,fname)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%only read dcb for GPS or GLONASS(CODE procuct)
%convert CODE product to consistent with CAS product
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global glc
type=0; DCBPAIR=glc.DCBPAIR;%导入dcbpair的参数
idx=strfind(fname,glc.sep); fname0=fname(idx(end)+1:end);
fprintf('Info:reading DCB file %s...',fname0);%读到该文件的时候,在命令行输出“Info:reading DCB file %s...”
fid=fopen(fname);
%nav.cbias
%双差码偏差参数通常用于补偿不同频率的伪距观测之间的差异。将这个差异转化为以米为单位的数值,可以更准确地对卫星导航系统进行校正。
while ~feof(fid)
line=fgets(fid);
if ~isempty(strfind(line,'DIFFERENTIAL (P1-P2) CODE BIASES')),type=1;
elseif ~isempty(strfind(line,'DIFFERENTIAL (P1-C1) CODE BIASES')),type=2;
elseif ~isempty(strfind(line,'DIFFERENTIAL (P2-C2) CODE BIASES')),type=3;
end%将三个文件p1p2、p1c1与p2c2文件,将三种文件对应不同的类型,第一种类型是type=1;第二种类型是type=2,第三种类型是type=3.
if type==0,continue;end
if isempty(line),continue;end
if line(1)~='G'&&line(1)~='R',continue;end
buff=strsplit(line);%strsplit() 是一个字符串操作函数,通常用于将一个字符串分割成多个子字符串,并将这些子字符串存储在一个数组中。
%buff 一般指的是用于存储或处理数据的缓冲区或缓冲器。
cbias=str2double(line(27:35));%将第27-35列的值以str2double的形式赋值给cbias(abias表示的是载波偏差)。
%对于接收机来说
if strcmp(buff(1),'G')||strcmp(buff(1),'R')%如果存储的系统的类型为G或者R的时候。
if ~isfield(obs,'sta'),continue;end
if strcmp(buff(1),'G'),j=1;else,j=2;end
if strcmpi(obs.sta.name,char(buff(2)))
nav.rbias(j,type)=cbias*1E-9*glc.CLIGHT;%nav.rbias(j,type) = cbias * 1E-9 * glc.CLIGHT:如果上述条件满足,将计算得到的值 cbias 乘以 1E-9(表示1倍 1−910 −9,即纳秒转为秒的缩放因子),再乘以 glc.CLIGHT(光速)并赋值给 nav.rbias(j,type)。
end%rbias 是一个术语,表示接收机的接收机钟差(Receiver Clock Bias)。
%这段代码片段的功能是,如果 obs 中有 sta 字段,然后根据 buff 字符串中的第一个字符确定变量 j 的值,最后通过比较字符串 obs.sta.name 与 buff(2) 来决定是否给 nav.rbias(j,type) 赋值。
else
%对于卫星来说
sat=satid2no(char(buff(1)));%调用的是satid2no()函数,该函数用于将卫星的标识符(satid)转换为卫星编号(sats)。
if sat~=0
[sys,~]=satsys(sat);%调用的是satsys()函数,输入参数主要包括sat。
if type==1%如果是第一类型,p1p2的文件。
dcb_pair='C1W-C2W'; %P1-P2
%“dcb_pair” 是一个表示双差码偏差(Differential Code Bias)的字符串,例如 ‘C1W-C2W’。
%在这个字符串中,‘C1W’ 表示的是 P1 码,‘C2W’ 表示的是 P2 码。P1 码和 P2 码都是 GPS 卫星导航的伪距测量中的两个频率。
if sys==glc.SYS_GPS%如果系统是GPS系统的话
for i=1:glc.MAXDCBPAIR%使用for循环遍历双差码偏差参数数组DCBPAIR的每个元素。
if strcmp(dcb_pair,DCBPAIR(sys,i)),break;end%在每次循环中,通过使用strcmp函数将dcb_pair与DCBPAIR(sys, i)进行比较。如果它们相等,则执行break语句跳出循环。
end
%计算导航系统为GPS的卫星的双差码偏差。
nav.cbias(sat,i)=cbias*1E-9*glc.CLIGHT;%对于给定的卫星编号sat和双差码偏差参数索引i,将卫星卫星导航系统的dcbs数组中的元素设置为cbias乘以1E-9(表示以纳秒为单位)乘以光速glc.CLIGHT。
elseif sys==glc.SYS_GLO%如果系统是GLO系统话
dcb_pair='C1P-C2P'; %P1-P2
%“dcb_pair” 是一个表示双差码偏差(Differential Code Bias)的字符串,例如 ‘C1P-C2P’。
for i=1:glc.MAXDCBPAIR%%使用for循环遍历双差码偏差参数数组DCBPAIR的每个元素。
if strcmp(dcb_pair,DCBPAIR(sys,i)),break;end%在每次循环中,通过使用strcmp函数将dcb_pair与DCBPAIR(sys, i)进行比较。如果它们相等,则执行break语句跳出循环。
end
%计算导航系统GLO的卫星的双差码偏差
nav.cbias(sat,i)=cbias*1E-9*glc.CLIGHT;%;%对于给定的卫星编号sat和双差码偏差参数索引i,将卫星卫星导航系统的dcbs数组中的元素设置为cbias乘以1E-9(表示以纳秒为单位)乘以光速glc.CLIGHT。
end
elseif type==2%如果是p1c1类型的文件的话
dcb_pair='C1C-C1W'; %C1-P1
%“dcb_pair” 是一个表示双差码偏差(Differential Code Bias)的字符串,例如 ‘C1C-C1W’。
if sys==glc.SYS_GPS%如果系统是GPS系统的话
for i=1:glc.MAXDCBPAIR%%%使用for循环遍历双差码偏差参数数组DCBPAIR的每个元素。
if strcmp(dcb_pair,DCBPAIR(sys,i)),break;end%在每次循环中,通过使用strcmp函数将dcb_pair与DCBPAIR(sys, i)进行比较。如果它们相等,则执行break语句跳出循环。
end
%计算导航系统为GPS的卫星的双差码偏差。
nav.cbias(sat,i)=-cbias*1E-9*glc.CLIGHT;
%nav.cbias(sat, i) = -cbias * 1E-9 * glc.CLIGHT; 这一行代码的作用是将导航系统为GPS的卫星的双差码偏差的负值赋值给nav.cbias(sat, i)。
elseif sys==glc.SYS_GLO%如果系统是GLO系统的话
dcb_pair='C1C-C1P'; %C1-P1
%“dcb_pair” 是一个表示双差码偏差(Differential Code Bias)的字符串,例如 ‘C1C-C1P’。
for i=1:glc.MAXDCBPAIR%%使用for循环遍历双差码偏差参数数组DCBPAIR的每个元素。
if strcmp(dcb_pair,DCBPAIR(sys,i)),break;end%%在每次循环中,通过使用strcmp函数将dcb_pair与DCBPAIR(sys, i)进行比较。如果它们相等,则执行break语句跳出循环。
end
%计算导航系统GLO的卫星的双差码偏差
nav.cbias(sat,i)=-cbias*1E-9*glc.CLIGHT; %note the minus sign !
%nav.cbias(sat, i) = -cbias * 1E-9 * glc.CLIGHT; 这一行代码的作用是将导航系统为GLO的卫星的双差码偏差的负值赋值给nav.cbias(sat, i)。
end
elseif type==3%如果是p2c2类型的文件的话
dcb_pair='C2W-C2L'; %P2-C2
%“dcb_pair” 是一个表示双差码偏差(Differential Code Bias)的字符串,例如 ‘C2W-C2L’。
if sys==glc.SYS_GPS%如果卫星系统为GPS的话
for i=1:glc.MAXDCBPAIR%使用for循环遍历双差码偏差参数数组DCBPAIR的每个元素。
if strcmp(dcb_pair,DCBPAIR(sys,i)),break;end%在每次循环中,通过使用strcmp函数将dcb_pair与DCBPAIR(sys, i)进行比较。如果它们相等,则执行break语句跳出循环。
end
%计算卫星系统GPS的双差码偏差
nav.cbias(sat,i)=cbias*1E-9*glc.CLIGHT;
%在这个代码中,将导航系统为GPS的卫星的双差码偏差赋值给nav.cbias(sat, i)。具体而言,将cbias乘以1E-9(表示以纳秒为单位)乘以光速glc.CLIGHT的值,然后将结果赋给nav.cbias(sat, i)。
elseif sys==glc.SYS_GLO%如果卫星系统为GLO的话
dcb_pair='C2C-C2P'; %C2-P2
%“dcb_pair” 是一个表示双差码偏差(Differential Code Bias)的字符串,例如 ‘C2C-C2P’。
for i=1:glc.MAXDCBPAIR%使用for循环遍历双差码偏差参数数组DCBPAIR的每个元素。
if strcmp(dcb_pair,DCBPAIR(sys,i)),break;end%在每次循环中,通过使用strcmp函数将dcb_pair与DCBPAIR(sys, i)进行比较。如果它们相等,则执行break语句跳出循环。
end
%计算导航系统GLO的卫星的双差码偏差
nav.cbias(sat,i)=-cbias*1E-9*glc.CLIGHT; %note the minus sign !
%将cbias(双差码偏差值)乘以1E-9(表示以纳秒为单位)乘以光速glc.CLIGHT的值,并将结果赋值为nav.cbias(sat, i)的相反数。这样做的目的可能是根据特定的需求或公式要求,将双差码偏差的值取负。
end
end
end
end
end
fclose(fid);
fprintf('over\n');
return
1.6多系统差分码偏差文件(.BSX)
function nav=readdcb_mgex(nav,opt,obsr,fname)%定义一个名叫readdcb_mgex()函数,输入参数为nav,opt,obsr,fname,输出参数为nav
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%read dcb for MGEX (include GPS GLONASS GALILEO BDS QZSS)
%if CODE_DCB can be used for GPS and GLONASS,don't use CAS_DCB for GPS and
%GLONASS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%在示例的时候,一般需要的文件是CAS0MGXRAP_20190870000_01D_01D_DCB.BSX
global glc
type=0; time=str2time(opt.ts); DCBPAIR=glc.DCBPAIR;
if time.time==0
time=obsr.data(1).time;
end
idx=strfind(fname,glc.sep); %这段代码的目的是从字符串fname中提取出文件名(fname0)。具体而言,使用strfind函数在fname中查找glc.sep(可能是分隔符)第一次出现的位置,并将其索引存储在idx中。
fname0=fname(idx(end)+1:end);%然后,通过idx(end)+1,提取出最后一个分隔符的后一个位置作为文件名在fname中的起始位置。最后,使用起始位置和字符串切片操作,将起始位置到字符串末尾的部分提取出来,并存储在fname0中,即提取出了文件名。
fprintf('Info:reading DCB_MGEX file %s...',fname0);%读取DCB_MGEX文件(读到文件的名字)后,在命令行窗口输出“Info:reading DCB_MGEX file %s...”
fid=fopen(fname);%打开DCB_MGEX文件
while ~feof(fid)
line=fgets(fid);
if ~isempty(strfind(line,'+BIAS/SOLUTION')),type=1;end %进行循环,如果出现标签“+BIAS/SOLUTION”的时候,type为1
if ~isempty(strfind(line,'-BIAS/SOLUTION')),break;end%如果出现标签“-BIAS/SOLUTION”的话,结束操作
if ~type,continue;end%如果没有类型的话,结束操作
if ~strcmp(line(2:4),'DSB')||~strcmp(line(16:19),' '),continue;end
%如果第2-4列中出现DSB的话,或者第16-19列为空的话(' '),继续操作。
year =fix(str2num(line(36:39)));%#ok
%将第36到39列的数据以str2num的形式取整赋值到year上
doy_s=fix(str2num(line(41:43)));%#ok
%将第41到43列的数据以str2num的形式取整赋值到doy_s上----开始的年积日
doy_e=fix(str2num(line(56:58)));%#ok
%将第56到58列的数据以str2num的形式取整赋值到doy_e上----结束的年积日
if year<=50,year=year+2000;end
%如果是年分小于等于50的话,将其值加上2000
time_s=ydoy2time(year,doy_s); %调用ydoy2time()函数,输入参数为年份以及开始的年积日,输出为开始的时间
time_e=ydoy2time(year,doy_e);%调用ydoy2time()函数,输入参数为年份以及结束的年积日,输出为结束的时间
if ~(timediff(time,time_s)>=0&&timediff(time,time_e)<0),continue;end
%这一行代码的含义是,如果时间 time 不在时间范围 time_s 和 time_e 内(即时间早于 time_s 或晚于等于 time_e),那么就执行 continue,即跳过当前迭代,进入下一次迭代。
%换句话说,这行代码的作用是忽略不符合时间范围要求的情况,继续处理下一个迭代。
satid=line(12:14);%satid表示的是卫星编号,从第12列到14列(PRN)
sat=satid2no(satid);%调用satid2no函数,输入参数为satid,输出为sat
[sys,prn]=satsys(sat);%调用的是satsys()函数
dcb_pair=[line(26:28),'-',line(31:33)];%计算双差码偏差,以第一个为例,是将(C1C-C1W)
cbias=str2num(line(72:92)); %#ok
%以第一行为例,将第72-92列以str2num的形式赋值给cbias(为-0.9690)。----cbias(双差码偏差值)
if sys==glc.SYS_GPS%在上面的两个标签之间,如果系统为GPS系统的话
if ~nav.no_CODE_DCB,continue;end
%Code Differential Code Bias (简称 CODE_DCB)是卫星导航系统中的一个校准参数,用于修正由于 GPS 和 GLONASS 等系统的信号传播路径中的大气影响引起的误差。
for i=1:glc.MAXDCBPAIR
if strcmp(dcb_pair,DCBPAIR(sys,i)),break;end
end%%在每次循环中,通过使用strcmp函数将dcb_pair与DCBPAIR(sys, i)进行比较。如果它们相等,则执行break语句跳出循环。
nav.cbias(sat,i)=cbias*1E-9*glc.CLIGHT;%这行代码将 cbias 乘以 1E-9(即十的负九次方,相当于除以十亿),然后再乘以 glc.CLIGHT,其中 glc.CLIGHT 表示光速。最终的结果赋值给 nav.cbias(sat,i)。
elseif sys==glc.SYS_GLO%如果系统为GLO的话
if ~nav.no_CODE_DCB,continue;end%如果 nav.no_CODE_DCB 的值为假(false),即表示码差分码偏差不可用,
for i=1:glc.MAXDCBPAIR%使用for循环遍历双差码偏差参数数组DCBPAIR的每个元素。
if strcmp(dcb_pair,DCBPAIR(sys,i)),break;end%在每次循环中,通过使用strcmp函数将dcb_pair与DCBPAIR(sys, i)进行比较。如果它们相等,则执行break语句跳出循环。
end
nav.cbias(sat,i)=cbias*1E-9*glc.CLIGHT;
%在这个代码中,将导航系统为GPS的卫星的双差码偏差赋值给nav.cbias(sat, i)。具体而言,将cbias乘以1E-9(表示以纳秒为单位)乘以光速glc.CLIGHT的值,然后将结果赋给nav.cbias(sat, i)。
elseif sys==glc.SYS_GAL%如果系统为GAL的话
for i=1:glc.MAXDCBPAIR%使用for循环遍历双差码偏差参数数组DCBPAIR的每个元素。
if strcmp(dcb_pair,DCBPAIR(sys,i)),break;end%在每次循环中,通过使用strcmp函数将dcb_pair与DCBPAIR(sys, i)进行比较。如果它们相等,则执行break语句跳出循环。
end
nav.cbias(sat,i)=cbias*1E-9*glc.CLIGHT;
%%在这个代码中,将导航系统为GPS的卫星的双差码偏差赋值给nav.cbias(sat, i)。具体而言,将cbias乘以1E-9(表示以纳秒为单位)乘以光速glc.CLIGHT的值,然后将结果赋给nav.cbias(sat, i)。
elseif sys==glc.SYS_BDS%如果卫星系统是BSD的话
for i=1:glc.MAXDCBPAIR%使用for循环遍历双差码偏差参数数组DCBPAIR的每个元素。
if prn>18 %BDS-3%如果prn大于18的话
if strcmp(dcb_pair,DCBPAIR(sys,i))
if i==2,i=glc.BD3_C2IC6I;end %#ok
break;%如果 dcb_pair 和 DCBPAIR(sys,i) 相等,则在 i 等于 2 的情况下,将 i 赋值为 glc.BD3_C2IC6I,然后跳出当前循环或循环结构
end
else %BDS-2
if strcmp(dcb_pair,DCBPAIR(sys,i)),break;end%在每次循环中,通过使用strcmp函数将dcb_pair与DCBPAIR(sys, i)进行比较。如果它们相等,则执行break语句跳出循环。
end
end
nav.cbias(sat,i)=cbias*1E-9*glc.CLIGHT;%cbias为纳秒
%这段代码将一个浮点数变量 cbias 乘以 1E-9(即 1 除以 10 的 9 次方)和光速值(由常数 glc.CLIGHT 表示),并将结果赋值给 nav.cbias(sat, i)。也就是说,这段代码用于将一个以纳秒为单位的钟差偏差(cbias)转换为以米为单位的钟差偏差,并将结果存储在 nav.cbias 数组的相应位置(sat 和 i 为数组的索引)中。
elseif sys==glc.SYS_QZS%如果系统是QZSS系统的话
for i=1:glc.MAXDCBPAIR%使用for循环遍历双差码偏差参数数组DCBPAIR的每个元素。
if strcmp(dcb_pair,DCBPAIR(sys,i)),break;end %在每次循环中,通过使用strcmp函数将dcb_pair与DCBPAIR(sys, i)进行比较。如果它们相等,则执行break语句跳出循环。
end
nav.cbias(sat,i)=cbias*1E-9*glc.CLIGHT;
%在这个代码中,将导航系统为GPS的卫星的双差码偏差赋值给nav.cbias(sat, i)。具体而言,将cbias乘以1E-9(表示以纳秒为单位)乘以光速glc.CLIGHT的值,然后将结果赋给nav.cbias(sat, i)。
end
end
fclose(fid);
fprintf('over\n');
return
1.7地球参考框架文件(.ERP)
function nav=readerp(nav,fname)
%定义可readerp()函数,输入的参数为nav与fname,输出的参数书为nav
global glc gls%代码声明了三个全局变量,即 glc、gls
flag=0; n=0; erp=gls.erp;
idx=strfind(fname,glc.sep); fname0=fname(idx(end)+1:end);%在索引的时候发现.erp的文件名字的时候
fprintf('Info:reading erp file %s...',fname0);%在命令行输出“Info:reading erp file %s...”
fid=fopen(fname);%使用 fopen 函数打开文件.erp(后缀)。
% read data
while ~feof(fid)
line=fgets(fid,266);
if ~isempty(strfind(line,'MJD')),flag=1;end%if ~isempty(strfind(line,'MJD')),flag=1;end 判断 line 中是否包含字符串 ‘MJD’。如果包含,则将 flag 置为 1。
if flag==0,continue;end
if flag==1
line=fgets(fid); %#ok<NASGU>
flag=2; continue;
end
if flag==2
value=str2double(strsplit(line));%将每一列的元素进行存储,村相互到value里面
%value=str2double(strsplit(line)); 将字符串 line 按空格进行分割,并将每个分割后的字符串转换为数值,并将结果存储在 value 数组中。
erp.data(n+1).mjd =value(1);%.mjd:Modified Julian Date(修正朱利安日期),用来表示观测时间。
% 将 value 数组中的第一个元素赋值给 erp.data 结构体数组的成员 mjd,其中 n+1 为数组的索引。
erp.data(n+1).xp =value(2)*1e-6*glc.AS2R;%.xp:地球自转的 X 方向极移,即地球自转轴在国际参考坐标系下相对于地理坐标系的偏移。
%glc.AS2R 是一个常量,代表角秒转弧度的比例因子
erp.data(n+1).yp =value(3)*1e-6*glc.AS2R;%.yp:地球自转的 Y 方向极移,即地球自转轴在国际参考坐标系下相对于地理坐标系的偏移。
erp.data(n+1).ut1_utc=value(4)*1e-7;%.ut1_utc:UT1 与 UTC 之间的差距,用来校正地球自传的时间单位。
erp.data(n+1).lod =value(5)*1e-7;%.lod:地球自转率(Length of Day),表示地球自转的角速度。
erp.data(n+1).xpr =value(13)*1e-6*glc.AS2R;%.xpr:极移预测的 X 方向误差。
erp.data(n+1).ypr =value(14)*1e-6*glc.AS2R;%.ypr:极移预测的 Y 方向误差。
n=n+1;
end
end
fclose(fid);
fprintf('over\n');
erp.n=n;
erp.nmax=n;
nav.erp=erp;
return
1.8坐标时序文件(.blq)
function nav=readblq(nav,obsr,obsb,fname)%定义一个readblq()函数,输入参数有基准站的观测数据、流动站的观测数据、文件名字以及与观测文件相关的一个结构体nav,输出的参数包括nav
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%read ocean tide loading parameter(otlp) for rover and base
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global glc
idx=strfind(fname,glc.sep); fname0=fname(idx(end)+1:end);
fprintf('Info:reading blq file %s...',fname0);%在命令行输出“Info:reading blq file ocnload.blq...”
ns=0;%ns表示的是卫星的数量
if obsr.n>0,sta(1)=obsr.sta;ns=1;end
if obsb.n>0,sta(2)=obsb.sta;ns=2;end
if ns==0,return;end
for i=1:ns
fid=fopen(fname);%打开(.blq)文件
sta_name=upper(sta(i).name);%将测站的名字都转换成大写的形式,赋值给sta_name。
%upper()函数将字符串中的所有字符转换为大写字母形式,并返回转换后的结果。
while ~feof(fid)
line=fgets(fid);
%line = fgets(fid); 表示从打开的文件中读取一行文本,并将其存储在 line 变量中。
%每次调用 fgets 函数,它会读取文件中的下一行,并将结果作为一个字符串返回。
%如果文件中已经没有更多的行可读取,fgets 函数将返回一个空指针或空字符串。
if ~isempty(strfind(line,'-cmc')),break;end%如果文件中找到“-cmc”的话,结束操作
end
while ~feof(fid)
line=fgets(fid);
%line = fgets(fid); 表示从打开的文件中读取一行文本,并将其存储在 line 变量中。
%每次调用 fgets 函数,它会读取文件中的下一行,并将结果作为一个字符串返回。
%如果文件中已经没有更多的行可读取,fgets 函数将返回一个空指针或空字符串。
if ~isempty(strfind(line,'$$'))||size(line,2)<=2,continue;end %如果在文件不为空,且找到“$$”的话
name=upper(line(3:6));%以第一个为例,主要将CEBR转为大写字母
if isempty(strfind(sta_name,name)),continue;end%如果在文件中找到测站的话
[otlp,fid,stat]=decode_blqdata(fid);%调用decode_blqdata()函数,输入参数包括fid,输出的为otlp、fid与stat
if stat==1%使用条件语句,判断变量 stat 是否等于1。如果条件成立,执行接下来的代码块,否则跳过。
nav.otlp(:,:,i)=otlp';%将 otlp 的转置值赋给 nav.otlp(:,:,i)。这段代码假设 nav.otlp 是一个三维数组,将 otlp 的转置值放入 nav.otlp 的第 i 个维度。
end
end
fclose(fid);%关闭fid
end
fprintf('over\n');
return
在进行读取数据的时候调用的是decode_blqdata()函数,在该函数中,对将数据以矩阵的形式存储下来。erer
function [otlp,fid,stat]=decode_blqdata(fid)%定义了一个名叫decode_blqdata()函数,主要是用以解决blq文件数据存储问题,输入参数包括fid,输出参数为otlp,fid,stat
stat=0; n=0;
otlp=zeros(6,11);
%将变量 stat 和 n 初始化为0;创建一个大小为 6x11 的全零矩阵,并将其赋给变量 otlp。
while ~feof(fid)
line=fgets(fid);
if ~isempty(strfind(line,'$$'))||size(line,2)<=2,continue;end%if ~isempty(strfind(line, '$$')) || size(line, 2) <= 2:如果 line 中包含 $$ 或者 line 的长度小于等于2,则执行 continue;,跳过当前循环的剩余代码。
v=str2double(strsplit(line));%v = str2double(strsplit(line));:将 line 按空格分割,并将得到的子字符串转换为数值,并存储在向量 v 中
idx=~isnan(v);
value=v(idx);
for i=1:11%刚好是文件的一个测站的数据
otlp(n+1,i)=value(i);
end
n=n+1;
if n==6
stat=1;
return;
end
end
return
2.设置文件输出
2.1模式设置(PPP)
function rtk=initoutfile(rtk,opt,file,obsr)%定义了一个initoutfile()函数,输入参数为rtk、opt(结构体)、文件 以及观测数据,输出参数为rtk
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Copyright(c) 2020-2025, by Kai Chen, All rights reserved.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global glc
if opt.ins.mode==glc.GIMODE_OFF
if opt.mode==glc.PMODE_SPP
mode='SPP';
elseif opt.mode==glc.PMODE_PPP_KINEMA
mode='PPP_K';
end
end
%根据模型是什么,选择相应的模式----PPP
2.2路径及结果文件设置
fullname=char(file.obsr);%glc.sep为/分隔符 输出的结构文件的名字是根据观测文件的名字
idx1=find(fullname==glc.sep); idx2=find(fullname=='.');
outfilename=[fullname(idx1(end)+1:idx2(end)-1),'_',mode,'.pos'];%设置输出文件名字
fullpath=which('main.m');%result的路径是与main.m在同一个文件夹里面
[path,~,~]=fileparts(fullpath);
pathname=[path,glc.sep,'result',glc.sep];
rtk.outfile=[pathname,outfilename];%完整的路径
3.设置卫星和接收者的天线参数
% 设置对于卫星和接收者的天线参数
if opt.mode~=glc.PMODE_SPP%首先判断优化模式是否为SPP(Single Point Positioning),如果不是,则进行以下操作。
nav=L2_L5pcv_copy(nav);%将导航数据(nav)中的L2-L5频段的天线相位中心变量(pcv)拷贝一份(L2_L5pcv_copy(nav)),并同时将观测数据(obsr、obsb)中的站点信息(sta)提取出来。
if isfield(obsr,'sta'),stas(1)=obsr.sta;end
if isfield(obsb,'sta'),stas(2)=obsb.sta;end
time0.time=obsr.data(1,1);time0.sec=obsr.data(1,2);%将观测数据的第一条数据的时间信息(time0)提取出来,包括时间(time0.time)和秒(time0.sec)。
[nav,opt]=setpcv(time0,opt,nav,nav.ant_para,nav.ant_para,stas);%调用setpcv函数来设置这些天线参数。setpcv函数会在指定的时间(time0)和导航数据(nav)中,根据天线参数(nav.ant_para)设置卫星和接收者的天线参数,并将设置后的天线参数保存在opt结构体中。
end
3.1调用L2_L5pcv_copy()函数
nav=L2_L5pcv_copy(nav)
这一行代码的作用是创建一个导航数据(nav)的副本,并将副本命名为 L2_L5pcv;
L2_L5pcv_copy
函数是 RTKLIB 中的一个函数,它用于复制导航数据的 L2-L5 频段的天线相位中心变量(PCV)参数。这些 PCV 参数描述了卫星和接收机的天线特性,包括天线相位中心、天线相位中心变化随着观测角度的变化等信息。通过调用 L2_L5pcv_copy
函数,可以将导航数据中的 L2-L5 频段的 PCV 参数复制到新的导航数据副本 L2_L5pcv
中。这样做的目的是为了在设置卫星和接收机的天线参数时,使用 L2-L5 频段的 PCV 参数进行计算和处理。
总结起来,nav=L2_L5pcv_copy(nav)
这行代码的作用是复制导航数据的 L2-L5 频段的 PCV 参数,并将副本保存在 L2_L5pcv
中,以备后续使用。
function nav=L2_L5pcv_copy(nav)%定义了一个名叫L2_L5pcv_copy()函数,输入参数为nav,输出参数为nav
%L2_L5pcv_copy() 函数是用于复制 L2 和 L5 相位中心变化 (PCV) 数据的函数。
if nav.ant_para.n==0,return;end%L2_L5pcv_copy() 函数的作用是将 L2 和 L5 频率的相位中心变化(PCV)数据从一个数据结构复制到另一个数据结构中。
for i=1:nav.ant_para.n%定义了一个循环,从 1 到 nav.ant_para.n 进行迭代。这里的 nav.ant_para.n 可能表示接收机天线的个数或者与 L2 和 L5 相位中心变化数据相关的其他参数。
if norm(nav.ant_para.pcv(i).off(3,:))>0,continue;end%判断第 i 个天线的 off 字段的第 3 行是否大于零。off 可能表示相位中心变化的偏移量。如果大于零,则跳过当前迭代,继续下一个循环迭代
nav.ant_para.pcv(i).off(3,:)=nav.ant_para.pcv(i).off(2,:);%将第 i 个天线的 off 字段的第 2 行的值赋给第 3 行。这个赋值操作可能是复制 L2 相位中心变化到 L5 相位中心变化。
nav.ant_para.pcv(i).var(3,:)=nav.ant_para.pcv(i).var(2,:);%将第 i 个天线的 var 字段的第 2 行的值赋给第 3 行。这个赋值操作可能是复制 L2 相位中心变化的方差到 L5 相位中心变化的方差。
end
return
3.2调用setpcv()函数
[nav,opt]=setpcv(time0,opt,nav,nav.ant_para,nav.ant_para,stas)
这行代码的作用是调用 setpcv
函数来设置卫星和接收机的天线参数,并将结果保存在 nav
和 opt
变量中。
opt.pcvr=repmat(gls.pcv,mode,1);%这行代码的作用是将一个长度为 mode 的向量 gls.pcv 在第一个维度上复制多次,并将结果存储在 opt.pcvr 变量中。
repmat
是 MATLAB 中的函数,用于将输入向量或矩阵在指定的维度上进行复制。gls.pcv
是一个原始的向量。mode
是一个指定复制次数的整数值。
通过调用 repmat(gls.pcv,mode,1)
,我们将 gls.pcv
向量在第一个维度上复制了 mode
次,也就是将其复制了 mode
行。第二个参数 1
表示在第二个维度上不进行复制,即保持原始向量的列数不变。
function mask=set_sysmask(opt)%设置set_sysmask()函数,输入参数opt。
mask(1:5)=0;
%这个函数可以用于根据输入参数的不同设置系统屏蔽状态,进而控制系统的行为。
if ~isempty(strfind(opt,'G')),mask(1)=1;end%如果输入参数 opt 中包含字母 “G”,则将 mask 中的第一个元素设置为 1。
if ~isempty(strfind(opt,'R')),mask(2)=1;end%如果输入参数 opt 中包含字母 “R”,则将 mask 中的第一个元素设置为 2。
if ~isempty(strfind(opt,'E')),mask(3)=1;end%如果输入参数 opt 中包含字母 “E”,则将 mask 中的第一个元素设置为 3。
if ~isempty(strfind(opt,'C')),mask(4)=1;end%如果输入参数 opt 中包含字母 “C”,则将 mask 中的第一个元素设置为 4。
if ~isempty(strfind(opt,'J')),mask(5)=1;end%如果输入参数 opt 中包含字母 “J”,则将 mask 中的第一个元素设置为 5。
return
获取卫星的编号,卫星的数量,卫星系统的pcvs时间等。
if pcvs.n>0
for i=1:glc.MAXSAT%最大的卫星数量(1:150)
[sys,~]=satsys(i);%使用函数 satsys(i) 获取卫星编号 i 对应的系统(如 GPS、GLONASS、Galileo 等)。
if mask(sys)==0,continue;end
[pcv,stat]=searchpcv(i,'',time,pcvs);%使用函数 searchpcv 根据卫星编号 i、空字符串 ''、时间 time 和 pcvs(可能是存储卫星天线参数的结构体或变量)来搜索卫星的天线参数。所得的天线参数存储在变量 pcv 中,同时还会返回一个状态 stat。
if stat==0%如果 stat 的值为 0,表示未找到卫星的天线参数,则执行以下操作:
id=satno2id(i);%使用函数 satno2id(i) 将卫星编号 i 转换为卫星标识符(通常是卫星的名称或编号)
nav.pcvs(i)=gls.pcv;%将 gls.pcv 赋值给导航数据 nav 中的 pcvs(表示卫星天线参数),以备后续使用。
fprintf('Warning:%s,have no antenna parameters!\n',id);%打印警告信息,提示该卫星缺少天线参数。
continue;
end
nav.pcvs(i)=pcv;
调用satsys()函数
function [sys,prn]=satsys(sat)%定义了一个satsys()函数,该函数输入参数包括sat,输出的参数包括sys,prn。
global glc
sys=0;
if sat<=0||sat>glc.MAXSAT
sat=0;%如果卫星编号小于等于0或大于最大卫星数glc.MAXSAT,则将sys设置为0,sat设置为0。
elseif sat<=glc.NSATGPS%%%%%%在进行PPP的时候主要是进行的是这一块
sys=glc.SYS_GPS; %如果卫星编号小于等于GLONASS卫星数glc.NSATGPS,则将sys设置为GPS,sat设置为卫星编号加上最小PRN编号减1。
sat=sat+glc.MINPRNGPS-1;
elseif (sat-glc.NSATGPS)<=glc.NSATGLO
sys=glc.SYS_GLO; %如果卫星编号减去GPS卫星数glc.NSATGPS后小于等于GLONASS卫星数glc.NSATGLO,则将sys设置为GLONASS,sat设置为卫星编号减去GPS卫星数再加上最小PRN编号减1。
sat=(sat-glc.NSATGPS)+glc.MINPRNGLO-1;
elseif (sat-glc.NSATGPS-glc.NSATGLO)<=glc.NSATGAL
sys=glc.SYS_GAL; %如果卫星编号减去GPS卫星数glc.NSATGPS再减去GLONASS卫星数glc.NSATGLO后小于等于GALILEO卫星数glc.NSATGAL,则将sys设置为GALILEO,sat设置为卫星编号减去GPS卫星数再减去GLONASS卫星数再加上最小PRN编号减1。
sat=(sat-glc.NSATGPS-glc.NSATGLO)+glc.MINPRNGAL-1;
elseif (sat-glc.NSATGPS-glc.NSATGLO-glc.NSATGAL)<=glc.NSATBDS
sys=glc.SYS_BDS; %如果卫星编号减去GPS卫星数glc.NSATGPS再减去GLONASS卫星数glc.NSATGLO再减去GALILEO卫星数glc.NSATGAL后小于等于BDS卫星数glc.NSATBDS,则将sys设置为BDS,sat设置为卫星编号减去GPS卫星数再减去GLONASS卫星数再减去GALILEO卫星数再加上最小PRN编号减1。
sat=(sat-glc.NSATGPS-glc.NSATGLO-glc.NSATGAL)+glc.MINPRNBDS-1;
elseif (sat-glc.NSATGPS-glc.NSATGLO-glc.NSATGAL-glc.NSATBDS)<=glc.NSATQZS
sys=glc.SYS_QZS; %如果卫星编号减去GPS卫星数glc.NSATGPS再减去GLONASS卫星数glc.NSATGLO再减去GALILEO卫星数glc.NSATGAL再减去BDS卫星数glc.NSATBDS后小于等于QZSS卫星数glc.NSATQZS,则将sys设置为QZSS,sat设置为卫星编号减去GPS卫星数再减去GLONASS卫星数再减去GALILEO卫星数再减去BDS卫星数再加上最小PRN编号减1。
sat=(sat-glc.NSATGPS-glc.NSATGLO-glc.NSATGAL-glc.NSATBDS)+glc.MINPRNQZS-1;
else
sat=0;
end
prn=sat;%将sat值赋值给prn(每颗卫星唯一的编码值)。
return
4.数据处理过程
调用gnss_processor()函数
if opt.ins.mode==glc.GIMODE_OFF%根据所选择的模式将进行相应的处理过程
% gnss
gnss_processor(rtk,opt,obsr,obsb,nav);%调用 gnss_processor 函数对 GNSS 数据进行处理。
end
4.1设置输出的进条度
ti=0;
%创建一个名为"GINav"的进度条,并在创建进度条时设置进度条窗口的位置和大小。
hbar=waitbar(0,'Processing...','Name','GINav', 'CreateCancelBtn', 'delete(gcbf);');%代码中的waitbar函数用于创建进度条窗口,第一个参数"0"表示初始进度为0,第二个参数"‘Processing…’“表示在进度条窗口中显示的文本信息为"Processing…”,第三个参数"‘Name’,‘GINav’“表示设置进度条窗口的名称为"GINav”,第四个参数"‘CreateCancelBtn’, ‘delete(gcbf);’"表示为进度条窗口创建一个取消按钮,并设置点击取消按钮时关闭进度条窗口。
H=get(0,'ScreenSize'); w=600; h=450; x=H(3)/2-w/2; y=H(4)/2-h/2; %通过函数get获取屏幕的尺寸,然后计算出进度条窗口的位置和大小,使用figure函数创建一个新的图形窗口,并将图形窗口的位置和大小设置为之前计算得到的值。
hfig=figure;set(gcf,'Position',[x y w h]);%使用set函数将当前图形窗口设置为gcf (get current figure),并将位置和大小设置为之前计算得到的值。
4.2初始化rtk结构
调用initrtk()函数
rtk=initrtk(rtk,opt);
function rtk=initrtk(rtk,opt)%定义一个initrtk()函数,输入参数为rtk与opt输出参数为rtk
global glc
rtk.opt=opt;%将结构体opt赋值给rtk.opt
rtk.opt.ts=str2time(rtk.opt.ts); %将rtk.opt.ts(表示的是开始的时间)以str2time的形式赋值给rtk.opt.ts
rtk.opt.te=str2time(rtk.opt.te);%rtk.opt.te(表示的是结束的时间)以str2time的形式赋值给rtk.opt.te
rtk.mask=set_sysmask(opt.navsys);%设置卫星系统的掩码,用以确定对应卫星系统
%在进行PPP的时候不进行这一快
if opt.mode<=glc.PMODE_STATIC
if opt.ionoopt==glc.IONOOPT_IFLC,NF=1;else,NF=opt.nf;end
if opt.dynamics==1,NP=9;else,NP=3;end
if opt.ionoopt==glc.IONOOPT_EST,NI=glc.MAXSAT;else,NI=0;end
if opt.tropopt<glc.TROPOPT_EST,NT=0;elseif opt.tropopt==glc.TROPOPT_EST,NT=2;else,NT=6;end
if opt.glomodear~=2,NL=0;else,NL=2;end
if opt.mode<=glc.PMODE_DGNSS,NB=0;else,NB=NF*glc.MAXSAT;end
% 参数的数量
NR=NP+NI+NT+NL; NX=NR+NB;
rtk.NF=NF;rtk.NP=NP;rtk.NI=NI;rtk.NT=NT;rtk.NL=NL;rtk.NB=NB;
% index of parameter
II=NP;ITR=NP+NI;ITB=NP+NI+NT/2;IL=NP+NI+NT;IB=NR;
rtk.ii=II; rtk.itr=ITR; rtk.itb=ITB; rtk.il=IL; rtk.ib=IB;
else%IFLC 方法是在两个或多个频率的载波相位观测值之间构建线性组合,利用它们之间的差异来消除电离层延迟的影响。
if opt.ionoopt==glc.IONOOPT_IFLC,NF=1;else,NF=opt.nf;end%“NF” 是指导航系统中的频率数量或频段的数量。
if opt.dynamics==1,NP=9;else,NP=3;end%opt.dynamics 是一个选项变量,用于设置动态模型的处理方式。
NC=glc.NSYS;%表示支持的卫星数量
% GLONASS
if opt.gloicb==0, NICB=0;%这是一个变量,用于存储根据 opt.gloicb 的值确定的国际代码偏差类型
elseif opt.gloicb==glc.GLOICB_LNF,NICB=1;%表示伽利略的国际代码偏差(inter-frequency code bias)类型为线性(linear)。
else, NICB=2;
end
if opt.tropopt<glc.TROPOPT_EST,NT=0;elseif opt.tropopt==glc.TROPOPT_EST,NT=1;else,NT=3;end%glc.TROPOPT_EST用以在GLONASS系统中估计对流层延迟
if opt.ionoopt==glc.IONOOPT_EST,NI=glc.MAXSAT;else,NI=0;end
if opt.nf>=3,ND=1;else,ND=0;end
% 参数的数量
NR=NP+NC+NICB+NT+NI+ND; NB=NF*glc.MAXSAT; NX=NR+NB;
rtk.NF=NF;rtk.NP=NP;rtk.NC=NC;rtk.NICB=NICB;rtk.NT=NT;rtk.NI=NI;rtk.ND=ND;rtk.NB=NB;
% 参数的索引
IC=NP; IICB=NP+NC; IT=NP+NC+NICB; II=NP+NC+NICB+NT; ID=NP+NC+NICB+NT+NI; IB=NR;
rtk.ic=IC; rtk.iicb=IICB; rtk.it=IT; rtk.ii=II; rtk.id=ID; rtk.ib=IB;
end
rtk.nx=NX; rtk.na=NR;%“rtk.nx” 表示状态向量的维度。它用于确定 RTK 过程中所估计的状态量的数量。这些状态量可能包括位置、速度、钟差等信息。 “rtk.na” 表示滤波状态向量的维度。它用于确定滤波器的状态估计向量的尺寸。在滤波过程中,滤波器通过更新状态向量来估计和更新系统的状态,而状态向量是通过测量和预测进行更新的。
rtk.x=zeros(rtk.nx,1);%通过将rtk.x初始化为大小为rtk.nx行1列的零向量,将rtk对象中的x变量初始化为大小为rtk.nx的变量向量。---“x” 是状态估计向量,它包含了定位过程中估计的各个状态量的值。
rtk.P=zeros(rtk.nx,rtk.nx);%通过将rtk.P初始化为大小为rtk.nx行rtk.nx列的零矩阵,将rtk对象中的P变量初始化为大小为rtk.nx的协方差矩阵。---P是状态估计的协方差矩阵,表示各个状态量之间的误差协方差。
rtk.xa=zeros(rtk.na,1);%通过将rtk.xa初始化为大小为rtk.na行1列的零向量,将rtk对象中的xa变量初始化为大小为rtk.na的变量向量。---“xa” 是滤波中的滤波状态向量(滤波状态估计),表示滤波器输出的最优估计状态。
rtk.Pa=zeros(rtk.na,rtk.na);%通过将rtk.Pa初始化为大小为rtk.na行rtk.na列的零矩阵,将rtk对象中的Pa变量初始化为大小为rtk.na的协方差矩阵。---“Pa” 是滤波中的滤波协方差矩阵,表示滤波器输出的状态估计的协方差。
return