（四）GNSS定位中的距离测量（模糊度固定LAMBDA算法）

定位的手段源自距离交会，因此首先要搞懂如何使用测距码，载波相位的观测值测定卫星至接收机的距离。为获得完整的载波相位观测值，还需要进行周跳的探测，整周模糊度的确定等相关问题。这里根据我自己的理解简单的叙述一下该部分内容。

测距码

测距码是一串二进制码序列。卫星发射某一结构的测距码，接收机复制出相同的测距码，当接收机接收到卫星的测距码后进行信号比对，调整延迟时间（搜索卫星信号阶段）使得信号完全对齐（两组信号的相关系数为1）。求相关系数采用积分实现，积分间隔T内将信号分割很多小份，复制码和测距码相乘，只有完全对其，积分和才等于T（-1 × -1 = 1，1 × 1 = 1,若没对齐会出现1 × -1的情况）。而当相关系数 = T 时，称为卫星信号已锁定。延迟时间即卫星信号的传播时间。伪距即速度乘时间。上述工作在接收机内在极短的时间内完成，因此可以认为伪距观测值是瞬时观测值。

使用测距码的原因：

1. 易于提取微弱的卫星信号。卫星发射功率只有20w，传到地面信号强度甚至只有噪声的万分之一，只有依赖测距码的结构，才能将其从噪声中分离出来。
2. 提高精度。由于存在误差（钟误差，传播中信号畸变），因此R不可能严格等于1，实际中采用RMax时刻进行伪距测量。因此本质是两个码在对齐的最好的状态下进行伪距测量，这样获得的伪距观测值可以认为用积分中每个码分别测距，得到的均值。精度一般可达码元宽度的1/100左右。
3. 利用码分多址对信号进行识别处理。采用码分多址技术所有信号都是相同的频率，接收机采用全向型天线，信号经过放大器进入接受通道，如指定第一通道观测3号卫星，该通道只产生3号卫星的复制码，而其他卫星的码是相互正交的，噪声因为存在随机性相关性系数也趋为0。只有3号卫星的测距码可以起作用。视场中存在n颗卫星，接收机只需要拿出n个通道进行处理，就能同时对n颗卫星进行伪距测量
4. 便于对系统进行控制和处理。由上可知，对码结构的公开和保密便可以控制用户对卫星系统的使用。GPS设计初期，也有方案是采用10种频率发射信号，每种相差一个数量级，分别测定距离中的某一位数。但毕竟GPS不是民用系统，而是军用系统，该方案不可能通过。

伪距观测方程

在伪距量测中，观测值P是信号到达接收机的时间和离开卫星时间之差乘以光速获得：$$P=c\cdot (t_R-t^S)$$由于接收机和卫星的钟均有误差，所以由卫星钟和接收机钟给出的时间和真实到达的时间有如下关系：$t^S+V_{t^S}={\tau }_a$，$t_R+V_{t_R}={\tau }_b$ 代入得：
$$P=c\cdot ({\tau }_a-{\tau }_b)+c{V}_{t^S}-c{V}_{t_R}$$但由于信号穿过电离层和对流层时的速度不是光速，因此需要加上电离层延迟改正和对流层延迟改正，才能得到真正的几何距离：
$$\rho =c\cdot ({\tau }_a-{\tau }_b)+I+T$$代入得：$P=\rho -I-T+c{V}_{t^S}-c{V}_{t_R}$
则卫星在观测瞬间在空间位置和接收机在的空间位置的几何距离为：$${{\rho }_i}^{\prime }=\sqrt{{({X^i}^{\prime }-X)}^2+{({Y^i}^{\prime }-Y)}^2+{({Z^i}^{\prime }-Z)}^2}$$由于卫星真正的空间位置和卫星星历给出的位置有误差，误差在传播路径上的投影为$\delta {\rho }^{\prime }$，多路径效应为$\delta {\rho }_{mul}$，测量噪声为${\epsilon }_i$ :
$$P=\sqrt{{(X^i-X)}^2+{(Y^i-Y)}^2+{(Z^i-Z)}^2}+c{V}_{t^S}-c{V}_{t_R}-I-T+\delta {\rho }^{\prime }+\delta {\rho }_{mul}+{\epsilon }_i$$卫星星历误差对测距的影响时未知的，通常采用足够精密的星历，使其影响可以忽略不计，多路径效应也无法模型化，只能选择接收机和观测方式进行优化，因此也忽略不计，测量噪声通常当作随机误差，无法计算具体数值，因此在伪距计算的时候，观测方程为：$$P=\sqrt{{(X^i-X)}^2+{(Y^i-Y)}^2+{(Z^i-Z)}^2}+c{V}_{t^S}-c{V}_{t_R}-I-T$$

