一种对流层长距离亚厘米级实时动态卫星导航定位方法

摘要

本发明涉及一种对流层长距离亚厘米级实时动态卫星导航定位方法，该方法包括：首先，固定L1和L2载波的双差模糊度，固定模糊度之后，建立含对流层参数的电离层加权双差观测方程；其次，确定降相关参数，在计算降相关参数过程中，需要用到未知数的真值，采用不含对流层参数的最小二乘位置参数的协因数矩阵代替求解降相关参数过程中的参数真值的二次型矩阵；最后，根据降相关参数求得最优降相关解，从而得到可靠的亚厘米级坐标估计值。与现有技术相比，本发明具有精度高、稳定性强、适用范围广等优点。

说明书

技术领域

本发明涉及用于卫星导航的定位技术，特别是涉及一种对流层长距离亚厘米级实时动态卫星导航定位方法。 

背景技术

在过去的二十年，GNSS RTK精密定位技术有了长足地发展，并被广泛应用于很多领域，如测地学、工程学以及仪器自动化等。在地学应用中，高频GNSS数据可以通过一系列的频率和振幅捕捉到同震位移，从某种角度讲，这种测量的广度要比地震传感器更广。它能够检测来自几百公里甚至几千公里外的任意大小的地面移动和同震平面波。另一方面，也可以利用GNSS RTK确定当地大地水准面模型。然而，制约现有RTK技术在地学研究中应用的主要因素有两个：(1)现在RTK技术可以达到厘米级的水平，但在高度方向上的精度要比在水平方向上差1.5至2倍。(2)RTK技术的定位精度随着基准站与流动站的距离的增加迅速衰减。 

在载波相位模糊度解算和点位精度计算中，电离层和对流层的影响是误差的最主要来源。如果模糊度能正确固定，利用双频GNSS信号，采用无电离层组合，可以将电离层对点位计算的影响基本消除。然而，即使利用多频信号，对流层延迟的影响也无法削弱。且对流层延迟对点位误差的影响与基准站至流动站的距离有关，通常数十公里有几厘米的误差。在未来的GNSS系统中，多频GNSS信号将可以民用。对于数百公里的距离，所有载波相位的整周模糊度能可靠地固定，而此时，对流层对点位的影响可以达到几分米。如果精密GNSS轨道代替了广播星历，那么，对于定位，对流层延迟将会是最重要的因素。 

为了解决对流层影响，已发展了许多方法。 

(1)对于短基线而言，最简单的方法是忽略对流层残差的影响。因为对于这种情况，假设两个接收机是相关的，那么，在双差之后，它们的影响可以消除。例如，对于传统的单基站RTK，基准站与流动站的距离限制在20km以内，是因为这 样能最大限度消除大气影响。 

(2)网络RTK则是另一种减小对流层影响的方法，这种方法可以扩大服务范围。在网络覆盖范围内，流动站与基准站的对流层延迟是通过对多基准站相关的对流层延迟内插得到。虽然网络RTK已经在商业服务中广泛应用，但是随着站间的距离增加，两个站间的对流层延迟相关性减弱，导致对流层残余误差增加。因此，现有的基于网络RTK的双频系统的有效距离在50-70公里。即使三频信号的加入，这个距离也只是在100公里-140公里。另外，在大气情况不利的条件下，对于局部的异常情况，流动站的内插改正显得捉襟见肘。 

(3)引入对流层参数，例如引入天顶对流层延迟参数来吸收对流层残差。这种方法能克服距离的限制，对于长基线，不管是单基站还是网络RTK，都基本能实现良好的应用。但是，由于对流层参数和高程参数的强相关性，一旦将对流层参数和坐标参数一起平差处理，将导致平差模型严重病态，通常需要积累足够的观测值以克服该病态问题，因此，通常采用滤波技术逐渐地引入新的观测值，序贯计算定位结果。然而，采用滤波方法求解对流层参数充其量是多个历元对流层延迟的平均值，不能很好地反映对流层的实时变化，尤其是流动站在复杂的大气环境和高速移动的情况下，所以本发明为解决这一问题提出了一种可行的方法。 

