CN103197335B - 采用改进正则化方法抑制dgps整周模糊度病态性的方法 - Google Patents

Abstract

本发明是一种采用改进正则化方法抑制DGPS整周模糊度病态性的方法，首先采集GPS载波相位的观测数据，建立DGPS载波相位双差观测方程；再根据观测方程，基于最小二乘方法，获得DGPS整周模糊度的浮点解和相应的方差-协方差矩阵；然后采用两步解法构造Tikhonov正则化算法中的正则化矩阵，根据DFP拟牛顿法求得相应的正则化参数，利用得到的Tikhonov正则化算法对方差-协方差矩阵进行处理，抑制DGPS整周模糊度的病态性，最后获得比较准确的整周模糊度。本发明方法采用改进的Tikhonov正则化算法来抑制DGPS整周模糊度中的病态性问题，有利于得到准确的整周模糊度，实现DGPS高精度定位和姿态测量。

Description

采用改进正则化方法抑制DGPS整周模糊度病态性的方法

技术领域

本发明涉及一种采用改进正则化方法抑制DGPS（Difference Global Positioning System，差分全球定位系统）整周模糊度病态性的方法，属于抑制DGPS整周模糊度求解过程中存在的病态问题的技术领域。

背景技术

在利用DGPS进行高精度定位和姿态测量时，最为重要的就是整周模糊度的求解。在DGPS进行高精度定位的数据处理中，一般都是基于最小二乘理论（Least Square Method）进行参数估计的，传统的处理方法是首先进行最小二乘估计，得到模糊度的浮点解，然后结合各种搜索模糊度的方法来确定模糊度。

在众多的求解整周模糊度的方法中，LAMBDA（Least-squares AMBiguity DecorrelationAdjustment）方法应用最为广泛。采用LAMBDA方法准确地确定整周模糊度需要两个前提，一是较准确地模糊度浮点解；二是较为合理的方差-协方差阵。但是对于DGPS的快速高精度定位而言，观测量之间具有较强的相关性，会导致DGPS求解中的法方程严重病态。当法方程严重病态时，法方程的求逆会不稳定，原始数据的微小变化都会造成解的较大变化，导致结果质量降低，又由于观测噪声不可避免，也将会导致DGPS整周模糊度的浮点解与准确值的差值很大，这个时候采用LAMBDA方法也难以准确求解整周模糊度，所以以往只有通过观测时间来改善法矩阵的状态和模糊度的浮点解准确性，这样也就降低了快速定位的效率，无法满足快速定位的需求。因此在DGPS整周模糊度求解时必须考虑其存在的病态问题并进行改善和解决。

在改善法矩阵病态性研究方面，主要方法有岭估计法、截断奇异值法、QR分解算法、Tikhonov正则化法和基于信噪比的正则化方法等，其中Tikhonov正则化法应用较为广泛。正则化的目的是对病态问题进行求解并得到稳定的解。在正则化方法中，其核心步骤就是选择合适的正则化矩阵和相应的参数，在标准的Tikhonov正则化方法中，采用的是单位矩阵作为正则化矩阵，而正则化参数选择为1。近些年，对于正则化方法研究中，也出现了许多正则化矩阵和参数的选取方法，如岭迹法、L-curve曲线法、广义交叉检核法（GCV）等，但是这些方法也相应存在着各自的缺陷。岭迹法是将各分量画在同一张图上，选择正则化参数使得正则化方程值大体趋于稳定，这种方法的缺点是缺乏严格的理论根据，正则化参数的选取具有一定的主观随意性；L-curve曲线法是通过选择不同的正则化参数值，经过作图拟合，得到合理的正则化参数，这种方法较为直观，但是至今也没有完整正确的理论对其进行系统的解释和说明；GCV不依赖任何先验信息和假设，仅从实际观测中提取信息，它不需要准确知道观测信息的方差因子，但利用GCV求解正则化参数时经常会出现发散现象，所求得的值也不稳定。

整周模糊度一旦出现问题，就将导致大量系统性的粗差，从而严重影响定位的精度和可靠性。

发明内容

本发明的目的是为了改善和解决DGPS整周模糊度求解过程中存在的病态性问题，提出了一种采用改进正则化方法抑制DGPS整周模糊度病态性的方法，本发明解决了DGPS整周模糊度求解中的病态问题，从而能得到稳定的整周模糊度解，实现DGPS的高精度定位和姿态测量。

本发明提供的一种采用改进正则化方法抑制DGPS整周模糊度病态性的方法，包括以下几个步骤：

步骤一：采集GPS载波相位的观测数据，建立DGPS载波相位双差观测方程；

步骤二：根据步骤一得出的DGPS载波相位双差观测方程，利用最小二乘方法求解DGPS整周模糊度的浮点解和相应的方差-协方差矩阵；

步骤三：利用Tikhonov正则化方法对方差-协方差矩阵进行处理，采用两步解法构造Tikhonov正则化方法中的正则化矩阵；

步骤四：采用最小均方误差作为标准来确定两步解法中的正则化参数α₁和α₂，相应确定正则化矩阵R₁和R₂，最终得到本发明的Tikhonov正则化方法；

步骤五：由步骤三和步骤四获得的Tikhonov正则化方法抑制整周模糊度向量的病态性之后，再利用整数最小二乘法来确定最终的整周模糊度，以应用于DGPS中。

步骤一中所述的DGPS载波相位双差观测方程的建立如下：

两个不同的GPS接收机在同一时刻i跟踪到(n+1)颗卫星，则建立的相应的GPS载波相位双差方程为：

$y_i=A_{i}x+\mathrm{\λb}+e_i;{\mathrm{\σ}}_0^2{W}_0^{-1}$

其中，i表示星历历元，y_i表示星历历元i时刻的双差载波相位观测值，为n维列向量；A_i为星历历元i时刻的包含了接收机与卫星之间视线矢量的矩阵；x为位置参数，为3维列向量；λ为载波L₁的波长；b表示整周模糊度，为n维列向量；e_i为星历历元i时刻的观测噪声；表示观测值y_i的观测方差矩阵；表示观测方差矩阵的单位权重；W₀表示正定的权矩阵；

