CN113032713B - 一种基于加权多通道奇异谱的数据序列信号提取方法 - Google Patents

一种基于加权多通道奇异谱的数据序列信号提取方法 Download PDF

Info

Publication number
CN113032713B
CN113032713B CN202110212835.3A CN202110212835A CN113032713B CN 113032713 B CN113032713 B CN 113032713B CN 202110212835 A CN202110212835 A CN 202110212835A CN 113032713 B CN113032713 B CN 113032713B
Authority
CN
China
Prior art keywords
sequence
signal
data sequence
matrix
mssa
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.)
Active
Application number
CN202110212835.3A
Other languages
English (en)
Other versions
CN113032713A (zh
Inventor
沈云中
王奉伟
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Tongji University
Original Assignee
Tongji University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Tongji University filed Critical Tongji University
Priority to CN202110212835.3A priority Critical patent/CN113032713B/zh
Publication of CN113032713A publication Critical patent/CN113032713A/zh
Application granted granted Critical
Publication of CN113032713B publication Critical patent/CN113032713B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • 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
    • 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
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02D—CLIMATE CHANGE MITIGATION TECHNOLOGIES IN INFORMATION AND COMMUNICATION TECHNOLOGIES [ICT], I.E. INFORMATION AND COMMUNICATION TECHNOLOGIES AIMING AT THE REDUCTION OF THEIR OWN ENERGY USE
    • Y02D30/00—Reducing energy consumption in communication networks
    • Y02D30/70—Reducing energy consumption in communication networks in wireless communication networks

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • General Physics & Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Computing Systems (AREA)
  • Complex Calculations (AREA)

Abstract

本发明涉及一种基于加权多通道奇异谱的数据序列信号提取方法，1、获取多变量数据序列及先验误差序列，检验数据序列，判断其是否存在缺值，并进行中心化处理；2、对处理后的多变量数据序列进行信号估计；3、在估计信号的基础上获取多变量数据序列的噪声；4、根据先验误差序列构建权因子，结合步骤2得到的估计信号，生成新的数据序列；5、对新的数据序列再次进行信号估计，重复步骤2～4，提取新的信号；6、计算相邻两次提取信号的差值序列，判断其是否满足收敛条件，若满足，输出最后一次估计信号，否则，继续进行循环迭代直到信号结果收敛为止。与现有技术相比，本发明具有不破坏信号结构、加权方式效果更好且可靠性更高等优点。

Description

一种基于加权多通道奇异谱的数据序列信号提取方法
技术领域
本发明涉及时间序列数据分析及处理领域，尤其是涉及一种基于加权多通道奇异谱的数据序列信号提取方法。
背景技术
对于多变量观测数据序列xl(t),(t＝1,2,…,N；l＝1,2,…,L)，其中N表示序列长度，L为通道数，通常将观测数据简单的划分为信号与噪声两部分。目前针对噪声的滤除方法有很多种，如多元小波分析、主成分分析、以及多通道奇异谱分析(MultichannelSingular Spectrum Analysis，MSSA)等一些方法。MSSA与其他方法相比，是基于单变量SSA方法的拓展，其不需要先验信息，能够准确识别数据序列中的周期信号，并将序列分为趋势项、周期项和噪声项三部分，对于处理多变量时间序列具有更明显的优势，被广泛应用于时间序列分析滤波处理和信号提取。
对于不同的研究领域，如何从大量观测时间序列数据中更为准确的滤除噪声并提取有用信号成为诸多领域研究的关键问题，如GRACE时变重力场模型滤波处理。MSSA被成功应用于GRACE时变重力场模型高频噪声及条带误差的滤波，可较为准确的提取有用地球物理信号。通常处理GRACE观测时间序列时，均假定时间序列中的精度是相同的，但是通常观测时间序列的精度是不等的，但对于非完整时间序列而言，需要进行数据插值或多次迭代处理才能进行信号提取，导致滤波处理效果不够好，容易造成信号结构的破坏。
发明内容
本发明的目的就是为了克服上述现有技术存在的缺陷而提供一种基于加权多通道奇异谱的数据序列信号提取方法。
本发明的目的可以通过以下技术方案来实现：
一种基于加权多通道奇异谱的数据序列信号提取方法，该方法的具体步骤如下：
S1：获取多变量数据序列及先验误差序列，对多变量数据序列进行检验，判断其是否存在缺值，并进行中心化处理。
S2：对S1处理后的多变量数据序列进行信号估计。
S3：在S2估计信号的基础上，获取多变量数据序列的噪声项。
S4：以能量保持守恒且不破坏信号结构为基础，根据先验误差序列构建权因子，结合S2得到的估计信号，生成新的数据序列。
S5：对新的数据序列再次进行信号估计，重复步骤S2～S4，提取新的信号。
S6：计算相邻两次提取信号的差值序列，判断其是否满足收敛条件，若满足，输出最后一次估计信号，否则，继续进行循环迭代直到信号结果收敛为止。
对S1处理后的多变量数据序列进行信号估计的具体内容为：
若S1中判断多变量数据序列不存在缺值，则采用传统MSSA方法进行信号估计；若S1中判断多变量数据序列存在缺值，则采用改进的MSSA方法进行信号估计。
采用传统MSSA方法进行信号估计的具体内容为：
211)选定滞后窗口，生成轨迹矩阵，构建滞后协方差阵；
212)对协方差阵进行特征值分解，获取相应特征值和特征向量，并按照特征值大小降序排列；
213)采用W-correlations方法选取d个重构分量，作为前d个主成分；
214)利用确定的前d个主成分重构信号，其余成分归为噪声。
采用改进的MSSA方法进行信号估计的具体内容为：
221)选定滞后窗口，利用窗口内的有效历元构建相应的滞后协方差阵；
222)对滞后协方差阵进行特征值分解，获取相应特征值和特征向量，并按照特征值大小降序排列；
223)基于改进的MSSA方法构建主成分解算方程，采用最小范数约束进行主成分求解；
224)确定重构分量个数d，作为前d个主成分；
225)利用确定的前d个主成分重构信号，将其余的分量归为噪声。
S4中，为保持加权转换前后序列能量不变，构建权因子需满足：
Figure GDA0003632144800000021
式中：σ0l为单位权中误差，σl(t)为先验误差序列，N为有效历元的总数，l为单个通道，通过上式解得σ0l为：
Figure GDA0003632144800000031
则构建的权因子pl(t)的表达式为：
Figure GDA0003632144800000032
生成的新的数据序列x′l(t)的表达式为：
Figure GDA0003632144800000033
式中，sl(t)为步骤2)得到的估计信号，el(t)为步骤3)得到的噪声。
S6中，所述收敛条件为所有历元差值的绝对值小于收敛限差，即|Δsl(t)|＜μ，其中Δsl(t)为所有历元差值，μ为收敛限差。
进一步地，对S1处理后的多变量数据序列进行信号估计时，对于完整数据序列，其迭代过程的具体内容为：
a)对于多变量数据序列xl(t),(t＝1,2,…,N；l＝1,2,…,L)，其中N为有效历元的总数，L为通道数总量，选择滞后窗口M构建大小为的(K×LM)轨迹矩阵X，其中K＝N-M+1为轨迹矩阵列数，轨迹矩阵具体写为X＝[X1；X2；…；XL]T，则单个通道l序列的轨迹矩阵构建如下：
Figure GDA0003632144800000034
b)获取相应滞后协方差阵C＝(1/LM)XXT，其中每个元素计算公式为：
Figure GDA0003632144800000035
式中，xl(i+m-1)和xl(j+m-1)为有效观测值，i、j分别为每个元素在滞后协方差阵中的行数、列数，m为轨迹矩阵X中相应通道l滞后窗口M内的元素列数；
c)对协方差阵C进行特征值分解：
C＝VΛVT
式中，Λ为特征值λk构成的对角阵，λk为C的特征值并按降序排列，vk为正交矩阵V的行向量，按下式计算主成分矩阵A：
A＝VX＝V[X1,X2,…,XL]
相应单个通道l序列的主成分矩阵Al计算如下：
Al＝VXl,(1≤l≤L)
Al的第k行即为第k个主成分，其中每个数值Al(k,i)为：
Figure GDA0003632144800000041
其中v(j,k)是vk的第j个元素，对于单个通道l数据序列第k个重构分量计算如下：
Figure GDA0003632144800000042
d)根据前d个主成分重构获取单个通道l序列的信号sl(t)：
Figure GDA0003632144800000043
对S1处理后的多变量数据序列进行信号估计时，对于缺失数据序列，其迭代过程的具体内容为：
a)对于多变量数据序列xl(t),(t＝1,2,…,N；l＝1,2,…,L)，其中N为有效历元的总数，L为通道数总量，选择滞后窗口M构建(K×LM),K＝N-M+1，构建轨迹矩阵X＝[X1；X2；…；XL]T，则单个通道l序列的轨迹矩阵构建如下：
Figure GDA0003632144800000044
b)获取相应滞后协方差阵C＝(1/LM)XXT，利用改进的MSSA方法的获取滞后协方差阵元素：
Figure GDA0003632144800000045
式中，xl(i+m-1)和xl(j+m-1)均为有效观测值而非缺失值，Nij为有效历元个数，i、j分别为每个元素在滞后协方差阵中的行数、列数，m为轨迹矩阵X中相应通道l滞后窗口M内的元素列数；
c)对协方差阵C进行特征值分解：
C＝VΛVT
式中，Λ为特征值λk构成的对角阵，λk为C的特征值并按降序排列，vk为正交矩阵V的行向量，按下式计算主成分矩阵A：
A＝VX＝V[X1,X2,…,XL]
相应单个通道l序列的主成分矩阵Al计算如下：
Al＝VXl,(1≤l≤L)
Al的第k行即为第k个主成分，对于单个通道l序列，每个数值Al(k,i)为：
Figure GDA0003632144800000051
式中，v(j,k)为vk的第j个元素，Si和
Figure GDA0003632144800000052
分别表示滞后窗口[i,i+K+1]内有效数据集和缺失数据集；
d)利用主成分重构的序列代替上式中的缺值部分，可得：
Figure GDA0003632144800000053
e)将代替后的xl(i+j-1)代入式Al(k,i)，构建解算主成分的方程；
f)依次对每个通道l序列采用W-correlations法进行重构；
g)根据前d个主成分重构获取单个通道l序列的信号sl(t)：
Figure GDA0003632144800000054
进一步地，所述滞后窗口的大小选择为数据序列中包含的周期信号所对应周期的整数倍，且不超过序列总长度的一半。
本发明提供的数据序列信号提取方法，相较于现有技术至少包括如下有益效果：
1)本发明基于先验误差构建权因子，考虑到先验误差主要与噪声项有关，因此只对噪声项加权生成新的时间序列，从而不破坏信号结构；
2)本发明方法由于信号和噪声的初值不确定，因此需要通过迭代处理估计信号直到满足收敛条件；本发明在处理已知先验信息的时间序列数据时，根据先验误差序列构建权因子，仅对噪声加权，在不破坏信号结构的基础上将不等精度序列化为等精度序列，因此相比于传统的MSSA方法在理论上更为合理，加权方式效果更好且可靠性更高。
附图说明
图1为实施例中基于加权多通道奇异谱的数据序列信号提取方法的流程示意图。
图2为实施例中GRACE球谐系数对应的先验误差，其中，子图(a)为2002年4月的先验误差数据，子图(b)为2015年2月的先验误差数据。
图3为实施例中利用Improved MSSA算法得到的前30个主成分的W-correlations结果。
图4为实施例中两种算法提取信号的精度对比及相对提升百分比，其中，子图(a)为Improved MSSA不同阶次球谐系数拟合误差；子图(b)为Weighted MSSA不同阶次球谐系数拟合误差；子图(c)为两种方法拟合误差差异(Improved MSSA-Weighted MSSA)；子图(d)为拟合误差的相对提升百分比。
图5为实施例中两种算法全球质量变化信号对比图，其中，子图(a)、(b)为未做任何滤波情况下的全球质量变化信号，子图(c)、(d)为Improved MSSA提取的全球质量变化信号，子图(e)、(f)为Weighted MSSA提取的全球质量变化信号，子图(a)、(c)、(e)为2002年4月的全球质量变化信号，子图(b)、(d)、(f)为2015年2月的全球质量变化信号。
图6为实施例中两种算法全球陆海信噪比对比图。
图7为实施例中涉及到的完整序列传统MSSA方法流程图。
图8为实施例中涉及到的缺失序列Improved MSSA方法流程图。
具体实施方式
下面结合附图和具体实施例对本发明进行详细说明。显然，所描述的实施例是本发明的一部分实施例，而不是全部实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都应属于本发明保护的范围。
实施例
有学者提出了考虑先验误差的SSA方法并用于提取全球平均海平面(GMSL)变化信号。与SSA类似，也有学者提出了一种考虑先验误差的加权PCA(wPCA)。由于先验误差主要与噪声有关，与信号无关，因此以上两种加权方法，对信号和噪声项均进行加权，在一定程度上会破坏信号结构。因此，考虑到观测的先验误差主要与噪声有关，与信号无关，为避免破坏信号结构，Li&Shen(2018)提出了一种改进的PCA(IPCA)方法，用于GNSS坐标时间序列分析，实验结果证明了IPCA的优越性，验证了考虑先验误差的必要性。考虑到MSSA方法要求时间序列必须是完整的，因此对于完整数据序列，本发明采用传统MSSA方法；对于存在缺失数据的时间序列，采用已有研究提出的能够更好的直接处理不完整多变量时间序列且无需进行任何数据插值或迭代处理的Improved MSSA方法。即基于Li&Shen(2018)和Wang et al.(2020)的基础上，本发明提出了一种基于加权多通道奇异谱的数据序列信号提取方法，该方法顾及先验信息，即基于先验误差构建权因子，仅对噪声项进行加权并通过迭代算法估计信号直到收敛。该算法处理时间序列时由于只对噪声项加权，所以未破坏信号结构，从而能高质量、高可靠性地提取时间序列数据中的信号。因此，本发明在时间序列数据分析和信号处理相关的诸多领域，有潜在的应用推广价值。上述Li&Shen(2018)的具体文献为WeiweiLi、Yunzhong Shen于2018年发表的《The Consideration of Formal Errors inSpatiotemporal Filtering Using Principal Component Analysis for Regional GNSSPosition Time Series》。上述Wang et al.(2020)具体文献为F Wang，Y Shen，T Chen，QChen，W Li于2020年发表的《Improved multi-channel singular spectrum analysis forpost-processing GRACE monthly gravity field models》。
如图1所示，本发明涉及的基于加权多通道奇异谱的数据序列信号提取方法，包括如下步骤内容：
步骤一、获取多变量数据序列xl(t),(t＝1,2,…,N；l＝1,2,…,L)及先验误差序列σl(t)，对多变量数据序列进行检验，判断其是否存在缺值，并进行中心化处理。N为有效历元的总数，L为通道数总量。
步骤二、采用传统MSSA方法进行信号估计。
若步骤一的判断结果为“否”，即数据序列完整，则采用传统MSSA方法进行信号估计；具体地：选定滞后窗口M，生成轨迹矩阵，构建滞后协方差阵；对协方差阵进行特征值分解，得到相应特征值和特征向量，并按照特征值大小降序排列；选取重构分量个数，假设选取的个数为d，采用W-correlations方法进行选择；利用确定的前d个主成分重构信号sl(t)，其余剩下的则归为噪声el(t)。
若步骤一的判断结果为“是”，即数据序列存在缺值，则采用改进的MSSA方法进行信号估计，可采用Wang et al.(2020)提出的Improved MSSA方法进行信号估计；具体地：选定滞后窗口M，利用窗口M内的有效历元构建相应的滞后协方差阵；对滞后协方差阵进行特征值分解，得到相应特征值和特征向量，并按照特征值大小降序排列；基于改进的MSSA方法构建主成分解算方程，存在缺值时方程存在秩亏问题，采用最小范数约束进行求解主成分；确定重构分量个数，假设确定的个数为d，；利用确定的前d个主成分重构信号sl(t)，其余的分量则归为噪声el(t)。
步骤三、将观测序列视为xl(t)＝sl(t)+el(t)，在步骤二估计信号sl(t)的基础上，获取噪声项el(t)＝xl(t)-sl(t)。
步骤四、在能量保持守恒且不破坏信号结构的基础上，构建权因子。为保持加权转换前后序列能量不变，定权需满足，
Figure GDA0003632144800000081
式中，σ0l为单位权中误差，σl(t)为先验误差序列，N为有效历元的总数，l为单个通道。然后可得：
Figure GDA0003632144800000082
设计构建的权因子pl(t)的表达式为：
Figure GDA0003632144800000083
最后生成新的数据序列
Figure GDA0003632144800000084
步骤五、在完成步骤四生成新的数据序列x′l(t)后，需要对新序列再次进行传统或缺失MSSA方法，即重复步骤二～四，提取得到新的信号。
步骤六、计算相邻两次提取信号的差值序列Δsl(t)，判断是否满足收敛条件即所有历元差值绝对值|Δsl(t)|＜μ，其中μ为收敛限差；若判断结果为“是”，则实验终止，输出最后一次估计信号，若判断结果为“否”，则继续进行循环迭代直到信号结果收敛为止。
具体地，迭代过程中涉及的步骤二包括：
对于完整数据序列，采用传统多通道奇异谱分析方法，具体步骤如下：
对于多变量数据序列xl(t),(t＝1,2,…,N；l＝1,2,…,L)，其中，N为有效历元的总数，L为通道数总量。选择滞后窗口M构建大小为(K×LM)的轨迹矩阵X，其中K＝N-M+1为通道l轨迹矩阵列数，轨迹矩阵具体可写为X＝[X1；X2；…；XL]T，则单个通道l序列的轨迹矩阵构建如下：
Figure GDA0003632144800000091
相应滞后协方差阵C＝(1/LM)XXT，其中每个元素计算公式为：
Figure GDA0003632144800000092
式中，式中，xl(i+m-1)和xl(j+m-1)均为有效观测值而非缺失值，i、j对应于每个元素在滞后协方差阵中的行数、列数，m为轨迹矩阵X中相应通道l滞后窗口M内的元素列数。
对协方差阵C进行特征值分解：
C＝VΛVT (3)
式中，Λ是特征值λk构成的对角阵，λk是C的特征值并按降序排列，vk是正交矩阵V的行向量。主成分矩阵A按下式计算：
A＝VX＝V[X1,X2,…,XL] (4)
相应单个通道l序列的主成分矩阵Al计算如下：
Al＝VXl,(1≤l≤L) (5)
Al的第k行即为第k个主成分。其中每个数值Al(k,i)为：
Figure GDA0003632144800000093
其中v(j,k)是vk的第j个元素。对于单通道l数据序列第k个重构分量计算如下：
Figure GDA0003632144800000094
由于式(3)中的特征值是按降序排列的，因此时间序列可由前几个主成分分量重构，其余则视为噪声，即单个通道l序列的信号sl(t)可以由前d个主成分重构：
Figure GDA0003632144800000095
对于缺值数据序列，则利用缺值多通道奇异谱分析方法，具体步骤如下：对于缺值数据序列，Wang et al.(2020)针对性提出了缺值MSSA(Improved MSSA)方法，滞后协方差阵元素计算公式为：
Figure GDA0003632144800000101
其中xl(i+m-1)和xl(j+m-1)均为有效观测值而非缺失值，Nij表示有效历元个数。对于单个通道l序列，式(6)可重新表示为：
Figure GDA0003632144800000102
其中Si和
Figure GDA0003632144800000103
分别表示窗口[i,i+K+1]内有效数据集和缺失数据集。利用主成分重构的序列代替式(10)中的缺值部分：
Figure GDA0003632144800000104
将式(11)代入式(10)，构建解算主成分的方程。由于方程是秩亏的因此引入最小范数准则求解主成分。主成分计算完成以后，后续过程则采用与传统MSSA算法同样过程依次对每个通道l序列进行重构。对于实施过程中滞后窗口M和重构分量个数d选取问题，为了提高序列中不同尺度信号的可区分性，窗口大小通常优选为数据序列中含有的周期信号的周期的整数倍，且不超过序列总长度的一半，重构阶数则可依据W-correlations法进行确定。
具体地，迭代过程中涉及的步骤二～步骤四主要包括：
对于多通道时间序列xl(t),(t＝1,2,…,N；l＝1,2,…,L)通常可写为：
xl(t)＝sl(t)+el(t) (12)
其中sl(t)为信号，el(t)视为噪声。鉴于序列xl(t)是不等精度序列，因此根据先验误差σl(t)构建权重因子来处理式(12)的数据序列，相当于直接处理生成的具有等权的数据序列x′l(t)，考虑到先验误差σl(t)主要与噪声项有关，因此仅对噪声项el(t)进行加权，从而不破坏信号结构。得到的数据序列x′l(t)的表达式为：
Figure GDA0003632144800000105
其中pl(t)是在保证数据序列变换前后能量不变的基础上提出的，即所有有效历元权重和应等于有效历元的总数N。
Figure GDA0003632144800000106
由式(14)可推导出相应单位权中误差为：
Figure GDA0003632144800000111
从而可得设计构建权因子为：
Figure GDA0003632144800000112
通常情况下信号sl(t)和噪声el(t)都是不确定的，因此需进行迭代计算，信号和噪声的迭代初值可利用传统MSSA算法(步骤二中的内容)进行估计，更新信号sl(t)和噪声el(t)，依据式(13)生成新的数据序列，需要对新序列x′l(t)再次进行传统或缺值多通道奇异谱分析即重复步骤二～四，进行迭代实验直到提取信号序列收敛为止。
本实施例以文献Wang et al.(2020)中的GRACE时变重力场模型数据算例作为实施算例，由于GRACE数据本身存在缺失，为了更好的比较本发明从GRACE数据序列提取信号的效果，实施例中将利用缺失MSSA方法(Improved MSSA)与本发明加权MSSA(WeightedMSSA)算法进行对比。采用GRACE CSR-RL06球谐数据序列；时间跨度为2002.04～2016.08；时间分辨率为月，截断阶次为60。以2002.04与2005.02两个月为例展示了GRACE的先验误差具体见图2，可以发现GRACE阶数越高误差越大，阶数越低精度越高，且不同阶次的先验误差是不同的。
首先进行步骤S1：对GRACE数据序列进行检验判断是否存在缺值并进行中心化处理，首先实施Improved MSSA对含噪不完整GRACE球谐数据序列进行处理。
步骤S2：利用Improved MSSA算法对含噪序列进行信号提取，通过实验对比优选窗口长度M＝60；前10个主成分(PC)能量占比为80.19％，其中PCs 1+4主要为趋势信号，能量占比约为29.42％，PCs 2+3代表周年，PCs 8-10代表半年，PCs5-7代表其他长周期信号。根据W-corrections确定重构分量个数取10(见图3)。对含噪序列进行估计信号sl(t)作为迭代实验的信号和噪声初值。
步骤S3：在步骤S2估计信号sl(t)的基础上，获取噪声项el(t)＝xl(t)-sl(t)；
步骤S4：在能量保持守恒且不破坏信号结构的基础上，基于先验误差序列构建权因子
Figure GDA0003632144800000113
生成新的数据序列
Figure GDA0003632144800000114
步骤S5：在完成步骤S4生成新的数据序列后，需要对新序列x′l(t)再次进行缺值多通道奇异谱分析(Improved MSSA)即重复步骤S2～S4，提取得到新的信号。
步骤S6：计算相邻两次提取信号的差值序列Δsl(t)，判断是否满足收敛条件即所有历元差值绝对值
Figure GDA0003632144800000121
若判断结果为“是”，则实验终止，输出最后一次估计信号，若判断结果为“否”，则继续进行循环迭代直到信号结果收敛为止。为验证两种算法提取信号的效果，采用单位权中误差进行评价，公式如下：
Figure GDA0003632144800000122
Figure GDA0003632144800000123
其中xl(t)为原始序列，
Figure GDA0003632144800000124
为重构信号(Improved MSSA)，
Figure GDA0003632144800000125
为重构信号(Weighted MSSA)，pl(t)为权因子。图4给出了两种算法提取GRACE球谐系数信号的精度统计结果，可以看出考虑先验误差的Weighted MSSA方法的所有阶次的单位权中误差均小于Improved MSSA，表明其提取信号更多；相对提升百分比结果表明球谐系数随着阶数的增高，提取信号效果提升越明显，该结果与GRACE时变重力场阶数越高噪声越大，考虑先验误差的Weighted MSSA信号提取效果就越好相一致。
为进一步验证两种算法提取信号在空间质量变化上的差异，将两种算法提取的球谐系数按照时变重力场反演公式计算出相应的质量变化，具体计算公式采用Wahr et al.,(1998)中的计算公式，Wahr et al.,(1998)的具体文献为Wahr J，Molenaar M，Bryan F于1998发表的《Time variability of the Earth's gravity field:Hydrological andoceanic effects and their possible detection using GRACE》。图5以2002.04和2015.02为例展示两种算法的信号提取对比效果。
图5结果表明两种算法滤波效果均较为明显，通过与Improved MSSA对比发现，考虑先验误差的Weighted MSSA的结果中条带噪声更少，而前面结果显示其提取信号更多，表明该算法能够提取更多的地球物理信号且滤波效果更好。为进一步确认Weighted MSSA滤波效果的优越性，利用文献Chen et al.(2006)提出的全球陆海信噪比RMS_Ratio指标对两种算法滤波效果再次进行评估，计算公式如下：
Figure GDA0003632144800000126
其中MASSland和MASSocean分别表示陆地和海洋上的信号，Err为相应噪声。上述Chenet al.(2006)具体文献为Chen J.L.,Wilson C.R.,K.-W.Seo.于2006年发表的《Optimizedsmoothing of gravity recovery and climate experiment(GRACE)time-variablegravity observations》。具体计算信噪比结果见图6，结果显示考虑先验误差的WeightedMSSA的所有月份信噪比均高于Improved MSSA，充分证明了Weighted MSSA在信号提取滤除噪声方面的优越性，也表明滤波过程中考虑先验误差的必要性。
以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的工作人员在本发明揭露的技术范围内，可轻易想到各种等效的修改或替换，这些修改或替换都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以权利要求的保护范围为准。

