CN110531386A - 基于wmad的gnss卫星钟差数据预处理方法 - Google Patents

Abstract

本发明的目的是提供一种基于WMAD的GNSS卫星钟差数据预处理方法。本发明有效解决了经典粗差探测方法不能有效剔除钟差中量级较小的误差问题，本发明的方法具体是首先对GNSS卫星钟差进行一次差计算，得到钟差一次差分数据；然后对钟差一次差分数据进行小波分解，得到分解后小波系数；再利用MAD方法剔除较大的小波系数；计算各层阈值，进行阈值化处理；最后小波重构得到预处理后的钟差一次差分数据。本发明为卫星钟差数据预处理提供一种有效方法，提高钟差数据的质量。

Description

基于WMAD的GNSS卫星钟差数据预处理方法

技术领域

本发明涉及一种GNSS卫星钟差数据处理方法，属于大地测量与导航测绘技术领域。

背景技术

GNSS星载原子钟在长期的运行中会出现调频、钟切换等操作，同时还会受到各种不确定因素的影响，钟差数据不可避免地含有一些粗差，如果用含有粗差的钟差数据进行数据分析及建模预报，会大大影响数据分析及建模预报的质量。因此，国内外学者在卫星钟差粗差剔除方面进行了大量的研究，对于钟差数据中的粗差处理，常用的方法有中位数(Median Absolute Deviation,MAD)方法、抗差估计方法以及Bayesian方法等。其中，MAD方法具有原理简单、计算效率高等优点，是目前钟差数据预处理中普遍采用的一种方法。

通过上述方法可以有效剔除钟差数据中含有的粗差，但在实际的钟差处理过程中，发现剔除完粗差的钟差数据中常会存在量级较小的误差(异常值)，而经典的粗差探测的方法对于这些异常值很难进行有效的剔除，而这些异常值也会对影响钟差数据分析及建模预报的质量，所以在使用钟差数据之前，对钟差数据中量级较小的误差进行合理有效的预处理是非常必要的步骤。

发明内容

本发明的目的是提出一种基于WMAD的GNSS卫星钟差数据预处理方法，本发明分析了小波分析和中位数方法对应含义，并将两种方法融合引用到GNSS卫星钟差数据预处理过程中，提高GNSS卫星钟差数据预处理的质量。

为了实现上述目的，本发明采用下述技术方案：

基于WMAD的GNSS卫星钟差数据预处理方法，包括以下步骤：

A：对GNSS卫星钟差进行一次差计算，得到钟差一次差分数据；

B：对GNSS卫星钟差一次差分数据进行小波分解，得到分解后小波系数；

C：利用MAD方法剔除较大的小波系数；

D：计算各层阈值，进行阈值化处理；

E：小波重构得到预处理后的钟差一次差分数据。

进一步地，所述的步骤A包括以下步骤：

a：首先设L＝{l(i),i＝1,2,…,n}为一组不同历元时刻的钟差值；

b：计算卫星钟差一次差分数据，即

ΔL(j)＝l(i+1)-l(i),j＝1,2,…n-1

其中，ΔL(j)为钟差一次差分数据。

进一步地，所述的步骤B包括以下步骤：

c：确定小波函数；

分别利用不同小波函数对钟差一次差分数据进行小波分解，得到分解后的小波系数，即

其中，WT_f(a,b)为小波变换系数，f(g)为卫星钟差一次差分数据，ψ_a,b(g)为小波母函数，a为伸缩因子，b平移因子(g＝0，1,2，…T-1,T为卫星钟差一次差数据个数)；

计算小波变换后小波系数的能量值，即

其中，E_s(W)为提取信号的能量值，N为小波系数的个数，W(s,t)为小波系数，s表示尺度参数，t表示信号幅值；获得最大能量值的小波函数为选择的最优小波函数；

d：确定小波分解层数；

根据尺度与频率的对应关系，即：

其中f_s ^ps是尺度s对应的伪频率ps(pseudo frequency),f_c是小波对应的中心频率，Δt是采用周期；伪频率应覆盖有用信号f_sig的全部频率范围，即：

因此有用信号的最小频率(minf_sig)应该不小于尺度s为L的伪频率，也就是：

离散小波变换的分解层数是对尺度s进行二进离散采样，因此L对应的分解层数j为：

2^j-1≤L≤2^j

e：对钟差数据进行小波分解得出小波系数。