若观测数据中有m个星历历元，则相应的DGPS载波相位双差观测方程为：

$y=\mathrm{Ax}+\mathrm{Bb}+e;{\mathrm{\σ}}_0^2{W}^{-1}$

其中，I_m表示m维的单位矩阵，双差模式下的观测值向量y＝[y₁ y₂ … y_m-1 y_m]^T，位置参数向量系数矩阵A＝[A₁ A₂ … A_m-1 A_m]^T，整周模糊度系数矩阵B＝[B₁ B₂ … B_m-1 B_m]^T；双差模式下的观测误差向量e=[e₁ e₂ … e_m-1 e_m]^T。

步骤二中所述的的整周模糊度的浮点解为：其中，方差-协方差矩阵Q＝W-WA(A^TWA)^-1A^TW，正规矩阵P=B^TQB。

步骤三中所述的Tikhonov正则化方法的函数表达式F^α(b)为：

$\min :F^{\mathrm{\α}}\left(b\right)={\left|\right|\tilde{y}-\tilde{B}\hat{b}\left|\right|}^2+\mathrm{\α}{\hat{b}}^{T}R\hat{b}$

其中，α表示非负的正则化参数，||·||表示正则矩阵的2-范数，R表示正定或半正定的正则化矩阵。

步骤三中采用两步解法构造Tikhonov正则化方法中的正则化矩阵，设R₁、R₂分别表示第一步解法和第二步解法中所对应的正则化矩阵，α₁、α₂分别表示第一步解法和第二步解法中所对应的正则化参数；

第一步解法中，选择正则化矩阵R＝R₁=I，得到相应的正则化解和偏差分别为：

${\hat{b}}_{{\mathrm{\α}}_1}={P_{{\mathrm{\α}}_1}}^{-1}{\tilde{B}}^T\tilde{y}$

$\mathrm{bias}\left({\hat{b}}_{{\mathrm{\α}}_1}\right)=E({\hat{b}}_{{\mathrm{\α}}_1}-\hat{b})=-{\mathrm{\α}}_1{P_{{\mathrm{\α}}_1}}^{-1}\hat{b}$

正规矩阵 $P_{{\mathrm{\α}}_1}={\tilde{B}}^T\tilde{B}+{\mathrm{\α}}_{1}I,$ b利用整数最小二乘方法求解：Q＝W-WA(A^TWA)^-1A^TW，为相应的方差-协方差矩阵；Zⁿ表示n维的整数向量，矩阵 $P_{\hat{b}}={\mathrm{\σ}}_0^{2}P;$

进一步，根据求得正则化解的方差矩阵 $D\left({\hat{b}}_{{\mathrm{\α}}_1}\right)$ 为： $D\left({\hat{b}}_{{\mathrm{\α}}_1}\right)={\mathrm{\σ}}_0^2{P_{{\mathrm{\α}}_1}}^{-1}{\tilde{B}}^T\tilde{B}{P_{{\mathrm{\α}}_1}}^{-1};$

正则化解的均方误差为： $\begin{array}{c}\mathrm{MSE}\left({\hat{b}}_{{\mathrm{\α}}_1}\right)=D\left({\hat{b}}_{{\mathrm{\α}}_1}\right)+\mathrm{bias}\left({\hat{b}}_{{\mathrm{\α}}_1}\right){\left[\mathrm{bias}\left({\hat{b}}_{{\mathrm{\α}}_1}\right)\right]}^T\\ ={\mathrm{\σ}}_0^2{P_{{\mathrm{\α}}_1}}^{-1}{\tilde{B}}^T\tilde{B}{P_{{\mathrm{\α}}_1}}^{-1}+{\mathrm{\α}}_1^2{P_{{\mathrm{\α}}_1}}^{-1}\hat{b}\·{\hat{b}}^T{P_{{\mathrm{\α}}_1}}^{-1}\end{array};$ 正则化参数α₁采用步骤四中方法进行确定；

第二步解法中，取正则化矩阵 $R=R_2=\mathrm{diag}{\left(\mathrm{MSE}\left({\hat{b}}_{{\mathrm{\α}}_1}\right)\right)}^{-1},$ 则得到正则化解及对应的均方误差分别为：

${\hat{b}}_{{\mathrm{\α}}_2}={P_{{\mathrm{\α}}_2}}^{-1}{\tilde{B}}^T\tilde{y}$

$\mathrm{MSE}\left({\hat{b}}_{{\mathrm{\α}}_2}\right)=D\left({\hat{b}}_{{\mathrm{\α}}_2}\right)+\mathrm{bias}\left({\hat{b}}_{{\mathrm{\α}}_2}\right){\left[\mathrm{bias}\left({\hat{b}}_{{\mathrm{\α}}_2}\right)\right]}^T$

$={\mathrm{\σ}}_0^2{P_{{\mathrm{\α}}_2}}^{-1}{\tilde{B}}^T\tilde{B}{P_{{\mathrm{\α}}_2}}^{-1}+{\mathrm{\α}}_2^2{P_{{\mathrm{\α}}_2}}^{-1}\hat{b}\·{\hat{b}}^T{P_{{\mathrm{\α}}_2}}^{-1}$

第二步解法中的正规矩阵 $P_{{\mathrm{\α}}_2}={\tilde{B}}^T\tilde{B}+{\mathrm{\α}}_2{R}_2,$ 正则化参数α₂采用步骤四中方法进行确定。

步骤五中,首先，通过Tikhonov正则化方法对方差-协方差矩阵进行处理，并结合整周模糊度的浮点解的公式，得到正规矩阵 $\tilde{P}:\tilde{P}={\tilde{B}}^T\tilde{B}+{\mathrm{\α}}_2{R}_2;$

进一步，得到整周模糊度均方误差矩阵 ${\tilde{P}}_{\hat{b}}:{\tilde{P}}_{\hat{b}}={\mathrm{\σ}}_0^2{\tilde{P}}^{-1}={\mathrm{\σ}}_0^2{({\tilde{B}}^T\tilde{B}+{\mathrm{\α}}_2{R}_2)}^{-1};$