Claims (7)

1.一种基于加权多通道奇异谱的数据序列信号提取方法，其特征在于，包括下列步骤：
1)获取多变量数据序列及先验误差序列，对多变量数据序列进行检验，判断其是否存在缺值，并进行中心化处理；
2)对步骤1)处理后的多变量数据序列进行信号估计；若步骤1)中判断多变量数据序列不存在缺值，则采用传统MSSA方法进行信号估计，采用传统MSSA方法进行信号估计的具体内容为：
211)选定滞后窗口，生成轨迹矩阵，构建滞后协方差阵；
212)对协方差阵进行特征值分解，获取相应特征值和特征向量，并按照特征值大小降序排列；
213)采用W-correlations方法选取d个重构分量，作为前d个主成分；
214)利用确定的前d个主成分重构信号，其余成分归为噪声；
若步骤1)中判断多变量数据序列存在缺值，则采用改进的MSSA方法进行信号估计；采用改进的MSSA方法进行信号估计的具体内容为：
221)选定滞后窗口，利用窗口内的有效历元构建相应的滞后协方差阵；
222)对滞后协方差阵进行特征值分解，获取相应特征值和特征向量，并按照特征值大小降序排列；
223)基于改进的MSSA方法构建主成分解算方程，采用最小范数约束进行主成分求解；
224)确定重构分量个数d，作为前d个主成分；
225)利用确定的前d个主成分重构信号，将其余的分量归为噪声；
3)在步骤2)估计信号的基础上，获取多变量数据序列的噪声项；
4)以能量保持守恒且不破坏信号结构为基础，根据先验误差序列构建权因子，结合步骤2)得到的估计信号，生成新的数据序列；
5)对新的数据序列再次进行信号估计，重复步骤2)～4)，提取新的信号；
6)计算相邻两次提取信号的差值序列，判断其是否满足收敛条件，若满足，输出最后一次估计信号，否则，继续进行循环迭代直到信号结果收敛为止。
2.根据权利要求1所述的基于加权多通道奇异谱的数据序列信号提取方法，其特征在于，步骤4)中，为保持加权转换前后序列能量不变，构建权因子需满足：
Figure FDA0003632144790000021
式中：σ0l为单位权中误差，σl(t)为先验误差序列，N为有效历元的总数，l为单个通道，通过上式解得σ0l为：
Figure FDA0003632144790000022
则构建的权因子pl(t)的表达式为：
Figure FDA0003632144790000023
3.根据权利要求2所述的基于加权多通道奇异谱的数据序列信号提取方法，其特征在于，步骤4)中，生成的新的数据序列x′l(t)的表达式为：
Figure FDA0003632144790000024
式中，sl(t)为步骤2)得到的估计信号，el(t)为步骤3)得到的噪声。
4.根据权利要求1所述的基于加权多通道奇异谱的数据序列信号提取方法，其特征在于，步骤6)中，所述收敛条件为所有历元差值的绝对值小于收敛限差，即|Δsl(t)|＜μ，其中Δsl(t)为所有历元差值，μ为收敛限差。
5.根据权利要求1所述的基于加权多通道奇异谱的数据序列信号提取方法，其特征在于，对步骤1)处理后的多变量数据序列进行信号估计时，对于完整数据序列，其迭代过程的具体内容为：
a)对于多变量数据序列xl(t),(t＝1,2,…,N；l＝1,2,…,L)，其中N为有效历元的总数，L为通道数总量，选择滞后窗口M构建大小为的(K×LM)轨迹矩阵X，其中K＝N-M+1为轨迹矩阵列数，轨迹矩阵具体写为X＝[X1；X2；…；XL]T，则单个通道l序列的轨迹矩阵构建如下：
Figure FDA0003632144790000025
b)获取相应滞后协方差阵C＝(1/LM)XXT，其中每个元素计算公式为：
Figure FDA0003632144790000026
式中，xl(i+m-1)和xl(j+m-1)为有效观测值，i、j分别为每个元素在滞后协方差阵中的行数、列数，m为轨迹矩阵X中相应通道l滞后窗口M内的元素列数；
c)对协方差阵C进行特征值分解：
C＝VΛVT
式中，Λ为特征值λk构成的对角阵，λk为C的特征值并按降序排列，vk为正交矩阵V的行向量，按下式计算主成分矩阵A：
A＝VX＝V[X1,X2,…,XL]
相应单个通道l序列的主成分矩阵Al计算如下：
Al＝VXl,(1≤l≤L)
Al的第k行即为第k个主成分，其中每个数值Al(k,i)为：
Figure FDA0003632144790000031
其中v(j,k)是vk的第j个元素，对于单个通道l数据序列第k个重构分量计算如下：
Figure FDA0003632144790000032
d)根据前d个主成分重构获取单个通道l序列的信号sl(t)：
Figure FDA0003632144790000033
6.根据权利要求1所述的基于加权多通道奇异谱的数据序列信号提取方法，其特征在于，对步骤1)处理后的多变量数据序列进行信号估计时，对于缺失数据序列，其迭代过程的具体内容为：
a)对于多变量数据序列xl(t),(t＝1,2,…,N；l＝1,2,…,L)，其中N为有效历元的总数，L为通道数总量，选择滞后窗口M构建(K×LM),K＝N-M+1，构建轨迹矩阵X＝[X1；X2；…；XL]T，则单个通道l序列的轨迹矩阵构建如下：
Figure FDA0003632144790000034
b)获取相应滞后协方差阵C＝(1/LM)XXT，利用改进的MSSA方法的获取滞后协方差阵元素：
Figure FDA0003632144790000041
式中，xl(i+m-1)和xl(j+m-1)均为有效观测值而非缺失值，Nij为有效历元个数，i、j分别为每个元素在滞后协方差阵中的行数、列数，m为轨迹矩阵X中相应通道l滞后窗口M内的元素列数；
c)对协方差阵C进行特征值分解：
C＝VΛVT
式中，Λ为特征值λk构成的对角阵，λk为C的特征值并按降序排列，vk为正交矩阵V的行向量，按下式计算主成分矩阵A：
A＝VX＝V[X1,X2,…,XL]
相应单个通道l序列的主成分矩阵Al计算如下：
Al＝VXl,(1≤l≤L)
Al的第k行即为第k个主成分，对于单个通道l序列，每个数值Al(k,i)为：
Figure FDA0003632144790000042
式中，v(j,k)为vk的第j个元素，Si和
Figure FDA0003632144790000043
分别表示滞后窗口[i,i+K+1]内有效数据集和缺失数据集；
d)利用主成分重构的序列代替上式中的缺值部分，可得：
Figure FDA0003632144790000044
e)将代替后的xl(i+j-1)代入式Al(k,i)，构建解算主成分的方程；
f)依次对每个通道l序列采用W-correlations法进行重构；
g)根据前d个主成分重构获取单个通道l序列的信号sl(t)：
Figure FDA0003632144790000045
7.根据权利要求6所述的基于加权多通道奇异谱的数据序列信号提取方法，其特征在于，所述滞后窗口的长度选择为数据序列中包含的周期信号所对应周期的整数倍，且不超过序列总长度的一半。
CN202110212835.3A 2021-02-25 2021-02-25 一种基于加权多通道奇异谱的数据序列信号提取方法 Active CN113032713B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110212835.3A CN113032713B (zh) 2021-02-25 2021-02-25 一种基于加权多通道奇异谱的数据序列信号提取方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110212835.3A CN113032713B (zh) 2021-02-25 2021-02-25 一种基于加权多通道奇异谱的数据序列信号提取方法

