刘基余，陈小明

（武汉大学测绘学院，武汉 430079）

利用卡尔漫滤波处理GNSS动态测量数据，需要首先建立滤波的动态模型（状态方程）和观测模型（观测方程）。对GNSS动态定位而言，动态模型难以用精确的数学公式来表达，实用中一般在精度损失较小的前提下对动态模型进行简化。最常用的动态模型为常速模型或常加速模型。常速模型或常加速模型的选择依赖于运动载体的运动状态以及数据采样率的高低，在高精度GNSS动态定位中数据采样率一般为1秒或更高，此时可采用常速模型。动态噪声可假定为零均值高斯噪声。在整周模糊度已知时（利用基线初始化或模糊度在航解算OTF求得），单纯利用双差载波相位观测量就可以获得极高的动态定位精度。

由于载波相位测量的整周模糊度特性，在载波相位测量出现周跳后，整周模糊度会发生改变。若不能及时修正，将严重影响动态定位的精度。在滤波设计时，需要考虑对周跳的探测与修复。另一个方面要考虑到动态接收机机动情形，在滤波设计时，更需考虑对机动加速度的补偿。中国魏明博士等（1990）和卢刚博士等（1991）分别提出了用统计检验的方法来探测周跳和机动加速度，并利用两步卡尔曼滤波来估计周跳的方法；卢刚博士等（1990）称之为GPS动态定位的统计质量控制（Statistical quantity control）。为了研究质量控制的有效性，本文提出利用可靠性理论来分析周跳探测，机动探测的能力以及周跳与机动加速度的可区分性问题，并在此基础上提出一种更为简单实用的高精度GNSS动态定位模型。

1 基本滤波模型及其递推方程

状态方程：

观测方程：

系统动态噪声和量测噪声的统计特性假设如下：

式（1a）中动态模型采用常速模型。状态转移矩阵Φk，k-1和干扰矩阵Γk-1定义如下：

式中，I3为3×3维单位阵；ΔT为采样时间间隔。

式（1b）为采用双差载波相位观测量的观测方程。双差载波相位观测值可选为L1、L2单频载波相位或双频组合观测值，其频率为f，波长为λ=C//f。该式为非线性方程。令式（1b）围绕参考估计（一般选为当前时刻的预报值Xk，K-1）线性化，可得到线性化观测方程

或记为

式中，YK为双差载波相位观测向量；为相位的双差整周模糊度向量。

卡尔曼滤波的初值为X0/0，P0/0。

2 GPS动态定位的质量控制

在前述的基本滤波模型中未考虑运动载体的机动以及观测值中可能存在的粗差（周跳）。由于机动加速度和周跳可以常值为差来模拟，在出现机动或发生周跳时可在滤波状方程和观测方程中加入相应的修正量，此时滤波模型为：

式中，b，d分别表示机动加速度和双差周跳。当不存在机动加速度时，b=0；当不存在周跳时，d=0；当b，d同时为零时，模型（4a）、（4b）退化为（1a），（2b）。BK，DK分别为机动加速度、双差周跳对应的系数阵。BK定义如下：

DK定义如下：对应第i个双差周跳，有

若已知机动加速度b和双差周跳d，则时刻k动态接收机天线的位置、速度预报值和滤波分别为：

当存在机动加速度和（或）双差周跳，但仍采用模型式（1a）和式（3b），而利用式（3a）～（3f）进行递推计算值时，递推估值，将不是最优无偏估计而是有偏的。当观测方程是非线性方程时，有偏估计与无偏估计之间的偏差关系也是非线性的。可以证明，当忽略线性化误差时，无偏估计和有偏估计间的偏差值可用b、d的线性关系式表达。下面进行详细推导。

设有偏估计和无偏估计差值可用下列线性关系式表达

比较式（7）两端系数，可得到递推关系

将观测方程围绕参考估计线性化，有

比较式（9）两端系数可得到

式（8a）、式（10a）、式（8b）和式（10b）分别构成递推计算公式。在b、d已知时，可将无偏差参数模型式（1a）和式（2b）计算得到的有偏估计修正为最优无偏估计。b、d可由两步卡尔曼滤波的方法来确定。

