一种缩短整周模糊度求解时间的方法与流程

本发明涉及卫星导航定位领域中的整周模糊度求解方法，具体涉及一种缩短整周模糊度求解时间的方法。

背景技术

正确使用载波相位观测测量载体姿态的关键是实时、精确的固定载波相位的整周模糊度，只有精确的整周模糊度固定解才能为精密观测系统以及航天航空等制导系统提供精密服务成为可能；lambda算法是固定整周模糊度算法中较为常见的算法，但是lambda算法存在固有的缺陷；首先，模糊度固定过程中置换发生的次数与整个解算过程所消耗的时间成正比，而置换过程存在一定的冗余；即置换过程的计算量与计算时间占比整个置换过程较高；其次，对整周模糊度进行搜索求解时，始终用固定的超椭球体尺寸对搜索空间进行限制，这样每次的搜索都重复搜索同一个空间内的点，这些点也包括不适合做模糊度的点；但是单更需要得到两个及以上的最小二乘估计的最佳值时，这就可能造成在对其他的最小二乘值进行求解的搜索空间相对第一次求解时的搜索空间较“臃肿”的现象。

技术实现要素：

本发明提供一种针对现有lambda算法中存在的不足，在保证整周模糊度解算正确率的条件下可缩短整周模糊度解算时间的方法。

本发明采用的技术方案是：一种缩短整周模糊度求解时间的方法，包括以下步骤：

步骤1：建立双差方程，采用最小二乘法获得整周模糊度的浮点解a及其协方差矩阵

步骤2：根据协方差矩阵构造转换矩阵z，对浮点解和协方差矩阵进行转换；

步骤3：根据转换后浮点解和协方差矩阵进行整周模糊度的搜索得到整周解；

步骤4：对整周解进行z的反变换即可得到整周模糊度解；

步骤2中通过z变换对协方差矩阵进行变换分解过程中，将对角矩阵中最小的元素变换到n的位置上，然后求出相关的变换矩阵，最后采用相关矩阵更新三角矩阵和对角矩阵，具体过程如下：

s1：协方差矩阵的变换形式如下：

式中：p为对角矩阵中最小元素变换到n位置得到的变换矩阵，p1为第一次变换后得到的矩阵，u为三角矩阵，u^t为u的转置矩阵，qn为中的最小对角元素；为为(n-1)*(n-1)维矩阵，为1*(n-1)维向量，为(n-1)维向量，d为对角矩阵；

其中：

式中：d为n处的对角元素，d＝qn；l为三角矩阵，o为0矩阵，为(n-1)*(n-1)矩阵，为(n-1)*(n-1)矩阵，为单位矩阵；

根据式(1)得到变换矩阵p1；

s2；重复步骤s1，即可得到整个变换矩阵：

式中：pi为变换矩阵，i＝1，2，...，n-1；

根据变换矩阵得到更新后的对角矩阵d；

s3：当更新后的对角矩阵d中位置k满足以下条件时，进行置换；如果满足条件的k不存在，则不进行置换；

进一步的，所述步骤3的求解公式如下：

式中：为浮点解向量经过z变换得到的向量，z为，为协方差矩阵经z变换得到的向量，χ²为搜索空间尺寸；

搜索空间式为：

根据式(5)求得其中一组解为z1，将其重新带入式(5)得到s(z1)；

将求得的s(z1)值替换式(4)中的χ²；

重复上述求解步骤，得到所有的整周解。

进一步的，所述整周解求解过程如下：

s11：假定χ²的值趋于正无穷，根据式(4)和(5)得到第一个整数矢量式中为；

s12：根据步骤s11得到的整数矢量，取z2中的第一个元素和z1中与第一个元素相近的整数，剩余元素取z1向量中的元素；重复此步骤得到第p个整数矢量zp；

s13：将步骤s12中得到的整数矢量分别带入式(5)，得到关于s函数的解向量s(z)＝[s(z1),s(z2),...,s(zp)]；

s14：取χ²＝s(zp)，基于此进行空间搜索，重复步骤s11-s14，得到p组整周解。