Publications (2)

Publication Number Publication Date
CN113032713A CN113032713A (zh) 2021-06-25
CN113032713B true CN113032713B (zh) 2022-07-05

Family

ID=76462125

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110212835.3A Active CN113032713B (zh) 2021-02-25 2021-02-25 一种基于加权多通道奇异谱的数据序列信号提取方法

Country Status (1)

Country Link
CN (1) CN113032713B (zh)

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107449964A (zh) * 2017-07-20 2017-12-08 淮阴工学院 一种用于模式重建和预测的广义多元奇异谱分析方法
CN110618985A (zh) * 2019-08-28 2019-12-27 同济大学 一种基于改进奇异谱的包含乘性噪声数据序列处理方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106226407B (zh) * 2016-07-25 2018-12-28 中国电子科技集团公司第二十八研究所 一种基于奇异谱分析的超声回波信号在线预处理方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107449964A (zh) * 2017-07-20 2017-12-08 淮阴工学院 一种用于模式重建和预测的广义多元奇异谱分析方法
CN110618985A (zh) * 2019-08-28 2019-12-27 同济大学 一种基于改进奇异谱的包含乘性噪声数据序列处理方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
利用奇异谱分析探测GPS坐标时间序列时变周期信号;刘磊;《江西测绘》;20200325(第01期);全文 *
基于MSSA的区域GPS站点季节性信号提取;汪浩等;《大地测量与地球动力学》;20190515(第05期);全文 *

