CN114510969A - 一种坐标时间序列的降噪方法 - Google Patents

Abstract

本发明提出一种坐标时间序列的降噪方法，首先获取原始的坐标时间序列，并对其进行预处理；再将预处理后的坐标时间序列转换成轨迹矩阵，然后基于轨迹矩阵对坐标时间序列进行奇异值分解，得到分解分量并计算各分解分量与时间之间的MIC值；最后选择MIC值较大的前P个分量进行重构，重构后的坐标时间序列作为降噪结果，以实现对原始坐标时间序列的降噪。本发明通过MIC自适应准则确定用于重构的分解分量，综合考虑了分解分量、重构序列以及残差序列的MIC值，在保证重构后的坐标时间序列具有较强的非线性关系的同时兼顾更多的分解分量和更少的信号残留，实现了科学有效的降噪。

Description

一种坐标时间序列的降噪方法

技术领域

本发明涉及一种坐标时间序列的降噪方法，属于坐标时间序列降噪技术领域。

背景技术

坐标时间序列中包含着许多与板块构造运动、地表质量迁移等地球动力学过程等有关的重要信息，这些信息成分往往是时间相关的，能够为地球动力学研究提供丰富的数据基础。坐标时间序列中的时间相关性成分主要由以下几部分组成：由构造运动引起的测站坐标的长期变化、由地球物理效应引起的测站坐标的季节性变化、地震引起的震后形变以及其他与时间相关的坐标变化。而这些能够反映地球物理学过程的成分往往与噪声混叠在一起，甚至被噪声所掩盖，坐标时间序列中的噪声主要来自于观测误差以及数据处理策略不完善带来的误差等。因此，对坐标时间序列进行有效的降噪，能够实现时间相关性成分的提取，继而为地球动力学研究提供准确数据，同时坐标时间序列也是地球参考框架建立和维持的关键基础数据，科学的降噪方法有助于基准站坐标与速度的精确估计和地球参考框架的高精度维持。

现有的坐标时间序列的降噪方法有许多。其中有采用局部均值分解对坐标时间序列数据进行分解，将认定噪声分量剔除，但这会导致部分有用信息丢失；还有通过奇异谱分析对时间序列进行分解和重构，该方法作为一种数据驱动的非参数方法，能够在不需要任何先验信息的前提下，通过时间序列的分解和重构从包含噪声的时间序列中提取出有用的信息，进而实现降噪处理。但是奇异谱分析中重构成分的确定对于降噪效果具有很大的影响，重构时目前采用是根据贡献率大小来确定重构成分，且重构成分的确定往往是基于人的主观判断，缺乏科学的自适应准则，会出现有用信息丢失的情况，无法保证降噪效果。

发明内容

本发明的目的是提供一种坐标时间序列的降噪方法，以解决在现有技术对坐标时间序列数据降噪过程中因缺乏科学的自适应准则而导致的有用信息丢失的问题。

本发明提供了一种坐标时间序列的降噪方法，该方法包括以下步骤：

1)获取原始坐标时间序列；

2)将获取的坐标时间序列转换成相应的轨迹矩阵，基于转换的轨迹矩阵对坐标时间序列进行奇异值分解，得到分解分量；

3)计算各分解分量与时间之间的MIC值，得到各分解分量的MIC值，按照各分解分量MIC值大小对所述分解分量进行排序，选取其中MIC值较大的P个分解分量进行重构，生成重构序列，该重构序列为降噪后的坐标时间序列，以实现对原始坐标时间序列的降噪；

所述P的选取依据为：前P个分解分量生成的重构序列的MIC值大于第一设定阈值而前P+1个分解分量生成的重构序列的MIC值小于第一设定阈值或者前P个分解分量生成的重构序列的MIC值大于第一设定阈值且残差序列的MIC值小于第二设定阈值，其中，残差序列是选取出P个分解分量后剩余的分解分量所构成的序列。

