共识卫星数量

残差阵,数量ns*nf*2(每颗卫星的每个频率×2是伪距和载波)

dz设计矩阵,数量ns*3

每个频点的伪距和载波频率一样而且伪距也用不到频率是载波求波长时用的

N2等模糊度

blbaseline(x,rtk-rb,dr);ecef2pos(x,posu);

Ri是参考卫星方差RJ是非参考卫星方差。

用于后面计算双差后的观测值方差--协方差阵

Rimat(ns*nf*22,1);

初始化伪距和载波相位残差为0每个卫星每个频点的用于以后记录双差残差

double

MAXSAT最大卫星数不管观没观测到都初始化。

这也造成了空间上的浪费

NFREQ最大频率

{rtk-ssat[i].resp[j]rtk-ssat[i].resc[j]0.0;}

对流层和电离层处理暂不考虑

{im[i](ionmapf(posu,azeliu[i]*2)ionmapf(posr,azelir[i]*2))/2.0;}if

{tropu[i]prectrop(rtk-sol.time,posu,0,azeliu[i]*2,opt,x,dtdxui*3);tropr[i]prectrop(rtk-sol.time,posr,1,azelir[i]*2,opt,x,dtdxri*3);}}

for

m0:GPS/SBS,1:GLO,2:GAL,3:BDS,4:QZS,5:IRN

二重循环根据所选择的定位模式进行相应的频点遍历。

nf*2是因为每个系统的每个频点(*2

包括伪距和载波)

是载波(动态定位模式及以后)那我们从0--2*nf因为同时使用伪距和载波残差

for

(fopt-modePMODE_DGPS?0:nf;fnf*2;f)

ns---共识卫星数量

这个就是找参考卫星的在这里是把高度角最高的那个作为参考卫星。

我们也可以改写代码设置其他的规则。

参考卫星下标

for

continue;//当前共识卫星的系统和当前遍历的系统是否一样if

(!validobs(iu[j],ir[j],f,nf,y))

continue;

(i0||azel[1iu[j]*2]azel[1iu[i]*2])

又是一个大循环遍历所有的共识卫星这次是要计算双差残差和设计矩阵。

组双差的话两台接收机是要都观测到相同的卫星的所以需要共识卫星处理。

for

continue;sysirtk-ssat[sat[i]-1].sys;//参考卫星系统sysjrtk-ssat[sat[j]-1].sys;//非参考卫星系统freqifreq[f%nfiu[i]*nf];freqjfreq[f%nfiu[j]*nf];if

(!test_sys(sysj,m))

(!validobs(iu[j],ir[j],f,nf,y))

(H)

(y[fiu[i]*nf*2]-y[fir[i]*nf*2])

v[nv](y[fiu[i]*nf*2]-y[fir[i]*nf*2])-(y[fiu[j]*nf*2]-y[fir[j]*nf*2]);

设计矩阵e[

因为在relpos()函数中站间单差[先是算基准站后是流动站]设计矩阵都填到e里面了。

所以大小如此

设计矩阵非参考卫星-参考卫星

我的猜想是可能一样因为由计算系数的公式可知这个与卫星位置和接收机概略位置有关

(H)

{Hi[k]-e[kiu[i]*3]e[kiu[j]*3];}}

Hi[

{v[nv]-CLIGHT/freqi*x[IB(sat[i],f,opt)]-CLIGHT/freqj*x[IB(sat[j],f,opt)];if

(H)

CLIGHT/freqi;Hi[IB(sat[j],f,opt)]-CLIGHT/freqj;}}else

{//根据站星双差方差还要减去(模糊度的derT)见理论。

而在y[]的计算中还没有考虑模糊度在这里给补上//x[]

模糊度v[nv]-x[IB(sat[i],f,opt)]-x[IB(sat[j],f,opt)];if

(H)

哪个卫星哪个频点的模糊度位置索引不懂Hi[IB(sat[j],f,opt)]-1.0;}}}

Hi[IB(sat[i],f,opt)]

这里是在补设计矩阵参考卫星系数为-1。

而这里为1后续应该在某个计算中会改掉

(fnf)