发明内容

本发明的目的就是为了克服上述现有技术存在的缺陷而提供一种精度高、稳定性强、适用范围广的对流层长距离亚厘米级实时动态卫星导航定位方法。本发明针对在模糊度正确固定的基础上，采用GNSS信号难以实现长距离亚厘米级实时定位的两个主要限制因素：即对流层延迟的有效抑制和引入对流层参数后模型的严重病态性，提出了采用降相关方法实现长距离亚厘米级实时定位算法。 

本发明的目的可以通过以下技术方案来实现： 

一种对流层长距离亚厘米级实时动态卫星导航定位方法，该方法包括以下步骤： 

(a)获取GNSS双频或多频实时观测数据； 

(b)对实时观测数进行数据预处理； 

(c)构建相位与伪距的差分观测方程，包括单差观测方程和双差观测方程； 

(d)采用最小二乘方法计算模糊度浮点解； 

(e)由于双差模糊度应该是一个整数，而采用最小二乘解得的模糊度为浮点解，所以采用LAMBDA方法对模糊度浮点解进行固定，得到固定后的双差模糊度，并采用Ps-LAMBDA检验模糊度的可靠性，当固定了双差模糊度后，就成为整周模糊度； 

(f)根据步骤(d)和(e)，得到双差模糊度改正后的相位与伪距的双差观测方程； 

(g)利用等价性原理，在步骤(f)获得的双差观测方程两边乘以一个变换矩阵，消除电离层延迟，并利用最小二乘法形成法方程矩阵，获得最小二乘解； 

(h)构建不含对流层参数的误差方程，并基于最小二乘准则计算最小二乘解及其相应的协方差矩阵； 

(i)采用不含对流层天顶延迟参数的最小二乘解的协方差矩阵代替求解降相关参数过程中涉及到的参数真值的二次型矩阵，从而计算得到降相关参数； 

(j)确定降相关参数后，采用降相关方法重新求解步骤(g)的病态问题，得到可靠的亚厘米级坐标估计值。 

所述的步骤(b)中，数据预处理包括卫星截止高度角设置、时标校正、相位观测值周跳探测与修复、粗差探测与处理以及卫星和接收机的天线相位中心修正。 

优选地，采用的高度角定权公式为

所述的步骤(c)中构建的双差观测方程如下： 

$\left(\begin{array}{c}{\Phi }_1=\rho +T-I_1-{\lambda }_1{N}_1+{ϵ}_{{\Phi }_1}\\ {\Phi }_2=\rho +T-f_1^2/f_2^2\times I_1-{\lambda }_2{N}_2+{ϵ}_{{\Phi }_2}\\ P_1=\rho +T+I_1+{ϵ}_{P_1}\\ P_2=\rho +T+f_1^2/f_2^2\times I_1+{ϵ}_{P_2}\end{array}\right)---\left(1\right)$

式中，下标“1”、“2”分别表示与L1和L2载波相位有关；Φ和P分别表示双差相位和伪距观测值；ρ是双差卫星-地球距离真值；I₁表示L1频率的双差电离层延迟；T是用UNB3对流层标准模型和Niell映射函数改正后的双差对流层延迟残差矩阵，T＝b×τ，τ是对流层天顶延迟参数，b为对应系数矩阵；ε是期望为0的正太分布的随机噪声；f₁和f₂表示L1和L2上的频率；λ₁和λ₂表示L1和L2的波长；N₁和N₂表示L1和L2上的整周模糊度。 

所述的步骤(f)中获得的双差观测方程如下： 

$\left(\begin{array}{c}{\Phi }_1=\rho +T-I_1+{ϵ}_{{\Phi }_1}\\ {\Phi }_2=\rho +T-f_1^2/f_2^2\times I_1+{ϵ}_{{\Phi }_2}\\ P_1=\rho +T+I_1+{ϵ}_{P_1}\\ P_2=\rho +T+f_1^2/f_2^2\times I_1+{ϵ}_{P_2}\end{array}\right)---\left(2\right)$

