一种坐标时间序列周期项的提取方法与流程

1.本发明涉及一种坐标时间序列周期项的提取方法，属于定位导航技术领域。


背景技术：

2.坐标时间序列中包含显著的非线性变化(尤其是垂直方向上)，研究表明，非线性变化 的振幅可以达到1～2cm。然而，基准站的坐标时间序列作为地球参考框架建立和维持的关键 基础数据，这种非线性变化会对地球参考框架的精度和稳定性产生影响，因此研究坐标时间 序列的非线性变化是进一步提升参考框架精度和稳定性的重要途径。
3.坐标时间序列中的非线性变化主要包括：1.由地球物理效应引起的基准站坐标真实的非 线性变化，比如环境负载(包括大气负载、水文负载、非潮汐海洋负载等)和热膨胀效应等； 2.由数据处理模型不完善或其他系统误差引起的虚假的非线性变化，比如天线相位中心误差 和光压模型误差等；3.其他未模型化的因素引起的形变以及观测噪声。从非线性变化的来源 来看，其主要表现形式是各种周期性变化，即不同的周期项。因此，研究坐标非线性变化无 法避开的一个问题就是坐标时间序列中周期项的探测识别和提取。
4.目前，坐标时间序列中周期项的提取方法主要有功率谱分析、小波分析、奇异值谱分析 (singular spectrum analysis,ssa)等。其中，功率谱分析是最直观、也是目前使用最多的 方法。采用功率谱分析进行周期项的探测识别的主要思路是通过求取坐标时间序列的功率谱 或功率谱密度，进而在频域对潜在的周期项进行探测。
5.传统的功率谱估计方法主要是基于傅里叶变换，包括：相关函数法、周期图法等，其中 相关函数法以时间序列的相关函数为媒介，通过求取相关函数的傅里叶变换的得到相应的功 率谱；周期图法则直接对时间序列做傅里叶变换求得功率谱。后续，有人提出通过对已有时 间序列进行建模来获取相应的功率谱，比较常用的是最大熵谱估计法。
6.无论是传统的功率谱估计方法还是后续提出的现代功率谱估计方法，在得到坐标时间序 列的功率谱后，只能通过人工或者设置阈值的方式对潜在的周期项进行探测识别，而坐标时 间序列中往往包含着噪声，噪声会对潜在的周期信号产生影响，甚至掩盖较小的周期信号， 进而导致潜在的周期信号无法被识别或错将噪声识别为周期信号，导致周期项提取的准确性 较低。


技术实现要素：