本发明的有益效果是：

(1)本发明通过采用路径最优算法，减少了置换次数，整个整周模糊度的求解时间进一步减少；

(2)本发明通过对searching过程进行优化，模糊度的搜索空间可以在每一个最佳的模糊度求解过程中都得到缩减，使得整周解在搜索空间上的搜索效率得到提升；

(3)本发明在保证解算整周模糊度解的准确度的条件下，缩短了解算时间，减少变换冗余，提升了搜索效率。

附图说明

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

图2为本发明实施例中采用本发明中对reduction过程优化后与现有lambda算法解算时间对比图。

图3为本发明实施例中采用本发明中对searching过程优化后与现有lambda算法解算时间对比图。

具体实施方式

下面结合附图和具体实施例对本发明做进一步说明。

本发明针对现有lambda存在的不足基于整周模糊度的浮点解及协方差矩阵，在去相关变换搜索求解的过程中，对reduction和searching两个过程分别进行优化；基于reduction过程的优化可以进一步划分为置换过程路径化和路径最优化两个基本步骤，而基于searching过程的优化目的在于减少搜索空间的冗余。

lambda算法是一种基于目标函数的整周模糊度的搜索解算算法，主要流程包括：

步骤1：建立双差方程，采用最小二乘法获得整周模糊度的浮点解a及其协方差矩阵双差方程如下：

y＝b·b+a·a+e(6)

步骤2：根据协方差矩阵构造转换矩阵z，对浮点解和协方差矩阵进行转换；利用z变换将浮点解矩阵和协方差矩阵变换成和达到去相关的目的；

步骤3：根据转换后浮点解和协方差矩阵进行整周模糊度的搜索得到整周解；求解公式如式(4)，求得解记为n；

步骤4：对整周解进行z的反变换即可得到整周模糊度解；一般将步骤2称为reduction过程，将步骤3记为searching过程。

reduction过程具体实现

1、z变换分解协方差矩阵

reduction过程通过z变换对协方差矩阵进行变换达到去相关效果；其具体实现是，通过z变换，将模糊度的浮点解向量及其相应的协方差矩阵转换为新的模糊度向量和对应的协方差矩阵如式4所示。

对浮点解和协方差矩阵进行z变换后，根据式(4)进行整周模糊度的搜索求解；下面将阐述z变换矩阵的具体变换过程。

将协方差矩阵进行cholesky分解可得其三角矩阵与对角矩阵相乘的形式，并将其代入式(7)可得如式(8)的表达式。

式(8)中，d为对角矩阵，l^-t单位三角矩阵。

在以上的变换过程中，需要满足以下两个条件：

1)变换后的协方差系数尽可能小，即单位三角矩阵l的非对角线元素尽可能的小，以此达到去相关性；

2)对角矩阵d的对角元素在矩阵中从(1，1)至(n，n)的排列位置应按照从大到小的顺序排列，即

以上进行z变换并满足所诉的两个条件的过程，是通过迭代的整数高斯变换(以下简称高斯变换)过程和置换过程完成的，高斯变换在一定程度上降低相关性，但此时并未解决频谱不连续的问题，因此需要通过置换过程满足条件2)以使得解决不连续的问题，提高模糊度搜索的效率。

对cholesky分解得到的l和d进行分块表示，可以分别表示成如式(9)所示的形式。

在式(9)中，l11是i-1维矩阵，l22是二阶矩阵，l33表示n-i-1维矩阵，其他的分块矩阵可分别确定阶数；则对应的第i和i+1维的z变换矩阵可以表示如式(10)的形式。

式(10)中，z22表示二阶满秩方阵；根据式(8)则可得到新的关于z的分解形式如式(11)所示。

此外，根据式(7)和式(10)，可以将变换后的协方差矩阵表示为式(12)的形式。

从式(11)可以看出，更新后的分解矩阵和中，只有和在变化过程中发生了变化，则结合式(10)，式(7)以及式(12)可得到式(13)表示的变换更新公式。