此时，整周模糊度的浮点解为：

最后，根据式子：确定最终的整周模糊度b。

本发明的优点和积极效果在于：

（1）本发明采用两步解法构造了正则化方法中的正则化矩阵，与利用单位矩阵作为正则化矩阵相比，更加适合DGPS整周模糊度中的病态性问题，有利于得到准确的整周模糊度，从而有利于实现DGPS的高精度定位和姿态测量；

（2）本发明将DFP拟牛顿法引入到正则化方法中对正则化参数进行选取，与标准的Tikhonov正则化方法中利用1作为正则化参数相比，可得到最优的正则化参数，较好地抑制了DGPS整周模糊度中的病态性问题；

（3）本发明利用改进的Tikhonov正则化方法对DGPS整周模糊度病态问题进行求解，可以大大改善其病态性，得到更为稳定的解，从而有利于实现DGPS的高精度定位和姿态测量。

附图说明

图1为本发明抑制DGPS整周模糊度病态性的改进正则化方法的步骤流程图。

具体实施方式

下面结合附图具体对本发明的技术方案进行说明。

本发明方法首先利用整数最小二乘求解观测方程的浮点解，再将改进的正则化方法引入到整周模糊度病态性的抑制中，从而使DGPS整周模糊度求解过程中的病态性得以改善，可得到稳定的整周模糊度的解，从而在利用DGPS进行高精度定位和姿态测量时更加准确。

本发明是一种采用改进正则化方法抑制DGPS整周模糊度病态性的方法，如图1所示，各步骤具体说明如下。

步骤一：采集GPS载波相位的观测数据，建立DGPS载波相位双差观测方程；

在利用GPS载波相位进行定位时，存在着很多的误差如：卫星和接收机钟差、电离层和对流层误差、卫星轨道误差、以及噪声误差等，因此在利用GPS进行高精度姿态测量和相对定位时，一般采用站星的双差模型来处理，用来消除这些存在的误差。若两个不同的GPS接收机在同一时刻i跟踪到(n+1)颗卫星，那么相应的GPS载波相位双差方程可以线性表示为：

$y_i=A_{i}x+\mathrm{\λb}+e_i;{\mathrm{\σ}}_0^2{W}_0^{-1}---\left(11\right)$

式中，i表示星历历元，y_i表示星历历元i时刻的双差载波相位观测值，为n维列向量；A_i为星历历元i时刻的包含了接收机与卫星之间视线矢量的矩阵；x为位置参数，为3维列向量；λ为载波L₁的波长，本发明实施例中设置为19cm；b表示整周模糊度，为n维列向量；e_i为星历历元i时刻的观测噪声；表示观测值y_i的观测方差矩阵；表示观测方差矩阵的单位权重；W₀表示正定的权矩阵。

若观测数据中有m个星历历元，则相应的DGPS载波相位方程为：

$y=\mathrm{Ax}+\mathrm{Bb}+e;{\mathrm{\σ}}_0^2{W}^{-1}---\left(2\right)$

式中：y表示双差模式下的观测值向量，y＝[y₁ y₂ … y_m-1 y_m]^T；A和B分别表示位置参数向量系数矩阵和整周模糊度系数矩阵，A＝[A₁ A₂ … A_m-1 A_m]^T，B＝[B₁ B₂ … B_m-1 B_m]^T；x与式（1）中x表示的含义相同，表示位置参数向量；b与式（1）中b表示的含义相同，表示整周模糊度，为n维向量；e表示双差模式下的观测误差（噪声）向量，e=[e₁ e₂ … e_m-1 e_m]^T；式（2）中的与式（1）中的表示的含义相同，都表示观测方差矩阵的单位权重；I_m表示m维的单位矩阵。

步骤二：根据步骤一得出的DGPS载波相位双差观测方程，利用最小二乘方法求解DGPS整周模糊度的浮点解和相应的方差-协方差矩阵。

通过方程（2）得到误差方程的法方程，用加权最小二乘法求位置向量的浮点解和整周模糊度向量的浮点解。求解函数表达式如下：

$\underset{x\∈R^3,b\∈Z^n}{\min }:F={(y-\mathrm{Ax}-\mathrm{Bb})}^{T}W(y-\mathrm{Ax}-\mathrm{Bb})---\left(3\right)$

其中，R³表示3维的实数向量；Zⁿ表示n维的整数向量；F表示加权最小二乘法求解的函数表达式；W表示DGPS载波相位双差观测方程的权矩阵，如式（2）中所表示。

首先，利用式（3）求得位置参数的浮点解

$\hat{x}={\left(A^T\mathrm{WA}\right)}^{-1}{A}^{T}W(y-\mathrm{Bb})---\left(4\right)$

然后，将式（4）带入式（3），得到加权最小二乘法求解的函数表达式F(b)如下：

$\min :F\left(b\right)={(y-B\hat{b})}^{T}Q(y-B\hat{b})---\left(5\right)$

式中，Q＝W-WA(A^TWA)-¹A^TW，为相应的方差-协方差矩阵。

因此，基于最小二乘理论所求得的整周模糊度向量的浮点解根据下式可得到：

$P\hat{b}=B^T\mathrm{Qy}---\left(6\right)$

式中，P=B^TQB，是正规矩阵。

再利用整数最小二乘方法进行整周模糊度固定值的求解：

$b=\arg \underset{b\∈Z^n}{\min }{(\hat{b}-b)}^T{P}_{\hat{b}}(\hat{b}-b)---\left(7\right)$

式中，整周模糊度均方误差矩阵 $P_{\hat{b}}={\mathrm{\σ}}_0^{2}P.$

步骤三：利用Tikhonov正则化方法对方差-协方差矩阵Q进行处理，以抑制DGPS整周模糊度的病态性。本步骤根据步骤二中得出的方差-协方差矩阵Q，采用两步解法构造出Tikhonov正则化方法中相应的正则化矩阵。

为了求得有用和稳定的解，采用Tikhonov正则化方法对式（5）中的方差-协方差矩阵Q进行处理，以抑制DGPS整周模糊度的病态性。

