CN107153208B - 一种gps载波相位周跳探测与修复的方法 - Google Patents
一种gps载波相位周跳探测与修复的方法 Download PDFInfo
- Publication number
- CN107153208B CN107153208B CN201710469425.0A CN201710469425A CN107153208B CN 107153208 B CN107153208 B CN 107153208B CN 201710469425 A CN201710469425 A CN 201710469425A CN 107153208 B CN107153208 B CN 107153208B
- Authority
- CN
- China
- Prior art keywords
- value
- exponential smoothing
- prediction
- cycle slip
- observation
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
- 238000000034 method Methods 0.000 title claims abstract description 46
- 238000000819 phase cycle Methods 0.000 title claims abstract description 10
- 238000001514 detection method Methods 0.000 title claims description 15
- 238000009499 grossing Methods 0.000 claims abstract description 53
- 238000000354 decomposition reaction Methods 0.000 claims abstract description 22
- 230000003044 adaptive effect Effects 0.000 claims abstract description 9
- 239000011159 matrix material Substances 0.000 claims description 39
- 238000004364 calculation method Methods 0.000 claims description 6
- 238000005457 optimization Methods 0.000 claims description 6
- 238000001228 spectrum Methods 0.000 claims description 4
- 238000001914 filtration Methods 0.000 claims description 2
- 230000005540 biological transmission Effects 0.000 description 9
- 238000011160 research Methods 0.000 description 6
- 238000010586 diagram Methods 0.000 description 3
- 108010076504 Protein Sorting Signals Proteins 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 230000000717 retained effect Effects 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000000903 blocking effect Effects 0.000 description 1
- 239000005433 ionosphere Substances 0.000 description 1
- 238000010183 spectrum analysis Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
- G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
- G01S19/13—Receivers
- G01S19/35—Constructional details or hardware or software details of the signal processing chain
- G01S19/37—Hardware or software details of the signal processing chain
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
- G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
- G01S19/13—Receivers
- G01S19/23—Testing, monitoring, correcting or calibrating of receiver elements
- G01S19/235—Calibration of receiver components
Landscapes
- Engineering & Computer Science (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Computer Networks & Wireless Communication (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Signal Processing (AREA)
- Position Fixing By Use Of Radio Waves (AREA)
Abstract
本发明公开一种利用小波奇异值分解和自适应指数平滑进行GPS载波相位周跳探测与修复的方法,利用GPS观测文件,采用小波奇异值分解方法,探测出GPS载波相位传递中,整周跳变的现象以及跳变发生的历元,并采用自适应指数平滑法以修复周跳,保证GPS载波相位时频传递的精度。
Description
技术领域
本发明属于时间频率传递的技术领域,尤其涉及一种利用小波奇异值分解和自适应指数平滑进行GPS载波相位周跳探测与修复的方法。
背景技术
时间频率的研究是基础科学研究中的一个重要分支,时间频率在科研、定位系统、电力系统、军事和国家安全等方面处于举足轻重的地位。世界各国大都建有自己的守时实验室,产生本国的原子时标,并参与国际比对。
中国计量科学研究院承担着向地方实验室发布标准时间与频率,并参与国际原子时比对的任务。因此,与国际原子时比对,以获得稳定度及准确度更高的本地原子时是至关重要的。中国计量科学研究院目前采用的时间比对方法为GPS共视法,但共视法要求两地的实验室在同一时刻需观测到同一颗卫星,且持续较长时间。若两地实验室距离较远,无法观测到同一颗卫星,则需第三个实验室作为中继。GPS载波相位法无需考虑观测卫星的问题,且时频传递精度更高。但GPS载波相位时频传递中有一个重要的问题,即周跳。周跳是指由于GPS接收机失锁,或接收载波由于其他原因被阻断,而产生间断的现象。周跳的发生对GPS载波相位时频传递影响较大。故而周跳探测与修复是非常重要的。本课题研究GPS载波相位周跳探测与修复方法,旨在保证时频传递的精度,并提高本地原子时的稳定度和准确度。
发明内容
本发明要解决的技术问题是,提供一种利用小波奇异值分解和自适应指数平滑进行GPS载波相位周跳探测与修复的方法,利用GPS观测文件,采用小波奇异值分解方法,探测出GPS载波相位传递中,整周跳变的现象以及跳变发生的历元,并采用自适应指数平滑法以修复周跳,保证GPS载波相位时频传递的精度。
为解决上述问题,本发明采用如下的技术方案:
一种GPS载波相位周跳探测与修复的方法,包括以下步骤:
第一步、利用Matlab读取Rinex型的GPS观测文件,得到一天内观测到的卫星号、GPS伪距观测值和相位观测值;
第二步、选择伪距观测值以及相位观测值,由于同一个历元,GPS接收机所接收到的观测值来自多颗卫星,且跟踪同一颗卫星的时间长短也各有不同,为提高时频传递的精度,需尽可能的得到持续历元较长的观测值,首先选择观测到的某一颗卫星,然后选择持观测历元的观测值,保留相位观测值和伪距观测值;
在历元间对MW组合量做差,即可得到:
第四步、构造出的MW组合差分量是一个一维序列,利用Matlab小波分解函数(wavedec函数)对此MW组合差分量进行小波分解,得到各层的细节分量,以四层小波分解为例,设得到细节分量为D1,D2,D3,D4,细节分量是与MW组合差分量长度相同的一维序列;
第五步、利用各层细节分量构造Hankel矩阵,对所述矩阵进行奇异值分解,得到它们的奇异值矩阵,对奇异值进行差分谱和能量比分析,选择合适的奇异值进行降秩重构,滤除分量中的噪声,即可判断出发生周跳的历元;
第六步、利用未发生周跳的多个历元,进行自适应指数平滑方法预测出未来部分历元的值,将已修复的历元加入历史历元中,参与预测,修复下一个发生整周跳变的历元。
作为优选,第六步中,自适应指数平滑法通过对预测目标历史统计序列进行逐层的平滑计算,找出预测目标的基本变化趋势并以此预测,
预测模型为
二次指数平滑即为对一次指数平滑的结果做指数平滑。
采用三次指数平滑,即对二次指数平滑的结果做指数平滑,计算公式为:
三次指数平滑法的预测模型为:
其中:
采用0.618优选法确定最终的预测平滑系数α,0<α<1,首先选择α的值为0.618,使用三次指数平滑模型进行预测,利用预测出的结果计算误差平方和,第二次选择(1-0.618+0=0.382)为α的值,同样的方法进行预测,计算误差平方和,若第二次预测的误差平方和小于第一次预测,则去除0.618以上的部分,反之去除0.618一下的部分,对得到的新的区间继续进行0.618优选法选择α的值,直到选择出最优的α。
本发明的特征如下:
(1)能够准确的判断出多个历元内是否有周跳发生,和周跳发生的历元。
(2)能够实时的进行周跳探测。
(3)可探测出一周的小周跳。
与现有技术相比,本发明具有以下有益效果:
本发明提出了一种基于小波奇异值分解的GPS载波相位周跳探测的方法,该方法与现行的电离层残差法、多项式拟合法和MW组合法等相比,可以准确探测出较小周跳的发生,漏检率大大降低,为GPS载波相位时频传递精度提供了保障。采用自适应指数平滑法预测发生周跳历元的值以修复周跳,结果表明,预测值与实际值误差较小,精度较高。
附图说明
图1 GPS载波相位周跳探测流程图;
图2 MW组合法探测周跳结果图;
图3奇异值能量比示意图;
图4奇异值差分谱图;
图5基于小波分解和奇异值分解探测周跳结果图;
图6自适应指数平滑法预测结果图。
具体实施方式
以下结合具体实施例,并参照附图,对本发明进一步详细说明。
如图1所示,本发明实施例提供一种利用小波奇异值分解和自适应指数平滑进行GPS载波相位周跳探测与修复的方法,包括以下步骤:
步骤1、对Rinex文件进行处理。
使用Matlab读取Rinex文件,得到一天内观测到的所有卫星和其观测值等信息。观测值一般包括相位观测值、C/A码伪距观测值、P码伪距观测值、多普勒频率等。
步骤2、从步骤1所得到的卫星中选择某一颗卫星,然后选择持续观测历元的观测值,保留相位观测值和P码伪距观测值。
步骤3、构造MW组合差分量。
在历元间对其做差,即可得到:
以此作为后续周跳探测的检测量。此检测量为一维序列,本发明以长度为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矩阵进行奇异值分解。
奇异值分解(Singular Value Decomposition,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历元值取代实际值,即可修复周跳。本发明采用自适应指数平滑法进行周跳修复。
指数平滑法通过对预测目标历史统计序列进行逐层的平滑计算,消除由于随机因素造成的影响,找出预测目标的基本变化趋势并以此预测。
预测模型为
二次指数平滑即为对一次指数平滑的结果做指数平滑。
本发明采用三次指数平滑,即对二次指数平滑的结果做指数平滑。计算公式为:
三次指数平滑法的预测模型为:
其中:
步骤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、利用自适应指数平滑法预测,以修复周跳。
以上实施例仅为本发明的示例性实施例,不用于限制本发明,本发明的保护范围由权利要求书限定。本领域技术人员可以在本发明的实质和保护范围内,对本发明做出各种修改或等同替换,这种修改或等同替换也应视为落在本发明的保护范围内。
Claims (2)
1.一种GPS载波相位周跳探测与修复的方法,其特征在于,包括以下步骤:
第一步、利用Matlab读取Rinex型的GPS观测文件,得到一天内观测到的卫星号、GPS伪距观测值和相位观测值;
第二步、选择观测到的某一颗卫星,然后选择观测历元的观测值,保留相位观测值和伪距观测值;
在历元间对MW组合量做差,即可得到:
第四步、构造出的MW组合差分量是一个一维序列,利用Matlab小波分解函数对此MW组合差分量进行四层小波分解,得到各层的细节分量,所述细节分量是与MW组合差分量长度相同的一维序列;
第五步、利用各层细节分量构造Hankel矩阵,对所述矩阵进行奇异值分解,得到它们的奇异值矩阵,对奇异值进行差分谱和能量比分析,选择合适的奇异值进行降秩重构,滤除分量中的噪声,即可判断出发生周跳的历元;
第六步、利用未发生周跳的多个历元,进行自适应指数平滑方法预测出未来部分历元的值,将已修复的历元加入历史历元中,参与预测,修复下一个发生整周跳变的历元。
2.如权利要求1所述的GPS载波相位周跳探测与修复的方法,其特征在于,第六步中,自适应指数平滑法通过对预测目标历史统计序列进行逐层的平滑计算,找出预测目标的基本变化趋势并以此预测,
式中:Vt (1)为第t周期的一次指数平滑值,α为平滑系数,0<α<1,
预测模型为
二次指数平滑即为对一次指数平滑的结果做指数平滑,
采用三次指数平滑,即对二次指数平滑的结果做指数平滑,计算公式为:
三次指数平滑法的预测模型为:
其中:
αt=3Vt (1)-3Vt (2)+Vt (3)
采用0.618优选法确定最终的预测平滑系数α,0<α<1,首先选择α的值为0.618,使用三次指数平滑模型进行预测,利用预测出的结果计算误差平方和,第二次选择(1-0.618+0=0.382)为α的值,同样的方法进行预测,计算误差平方和,若第二次预测的误差平方和小于第一次预测,则去除0.618以上的部分,反之去除0.618以下的部分,对得到的新的区间继续进行0.618优选法选择α的值,直到选择出最优的α。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201710469425.0A CN107153208B (zh) | 2017-06-20 | 2017-06-20 | 一种gps载波相位周跳探测与修复的方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201710469425.0A CN107153208B (zh) | 2017-06-20 | 2017-06-20 | 一种gps载波相位周跳探测与修复的方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN107153208A CN107153208A (zh) | 2017-09-12 |
CN107153208B true CN107153208B (zh) | 2020-06-19 |
Family
ID=59796380
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201710469425.0A Expired - Fee Related CN107153208B (zh) | 2017-06-20 | 2017-06-20 | 一种gps载波相位周跳探测与修复的方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN107153208B (zh) |
Families Citing this family (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN107728168B (zh) * | 2017-11-09 | 2021-08-20 | 昆明理工大学 | 一种基于形态滤波和奇异值分解的周跳检测方法 |
CN109358350B (zh) * | 2018-10-08 | 2021-02-05 | 中国人民解放军战略支援部队信息工程大学 | 一种北斗三频周跳探测方法与装置 |
CN110505009A (zh) * | 2019-09-12 | 2019-11-26 | 国家电网有限公司 | 一种基于相干光时域反射计的电力光缆监控装置及方法 |
CN110826017A (zh) * | 2019-09-25 | 2020-02-21 | 中国地质大学(武汉) | 一种基于参数优化Hankel矩阵和奇异值分解的信号去噪方法 |
CN110967717A (zh) * | 2019-12-23 | 2020-04-07 | 合肥工业大学 | 一种基于小波变换法的周跳探测和修复方法 |
CN111239786B (zh) * | 2020-02-28 | 2021-04-13 | 同济大学 | 一种无人驾驶定位测姿的载波相位整周模糊度测定方法 |
CN111680398B (zh) * | 2020-05-06 | 2023-06-27 | 北京航空航天大学 | 一种基于Holt-Winters模型的单机性能退化预测方法 |
Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN104570013A (zh) * | 2014-12-30 | 2015-04-29 | 北京无线电计量测试研究所 | 一种用于频率驯服的实时gps载波相位周跳的探测方法 |
CN105137459A (zh) * | 2015-07-29 | 2015-12-09 | 昆明理工大学 | 一种北斗单频周跳探测方法 |
CN106125107A (zh) * | 2016-06-14 | 2016-11-16 | 昆明理工大学 | 一种利用mw与小波变换探测北斗周跳的方法 |
-
2017
- 2017-06-20 CN CN201710469425.0A patent/CN107153208B/zh not_active Expired - Fee Related
Patent Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN104570013A (zh) * | 2014-12-30 | 2015-04-29 | 北京无线电计量测试研究所 | 一种用于频率驯服的实时gps载波相位周跳的探测方法 |
CN105137459A (zh) * | 2015-07-29 | 2015-12-09 | 昆明理工大学 | 一种北斗单频周跳探测方法 |
CN106125107A (zh) * | 2016-06-14 | 2016-11-16 | 昆明理工大学 | 一种利用mw与小波变换探测北斗周跳的方法 |
Non-Patent Citations (4)
Title |
---|
GPS周跳探测及修复的小波变换法;蔡昌盛 等;《武汉大学学报》;20070131;第39-42页 * |
Small cycle-slip detection of single-frequency in BDS based on SVD;Yang Gao et al.;《The 27th Chinese Control and Decision Conference (2015 CCDC)》;20150720;第3627-3631页 * |
一种非差相位观测值的周跳探测与修复方法;夏朋飞 等;《测绘》;20111231;第243-245页 * |
基于小波变换的GPS精密单点定位中的周跳探测;黄兵杰 等;《武汉大学学报》;20060630;第512-515页 * |
Also Published As
Publication number | Publication date |
---|---|
CN107153208A (zh) | 2017-09-12 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN107153208B (zh) | 一种gps载波相位周跳探测与修复的方法 | |
CN102597701B (zh) | 用于补偿错误测量的系统和方法 | |
Li | Cycle slip detection and ambiguity resolution algorithms for dual-frequency GPS data processing | |
KR100964935B1 (ko) | 이력 상관 데이터를 이용하여 신호 상관을 수행하는 방법및 장치 | |
CN104714244A (zh) | 一种基于抗差自适应Kalman滤波的多系统动态PPP解算方法 | |
CN109752739B (zh) | 一种观测数据处理方法、装置、终端设备及存储介质 | |
CN111522032B (zh) | 一种北斗三号系统用户完好性处理的优化方法及优化装置 | |
CN116774264B (zh) | 基于低轨卫星机会信号多普勒的运动目标定位方法 | |
CN111123322B (zh) | 星载gnss接收机的观测值实时数据预处理方法、系统、介质及设备 | |
CN117724125B (zh) | 一种基于一致性的观测数据的质量控制方法及装置 | |
CN114280640A (zh) | 高精度定位方法、计算装置和计算机可读介质 | |
CN117388896A (zh) | 利用载波相位残差模型的gnss测量处理 | |
CN113671551A (zh) | Rtk定位解算方法 | |
CN111999750B (zh) | 针对杆臂不准的实时单站周跳探测改进方法 | |
Kaddour et al. | Fault detection and exclusion for GNSS measurements using observations projection on information space | |
CN105068092B (zh) | 一种应用于星基增强系统机载接收机的周跳检测与修复方法 | |
Kim et al. | Predicting IGS RTS corrections using ARMA neural networks | |
Di Lecce et al. | Neural technologies for increasing the GPS position accuracy | |
Tabatabaei et al. | MP Mitigation in Urban Canyons using GPS‐combined‐GLONASS Weighted Vectorized Receiver | |
CN112415547B (zh) | 卫星信号的周跳计算方法及装置 | |
CN117724128B (zh) | 一种低轨卫星轨道预报方法、系统、终端及介质 | |
CN116794700B (zh) | 一种用于船载北斗一体机的卫星故障检测方法 | |
Jwo et al. | Navigation Performance Enhancement Using IMM Filtering for Time Varying Satellite Signal Quality | |
De Weerdt et al. | New approach for integer ambiguity resolution using interval analysis | |
CN116893435A (zh) | 一种基于注意力机制增强时间卷积网络的三星定位方法 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant | ||
CF01 | Termination of patent right due to non-payment of annual fee | ||
CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20200619 |