从式(12)可以看出，式(7)表示的协方差矩阵的变换形式可以由原协方差矩阵经过cholesky分解后的分解矩阵得到；此外，通过不断地更新单位下三角矩阵l和对角矩阵d可以重复的进行二维变换，并以此实现对原协方差矩阵的变换。

2、高斯变换去相关

在式(10)中，将转换矩阵限制为单高斯转换矩阵，则可将其转变为式(14)所示的形式。

则由式(13)可得只发生了变化的单位三角矩阵l和对角矩阵d的子矩阵的更新形式分别如式(15)和(16)所示。

由式(15)和式(16)明显可以看出，高斯变换过后，只有对角矩阵l发送变化，而对角矩阵d并未发生变化，即协方差和条件方差会随着高斯变化的过程而变化，但方差即对角矩阵不会改变；通常为了使去相关的效果最好，一般会使α取离li+1,i最近的整数的相反数的数值，这样，就可以使l中的任一个非对角元素均小于0.5。

3、置换重排序过程

经过上诉高斯变换后，reduction过程满足了条件一，但此时可能会因为频谱之间存在连续性等的问题而造成模糊度搜索效率低的情况，而此时对模糊度进行置换重排列就可以解决这些问题。

如式(13)，将置换矩阵替代z22，从而可以使得第i个和第i+1个模糊度置换，具体如式(17)所示。

则可分别解得发生变化后的三角矩阵和对角矩阵的子矩阵如式(18)和式(19)所示；

对式(18)和式(19)进行简化，则可得到如式(20)的形式。

当第i个对角元素比第i+1个对角元素小时发生置换；将式(20)，式(19)与式(11)和式(10)对比可知，经过式(20)的置换，三角矩阵l和对角矩阵d都得到了更新，第i和第i+1模糊度也进行了重排列。

searching过程解析

采用最小二乘估计法求解整周模糊度等价于基于浮点解和协方差矩阵求解式(21)的最小化问题。

当采用lambda算法固定整周模糊度时，全局的搜索空间zⁿ由局部的搜索空间代替，即给式(21)限定搜索空间范围，如式(22)所示。

式(21)表示的是一个以a为中心的超椭球体，其形状由协方差矩阵控制，球体尺寸由χ²控制，球体形状和尺寸对整周模糊度的搜索效率有着很大的影响；此时的超椭球体狭长，搜索效率较低，为了提高搜索效率，对协方差矩阵和浮点解进行reduction过程后，式(22)转换为式(23)的形式。

式(23)中，为了保持整周模糊度的正确性，仅对椭球体的形状做了改变，其大小尺寸并不发生改变；将式(8)代入(23)可得式(24)的展开形式。

在式(24)中，

结合式(24)以及(25)，可得到式(26)所示的分步搜索环的形式。

搜索空间尺寸χ²的选择

由于reduction过程解决了椭球空间形状狭长的问题，因此，在searching过程中，实际控制搜索空间大小的li+1,i最影响搜索效率；在搜索过程中，不能选择太大的li+1,i，否则会在搜索空间中出现不必要的网格点，影响搜索效率，但如果选择较小的li+1,i又可能造成整周模糊度丢失的现象。

在lambda算法的搜索过程中，χ²的选取与整周模糊度度解的候选个数有关，通常可以通过超椭球空间的容积和计算数值的方法来选取χ²，下面就简单介绍如何通过超椭球空间的容积选取模糊度空间的尺寸。

在超椭球空间中，其容积可用如式(27)表示。

在式(27)中，vn表示单位超球体的容积，可以通过迭代的方式求得；而表示方差阵的行列式；在搜索过程中，候选整周模糊度与容积之间的关系如下式。

p＝[en](28)

式(28)中，p表示候选个数，[*]表示不超过*的离得最近的整数；根据式(28)和式(27)两式即可确定χ²的大小；利用这种方法确定的搜索空间尺寸仅和候选整周模糊度相关，这样就避免了搜索尺寸不合适造成的搜索效率低或缺失整周模糊度的问题。

本发明针对lambda算法基于reduction过程和searching两个过程进行优化，其中基于reduction过程的优化可进一步划分为置换过程路径化和路径最优化两个基本步骤；基于searching过程的优化目的在于减少搜索空间的冗余。

