一种基于卡方检验的AWPSO姿态解算方法与流程

一种基于卡方检验的awpso姿态解算方法
技术领域
1.本发明涉及全球卫星导航技术领域，具体涉及一种基于卡方检验的awpso 姿态解算方法。


背景技术：

2.姿态测量是现代导航系统的基础，基于卫星载波相位的姿态测量具有低成 本，高精度，无惯性元件的累积误差等优点。
3.目前基于载波相位的姿态测量可分为姿态角直接解算法和姿态矩阵估计法。 姿态角直接解算法利用已知的天线位置与旋转姿态角的关系直接解算姿态角， 包括使一条基线设置为载体的主轴方向，求解出航向角和俯仰角，利用第二条 基线的旋转关系求解出横滚角的两天线测姿法和利用两天线测姿公式确定航向 角和俯仰角，将另外天线经两次转动得到横滚角的多天线测姿法。整周模糊度 的精度决定了直接解算的角度精度。常见的算法有teunissen提出的基于模糊度 空间的lambda算法和基于位置域的afm算法；姿态矩阵估计法利用多历元 多颗卫星建立载波相位观测方程求解基线坐标，通过基线矢量在不同坐标系下 的转换关系估计出姿态矩阵，进而求得姿态角度。姿态矩阵的估计常采用最小 二乘法进行估计，姿态求解过程中包括了基线坐标估计误差和姿态矩阵的转换 估计误差，大大降低了姿态角的解算精度。
4.然而，传统姿态角直接解算法中存在整周模糊度求解复杂，解算时间长的 问题。


技术实现要素：

