去除GPS坐标时间序列中有色噪声的方法与流程

本发明涉及gps数据处理领域，尤其涉及一种去除gps坐标时间序列中有色噪声的方法。

背景技术：

gps基准站坐标时间序列被广泛地应用于大地测量和地球动力学研究。它可以为研究地壳运动的时空变化、结构变形特征和其他科学问题提供基础数据。由于各种因素包括地球物理环境以及与gps技术相关的系统误差，导致gps站坐标时间序列中含有白噪声以及有色噪声(误差)。由于gps坐标时间序列中的噪声影响了测站坐标的稳定性和可靠性，导致对一些地球物理现象的错误解释，因此有必要采取措施削弱gps坐标时间序列的噪声强度。然而，目前的去噪方法还不能剔除存在于低频中的有色噪声，并且即使削弱了噪声的强度，依然不能准确获得噪声振幅的估值，尤其是残留的有色噪声。

一些学者提出基于小波变换去除gps坐标时间序列中的噪声。常用的方法为离散小波变换，该方法基于mallat算法对gps坐标时间序列进行多分辨率分解，然后对获得的高频的小波变换系数进行阈值量化，最后通过对处理后的小波系数进行重构获得去噪后的坐标时间序列。然而该方法在设置阈值时仅考虑了白噪声，对于削弱坐标时间序列中的高频噪声(白噪声)具有较好的效果，对于存在于坐标时间序列中的低频噪声如闪烁噪声、随机漫步噪声，去噪效果并不理想。

[wu,h.,li,k.,shi,w.,clarke,k.c.,zhang,j.,&li,h.(2015).awavelet-basedhybridapproachtoremovetheflickernoiseandthewhitenoisefromgpscoordinatetimeseries.gpssolutions,19(4),511-523]提出基于小波熵和小波变换系数的模极大值去噪。该方法根据模拟的白噪声的小波系数以及坐标时间序列的小波系数的小波熵，获得闪烁噪声的小波系数，然后将其重构为闪烁噪声并从原始坐标时间序列中剔除。最后利用一般的离散小波变换削弱剩余的白噪声。该方法需要提前已知白噪声与余下信号的信噪比，以此获得白噪声的方差，在此基础上才能模拟出白噪声序列并获得其小波变换系数；此外，该方法通过白噪声的小波熵获得有色噪声的小波熵时，忽略了时间序列中周期信号的小波系数模极大值的影响以及存在于低频中的白噪声。因此该方法只适用于较短的坐标时间序列，否则很可能会对周期信号造成影响。

[张恒璟,程鹏飞.(2014).基于eemd的gps高程时间序列噪声识别与提取.大地测量与地球动力学,34(2),79-83]和[张双成,何月帆,李振宇,侯晓伟,瞿伟,南阳.(2017).emd用于gps时间序列降噪分析.大地测量与地球动力学,37(12),1248-1252]提出基于经验模式分解(empiricalmodedecomposition，emd)法去除gps坐标时间序列中的噪声。该方法类似于离散小波变换，将原始信号分解为一系列频率由高到低的本征函数模态分量(intrinsicmodefunction，imf)，残差项与邻近的低频imf分量合并，通过imf分量方差贡献率等指标，判断序列的非线性趋势项。接近趋势项的中低频imf分量是明显的信号周期运动项，高频部分可能是暂时无法识别的信号高频周期项，或者是信号中包含的噪声信息。利用该方法可以有效地提取gps坐标时间序列中的非线性信号。然而，有色噪声的功率集中于低频，而白噪声的功率在整个功率谱中都是维持在相同的水平，因此，利用emd方法同样无法剔除存在于低频中的白噪声以及有色噪声。

[王解先,连丽珍,沈云中.(2013).奇异谱分析在gps站坐标监测序列分析中的应用.同济大学学报:自然科学版,(2),282-288.]利用奇异谱分析法对gps测站坐标时间序列进行处理。该方法根据所观测到的时间序列构造出轨迹矩阵，并对轨迹矩阵进行分解、重构，从而提取出代表原时间序列不同成分的信号，如长期趋势信号、周期信号、噪声信号等，从而对时间序列的结构进行分析。实验结果表明，通过奇异谱分析法重构后的序列比原始观测序列更为光滑，有明显的降噪效果。然而该方法依然不能将低频的有色噪声与周期信号有效分离。