将其写成矩阵形式如下： 

$\left(\begin{array}{c}{\Phi }_1\\ {\Phi }_2\\ P_1\\ P_2\end{array}\right)=\left(\begin{array}{ccc}A& b& -E_m\\ A& b& -f_1^2/f_2^2{E}_m\\ A& b& E_m\\ A& b& f_1^2/f_2^2{E}_m\end{array}\right)\left(\begin{array}{c}x\\ \tau \\ I_1\end{array}\right)+\left(\begin{array}{c}{ϵ}_{{\Phi }_1}\\ {ϵ}_{{\Phi }_2}\\ {ϵ}_{P_1}\\ {ϵ}_{P_2}\end{array}\right)---\left(3\right)$

式中，A是m×3维有关坐标[x y z]^T在坐标系下的系数矩阵，E_m是m×m单位矩阵。 

所述的步骤(g)中，在双差观测方程两边乘以一个变换矩阵后得到如下误差方程： 

y_IF＝Hξ+ε_IF  (4) 

其中，$H=\left(\begin{array}{cc}{e}_3\otimes A& e_3\otimes b\end{array}\right),\xi =\left(\begin{array}{c}x\\ \tau \end{array}\right),$e₃为3×3单位矩阵，ε_IF为随机噪声，上述误差方程的最小二乘解为： 

${\hat{\xi }}_{\mathrm{LS}}=N_{\xi }^{-1}{w}_{\xi }---\left(5\right)$

其中， 

$N_{\xi }=H^T{Q}_{\mathrm{IF}}^{-1}H=e_3^T{\left(\mathrm{T\Lambda T}\right)}^{-1}{e}_3\otimes \left(\begin{array}{cc}{A}^T{Q}_{\mathrm{DD}}^{-1}A& A^T{Q}_{\mathrm{DD}}^{-1}b\\ b^T{Q}_{\mathrm{DD}}^{-1}A& b^T{Q}_{\mathrm{DD}}^{-1}b\end{array}\right)---\left(6\right)$

$w_{\xi }=H^T{Q}_{\mathrm{IF}}^{-1}{y}_{\mathrm{IF}}=\left(\begin{array}{cc}\left(\begin{array}{c}{e}_3^T{\left(\mathrm{T\Lambda T}\right)}^{-1}\otimes A^T{Q}_{\mathrm{DD}}^{-1}\\ e_3^T{\left(\mathrm{T\Lambda T}\right)}^{-1}\otimes b^T{Q}_{\mathrm{DD}}^{-1}\end{array}\right)& Q_{\mathrm{IF}}\end{array},=(T{\mathrm{\Lambda T}}^T\otimes Q_{\mathrm{DD}})\right)---\left(7\right)$

Q_DD是单历元对应于双差观测值的矩阵。 

所述的步骤(h)中，不含对流层参数的误差方程如下： 

$\left(\begin{array}{c}{\Phi }_1\\ {\Phi }_2\\ P_1\\ P_2\end{array}\right)=\left(\begin{array}{cc}A& -E_m\\ A& -f_1^2/f_2^2{E}_m\\ A& E_m\\ A& f_1^2/f_2^2{E}_m\end{array}\right)\left(\begin{array}{c}x\\ I_1\end{array}\right)+\left(\begin{array}{c}{ϵ}_{{\Phi }_1}\\ {ϵ}_{{\Phi }_2}\\ {ϵ}_{P_1}\\ {ϵ}_{P_2}\end{array}\right)---\left(3\right)$

其相应的最小二乘解为： 

${\hat{x}}_{\mathrm{LS}}=N_x^{-1}{w}_x---\left(9\right)$

其中， 