基于reduction过程的优化：

使用路径最优化算法减少置换次数，这个过程进一步分为置换过程路径化和路径最优化两个过程；置换过程路径化指的是：在reduction过程中用z变换对协方差矩阵进行变换分解时，将每一步变换都以“路径”的形式进行转换；首先将对角矩阵中最小的元素变换到n的位置上，然后求出相关的变换矩阵p，最后用相关的变换矩阵更新三角矩阵l和对角矩阵d。

具体过程如下：

s1：假定在协方差矩阵的变换中，将对角矩阵的最小元素置换到位置n处，其变换矩阵用p1表示，则对协方差矩阵的变换形式如下：

式中：p为对角矩阵中最小元素变换到n位置得到的变换矩阵，p1为第一次变换后得到的矩阵，u为三角矩阵，u^t为u的转置矩阵，qn为中的最小对角元素；为(n-1)*(n-1)维矩阵，为1*(n-1)维向量，为(n-1)维向量，d为对角矩阵；

其中：

式中：d为n处的对角元素，d＝qn；l为三角矩阵，o为0矩阵，为(n-1)*(n-1)矩阵，为(n-1)*(n-1)矩阵，为单位矩阵；

根据式(1)得到变换矩阵p1。

s2；进一步对以同样的方式进行变换，假定其变换矩阵为p2，重复步骤s1，直到上述变换至n-1处变换停止，且其变换矩阵为pn-1即可得到整个变换矩阵：

式中：pi为变换矩阵，i＝1，2，...，n-1。

在置换过程路径化的过程中，逐个求出各个“路径化”的变换矩阵pi为变换矩阵，i＝1，2，...，n-1，且按对角矩阵元素的从大到小的顺序进行路径化，因此变换后的对角矩阵符合从小到大的顺序排列。

而路径最优化过程则指，每一次置换都选择使对角矩阵产生变化最大的位置进行变换，以此达到最优化；由于每一次置换都选择使对角矩阵上的变化率最大的位置进行置换，很大程度上减少了置换次数。

s3：当更新后的对角矩阵d中位置k满足以下条件时，进行置换；如果满足条件的k不存在，则不进行置换；

在对置换过程完成路径化之后，具体的置换将不再按照从右到左(i＝n,n-1,...,1)的顺序置换，而是选择路径化后对角矩阵d中的对角元素dk产生变化最大的位置进行置换，即位置k应当满足的条件如式(3)所示：

如果满足式(3)放入k不存在，则说明路径已最优，不需要进行置换；由于置换次数得到了相应的减少，整个整周模糊度的求解时间也将得到进一步减少。

基于searching过程的优化算法：

步骤3的求解公式如下：

式中：为浮点解向量经过z变换得到的向量，z为，为协方差矩阵经z变换得到的向量，χ²为搜索空间尺寸；

假定z变化后的搜索空间式如下所示：

当需要求解n个(n大于或等于2)整周模糊度时，假定根据式(5)求得其中一组解为z1，将其重新带入式(5)得到s(z1)；

将求得的s(z1)值替换式(4)中的χ²；此时搜索空间将比第一次搜索模糊度解时的搜索空间更小，在固定模糊度求解的过程中，不断地以这种方式进行求解，模糊度的搜索空间可以在每一个最佳的模糊度求解的过程中都得到缩减。

具体实现过程如下：

s11：假定χ²的值趋于正无穷，根据式(4)和(5)得到第一个整数矢量式中为；

s12：根据步骤s11得到的整数矢量，取z2中的第一个元素与z1中的第一个元素相近的整数，剩余的元素取与z1向量中剩余的元素；重复此步骤得到第p个整数矢量zp

s13：将步骤s12中得到的整数矢量分别带入式(5)，得到关于s函数的解向量s(z)＝[s(z1),s(z2),...,s(zp)]；