进一步地，所述的步骤C包括以下步骤：

f：然后对每层的小波系数WT_f(a,b)与小波系数序列的中数(MED)m加上中位数(MAD)的n倍之和相比较，即若小波系数满足

|WT_f(a,b)|＞(m+n×MAD)

就认为该小波系数对应的就是较大异常值，同时剔除该小波系数。

进一步地，所述的步骤D包括以下步骤：

g：计算每层的通用阈值，即：

其中。λ为通用阈值，σ为噪声信号的标准差，N为小波系数的个数。

h：对高频系数进行软阈值处理，软阈值函数表达式为：

其中WT_f(a,b)为小波系数，为估计小波系数，λ为阈值。

进一步地，所述的步骤E包括以下步骤：

i：小波重构恢复去噪后的钟差一次差分数据，重构的计算公式为：

其中，Sf(a,b)为低频系数，WT_f(a,b)为小波系数，和分别对应于重构低通滤波器和重构高通滤波器。

与现有技术相比，本发明的有益效果在于：本发明提出一种基于WMAD的GNSS钟差数据预处理的方法，小波分析(Wavelet)由于具有多分辨率分析的特点，可以有效的除去数据中由白噪声引起的较小异常值，而中位数方法可以通过调节MAD的倍数n值来有效剔除数据中较大异常值。本发明分析了小波分析和中位数方法对应含义，并将两种方法融合(即WMAD)引用到GNSS卫星钟差数据预处理过程中，提高GNSS卫星钟差数据预处理的质量。首先，利用小波函数对钟差数据进行多层小波分解，得到分解后的小波系数。其次，利用MAD方法剔除较大的小波系数，再利用处理后的小波系数来计算各层的阈值并进行阈值化处理，最后，将处理后的小波系数进行小波重构得到预处理后钟差数据。因此，本发明结合两种方法的优点，提出一种小波分析与中位数法结合的钟差预处理方法，对钟差数据进行预处理，从而提高钟差数据的质量。

附图说明

图1为本发明所述的基于WMAD的GNSS卫星钟差数据预处理方法流程图；

图2是具体实施方式中PRN02钟差一次差分图；

图3是MAD预处理后一次差分数据显示图；

图4是小波阈值算法预处理后的一次差分数据显示图；

图5是WMAD预处理后一次差分数据显示图；

图6是MAD预处理、小波阈值算法以及本发明WMAD预处理三种方法预处理后的预报误差图。

具体实施方式

下面将结合实施例，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

参照图1，基于WMAD的GNSS卫星钟差数据预处理方法，包括以下步骤：

A：对卫星钟差进行一次差计算，得到钟差一次差分数据；

B：对GNSS卫星钟差一次差分数据进行小波分解，得到分解后小波系数；

C：利用MAD方法剔除较大的小波系数；

D：计算各层阈值，进行阈值化处理；

E：小波重构得到预处理后的钟差数据。

进一步地，所述的步骤A包括以下步骤：

a：首先设L＝{l(i),i＝1,2,…,n}为一组不同历元时刻的钟差值；

b：计算卫星钟差一次差分数据，即

ΔL(j)＝l(i+1)-l(i),j＝1,2,…n-1

其中，ΔL(j)为钟差一次差分数据。

进一步地，所述的步骤B包括以下步骤：

c：确定小波函数；

分别利用不同小波函数对钟差一次差分数据进行小波分解，得到分解后的小波系数，即

其中，WT_f(a,b)为小波变换系数，f(g)为卫星钟差一次差分数据，ψ_a,b(g)为小波母函数，a为伸缩因子，b平移因子(g＝0，1,2，…T-1,T为卫星钟差一次差数据个数)。

计算小波变换后小波系数的能量值，即

其中，E_s(W)为提取信号的能量值，N为小波系数的数量，W(s,t)为小波系数，s表示尺度参数，t表示信号幅值。获得最大能量值的小波函数为选择的最优小波函数

d：确定小波分解层数；

根据尺度与频率的对应关系，即：

其中f_s ^ps是尺度s对应的伪频率(pseudo frequency),f_c是小波对应的中心频率，Δt是采用周期。伪频率应覆盖有用信号f_sig的全部频率范围，即：

因此有用信号的最小频率(minf_sig)应该不小于尺度为L的伪频率，也就是：

离散小波变换的分解层数是对尺度s进行二进离散采样，因此L对应的分解层数j为：