载波相位

如果把载波当作测距信号使用，对载波进行相位测量，理论上可以获得很高的精度。但载波是没有任何标记的余弦波，因此只能测定不足一周的部分，产生整周数不确定的问题。整周计数也可能产生跳变等问题，因此数据处理过程较为复杂，也是高精度定位的代价。

重建载波

由于相位上调制了测距码和导航电文，因此需要解调，恢复载波。

方法：

1. 码相关法，使用复制码对卫星信号进行二进制相位调制，即可去掉测距码，留下载波和导航电文，而两者频率悬殊，容易用滤波器分离开。缺点：必须知道测距码的结构，若未公开则无法使用。
2. 平方法，将接受到的卫星信号自乘，得到的是载波的二次谐波，波长为原来的一半。缺点：去掉了测距码和导航电文信息，同时信噪比很差。
   GPS现代化后，已知结构的测距码越来越多，用码相关法即可得到多个频率的伪距观测值。

测量原理

某一瞬间，卫星信号在接收机R处相位为${\phi }_R$，在卫星S处相位为${\phi }_S$，以周为单位，则卫地距为：$$\rho =\lambda ({\phi }_S-{\phi }_R)$$式中，$\lambda$为载波的波长，其中包含整周计数。但卫星端并不量测载波相位，所以需要接收机产生一组与卫星载波频率和初相完全相同的基准信号（复制载波），这样相位之差即接收机产生的基准信号相位与接收到的来自卫星的载波相位之差。如果能求得完整的相位差，便可求出卫星至接收机的精确距离。

首次量测值为相位不足一周的部分，随着卫星的运动，相位差也会不断的变化，接收机可以记录下信号变化，当360°到0°，计数加一，因此量测值除了不足一周的部分，还有整周计数，因此完整的载波相位观测值由三部分给出:$$\tilde{\varphi }=\tilde{\phi }+N=\mathrm{int}(\phi )+Fr(\phi )+N$$其中N是位置的，需要通过其他途径求得。

观测方程

与伪距观测方程类似，可得：$$\lambda \tilde{\phi }=\sqrt{{(X^i-X)}^2+{(Y^i-Y)}^2+{(Z^i-Z)}^2}+c{V}_{t^S}-c{V}_{t_R}-\lambda N-I-T$$（式中忽略了星历误差，多路径效应和噪声）
与伪距观测方程相比，其中增加了一个新的未知数，整周模糊度N，同时如果广播星历求得的钟差改正数的误差可超过1m，超过载波相位测量的误差，因此需要采用精密钟差产品内插求得，或在导航电文给出的钟差基础上增加新的钟差参数，进行平差计算。或求差，必要时对对流层改正也可做类似处理。

组合观测值: 一般单差是做卫星间的观测值差，双差是在星间单差的基础上在测站间做差，三差是在双差的基础上再在历元间做差。常用的是双差，三差解几何强度很低。求差是在同一频率同一类型的观测值间做线性组合，而组合观测值则是在同一类型不同频率或者不同类型的双频观测值进行线性组合。常用的线性组合有宽项、窄项、无电离层组合、MW组合。注意组合观测值一定要保证整周模糊度的整数特性。

周跳的探测与修复

完整的载波相位由$\mathrm{int}(\phi )$、$Fr(\phi )$、$N$三部分组成。其中$Fr(\phi )$可以以极高的精度测定，但如果在观测时间计数器中置了累加，会导致$\mathrm{int}(\phi )$产生错误，这种现象称为整周跳变。导致周跳的原因有很多，卫星被遮挡、外界干扰导致信号失锁。如果求出丢失的整周数则可以恢复观测值，但如果长时间的中断或者接收机重启，信号丢失了连续性，整周模糊度则需要重新估计。

