Chapter.12 GNSS定位算法——模糊度固定LAMBDA算法
作者：齐花Guyc(CAUC)
文章目录
- Chapter.12 GNSS定位算法——模糊度固定LAMBDA算法
- 一.整周模糊度理论
- 1.LAMBDA算法干了一件什么事情？
- 2.LAMBDA算法步骤
- （1）去相关（Z变换）
- （2）寻找整数解（整数估计）
- （3）验证结果（比率验证）
一.整周模糊度理论
整周模糊度是全球导航卫星系统高精度定位中的核心概念，指的是载波相位观测值中未知的整周周期数。GNSS载波相位测量以波长为单位，但由于初始相位未知，观测值包含一个整数倍的波长和一个小数相位。准确解算整周模糊度并将其固定为整数是实现厘米级定位精度的关键。
在LAMBDA算法中，整周模糊度的解算基于双差载波相位数据，通过最小二乘方法和整数约束来求解。
1.LAMBDA算法干了一件什么事情？
高精度定位需要用到GNSS的载波相位数据，但这些数据里有个未知数——整周模糊度。这个模糊度是载波信号的整周周期数，类似于信号跑了几圈才到你这里，但它一开始是未知的。
LAMBDA算法就像一个“解谜大师”，专门用来破解这个模糊度，让它从一个模糊的“小数”（浮点解）变成一个确定的“整数”（整数解）。一旦模糊度被正确固定，定位精度就能从几米提升到厘米级。
2.LAMBDA算法步骤
LAMBDA算法主要分三个步骤：
- 去相关
- 寻找整数解
- 验证结果
（1）去相关（Z变换）
GNSS的模糊度数据（浮点解）就像一堆缠在一起的毛线，彼此高度相关，很难直接处理。Z变换把这些毛线解开，变成一堆独立的线（去相关后的模糊度）。这就像把一个复杂的迷宫变成规则的方格，搜索起来更省力。
setp1：获取原始数据
假设浮点解为 [1.23,4.56,−2.78]T[1.23, 4.56, -2.78]^T[1.23,4.56,−2.78]T，协方差矩阵为
Qa^a^=[1.00.80.30.81.50.40.30.41.2]Q_{\hat{a}\hat{a}} = \begin{bmatrix} 1.0 & 0.8 & 0.3 \\ 0.8 & 1.5 & 0.4 \\ 0.3 & 0.4 & 1.2 \end{bmatrix}Qa^a^=1.00.80.30.81.50.40.30.41.2
此时手上拿到一堆缠在一起的毛线和一张说明书，说明每根毛线有多粗（误差大小）和怎么缠在一起（相关性）。
setp2：分析误差表格
为了造出Z矩阵，先要通过 LTDLL^T D LLTDL 分解把误差表格拆开，弄清楚模糊度之间的缠法。
Qa^a^=LTDLQ_{\hat{a}\hat{a}} = L^T D LQa^a^=LTDL
比如分解后为：
D=[0.50000.70000.6],L=[100l2110l31l321]D = \begin{bmatrix} 0.5 & 0 & 0 \\ 0 & 0.7 & 0 \\ 0 & 0 & 0.6 \end{bmatrix}, \quad L = \begin{bmatrix} 1 & 0 & 0 \\ l_{21} & 1 & 0 \\ l_{31} & l_{32} & 1 \end{bmatrix}D=0.50000.70000.6,L=1l21l3101l32001
此时相当于把毛线团拆开，记录下每根线的长度（D）和怎么交叉缠绕（L），为打造整理机（Z矩阵）做准备。
setp3：构造Z矩阵
使用整数高斯变换构造Z矩阵，通过一系列整数运算，逐步构造Z矩阵，目标是让 Qz^z^Q_{\hat{z}\hat{z}}Qz^z^ 的非对角元素尽量小。
Z变换把高相关的模糊度变成低相关的，搜索空间从细长椭圆变成接近圆形。
整数高斯变换步骤
- 初始化Z矩阵
一开始，Z矩阵设为单位矩阵。Z=[100010001]Z = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}Z=100010001 - 逐对降低相关性
检查矩阵Qa^a^Q_{\hat{a}\hat{a}}Qa^a^的非对角元素，找到相关性最大的元素（Qa^a^(1,2)=0.8Q_{\hat{a}\hat{a}}(1,2) = 0.8Qa^a^(1,2)=0.8）。用整数高斯变换调整Z矩阵，减少这对模糊度的相关性。方法是构造一个整数变换矩阵，比如：
T=[1−k0010001]T = \begin{bmatrix} 1 & -k & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}T=100−k10001
其中 kkk 是整数（可以通过四舍五入 Qa^a^(1,2)/Qa^a^(2,2)Q_{\hat{a}\hat{a}}(1,2)/Q_{\hat{a}\hat{a}}(2,2)Qa^a^(1,2)/Qa^a^(2,2) 得到）。
更新Z矩阵：Z=Z⋅TZ = Z \cdot TZ=Z⋅T，同时更新逆矩阵：Z−1=T−1⋅Z−1Z^{-1} = T^{-1} \cdot Z^{-1}Z−1=T−1⋅Z−1。 - 迭代优化
重复检查非对角元素，找下一个相关性最大的元素，继续重复步骤2。
每次调整后，计算新的协方差矩阵。
Qz^z^=Z−1Qa^a^Z−TQ_{\hat{z}\hat{z}} = Z^{-1} Q_{\hat{a}\hat{a}} Z^{-T}Qz^z^=Z−1Qa^a^Z−T
直到非对角元素足够小（接近对角矩阵）为止。
setp4：变换数据
用Z矩阵变换模糊度和协方差矩阵
用Z矩阵的逆把原始模糊度变成新模糊度：z^=Z−1a^\hat{z} = Z^{-1} \hat{a}z^=Z−1a^
用Z矩阵把原始协方差变成新协方差：Qz^z^=Z−1Qa^a^Z−TQ_{\hat{z}\hat{z}} = Z^{-1} Q_{\hat{a}\hat{a}} Z^{-T}Qz^z^=Z−1Qa^a^Z−T
Qz^z^=[0.50.10.050.10.70.030.050.030.6]Q_{\hat{z}\hat{z}} = \begin{bmatrix} 0.5 & 0.1 & 0.05 \\ 0.1 & 0.7 & 0.03 \\ 0.05 & 0.03 & 0.6 \end{bmatrix}Qz^z^=0.50.10.050.10.70.030.050.030.6
非对角元素的值变小，说明相关性降低。
（2）寻找整数解（整数估计）
在变换后的模糊度中，寻找最优整数解，用Z反变换回到原始模糊度的整数解。
整数估计的方法有三种：
- 四舍五入
- 整数自举法
- 整数最小二乘（ILS）
四舍五入
没考虑模糊度间的关系，成功率低；随便猜，简单但容易错。
整数自举法
从最精确的模糊度开始取整，根据相关性调整其他模糊度，再取整；考虑部分相关性，成功率比四舍五入高。跟着线索一步步找，比瞎猜好，但可能不是最佳位置。
整数最小二乘（ILS）
寻找让误差最小的整数解，数学上：
a^=argmina∈Zn(a^−a)TQa^a^−1(a^−a)\hat{a} = \arg \min_{a \in \mathbb{Z}^n} (\hat{a} - a)^T Q_{\hat{a}\hat{a}}^{-1} (\hat{a} - a)a^=arga∈Znmin(a^−a)TQa^a^−1(a^−a)
a^\hat{a}a^ 是原始模糊度(浮点解)；aaa 是可能的整数解;Qa^a^Q_{\hat{a}\hat{a}}Qa^a^ 是协方差矩阵，记录模糊度的可靠性和相关性。(a^−a)TQa^a^−1(a^−a)(\hat{a} - a)^T Q_{\hat{a}\hat{a}}^{-1} (\hat{a} - a)(a^−a)TQa^a^−1(a^−a)：衡量整数解离浮点解有多远，考虑了数据的可靠性和相关性。
步骤：
- 定义搜索空间
在去相关空间中，找到误差最小的整数解，目标函数是：
F(z)=(z^−z)TQz^z^−1(z^−z)≤χ2F(z) = (\hat{z} - z)^T Q_{\hat{z}\hat{z}}^{-1} (\hat{z} - z) \leq \chi^2F(z)=(z^−z)TQz^z^−1(z^−z)≤χ2
定义了一个超椭球（搜索空间），包含所有足够接近 z^\hat{z}z^ 的整数点 zzz；χ2\chi^2χ2 是阈值，控制椭球大小。 - 搜索整数解
在超椭球内，找到让 F(z)F(z)F(z) 最小的整数解（最优解），并记录次优解。
有两种搜索方式：枚举搜索和收缩搜索
枚举搜索检查椭球内所有整数点，计算每个点的 F(z)F(z)F(z)，找出最小和次小值。全面不会漏解，但耗时长。
收缩搜索动态缩小椭球，优先检查最可能的整数点。
初始化搜索空间后，分解Qz^z^Q_{\hat{z}\hat{z}}Qz^z^，用Cholesky分解将其分解为Qz^z^=LLTQ_{\hat{z}\hat{z}} = L L^TQz^z^=LLT，令y=L−1(z^−z)y = L^{-1} (\hat{z} - z)y=L−1(z^−z)，则变换目标函数为F(z)=(z^−z)TQz^z^−1(z^−z)=yTy=∑i=1nyi2F(z) = (\hat{z} - z)^T Q_{\hat{z}\hat{z}}^{-1} (\hat{z} - z) = y^T y = \sum_{i=1}^n y_i^2F(z)=(z^−z)TQz^z^−1(z^−z)=yTy=i=1∑nyi2
椭球变成：∑i=1nyi2≤χ2\sum_{i=1}^n y_i^2 \leq \chi^2∑i=1nyi2≤χ2 更简单的球形约束。
假设前 i−1i-1i−1 个整数 z1,z2,…,zi−1z_1, z_2, \dots, z_{i-1}z1,z2,…,zi−1 已确定，计算第 iii 个整数 ziz_izi 的条件范围：
zi∈[⌊z^i∣i−1−χ2−∑j=1i−1yj2lii⌉,⌈z^i∣i−1+χ2−∑j=1i−1yj2lii⌉]z_i \in \left[ \left\lfloor \hat{z}_{i|i-1} - \frac{\sqrt{\chi^2 - \sum_{j=1}^{i-1} y_j^2}}{l_{ii}} \right\rceil, \left\lceil \hat{z}_{i|i-1} + \frac{\sqrt{\chi^2 - \sum_{j=1}^{i-1} y_j^2}}{l_{ii}} \right\rceil \right]zi∈z^i∣i−1−liiχ2−∑j=1i−1yj2,z^i∣i−1+liiχ2−∑j=1i−1yj2
z^i∣i−1\hat{z}_{i|i-1}z^i∣i−1：基于前 i−1i-1i−1 个整数的条件估计；liil_{ii}lii：LLL 的对角元素。
从初始椭球开始，逐维度搜索整数点，找到更小误差的解，更新 χ2\chi^2χ2，缩小椭球。
复习一下卡方分布χ2\chi^2χ2的概念：
想象在玩投飞镖游戏，每次投的飞镖偏离靶心的距离是一个正态分布（随机但围绕中心）。如果你投了很多次，把每次偏离距离的平方加起来，这个总和的表现规律就是卡方分布。告诉你，误差平方的总和大概会落在什么范围内。
如果有 kkk 个独立的标准正态分布变量（均值为0，方差为1），记为 Z1,Z2,…,ZkZ_1, Z_2, \dots, Z_kZ1,Z2,…,Zk，那么它们的平方和：
X=Z12+Z22+⋯+Zk2X = Z_1^2 + Z_2^2 + \dots + Z_k^2X=Z12+Z22+⋯+Zk2
服从卡方分布，自由度为 kkk，记为 X∼χ2(k)X \sim \chi^2(k)X∼χ2(k)，自由度 kkk 表示有多少个独立变量，决定了卡方分布的形状。
在此参考：张镇虎的博客
自由度 kkk：表示独立正态变量的个数；显著性水平 α\alphaα：表示误差平方和超过 χ2\chi^2χ2 的概率，置信水平 1−α1-\alpha1−α 表示包含在 χ2\chi^2χ2 内的概率。
例子：自由度 k=3k = 3k=3，置信水平95%（α=0.05\alpha = 0.05α=0.05），查表得 χ2=7.815\chi^2 = 7.815χ2=7.815，表示95%的概率误差平方和小于7.815。
（3）验证结果（比率验证）
用比率检验来检查ILS找到的整数模糊度解是否可靠，确保最佳整数解是正确的，避免选错导致定位误差。ILS会输出一个最优解和一个次优解，比率检验通过比较这两者的误差大小，判断最佳解是否明显优于其他解，从而决定是否接受。
比较这两个整数解离浮点解有多远。如果最优解明显比次优解近很多，就放心取；如果两者差不多近，就需要保持怀疑，可能选择保留浮点解。
比率=F(a^)F(a^′)≤μ\text{比率} = \frac{F(\hat{a})}{F(\hat{a}')} \leq \mu比率=F(a^′)F(a^)≤μ
F(a^)=(a^−a^)TQa^a^−1(a^−a^)F(\hat{a}) = (\hat{a} - \hat{a})^T Q_{\hat{a}\hat{a}}^{-1} (\hat{a} - \hat{a})F(a^)=(a^−a^)TQa^a^−1(a^−a^): 最优解的加权平方误差；F(a^′)=(a^−a^′)TQa^a^−1(a^−a^′)F(\hat{a}') = (\hat{a} - \hat{a}')^T Q_{\hat{a}\hat{a}}^{-1} (\hat{a} - \hat{a}')F(a^′)=(a^−a^′)TQa^a^−1(a^−a^′): 次优解的加权平方误差。
Qa^a^Q_{\hat{a}\hat{a}}Qa^a^: 原始模糊度的协方差矩阵；
μ\muμ: 阈值，决定接受最佳解的严格程度，μ\muμ的取值可以根据设置可接受的失败率 PfP_fPf，比如0.001（0.1%）计算，也可以固定阈值（经验值）。
如果比率 ≤μ\leq \mu≤μ，接受最佳解 a^\hat{a}a^; 否则拒绝，保留浮点解 a^\hat{a}a^。
)