技术实现要素：

本发明的目的在于提供一种去除gps坐标时间序列中有色噪声的方法，旨在用于解决现有的去除gps坐标时间序列中的噪声的方法并不能将有色噪声从gps坐标时间序列中有效去除的问题。

本发明是这样实现的：

本发明提供一种去除gps坐标时间序列中有色噪声的方法，包括以下步骤：

s1，获得gps基准站坐标时间序列的残差序列，估计出残差序列中白噪声的方差；

s2，根据白噪声的方差模拟出白噪声；

s3，将模拟白噪声以及残差序列同时进行小波分解；

s4，获取模拟白噪声的小波系数的信息熵以及残差序列的小波系数的信息熵；

s5，根据模拟白噪声的小波系数的信息熵以及残差序列的小波系数的信息熵，获得有色噪声的小波系数的信息熵：

s6，根据有色噪声的小波系数的信息熵获得有色噪声小波系数；

s7，将有色噪声的小波系数重构，获得有色噪声；

s8，将gps坐标时间序列减去获得的有色噪声，最终得到不受有色噪声影响的时间序列。

进一步地，所述步骤s1具体包括：

获得一组gps基准站干净的坐标时间序列，根据基准站运动模型以及坐标时间序列获得坐标残差序列，利用极大似然估计法分析坐标时间序列的最佳噪声组合模型，并估计出白噪声的方差。

进一步地，所述步骤s3中，分别利用离散小波变换和小波包变换同时对模拟白噪声以及残差序列进行小波分解。

进一步地，所述步骤s4具体包括：

利用离散小波变换时，根据以下公式分别计算模拟白噪声和坐标残差序列的小波高频系数以及最后一层低频系数的香农熵；利用小波包变换时，根据以下公式分别计算最后一层每个频段模拟白噪声以及坐标残差序列的小波系数的香农熵：

其中，si代表信号s在一个正交小波包基上的投影系数，e表示熵。

进一步地，所述步骤s5具体包括：

利用离散小波变换时，根据以下公式获得有色噪声小波高频系数以及最后一层低频系数的香农熵；利用小波包变换时，根据以下公式获得有色噪声在最后一层每个频段小波系数的香农熵：

其中，表示有色噪声小波系数香农熵，表示坐标残差序列的小波系数香农熵，表示白噪声小波系数香农熵。

进一步地，所述步骤s6具体包括：

利用离散小波变换时，根据以下公式获得有色噪声的高频以及最后一层低频小波系数；利用小波包变换时，根据以下公式获得有色噪声在最后一层每个频段的小波系数：

其中，和分别为坐标残差序列、有色噪声以及白噪声的小波系数；j表示小波系数所在的层数，k表示平移时间。

进一步地，所述步骤s7具体包括：

利用离散小波变换时，根据有色噪声的小波系数利用离散小波逆变换获得有色噪声；利用小波包变换时，根据有色噪声的小波系数利用小波包逆变换获得有色噪声，计算利用离散小波变换和小波包变换获得的有色噪声在每个时间点的均值，获得新的有色噪声。

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

本发明提供的这种去除gps坐标时间序列中有色噪声的方法，第一，该方法有效去除了gps坐标时间序列中有色噪声的影响。对于剩余的仅受白噪声影响的坐标时间序列，利用普通最小二乘即可获得准确的测站运动参数。第二，该方法仅对gps坐标残差序列进行处理，保留了gps坐标时间序列的非线性运动特征。第三，该方法同时采用离散小波变换以及小波包变换对gps坐标残差序列进行分解，获得的有色噪声顾及了其在不同频带(包括高频)上的特征，与原始的有色噪声在时域以及频域上更加接近。

附图说明

图1为本发明实施例提供的一种去除gps坐标时间序列中有色噪声的方法的流程图；

图2为离散小波分解树(左)及系小波系数所占频带宽度(右)；

图3为信号4层小波包分解的小波树(a)及最后一层小波系数频带宽度(b)；

图4为去除gps坐标时间序列中有色噪声流程图；

图5为kmin站坐标时间序列，从上到下依次为n、e、u分量的时间序列；

图6为2004年12月-2008年5月kmin站垂向残差序列及其双对数功率谱；

图7为2004年12月-2008年5月去除有色噪声后的kmin站垂向残差序列的双对数功率谱。