关键的问题在于如何探测何时出现机动、周跳，以及机动出现在哪个方向上、周跳出现在哪个双差观测值上。机动以及周跳的探测可以通过对预报残差进行统计检验来完成。

根据卡尔曼滤波理论，在滤波模型存在偏差或对偏差进行了修正时，预报残偏差应该是零均值的白噪声。若实际系统存在偏差，但在滤波模型中未进行修正，则预报残差将不再是零均值的白噪声。当存在偏差b、d时，采用滤波模型式（4a）和式（4b）的预报残差为：

若采用滤波模型式（6a）和式（6b），则预报残差为：

由于b、d是常值偏差，因而rK服从均值为GKU，方差为的正态分布

而不存在偏差即b=0，d=0时预报残差则服从零均值，方差为的正态分布

根据采用无偏差参数模型式（2a）和式（2b）计算得到的新息序列在有、无偏差时的统计特性，可以构造如下原假设H0和备选假设Ha：

上述检验可构造如下统计检验量

式中，nK为偏差向量U的维数；λK为非中心化参数，定义为

上式所构成的偏差统计量可分为两种情况；当l=k时，仅利用当前观测时元的预报残差构成统计量，称为局部模型检验 LMT（Local Model Test）；

当K＞l时，则利用多个时元的预报残差构成统计量，称为全局模型检验GMT（Global Model Test）。全局模型检验探测系统模型偏差的能力比局部检验强，但存在检验延迟的问题，且计算量比局部检验大。

根据假设检验理论，选定所需检测的偏差，构造相应系数阵GK并选定某一显著性水平a，并取F=（nk）可以得出如下判断标准：当TG＜F时，接受原假设H0；当TG≥F时，接受备选假设Ha。

从理论上讲，这种检验方法是成立的。但从实际应用来看，此方法很难实行，因为在滤波计算时我们并不知道哪个方向发生机动或哪个观测值有周跳。因此，实用中常采用以下方法。

首先选取偏差向量的维数等于预报残差的维数，此时Gi阵为可逆方阵，此时统计量TG变为如下形式

在H0下有

式中，mi为第i个观测时元预报残差的维数。上式实际上只是一个预警检验统计量，它仅能用于判断系统是否出现故障。当k=l时，即仅利用当前观测的预报残差，称为局部预警检验。此时在H0下有

选择一显著性水平a，并取接受H0，即认为既无机动也无周跳，当≥F1时，接受Ha，即认为可能存在机动或周跳。此时应对故障原因进行进一步诊断，从而对偏差源进行定位。中国魏明博士和卢刚博士建议采用类似巴尔达的数据探测（Data Snooping）的方法进行，即对每假定的偏差源分别采用一维的检验一统计量进行定位检验，直到检测出全部可能的偏差。

当检测机动加速度时，则有

或者写成

式（23a）和式（23b）为检验机动加速度的全局一维定位检验量，若令k=1，则形成局部一维定位检验量，此时有

式（24a）和（24b）为检验周跳的全局一维定位检验量，若令k=1则形式相应的局部一维定位检验量，此时有

卢刚博士（1991）认为局部预警检验和局部定位检验对GPS相对定位中观测值进行偏差的检验和定位已足够，因为最常见的偏差如周跳等均影响当前时刻的滤波值。

在实用中，我们发现上述局部统计量对周跳很敏感而对机动加速度则不甚敏感，往往在机动发生若干时元之后才能探测到。另外，当QK选择大于某一数值时，或者选得更大，即使不对机动加速度引起的滤波估计偏差进行补偿也能获得很好的滤波结果。下面我们用可靠性理论对上述现象进行研究分析。

3 可靠性分析

可靠性研究建立在数理统计假设检验的基础上。经典的假设检验理论是1955年由莱曼和皮尔孙提出的。在测量平范畴内，可靠性研究理论是由荷兰巴尔达教授在1967～1968年提出的。巴尔达的可靠性理论是从单个一维备选假设出发，研究平差系统发现单个模型误差的能力和不可发现的模型误差对平差结果的影响。前者称为内部可靠性，后者称为外部可靠性。这里的模型误差包括粗差和系统误差。此外，从已知单位权方差出发，巴尔达教授还导出了检验粗差的数据探测法（Data Snooping），即以服从正态分布的标准化残差作为统计量。

