当前位置：首页 > news >正文

多模卫星导航定位与应用-原理与实践（RTKLib）6

news 2025/10/8 5:35:26

若基站坐标不准，

基站的坐标误差，

伪距和相位双差

还需要广播星历，知道卫星的坐标。

第1个问题，共视卫星变了，双差模糊度也会变，法方程的维护问题

rtkpos.c

/* relpos()relative positioning ------------------------------------------------------
* args: rtk IO gps solution structure
obs I satellite observations 观测值
nu I # of user observations (rover)流动站
nr I # of ref observations (base)基站
nav I satellite navigation data 广播星历
*/
static int relpos(rtk_t *rtk, const obsd_t *obs, int nu, int nr,
const nav_t *nav)
{
prcopt_t *opt=&rtk->opt;
gtime_t time=obs[0].time;
double *rs,*dts,*var,*y,*e,*azel,*freq,*v,*H,*R,*xp,*Pp,*xa,*bias,dt;
int i,j,f,n=nu+nr,ns,ny,nv,sat[MAXSAT],iu[MAXSAT],ir[MAXSAT];
int info,vflg[MAXOBS*NFREQ*2+1],svh[MAXOBS*2];
int stat=rtk->opt.mode<=PMODE_DGPS?SOLQ_DGPS:SOLQ_FLOAT;
int nf=opt->ionoopt==IONOOPT_IFLC?1:opt->nf;

trace(3,"relpos : nu=%d nr=%d\n",nu,nr);

/* define local matrices, n=total observations, base + rover */
rs=mat(6,n); /* range to satellites */
dts=mat(2,n); /* satellite clock biases */
var=mat(1,n);
y=mat(nf*2,n);
e=mat(3,n);
azel=zeros(2,n); /* [az, el] */
freq=zeros(nf,n);

/* init satellite status arrays */
for (i=0;i<MAXSAT;i++) {
rtk->ssat[i].sys=satsys(i+1,NULL); /* gnss system */
for (j=0;j<NFREQ;j++) {
rtk->ssat[i].vsat[j]=0; /* valid satellite */
rtk->ssat[i].snr_rover[j]=0;
rtk->ssat[i].snr_base[j] =0;
}
}

//计算卫星位置、速度和钟差
/* compute satellite positions, velocities and clocks for base and rover */
satposs(time,obs,n,nav,opt->sateph,rs,dts,var,svh);

/* calculate [range - measured pseudorange] for base station (phase and code)
output is in y[nu:nu+nr], see call for rover below for more details */
trace(3,"base station:\n");

//计算基准站非差残差
if (!zdres(1,obs+nu,nr,rs+nu*6,dts+nu*2,var+nu,svh+nu,nav,rtk->rb,opt,
y+nu*nf*2,e+nu*3,azel+nu*2,freq+nu*nf)) {
errmsg(rtk,"initial base station position error\n");

free(rs); free(dts); free(var); free(y); free(e); free(azel); free(freq);
return 0;
}

//计算基准站和流动站的延迟
/* time diff between base and rover observations */
if (opt->intpref) {
/* time-interpolation of base residuals */
dt=intpres(time,obs+nu,nr,nav,rtk,y+nu*nf*2);
} else dt=timediff(time,obs[nu].time);
trace(3,"relpos : dt=%.3f\n",dt);

if (opt->mode!=PMODE_MOVEB) {
/* check if exceeded max age of differential */
rtk->sol.age=dt;
if (fabs(rtk->sol.age)>opt->maxtdiff) {
errmsg(rtk,"age of differential error (age=%.1f)\n",rtk->sol.age);
free(rs); free(dts); free(var); free(y); free(e); free(azel); free(freq);
return 1;
}
}

//选择共视卫星
/* select common satellites between rover and base-station */
if ((ns=selsat(obs,azel,nu,nr,opt,sat,iu,ir))<=0) {
errmsg(rtk,"no common satellite\n");

free(rs); free(dts); free(var); free(y); free(e); free(azel); free(freq);
return 0;
}

//卡尔曼滤波
/* update kalman filter states (pos,vel,acc,ionosp, troposp, sat phase biases) */
trace(4,"before udstate: x="); tracemat(4,rtk->x,1,NR(opt),13,4);
udstate(rtk,obs,sat,iu,ir,ns,nav);
trace(4,"after udstate x="); tracemat(4,rtk->x,1,NR(opt),13,4);

for (i=0;i<ns;i++) for (j=0;j<nf;j++) {

/* snr of base and rover receiver */
rtk->ssat[sat[i]-1].snr_rover[j]=obs[iu[i]].SNR[j];
rtk->ssat[sat[i]-1].snr_base[j] =obs[ir[i]].SNR[j];
}

/* initialize Pp,xa to zero, xp to rtk->x */
xp=mat(rtk->nx,1); Pp=zeros(rtk->nx,rtk->nx); xa=mat(rtk->nx,1);
matcpy(xp,rtk->x,rtk->nx,1);
matcpy(Pp,rtk->P,rtk->nx,rtk->nx);

ny=ns*nf*2+2;
v=mat(ny,1); H=zeros(rtk->nx,ny); R=mat(ny,ny); bias=mat(rtk->nx,1);

trace(3,"rover: dt=%.3f\n",dt);
for (i=0;i<opt->niter;i++) {
/* calculate zero diff residuals [range - measured pseudorange] for rover (phase and code)
output is in y[0:nu-1], only shared input with base is nav
obs = sat observations
nu = # of sats
rs = range to sats
dts = sat clock biases (rover)
svh = sat health flags
nav = sat nav data
xp = kalman states
opt = options
y = zero diff residuals (code and phase)
e = line of sight unit vectors to sats
azel = [az, el] to sats */

//流动站非差残差
if (!zdres(0,obs,nu,rs,dts,var,svh,nav,xp,opt,y,e,azel,freq)) {
errmsg(rtk,"rover initial position error\n");
stat=SOLQ_NONE;
break;
}
/* calculate double-differenced residuals and create state matrix from sat angles
O rtk->ssat[i].resp[j] = residual pseudorange error
O rtk->ssat[i].resc[j] = residual carrier phase error
I dt = time diff between base and rover observations
I Pp = covariance matrix of float solution
I sat = list of common sats
I iu,ir = user and ref indices to sats
I ns = # of sats
O v = double diff residuals (phase and code)
O H = partial derivatives
O R = double diff measurement error covariances
O vflg = list of sats used for dd */

//双差残差
if ((nv=ddres(rtk,obs,dt,xp,Pp,sat,y,e,azel,freq,iu,ir,ns,v,H,R,vflg))<4) {
errmsg(rtk,"not enough double-differenced residual, n=%d\n", nv);
stat=SOLQ_NONE;
break;
}
/* kalman filter measurement update, updates x,y,z,sat phase biases, etc
K=P*H*(H'*P*H+R)^-1
xp=x+K*v
Pp=(I-K*H')*P */
trace(3,"before filter x=");tracemat(3,rtk->x,1,9,13,6);
if ((info=filter(xp,Pp,H,v,R,rtk->nx,nv))) {
errmsg(rtk,"filter error (info=%d)\n",info);
stat=SOLQ_NONE;
break;
}
trace(3,"after filter x=");tracemat(3,xp,1,9,13,6);
trace(4,"x(%d)=",i+1); tracemat(4,xp,1,NR(opt),13,4);
}
/* calc zero diff residuals again after kalman filter update */
if (stat!=SOLQ_NONE&&zdres(0,obs,nu,rs,dts,var,svh,nav,xp,opt,y,e,azel,freq)) {

/* calc double diff residuals again after kalman filter update for float solution */
nv=ddres(rtk,obs,dt,xp,Pp,sat,y,e,azel,freq,iu,ir,ns,v,NULL,R,vflg);

/* validation of float solution, always returns 1, msg to trace file if large residual */
if (valpos(rtk,v,R,vflg,nv,4.0)) {

/* copy states */
matcpy(rtk->x,xp,rtk->nx,1);
matcpy(rtk->P,Pp,rtk->nx,rtk->nx);

/* update valid satellite status for ambiguity control */
rtk->sol.ns=0;
for (i=0;i<ns;i++) for (f=0;f<nf;f++) {
if (!rtk->ssat[sat[i]-1].vsat[f]) continue;
rtk->ssat[sat[i]-1].outc[f]=0;
if (f==0) rtk->sol.ns++; /* valid satellite count by L1 */
}
/* too few valid phases */
if (rtk->sol.ns<4) stat=SOLQ_DGPS;
}
else stat=SOLQ_NONE;
}
/* resolve integer ambiguity by LAMBDA */
if (stat==SOLQ_FLOAT) {
/* if valid fixed solution, process it */
if (manage_amb_LAMBDA(rtk,bias,xa,sat,nf,ns)>1) {

/* find zero-diff residuals for fixed solution */
if (zdres(0,obs,nu,rs,dts,var,svh,nav,xa,opt,y,e,azel,freq)) {

/* post-fit residuals for fixed solution (xa includes fixed phase biases, rtk->xa does not) */
nv=ddres(rtk,obs,dt,xa,Pp,sat,y,e,azel,freq,iu,ir,ns,v,NULL,R,vflg);

/* validation of fixed solution, always returns valid */
if (valpos(rtk,v,R,vflg,nv,4.0)) {

/* hold integer ambiguity if meet minfix count */
if (++rtk->nfix>=rtk->opt.minfix) {
// Note that the modear needs to be fix-and-hold in
// order for glomodear fix-and-hold to be applied here,
// to prevent glomodear alone forcing fix-and-hold.
if (rtk->opt.modear==ARMODE_FIXHOLD)
holdamb(rtk,xa);
/* switch to kinematic after qualify for hold if in static-start mode */
if (rtk->opt.mode==PMODE_STATIC_START) {
rtk->opt.mode=PMODE_KINEMA;
trace(3,"Fix and hold complete: switch to kinematic mode\n");
}
}
stat=SOLQ_FIX;
}
}
}
}

/* save solution status (fixed or float) */
if (stat==SOLQ_FIX) {
for (i=0;i<3;i++) {
rtk->sol.rr[i]=rtk->xa[i];
rtk->sol.qr[i]=(float)rtk->Pa[i+i*rtk->na];
}
rtk->sol.qr[3]=(float)rtk->Pa[1];
rtk->sol.qr[4]=(float)rtk->Pa[1+2*rtk->na];
rtk->sol.qr[5]=(float)rtk->Pa[2];

if (rtk->opt.dynamics) { /* velocity and covariance */
for (i=3;i<6;i++) {
rtk->sol.rr[i]=rtk->xa[i];
rtk->sol.qv[i-3]=(float)rtk->Pa[i+i*rtk->na];
}
rtk->sol.qv[3]=(float)rtk->Pa[4+3*rtk->na];
rtk->sol.qv[4]=(float)rtk->Pa[5+4*rtk->na];
rtk->sol.qv[5]=(float)rtk->Pa[5+3*rtk->na];
}
}
else { /* float solution */
for (i=0;i<3;i++) {
rtk->sol.rr[i]=rtk->x[i];
rtk->sol.qr[i]=(float)rtk->P[i+i*rtk->nx];
}
rtk->sol.qr[3]=(float)rtk->P[1];
rtk->sol.qr[4]=(float)rtk->P[1+2*rtk->nx];
rtk->sol.qr[5]=(float)rtk->P[2];

if (rtk->opt.dynamics) { /* velocity and covariance */
for (i=3;i<6;i++) {
rtk->sol.rr[i]=rtk->x[i];
rtk->sol.qv[i-3]=(float)rtk->P[i+i*rtk->nx];
}
rtk->sol.qv[3]=(float)rtk->P[4+3*rtk->nx];
rtk->sol.qv[4]=(float)rtk->P[5+4*rtk->nx];
rtk->sol.qv[5]=(float)rtk->P[5+3*rtk->nx];
}
rtk->nfix=0;
}
trace(3,"sol_rr= ");tracemat(3,rtk->sol.rr,1,6,15,3);
/* save phase measurements */
for (i=0;i<n;i++) for (j=0;j<nf;j++) {
if (obs[i].L[j]==0.0) continue;
rtk->ssat[obs[i].sat-1].pt[obs[i].rcv-1][j]=obs[i].time;
rtk->ssat[obs[i].sat-1].ph[obs[i].rcv-1][j]=obs[i].L[j];
}
for (i=0;i<MAXSAT;i++) for (j=0;j<nf;j++) {
/* Don't lose track of which sats were used to try and resolve the ambiguities */
/* if (rtk->ssat[i].fix[j]==2&&stat!=SOLQ_FIX) rtk->ssat[i].fix[j]=1; */
if (rtk->ssat[i].slip[j]&LLI_SLIP) rtk->ssat[i].slipc[j]++;
/* Inc lock count if this sat used for good fix */
if (!rtk->ssat[i].vsat[j]) continue;
if (rtk->ssat[i].lock[j]<0||(rtk->nfix>0&&rtk->ssat[i].fix[j]>=2))
rtk->ssat[i].lock[j]++;
}
free(rs); free(dts); free(var); free(y); free(e); free(azel); free(freq);
free(xp); free(Pp); free(xa); free(v); free(H); free(R); free(bias);

if (stat!=SOLQ_NONE) rtk->sol.stat=stat;

return stat!=SOLQ_NONE;
}

