CN114510969A - 一种坐标时间序列的降噪方法 - Google Patents
一种坐标时间序列的降噪方法 Download PDFInfo
- Publication number
- CN114510969A CN114510969A CN202210083780.5A CN202210083780A CN114510969A CN 114510969 A CN114510969 A CN 114510969A CN 202210083780 A CN202210083780 A CN 202210083780A CN 114510969 A CN114510969 A CN 114510969A
- Authority
- CN
- China
- Prior art keywords
- sequence
- decomposition
- coordinate time
- value
- time sequence
- 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.)
- Pending
Links
- 238000000034 method Methods 0.000 title claims abstract description 50
- 238000000354 decomposition reaction Methods 0.000 claims abstract description 61
- 239000011159 matrix material Substances 0.000 claims abstract description 26
- 238000007781 pre-processing Methods 0.000 claims abstract description 11
- 238000012935 Averaging Methods 0.000 claims description 2
- 238000001514 detection method Methods 0.000 claims description 2
- 230000008030 elimination Effects 0.000 claims description 2
- 238000003379 elimination reaction Methods 0.000 claims description 2
- 238000012163 sequencing technique Methods 0.000 claims description 2
- 238000004088 simulation Methods 0.000 description 11
- 230000000694 effects Effects 0.000 description 10
- 238000012545 processing Methods 0.000 description 9
- 238000004364 calculation method Methods 0.000 description 8
- 238000002474 experimental method Methods 0.000 description 6
- 238000010586 diagram Methods 0.000 description 2
- 238000005259 measurement Methods 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 238000005070 sampling Methods 0.000 description 2
- 230000002123 temporal effect Effects 0.000 description 2
- 230000036962 time dependent Effects 0.000 description 2
- 102000005717 Myeloma Proteins Human genes 0.000 description 1
- 108010045503 Myeloma Proteins Proteins 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 238000000605 extraction Methods 0.000 description 1
- 238000001914 filtration Methods 0.000 description 1
- 230000007774 longterm Effects 0.000 description 1
- 238000012423 maintenance Methods 0.000 description 1
- 238000013508 migration Methods 0.000 description 1
- 230000005012 migration Effects 0.000 description 1
- 229940037201 oris Drugs 0.000 description 1
- 230000000737 periodic effect Effects 0.000 description 1
- 230000001932 seasonal effect Effects 0.000 description 1
- 238000010183 spectrum analysis Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
- G06F2218/02—Preprocessing
- G06F2218/04—Denoising
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Mathematical Optimization (AREA)
- Pure & Applied Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Computational Mathematics (AREA)
- General Engineering & Computer Science (AREA)
- Mathematical Analysis (AREA)
- Signal Processing (AREA)
- Computing Systems (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Algebra (AREA)
- Artificial Intelligence (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- Complex Calculations (AREA)
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为预处理后的坐标时间序列长度。
进一步地,所述分解分量的计算公式为:
附图说明
图1是本发明坐标时间序列降噪的流程图;
图2(a)是采用本发明方法对仿真时间序列数据降噪处理后的重构序列结果图;
图2(b)是采用本发明方法对仿真时间序列数据降噪处理后的重构成分结果图;
图3是采用本发明方法对实测数据进行降噪处理后的结果图。
具体实施方式
下面结合附图对本发明的具体实施方式作进一步地说明。
本发明提出一种坐标时间序列的降噪方法,其具体流程如图1所示。首先获取原始的坐标时间序列,并对其进行预处理;再将预处理后的坐标时间序列转换成轨迹矩阵,然后基于轨迹矩阵对坐标时间序列进行奇异值分解,得到分解分量;最后根据MIC准则确定重构分解分量,生成重构后的坐标时间序列作为降噪结果。本发明通过MIC准则确定重构分量,综合考虑了重构序列以及残差序列的MIC值,以保证重构后的坐标时间序列具有较强的非线性关系,同时兼顾更多的分解分量和更少的信号残留,实现科学有效的降噪。
步骤1.数据获取及预处理
本发明首先获取原始坐标时间序列数据,主要为GNSS、VLBI、SLR和DORIS等技术产生的坐标时间序列数据。
以GNSS为例,由于在卫星信号接收过程中会受到电离层延迟、多路径效应、卫星周跳、接收机误差、环境干扰信号等多方面的影响,所获取的原始坐标时间序列会存在粗差、阶跃及数据缺失的情况。为了保证采样的均匀性,当存在数据缺失时,可以通过插值得到缺失时刻的数据,例如采用最近邻法、三次多项式插值、奇异值迭代差值等方法实现缺失值内插;由于数据中的粗差、阶跃会影响后续处理精度,为了尽可能避免这一问题,需对数据中粗差和阶跃进行探测并剔除,例如可以通过最小二乘残差法、IQR稳健估计等方法进行粗差和阶跃的探测与剔除,这样可以事先剔除一部分噪声,更有利于后续处理;此外,由于奇异值分解对不同成分设置的窗口长度有所差异,因此为了简化坐标时间序列的成分,提高分解的准确度,需要对其进行了零均值化,可以通过线性拟合消除其线性趋势项,使得处理后的数据在0值上下波动、均匀分布。
步骤2.进行奇异值分解,求解分解分量
设经过预处理后的坐标时间序列为X=(x1,x2,x3,…,xN),选择一个合适的窗口长度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为ak构成的系数矩阵,ak为对应特征向量Ek的系数,Ek为轨迹矩阵的协方差阵中第k个特征值对应的特征向量,轨迹矩阵D的协方差矩阵Tx如公式(3)所示:
同时ak也可以看作是X经过滤波得到的结果,其第i个元素可以通过公式(4)计算得到:
通过上述公式(2)-(4),可以得到第k个分解分量的第i个值为:
步骤3.数据重构
本发明利用MIC准则确定用于数据重构的分量数量。对分解得到的M个分量分别计算其与时间之间的MIC值,具体计算方法如下:
将历元t和坐标时间序列数据X形成的数据点的集合(ti,xi)分布在二维空间(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),将剩余的分解分量构成残差序列;计算重构序列XP的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为预处理后的坐标时间序列的长度。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202210083780.5A CN114510969A (zh) | 2022-01-20 | 2022-01-20 | 一种坐标时间序列的降噪方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202210083780.5A CN114510969A (zh) | 2022-01-20 | 2022-01-20 | 一种坐标时间序列的降噪方法 |
Publications (1)
Publication Number | Publication Date |
---|---|
CN114510969A true CN114510969A (zh) | 2022-05-17 |
Family
ID=81549722
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202210083780.5A Pending CN114510969A (zh) | 2022-01-20 | 2022-01-20 | 一种坐标时间序列的降噪方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN114510969A (zh) |
Cited By (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN115272301A (zh) * | 2022-09-20 | 2022-11-01 | 江苏新世嘉家纺高新科技股份有限公司 | 一种基于机器人的筒子纱缺陷自动检测方法 |
CN116450711A (zh) * | 2023-06-20 | 2023-07-18 | 山东科技大学 | Gnss坐标时间序列数据流匹配方法 |
-
2022
- 2022-01-20 CN CN202210083780.5A patent/CN114510969A/zh active Pending
Cited By (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN115272301A (zh) * | 2022-09-20 | 2022-11-01 | 江苏新世嘉家纺高新科技股份有限公司 | 一种基于机器人的筒子纱缺陷自动检测方法 |
CN116450711A (zh) * | 2023-06-20 | 2023-07-18 | 山东科技大学 | Gnss坐标时间序列数据流匹配方法 |
CN116450711B (zh) * | 2023-06-20 | 2023-08-18 | 山东科技大学 | Gnss坐标时间序列数据流匹配方法 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN109243483B (zh) | 一种含噪频域卷积盲源分离方法 | |
CN114510969A (zh) | 一种坐标时间序列的降噪方法 | |
North et al. | Detection of forced climate signals. Part 1: Filter theory | |
CN110619296B (zh) | 一种基于奇异分解的信号降噪方法 | |
CN109523486B (zh) | 噪声环境下基于鲁棒压缩感知的多通道脑电信号重构方法 | |
CN109343060B (zh) | 基于深度学习时频分析的isar成像方法及*** | |
CN114253962B (zh) | 一种考虑非线性因素的区域格网速度场构建方法及*** | |
US20120109563A1 (en) | Method and apparatus for quantifying a best match between series of time uncertain measurements | |
CN111160317A (zh) | 微弱信号盲提取方法 | |
CN111079591B (zh) | 基于改进多尺度主成分分析的不良数据修复方法及*** | |
CN110118958B (zh) | 基于变分编码-解码网络的宽带雷达复回波去噪方法 | |
CN111142134B (zh) | 一种坐标时间序列处理方法及装置 | |
CN113655534B (zh) | 基于多线性奇异值张量分解核磁共振fid信号噪声抑制方法 | |
CN112817056B (zh) | 一种大地电磁信号去噪方法及*** | |
CN113139918B (zh) | 一种基于决策灰狼优化字典学习的图像重构方法 | |
CN111398912B (zh) | 基于张量低秩逼近的合成孔径雷达干扰抑制方法 | |
CN109272054B (zh) | 一种基于独立性的振动信号去噪方法及*** | |
Ge et al. | A new MCMC algorithm for blind Bernoulli-Gaussian deconvolution | |
CN108020324B (zh) | 单束激光叠加脉冲信号的滤波检测方法 | |
CN109188473B (zh) | 基于盲分离技术的北斗卫星微弱信号高精度快速捕获方法 | |
CN110687605A (zh) | 基于改进的k-svd算法在地震信号处理的算法分析应用 | |
CN110764135A (zh) | 不规则地震数据全频带重建方法 | |
CN117250585A (zh) | 基于不完全谱重构的多分量线性调频信号参数估计方法 | |
CN112711074B (zh) | 一种地震初至波的去噪方法及装置 | |
Kharkov et al. | Blind estimation of noise variance for 1D signal denoising |
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 |