随后，由Fotstner和Koch等人将该理论推广到单个多维备选假设，从而研究系统发现多个模型误差的能力。1983年Fotstner第一次提出了模型误差的区分可能性，并从两个一维备选假设出发，导出了区分可能性本质上取决于检验量之间相关系数的结论。李德仁院士在他的博士论文中从高斯-马尔可夫模型含两个多维备选假设出发，提出了平差系统的可区分性和可靠性理论。该理论可研究全系统发现并区分不同模型误差的区分和定位提供了研究的基本理论和定量的尺度。

近年来，可靠性研究理论已在大地测量、摄影测量、工程测量及形变测量中取得了广泛的应用。这一理论也已被引入到集成导航系统的设计之中。下面应用可靠性理论来研究卡尔曼滤波处理高精度GPS动态定位数据的内外部可靠性。

3.1 内部可靠性分析

在单个备选假设下，内部可靠性研究的是若检验以一定的显著水平进行时能以什么样的把握（检验功效βP）发现模型误差。在多数情况下，由于模型误差的大小是未知的，内部可靠性主要研究至少要出现多大的模型误差，它才能在所规定的检验功效β0在显著性水平为α0的检验中发现。在两个备选假设下，可靠性理论主要研究的是不同备选假设下模型误差的可区分性（可区分性理论）。弗斯特勒尔（1983）指出两个检验量之间的相关系数可作为衡量可区分性的指标。李德仁教授在其博士论文中从两个多维备选假设的情况下，导出了检验量之间的总体相关和最大相关，并进而导出了可发现且可与它种模型误差相区分的模型误差下界值。

鉴于在前述质量控制模型中以类似巴尔达的数据探测法来探测周跳和机动加速度，因而在下面的可靠性研究中主要考虑备选假设为一维时的内部可靠性分析。

首先考虑单个一维备选假设的情形。

对于单个方向机动加速度，有零假设H0：不存在机动加速度

备选假设Ha：存在一个方向的机动加速度

根据内部可靠性理论，可发现的机动加速度下界定义为

式中，δ0为非中心化参数，它可在给定显著水平α0和检验功效β0下由标准化正态分布表求出。同样对于单个观测值粗差，有

零假设H0：不存在观测值粗差

备选假设Ha：存在一个观测值粗差可发现观测值粗差的下界定义为：

当仅选择当前时元的预报残差进行检验，即采用局部一维定位检测时，k=l，b，d分别定义为：

对于两个一维备选假设，即

H0：无机动，无观测值粗差

Ha1：存在一个方向机动加速度

Ha2：存在一个观测值粗差

此时，内部可靠性要研究发现并区分两类模型误差的能力。可区分性取决于两个备选假设之间的相关系数。单个方向机动加速度和单个观测值粗差检验统计量的相关系数定义为

在k=l时，即选用统计量为局部一维定位统计量时

相关系数越大，则模型误差的可区分性越差，当ρbd=1时，机动加速度与观测值粗差不可区分。根据相关系数ρbd可以查表得到可区分性放大倍数Kρbd。该值表示由于相关引起的非中心化参数的放大倍数，更具体地讲，它代表可区分且可发现的模型误差下界为可发现模型误差下界值的多少倍。

可发现且可与单个机动加速度区分的观测值粗差下界为

可发现且可与单个观测值粗差区分的机动加速度下界为

研究中发现，对于本节的滤波模型，内部可靠性主要与动态噪声QK和观测噪声RK有关，在RK一定的情况下，随着给定QK的增大可发现机动加速度和观测值粗差的下界呈增大趋势，而机动加速度与观测值粗差的相关系数则呈减小趋势。对于相同的QK，当观测值精度较高时，可发现机动加速度和观测值差的下界值较小。

3.2 可靠性分析实践