具体实施方式

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

如图1所示，本发明实施例提供一种去除gps坐标时间序列中有色噪声的方法，其特征在于，包括以下步骤：

s1，获得gps基准站坐标时间序列的残差序列，估计出残差序列中白噪声的方差(或振幅)；

s2，根据白噪声的方差模拟出白噪声；

s3，将模拟白噪声以及残差序列同时进行小波分解；

s4，获取模拟白噪声的小波系数的信息熵以及残差序列的小波系数的信息熵；

s5，根据模拟白噪声的小波系数的信息熵以及残差序列的小波系数的信息熵，获得有色噪声的小波系数的信息熵：

s6，根据有色噪声的小波系数的信息熵获得有色噪声小波系数；

s7，将有色噪声的小波系数重构，获得有色噪声；

s8，将gps坐标时间序列减去获得的有色噪声，最终得到不受有色噪声影响的时间序列。

下面对上述各步骤进行详细说明。

作为实施方式之一，所述步骤s1具体包括：获得一组gps基准站干净的坐标时间序列，根据基准站运动模型以及坐标时间序列获得坐标残差序列，利用极大似然估计法分析坐标时间序列的最佳噪声组合模型，并估计出白噪声的方差。

作为优选实施例，所述步骤s3中，为使gps坐标残差序列在频率上能够进行更加精细的分解，本发明采用了小波包变换，而且，为了使提取出的有色噪声与原有色噪声更加接近，本发明还采用了小波离散变换，即，本发明实施例分别利用小波包变换以及小波离散变换同时对模拟白噪声以及残差序列进行小波分解。一般分解为4或5层，分解层数视残差序列的时间跨度而定。

下面简要介绍这两种分解方法。

对小波母函数进行平移和伸缩可以得到一个小波基函数集，

式(1)中，为小波母函数，为小波基函数，其中a称为尺度伸缩因子，简称为尺度因子；b称为时间平移因子，简称平移因子，r为实数集。若尺度伸缩因子a和时间平移因子b连续变化，称为连续小波函数。通常取：将其代入(1)式可得：

此时离散小波变换公式为：

其中fdwt(a,b)为离散小波变换系数，*表示共轭。一维离散小波变换实现的算法一般是mallat算法，即先对较大的信号进行小波变换，获得低频系数和高频系数，再选取低频系数在原尺度的1/2尺度上再进行小波变换。

图2为对信号进行四层离散小波分解的树形图(左图)以及各层小波系数所占频带宽度(右图)。可以看出，对该信号进行4层离散小波分解，可以得到5个小波分解系数(4个高频小波系数和一个低频小波系数)。从第一层开始，每层小波系数所占频带宽度呈指数递减。因此，离散小波变换虽然可以对信号进行有效的时频分解，但其尺度按照二进制变化，即对信号频带进行指数等间隔划分，在高频频率分辨率较差。由于实际上每个频率所对应的有色噪声与含噪信号的小波熵的比不同，并且有色噪声的功率集中于低频，因此仅仅基于离散小波变换的高频小波系数获得的有色噪声仍与含噪信号中实际的有色噪声有较大差距。

从滤波器的角度，小波包变换的实现和离散小波变换没有本质区别，只是在原有的基础上按同样的方法对低频小波系数进行分解。所以其实现的方法与离散小波变换相同，但比离散小波变换分析的更为精细。

图3为信号4层小波包分解的小波树(a)及最后一层小波系数频带宽度(b)。由图3可以看出，对信号做n层分解，可以获得包括低频以及高频系数在内的2n个频带宽度相同的小波系数。

对gps残差序列进行小波变换后，本发明通过信息熵度量有色噪声的小波系数与残差序列的小波系数的比，在此基础上获得有色噪声的小波系数，进而通过小波逆变换获得有色噪声。

在信息论中，熵是用来度量信息规律性的概念，熵越小，则信息的规律性越强。熵具备可加性，即si代表信号s在一个正交小波包基上的投影系数，e表示熵，则e是一个由每个正交基上的系数的某种变换叠加起来的值，即且e(0)＝0。常用的熵有香农熵、p范数熵、对数能量熵、阈值熵、sure熵等。本发明实施例主要采用香农熵，其定义如下，约定0log(0)＝0：