方法：

1. 高次差法，由于卫星轨迹是连续的，周跳会破坏这种连续性，整周计数每秒可变化数千周。对相邻观测值做四次差后，基本为0，如果有明显偏大的数，则发生周跳。
2. 多项式拟合法，本质上可高次差法一样，但更适合计算机运算。
3. 星间单差观测值探测修复周跳，星间单差削去了接收机钟差，对较小的周跳可以发现，但如果不同卫星同时发生周跳，也无法识别，因此探测周跳往往通过多种方法结合。
4. 双频P码伪距观测值探测修复周跳，
5. 三差观测探测修复周跳，三差消除了大部分的误差，探测出来的周跳为双差虚拟观测值的周跳，可以在双差观测值计算的时候直接使用。

整周模糊度的确定

正确快速的确定整周模糊度是载波相位测量中关键的一步。理论上在考虑接收机和卫星钟差的情况下，各种组合观测值的整周模糊度应为整数，但由于噪声和残余误差的存在，导致模糊度参数不为整数，以非整数解为基础确定模糊度的真值，对定位非常的关键。实际上，定位所需的时间相当于整周模糊度确定的时间。
解决模糊度问题一般有三种方法：
（1）用伪距观测值确定
用测距码测定的卫地距来推算载波相位中的整周模糊度，这是一种在观测值域中解决问题的方法，需要较精确的伪距观测值和电离层延迟修正方法。
从公式看，如果伪距没有误差，那么伪距除以波长减去接收机跟踪的相位部分即是整周模糊度，但是由于伪距观测值精度较低，一个历元很难准确确定，对于不发生周跳的情况，整周模糊度是不变的。可以通过多个历元平均值求N。采用这种方法对于伪距和相位测量误差相同的量显然不会有影响，但对于不同误差（如电离层）则会影响模糊度的正确确定，必须消除。CODE提高的电离层格网产品存在2-8个TECU的平均误差，局部区域甚至达到15 TECU（1TECU≈0.16m），因此精度不足，仍需要采用双频改正模型
$$N_{IF}=\frac{f_1^2}{f_1^2-f_2^2}{N}_{WL}-\frac{f_1^2}{f_1^2-f_2^2}{N}_{NL}$$右边分别为宽项和窄项的模糊度。

（2）用精确的卫星星历和先验站坐标确定
精密星历给出的是卫星质心在空间的位置，因此需要卫星中心偏差改正，接收机同理。相应的对相位观测值也要进行电离层对流层等改正。主要用于中短基线的相对定位和静态定位。如模糊函数法，不显式求解模糊度参数，直接通过搜索算法求出待定点坐标的最优解。对近似坐标的精度要求比较高。

（3）通过平差方法确定
①　固定模糊的基本原理和方法
初始解得到的模糊度参数N，及其中误差，可以根据概率密度分布得到置信区间，从数理统计分析，当置信度取得99.9％时，正确解不应落在区间外，当区间内只有一个整数时，则可认为其是正确的模糊度。由于区间内可能存在多个整数，因此一开始只能将部分模糊度固定为整数，将其代回法方程，增加解的强度，往往可以固定其余的一些模糊度，如此循环直至全部固定。
由于短时间内卫星位置变化不大，因此同一卫星在不同历元建立的方程误差系数和常数项几乎相同，这意味着方程几乎线性相关，即对解算提供的数据很少。在这种情况下，模糊度实数解的中误差必然很大，置信区间也很大。例如当中误差为±1周时，置信度取99.9％，置信区间为6.56周。置信区间内的6-7个整数解称为模糊度参数的备选解。所有卫星的模糊度参数备选组数量非常的大，如果将组合一一代入方程，观测值残差平方和最小的组合便是最终解。当整周模糊度都正确时，残差一般只有1cm。因此需要用搜索算法将正确解挑选出来。
②　快速静态定位中的FARA法：
静态定位需要长时间的静态观测，并且基线往往较长。快速定位中，基线向量一般在数公里内，误差相关性好，相对定位中容易消除误差，处理的好几个历元便可得到高精度结果。而快速定位的时间实际就是快速确定整周模糊度的时间。
FARA法，Fast Ambiguity Resolution Approach. 该方法本质是对备选组进行数理统计检验，把大量的不合理的组合剔除。统计检验的标准：任意两个整周模糊度之差是否在模糊度差值的置信区间内。FARA法充分利用了初始解协因数阵的非对角线元素所提供的模糊度间的互相关信息，对参数进行进一步的数理统计检验。对于双频观测还可用L1和L2的模糊度进行进一步统计检验。
用FARA法求出后还需要进行解的确认，即确定最优解进行三项统计检验，由于FARA是基于概率论，所以还需进行统计检验确认固定解的正确性。a. 整数解和初始解的基线向量一致性检验。 b. 整数解和初始解的单位权中误差的一致性检验。 c. 整数解中最小单位权中误差与次最小单位权中误差的显著性检验，也称ratio检验，通常要求ratio≥3。
③　LAMBDA法
Least-square AMBiguity Decorrelation Adjustment. 最小模糊度降相关平差法。该方法可以缩小搜索范围，加快搜索过程，是最成功的一种模糊度搜索算法。步骤如下：

