一种用于削弱GNSS多路径误差的组合滤波方法

一种用于削弱gnss多路径误差的组合滤波方法
技术领域
1.本发明属于测绘领域，涉及一种用于削弱gnss多路径误差的组合滤波方法.


背景技术：

2.gnss变形监测技术具有采样高频率、实时高精度、无需通视、全自动、全天候等优点，可应用于土木工程结构与自然灾害的实时监测。gnss变形监测通常以短基线静态测量的方式进行，采用载波相位差分技术能基本消除卫星与接收机钟差、星历误差、大气延迟等相关性较强的误差，但变形监测环境往往比较复杂，由测站附近反射物造成的多路径效应不仅无法避免，而且难以通过差分削弱，直接导致gnss定位精度下降。另外，由测站接收机自生产生的噪声频域分布范围宽、高斯白噪声特性明显，其与多路径误差混叠，进一步增加多路径误差处理难度。
3.gnss测站周边环境变化较小或是未发生改变时，多路径效应呈现出较强的周日重复性，基于此特性利用恒星日滤波对连续多天gnss监测数据进行多路径误差提取、重复性建模能有效削弱多路径误差，处理效率快、性价比高。其中，多路径误差提取在整个处理过程十分关键，现主要的提取方法有小波分析(wavelet)系列、经验模态分解(emd)系列等。小波分析(wavelet)凭其多分辨率分析、时频域同时局部化的特点，在多路径误差处理上展现出明显优势，但其效果受限于小波基、分解层数、阈值及其估计方法等参数确定，参数选择没有统一的支撑理论；集合经验模态分解(eemd)是emd系列中的改进方法，它通过添加高斯白噪声、多次emd、均值处理有效改善了模态混叠问题，并且继承了能把信号自适应分解为若干个本征模态函数(imf)项与1个余项的特点，能够有效、快速地提取到信号中的有用成分，在多路径误差处理上展现出独特的优势，但其效果受限于白噪声添加量，白噪声量过多过少都会导致最终效果不理想。
4.作为gnss变形监测的主要误差源，多路径效应严重影响变形监测准确度。随着gnss技术在变形监测中的应用普及，迫切需要一种便于计算机运算、理论严密、性价比高的多路径误差处理方法，实现gnss实时高精度监测。


技术实现要素：