$N_x=e_3^T{\left(\mathrm{T\Lambda T}\right)}^{-1}{e}_3\otimes A^T{Q}_{\mathrm{DD}}^{-1}A---\left(10\right)$

$w_{\xi }=[e_3^T{\left(\mathrm{T\Lambda T}\right)}^{-1}\otimes A^T{Q}_{\mathrm{DD}}^{-1}]y_{\mathrm{IF}}---\left(11\right)$

相应的协方差矩阵为： 

$\Sigma {\hat{x}}_{\mathrm{LS}}={\sigma }_0^2{N}_x^{-1}---\left(12\right)$

其中，为方差分量单位权。 

所述的步骤(i)中，降相关参数α根据以下MSE最小准则求得： 

α＝arg min trace(M_R)  (13) 

$M_R={\sigma }_0^2{N}_R^{-1}{N}_{\xi }{N}_R^{-1}+{\alpha }^2{N}_R^{-1}S{\overline{\mathrm{\xi \xi }}}^{T}S{N}_R^{-1}---\left(\text{14}\right)$

其中，S为降相关矩阵，N_R＝N_ξ+αS； 

采用步骤(h)获得协方差矩阵代替式(14)中的二次型矩阵，则 

$S{\overline{\mathrm{\xi \xi }}}^{T}S=\left(\begin{array}{cc}\Sigma {\hat{x}}_{\mathrm{LS}}& 0\\ 0& 0\end{array}\right)={\sigma }_0^2\left(\begin{array}{cc}{N}_x^{-1}& 0\\ 0& 0\end{array}\right)={\sigma }_0^2{Q}_x---\left(15\right)$

则 

$M_R={\sigma }_0^2(N_R^{-1}{N}_{\xi }{N}_R^{-1}+{\alpha }^2{N}_R^{-1}{Q}_x{N}_R^{-1})={\sigma }_0^2{Q}_R---\left(16\right)$

α＝arg min_α＞0trace(Q_R)  (17) 

根据式(17)计算获得降相关参数α。 

所述的步骤(j)中，根据步骤(i)获得的降相关参数α，计算降相关解： 

${\hat{\xi }}_R=N_R^{-1}{w}_{\xi }---\left(14\right)$

相应的协方差矩阵为： 

${\Sigma }_{\stackrel{·}{\xi }R}={\sigma }_0^2{N}_R^{-1}{N}_{\xi }{N}_R^{-1}---\left(15\right)$

以求得的降相关解作为最终的亚厘米级坐标估计值。 

与现有技术相比，本发明具有以下优点： 

(1)引入对流层参数来吸收对流层残差，克服了距离的限制。 

(2)采用电离层加权模型将电离层的影响加以处理，能有效消除电离层误差，且保持定位模型的强度。根据等价性原理，在观测方程两边乘以一个合适的变换矩阵，既能消除电离层影响而又对其余参数的估计没有任何影响，此外，还可以对电离层进行适当的加权，从而得到更好的结果。 

(3)本发明采用LAMBDA方法求解双差模糊度，速度快，且能充分利用模糊度之间的相关性，采用Ps-LAMBDA方法能有效控制模糊度固定的可靠性。 

(4)引入对流层参数后，模型病态性很强，采用降相关方法逐个历元独立计算获得稳定的亚厘米级定位解，特别是高程定位的精度达到了毫米级。 

(5)引入对流层天顶延迟参数，能有效吸收经对流层基本模型改正后的残余对流层误差。 

(6)在计算降相关参数时，采用不含对流层参数的最小二乘解的协方差阵代替参数真值的二次型矩阵。由此，降相关参数的解算与观测几何的变化有关，显著降低了对流层参数与高程之间的相关性，从而明显降低了含对流层参数的误差方程的病态问题。 

(7)由于该发明实现单历元单独计算位置解，因此，每个历元求得的对流层参数具有实时性，能够真实地表征对流层的实时变化，避免了前后历元错误观测对当前历元解的影响。 

附图说明

图1为本发明的流程示意图。 