Tikhonov正则化方法的一般形式为：

$\min :F^{\mathrm{\α}}\left(b\right)={\left|\right|y-B\hat{b}\left|\right|}_Q^2+\mathrm{\α\Ω}\left(\hat{b}\right)$              $\left(8\right)$

$={\left|\right|y-B\hat{b}\left|\right|}_Q^2+\mathrm{\α}{\hat{b}}^{T}R\hat{b}$

式中，||·||表示正则矩阵的2-范数，下角标Q表示式的方差-协方差矩阵；α表示非负的正则化参数，在标准的Tikhonov正则化方法中选取为1；R表示正定（或半正定）的正则化矩阵，在标准的Tikhonov正则化方法中选取为单位矩阵；F^α(b)表示Tikhonov正则化方法求解的函数表达式；Ω(.)表示稳定泛函。因此上式（8）变为标准的Tikhonov正则化方法表达式如下所示：

$\min :F^{\mathrm{\α}}\left(b\right)={\left|\right|y-B\hat{b}\left|\right|}_Q^2+{\hat{b}}^T\hat{b}---\left(9\right)$

当采用正则化方法来解决病态问题时，第一件事就是确定合适的正则化矩阵，本发明采用一种新的方法来构造正则化矩阵，而非采用式（9）所示的正则化方法来解决病态性。

将式（5）中的矩阵Q进行单位化，得Q＝MM^T，那么式（5）可变为：

$\min :F\left(b\right)={(\tilde{y}-\tilde{B}b)}^T(\tilde{y}-\tilde{B}b)---\left(10\right)$

式中，矩阵 $\tilde{y}=M^{T}y,$ 矩阵 $\tilde{B}=M^{T}B.$

那么方程（8）可表示为：

$\min :F^{\mathrm{\α}}\left(b\right)={\left|\right|\tilde{y}-\tilde{B}\hat{b}\left|\right|}^2+\mathrm{\α}{\hat{b}}^{T}R\hat{b}---\left(11\right)$

为了求得合适的正则化矩阵，采用新的两步解法构造正则化矩阵，设R₁、R₂分别表示第一步解法和第二步解法中所对应的正则化矩阵，α₁、α₂分别表示第一步解法和第二步解法中所对应的正则化参数。首先，对式（11）中选择R＝R₁=I，可以得到相应的正则化解和偏差

${\hat{b}}_{{\mathrm{\α}}_1}={P_{{\mathrm{\α}}_1}}^{-1}{\tilde{B}}^T\tilde{y}---\left(12\right)$

$\mathrm{bias}\left({\hat{b}}_{{\mathrm{\α}}_1}\right)=E({\hat{b}}_{{\mathrm{\α}}_1}-\hat{b})=-{\mathrm{\α}}_1{P_{{\mathrm{\α}}_1}}^{-1}\hat{b}---\left(13\right)$

式中，与式（6）中P类似，为经过正则化变换之后的正规矩阵；α₁表示第一步解法中的正则化参数；是未知整周模糊度向量的值，由式（6）求得。

根据式（12）可求得正则化解的方差矩阵 $D\left({\hat{b}}_{{\mathrm{\α}}_1}\right)$ 为：

$D\left({\hat{b}}_{{\mathrm{\α}}_1}\right)={\mathrm{\σ}}_0^2{P_{{\mathrm{\α}}_1}}^{-1}{\tilde{B}}^T\tilde{B}{P_{{\mathrm{\α}}_1}}^{-1}---\left(14\right)$

那么正则化解的均方误差可表示为：

$\begin{array}{c}\mathrm{MSE}\left({\hat{b}}_{{\mathrm{\α}}_1}\right)=D\left({\hat{b}}_{{\mathrm{\α}}_1}\right)+\mathrm{bias}\left({\hat{b}}_{{\mathrm{\α}}_1}\right){\left[\mathrm{bias}\left({\hat{b}}_{{\mathrm{\α}}_1}\right)\right]}^T\end{array}$            $\left(15\right)$

$\begin{array}{c}={\mathrm{\σ}}_0^2{P_{{\mathrm{\α}}_1}}^{-1}{\tilde{B}}^T\tilde{B}{P_{{\mathrm{\α}}_1}}^{-1}+{\mathrm{\α}}_1^2{P_{{\mathrm{\α}}_1}}^{-1}\hat{b}\·{\hat{b}}^T{P_{{\mathrm{\α}}_1}}^{-1}\end{array}$

利用步骤四中的方法确定此时的正则化参数α₁，再将α₁带入（15）式中求得此时的均方误差。接下来再取式（11）中的正则化矩阵 $R=R_2=\mathrm{diag}{\left(\mathrm{MSE}\left({\hat{b}}_{{\mathrm{\α}}_1}\right)\right)}^{-1},$ 则相应的式（11）的正则化解及均方误差表示如下：

${\hat{b}}_{{\mathrm{\α}}_2}={P_{{\mathrm{\α}}_2}}^{-1}{\tilde{B}}^T\tilde{y}---\left(16\right)$

$\mathrm{MSE}\left({\hat{b}}_{{\mathrm{\α}}_2}\right)=D\left({\hat{b}}_{{\mathrm{\α}}_2}\right)+\mathrm{bias}\left({\hat{b}}_{{\mathrm{\α}}_2}\right){\left[\mathrm{bias}\left({\hat{b}}_{{\mathrm{\α}}_2}\right)\right]}^T$            $\left(17\right)$

$={\mathrm{\σ}}_0^2{P_{{\mathrm{\α}}_2}}^{-1}{\tilde{B}}^T\tilde{B}{P_{{\mathrm{\α}}_2}}^{-1}+{\mathrm{\α}}_2^2{P_{{\mathrm{\α}}_2}}^{-1}\hat{b},{\hat{b}}^T{P_{{\mathrm{\α}}_2}}^{-1}$

式中，与（6）式中P类似，为经过正则化变换之后的正规矩阵。α₂表示第二步解法中的正则化参数。