//双差残差

/* double-differenced residuals and partial derivatives -----------------------------------
O rtk->ssat[i].resp[j] = residual pseudorange error
O rtk->ssat[i].resc[j] = residual carrier phase error
I rtk->rb= base location
I dt = time diff between base and rover observations
I x = rover pos & vel and sat phase biases (float solution)
I P = error covariance matrix of float states
I sat = list of common sats
I y = zero diff residuals (code and phase, base and rover)
I e = line of sight unit vectors to sats
I azel = [az, el] to sats
I iu,ir = user and ref indices to sats
I ns = # of sats
O v = double diff innovations (measurement-model) (phase and code)
O H = linearized translation from innovations to states (az/el to sats)
O R = measurement error covariances
O vflg = bit encoded list of sats used for each double diff */
static int ddres(rtk_t *rtk, const obsd_t *obs, double dt, const double *x,
const double *P, const int *sat, double *y, double *e,
double *azel, double *freq, const int *iu, const int *ir,
int ns, double *v, double *H, double *R, int *vflg)
{
prcopt_t *opt=&rtk->opt;
double bl,dr[3],posu[3],posr[3],didxi=0.0,didxj=0.0,*im;
double *tropr,*tropu,*dtdxr,*dtdxu,*Ri,*Rj,freqi,freqj,*Hi=NULL,df;
int i,j,k,m,f,nv=0,nb[NFREQ*NSYS*2+2]={0},b=0,sysi,sysj,nf=NF(opt);
int frq,code;

trace(3,"ddres : dt=%.4f ns=%d\n",dt,ns);

/* bl=distance from base to rover, dr=x,y,z components */
bl=baseline(x,rtk->rb,dr);
/* translate ecef pos to geodetic pos */
ecef2pos(x,posu); ecef2pos(rtk->rb,posr);

Ri=mat(ns*nf*2+2,1); Rj=mat(ns*nf*2+2,1); im=mat(ns,1);
tropu=mat(ns,1); tropr=mat(ns,1); dtdxu=mat(ns,3); dtdxr=mat(ns,3);

/* zero out residual phase and code biases for all satellites */
for (i=0;i<MAXSAT;i++) for (j=0;j<NFREQ;j++) {
rtk->ssat[i].resp[j]=rtk->ssat[i].resc[j]=0.0;
}
/* compute factors of ionospheric and tropospheric delay
- only used if kalman filter contains states for ION and TROP delays
usually insignificant for short baselines (<10km)*/
for (i=0;i<ns;i++) {
if (opt->ionoopt==IONOOPT_EST) {
im[i]=(ionmapf(posu,azel+iu[i]*2)+ionmapf(posr,azel+ir[i]*2))/2.0;
}
if (opt->tropopt>=TROPOPT_EST) {
tropu[i]=prectrop(rtk->sol.time,posu,0,azel+iu[i]*2,opt,x,dtdxu+i*3);
tropr[i]=prectrop(rtk->sol.time,posr,1,azel+ir[i]*2,opt,x,dtdxr+i*3);
}
}

//对不同导航系统，各自选择参考卫星，做双差
/* step through sat systems: m=0:gps/sbs,1:glo,2:gal,3:bds 4:qzs 5:irn*/
for (m=0;m<6;m++) {

/* step through phases/codes */
for (f=opt->mode>PMODE_DGPS?0:nf;f<nf*2;f++) {
frq=f%nf;code=f<nf?0:1;

/* find reference satellite with highest elevation, set to i */
for (i=-1,j=0;j<ns;j++) {
sysi=rtk->ssat[sat[j]-1].sys;
if (!test_sys(sysi,m) || sysi==SYS_SBS) continue;
if (!validobs(iu[j],ir[j],f,nf,y)) continue;
/* skip sat with slip unless no other valid sat */
if (i>=0&&rtk->ssat[sat[j]-1].slip[frq]&LLI_SLIP) continue;
if (i<0||azel[1+iu[j]*2]>=azel[1+iu[i]*2]) i=j;
}
if (i<0) continue;

//计算双差残差

/* calculate double differences of residuals (code/phase) for each sat */
for (j=0;j<ns;j++) {
if (i==j) continue; /* skip ref sat */
sysi=rtk->ssat[sat[i]-1].sys;
sysj=rtk->ssat[sat[j]-1].sys;
freqi=freq[frq+iu[i]*nf];
freqj=freq[frq+iu[j]*nf];
if (freqi<=0.0||freqj<=0.0) continue;
if (!test_sys(sysj,m)) continue;
if (!validobs(iu[j],ir[j],f,nf,y)) continue;

if (H) {
Hi=H+nv*rtk->nx;
for (k=0;k<rtk->nx;k++) Hi[k]=0.0;
}

//两个接收站对两个卫星的双差

/* double-differenced measurements from 2 receivers and 2 sats in meters */
v[nv]=(y[f+iu[i]*nf*2]-y[f+ir[i]*nf*2])-
(y[f+iu[j]*nf*2]-y[f+ir[j]*nf*2]);

//流动站对两个卫星位置偏导数

/* partial derivatives by rover position, combine unit vectors from two sats */
if (H) {
for (k=0;k<3;k++) {
Hi[k]=-e[k+iu[i]*3]+e[k+iu[j]*3]; /* translation of innovation to position states */
}
}

//为了维护法方程的稳定性，列出每颗卫星对流层和电离层改正值，与PPP 一样

if (opt->ionoopt==IONOOPT_EST) {
/* adjust double-differenced measurements by double-differenced ionospheric delay term */
didxi=(code?-1.0:1.0)*im[i]*SQR(FREQL1/freqi);
didxj=(code?-1.0:1.0)*im[j]*SQR(FREQL1/freqj);
v[nv]-=didxi*x[II(sat[i],opt)]-didxj*x[II(sat[j],opt)];
if (H) {
Hi[II(sat[i],opt)]= didxi;
Hi[II(sat[j],opt)]=-didxj;
}
}
if (opt->tropopt>=TROPOPT_EST) {
/* adjust double-differenced measurements by double-differenced tropospheric delay term */
v[nv]-=(tropu[i]-tropu[j])-(tropr[i]-tropr[j]);
for (k=0;k<(opt->tropopt<TROPOPT_ESTG?1:3);k++) {
if (!H) continue;
Hi[IT(0,opt)+k]= (dtdxu[k+i*3]-dtdxu[k+j*3]);
Hi[IT(1,opt)+k]=-(dtdxr[k+i*3]-dtdxr[k+j*3]);
}
}
if (opt->mode > PMODE_DGPS && !code) {

//无电离组合IF法求整周模糊度
int ii = IB(sat[i], frq, opt);
int jj = IB(sat[j], frq, opt);
/* Adjust phase residual by double-differenced phase-bias term,
IB=look up index by sat&freq */
if (opt->ionoopt!=IONOOPT_IFLC) {
/* Phase-bias states are single-differenced so need to difference them */
v[nv]-=CLIGHT/freqi*x[ii]-CLIGHT/freqj*x[jj];
if (H) {
Hi[ii]= CLIGHT/freqi;
Hi[jj]=-CLIGHT/freqj;
}
}
else {
v[nv]-=x[ii]-x[jj];
if (H) {
Hi[ii]= 1.0;
Hi[jj]=-1.0;
}
}
}

/* adjust double-difference for glonass sats */
if (sysi==SYS_GLO&&sysj==SYS_GLO) {
if (rtk->opt.glomodear==GLO_ARMODE_AUTOCAL && frq<NFREQGLO) {
/* auto-cal method */
df=(freqi-freqj)/(f==0?DFRQ1_GLO:DFRQ2_GLO);
v[nv]-=df*x[IL(frq,opt)];
if (H) Hi[IL(frq,opt)]=df;
}
else if (rtk->opt.glomodear==GLO_ARMODE_FIXHOLD && frq<NFREQGLO) {
/* fix-and-hold method */
double icb=rtk->ssat[sat[i]-1].icbias[frq]*CLIGHT/freqi - rtk->ssat[sat[j]-1].icbias[frq]*CLIGHT/freqj;
v[nv]-=icb;
}
}

/* adjust double-difference for sbas sats */
if (sysj==SYS_SBS&&sysi==SYS_GPS) {
if (rtk->opt.glomodear==GLO_ARMODE_FIXHOLD && frq<NFREQ) {
/* fix-and-hold method */
double icb=rtk->ssat[sat[i]-1].icbias[frq]*CLIGHT/freqi - rtk->ssat[sat[j]-1].icbias[frq]*CLIGHT/freqj;
v[nv]-=icb;
}
}

/* save residuals */
if (code) rtk->ssat[sat[j]-1].resp[frq]=v[nv]; /* pseudorange */
else rtk->ssat[sat[j]-1].resc[frq]=v[nv]; /* carrier phase */

/* open up outlier threshold if one of the phase biases was just initialized */
double threshadj = 1;
if (opt->mode > PMODE_DGPS) {
// Open up outlier threshold if one of the phase biases was just initialized.
int ii = IB(sat[i], frq, opt);
int jj = IB(sat[j], frq, opt);
if (P[ii + rtk->nx * ii] == SQR(rtk->opt.std[0]) ||
P[jj + rtk->nx * jj] == SQR(rtk->opt.std[0]))
threshadj = 10;
}
/* if residual too large, flag as outlier */
if (fabs(v[nv])>opt->maxinno[code]*threshadj) {
rtk->ssat[sat[j]-1].vsat[frq]=0;
rtk->ssat[sat[j]-1].rejc[frq]++;
errmsg(rtk,"outlier rejected (sat=%3d-%3d %s%d v=%.3f)\n",
sat[i],sat[j],code?"P":"L",frq+1,v[nv]);
continue;
}

//单差误差协方差

/* single-differenced measurement error variances (m) */
Ri[nv] = varerr(sat[i], sysi, azel[1+iu[i]*2],
rtk->ssat[sat[i]-1].snr_rover[frq],
rtk->ssat[sat[i]-1].snr_base[frq],
bl,dt,f,opt,&obs[iu[i]]);
Rj[nv] = varerr(sat[j], sysj, azel[1+iu[j]*2],
rtk->ssat[sat[j]-1].snr_rover[frq],
rtk->ssat[sat[j]-1].snr_base[frq],
bl,dt,f,opt,&obs[iu[j]]);
/* increase variance if half cycle flags set */
if (!code&&(obs[iu[i]].LLI[frq]&LLI_HALFC)) Ri[nv]+=0.01;
if (!code&&(obs[iu[j]].LLI[frq]&LLI_HALFC)) Rj[nv]+=0.01;

/* set valid data flags */
if (opt->mode>PMODE_DGPS) {
// Only valid for AR if there is phase data.
if (!code) rtk->ssat[sat[i]-1].vsat[frq]=rtk->ssat[sat[j]-1].vsat[frq]=1;
}
else {
rtk->ssat[sat[i]-1].vsat[frq]=rtk->ssat[sat[j]-1].vsat[frq]=1;
}

#ifdef TRACE
double icb;
if (rtk->opt.glomodear==GLO_ARMODE_AUTOCAL)
icb=x[IL(frq,opt)];
else
icb=rtk->ssat[sat[i]-1].icbias[frq]*CLIGHT/freqi -
rtk->ssat[sat[j]-1].icbias[frq]*CLIGHT/freqj;
double xjj = 0.0, Pjj = 0.0;
if (opt->mode>PMODE_DGPS) {
int jj = IB(sat[j], frq, &rtk->opt);
xjj = x[jj];
Pjj = P[jj + jj * rtk->nx];
}
trace(3,"sat=%3d-%3d %s%d v=%13.3f R=%9.6f %9.6f icb=%9.3f lock=%5d x=%9.3f P=%.3f\n",
sat[i],sat[j],code?"P":"L",frq+1,v[nv],Ri[nv],Rj[nv],icb,
rtk->ssat[sat[j]-1].lock[frq],xjj,Pjj);
#endif

vflg[nv++]=(sat[i]<<16)|(sat[j]<<8)|((code?1:0)<<4)|(frq);
nb[b]++;
}
b++;
}
} /* end of system loop */

/* baseline length constraint, for fixed distance between base and rover */
if (rtk->opt.baseline[0]>0.0&&constbl(rtk,x,P,v,H,Ri,Rj,nv)) {
vflg[nv++]=3<<4;
nb[b++]++;
}
if (H) {trace(5,"H=\n"); tracemat(5,H,rtk->nx,nv,7,4);}

/* double-differenced measurement error covariance */
ddcov(nb,b,Ri,Rj,nv,R);

free(Ri); free(Rj); free(im);
free(tropu); free(tropr); free(dtdxu); free(dtdxr);

return nv;
}