具体实施方式

下面结合附图和具体实施例对本发明进行详细说明。本实施例以本发明技术方案为前提进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。 

一种对流层长距离亚厘米级实时动态卫星导航定位方法，该方法实施步骤为： 

首先，固定L1和L2载波的双差模糊度，固定模糊度之后，建立含对流层参数的电离层加权双差观测方程。 

其次，确定降相关参数。求得降相关参数是求降相关解的前提。降相关参数是根据求最小均方根误差(MSE)得到。而在计算降相关参数过程中，需要用到未知数的真值，通常情况下，可以采用含对流层参数的最小二乘估值来代替未知数的真值，但是，这个模型本身具有严重病态性，这个解以及它的协方差矩阵本身就不稳定，用它来代替未知数的真值显然是不合理的。本发明中采用不含对流层参数的最小二乘位置参数的协因数矩阵代替求解降相关参数过程中的参数真值的二次型矩阵。 

最后，根据降相关参数求得最优降相关解。 

本发明方法在实际操作中应注意的具体问题如下： 

需要设置参考站，参考站的设置应该尽量位于测区覆盖范围的中心部位，要求周围观测环境良好，视野开阔，视场内周围障碍物的高度角一般应小于10度；点位应远离大功率无线电发射源(如电台，微波站及微波通道等)及高压线，以避免周围磁场对信号的干扰；点位周围不应有对电磁波反射(或吸收)强烈的物体，以 减弱多路径效应的影响，最好使用扼流圈天线以有效抑制多路径效应。 

如图1所示，本发明方法具体包括以下步骤： 

在步骤a中，获取GNSS双频或多频实时观测数据。 

在步骤b中，对实时观测数进行数据预处理，包括卫星截止高度角设置、时标校正、相位观测值周跳探测与修复、粗差探测与处理、卫星和接收机的天线相位中心修正以及其他必要的改正等。其中，通过设置卫星截止高度角剔除受电离层和对流层影响较大的低高度角观测数据；通过单点定位的方法进行时标校正，实现流动站与参考站的时标统一；采用GF-MW组合方法进行周跳探测与修复。 

在步骤c中，构建相位与伪距的差分观测方程，包括单差、双差、三差观测方程，构建的双差观测方程如下： 

$\left(\begin{array}{c}{\Phi }_1=\rho +T-I_1-{\lambda }_1{N}_1+{ϵ}_{{\Phi }_1}\\ {\Phi }_2=\rho +T-f_1^2/f_2^2\times I_1-{\lambda }_2{N}_2+{ϵ}_{{\Phi }_2}\\ P_1=\rho +T+I_1+{ϵ}_{P_1}\\ P_2=\rho +T+f_1^2/f_2^2\times I_1+{ϵ}_{P_2}\end{array}\right)---\left(1\right)$

式中，下标“1”、“2”分别表示与L1和L2载波相位有关；Φ和P分别表示双差相位和伪距观测值；ρ是双差卫星-地球距离真值：I₁表示L1频率的双差电离层延迟；T是用UNB3对流层标准模型和Niell映射函数改正后的双差对流层延迟残差矩阵，T＝b×τ，τ是对流层天顶延迟参数，b为对应系数矩阵；ε是期望为0的正太分布的随机噪声；f₁和f₂表示L1和L2上的频率；λ₁和λ₂表示L1和L2的波长；N₁和N₂表示L1和L2上的整周模糊度。 

在步骤d中，采用最小二乘方法计算模糊度浮点解。 

在步骤e中，由于双差模糊度应该是一个整数，而采用最小二乘解得的模糊度为浮点解，所以采用LAMBDA方法对模糊度浮点解进行固定，得到固定后的双差模糊度，并采用Ps-LAMBDA检验模糊度的可靠性，当固定了双差模糊度后，就成为整周模糊度。 

在步骤f中，根据步骤d和e，得到双差模糊度改正后的相位与伪距的双差观测方程，如式(2)所示： 