再次利用步骤四中的方法确定此时的正则化参数α₂。利用本两步方法求得的正则化矩阵与步骤四中正则化参数的求解是紧密相关的，一旦正则化参数能够准确求得，那么相应的正则化矩阵也能相应地确定，所以接下来最为重要的就是选取合适的正则化参数了。

步骤四：根据DFP拟牛顿算法求得相应的正则化参数，正则化参数采用最小均方误差作为准则来确定。确定正则化参数α₁和α₂后，相应的正则化矩阵R₁和R₂也能相应地确定，根据式（11）得到本发明改进的正则化方法，可利用本发明改进的Tikhonov正则化方法对方差-协方差矩阵进行处理，以抑制DGPS整周模糊度的病态性。

本发明中采用最小均方误差作为标准来确定正则化参数，确定正则化参数的函数表达式f(α)表示为：

$\min :f\left(\mathrm{\α}\right)=\mathrm{tr}\left[\mathrm{MSE}\left({\hat{b}}_{\mathrm{\α}}\right)\right]$           $\left(18\right)$

$=\mathrm{tr}\left[D\left({\hat{b}}_{\mathrm{\α}}\right)\right]+{\left[\mathrm{bias}\left({\hat{b}}_{\mathrm{\α}}\right)\right]}^T\mathrm{bias}\left({\hat{b}}_{\mathrm{\α}}\right)$

其中，tr[.]表示求矩阵的迹运算。

那么分别结合式（15）和式（17），相应的（18）式可变为：

$\min :f\left({\mathrm{\α}}_1\right)={\mathrm{\σ}}_0^2\mathrm{tr}\left\{{P_{{\mathrm{\α}}_1}}^{-1}{\tilde{B}}^T\tilde{B}{P_{{\mathrm{\α}}_1}}^{-1}\right\}+{\mathrm{\α}}_1^2{\hat{b}}^T{P}_{{\mathrm{\α}}_1}^{-2}\hat{b}---\left(19a\right)$

$\min :f\left({\mathrm{\α}}_2\right)={\mathrm{\σ}}_0^2\mathrm{tr}\left\{{P_{{\mathrm{\α}}_2}}^{-1}{\tilde{B}}^T\tilde{B}{P_{{\mathrm{\α}}_2}}^{-1}\right\}+{\mathrm{\α}}_2^2{\hat{b}}^T{P}_{{\mathrm{\α}}_2}^{-2}\hat{b}---\left(19b\right)$

对式（19a）和式（19b）中的f(α₁)、f(α₂)分别对α₁、α₂求偏导：

$\▿f\left({\mathrm{\α}}_1\right)=\frac{\∂f\left({\mathrm{\α}}_1\right)}{\∂{\mathrm{\α}}_1}=-2{\mathrm{\σ}}_0^2\mathrm{tr}\left[(I-{\mathrm{\α}}_1{P}_{{\mathrm{\α}}_1}^{-1})P_{{\mathrm{\α}}_1}^{-2}\right]+2{\mathrm{\α}}_1{b}^T{P}_{{\mathrm{\α}}_1}^{-2}\hat{b}-2{\mathrm{\α}}_1^2{\hat{b}}^T{P}_{{\mathrm{\α}}_1}^{-3}\hat{b}---\left(20a\right)$

$\▿f\left({\mathrm{\α}}_2\right)=\frac{\∂f\left({\mathrm{\α}}_2\right)}{{\∂\mathrm{\α}}_2}$

$=-2{\mathrm{\σ}}_0^2\left[(R_2-{\mathrm{\α}}_2{R}_2{P}_{{\mathrm{\α}}_2}^{-1}{R}_2)P_{{\mathrm{\α}}_2}^{-2}\right]+{2\mathrm{\α}}_2{\hat{b}}^T{R}_2{P}_{{\mathrm{\α}}_1}^{-2}{R}_2\hat{b}---\left(20b\right)$

$-{\mathrm{\α}}_2^2{\hat{b}}^T{R}_2{P}_{{\mathrm{\α}}_2}^{-1}\left(R_2{P}_{{\mathrm{\α}}_2}^{-1}{+P}_{{\mathrm{\α}}_2}^{-1}{R}_2\right)P_{{\mathrm{\α}}_2}^{-1}{R}_2\hat{b}$

为了求得上式的最优解，采用DFP拟牛顿算法求解正则化参数的最优解，具体将α＝α₁与α＝α₂分别代入下面求解过程的式子中，来获取α₁和α₂的最优解，采用DFP拟牛顿法求解正则化参数的最优解过程如下：

α_(k+1)＝α_(k)+μ_(k)q_(k)      （21）

式中，k表示迭代次数；α_(k)表示在第k次迭代过程中的正则化参数；μ_(k)为第k次迭代过程中的正数常值，决定了每次迭代的步长；参数q_(k)决定了第k+1次迭代的方向。

步长μ_(k)由可由一维搜索确定：

$f({\mathrm{\α}}_{\left(k\right)}+{\mathrm{\μ}}_{\left(k\right)}{q}_{\left(k\right)})=\underset{\mathrm{\μ}>0}{\min }f({\mathrm{\α}}_{\left(k\right)}+{\mathrm{\μq}}_{\left(k\right)})---\left(22\right)$

也就是说，步长μ_(k)等于当f'(α_(k)+μq_(k))＝0所对应的值，μ表示每次一维搜索的步长。

式（21）中的q_(k)决定了下次搜索迭代的方向，它的迭代表达式为：

$q_{\left(k\right)}=-H_{\left(k\right)}\▿f\left({\mathrm{\α}}_{\left(k\right)}\right)---\left(23\right)$

式中，H_(k)为第k次迭代过程中的正规矩阵，H_(k)的表达式如下：

$H_{(k+1)}=H_{\left(k\right)}-\frac{H_{\left(k\right)}{t}_{\left(k\right)}{t}_{\left(k\right)}^T{H}_{\left(k\right)}}{t_{\left(k\right)}^T{H}_{\left(k\right)}{t}_{\left(k\right)}}+\frac{s_{\left(k\right)}{s}_{\left(k\right)}^T}{s_{\left(k\right)}^T{t}_{\left(k\right)}}---\left(24a\right)$