式中对数一般取2为底，单位为比特。假设r(t)＝w(t)+x(t),其中w(t)为白噪声，x(t)为有色噪声，r(t)为观测信号。由于小波变换是线性的，因此有：

其中，和分别为观测值、有色噪声以及白噪声的小波系数；j表示小波系数所在的层数，k表示平移时间。假设w(t)为方差δ^w的白噪声，则其小波系数绝对值平方的均值为：

可以看出，随着分解层数的增加而减小。根据香农熵的定义，可得白噪声的小波系数的香农熵：

可以看出，白噪声小波系数的香农熵只与小波基方差δ^w以及小波系数所在层数j有关。因此，虽然不能获得观测序列中的白噪声，但可以模拟出具有相同方差的白噪声。然后将模拟白噪声在各层的小波系数香农熵替换实际白噪声的香农熵。同理，有色噪声的小波系数的香农熵公式为：

观测序列r的小波系数的香农熵为：

由于白噪声和有色噪声并不相关，因此，因此，上式可变为：

进而根据模拟白噪声的小波熵以及观测序列的小波熵，获得有色噪声的小波熵：

根据公式(5)和(11)可以获得有色噪声小波系数：

通过对有色噪声的小波系数的重构，获得有色噪声x：

基于上述原理，作为实施方式之一，所述步骤s4具体包括：

利用离散小波变换时，根据以下公式(4)分别计算模拟白噪声和坐标残差序列的小波高频系数以及最后一层低频系数的香农熵，即在离散小波变换下模拟白噪声和坐标残差序列的小波系数香农熵；利用小波包变换时，根据以下公式分别计算最后一层每个频段模拟白噪声以及坐标残差序列的小波系数的香农熵，即在小波包变换下模拟白噪声和坐标残差序列的小波系数香农熵：

其中，si代表信号s在一个正交小波包基上的投影系数，e表示熵，约定0log(0)＝0。

作为实施方式之一，所述步骤s5具体包括：

利用离散小波变换时，根据以下公式(11)获得有色噪声小波高频系数以及最后一层低频系数的香农熵，即在离散小波变换下有色噪声的小波系数香农熵；利用小波包变换时，根据以下公式获得有色噪声在最后一层每个频段小波系数的香农熵，即在小波包变换下有色噪声的小波系数香农熵：

其中，表示有色噪声小波系数香农熵，表示坐标残差序列的小波系数香农熵，表示白噪声小波系数香农熵。

作为实施方式之一，所述步骤s6具体包括：

利用离散小波变换时，根据以下公式(12)获得有色噪声的高频以及最后一层低频小波系数，即在离散小波变换下有色噪声的小波系数；利用小波包变换时，根据以下公式获得有色噪声在最后一层每个频段的小波系数，即在小波包变换下有色噪声的小波系数：

其中，和分别为坐标残差序列、有色噪声以及白噪声的小波系数；j表示小波系数所在的层数，k表示平移时间。

作为实施方式之一，所述步骤s7具体包括：

利用离散小波变换时，根据有色噪声的小波系数利用离散小波逆变换获得有色噪声；利用小波包变换时，根据有色噪声的小波系数利用小波包逆变换获得有色噪声，计算利用离散小波变换和小波包变换获得的有色噪声在每个时间点的均值，获得新的有色噪声。

进一步地，具体可根据以下公式(13)来获得有色噪声x：

本发明实施例分别利用离散小波变换、小波包变换对gps坐标残差序列进行分解，获得其在不同频带上的小波系数。根据每个频带上有色噪声与残差序列的小波系数香农熵之比，获取有色噪声的小波系数，进而在时域上通过小波逆变换获得有色噪声的序列。对于基于这两种方法获取的有色噪声，对其取均值，获得与原有色噪声在频域和时域上较为相似的有色噪声，并将其剔除，最终得到仅受白噪声影响的时间序列。

与以往的gps坐标时间序列去噪方法相比，该方法具有以下优点：第一，该方法有效去除了gps坐标时间序列中有色噪声的影响。对于剩余的仅受白噪声影响的坐标时间序列，利用普通最小二乘即可获得准确的测站运动参数。第二，该方法仅对gps坐标残差序列进行处理，保留了gps坐标时间序列的非线性运动特征。第三，该方法同时采用离散小波变换以及小波包变换对gps坐标残差序列进行分解，获得的有色噪声顾及了其在不同频带(包括高频)上的特征，与原始的有色噪声在时域以及频域上更加接近。