$\left(\begin{array}{c}{\Phi }_1=\rho +T-I_1+{ϵ}_{{\Phi }_1}\\ {\Phi }_2=\rho +T-f_1^2/f_2^2\times I_1+{ϵ}_{{\Phi }_2}\\ P_1=\rho +T+I_1+{ϵ}_{P_1}\\ P_2=\rho +T+f_1^2/f_2^2\times I_1+{ϵ}_{P_2}\end{array}\right)---\left(2\right)$

将其写成矩阵形式如下： 

$\left(\begin{array}{c}{\Phi }_1\\ {\Phi }_2\\ P_1\\ P_2\end{array}\right)=\left(\begin{array}{ccc}A& b& -E_m\\ A& b& -f_1^2/f_2^2{E}_m\\ A& b& E_m\\ A& b& f_1^2/f_2^2{E}_m\end{array}\right)\left(\begin{array}{c}x\\ \tau \\ I_1\end{array}\right)+\left(\begin{array}{c}{ϵ}_{{\Phi }_1}\\ {ϵ}_{{\Phi }_2}\\ {ϵ}_{P_1}\\ {ϵ}_{P_2}\end{array}\right)---\left(3\right)$

式中，A是m×3维有关坐标[x y z]^T在坐标系下的系数矩阵，E_m是m×m单位矩阵，和是(m×1)维L1双差相位观测值和它们的噪 

声，其他的符号类似。 

可将式(2)简写为： 

$y_{\mathrm{DD}}={\left(\begin{array}{cccc}{\Phi }_1^T& {\Phi }_2^T& P_1^T& P_2^T\end{array}\right)}^T---\left(4\right)$

${ϵ}_{\mathrm{DD}}={\left(\begin{array}{cccc}{ϵ}_{{\Phi }_1}^T& {ϵ}_{{\Phi }_2}^T& {ϵ}_{P_1}^T& {ϵ}_{P_2}^T\end{array}\right)}^T---\left(5\right)$

则单历元双差观测值的随机模型为： 

$\Sigma y_{\mathrm{DD}}={\sigma }_{{\Phi }_1}^2(\Lambda \otimes Q_{\mathrm{DD}})---\left(6\right)$

其中，对角阵为$\Lambda =\mathrm{diag}\left(\begin{array}{cccc}1& s_{{\Phi }_2}& s_{P_1}& s_{P_2}\end{array}\right),s_{{\Phi }_2}={\sigma }_{{\Phi }_2}^2/{\sigma }_{{\Phi }_1}^2,s_{P_1}={\sigma }_{P_1}^2/{\sigma }_{{\Phi }_1}^2,s_{P_2}={\sigma }_{P_2}^2/{\sigma }_{{\Phi }_1}^2$是Φ₂、P₁和P₂相对于Φ₁的方差系数；和是不同观测类型的非差观测值的方差分量。Q_DD是单历元对应于双差观测值的矩阵，如果所有观测值独立时，那么它们有相同的方差，其中，e_m是m维元素都为1的列向量。在本发明中，采用作为方差分量的单位权，并对原始观测值采用高度角定权的方法，采用的高度角定权公式为

在步骤g中，利用等价性原理，在步骤(f)获得的双差观测方程两边乘以一个变换矩阵消除电离层延迟，并利用最小二乘法形成法方程矩阵，获得最小二乘解，具体如下： 

$T=\left(\begin{array}{cccc}\beta & \gamma & 0& 0\\ 0& \gamma /(\gamma -\beta )& \beta /(\beta -\gamma )& 0\\ 0& 0& \beta & \gamma \end{array}\right),\beta =f_1^2/(f_1^2-f_2^2),\gamma =-f_2^2/(f_1^2-f_2^2).$

在双差观测方程两边乘以一个变换矩阵后得到如下误差方程： 

y_IF＝Hξ+ε_IF  (7) 