2^j-1≤L≤2^j

e：对钟差数据进行小波分解得出小波系数。

进一步地，所述的步骤C包括以下步骤：

f：然后对每层的小波系数WT_f(a,b)与小波系数序列的中数(MED)m加上中位数(MAD)的n倍之和相比较，即若高频系数满足

|WT_f(a,b)|＞(m+n×MAD)

就认为该小波系数对应的就是较大异常值，同时剔除该小波系数。

进一步地，所述的步骤D包括以下步骤：

g：计算每层的通用阈值，即：

其中。λ为通用阈值，σ为噪声信号的标准差，N为小波系数的个数。

h：对高频系数进行软阈值处理，软阈值函数表达式为：

其中WT_f(a,b)为小波系数，为估计小波系数，λ为阈值。

进一步地，所述的步骤E包括以下步骤：

i：小波重构恢复去噪后的钟差一次差分数据，重构的计算公式为：

其中，Sf(a,b)为低频系数，WT_f(a,b)为小波系数，和分别对应于重构低通滤波器和重构高通滤波器。

为了验证本发明提出钟差预处理方法的有效性，本实施例使用IGS网站提供的GPS系统最终精密星历中30s间隔的精密卫星钟差数据进行试算分析。以GPS week1584第一天(2010年5月17号)的PRN02卫星钟差数据为例(其他卫星数据均可)，首先将钟差数据做一次差分，PRN02钟差一次差分如图2所示；然后在将一次差分数据分别利用MAD、小波阈值算法以及本发明提出的WMAD预处理策略进行预处理，用预处理后的数据分别通过BP神经网络建模预报卫星钟差。

星载原子钟正常情况是比较稳定的，相邻历元间钟差数据的数值变化相对较小，从图2上看PRN02卫星钟差数据受到噪声的影响相邻两个历元的钟差数据变化较大，进行一次差分后出现了很多峰值，如果用含有大量异常值的数据进行建模预报钟差数据，会大大影响预报的精度，所以在建模之前对钟差一次差分数据进行合理的预处理是很有必要的。

图3、图4和图5分别是MAD预处理、小波阈值算法以及本发明WMAD预处理三种方法对PRN02卫星钟差一次差分预处理后的效果图，从图3中可以很明显的看出，MAD方法对剔除相对较大的异常值效果明显，但是对相对较小的异常值却达不到满意的效果，这是因为钟差进行一次差分后变化量较小且数值变小，用一次差分数据的中位数作为阈值，往往相对较小的异常值不能被阈值发现，因此不能很好的剔除相对较小的异常值；从图4中可以明显看出，小波阈值算法经过频率域处理能够很好的将相对较小的异常值处理掉，但是不能准确的剔除相对较大异常值，这是因为由于较大异常值得影响，会使每层计算的噪声方差变小，从而计算的阈值会变小，不能探测出一些代表较大异常值得小波系数；而从图5中可以看出，本发明的WMAD预处理方法弥补了上面两种方法的不足，即有效剔除较大异常值又能很好的抑制较小的噪声误差，更好的净化数据。

为了更好的验证本发明提出方法的有效性，将三种方法预处理后的一次差数据通过BP神经网络建模对钟差进行预测，利用预报后的一次差数据与已知钟差数据对于相加便可得到所需时刻的预报钟差，最后，以IGS相应的精密钟差值作为基准，用均方根误差(RMS)和平均误差(Mean)评价预报的结果的精度；利用最大、最小误差之差的绝对值(Range)评价算法进行预报的稳定性。其中均方根误差计算公式为：

式中，t_i是IGS精密钟差；为钟差预报值。

为了验证WMAD预处理方法的有效性，使用没有进行预处理的一次差数据和MAD预处理、小波阈值算法以及本发明WMAD预处理三种方法预处理后的一次差数据前1440个历元(12h)进行BP神经网络建模,预报接下来的1440历元12个小时的钟差，比较分析预报的效果。表1给出了没有进行预处理及利用三种方法预处理后的预报统计值，预报误差图如图6所示。

表1预报结果统计表/ns