本发明利用奇异值分解和分解分量重构实现了坐标时间序列的降噪处理，在进行分解分量重构时，通过引入MIC准则，综合考虑分解分量、重构序列以及残差序列的MIC值，保证了重构后的坐标时间序列具有较强的非线性关系，同时兼顾更多的分解分量和更少的信号残留，实现了科学有效的降噪。

进一步地，所述步骤3)中的第一设定阈值的取值范围为[0.8,1]；第二设定阈值的取值范围为[0,0.3]。

通过上述这种方式，可以保证重构序列的强时间相关性和残差序列的弱时间相关性。

进一步地，所述步骤2)奇异值分解前需对原始坐标时间序列数据进行预处理，预处理包括粗差和阶跃的探测与剔除以及缺失值的内插。

通过对原始坐标时间序列的粗差和阶跃探测和剔除，可以减少一部分由测量误差造成的噪声，能够提高数据质量，保证后续重构序列的精度；通过对缺失值进行内插，可以保证采样的均匀性。

进一步地，所述预处理还包括对原始坐标时间序列的零均值化处理。

为满足后续计算过程中对数据零均值性的要求，通过对原始坐标时间序列进行零均值化处理，使得处理后数据在0值处上下波动且均匀分布，消除线性趋势项，方便后续对坐标时间序列进行分析，保障后续运算顺利进行。

进一步地，所述轨迹矩阵为M行L列的矩阵，其中1＜M＜N/2，L＝N-M+1，N为预处理后的坐标时间序列长度。

进一步地，所述分解分量的计算公式为：

其中，

为第k个重建分量的第i个值，x_i为经过预处理后坐标时间序列中的第i个元素，D_i为轨迹矩阵，E^k为轨迹矩阵的协方差阵中第k个特征值对应的特征向量，a^k为对应特征向量E^k的系数，

是E^k中的第j个元素，

是a^k的第i个元素。

附图说明

图1是本发明坐标时间序列降噪的流程图；

图2(a)是采用本发明方法对仿真时间序列数据降噪处理后的重构序列结果图；

图2(b)是采用本发明方法对仿真时间序列数据降噪处理后的重构成分结果图；

图3是采用本发明方法对实测数据进行降噪处理后的结果图。

具体实施方式

下面结合附图对本发明的具体实施方式作进一步地说明。

本发明提出一种坐标时间序列的降噪方法，其具体流程如图1所示。首先获取原始的坐标时间序列，并对其进行预处理；再将预处理后的坐标时间序列转换成轨迹矩阵，然后基于轨迹矩阵对坐标时间序列进行奇异值分解，得到分解分量；最后根据MIC准则确定重构分解分量，生成重构后的坐标时间序列作为降噪结果。本发明通过MIC准则确定重构分量，综合考虑了重构序列以及残差序列的MIC值，以保证重构后的坐标时间序列具有较强的非线性关系，同时兼顾更多的分解分量和更少的信号残留，实现科学有效的降噪。

步骤1.数据获取及预处理

本发明首先获取原始坐标时间序列数据，主要为GNSS、VLBI、SLR和DORIS等技术产生的坐标时间序列数据。

以GNSS为例，由于在卫星信号接收过程中会受到电离层延迟、多路径效应、卫星周跳、接收机误差、环境干扰信号等多方面的影响，所获取的原始坐标时间序列会存在粗差、阶跃及数据缺失的情况。为了保证采样的均匀性，当存在数据缺失时，可以通过插值得到缺失时刻的数据，例如采用最近邻法、三次多项式插值、奇异值迭代差值等方法实现缺失值内插；由于数据中的粗差、阶跃会影响后续处理精度，为了尽可能避免这一问题，需对数据中粗差和阶跃进行探测并剔除，例如可以通过最小二乘残差法、IQR稳健估计等方法进行粗差和阶跃的探测与剔除，这样可以事先剔除一部分噪声，更有利于后续处理；此外，由于奇异值分解对不同成分设置的窗口长度有所差异，因此为了简化坐标时间序列的成分，提高分解的准确度，需要对其进行了零均值化，可以通过线性拟合消除其线性趋势项，使得处理后的数据在0值上下波动、均匀分布。