其中，$H=\left(\begin{array}{cc}{e}_3\otimes A& e_3\otimes b\end{array}\right),\xi =\left(\begin{array}{c}x\\ \tau \end{array}\right),$e₃为3×3单位矩阵，ε_IF为随机噪声，上述误差方程的最小二乘解为： 

${\hat{\xi }}_{\mathrm{LS}}=N_{\xi }^{-1}{w}_{\xi }---\left(8\right)$

其中， 

$N_{\xi }=H^T{Q}_{\mathrm{IF}}^{-1}H=e_3^T{\left(\mathrm{T\Lambda T}\right)}^{-1}{e}_3\otimes \left(\begin{array}{cc}{A}^T{Q}_{\mathrm{DD}}^{-1}A& A^T{Q}_{\mathrm{DD}}^{-1}b\\ b^T{Q}_{\mathrm{DD}}^{-1}A& b^T{Q}_{\mathrm{DD}}^{-1}b\end{array}\right)---\left(9\right)$

$w_{\xi }=H^T{Q}_{\mathrm{IF}}^{-1}{y}_{\mathrm{IF}}=\left(\begin{array}{cc}\left(\begin{array}{c}{e}_3^T{\left(\mathrm{T\Lambda T}\right)}^{-1}\otimes A^T{Q}_{\mathrm{DD}}^{-1}\\ e_3^T{\left(\mathrm{T\Lambda T}\right)}^{-1}\otimes b^T{Q}_{\mathrm{DD}}^{-1}\end{array}\right)& Q_{\mathrm{IF}}\end{array},=(T{\mathrm{\Lambda T}}^T\otimes Q_{\mathrm{DD}})\right)---\left(10\right)$

最小二乘估值的协方差矩阵为 

$\Sigma {\hat{x}}_{\mathrm{LS}}={\sigma }_0^2{N}_x^{-1}---\left(11\right)$

在步骤h中，构建不含对流层参数的误差方程： 

其最小二乘解为： 

${\hat{x}}_{\mathrm{LS}}=N_x^{-1}{w}_x---\left(13\right)$

其中， 

$N_x=e_3^T{\left(\mathrm{T\Lambda T}\right)}^{-1}{e}_3\otimes A^T{Q}_{\mathrm{DD}}^{-1}A---\left(14\right)$

$w_{\xi }[e_3^T{\left(\mathrm{T\Lambda T}\right)}^{-1}\otimes A^T{Q}_{\mathrm{DD}}^{-1}]y_{\mathrm{IF}}---\left(15\right)$

相应的协方差矩阵为： 

$\Sigma {\hat{x}}_{\mathrm{LS}}={\sigma }_0^2{N}_x^{-1}---\left(16\right)$

在步骤i中，采用不含对流层天顶延迟参数的最小二乘解的协方差矩阵代替求解降相关参数过程中涉及到的参数真值的二次型矩阵，从而计算得到降相关参数，具体为： 

采用降相关的准则为： 

${(H{\hat{\xi }}_R-y_{\mathrm{IF}})}^T{Q}_{\mathrm{IF}}^{-1}(H{\hat{\xi }}_R-y_{\mathrm{IF}})+\alpha {\hat{\xi }}_R^{T}S{\hat{\xi }}_R^T=\min ---\left(17\right)$

其中，下标“R”表示降相关解。正常数α和正定矩阵S分别是降相关参数和降相关矩阵。这两项在降相关过程中都十分重要。等式(7)的病态性是由于对流层参数和高程的强相关性造成的。因此本发明选择降相关矩阵为 

$S=\left(\begin{array}{cc}{E}_3& 0\\ 0& 0\end{array}\right)---\left(18\right)$

E₃为3×3单位矩阵。 

在准则(17)下，降相关解为 

${\hat{\xi }}_R=N_R^{-1}{w}_{\xi }---\left(19\right)$

其中N_R＝N_ξ+αS。根据误差传播定律，得出的协方差矩阵为 