7.本申请的目的在于提供一种坐标时间序列周期项的提取方法，用以解决现有周期项提取 准确性低的问题。
8.为实现上述目的，本申请提出了一种坐标时间序列周期项的提取方法的技术方案，包括 以下步骤：
9.1)获取原始坐标时间序列，计算原始坐标时间序列的功率谱，称为原始功率谱；所述 原始功率谱包括频点、以及频点对应的原始功率谱密度值；
10.2)将原始坐标时间序列进行r次随机重排，计算每次随机重排后的坐标时间序列
的功 率谱，称为重排功率谱；所述重排功率谱中包括频点、以及频点对应的重排功率谱密度值； r大于设定次数阈值；
11.3)根据每次重排后，每个频点的原始功率谱密度值、以及重排功率谱密度值计算每个 频点对应的置信水平；所述置信水平越大，该频点为噪声的概率越低；置信水平越小，该频 点为噪声的概率越大；
12.4)在原始功率谱中，挑选出置信水平大于置信阈值的频点，得到去噪后的功率谱；
13.5)在去噪后的功率谱中进行周期项的提取。
14.本发明中坐标时间序列周期项的提取方法的技术方案的有益效果是：本发明通过对原始 坐标时间序列和随机重排后的坐标时间序列求取功率谱，进而得到每个频点的置信水平，通 过置信水平判断频点是否为噪声点：置信水平大，频点为噪声的概率低；置信水平小，频点 为噪声的概率大，因此在原始功率谱中，将置信水平较小的频点进行剔除，得到去噪的功率 谱，进而在去噪的功率谱中提取周期项。本发明通过随机重排计算每个频点的置信水平，去 除了噪声对周期项提取的影响，提高了周期性提取的科学性和准确性。
15.进一步的，置信水平的计算过程如下：
[0016][0017][0018]
其中，cl(f
k
)为频率f
k
的置信水平；r为随机重排的次数；p
x
(f
k
)为频率f
k
对应的 原始功率谱密度值；为第i次重排后，频率f
k
对应的重排功率谱密度值。
[0019]
进一步的，为了提高周期项提取的准确性，所述步骤5)中提取周期项的方法为：在去 噪后的功率谱中确定若干个主瓣区域；找出每个主瓣区域中最大原始功率谱密度值对应的频 点，该频点的倒数即为该主瓣区域的周期项。
[0020]
进一步的，为了提高主瓣区域确定的准确性，利用分段原则将去噪后的功率谱划分为多 个主频点辐射区域，每个主频点辐射区域包括一个主瓣区域及若干旁瓣区域；所述分段原则 为：相邻候选频点的频率差值大于2ofac
×
t
s
，或候选频点的区间峰值出现先减小后增大的情 形；所述候选频点为挑选出的置信水平大于置信阈值对应的频点；所述候选频点的区间峰值 为主瓣或者旁瓣的峰值；ofac为过采样率，t
s
为原始坐标时间序列的采样间隔。
[0021]
进一步的，为了提高周期项提取的准确性，通过第一过采样率计算原始功率谱，通过第 二过采样率计算最大原始功率谱密度值，用于确定最终的周期项对应的频点，第二过采样率 大于第一过采样率。
[0022]
进一步的，通过lomb
–
scargle周期图法计算功率谱。
附图说明
[0023]
图1是本发明坐标时间序列周期项的提取方法的流程图；
[0024]
图2是本发明去噪后功率谱根据第一个分段原则进行主频点辐射区域划分的示意图；
[0025]
图3是本发明去噪后功率谱根据第二个分段原则进行主频点辐射区域划分的示意图。
具体实施方式
[0026]
坐标时间序列周期项的提取方法实施例：
[0027]
本发明的主要构思在于，针对坐标时间序列中噪声的影响，导致周期项提取准确性低的 问题，本发明首先得到原始坐标时间序列的功率谱、以及原始坐标时间序列经过r次随机 重排后的功率谱，接着通过原始功率谱和重排功率谱得到每个频点下原始功率和重排功率的 大小关系，进而得到每个频点的置信水平，通过置信水平的阈值判断将原始坐标时间序列中 噪声的影响去除或减弱，提高周期性提取的准确性。
[0028]
具体的，坐标时间序列周期项的提取方法如图1所示，包括以下步骤：
[0029]
1)获取原始坐标时间序列，原始坐标时间序列x(t
j
)包括n个采样点，其中t
j
为第j个 采样时刻，j＝1,2,
…
,n。
[0030]
2)通过lomb
–
scargle周期图法求取原始坐标时间序列的功率谱，称为原始功率谱，原 始功率谱包括频点、以及频点对应的原始功率谱密度值。
[0031]
lomb
–
scargle周期图法是由scargle在lomb和barning工作的基础上进一步完善得到 的一种功率谱估计方法，对于原始坐标时间序列x(t
j
)，其正则化的原始功率谱密度值为：
[0032][0033]
其中，ω为角频率，根据公式ω＝2πf可将角频率ω转换至频率f；p
x
(ω)为角频率ω 对应的原始功率谱密度值；σ2为原始坐标时间序列的方差；τ为相移因子，其计算公式为：
[0034][0035]
本实施例中，为了进一步提高周期性提取的准确性和科学性，采用lomb
–
scargle周期 图法进行功率谱的计算，当然，本发明对功率谱计算的方法并不做限制。lomb
–
scargle周期 图法相对于其他功率谱估计方法的最大优势在于其能够处理非均匀采样的时间序列数据，而 gnss坐标时间序列中通常包含着一定数量的缺失值，在选用其他方法进行周期项提取时必 须先进行缺失值的内插，内插算法的误差通常也会对周期项的探测产生影响。
[0036]
3)对步骤1)中的原始坐标时间序列x(t
j
)进行r次(r＞设定次数阈值，设定次数阈 值为1000或者更高，对此不做限制)的随机重排，第i次重排后得到随机重排后的坐标时 间序列并通过lomb
–
scargle周期图法计算每次随机重排后的功率谱，称为重排功率 谱，重排功率谱中包括频点、以及频点对应的重排功率谱密度值。
[0037]
4)通过步骤2)的原始功率谱和步骤3)的重排功率谱得到每个频点对应的置信水
平。
[0038]
置信水平的计算过程如下：
[0039][0040][0041]
其中，cl(f
k
)为频率f
k
的置信水平，能够充分反映原始功率谱中信号的真实性； p
x
(f
k
)为频率f
k
对应的原始功率谱密度值；为第i次重排后，频率f
k
对应的重排功 率谱密度值。
[0042]
5)在原始功率谱中，挑选出置信水平cl(f
k
)大于置信阈值s的频点f
k
，也即挑选出噪 声概率小的频点，得到去噪后的功率谱；所挑选出的频点f
k
，构成一个候选频点集合f
candi
， f
candi
＝[f1,f2,
…
,f
n
]，f1＜f2＜
…
＜f
n
。
[0043]
频点集合f
candi
中理论上包含了所有真实的周期信号对应的频点，得到去除噪声后的如 图2所示的频点
‑
原始功率谱密度值的功率谱。
[0044]
6)在去噪后的功率谱中进行周期项的提取。
[0045]
本实施例中，利用分段原则将去噪后的功率谱划分为多个主频点辐射区域，如图2和图3所示，每个主频点辐射区域包括一个主瓣区域及若干旁瓣区域，确定主瓣区域后，找出每 个主瓣区域中最大原始功率谱密度值对应的频点，该频点的倒数即为该主瓣区域的周期项。
[0046]
分段原则的目的是为了将真实的频点提取出来，也即将主瓣区域提取出来，分段原则的 分段原理如下：
[0047]
lomb
–
scargle周期图法中频率的采样间隔δf为：
[0048][0049]
其中，ofac为过采样率，t
s
为原始坐标时间序列的采样间隔，而根据傅里叶变换的原理： 求解有限长坐标时间序列的功率谱，相当于给无限长的时间序列加了一个跨度为l的矩形窗， 因此，在时间序列的功率谱中会出现以主瓣为中心的对称分布并且逐级递减的旁瓣。由矩形 窗的性质可知，主瓣的宽度为2/l,旁瓣的宽度1/l，因此主瓣峰值与相邻旁瓣峰值以及相邻 的两个旁瓣峰值之间的间隔通常小于2/l，即小于2ofac
×
t
s
。
[0050]
而且由一个主频点辐射区域(主瓣和旁瓣形成的区域)过渡到另一个主频点辐射区域必 定会出现功率谱密度的区间峰值先减小后增大的情况。因此，分段原则为：
[0051]
①
相邻候选频点的频率差值大于2ofac
×
t
s
。这里的候选频点为挑选出的置信水平大于置 信阈值的频点，即去噪后的功率谱中的频点。原始功率谱进行去噪时，可能会出现连续去除 多个频点的现象，导致去噪后的功率谱中出现大片空白的现象，因此，这里将候选频点中相 邻候选频点的频率差值大于2ofac
×
t
s
进行直接划分。如图2所示：相邻候选频点f
m
和f
m 1
之 间的频率差值大于2ofac
×
t
s
，则候选频点f
m
属于主频点辐射区域1，而候选频点f
m 1
属于主 频点辐射区域2。
[0052]
②
如图3所示，候选频点的区间峰值出现先减小后增大的情形，也即功率谱密度的
区间 峰值出现先减小后增大的情形，候选频点的区间峰值为主瓣或者旁瓣的峰值。
[0053]
如果满足以上两个分段原则中的一个，则表明在这两个候选频点属于不同的主频点辐射 区域，进而在这两个频点之间进行主频点辐射区域的划分。
[0054]
主频点辐射区域划分之后，每个主频点辐射区域中，最大原始功率谱密度值对应的频点 即为最终的周期项对应的频点，通过t＝1/f的转换得到最终的周期项。
[0055]
在实际处理过程中，为了提升计算速度，原始功率谱和重排功率谱的频率采样间隔可适 当采用较小的过采样率，为第一过采样率，但为了进一步避免频率采样间隔过大造成的峰值 位置的偏差，通常会以更高的过采样率(即第二过采样率，第二过采样率大于第一过采样率) 在主频点辐射区域内最大原始功率谱密度值对应的频点附近的设定频率范围内重新计算原始 坐标时间序列的功率谱密度值，以获得更加准确的周期项。当然，作为其他实施方式，在对 周期项精度要求不高的情况下，也可以只通过第一过采样率确定周期。
[0056]
本实施例中，周期对应的频率的不确定度(也即周期对应的频率的误差)表示为：
[0057][0058]
其中，为去除周期分量后噪声的方差，σ
n
为去除周期分量后噪声的标准差；n
t
为原 始坐标时间序列的总长度，a为对应周期信号的振幅。
[0059]
进而得到周期d的不确定度表示为：
[0060][0061]
上述实施例中，为了提高周期项提取的准确性，通过分段原则进行主频点辐射区域的划 分，进而确定真实频点，得到周期项，作为其他实施方式，在去除噪声之后，也可以通过人 眼识别去噪后的功率谱，直接读取各个主频点辐射区域的最大原始功率谱密度值对应的频点 确定周期项，本发明对此不做限制。
[0062]
本发明通过随机重排能够使原始坐标时间序列变成杂乱无序的噪声，并且每一次重排后 得到的相同频点上的功率谱密度值也在发生变化，通过大量的重排统计得到的置信水平能够 充分反映原始坐标时间序列功率谱中信号的真实性，进而对噪声信号进行剔除，减弱或避免 噪声的影响，大大提高周期项提取的科学性和准确性。

