一种GPS载波相位周跳探测与修复的方法与流程

文档序号:11233166阅读:913来源:国知局
一种GPS载波相位周跳探测与修复的方法与流程

本发明属于时间频率传递的技术领域,尤其涉及一种利用小波奇异值分解和自适应指数平滑进行gps载波相位周跳探测与修复的方法。



背景技术:

时间频率的研究是基础科学研究中的一个重要分支,时间频率在科研、定位系统、电力系统、军事和国家安全等方面处于举足轻重的地位。世界各国大都建有自己的守时实验室,产生本国的原子时标,并参与国际比对。

中国计量科学研究院承担着向地方实验室发布标准时间与频率,并参与国际原子时比对的任务。因此,与国际原子时比对,以获得稳定度及准确度更高的本地原子时是至关重要的。中国计量科学研究院目前采用的时间比对方法为gps共视法,但共视法要求两地的实验室在同一时刻需观测到同一颗卫星,且持续较长时间。若两地实验室距离较远,无法观测到同一颗卫星,则需第三个实验室作为中继。gps载波相位法无需考虑观测卫星的问题,且时频传递精度更高。但gps载波相位时频传递中有一个重要的问题,即周跳。周跳是指由于gps接收机失锁,或接收载波由于其他原因被阻断,而产生间断的现象。周跳的发生对gps载波相位时频传递影响较大。故而周跳探测与修复是非常重要的。本课题研究gps载波相位周跳探测与修复方法,旨在保证时频传递的精度,并提高本地原子时的稳定度和准确度。



技术实现要素:

本发明要解决的技术问题是,提供一种利用小波奇异值分解和自适应指数平滑进行gps载波相位周跳探测与修复的方法,利用gps观测文件,采用小波奇异值分解方法,探测出gps载波相位传递中,整周跳变的现象以及跳变发生的历元,并采用自适应指数平滑法以修复周跳,保证gps载波相位时频传递的精度。

为解决上述问题,本发明采用如下的技术方案:

一种gps载波相位周跳探测与修复的方法,包括以下步骤:

第一步、利用matlab读取rinex型的gps观测文件,得到一天内观测到的卫星号、gps伪距观测值和相位观测值;

第二步、选择伪距观测值以及相位观测值,由于同一个历元,gps接收机所接收到的观测值来自多颗卫星,且跟踪同一颗卫星的时间长短也各有不同,为提高时频传递的精度,需尽可能的得到持续历元较长的观测值,首先选择观测到的某一颗卫星,然后选择持观测历元的观测值,保留相位观测值和伪距观测值;

第三步、采用mw组合法(melboutne-wubbena)构造初始检测量,gps信号是调制在频带为l1,l2上的,设第i历元的伪距观测值为p1,p2,相位观测值为宽巷模糊度为那么:

其中,λ1=19.06cm,λ2=24.45cm,为l1,l2的波长,f1=1575.42mhz,f2=1227.60mhz,为l1,l2的频率,为宽巷波长,约为86.2cm,c为光速。

在历元间对mw组合量做差,即可得到:

第四步、构造出的mw组合差分量是一个一维序列,利用matlab小波分解函数(wavedec函数)对此mw组合差分量进行小波分解,得到各层的细节分量,以四层小波分解为例,设得到细节分量为d1,d2,d3,d4,细节分量是与mw组合差分量长度相同的一维序列;

第五步、利用各层细节分量构造hankel矩阵,对所述矩阵进行奇异值分解,得到它们的奇异值矩阵,对奇异值进行差分谱和能量比分析,选择合适的奇异值进行降秩重构,滤除分量中的噪声,即可判断出发生周跳的历元;

第六步、利用未发生周跳的多个历元,进行自适应指数平滑方法预测出未来部分历元的值,将已修复的历元加入历史历元中,参与预测,修复下一个发生整周跳变的历元。

作为优选,第六步中,自适应指数平滑法通过对预测目标历史统计序列进行逐层的平滑计算,找出预测目标的基本变化趋势并以此预测,

mw组合差分量用式:计算,设选择持续观测历元数为t,组合差分量值用v1,v2,v3,...,vt表示,则一次指数平滑公式为:

式中:为第t周期的一次指数平滑值,α为加权系数,0<α<1,

预测模型为

为预测值,

二次指数平滑即为对一次指数平滑的结果做指数平滑。

采用三次指数平滑,即对二次指数平滑的结果做指数平滑,计算公式为:

三次指数平滑法的预测模型为:

其中:

使用误差平方和:表征预测精度,

采用0.618优选法确定最终的预测平滑系数α,0<α<1,首先选择α的值为0.618,使用三次指数平滑模型进行预测,利用预测出的结果计算误差平方和,第二次选择(1-0.618+0=0.382)为α的值,同样的方法进行预测,计算误差平方和,若第二次预测的误差平方和小于第一次预测,则去除0.618以上的部分,反之去除0.618一下的部分,对得到的新的区间继续进行0.618优选法选择α的值,直到选择出最优的α。

本发明的特征如下:

(1)能够准确的判断出多个历元内是否有周跳发生,和周跳发生的历元。

(2)能够实时的进行周跳探测。

(3)可探测出一周的小周跳。

与现有技术相比,本发明具有以下有益效果:

本发明提出了一种基于小波奇异值分解的gps载波相位周跳探测的方法,该方法与现行的电离层残差法、多项式拟合法和mw组合法等相比,可以准确探测出较小周跳的发生,漏检率大大降低,为gps载波相位时频传递精度提供了保障。采用自适应指数平滑法预测发生周跳历元的值以修复周跳,结果表明,预测值与实际值误差较小,精度较高。

附图说明

图1gps载波相位周跳探测流程图;

图2mw组合法探测周跳结果图;

图3奇异值能量比示意图;

图4奇异值差分谱图;

图5基于小波分解和奇异值分解探测周跳结果图;

图6自适应指数平滑法预测结果图。

具体实施方式

以下结合具体实施例,并参照附图,对本发明进一步详细说明。

如图1所示,本发明实施例提供一种利用小波奇异值分解和自适应指数平滑进行gps载波相位周跳探测与修复的方法,包括以下步骤:

步骤1、对rinex文件进行处理。

使用matlab读取rinex文件,得到一天内观测到的所有卫星和其观测值等信息。观测值一般包括相位观测值、c/a码伪距观测值、p码伪距观测值、多普勒频率等。

步骤2、从步骤1所得到的卫星中选择某一颗卫星,然后选择持续观测历元的观测值,保留相位观测值和p码伪距观测值。

步骤3、构造mw组合差分量。

gps信号是调制在频带为l1,l2上的,设第i历元的伪距观测值为p1,p2,相位观测值为第i个历元的宽巷模糊度可用下式表示:

其中,λ1=19.06cm,λ2=24.45cm,为l1,l2的波长,f1=1575.42mhz,f2=1227.60mhz,为l1,l2的频率,为宽巷波长,约为86.2cm,c为光速。

在历元间对其做差,即可得到:

以此作为后续周跳探测的检测量。此检测量为一维序列,本发明以长度为400的序列为例。

步骤4、利用matlab函数(wavedec)对mw组合差分量进行四层小波分解。得到四层细节分量,分别为d1,d2,d3,d4,四层细节分量长度都为399。设d1=(x1,x2,...,x399)

步骤5、利用细节分量构造hankel矩阵。

对于一个一维的信号序列,为对其进行奇异值分解处理,必须首先构造一个矩阵。本文采用hankel矩阵的生成方式。hankel矩阵的构造形式如下:

hankel矩阵的特点为,其反对角线上的元素相同。若利用时间序列构造hankel矩阵,则下一行矢量元素仅比上一行元素滞后一个时间点。

对于mw组合差分量进行小波分解重构得到细节分量,其长度为n,利用此序列构造hankel矩阵,若n为偶数,则令此矩阵的行数m=n/2+1,列数n=n/2+1;若n为奇数,则令矩阵行数m=(n+1)/2,列数n=(n+1)/2。

本发明以长度400的序列为例,细节分量长度为399。构造的hankel矩阵行数为200,列数为200。d1,d2,d3,d4构造的矩阵为m1,m2,m3,m4。以第一层分量d1=(x1,x2,...,x399)为例,构造出的m1为:

步骤6、对构造成的hankel矩阵进行奇异值分解。

奇异值分解(singularvaluedecomposition,svd)是一种矩阵的正交化分解的方法。设m是一个m×n阶的实矩阵,则必定存在一个正交矩阵u∈rm×m和另一个正交矩阵v∈rn×n,使得

m=udvt(4)

式中,d是一个半正定对角矩阵,且d∈rm×n,称为奇异值矩阵,矩阵d可表示为:

式中,矩阵s=diag(σ1,σ2,σ3,…,σp),其中p=min(m,n),且σ1≥σ2≥σ3≥…≥σp,对角线上的元素为矩阵m的奇异值,奇异值的数目及为此矩阵的秩。

本发明中进行奇异值分解的矩阵为200×200的矩阵,得到其奇异值矩阵d大小为200×200,即得到200个奇异值。

步骤7、奇异值能量比和差分谱分析。

将奇异值由大到小排列成一个序列,每个奇异值的能量比和差分定义为:

δσi=σi+1-σi(6)

步骤8、降秩重构

结合差分谱和能量比,保留快速下降及之前的奇异值进行降秩重构,即可去除噪声。

以第一层分量为例,如图3、4所示,本例保留前60个奇异值,原奇异值矩阵d对角线上有200个不为零的数据,即矩阵m有200个奇异值。现保留前60个较大的奇异值重新构造奇异值矩阵d,奇异值依旧排列在对角线上,对角线的后140个值为零。再利用式:m=udvt重构出新的矩阵m。再根据hankel矩阵m的特点,还原出长度为399的细节分量。

重构后的细节分量如图5所示,100、200、300历元有峰值,即可判断出周跳发生在此三个历元上。

步骤9、指数平滑法修复周跳。

本发明中,以第一个周跳为例,100历元发生周跳,则利用前99个数据进行预测,预测出的100历元值取代实际值,即可修复周跳。本发明采用自适应指数平滑法进行周跳修复。

指数平滑法通过对预测目标历史统计序列进行逐层的平滑计算,消除由于随机因素造成的影响,找出预测目标的基本变化趋势并以此预测。

mw组合差分量用式:计算,设选择持续观测历元数为t,组合差分量值用v1,v2,v3,...,vt表示,则一次指数平滑公式为:

式中:为第t周期的一次指数平滑值,α为加权系数,0<α<1。

预测模型为

为预测值。

二次指数平滑即为对一次指数平滑的结果做指数平滑。

本发明采用三次指数平滑,即对二次指数平滑的结果做指数平滑。计算公式为:

三次指数平滑法的预测模型为:

其中:

本发明使用误差平方和:表征预测精度。

步骤10、确定最终的平滑系数。

本发明采用0.618优选法确定最终的预测平滑系数α。0<α<1,首先选择α的值为0.618,使用三次指数平滑模型进行预测,利用预测出的结果计算误差平方和。第二次选择(1-0.618+0=0.382)为α的值,同样的方法进行预测,计算误差平方和,若第二次预测的误差平方和小于第一次预测,则去除0.618以上的部分,反之去除0.618一下的部分,对得到的新的区间继续进行0.618优选法选择α的值,直到选择出最优的α。

利用最优的α的值进行三次指数平滑预测,预测到的第100个值取代原始值,即可修复周跳,如图6所示。

本发明的特点在于:

1、利用小波分解,得到mw组合差分量的各层细节分量更好的反映出周跳发生的现象。

2、利用svd分解对细节分量进行降秩重构,一定程度上减少了噪声对于周跳探测的影响,大大降低了漏检率。

3、利用自适应指数平滑法预测,以修复周跳。

以上实施例仅为本发明的示例性实施例,不用于限制本发明,本发明的保护范围由权利要求书限定。本领域技术人员可以在本发明的实质和保护范围内,对本发明做出各种修改或等同替换,这种修改或等同替换也应视为落在本发明的保护范围内。

当前第1页1 2 
网友询问留言 已有0条留言
  • 还没有人留言评论。精彩留言会获得点赞!
1