步骤2.进行奇异值分解，求解分解分量

设经过预处理后的坐标时间序列为X＝(x₁,x₂,x₃,…,x_N)，选择一个合适的窗口长度M，将X转换为M行L列的轨迹矩阵D，如公式(1)所示，其中，M应当满足1＜M＜N/2，但在一般情况下M不宜超过时间序列长度的1/3；L＝N-M+1，N为预处理后的坐标时间序列长度。

基于轨迹矩阵D对预处理后的坐标时间序列进行奇异值分解，即可得到分解分量，具体计算过程如下：

由式(1)可知D的行和列是X的子序列，因此其可以在一组标准正交基下展开：得到如公式(2)所示：

式中，A为a^k构成的系数矩阵，a^k为对应特征向量E^k的系数，E^k为轨迹矩阵的协方差阵中第k个特征值对应的特征向量，轨迹矩阵D的协方差矩阵T_x如公式(3)所示：

式中，

求解T_x可得到特征值和特征向量E^k。

同时a^k也可以看作是X经过滤波得到的结果，其第i个元素可以通过公式(4)计算得到：

式中，

是E^k中的第j个元素，

是a^k的第i个元素。

通过上述公式(2)-(4)，可以得到第k个分解分量的第i个值为：

其中，

为第k个重建分量的第i个值，D_i为轨迹矩阵。

步骤3.数据重构

本发明利用MIC准则确定用于数据重构的分量数量。对分解得到的M个分量分别计算其与时间之间的MIC值，具体计算方法如下：

将历元t和坐标时间序列数据X形成的数据点的集合(t_i,x_i)分布在二维空间(R,S)中，使用m乘以n的网格划分数据空间，将落在第(r,s)格子中的数据点的频率作为P(r,s)的估计，即

然后利用公式(7)计算得到时间序列X的MIC值，

其中网格的分辨率限制为m×n＜B，其中B设置为数据量的0.6次方。

在计算得到M个分解分量的MIC值后，将分解分量按照MIC值从大到小进行排列，然后选取MIC值最大的前P个分解分量进行累加得到重构序列

其中P为重构的阶次(0＜P≤M)，将剩余的分解分量构成残差序列；计算重构序列X_P的MIC值

和残差序列的MIC值

其中重构序列的MIC值

和残差序列的MIC值

的计算方法与上述分解分量的MIC值计算方法一致，只是将输入数据转成重构序列和残差序列。

根据上式计算得到的

和

可以确定P的个数，具体P的选取依据为：

大于第一设定阈值而

小于第一设定阈值或者

大于第一设定阈值且

小于第二设定阈值；其中，为了保证重构序列的强时间相关性和残差序列的弱时间相关性，第一设定阈值的取值范围为[0.8,1]，第二设定阈值的取值范围为[0,0.3]，作为其他实施方式，第一阈值和第二阈值可以根据提取分量相关性强弱的需求及获取的坐标时间序列数据的质量进行设定。

根据上述P的选取依据确定了P的个数后，按照分解分量MIC值的大小对较大的前P个分解分量进行重构，得到的重构序列就是坐标时间序列的降噪后序列，本发明通过上述方式可以实现对原始坐标时间序列的降噪处理。

为了进一步说明本发明方法对坐标时间序列的降噪效果，分别利用仿真数据和实测数据对本发明方法进行了验证。其中仿真数据由若干个周期项(见表1)以及不同量级的白噪声(White Noise,WH)和闪烁噪声(Flicker Noise,FL)构成。

表1 仿真时间序列中的周期项及其参数

仿真数据实验：