$s_{\left(k\right)}={\mathrm{\α}}_{(k+1)}-{\mathrm{\α}}_{\left(k\right)}---\left(24b\right)$

$t_{\left(k\right)}=\▿f\left({\mathrm{\α}}_{(k+1)}\right)-\▿f\left(a_{\left(k\right)}\right)---\left(24c\right)$

s_(k)表示相邻迭代过程中所对应的正则化参数之差；t_(k)表示相邻迭代过程中所对应的偏导之差。

在经过DFP拟牛顿方法求得相应的正则化参数之后，在与步骤三相结合，便可得到改进的正则化方法，可以解决DGPS整周模糊度病态性问题。

步骤五：确定正则化参数α₁和α₂后，正则化矩阵R₁和R₂也能相应地确定，则利用本发明改进的Tikhonov正则化方法，抑制整周模糊度的病态性，再利用整数最小二乘法确定准确的整周模糊度，将准确的整周模糊度应用于DGPS中，以实现DGPS的高精度定位和姿态测量。

通过本发明改进的Tikhonov正则化方法对式（6）进行求解，得到正规矩阵为：

$\tilde{P}={\tilde{B}}^T\tilde{B}+{\mathrm{\α}}_2{R}_2---\left(25\right)$

进一步得到相应的式（7）中的整周模糊度均方误差矩阵为：

${\tilde{P}}_{\hat{b}}={\mathrm{\σ}}_0^2\tilde{P}={\mathrm{\σ}}_0^2({\tilde{B}}^T\tilde{B}+{\mathrm{\α}}_2{R}_2)---\left(26\right)$

将矩阵分别与式（6）和（7）中的P和相比，由于正则化矩阵的参与以及合理的正则化参数，改善了式（6）中矩阵P的状态，使得矩阵P的求逆变得正常，不会出现病态的问题，因而也能够得到较为可靠的估值。因此整周模糊度的浮点解可以利用下式求得：

$\hat{b}={\tilde{P}}^{-1}B\tilde{y}---\left(27\right)$

新解法在利用式（27）得到较为准确模糊度浮点解的同时，应用式（26）的均方误差矩阵代替方差-协方差矩阵来确定整周模糊度向量的搜索范围，来确定整周模糊度向量。那么相应的式（7）可变为：

$b=\arg \underset{b\∈Z^n}{\min }{(\hat{b}-b)}^T{\tilde{P}}_{\hat{b}}(\hat{b}-b)---\left(28\right)$

此时利用式（28）求解所得的整周模糊度向量才是可靠的较为准确的整周模糊度向量。

Claims (4)

1.一种采用改进正则化方法抑制DGPS整周模糊度病态性的方法，其特征在于，包括以下几个步骤：

步骤一：采集GPS载波相位的观测数据，建立DGPS载波相位双差观测方程；

两个不同的GPS接收机在同一时刻i跟踪到(n+1)颗卫星，则建立的相应的GPS载波相位双差方程为：

$y_i=A_{i}x+\mathrm{\λb}+e_i;{\mathrm{\σ}}_0^2{W}_0^{-1}---\left(1\right)$

其中，i表示星历历元，y_i表示i时刻的双差载波相位观测值，为n维列向量；A_i为星历历元i时刻的包含了接收机与卫星之间视线矢量的矩阵；x为位置参数向量，为3维列向量；λ为载波L₁的波长；b表示整周模糊度，为n维列向量；e_i为星历历元i时刻的观测噪声；表示观测值y_i的观测方差矩阵；表示观测方差矩阵的单位权重；W₀表示正定的权矩阵；

若观测数据中有m个星历历元，则相应的DGPS载波相位双差观测方程为：

$y=Ax+\mathrm{Bb}+e;{\mathrm{\σ}}_0^2{W}^{-1}---\left(2\right)$

其中，W表示DGPS载波相位双差观测方程的权矩阵，I_m表示m维的单位矩阵，双差模式下的观测值向量y＝[y₁ y₂ … y_m-1 y_m]^T，位置参数向量系数矩阵A＝[A₁ A₂ … A_m-1 A_m]^T，整周模糊度系数矩阵B＝[B₁ B₂ … B_m-1 B_m]^T；双差模式下的观测误差向量e＝[e₁ e₂ … e_m-1 e_m]^T；

步骤二：根据步骤一得出的DGPS载波相位双差观测方程，基于最小二乘方法，获得DGPS整周模糊度的浮点解和相应的方差-协方差矩阵Q；

所述的整周模糊度的浮点解根据下式获得：

$P\hat{b}=B^T\mathrm{Qy}---\left(3\right)$

方差-协方差矩阵Q＝W-WA(A^TWA)^-1A^TW，正规矩阵P＝B^TQB；

利用整数最小二乘方法进行整周模糊度b的固定值的求解：

$b=\arg \underset{b\∈Z^n}{\min }{(\hat{b}-b)}^T{P}_{\hat{b}}(\hat{b}-b)---\left(4\right)$

Zⁿ表示n维的整数向量，整周模糊度均方误差矩阵

步骤三：采用Tikhonov正则化方法对方差-协方差矩阵进行处理，采用两步解法构造Tikhonov正则化方法中的正则化矩阵；

所述的Tikhonov正则化方法的函数表达式F^α(b)为：

$\min :F^{\mathrm{\α}}\left(b\right)={\left|\right|\tilde{y}-\tilde{B}\hat{b}\left|\right|}^2+\mathrm{\α}{\hat{b}}^{T}R\hat{b}---\left(5\right)$

其中，α表示非负的正则化参数，表示矩阵的2-范数，R表示正定或半正定的正则化矩阵，矩阵矩阵M获取方法：将方差-协方差矩阵Q进行单位化得到Q＝MM^T；

设R₁、R₂分别表示第一步解法和第二步解法中所对应的正则化矩阵，α₁、α₂分别表示第一步解法和第二步解法中所对应的正则化参数；

第一步解法中，式(5)中选择正则化矩阵R＝R₁＝I，得到相应的正则化解和偏差分别为：