• 整数变换
  经典搜索算法的主要问题是在观测时间较短时，初始解的实数模糊度参数精度很低，参数的相关性又很强。相关性高会导致一个模糊度参数变化影响其他模糊度的搜索，使得搜索算法的计算量巨大。如果能够设法降低模糊度参数之间的相关性，使得某一模糊度的变化对其他模糊度的取值的影响尽可能小，就能大大加快模糊度的搜索过程。而LAMBDA算法中不对整数模糊度参数搜索，而先对初始解的实数模糊度参数及其协因数阵进行整数变换，利用此来降相关。

• 为什么要进行降相关
  当模糊度参数强相关时，对某一模糊度做出微小的调整时，其余的模糊度也会互相响应，导致很难快速收敛。若第i个模糊度参数的备选组 $n_i$，如果直接搜索，模糊度组合的总和可以达到 $N=n_1\times n_2\times \cdots \times n_n$ 。每增加一个模糊度参数，搜索的数量成倍的增加。而相反，如果各个模糊度完全的不相关，任何模糊度的搜索都不会影响其余参数的取值，那便可以在 $[PVV]=\min$ 的条件下分别对每个参数进行搜索，这样总的搜索数便为$N=n_1+n_2+\cdots +n_n$，数量大大降低。事实上模糊度参数不可能完全不相关，但降相关可以大幅降低搜索难度。降相关也可以从几何上进行解释，如下：
  
  左图表示a1和a2具有较大相关性,当模糊度a1的整数选择范围如图所示时,由于椭圆的约束, a2的整数选择范围就变得很小,造成的结果是:模糊度a1引入大量不必要的候选整数,而模糊度a2的整数选择范围很小,有可能包含不了一个整数,造成搜索过程的失败。而右图表示模糊度去相关处理后的搜索椭圆,由于模糊度的整数约束条件,通常很难完全消除模糊度的相关性,使得椭圆主轴完全平行于网格主轴,但去相关处理之后,可以达到了近似平行,同时去相关处理之后搜索椭圆变得更园。当搜索范围确定了以后,搜索椭圆中所包含的格点数目就确定了,上面两图均包含11个整数候选格点,显然右图整周模糊度搜索效率要高于左图,同时右图所对应的搜索范围可以选择的更小,而左图不太可能。
  由于模糊度转换矩阵Z的整数约束条件,双差模糊度去相关处理,需要由多次整数Z变换构成,同时也很难彻底消除双差模糊度相关性,只能降低双差模糊度相关性。从几何学原理可以证明,去相关处理并没有改变搜索空间的体积,基于这个结果,转换后的搜索空间变得更加圆,模糊度变得更加不相关。

• 搜索算法
  即求新参数z的整数最小二乘解，首先要采用搜索算法从备选组中挑选出来。由于变换后的新参数的方差及参数间的互相关性都大大减小，搜索更为迅速。求得最优整数组合z后，再进行逆变换，逆变换的结果即整数模糊度向量。

  对于模糊度搜索的多维椭球区域,其搜索空间的大小完全取决于范围的选择,为了获得合理的范围值,存在多种选择的办法和手段,这里本文仅讨论在实际应用中经常考虑的两种方式。第一种方式是确保搜索空间内至少存在一个整数矢量,第二种方式是确保搜索空间内至少存在两个整数矢量。后一种方式对于实际应用非常重要,采用这种方式可以提高了模糊度搜索的可靠性,这是因为,在搜索空间中不仅包含了到实数模糊度估计矢量距离最短整数矢量,并且一定还包含了到了实数模糊度估计矢量距离次最短整数矢量。