下面以1996年8月14日在哈尔滨进行的一次秒采样率机载GPS动态定位为例予以说明（如表1～表6所示）。该次定位中所采用的机载GPS信号接收机和基准接收机均为Trimble 4000SSE双频接收机。表1～表5为采用局部检验量时的内部可靠性分析，其中各量和计算均只采用滤波递推估计10个时元后的一个时元（GPS时间288880s）的数据（该时元PDOP=2.8，共观测下列6颗卫星：PRN26，06，27，16，17，18。其中PRN26号卫星高度角最大，其他卫星分别与它构成双差观测值）。表3～表4为采用宽巷观测值进行滤波计算时的内部可靠性分析数据，表4～表5为采用L1单频观测值进行滤波计算时的内部可靠性分析数据。计算中假定L1，L2载波相位量测精度等于0.025周，因而对于L1双差观测值，有协方差阵

对于宽巷双差载波相位观测值有

式中，λL1，λWL分别为L1载波相位和宽巷载波相位的波长（单位：cm）。

表1 QK与可发现机动加速度关系（宽巷）

表2 QK与可发现观测值粗差的关系（宽巷）

动态噪声协方差阵定义为

Q为标量，在各种QK阵设计中分别选取Q=le-8，le-4，le-2，0.1，1，10，100。各表计算中选取δ0为显著性水平a0=0.1%，检验功效为80%时的非中心化参数δ0=4.13。

由表1可知，当RK一定时（宽巷观测值），随QK设定值的增大可发现机动加速度下界值呈增大趋势。当Q=10时，可发现机动加速度下界约为重力加速度的2.6倍，超过了绝大多数运动载体的最大机动能力。而当Q=le-8时，所能发现的机动加速度下界约为0.2～0.4m。这是由于宽巷观测量测量噪声较大引起的。当选择观测噪声较小的观测量L1载波相位时，对应相同的QK可发现的机动加速度下界相对较小。

表3 QK与相关系数

表4 QK与可发现机动加速度关系（L1）

表5 QK与可发现观测值粗差（L1）

从整体上看（尤其是QK较大时），局部检验量对机动加速度不敏感，这是前面提到的采用局部检验量时会产生机动探测延迟的主要原因，为了较好地探测机动，应考虑采用全局检验量。

由表2可知，RK一定QK增大时，可发现的载波相位粗差的下界也呈增大的趋势，但增大速度明显比可发现的机动加速度下界要慢，而且呈现出一种饱和趋势。在Q≥0.1增大速度明显变慢，对于不同的双差载波相位，可发现观测值粗差的下界值不同。更值得注意的是对于各种Q值可发现的粗差下界均小于1周，即1周大小的周跳是完全可以发现的。而此时我们所采用的仅是局部定位检验统计量，这也是卢刚（1991）认为仅利用局部检验量就可以探测出小周跳的原因。

比较表5和表2还可以发现，当选择观测精度较高的L1观测值时，可发现粗差的下界值更小，随QK的增大趋于饱和的速度更快，因而探测周跳的能力也更强。

表3中未列出机动加速与所有双差载波相位的相关系数，仅列出了X、Y、Z方向与双差载波相位相关系数值的最大值。表中第一列的“/”上为机动加速度方向，“/”下为对应最大相关系数的双差载波相位卫星号。从表中可知，当Q较小时相关系数很大，甚至大到0.977。而随着Q的增大，相关系数呈现逐渐减小的趋势。而与可发现的机动加速度和粗差的趋势正好相反。在由相应的相关系数查表得到可区分性放大倍数并乘以表3中可发现粗差下界值后可得到可与机动加速度区分且可发现的粗差下界值（见表2中“-”下方）。相应的区分可能性1-γ0选为95%。这一可区分且可发现的粗差下界值仍小于1 周，这表明即使对于1周的周跳也是可以发现并与机动加速度相区分。

以上实际是处理大量高精度GPS动态定位内部可靠性分析的基本特点。根据上述实例分析，可以得到以下几个基本结论：

（1）对于周跳的探测，可仅采用局部检验量。由于载波相位量测精度高，即使动态噪声选择得较大，仅采用局部检验量，周跳均可发现且与机动加速度相区分，观测量精度越高（RK越小）周跳的探测能力越强。