Also Published As

Publication number Publication date
CN113032713A (zh) 2021-06-25

Similar Documents

Publication Publication Date Title
Fessler Optimization methods for magnetic resonance image reconstruction: Key models and optimization algorithms
CN110361778B (zh) 一种基于生成对抗网络的地震数据重建方法
Elsner et al. Singular spectrum analysis: a new tool in time series analysis
CN109740453B (zh) 一种基于小波变换的卫星磁场数据地震前兆异常提取方法
US20170108604A1 (en) Denoising seismic data
CN112069919A (zh) 基于非凸低秩矩阵逼近和全变分正则化的高光谱图像去噪方法
CN108303740B (zh) 一种航空电磁数据噪声压制方法及装置
CN112084847B (zh) 基于多通道截断核范数及全变差正则化的高光谱图像去噪方法
Francis et al. Nonlinear prediction of the ionospheric parameter ƒ o F 2 on hourly, daily, and monthly timescales
CN110244353B (zh) 一种基于稀疏范数优化算法的地震数据规则化方法
CN114091538B (zh) 一种基于信号特征的判别损失卷积神经网络智能降噪方法
Wang et al. Data and model dual-driven seismic deconvolution via error-constrained joint sparse representation
Kumar et al. Enabling uncertainty quantification for seismic data preprocessing using normalizing flows (NF)—An interpolation example
CN114372239A (zh) 一种去除钻孔应变数据环境影响因素的方法
CN113421198B (zh) 一种基于子空间的非局部低秩张量分解的高光谱图像去噪方法
CN113032713B (zh) 一种基于加权多通道奇异谱的数据序列信号提取方法
CN107944474B (zh) 基于局部自适应字典的多尺度协作表达高光谱分类方法
CN114428278A (zh) 基于压缩感知地震数据重建方法、装置、电子设备及介质
CN113158830A (zh) 一种剩余重力异常场分离方法
CN110618985B (zh) 一种基于改进奇异谱的包含乘性噪声数据序列处理方法
CN113009564B (zh) 地震数据处理方法和装置
CN113640891B (zh) 一种基于奇异谱分析的瞬变电磁探测数据噪声滤除方法
CN115034978A (zh) 基于子空间表示结合图嵌入约束的高光谱图像去噪方法
CN111856559B (zh) 基于稀疏贝叶斯学习理论的多道地震谱反演方法及系统
CN108387928B (zh) 基于地震特征变换空间的数据构建方法

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
GR01 Patent grant
GR01 Patent grant