下面结合一个具体的例子对上述方法进行说明。

本实施例是剔除中国大陆构造环境监测网络(简称陆态网)kmin站u(up)分量时间序列中的有色噪声。图4为该方法的流程图，该方法的具体步骤如下：

(1)从中国地震局gnss产品服务平台获取kmin站1999-2018年的坐标时间序列(时间序列图图5所示)。可以看出，在2004年末和2008年5月有地震发生，致使地震前后的时间序列有明显的偏移。根据地震发生的时间，将kmin站时间序列分为三段：1999年3月-2004年12月,2004年12月-2008年5月,2008年5月-2018年3月。由图可以看出，2004年末发生的地震使kmin站垂直分量的时间序列产生的偏移较为显著，为避免噪声分析结果受地震的影响，本发明对该站2004年12月-2008年5月垂直分量的时间序列进行分析。估计并去除这段时间序列中由于更换天线或者地震所引起的明显的“跳变”。然后，利用被称为三倍中误差准则剔除粗差，最终获得干净的垂向时间序列。

(2)根据gps单站单分量的运动模型(公式(14))及其坐标时间序列，在仅考虑白噪声的情况下，利用一般最小二乘分别估计了测站在这段时间内u分量的运动速度、周年以及半周年信号的振幅以及白噪声的振幅(标准差)。

公式(14)中，y(ti)为gps单站单分量的坐标时间序列，x⁽¹⁾和x⁽²⁾分别为初始位置和速率，q-1为时间序列中包含的周期信号数，x^(2k-1)和x^(2k)是调和函数的系数，用于描述周期信号的振幅，vi为白噪声与有色噪声的线性组合，为误差。由于噪声模型并不会影响参数的估值，因此在该步骤中假设vi仅为白噪声。将运动参数的估值以及坐标时间序列代入公式(14)，获得这段时间序列的残差序列。该段时间序列的残差序列的功率谱如图6所示。可以看出，去除有色噪声之前，残差序列的功率集中在低频，这是因为有色噪声的功率也集中在低频，且闪烁噪声的噪声振幅估值远大于白噪声。

(3)根据公式(14)利用极大似然估计法计算分析这段时间序列的最佳噪声组合模型及其噪声振幅，结果为白噪声+闪烁噪声。

(4)根据白噪声的振幅估值，利用matlab软件生成具有相同噪声振幅的高斯白噪声。

(5)分别利用离散小波以及小波包对该站垂向的残差序列进行小波变换，并且用相同的方法对生成的白噪声进行小波变换。

(6)分别计算在离散小波变换以及小波包变换下，小波变换系数的香农熵。对于离散小波变换，根据公式(4)计算白噪声和残差序列的小波高频系数以及最后一层低频系数的香农熵；对于小波包变换，根据公式(4)分别计算最后一层每个频段白噪声以及残差序列的小波系数的香农熵。

(7)计算闪烁噪声的香农熵。根据公式(11)，利用白噪声以及残差序列的香农熵计算出闪烁噪声的香农熵。对于离散小波变换，获得闪烁噪声在小波高频系数以及最后一层低频系数的香农熵；对于小波包变换，获得闪烁噪声在最后一层每个频段小波系数的香农熵。

(8)计算闪烁噪声的小波系数。根据公式(12)，利用闪烁噪声小波系数的香农熵与残差序列的香农熵的比，计算出闪烁噪声的小波系数。对于离散小波变换，获得闪烁噪声的高频以及最后一层低频小波系数；对于小波包变换，获得闪烁噪声在最后一层每个频段的小波系数。

(9)根据闪烁噪声的小波系数，利用离散小波逆变换和小波包逆变换，分别获得这两种方法下的闪烁噪声。

(10)获得计算这两种方法下获得的闪烁噪声在每个时间点的均值，将其构成与残差序列中原有的闪烁噪声较为一致的闪烁噪声。

(11)将(10)中获得的闪烁噪声从残差序列中剔除，最终获得不受闪烁噪声影响的时间序列。剔除闪烁噪声后的残差序列的其功率谱，如图7所示。可以看出，去除有色噪声之后，残差序列的功率集中在高频。在频率大于7cpy的范围内，功率趋于平稳，而频率低于7cpy的功率基本为0。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。