（2）对于机动加速度的检测，应采用全局检验量。局部检验量对机动加速度的检测不甚敏感。当采用多个时元构造全局检验时，随着选择时元数的增多，机动加速度检测能力逐步增强。

（3）当QK较大时，可发现机动加速度下界值往往会超过运动载体的最大机动能力，然而此时对周跳探测能力仍优于1周。此时，可不必考虑对机动加速度的检测，而仅用数据探测法检测周跳。

表6 采用全局检验量时QK与可发现机动加速度关系（宽巷）

表7 SbKK相应元素（宽巷）

3.3 外部可靠性研究

外部可靠性研究的是不可发现的模型误差对平差结果的影响。在对周跳进行修正后，仅存在机动加速度，此时机动加速度时滤波值的影响为：

由式（31）可知机动加速度对滤波值的影响是由机动加速度与中相应元素相乘得到的，由于直角坐标系轴垂直的，因而沿X、Y或Z轴方向的单方向机动加速度，仅影响该方向的坐标和速度。

仍采用与检验内部可靠性相同的数据，采用宽巷观测值，利用式（31）进行了外部可靠展性分析（表7）

表7中列出了与X方向位置、速度滤波估计偏差有关的中相应的元素（K=1，2，3，4）。由表列数据乘以X方·向机动加速度，可以得到这几个时元的X，X滤波估计偏差。由表列数据可知，当给定QK较小时（Q＜le-2）时，随着时元数的增加，由机动加速度引起的位置、速度偏差呈递增趋势，机动加速度会引起滤波估计的变坏甚至发散。而QK较大时（Q＞0.1）时，随着时元数的增加，由机动加速度引起的位置估计误差趋于常值。此常值随QK的增大逐渐趋于0。速度估计偏差则趋于机动加速度的1/2。此时滤波得到的速度估值实际上是该时元与前一时元速度的平均值，即存在机动加速度时，采用无偏差参数模型，速度滤波估计随QK增大逐渐趋于平均速度。

以上实例可以作如下解释：RK一定，当QK较小时，滤波增益也小，即过去观测在滤波中的加权较大，此时滤波估计对过去观测的依赖性大，在出现机动后，会引起估计偏差的增大。当QK较大时滤波对于过去观测的依赖性或者说受过去观测的影响也相对较小。在出现机动后滤波估计主要取决于当前的观测，滤波估计偏差随时元数增加的影响也较小，而且趋于较小的稳定值（滤波位置估计偏差趋于0）。另外，由于双差载波相位观测对速度的偏导数为0，即载波相位观测量不直接对速度分量进行观测。QK较大时，速度估计基本上由在时元和上一时元的位置求平均得到，此时的速度估计为平均速度，而不是当前时刻的瞬时速度

4 简单、实用的高精度GPS动态定位滤波处理模型

由前述的内部可靠性研究可知，采用局部检验量周跳均可发现且与机动加速度相区分。观测量精度越高，周跳探测的能力也越强。且当QK较大时，可发现机动加速度下界值往往会超过运动载体的最大机动能力，机动加速度的检测应采用全局检验量。

由外部可靠性研究可知，在QK较小时不可发现的机动加速度会引起滤波估计变坏甚至发散，而QK较大时，随时元数的增多，由机动加速度引起的位置估计误差趋于常值，此常值随QK的增大趋于零。随QK的增大，速度估计偏差趋于机动加速度的1/2，即速度滤波估值趋于平均速度。

当我们不以测速为目的，而仅仅以确定各时元动态点位为目的时，可以不顾及机动加速度的影响，即可以不采用机动识辩动态模型而采用常速模型并假定较大的动态噪声将机动加速度引起的速度扰动归入动态噪声之中。由于载波相位具有很高的量测精度，此时滤波位置估计仍具有很高的精度，且滤波不会发散。在多种机载GPS动态定位数据处理中发现，当选择Q≥0.4，L1、L2载波相位量测精度分别假定为0.025周时，各次滤波结果中坐标分量相差仅在毫米级，图7.2.1为一次秒采样率机载GPS动态定位中，QK中元素Q分别选择为0.4、1、4、10 时滤波结果中X坐标分量的互差。