5.为了解决上述背景技术提到的技术问题，本发明一种用于削弱gnss多路径误差的组合滤波方法，针对wavelet、eemd的多路径误差处理上存在的不足实现更有效地削弱多路径误差。
6.为了实现上述技术目的，本发明的技术方案为：
7.一种用于削弱gnss多路径误差的组合滤波方法，包括以下步骤：
8.步骤一、读取gnss观测数据的原始坐标序列，读取完后开始解算；
9.步骤二、利用集合经验模态分解eemd将原始坐标序列分解为一系列本征模态函数项imf和余项；
10.步骤三、利用分类指标k1和k2将imf项与余项细分为噪声项、过渡项和有用项；
11.步骤四、舍弃噪声项，对过渡项进行小波滤波wavelet处理，得到滤波后的过渡项；
12.步骤五、对有用项和滤波后的过渡项进行重构，得到多路径误差序列；
13.步骤六、将第一天的多路径误差序列作为参考信号序列，计算参考信号序列与其余天多路径误差序列间的幅值比例系数；
14.步骤七、恢复其余天的多路径误差序列幅值，并以此作为误差改正模型，与对应的原始坐标序列相减，得到改正后的坐标序列；
15.步骤八、结束解算，输出多路径误差改正后的坐标序列。
16.优选地，所述步骤二中，eemd处理过程步骤如下：
17.21)对原始坐标序列y(t)加入高斯白噪声n(t)得到待处理坐标序列ym(t)，加入的高斯白噪声n(t)大小由其与原始坐标序列y(t)间的标准差比值σ决定；
18.22)对ym(t)进行emd处理，得到一组imf项与余项，emd处理步骤如下：
19.(221)令ym(t)＝r
i-1
(t)，i＝1；
20.(222)寻找r
i-1
(t)上所有的极值点，利用三次样条插值分别拟合极小值点、极大值点序列，获得r
i-1
(t)上、下包络线；
21.(223)计算上、下包络线的均值序列mi(t)，并与r
i-1
(t)做差hi(t)＝r
i-1
(t)-mi(t)；
22.(224)判断hi(t)是否满足：
①hi
(t)极值点与零点数相等或仅差1；
②hi
(t)上、下包络线关于时间轴局部对称；若满足，hi(t)＝imfi(t)，否则，令r
i-1
(t)＝hi(t)，跳至步骤(222)直至满足条件；
23.(225)计算ri(t)＝r
i-1
(t)-imfi(t)，判断ri(t)是否为单调或其极值点个数是否小于2，若ri(t)为单调或其极值点个数小于2，则分解结束，令n＝i，否则，令i＝i 1，跳至步骤(222)，直至满足条件，得到一组由n个imf项与1个余项组成的分解结果，表达式如下：
[0024][0025]
式中，为由emd得到的第i个imf项，为余项；i＝1,2,
…
,n，n为imf项的总个数；
[0026]
23)重复步骤21)与22)n次，得到n组imf项与余项；
[0027]
24)对n组imf项与余项分别求均值，得到的n个imf项均值与1个余项均值作为最终结果，即由eemd得到的n个imf项均值与1个余项均值，如下式所示：
[0028][0029]
式中，j＝1,2,
…
,n，n为emd重复总次数；imfi(t)为第i个imf项均值，rn(t)为余项均值。
[0030]
优选地，所述步骤三中，分类指标k1的确定步骤如下：
[0031]
31)通过eemd处理，从第k个imf项均值开始重构的坐标序列表示为如下：
[0032][0033]
32)利用连续均方根误差准则计算2个连续重构的坐标序列间的平方欧氏距离，计算如下：
[0034][0035]
式中，k＝1,2,3,....,n-1；l为坐标序列的长度，即历元数；
[0036]
33)根据k取值范围计算出cmse值，当cmse值首次出现极小值处的k值作为过渡项的起点，即k＝k1：
[0037][0038]
式中，argfirstlocalmin(
·
)为获取首次极小值时对应的k值；为不超过的最大整数，n为imf项均值的总个数。
[0039]
优选地，所述步骤三中，分类指标k2的确定步骤如下：
[0040]
41)根据eemd分解结果，计算第k个imf均值项的能量密度与其平均周期的乘积etk，公式如下：
[0041][0042]
式中，ek为第k个imf均值项能量密度，为第k个imf均值项的平均周期，计算公式分别如下：
[0043][0044][0045]
式中，k＝1,2,3,....,n，ok为第k个imf项的极值点数；l为第k个imf项的长度，即历元数；
[0046]
42)计算第k个imf项的有效系数rk，表达式如下：
[0047][0048]
式中，k＝2,3,....,n；
[0049]
43)根据k取值范围依此计算rk，将首次出现rk》3时的k值作为有用项的起点，即k＝k2。
[0050]
优选地，所述步骤四中，小波滤波的处理步骤如下：
[0051]
51)对过渡项序列进行小波分解，得到一组由多个高频、低频项组成的分解结果；
[0052]
52)对高频项序列采用软阈值化处理，低频项不作处理，软阈值函数soft(
·
)表达如下：
[0053]
[0054]
式中，x
t
为高频项序列幅值，t＝1,2,3,....,l，l为序列长度，即历元数；sgn(
·
)为阶跃函数，x
t
》0结果取1，x
t
《0结果取-1，x
t
＝0取0；λ为阈值，计算如下：
[0055]
53)重构低频项与软阈值化后的高频项，得到滤波后的过渡项。
[0056]
所述步骤五中，重构得到的gnss多路径误差序列表达如下:
[0057][0058]
式中，wavelet为小波滤波处理；imfi(t)为由eemd得到的第i个imf项，rn(t)为余项；k1与k2分别为区分噪声项与过渡项、过渡项与有用项的指标。优选地，所述步骤六中，幅值比例系数ar计算如下:
[0059][0060]
式中，average(
·
)为均值处理；r(t)为参考信号序列，为重构得到的多路径误差序列。
[0061]
优选地，所述步骤七中，多路径误差幅值恢复计算如下：
[0062][0063]
式中，s(t)为幅值恢复后的多路径误差；ar为幅值比例系数；为重构得到的多路径误差序列。
[0064]
采用上述技术方案带来的有益效果：
[0065]
本发明充分考虑了gnss变形监测中多路径效应的物理特性以及wavelet与eemd方法在多路径误差处理上的不足，在eemd-wavelet组合滤波的基础上引入一种组合分类准则。本方法强化了eemd与wavelet在多路径误差处理上的可靠性，不仅能更有效地提取、改正多路径误差，而且流程简单、算法统一，能够适应实际工程。
附图说明
[0066]
图1是一种用于削弱gnss多路径误差的组合滤波方法得实施流程图；
[0067]
图2是3天的高程方向坐标序列图；
[0068]
图3是eemd分解结果图；
[0069]
图4是滤波后的过渡项波形图；
[0070]
图5是第一天多路径误差序列图；
[0071]
图6是多路径误差改正前后的坐标序列图。
具体实施方式
[0072]
以下将结合附图，对本发明的技术方案进行详细说明。
[0073]
本发明提出了一种用于削弱gnss多路径误差的组合滤波方法，实施流程如图1所示。
[0074]
本实施案例数据来源于科廷大学超短基线gnss观测数据，基线长度为6.15m，采用的接收机为trimble-netr9，观测历元间隔为30s，截至高度角为10
°
，利用载波相位差分解算得到的原始坐标序列中主要为多路径误差与随机噪声。
[0075]
步骤一、计算机读取连续三天、同一时段、高程方向、长度为2500历元的原始坐标序列，如图2所示，取第一天的原始坐标序列对接下来的步骤进行阐述。
[0076]
步骤二，利用eemd对原始坐标序列进行分解，参数默认，得到10个imf项与1个余项，图3所示为eemd分解结果；
[0077]
步骤三，根据eemd分解结果，利用cmse准则与基于能量密度与平均周期定义的有效系数r依次计算分类指标k1与k2，计算结果如表1所示。
[0078]
表1.cmse值与能量有效系数r
[0079]
imf项数k12345678910cmse值11.383.752.993.391.852.291.370.451.280.03r值-0.320.571.730.280.903.460.9512.525.39
[0080]
表中，第3阶imf项的cmse值首次出现极小值，k1＝3，第6阶imf项的r值大于3，k2＝6，imf1～2项为噪声项、imf3～6项为过渡项，imf7～10与余项为有用项。
[0081]
步骤四，剔除噪声项，对过渡项进行小波滤波处理，参数默认，得到滤波后的过渡项；如图4所示。
[0082]
步骤五，对滤波后的过渡项与有用项重构得到多路径误差；
[0083]
步骤六，第一天的多路径误差作为参考信号，再利用提取到的往后两天多路径误差分别计算幅值比例系数ar，ar得值如图5和表2所示；
[0084]
表2.ar值
[0085]
编号day1-day2day1-day3ar0.960.97
[0086]
步骤七，恢复多路径误差幅值，以此作为误差改正模型，与原始坐标序列做减实现多路径误差改正，图6与表3分别为改正后的坐标序列及其rmse：
[0087]
表3.多路径误差改正后的坐标序列rmse
[0088][0089]
实施例仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明保护范围之内。