${\Sigma }_{\stackrel{·}{\xi }R}={\sigma }_0^2{N}_R^{-1}{N}_{\xi }{N}_R^{-1}---2\left(20\right)$

降相关解(19)比最小二乘解(13)要稳定得多，但是这个结果是有偏的。 

采用均方误差用来评估降相关解偏差影响。 

$M_R={\sigma }_0^2{N}_R^{-1}{N}_{\xi }{N}_R^{-1}+{\alpha }^2{N}_R^{-1}S{\overline{\mathrm{\xi \xi }}}^{T}S{N}_R^{-1}---\left(\text{21}\right)$

降相关参数可以根据以下准则求得： 

α＝arg min trace(M_R)  (22) 

事实上，为不同目的，可以有许多准则来求得降相关参数。本发明选择MSE最小准则是因为本发明关注高质量的定位结果。从频率角度看，降相关参数的引入是为了减少高频信号的影响，而对于低频信号，几乎没有减少。一旦降相关参数α确定，那么相对于最小二乘解更稳定的降相关解就可以根据等式(19)求得。 

方差分量十分重要，因为α对此变量很敏感。从理论角度讲，α是用来平衡观测值和降相关偏差的因子。如果很大，那么，α也会相应变大，相应会引起较大的偏差；如果很小，那么α也会变小，那就意味着病态条件数不能有效控制，降相关解会不稳定。幸运的是，在本发明中，在计算α中，不需要用到具体解释参看下文。因此，最重要的是为未知数的真值找到合理的替代。 

通常，未知数真值可以用式(8)估计。但问题是，在这个病态模型中，最小二乘解以及它的协因数矩阵本身就不稳定。而这个就是要运用降相关的原因。当然，也可以先初始化一个小的降相关参数来计算降相关解。然后用这个解来代替真值然后再根据此真值得到更好的α，用这个α，可以求得更好的降相关解，再来代替真值通过迭代求出最优降相关解。然而，在现实中，在不同观测阶段和不同的观测几何条件下，模型(7)的条件数变化很大。而且，迭代的过程对于坐标的计算压力也很大。 

在上面提到，虽然不含对流层参数的最小二乘解是有偏估计，但是它比等式(8)求得的解稳定。因此，用代替真值是合理的。类似地，真值越大，降相关参数则越小，反之亦然。在这种情况下，通常通过迭代求取如果只是简单地把得到的结果代入等式(22)求取α，那么得到的降相关参数α一定会很大。这样会很大程度上削弱观测值的贡献，得到结果相当于经过了平滑处理，不能反映真实的情 况。此时，协方差矩阵(16)是对的合理描述。因此，在等式(21)中，我们用协方差矩阵代替二次型矩阵即 

$S{\overline{\mathrm{\xi \xi }}}^{T}S=\left(\begin{array}{cc}\Sigma {\hat{x}}_{\mathrm{LS}}& 0\\ 0& 0\end{array}\right)={\sigma }_0^2\left(\begin{array}{cc}{N}_x^{-1}& 0\\ 0& 0\end{array}\right)={\sigma }_0^2{Q}_x---\left(23\right)$

则 

$M_R={\sigma }_0^2(N_R^{-1}{N}_{\xi }{N}_R^{-1}+{\alpha }^2{N}_R^{-1}{Q}_x{N}_R^{-1})={\sigma }_0^2{Q}_R---\left(24\right)$

则准则为： 

α＝arg min_α＞0trace(Q_R)  (25) 

此法中，方差分量可以不计算。要强调的是，用最小二乘准则估计的不含对流层参数的坐标参数的协方差矩阵具有观测几何特征。因此，降相关参数的计算与观测几何变化有关。当得到合理的α，对流层参数和高程之间的相关性就可以显著降低。 

在步骤j中，确定降相关参数后，采用降相关方法重新求解，即采用式(19)和(20)解得降相关解，解决步骤g的病态问题，以求得的降相关解作为最终的亚厘米级坐标估计值，从可得到可靠的亚厘米级定位坐标估值。