基本以上内外可靠性分析以及在已知整周模糊度已知时载波相位观测值可转换为高精度测相位伪距的特点，提出了一种简单实用的高精度GPS动态定位滤波处理模型。由于状态定位向量中不包含模糊度参数，称之为无模糊度参数滤波模型。

无模糊度参数滤波模型由式（35）定义：

图1 不同QK下滤波差值

系统动态噪声和量测噪声的统计特性假设如下：

式（35a）中动态模型采用常速模型。状态转移矩阵ΦK，K-1和干扰矩阵TK-1定义如下：

式中，I3为3×3维单位阵；ΔT为采样时间间隔。其中，状态向量为；xK，yK，zK为动态接收机天线在WGS-84地心坐标系中的三维位置，为动态接收天线的三维速度。在该模型中选取较大的模型噪声QK，根据内部可靠性分析可探测的机动加速度下界值也大，从而相应的局部检验量对机动加速度不敏感。根据外部可靠性分析，机动加速度对滤波位置影响很小，而动态定位所需要的也只是各观测时元的动态位置，因而选择较大的QK不会因机动加速度降低动态定位的精度。

在此模型下，关键的问题是周跳的探测与修复。周跳的探测和修复按如下步骤进行：

第一步：根据预处理中周跳探测结果，首先删除有周跳标记卫星所对应的双差载波相位观测值。

第二步：计算预报残差及其协方差阵，然后采用式（21）进行局部示警检验。当＜F时，认为剩余双差观测值没有周跳，转到第五步。当≥F时，则剩余双差观测值中仍存在周跳，但无法确定在哪个观测值上，为此需进行进一步检测。

第三步：对所剩余双差观测值分别采用式（24b）进行局部一维定位检测，并删除数值最大检验量所对应双差观测值。

第四步：重复第二、三步直到局部示警检验量TG1＜F。

第五步：利用剩余双差观测值进行滤波计算。

第六步：检查周跳标记，并根据滤波位置反求有周跳标记双差载波相位观测值的周跳数。

以上周跳的探测的修复基本上采用卢刚和魏明的周跳探测、修复方法，并基于内外部可靠性研究抛弃了并无多大作用而耗力的机动加速度检测部分，简化了处理模型。研究表明以上周跳探测和修复方法对于少量周跳（无周跳的双差载波相位观测值数量≥3）是十分有效的，若周跳卫星很多，连续跟踪的卫星数＜4时，以上方法则无法完全修复周跳，此时需利用在航模糊解算方法重新求各卫星整周模糊度，这和卢刚的研究结果也是一致的。

5 数值分析

5.1 周跳的探测与修复

根据上文所述模型，对无周跳数附加虚拟跳值后进行了周跳的探测和修复，所用数据为1996年8 月4日哈尔滨进行的一次秒采样率机载GPS动态定位（表8，表9）。表8中数据为飞机静止时的周跳探测与修复结果；表9中数据为飞机起飞后的周跳探测和修复结果。由于飞机在运动，因而各双关载波相位的预报残差较大，但此时仍有很好的周跳探测和修复能力，充分证明了上文模型的有效性。

5.2 与GPSur vey软件解算结果比较

Ttimble公司为其测量型接收机配备了高精度定位数据处理软件包GPSurvey，它可以处理静态和动态数据，为了验证本文算法的正确性，选择了一套动态定位数据分别采用文方法和GPSurvey软件进行了处理，其差值见图2（图中虚线表示基准站与流动站的距离，实线为解算差值）。从该图中可以看出，即使距离远达100km，本文方法与GPSurvey软件解算差值也不超过5cm，即不超过0.5ppm，这也证明了本文方法的正确性。

图2 无模糊度参数解算与GPSurvey软件解算成果比较

[1] 刘基余.GPS卫星导航定位原理与方法（第二版）.北京：北京科学出版社，2008.6.

[2] Lachapelle，G.，GPS Theory and Applications，University of Calgary，Fall 2000，PP.310.

[3] 刘基余，陈小明，李德仁等.GPS Kinematic Carrier Phase Measurements for Aerial Photogrammetry，ISPRS Journal of Photogrammetry and Remote Sensing，Vol.51，No.5，1996，P.230～242.