为了验证本发明方法在不同噪声大小和成分下的性能，共设置了两组实验，第一组实验控制噪声成分不变(WH和FL分别为40％和60％)，分别对信噪比为10dB、5dB、0dB、-5dB和-10dB的仿真时间序列进行降噪；第二组实验控制噪声强度不变(信噪比为0dB)，分别对不同噪声成分(WH+FL：100％+0，80％+20％，60％+40％，40％+60％，20％+80％，0+100％)的仿真时间序列进行降噪分析。对信噪比为5dB，噪声成分为40％WH和60％FL的仿真时间序列的降噪结果进行了展示，如图2(a)和2(b)所示，可以看到采用本发明方法很好地实现了仿真时间序列的降噪处理。同时还将本发明方法与传统SSA方法进行了对比，表2展示了相同噪声成分不同信噪比下不同方法对仿真时间序列的降噪效果，表3展示了相同信噪比不同噪声成分下不同方法对仿真时间序列的降噪效果。

表2：

表3：

从仿真实验的结果来看，相比于传统SSA通过设置截止特征值贡献率自动选取重构成分的降噪方法，本发明的降噪方法具有明显的优势，在不同的噪声成分和噪声强度下其降噪后信号与真实时间序列之间的相关性以及信噪比都得到了较大的提升。而本发明方法与传统SSA通过手动选取重构成分的降噪方法具有相当的降噪效果，同时能够实现自适应的降噪，避免了主观因素的影响，而且通过引入MIC准则使得降噪效果更具科学性。

实测数据实验：

在实测数据的实验中，选取了IGS测站中的ZIMM测站作为实验对象验证本发明方法的降噪效果。实验结果如图3所示(第一、第二阈值分别为0.999和0.01)。

从图3中可以看出，本发明方法实现了对实测坐标时间序列的降噪，并且经过计算，降噪后的时间序列MIC值为0.999，残差序列的MIC值为0.07，重构序列的特征值贡献率可以达到83％，在保证重构后的坐标时间序列具有较强的非线性关系的同时兼顾更多了的分解分量和更少的信号残留，实现了科学有效的降噪。

Claims (6)

1.一种坐标时间序列的降噪方法，其特征在于，该方法包括以下步骤：

1)获取原始坐标时间序列；

2)将获取的坐标时间序列转换成相应的轨迹矩阵，基于转换的轨迹矩阵对坐标时间序列进行奇异值分解，得到分解分量；

3)计算各分解分量与时间之间的MIC值，得到各分解分量的MIC值，按照各分解分量MIC值大小对所述分解分量进行排序，选取其中MIC值较大的P个分解分量进行重构，生成重构序列，该重构序列为降噪后的坐标时间序列，以实现对原始坐标时间序列的降噪；

所述P的选取依据为：前P个分解分量生成的重构序列的MIC值大于第一设定阈值而前P+1个分解分量生成的重构序列的MIC值小于第一设定阈值或者前P个分解分量生成的重构序列的MIC值大于第一设定阈值且残差序列的MIC值小于第二设定阈值，其中，残差序列是选取出P个分解分量后剩余的分解分量所构成的序列。

2.根据权利要求1所述的坐标时间序列的降噪方法，其特征在于，所述步骤3)中的第一设定阈值的取值范围为[0.8,1]；第二设定阈值的取值范围为[0,0.3]。

3.根据权利要求1所述的坐标时间序列的降噪方法，其特征在于，所述步骤2)奇异值分解前需对原始坐标时间序列数据进行预处理，预处理包括粗差和阶跃的探测与剔除以及缺失值的内插。

4.根据权利要求3所述的坐标时间序列的降噪方法，其特征在于，所述预处理还包括对原始坐标时间序列的零均值化处理。

5.根据权利要求3或4所述的坐标时间序列的降噪方法，其特征在于，所述轨迹矩阵为M行L列的矩阵，其中1＜M＜N/2，L＝N-M+1，N为预处理后的坐标时间序列的长度。

6.根据权利要求1或3或4所述的坐标时间序列的降噪方法，其特征在于，分解分量的计算公式为：

其中，

为第k个重建分量的第i个值，x_i为经过预处理后坐标时间序列中的第i个元素，D_i为轨迹矩阵，E^k为轨迹矩阵的协方差阵中第k个特征值对应的特征向量，a^k为对应特征向量E^k的系数，

是E^k中的第j个元素，

是a^k的第i个元素。

Images