5.本发明的目的在于提供一种基于卡方检验的awpso姿态解算方法，旨在 避免直接求解整周模糊度引起的复杂计算过程，同时提高姿态解算精度，降低 系统计算复杂度，提升解算实时性和鲁棒性。
6.为实现上述目的，本发明提供了一种基于卡方检验的awpso姿态解算方 法，包括下列步骤：
7.载入卫星原始观测量数据；
8.建立双差载波相位模型；
9.将基线矢量与姿态角关系引入所述双差载波相位模型；
10.建立姿态角的适应度函数；
11.改进粒子群算法，自适应更新粒子的惯性权重和学习因子；
12.当种群收敛超过30代，保存当前全局最优解至候选解序列，在迭代次数小 于总迭代次数时，重新初始化粒子位置及权重，重新搜索；
13.对候选解序列计算马氏距离，并进行卡方检验，剔除大误差解；
14.对通过检验的候选解进行加权计算，求解最终姿态角度值。
15.其中，所述双差载波相位模型由主天线a、从天线b、卫星i和卫星j组成， 其中所述主天线a和所述从天线b为固联在载体上的两个gps天线，所述卫星 i和所述卫星j到两个天
线的距离大于基线长度，同一颗卫星在所述主天线a和 所述从天线b处的高度角和方位角相同。
16.其中，在将基线矢量与姿态角关系引入所述双差载波相位模型的过程中， 通过基线矢量间的关系获得所述主天线a和所述从天线b分别对应卫星的单差 载波相位观测方程，再对同一时刻两颗卫星的单差载波相位观测方程求差，消 除天线的瞬时钟差，获得双差载波相位观测方程。
17.其中，在改进粒子群算法，自适应更新粒子的惯性权重和学习因子的过程 中，首先初始化粒子搜索空间，计算每个粒子的适应度值，对最佳适应度位置 进行更新，同时自适应更新粒子的惯性权重和学习因子。
18.其中，所述全局最优解由惯性权重和学习因子决定，通过重新初始化粒子 种群使粒子群跳出局部最优，并经过对惯性权重的迭代公式修改为自适应变化 使搜索空间随迭代动态减小，最终使粒子动态地向全局最优收敛。
19.其中，所述候选解序列的建立，包括下列步骤：
20.步骤1：在搜索空间中随机生成m个n维粒子xi；
21.步骤2：计算当前每个粒子的适应度值，并与粒子个体当前的最优位置p
best
比 较，若粒子当前的适应度值大于p
best
，则将p
best
更新为当前适应度值，否则不更 新；
22.步骤3：通过适应度函数计算粒子群的适应度值，得到当前的全局最优解g
best
23.步骤4：计算当前惯性权重和学习因子，同时，更新粒子的速度，位置：
24.vi＝w*vi+c1*rand*(p
i-xi)+c2*rand*(g
i-xi)
25.xi＝xi+vi26.ci＝c
max
+(c
min-c
max
)*i/t
27.w＝(w
max-w
min
)*i/t
28.其中xi＝(x
i1
,x
i2
,x
i3
,...,x
in
)为粒子当前位置，vi＝(v
i1
,v
i2
,v
i3
,...,v
in
)为粒子当前的速度， w为自适应惯性权重，ci为自适应学习因子，g为候选解序列，c
max
和c
min
分别为 最大和最小学习因子，w
max
和w
min
分别为最大和最小惯性权重，t为迭代次数；
29.步骤5：当粒子群中到目前为止的全局最优位置g
best
收敛保持30代，且迭代 次数小于预设总迭代次数，将当前适应度值及相应的粒子速度，位置保存至序 列g中，同时跳出收敛，重新初始化xi及自适应惯性权重w，重回步骤2：
30.xi＝rand*(x
max-x
min
)+xi31.w＝w-(w
max-w
min
)*i/t
32.步骤6：判断当前累积迭代次数是否大于总迭代次数，如果没有则返回步骤 2，否则迭代结束。
33.其中，所述卡方校验具体为服从正态分布的候选解其马氏距离的平方服从 卡方分布。
34.其中，对候选解序列计算马氏距离，并进行卡方检验，剔除大误差解过程， 具体为当累计迭代大于预设值时，退出迭代，对候选解序列进行马氏距离计算， 并对计算结果进行卡方校验。
35.其中，在对通过检验的候选解进行加权计算，求解最终姿态角度值过程中， 首先对不符合置信区间内的搜索结果，视为无效值将其剔除，再通过校验的候 选解进行加权运
算，求解出最终姿态角。
36.本发明提供了一种基于卡方检验的awpso姿态解算方法，利用姿态角与 基线矢量的关系，构造姿态角的单历元适应度函数模型，通过适应度函数进行 姿态角搜索，规避了整周模糊度求解带来的大运算量，同时采用自适应权重因 子和候选解的方式提高了姿态角的搜索范围，减少算法早熟，提高算法收敛性， 最后利用卡方检验对候选解进行有效性检测，并通过剔除无效解减小局部最优 解引起的误差，提高了解算的精度和鲁棒性。
附图说明
37.为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施 例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述 中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付 出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
38.图1是本发明的一种基于卡方检验的awpso姿态解算方法的流程示意图。
39.图2是本发明的具体实施例中基线与卫星i和j的空间几何关系图。
40.图3是本发明的具体实施例的首历元适应度函数三维图。
41.图4是本发明的具体实施例的姿态角(航向角，俯仰角)分布直方图。
42.图5是本发明的具体实施例中单历元候选解马氏距离卡方检验分布图。
43.图6是本发明的具体实施例中的适应度函数进化曲线图。
44.图7是本发明的具体实施例中不同算法分组对比实验中直接法姿态角误差示 意图。
45.图8是本发明的具体实施例中不同算法分组对比实验中最小二乘法姿态角误 差示意图。
46.图9是本发明的具体实施例中不同算法分组对比实验中传统pso算法姿态角 误差示意图。
47.图10是本发明的具体实施例中不同算法分组对比实验中本发明的姿态角误 差示意图。
具体实施方式
48.下面详细描述本发明的实施例，所述实施例的示例在附图中示出，其中自 始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。 下面通过参考附图描述的实施例是示例性的，旨在用于解释本发明，而不能理 解为对本发明的限制。
49.为了使本领域技术人员对本发明的技术方案更加清晰的理解，先对本发明 中的一些技术术语进行简要说明。
50.pso(particle swarm optimization)是粒子群优化算法的英文缩写，是一种 基于种群的随机优化技术，由eberhart和kennedy于1995年提出。粒子群算法 模仿昆虫、兽群、鸟群和鱼群等的群集行为，这些群体按照一种合作的方式寻 找食物，群体中的每个成员通过学习它自身的经验和其他成员的经验来不断改 变其搜索模式。
51.awpso(adaptive weight particle swarm optimization)是自适应粒子群优 化算法的英文缩写，在pso算法的基础上进行改进，通过自适应变化的学习因 子，惯性权重提
高算法的收敛性能。
52.卡方检验就是统计样本的实际观测值与理论推断值之间的偏离程度，实际 观测值与理论推断值之间的偏离程度就决定卡方值的大小，如果卡方值越大， 二者偏差程度越大；反之，二者偏差越小；若两个值完全相等时，卡方值就为0， 表明理论值完全符合。本发明通过卡方检验来检验解算结果与实际值的偏差程 度，予以校正，提高精度。
53.请参阅图1，本发明提出了一种基于卡方检验的awpso姿态解算方法，包 括下列步骤：
54.s1：载入卫星原始观测量数据；
55.s2：建立双差载波相位模型；
56.s3：将基线矢量与姿态角关系引入所述双差载波相位模型；
57.s4：建立姿态角的适应度函数；
58.s5：改进粒子群算法，自适应更新粒子的惯性权重和学习因子；
59.s6：当种群收敛超过30代，保存当前全局最优解至候选解序列，在迭代次 数小于总迭代次数时，重新初始化粒子位置及权重，重新搜索；
60.s7：对候选解序列计算马氏距离，并进行卡方检验，剔除大误差解；
61.s8：对通过检验的候选解进行加权计算，求解最终姿态角度值。
62.以下通过具体实施例结合实施步骤作进一步说明：
63.所述双差载波相位模型的建立：设固联在载体上的两gps天线分别为主天 线a和从天线b，两天线组成长为b的基线矢量b。基线与卫星i和j的空间几何 关系如图2所示。
64.其中，为天线a和b对卫星i的单差载波相位，为天线a和b对卫星 j的单差载波相位，p和y为基线矢量的俯仰角和航向角；α和ω为卫星j对基线 所处载波平面的高度角和方位角。由此可得天线a和b对卫星i单差载波相位观 测方程：
[0065][0066][0067]
式中，λ为载波长度；为卫星i到天线a的单位视线向量；为b与的夹角；c为光速；δt为天线a和b间的瞬时钟差；为卫星i到天线a和b的 单差整周模糊度；δi为观测噪声。
[0068]
对同一时刻两颗卫星的单差载波相位观测方程求差，消除天线的瞬时钟差， 可以得到如下双差载波相位观测方程：
[0069][0070]
由于卫星到两天线的距离远大于基线长度，载波信号可视为平面波，因此 同颗卫星在天线a与b处的高度角和方位角相同。将姿态角代入载波相位方程可 得:
[0071][0072]
式中为双差整周模糊度，远大于δ
ij
，因此即：
[0073][0074]
当p和y取值为真实航向角和俯仰角时，计算所得值到达极大值。此时，姿 态角的
求解转为给定区间的最优值搜索问题。选取最大仰角卫星l为参考卫星， 根据上式建立姿态角的适应度函数：
[0075][0076]
进一步的，粒子群的全局最优解由惯性权重和学习因子决定，而不同的模 型不同的应用场景，最佳的惯性权重和学习因子不尽相同，故标准pso容易陷 入局部最优解，难以搜索到精确的姿态角。具体如图3所示，姿态角适应度函 数在搜索区间内存在多个局部最优解，在搜索过程中需要排除局部最优解对全 局最优解的干扰，使用传统pso算法所得解会极易陷入局部最优。为此，需要 对算法进行改进。首先，pso搜索迭代收敛超过30代时，粒子可能陷入局部最 优的困境，通过下式重新初始化粒子种群，根据当前粒子位置生成新的粒子群， 使得粒子群能够跳出局部最优：
[0077]
xi＝rand*(x
max-x
min
)+xi[0078]
同时对惯性权重的迭代公式修改为自适应变化，使搜索空间随迭代动态减 小。
[0079]
w＝w-(w
max-w
min
)*i/t
[0080]
上述修改能使得粒子动态地向全局最优收敛。最后，我们在算法中引入了 候选解序列，保存每次收敛时的g
best
。并通过一定检验方法，进行有效性检验。
[0081]
awpso建立候选解序列的基本步骤如下：
[0082]
(1)在搜索空间中随机生成m个n维粒子xi[0083]
(2)计算当前每个粒子的适应度值，并与p
best
(p
best
为粒子个体当前的最优 位置)比较，若粒子当前的适应度值大于p
best
，则将p
best
更新为当前适应度值， 否则不更新
[0084]
(3)通过适应度函数计算粒子群的适应度值，得到当前的全局最优解g
best
[0085]
(4)计算当前惯性权重和学习因子，同时，更新粒子的速度，位置。
[0086]
vi＝w*vi+c1*rand*(p
i-xi)+c2*rand*(g
i-xi)
[0087]
xi＝xi+vi[0088]
ci＝c
max
+(c
min-c
max
)*i/t
[0089]
w＝(w
max-w
min
)*i/t
[0090]
其中xi＝(x
i1
,x
i2
,x
i3
,...,x
in
)为粒子当前位置，vi＝(v
i1
,v
i2
,v
i3
,...,v
in
)为粒子当前的速度， w为自适应惯性权重，ci为自适应学习因子，g为候选解序列，c
max
和c
min
分别为 最大和最小学习因子，w
max
和w
min
分别为最大和最小惯性权重，t为迭代次数；
[0091]
(5)当粒子群中到目前为止的全局最优位置g
best
收敛保持30代，且迭代次 数小于预设总迭代次数，将当前适应度值及相应的粒子速度，位置保存至序列g 中。同时跳出收敛，重新初始化xi及自适应惯性权重w，重回步骤(2)：
[0092]
xi＝rand*(x
max-x
min
)+xi[0093]
w＝w-(w
max-w
min
)*i/t
[0094]
(6)判断当前累积迭代次数是否大于总迭代次数，如果没有则返回步骤(2)， 否则迭代结束。
[0095]
图4是本发明的具体实施例的姿态角(航向角，俯仰角)分布直方图。
[0096]
从图4可知，awpso算法搜索结果呈正态分布状。利用此特性，可采用卡 方检验进行候选解的筛选。
[0097]
图5是本发明的具体实施例中单历元候选解马氏距离卡方检验分布图。
[0098]
由图5可知，通过本算法搜索所得解的马氏距离的卡方检验结果，超过95％ 的解的显著水平小于2.706，因此可取显著水平为2.706为检验阈值。
[0099]
进一步的，基于卡方检验的候选解马氏距离内容如下：
[0100]
卡方检验原理是服从正态分布的候选解其马氏距离的平方应服从卡方分布。 候选解gi的马氏距离为：
[0101][0102]
其中，μ为候选解的均值，s为候选解的协方差矩阵。d的平方服从卡方分 布：
[0103]
d2＝(g
i-μ)
t
s-1
(g
i-μ)～χ2[0104]
根据假设检验方法，设awpso搜索所得姿态角没有陷入局部最优，即假 设候选解服从正态分布。设未陷入局部最优的概率为α，则：
[0105][0106]
式中，为显著水平为α的卡方分布，t为检测阈值，显著性水平与阈值关 系如表1所示，α越小，则阈值越大，判定条件越宽松，则陷入局部的解与真实 值的误差越大。当d2＜t时，认为候选解未陷入局部最优，反之，则认为候选解 陷入局部最优。
[0107]
若判定候选解陷入局部最优时，可以将当前解视为无效值，将其剔除。经 卡方检验筛选后的候选解序列可认为最接近真实姿态角的全局最优解，对该候 选解序列做加权运算，即可得到最优姿态角。
[0108][0109]
表1性能分析
[0110][0111]
图6是本发明的具体实施例中的适应度函数进化曲线图。
[0112]
基于卡方检验的awpso算法历元间的适应度值曲线如图6所示，在总共 1000历元里，每历元的适应度值始终保持在0.8以上，说明算法解算出的角度 为最接近真实值的全局最优解。适应度值的波动范围不超过0.1，有效证明算法 鲁棒性强。
[0113]
本发明还就姿态角误差与传统姿态解算方法进行了对比，请参阅图7至图 10。
[0114]
由图7、8、9、10可知，四种算法中，直接法所求的姿态角误差最大，误差 值最大可达5
°
，最小二乘法误差小于直接法，最大不超过1.5
°
，pso算法误 差略小于最小二乘法，误差在
±1°
以内，awpso算法误差最小，相较pso， 误差缩小近0.5
°
，历元间姿态角度解算稳定性高。
[0115]
各算法的rms误差精度由高到低分别为awpso，pso，最小二乘和直接法。 在直接法下，两个姿态角的rms误差分别高达1.688
°
和1.046
°
，但在awpso 算法下，姿态角的误差减小至0.245
°
和0.243
°
。航向角精度较直接法提高近 85％，俯仰角精度提高约77％。相较最小二乘和pso算法，awpso算法的航 向角精度提高约53％，俯仰角精度分别提高将近29％和53％。
[0116]
以上所揭露的仅为本发明一种较佳实施例而已，当然不能以此来限定本发 明之权利范围，本领域普通技术人员可以理解实现上述实施例的全部或部分流 程，并依本发明权利要求所作的等同变化，仍属于发明所涵盖的范围。