RTKLIB中的lambda函数代码

/* lambda/mlambda integer least-square estimation ------------------------------
* integer least-square estimation. reduction is performed by lambda,
* and search by mlambda. 用lambda算法降相关，mlambda算法搜索
* args   : int    n      I  number of float parameters
*          int    m      I  number of fixed solutions
*          double *a     I  float parameters (n x 1)
*          double *Q     I  covariance matrix of float parameters (n x n)
*          double *F     O  fixed solutions (n x m)
*          double *s     O  sum of squared residulas of fixed solutions (1 x m)
* return : status (0:ok,other:error)
* notes  : matrix stored by column-major order (fortran convension)
*-----------------------------------------------------------------------------*/
extern int lambda(int n, int m, const double *a, const double *Q, double *F,
                  double *s)
{
    int info;
    double *L,*D,*Z,*z,*E;
    
    if (n<=0||m<=0) return -1;
    L=zeros(n,n); D=mat(n,1); Z=eye(n); z=mat(n,1); E=mat(n,m);
    
    /* LD factorization 对双差模糊度方差-协方差阵Q进行LDLT分解 */
    if (!(info=LD(n,Q,L,D))) {
        
        /* lambda reduction lambda降相关 */
        reduction(n,L,D,Z);
        matmul("TN",n,1,n,1.0,Z,a,0.0,z); /* z=Z'*a */
        
        /* mlambda search mlambda搜索 */
        if (!(info=search(n,m,L,D,z,E,s))) {

            /* 将在新空间中固定的模糊度逆变换回双差模糊度空间中 */
            info=solve("T",Z,E,n,m,F); /* F=Z'\E */
        }
    }
    free(L); free(D); free(Z); free(z); free(E);
    return info;
}


/* lambda reduction (z=Z'*a, Qz=Z'*Q*Z=L'*diag(D)*L) (ref.[1]) ---------------*/
static void reduction(int n, double *L, double *D, double *Z)
{
    int i,j,k;
    double del;
    
    j=n-2; k=n-2;
    //这里的调序变换类似插入排序的思路？
    while (j>=0) {
        //由于第k+1，k+2，...，n-2列都进行过降相关并且没有被上一次调序变换影响，
        //因此只需对第0,1，...，k-1，k列进行降相关
        if (j<=k) for (i=j+1;i<n;i++) gauss(n,L,Z,i,j);//从最后一列开始，各列非对角线元素从上往下依次降相关
        del=D[j]+L[j+1+j*n]*L[j+1+j*n]*D[j+1];
        if (del+1E-6<D[j+1]) { /* 检验条件，若不满足检验条件则开始进行调序变换 compared considering numerical error */
            perm(n,L,D,j,del,Z);//调序变换
            k=j; j=n-2;//完成调序变换后重新从最后一列开始进行降相关及排序，k记录最后一次进行过调序变换的列序号
        }
        else j--;
    }
}