rtk-ssat[sat[j]-1].resp[f-nf]v[nv];/*

test

(opt-maxinno0.0fabs(v[nv])opt-maxinno)

{if

{rtk-ssat[sat[i]-1].rejc[f];rtk-ssat[sat[j]-1].rejc[f];}errmsg(rtk,outlier

rejected

v%.3f)\n,sat[i],sat[j],fnf?L:P,f%nf1,v[nv]);continue;}

Ri[nv]varerr(sat[i],sysi,azel[1iu[i]*2],bl,dt,f,opt);Rj[nv]varerr(sat[j],sysj,azel[1iu[j]*2],bl,dt,f,opt);

(opt-modePMODE_DGPS)

rtk-ssat[sat[i]-1].vsat[f]rtk-ssat[sat[j]-1].vsat[f]1;}else

{rtk-ssat[sat[i]-1].vsat[f-nf]rtk-ssat[sat[j]-1].vsat[f-nf]1;}

独特的编码由此码可知双差的具体信息。

哪两个卫星作差载波还是伪距第几个频率

vflg[nv](sat[i]16)|(sat[j]8)|((fnf?0:1)4)|(f%nf);

至此共识卫星遍历结束

(opt-modePMODE_MOVEBconstbl(rtk,x,P,v,H,Ri,Rj,nv))

{vflg[nv]34;nb[b];}

{R[ki(kj)*nv]Ri[ki](ij?Rj[ki]:0.0);

//参考GNSS理论16.2讲。

非对角线参考卫星方差}}trace(5,R\n);

validobs()

有载波就不管伪距了解这是f(0--2*nf)的情况随着循环进行伪距最终会遍历到(nf--2*nf)检查y[fi*nf*2]!0

检查基准站共识卫星的载波残差fnf||(y[f-nfi*nf*2]!0.0y[f-nfj*nf*2]!0.0)

f-nf为真短路原则当fnf,伪距残差,载波也要有不然伪距不可用(注释)y[fi*nf*2]

载波观测

y[fi*nf*2]!0.0y[fj*nf*2]!0.0(fnf||(y[f-nfi*nf*2]!0.0y[f-nfj*nf*2]!0.0));

全部代码

bl,dr[3],posu[3],posr[3],didxi0.0,didxj0.0,*im;double

*tropr,*tropu,*dtdxr,*dtdxu,*Ri,*Rj,freqi,freqj,*HiNULL;int

i,j,k,m,f,nv0,nb[NFREQ*4*22]{0},b0,sysi,sysj,nfNF(opt);trace(3,ddres

dt%.1f

ns%d\n,dt,rtk-nx,ns);//计算基线长度blbaseline(x,rtk-rb,dr);ecef2pos(x,posu);

基准站和流动站方差计算以后协方差阵对一些中间变量进行初始化将双差伪距残差和双差载波相位残差初始化为0*/Rimat(ns*nf*22,1);

Rjmat(ns*nf*22,1);

将双差伪距残差和双差载波相位残差初始化为0每个卫星每个频点的为什么没有jNFREQ*2,因为伪距和载波分为resp和resc了*/for

(i0;iMAXSAT;i)

{rtk-ssat[i].resp[j]rtk-ssat[i].resc[j]0.0;}/*

compute

如果卡尔曼滤波中包含对流层状态量调用prectrop函数计算基站和移动站的对流层延迟湿分量存在tropu[i]和tropur[i]中如果卡尔曼滤波器包含电离层状态量调用ionmapf函数分别计算基站和移动站处的投影函数*/for

(i0;ins;i)

{im[i](ionmapf(posu,azeliu[i]*2)ionmapf(posr,azelir[i]*2))/2.0;}if

{tropu[i]prectrop(rtk-sol.time,posu,0,azeliu[i]*2,opt,x,dtdxui*3);tropr[i]prectrop(rtk-sol.time,posr,1,azelir[i]*2,opt,x,dtdxri*3);}}/*

对各个系统以及载波相位、伪距分别进行循环处理双重循环遍历每个系统的每个频点(*2

包括伪距和载波)哪种定位模式是伪距那我们从nf---2*nf因为伪距残差是在y[]后半是载波(动态定位模式及以后)那我们从0--2*nf因为同时使用伪距和载波*/for

(m0;m6;m)

m0:GPS/SBS,1:GLO,2:GAL,3:BDS,4:QZS,5:IRN

*/for

(fopt-modePMODE_DGPS?0:nf;fnf*2;f)