技术特征：
1.一种坐标时间序列周期项的提取方法，其特征在于，包括以下步骤：1)获取原始坐标时间序列，计算原始坐标时间序列的功率谱，称为原始功率谱；所述原始功率谱包括频点、以及频点对应的原始功率谱密度值；2)将原始坐标时间序列进行r次随机重排，计算每次随机重排后的坐标时间序列的功率谱，称为重排功率谱；所述重排功率谱中包括频点、以及频点对应的重排功率谱密度值；r大于设定次数阈值；3)根据每次重排后，每个频点的原始功率谱密度值、以及重排功率谱密度值计算每个频点对应的置信水平；所述置信水平越大，该频点为噪声的概率越低；置信水平越小，该频点为噪声的概率越大；4)在原始功率谱中，挑选出置信水平大于置信阈值的频点，得到去噪后的功率谱；5)在去噪后的功率谱中进行周期项的提取。2.根据权利要求1所述的坐标时间序列周期项的提取方法，其特征在于，置信水平的计算过程如下：算过程如下：其中，cl(f
k
)为频率f
k
的置信水平；r为随机重排的次数；p
x
(f
k
)为频率f
k
对应的原始功率谱密度值；为第i次重排后，频率f
k
对应的重排功率谱密度值。3.根据权利要求1所述的坐标时间序列周期项的提取方法，其特征在于，所述步骤5)中提取周期项的方法为：在去噪后的功率谱中确定若干个主瓣区域；找出每个主瓣区域中最大原始功率谱密度值对应的频点，该频点的倒数即为该主瓣区域的周期项。4.根据权利要求3所述的坐标时间序列周期项的提取方法，其特征在于，利用分段原则将去噪后的功率谱划分为多个主频点辐射区域，每个主频点辐射区域包括一个主瓣区域及若干旁瓣区域；所述分段原则为：相邻候选频点的频率差值大于2ofac
×
t
s
，或候选频点的区间峰值出现先减小后增大的情形；所述候选频点为挑选出的置信水平大于置信阈值的频点；所述候选频点的区间峰值为主瓣或者旁瓣的峰值；ofac为过采样率，t
s
为原始坐标时间序列的采样间隔。5.根据权利要求3所述的坐标时间序列周期项的提取方法，其特征在于，通过第一过采样率计算原始功率谱，通过第二过采样率计算最大原始功率谱密度值，用于确定最终的周期项对应的频点，第二过采样率大于第一过采样率。6.根据权利要求1所述的坐标时间序列周期项的提取方法，其特征在于，通过lomb
–
scargle周期图法计算功率谱。
技术总结
本发明涉及一种坐标时间序列周期项的提取方法，属于定位导航技术领域。方法包括：获取原始坐标时间序列，计算原始坐标时间序列的功率谱；将原始坐标时间序列进行R次随机重排，计算每次随机重排后的功率谱；根据每次重排后，每个频点的原始功率谱密度值、以及重排功率谱密度值确定每个频点的置信水平；置信水平越大，该频点为噪声的概率越低；置信水平越小，该频点为噪声的概率越大；在原始功率谱中，挑选出置信水平大于置信阈值的频点，得到去噪后的功率谱；在去噪后的功率谱中进行周期项的提取。本发明通过随机重排计算每个频点的置信水平，去除了噪声对周期项提取的影响，提高了周期性提取的科学性和准确性。期性提取的科学性和准确性。期性提取的科学性和准确性。


技术研发人员：贾彦锋 朱新慧 孙付平 肖凯 何劢航 柯能
受保护的技术使用者：中国人民解放军战略支援部队信息工程大学
技术研发日：2021.03.18
技术公布日：2021/6/29