static int search(int n, int m, const double *L, const double *D,
                  const double *zs, double *zn, double *s)
{
    int i,j,k,c,nn=0,imax=0;
    double newdist,maxdist=1E99,y;//maxdist，当前超椭圆半径
    double *S=zeros(n,n),*dist=mat(n,1),*zb=mat(n,1),*z=mat(n,1),*step=mat(n,1);
    
    k=n-1; dist[k]=0.0; //k表示当前层，从最后一层（n-1）开始计算
    zb[k]=zs[k];//即zn
    z[k]=ROUND(zb[k]); y=zb[k]-z[k]; step[k]=SGN(y);    //四舍五入取整；取整后的数与未取整的数作差；step记录z[k]是四舍还是五入
    for (c=0;c<LOOPMAX;c++) {
        newdist=dist[k]+y*y/D[k];
        if (newdist<maxdist) {  //如果当前累积目标函数计算值小于当前超椭圆半径
            if (k!=0) { //情况1：若还未计算至第一层，继续计算累积目标函数值
                dist[--k]=newdist;  //记录下当前层的累积目标函数值，dist[k]表示了第k,k+1,...,n-1层的目标函数计算和
                for (i=0;i<=k;i++)
                    S[k+i*n]=S[k+1+i*n]+(z[k+1]-zb[k+1])*L[k+1+i*n];//
                zb[k]=zs[k]+S[k+k*n];   //计算Zk，即第k个整数模糊度参数的备选组的中心
                z[k]=ROUND(zb[k]); y=zb[k]-z[k]; step[k]=SGN(y);    //四舍五入取整；取整后的数与未取整的数作差；记录是四舍还是五入
            }
            else {  //情况2：若已经计算至第一层，意味着所有层的累积目标函数值计算完毕
                //nn为当前候选解数，m为我们需要的固定解数，这里为2，表示需要一个最优解及一个次优解
                //s记录候选解的目标函数值，imax记录之前候选解中的最大目标函数值的坐标
                if (nn<m) { //若候选解数还没满
                    if (nn==0||newdist>s[imax]) imax=nn;//若当前解的目标函数值比之前最大的目标函数值都大，那么更新imax使s[imax]指向当前解中具有的最大目标函数值
                    for (i=0;i<n;i++) zn[i+nn*n]=z[i];//zn存放所有候选解
                    s[nn++]=newdist;//s记录当前目标函数值newdist，并加加当前候选解数nn
                }
                else {  //若候选解数已满（即当前zn中已经存了2个候选解）
                    if (newdist<s[imax]) {  //若 当前解的目标函数值 比 s中的最大目标函数值 小
                        for (i=0;i<n;i++) zn[i+imax*n]=z[i];    //用当前解替换zn中具有较大目标函数值的解
                        s[imax]=newdist;    //用当前解的目标函数值替换s中的最大目标函数值
                        for (i=imax=0;i<m;i++) if (s[imax]<s[i]) imax=i;    //更新imax保证imax始终指向s中的最大目标函数值
                    }
                    maxdist=s[imax];    //用当前最大的目标函数值更新超椭圆半径
                }
                z[0]+=step[0]; y=zb[0]-z[0]; step[0]=-step[0]-SGN(step[0]); //在第一层，取下一个有效的整数模糊度参数进行计算（若zb为5.3，则z取值顺序为5,6,4,7，...）
            }
        }
        else {  //情况3：如果当前累积目标函数计算值大于当前超椭圆半径
            if (k==n-1) break; //如果当前层为第n-1层，意味着后续目标函数各项的计算都会超出超椭圆半径，因此终止搜索
            else {  //若当前层不是第n-1层
                k++;    //退后一层，即从第k层退到第k+1层
                z[k]+=step[k]; y=zb[k]-z[k]; step[k]=-step[k]-SGN(step[k]); //计算退后一层后，当前层的下一个有效备选解
            }
        }
    }

    // 对s中的目标函数值及zn中的候选解进行排序（以s中目标函数值为排序标准，进行升序排序）
    // RTKLIB中最终可以得到一个最优解一个次优解，存在zn中，两解对应的目标函数值，存在s中
    for (i=0;i<m-1;i++) { /* sort by s */
        for (j=i+1;j<m;j++) {
            if (s[i]<s[j]) continue;
            SWAP(s[i],s[j]);
            for (k=0;k<n;k++) SWAP(zn[k+i*n],zn[k+j*n]);
        }
    }
    free(S); free(dist); free(zb); free(z); free(step);
    
    if (c>=LOOPMAX) {
        fprintf(stderr,"%s : search loop count overflow\n",__FILE__);
        return -1;
    }
    return 0;
}

其他方法
如Go and stop法，在初始化过程确定整周模糊度，在接收机迁站的过程中保持传递这些模糊度，到达新测站后，无需重新确定。初始化方法有：
（1）已知基线法，双差观测值解算基线向量方程中包含基线向量喝整周模糊度两类参数，若已有基线向量数据，将接收机假设到基线向量两端，可以快速确定模糊度。
（2）交换天线法，本质是求出基线向量。整周模糊度只和接受机天线和卫星有关。