${\hat{b}}_{{\mathrm{\α}}_1}={P_{{\mathrm{\α}}_1}}^{-1}{\tilde{B}}^T\tilde{y}---\left(6\right)$

$\mathrm{bias}\left({\hat{b}}_{{\mathrm{\α}}_1}\right)=E({\hat{b}}_{{\mathrm{\α}}_1}-\hat{b})=-{\mathrm{\α}}_1{P_{{\mathrm{\α}}_1}}^{-1}\hat{b}---\left(7\right)$

正规矩阵b根据式(4)求解得到；

进一步，根据式(6)求得正则化解的方差矩阵为：

$D\left({\hat{b}}_{{\mathrm{\α}}_1}\right)={\mathrm{\σ}}_0^2{P_{{\mathrm{\α}}_1}}^{-1}{\tilde{B}}^T\tilde{B}{P_{{\mathrm{\α}}_1}}^{-1}---\left(8\right)$

那么正则化解的均方误差为：

$\begin{array}{c}\mathrm{MSE}\left({\hat{b}}_{{\mathrm{\α}}_1}\right)=D\left({\hat{b}}_{{\mathrm{\α}}_1}\right)+\mathrm{bias}\left({\hat{b}}_{{\mathrm{\α}}_1}\right){\left[\mathrm{bias}\left({\hat{b}}_{{\mathrm{\α}}_1}\right)\right]}^T\\ ={\mathrm{\σ}}_0^2{P_{{\mathrm{\α}}_1}}^{-1}{\tilde{B}}^T\tilde{B}{P_{{\mathrm{\α}}_1}}^{-1}+{\mathrm{\α}}_1^2{P_{{\mathrm{\α}}_1}}^{-1}\hat{b}\·{\hat{b}}^T{P_{{\mathrm{\α}}_1}}^{-1}\end{array}---\left(9\right)$

正则化参数α₁采用步骤四中方法进行确定，将确定的α₁带入式(9)中求得此时的均方误差；第二步解法中，取式(5)中的正则化矩阵则第二步解法中正则化解及对应的均方误差分别为：

${\hat{b}}_{{\mathrm{\α}}_2}={P_{{\mathrm{\α}}_2}}^{-1}{\tilde{B}}^T\tilde{y}---\left(10\right)$

$\begin{array}{c}\mathrm{MSE}\left({\hat{b}}_{{\mathrm{\α}}_2}\right)=D\left({\hat{b}}_{{\mathrm{\α}}_2}\right)+\mathrm{bias}\left({\hat{b}}_{{\mathrm{\α}}_2}\right){\left[\mathrm{bias}\left({\hat{b}}_{{\mathrm{\α}}_2}\right)\right]}^T\\ ={\mathrm{\σ}}_0^2{P_{{\mathrm{\α}}_2}}^{-1}{\tilde{B}}^T\tilde{B}{P_{{\mathrm{\α}}_2}}^{-1}+{\mathrm{\α}}_2^2{P_{{\mathrm{\α}}_2}}^{-1}\hat{b}\·{\hat{b}}^T{P_{{\mathrm{\α}}_2}}^{-1}\end{array}---\left(11\right)$

第二步解法中的正规矩阵正则化参数α₂采用步骤四中方法进行确定；

步骤四：采用最小均方误差作为准则来确定正则化参数α₁和α₂，相应确定正则化矩阵R₁和R₂，最终根据式(5)得到Tikhonov正则化方法；

步骤五：通过Tikhonov正则化方法对方差-协方差矩阵进行处理后，再利用整数最小二乘法来确定最终的整周模糊度，以应用于DGPS中；

首先，通过Tikhonov正则化方法对方差-协方差矩阵进行处理，并结合式(3)得到正规矩阵 $\tilde{P}:\tilde{P}={\tilde{B}}^T\tilde{B}+{\mathrm{\α}}_2{R}_2;$

进一步，得到式(4)中的整周模糊度均方误差矩阵

此时，整周模糊度的浮点解为：

最后，根据式子：确定最终的整周模糊度b。

2.根据权利要求1所述的一种采用改进正则化方法抑制DGPS整周模糊度病态性的方法，其特征在于，所述的步骤二获取DGPS整周模糊度的浮点解的方法是：

通过方程(2)得到误差方程的法方程，用加权最小二乘法求解位置参数向量的浮点解和整周模糊度的浮点解函数表达式F如下：

$\underset{x\∈R^3,b\∈Z^n}{\min }:F={(y-\mathrm{Ax}-\mathrm{Bb})}^{T}W(y-\mathrm{Ax}-\mathrm{Bb})---\left(12\right)$

其中，R³表示3维的实数向量，Zⁿ表示n维的整数向量；

利用式(12)求得位置参数向量的浮点解

$\hat{x}={\left(A^T\mathrm{WA}\right)}^{-1}{A}^{T}W(y-\mathrm{Bb})---\left(13\right)$

然后，将式(13)带入式(12)，得到用加权最小二乘法求解整周模糊度的函数表达式F(b)：

$\min :F\left(b\right)={(y-B\hat{b})}^{T}Q(y-B\hat{b})---\left(14\right)$

最后，根据式(3)获得整周模糊度的浮点解

3.根据权利要求1或2所述的一种采用改进正则化方法抑制DGPS整周模糊度病态性的方法，其特征在于，步骤三中所述的Tikhonov正则化方法的函数表达式，根据下面过程得到：

Tikhonov正则化方法的形式为：

$\min :F^{\mathrm{\α}}\left(b\right)={\left|\right|y-B\hat{b}\left|\right|}_Q^2+\mathrm{\α}{\hat{b}}^{T}R\hat{b}---\left(15\right)$

其中，下角标Q表示式的方差-协方差矩阵；

将方差-协方差矩阵Q进行单位化，得Q＝MM^T，则用加权最小二乘法求解整周模糊度的函数表达式F(b)为：

$\min :F\left(b\right)={(\tilde{y}-\tilde{B}b)}^T(\tilde{y}-\tilde{B}b)---\left(16\right)$