{/*

elevation选择最大的高度角卫星作为参考卫星当然我们也可以自己设置规则ssat_t

ssat[MAXSAT];

所有的共识卫星(全系统)每次都会被遍历一遍我们选出所需的系统然后再*/for

(i-1,j0;jns;j)

continue;//当前共识卫星的系统和当前遍历的系统是否一样if

(!validobs(iu[j],ir[j],f,nf,y))

continue;

(i0||azel[1iu[j]*2]azel[1iu[i]*2])

(i0)

做双差遍历共识卫星一个历元一个历元处理输入流动站和基准站的数据流动站卫星残差

基准站卫星残差

每颗卫星每个频点的载波和伪距y[...............

................]站间单差两个接收机(a,b)观测同一颗卫星(p)作差(M1)。

卫星p肯定是它们的共识卫星星间双差参考卫星站间单差(M2),M1-M2*/for

(j0;jns;j)

continue;sysirtk-ssat[sat[i]-1].sys;//参考卫星系统sysjrtk-ssat[sat[j]-1].sys;//流动freqifreq[f%nfiu[i]*nf];freqjfreq[f%nfiu[j]*nf];if

(!test_sys(sysj,m))

(!validobs(iu[j],ir[j],f,nf,y))

continue;//检验数据有效yy[]//设计矩阵初始化nx---待解量的个数if

(H)

共识--参考卫星在基准站残差中的索引(y[fiu[i]*nf*2]-y[fir[i]*nf*2])

共识参考卫星站间单差*/v[nv](y[fiu[i]*nf*2]-y[fir[i]*nf*2])-(y[fiu[j]*nf*2]-y[fir[j]*nf*2]);/*

partial

因为在relpos()函数中站间单差[先是算基准站后是流动站]设计矩阵都填到e里面了。

所以大小如此设计矩阵非参考卫星-参考卫星因为我们求得是流动站位置所以只用iu[]即可至于不同频率,这里只有单频*/if

(H)

{Hi[k]-e[kiu[i]*3]e[kiu[j]*3];}}/*

ionospheric

{didxi(fnf?-1.0:1.0)*im[i]*SQR(FREQ1/freqi);didxj(fnf?-1.0:1.0)*im[j]*SQR(FREQ1/freqj);v[nv]-didxi*x[II(sat[i],opt)]-didxj*x[II(sat[j],opt)];if

(H)

didxi;Hi[II(sat[j],opt)]-didxj;}}/*

tropospheric

(opt-tropoptTROPOPT_EST||opt-tropoptTROPOPT_ESTG)

{v[nv]-(tropu[i]-tropu[j])-(tropr[i]-tropr[j]);for

(k0;k(opt-tropoptTROPOPT_ESTG?1:3);k)

{if

(dtdxu[ki*3]-dtdxu[kj*3]);Hi[IT(1,opt)k]-(dtdxr[ki*3]-dtdxr[kj*3]);}}/*

phase-bias

{v[nv]-CLIGHT/freqi*x[IB(sat[i],f,opt)]-CLIGHT/freqj*x[IB(sat[j],f,opt)];if

(H)

CLIGHT/freqi;Hi[IB(sat[j],f,opt)]-CLIGHT/freqj;}}else

{//根据站星双差方差还要减去(模糊度的derT)见理论。

而在y[]的计算中还没有考虑模糊度在这里给补上//x[]

模糊度v[nv]-x[IB(sat[i],f,opt)]-x[IB(sat[j],f,opt)];if

(H)

哪个卫星哪个频点的模糊度位置索引不懂Hi[IB(sat[j],f,opt)]-1.0;}}}//将计算好的残差存到ssat结构体if

(fnf)

rtk-ssat[sat[j]-1].resp[f-nf]v[nv];/*

test

(opt-maxinno0.0fabs(v[nv])opt-maxinno)

{if

{rtk-ssat[sat[i]-1].rejc[f];rtk-ssat[sat[j]-1].rejc[f];}errmsg(rtk,outlier

rejected

v%.3f)\n,sat[i],sat[j],fnf?L:P,f%nf1,v[nv]);continue;}/*

(single-differenced)

*/Ri[nv]varerr(sat[i],sysi,azel[1iu[i]*2],bl,dt,f,opt);Rj[nv]varerr(sat[j],sysj,azel[1iu[j]*2],bl,dt,f,opt);/*

set

rtk-ssat[sat[i]-1].vsat[f]rtk-ssat[sat[j]-1].vsat[f]1;}else

{rtk-ssat[sat[i]-1].vsat[f-nf]rtk-ssat[sat[j]-1].vsat[f-nf]1;}trace(4,sat%3d-%3d

%s%d

%8.6f\n,sat[i],sat[j],fnf?L:P,f%nf1,v[nv],Ri[nv],Rj[nv]);//标记一下双差的具体信息。

哪两个卫星作差载波还是伪距第几个频率vflg[nv](sat[i]16)|(sat[j]8)|((fnf?0:1)4)|(f%nf);nb[b];//每个小频率有几颗卫星(残差有多少个)}b;//每个小频率}/*

end

(opt-modePMODE_MOVEBconstbl(rtk,x,P,v,H,Ri,Rj,nv))

(H)

tracemat(5,H,rtk-nx,nv,7,4);}/*

measurement

*/ddcov(nb,b,Ri,Rj,nv,R);free(Ri);

free(Rj);