对比表1中使用没有进行预处理的一次差分数据和使用三种方法预处理一次差值钟差预报统计值，从预报结果的RMS和Mean值可以看出，经过三种方法预处理后的钟差预报精度均优于没有经过预处理钟差的预报精度，而根据预报结果的Range值则可看出，经过预处理后的数据预报稳定性也均优于没有预处理后钟差预报稳定性；经MAD和小波阈值算法两种方法预处理后钟差预报精度基本相当，预报稳定性上MAD方法优于小波阈值算法；而经本发明提出的WMAD预处理后的钟差预报精度和稳定性上都大大的优于前面的两种预处理方法。

从图6预报误差图也能看出，未经预处理一次差值预报的钟差误差最大，而且有较大的浮动；经MAD和小波阈值预处理后的钟差预报误差相比未经预处理一次差值预报的钟差误差减小，但MAD预处理后预报误差浮动较大，而小波阈值预处理后预报误差虽然浮动较小但出现相对较大误差稳定性不高；经WMAD预处理后的钟差误差值不仅较小，并且没有较大的浮动，具有很好的稳定性，验证了本发明方法的有效性。

Claims (6)

1.基于WMAD的GNSS卫星钟差数据预处理方法，其特征在于，包括以下步骤：

A：对GNSS卫星钟差进行一次差计算，得到钟差一次差分数据；

B：对钟差一次差分数据进行小波分解，得到分解后小波系数；

C：利用MAD方法剔除较大的小波系数；

D：计算各层阈值，进行阈值化处理；

E：小波重构得到预处理后的钟差一次差分数据。

2.如权利要求1所述的基于WMAD的GNSS卫星钟差数据预处理方法，其特征在于，所述的步骤A包括以下步骤：

a：首先设L＝{l(i),i＝1,2,…,n}为一组不同历元时刻的钟差值；

b：计算卫星钟差一次差分数据，即

ΔL(j)＝l(i+1)-l(i),j＝1,2,…n-1

其中，ΔL(j)为钟差一次差分数据。

3.如权利要求1所述的基于WMAD的GNSS卫星钟差数据预处理方法，其特征在于，所述的步骤B包括以下步骤：

c：确定小波函数；

分别利用不同小波函数对钟差一次差分数据进行小波分解，得到分解后的小波系数，即

其中，WT_f(a,b)为小波变换系数，f(g)为卫星钟差一次差分数据，ψ_a,b(g)为小波母函数，a为伸缩因子，b平移因子，其中g＝0，1,2，…T-1,T为卫星钟差一次差数据个数；

计算小波变换后小波系数的能量值，即

其中，E_s(W)为提取信号的能量值，N为小波系数的个数，W(s,t)为小波系数，s表示尺度参数，t表示信号幅值；获得最大能量值的小波函数为选择的最优小波函数；

d：确定小波分解层数；

根据尺度与频率的对应关系，即：

其中f_s ^ps是尺度s对应的伪频率ps(pseudo frequency),f_c是小波对应的中心频率，Δt是采用周期；伪频率应覆盖有用信号f_sig的全部频率范围，即：

因此有用信号的最小频率(minf_sig)应该不小于尺度s为L的伪频率，也就是：

离散小波变换的分解层数是对尺度s进行二进离散采样，因此L对应的分解层数j为：

2^j-1≤L≤2^j

e：对钟差数据进行小波分解得出小波系数。

4.如权利要求1所述的基于WMAD的GNSS卫星钟差数据预处理方法，其特征在于，所述的步骤C包括以下步骤：

f：然后对每层的小波系数WT_f(a,b)与小波系数序列的中数(MED)m加上中位数(MAD)的n倍之和相比较，即若小波系数满足

|WT_f(a,b)|＞(m+n×MAD)

就认为该小波系数对应的就是较大异常值，同时剔除该小波系数。

5.如权利要求1所述的基于WMAD的GNSS卫星钟差数据预处理方法，其特征在于，所述的步骤D包括以下步骤：

g：计算每层的通用阈值，即：

其中：λ为通用阈值，σ为噪声信号的标准差，N为小波系数的个数；

h：对高频系数进行软阈值处理，软阈值函数表达式为：

其中WT_f(a,b)为小波系数，为估计小波系数，λ为阈值。

6.如权利要求1所述的基于WMAD的GNSS卫星钟差数据预处理方法，其特征在于，所述的步骤E包括以下步骤：

i：小波重构恢复去噪后的钟差一次差分数据，重构的计算公式为：

其中，Sf(a,b)为低频系数，WT_f(a,b)为小波系数，和分别对应于重构低通滤波器和重构高通滤波器。