将式(15)和式(16)结合，得到式(5)。

4.根据权利要求1所述的一种采用改进正则化方法抑制DGPS整周模糊度病态性的方法，其特征在于，步骤四中所述的采用最小均方误差作为准则来确定正则化参数α₁和α₂，具体是：

首先，确定采用最小均方误差作为标准来确定正则化参数的函数f(α)的表达式：

$\begin{array}{c}\min :f\left(\mathrm{\α}\right)\mathrm{tr}\left[\mathrm{MSE}\left({\hat{b}}_{\mathrm{\α}}\right)\right]\\ =\mathrm{tr}\left[D\left({\hat{b}}_{\mathrm{\α}}\right)\right]+{\left[\mathrm{bias}\left({\hat{b}}_{\mathrm{\α}}\right)\right]}^T\mathrm{bias}\left({\hat{b}}_{\mathrm{\α}}\right)\end{array}---\left(17\right)$

其中，tr[.]表示求矩阵的迹运算；

将式(9)和式(17)结合，将式(11)和式(17)结合，分别得到如下两式：

$\min :f\left({\mathrm{\α}}_1\right)={\mathrm{\σ}}_0^2\mathrm{tr}\left\{{P_{{\mathrm{\α}}_1}}^{-1}{\tilde{B}}^T\tilde{B}{P_{{\mathrm{\α}}_1}}^{-1}\right\}+{\mathrm{\α}}_1^2{\hat{b}}^T{P}_{{\mathrm{\α}}_1}^{-2}\hat{b}---\left(18a\right)$

$\min :f\left({\mathrm{\α}}_2\right)={\mathrm{\σ}}_0^2\mathrm{tr}\left\{{P_{{\mathrm{\α}}_2}}^{-1}{\tilde{B}}^T\tilde{B}{P_{{\mathrm{\α}}_2}}^{-1}\right\}+{\mathrm{\α}}_2^2{\hat{b}}^T{P}_{{\mathrm{\α}}_2}^{-2}\hat{b}---\left(18b\right)$

式(18a)和式(18b)中的f(α₁)、f(α₂)分别对α₁、α₂求偏导，得到：

$\▿f\left({\mathrm{\α}}_1\right)=\frac{\∂f\left({\mathrm{\α}}_1\right)}{{\∂\mathrm{\α}}_1}=-{2\mathrm{\σ}}_0^2\mathrm{tr}\left[(I-{\mathrm{\α}}_1{P}_{{\mathrm{\α}}_1}^{-1})P_{{\mathrm{\α}}_1}^2\right]+{2\mathrm{\α}}_1{\hat{b}}^T{P}_{{\mathrm{\α}}_1}^{-2}\hat{b}-{2\mathrm{\α}}_1^2{\hat{b}}^T{P}_{{\mathrm{\α}}_1}^{-3}\hat{b}---\left(19a\right)$

$\begin{array}{c}\▿f\left({\mathrm{\α}}_2\right)=\frac{\∂f\left(\mathrm{\α}2\right)}{{\∂\mathrm{\α}}_2}\\ =-{2\mathrm{\σ}}_0^2\left[(R_2-{\mathrm{\α}}_2{R}_2{P}_{{\mathrm{\α}}_2}^{-1}{R}_2)P_{{\mathrm{\α}}_2}^{-2}\right]+{2\mathrm{\α}}_2{\hat{b}}^T{R}_2{P}_{{\mathrm{\α}}_1}^{-2}{R}_2\hat{b}\\ -{\mathrm{\α}}_2^2{\hat{b}}^T{R}_2{P}_{{\mathrm{\α}}_2}^{-1}(R_2{P}_{{\mathrm{\α}}_2}^{-1}+P_{{\mathrm{\α}}_2}^{-1}{R}_2)P_{{\mathrm{\α}}_2^{}}^{-1}{R}_2\hat{b}\end{array}---\left(19b\right)$

采用DFP拟牛顿法来确定正则化参数的最优解，具体将α＝α₁与α＝α₂分别代入下面求解过程的式子中，来获取α₁和α₂的最优解，采用DFP拟牛顿法求解正则化参数的最优解过程如下：

α_(k+1)＝α_(k)+μ_(k)q_(k)(20)

其中，k表示迭代次数；α_(k)表示在第k次迭代过程中的正则化参数；第k次迭代过程中的正数常值μ_(k)决定了每次迭代的步长；参数q_(k)决定了第k+1次迭代的方向；

步长μ_(k)由一维搜索确定：

$f({\mathrm{\α}}_{\left(k\right)}+{\mathrm{\μ}}_{\left(k\right)}{q}_{\left(k\right)})=\underset{\mathrm{\μ}>0}{\min }f({\mathrm{\α}}_{\left(k\right)}+{\mathrm{\μq}}_{\left(k\right)})---\left(21\right)$

步长μ_(k)等于当f'(α_(k)+μq_(k))＝0所对应的值，μ表示每次一维搜索的步长；

第k次迭代过程中的参数q_(k)通过下式获得：

$q_{\left(k\right)}=-H_{\left(k\right)}\▿f\left({\mathrm{\α}}_{\left(k\right)}\right)---\left(22\right)$

第k次迭代过程中的正规矩阵H_(k)为：

$H_{(k+1)}=H_{\left(k\right)}-\frac{H_{\left(k\right)}{t}_{\left(k\right)}{t}_{\left(k\right)}^T{H}_{\left(k\right)}}{t_{\left(k\right)}^T{H}_{\left(k\right)}{t}_{\left(k\right)}}+\frac{s_{\left(k\right)}{s}_{\left(k\right)}^T}{s_{\left(k\right)}^T{t}_{\left(k\right)}}---\left(23a\right)$

s_(k)＝α_(k+1)-α_(k)(23b)

$t_{\left(k\right)}=\▿f\left({\mathrm{\α}}_{(k+1)}\right)-\▿f\left({\mathrm{\α}}_{\left(k\right)}\right)---\left(23c\right)$

其中，s_(k)表示相邻迭代过程中的正则化参数之差；t_(k)表示相邻迭代过程中所对应的偏导之差。