s14：为了避免漏掉整数解的情况，取χ²＝s(zp)(此时的χ²相对于步骤s11中已经得到了相应的缩减)，基于此搜索空间再次进行搜索，重复步骤s11-s14，得到p组整周解；在这个过程中，搜索空间每一次的搜索都得到相应的缩减，这就使得下一个整周解在该搜索空间上的搜索效率得到了提升。

为了验证本发明方法的有效性，对本发明方法进行仿真实验并与现有lambda算法进行对比。

采用的操作系统是windows10，软件环境是matlab-2014a，硬件环境是i7-4790，3.6ghz。此外，姿态测量的实验基于后处理的方式，gps的观测数据来源于已经收集到的gps观测数据，这些数据均存储在rinex格式文件中。

通过仿真手段构建lambda算法及其优化算法所需的浮点解和协方差矩阵，并基于此构建算法的搜索解算空间；下文中or-过程算法为在lambda算法的基础上对reduction过程进行优化后的算法；os-过程算法为在lambda算法的基础上对searching过程进行优化后的算法。

协方差矩阵可通过cholesky分解成三角矩阵和对角矩阵相乘的形式，l是单位下三角矩阵，对于对角线以下的矩阵元素可通过matlab中自带的rand函数实现；d是对角矩阵，在实验中通过diag函数实现，对角元素由rand函数产生；此外，为了让仿真实验更具合理性，分别构建了2、4、6、8、10、12、14、16、18、20、22、24维共12组数据，且分别对每组数据进行了20次实验。

or-过程算法的优化仿真

将准备好的数据分别通过lambda算法与or-过程算法进行解算，并分别将这两个算法求解20次这12组数据的平均时间汇成对比结果图，其结果如图2所示；横轴表示的是整周模糊度的维度，纵轴表示的是各维度分别基于lambda算法和or-过程算法20次实验的平均解算时间；从图中的效果来看，lambda算法求解整周模糊度解的时间比or-过程算法所需的时间更多，表明基于reduction过程优化后的算法在一定程度上减少了变换冗余，提高了解算效率。

表1表示的是随机抽取20次实验中的3组基于16维浮点解用lambda算法和or-过程算法解算整周模糊度固定解的结果对比表格。

表1.lambda和or-过程算法的解算结果对比表格

表格数据记录部分总共有16行，表示16维数据，行内的大括号里有三列，表示20次实验里随机抽取的三次实验；组合的意思即为随机抽取的三次实验数据，每次数据都是16维整周模糊度浮点解向量；而“lambda算法结果”和“or-过程算法结果”关键字对应的数据则分别是对这三次随机抽取的16维数据的整数解。

从表1可以看出，在整周模糊度解维度为16维的情况下，lambda算法和or-过程解算的整周模糊度解的准确度并无差异；此外结合图2可以得出以下结论：在和lambda算法准确率相同的情况下，or-过程算法所需的解算时间比lambda算法所需的时间更少。

os-过程算法的优化仿真实验

针对os-过程算法的20组优化仿真实验中(取20组实验的时间平均值作为比较对象)用os-过程算法和lambda算法解算12组不同维度的整周模糊度浮点解，并分析二者解算整周模糊度的时间效率；如图3所示的是当lambda算法使用os-过程算法优化和使用lambda算法对整周模糊度进行解算的时间解算图。

其坐标轴含义与图2相同，由图3可以看出随着维度的增加，os-过程算法的运行时间比lambda算法运行的时间更少，这表明，基于searching过程的优化算法在一定程度上减少了空间冗余，提升了整周模糊度的搜索效率。

为了分析lambda算法和os-过程算法的解算结果对比，从20次实验中随机抽取了3组基于lambda算法解算和os-过程算法解算16维浮点解的实验结果，并绘制了如表2所示的结果对比表格。

表2.lambda和os-过程算法的解算结果对比表格

结合图3和表2可以看出，在整周模糊度固定解的解算正确率相同的情况下，用os-过程优化后的lambda算法可在一定程度上提升lambda算法的解算效率。

本发明方法能够在保证整周模糊度解算准确率的前提下减少一定的变换冗余，提高解算效率；并且searching过程的优化能提高整周模糊度的解算效率。