一种从GNSS时间序列中提取周日半周日海潮信号的方法与流程
本发明涉及高精度全球卫星定位系统定位及其应用领域，在利用高精度全球卫星定位系统提取海洋潮汐负荷形变领域中应用，具体地说是一种从gnss时间序列中提取周日半周日海潮信号的方法。
背景技术：
在月球和太阳引力的共同作用下固体地球表面会产生周期性位移及海洋质量周期性的迁移，称为固体地球潮和海潮负荷。这些潮汐会改变固体地球的旋转和平移状态，它被地球自转参数和地心参数吸收，也导致全球卫星定位系统持续观测台站在时间序列中出现周期性的位移变化。固体潮会导致低纬度区域的垂直位移超过几厘米，海潮会导致沿海地区台站产生厘米级的垂直位移。
通常固体潮和海潮导致的台站周期性位移会在全球卫星定位系统数据处理过程中利用模型模拟掉。其中固体潮和海潮长周期项(mf，mm和ssa)导致的台站位移模型拟合精度较高，约为1％左右，但是周日半周日的高频海潮信号模型拟合效果差，在单日解的时间序列中仍有较大残留，学者们提出各种方法利用空间观测技术来提取周日半周日海潮信号以便提高海潮模型精度。目前利用全球卫星定位系统来提取高频海潮信号的方法主要有以下三种：静态谐波估计方法(静态方法)，动态估计法和假频(alias)谐波分析法。
静态谐波估计法把主要的海洋潮汐负荷形变(通常8个)作为未知参数在求解单日解台站坐标的过程中解算，单日内解算这些海潮参数，其法方程矩阵肯定是奇异的，需要加入约束消除奇异性，并用多天解及协方差矩阵联合求解。动态估计法先求解高频的台站坐标，海洋潮汐负荷的信息包含在台站解中(解的频率必须高于周日和半日信号的尼奎斯nyquist频率)。然后用最小二乘拟合台站坐标时间序列中潮汐信号的振幅和相位。假频谐波分析法是指利用传统的谱分析方法求解单日解时间序列中没有全部平均掉的周日半周日海潮时间序列，这些高频的海潮信号会以低频信号的方式出现在单日解时间序列中。如果周日半周日假频海潮信号能被解出，可以由它转换算出原先的周日半周日潮汐信号。
静态和动态估计法都是直接用原始潮汐波频率解算，避免间接计算的模糊度和转换误差，解的稳定度和精度都较好。此外，动态估计法可用较短的序列，如不到1年序列。但是这两种方法都需要重新对所有的资料重新处理求解，工作量浩大。其中静态估计法还需要改动软件，求解过程中有秩亏，必须加约束，约束带来畸变，导致计算效率低，收敛慢；动态高频坐标解和大气天顶距延迟参数高度相关，如资料较短，季节性负荷影响很大。假频谐波分析法虽然可以直接使用单日解时间序列，但是主要用谱分析方法来检测海洋潮汐波残差的存在，对于周期较长的假频信号(例如s2,k1,p1和k2的假频信号周期为周年或者半周年)，如果没有长的时间序列则难以直接求解；此外，假频信息较弱，不唯一，容易和其他影响因素混淆，很难和其他负荷(如大气)的季节项形变区分开来。
技术实现要素：
本发明的目的是为了能简单直接使用各大数据分析中心提供的单日解时间序列提取周日半周日海潮信号，弥补假频谐波分析方法的不足而提出了另外一种提取方法，利用混频技术来从单日解时间序列中提取周日半周日海潮信号。该方法起源于电子通讯领域，然后引入到空间大地测量的vlbi数据处理领域，用于高频(周日和半周日)地球自转研究领域。
实现本发明目的的具体技术方案是：
一种从gnss时间序列中提取周日半周日海潮信号的方法，该方法包括以下具体步骤：步骤1：gnss单日解时间序列预处理
首先选择数据利用率至少90％的台站单日解原始时间序列进行预处理，探测中断，剔除粗差；然后从时间序列中去除常数项、线性项、同震跳跃、非地震跳跃、震后变形和局部多项式项，保留周年及半周年项；最后，用线性插值方法对时间序列分析结果进行插值，获得新的均匀间隔时间序列；
步骤2：周日半周日海潮信号混频
选取固定频率信号乘以原有时间序列，把原有时间序列转化成高频的频率和与低频的频率差两个信号之和；
公式(1)中a(t)为原始高频的时间序列，其中a是幅度，ω是频率，t即采样间隔，
公式(1)乘以一个固定频率信号cosω0t后变为新的、振幅相等和频率分量ω+ω0和差频分量ω-ω0的和，即公式(2)；
公式(2)中a′(t)为混频后新的时间序列，其中a是幅度，ω是频率，ω0是固定信号频率，t即采样间隔，
步骤3：周日半周日海潮信号滤波
选择低通滤波器将时间序列中含有周日半周日海潮信号的高频信号ω+ω0与低频信号ω-ω0分离；
步骤4：周日半周日海潮信号提取
用公式(3)所示的周期图谱分析法中welch谱分析法处理混频后含有周日半周日海潮信号的低频信号ω-ω0；
其中，jp(ω)表示分段数为p的周期函数，xp(n)表示长度为n的信号分成p段，每段m个数据，
步骤5：周日半周日海潮信号滤波振幅和相位拟合
使用最小二乘在检测到的差频信号的谱峰处分别拟合它们的振幅和相位；振幅拟合根据公式(5)计算振幅压缩效应；
其中，a'j为考虑压幅效应后差频信号j的振幅，fj是周期信号的频率，aj表示来利用最小二乘拟合的振幅，δt为数据处理时间间隔，π为圆周率。
本发明提出了一种从gnss时间序列中提取周日半周日海潮信号的方法，其优势有以下两点：一、通过利用固定频率信号乘以原有时间序列，使得原有时间序列变成为高频(频率和)和低频(频率差)两个信号的和，满足假频谐波分析法存在的，信号不唯一及容易受到其他信号干扰等弱点。二、可以直接利用各数据处理中心提供的单日解时间序列直接求解周日半周日海潮信号，数据处理简单，计算量小，效率低，便于使用，对于新手和非专业人士而言，学习和使用简便，有利于广泛推广。
附图说明
图1为本发明流程图；
图2为本发明的实施框图；
图3为本发明方法实施例chur台站单日解时间序列；
图4为利用假频谐波分析法提取chur台站坐标单日解时间序列中周日半周日海潮信号结果图；
图5为本发明方法提取chur台站坐标单日解时间序列中包含周日海潮信号结果图；
图6为本发明方法提取chur台站坐标单日解时间序列中包含半周日海潮信号结果图。
具体实施方式
以下结合附图1及附图2对本发明进行详细描述。
本发明包括以下具体步骤：
步骤1：gnss单日解时间序列预处理
首先需要对各大数据处理中心提供的单日解进行时间序列进行预处理。在台站选择时应该尽量选择数据可利用率至少90％的台站。数据来源该时间序列包含了中断、粗差以及观测网共同误差等，因此需要在时间序列分过程中首先进行预处理，探测中断，剔除粗差。然后需要从时间序列中去除常数项，线性项，同震跳跃，非地震跳跃，震后变形和局部多项式项，保留可能包含周日半周日假频海潮信号的周年及半周年项。由于微小的数据缺失也会影响谱分析结果，因此为了使影响最小化，还需要使用线性插值方法对时间序列处理结果进行插值，以获得新的均匀间隔时间序列。
步骤2：周日半周日海潮信号混频
混频就是选择一个固定频率信号乘以原有时间序列，使得原有时间序列转化成高频的频率和以及低频的频率差两个信号之和，使得后期处理可以更容易在高频信号或者低频信号上进行。
公式(1)中a(t)为原始高频的时间序列，其中a是幅度，ω是频率，t即采样间隔，
公式(1)乘以一个固定频率信号cosω0t后变为新的、振幅相等和频率分量ω+ω0和差频分量ω-ω0的和，即公式(2)；
公式(2)中a′(t)为混频后新的时间序列，其中a是幅度，ω是频率，ω0是固定信号频率，t即采样间隔，
步骤3：周日半周日海潮信号滤波
与原始时间序列相比，新序列的和频分量ω+ω0变为较高频信号，而差频分量ω-ω0变为低频信号。两个分量的振幅是原始时间序列振幅的一半。低频信号ω-ω0很容易利用低通滤波器将其和其他信号分离。然而，该方法存在一个基本问题，即镜像频率。例如，想要通过1mhz差频来检测大约20mhz的原始频率的时间序列，则有两个固定频率信号可产生1mhz差频，即为21mhz和19mhz，21mhz频率为需要的固定频率，而19mhz则为混频产生的镜像频率，可以通过选择固定信号频率的方法，降低镜像频率的影响。
步骤4：周日半周日的海潮信号提取
用公式(3)所示的周期图谱分析法中welch谱分析法处理混频后含有周日半周日海潮信号的低频信号ω-ω0；
其中，jp(ω)表示分段数为p的周期函数，xp(n)表示长度为n的信号分成p段，每段m个数据，
为了能很好的识别假频周期为1年左右周日半周日海潮信号，必须根据时间序列的长度选用合理的分段数，以便使每段的有效长度超过2.5年。
步骤5：周日半周日海潮信号振幅、相位拟合
welch谱分析法能够明确的检测出时间序列中谱峰对应的信号及其频率(差频)，进一步转换为原始频率及对应的周日半周日海潮信号。然后可以使用最小二乘在检测到的谱峰处(差频)分别拟合它们的振幅和相位。此外，由于全球卫星定位系统持续观测台站单日解是台站坐标在24小时时段的平均解，需要考虑了振幅压缩效应。周日半周日海潮信号振幅的压幅效应根据公式(6)来计算。
其中，a'j为考虑压幅效应后差频信号j的振幅，fj是周期信号的频率，aj表示来利用最小二乘拟合的振幅，δt为数据处理时间间隔，π为圆周率。
实施例
为验证本方法，以美国scrippsorbitandpermanentarraycenter(sopac)提供持续站chur的单日解时间序列为例，本实施例具体包括以下步骤：
1.gnss单日解时间序列预处理
首先需要利用时间序列分析法处理各大数据处理中心提供的单日解时间序列。在台站选择时应该尽量选择数据可利用率大于90％的台站。sopac提供持续站的单日解时间序列由sopac采用“st_filter”软件，将分别由美国nasajetpropulsionlaboratory(jpl)(采用gipsy软件)和sopac(采用gamit/globk软件)分析中心得到的松约束解进行组合，以消除由软件及处理策略导致的误差,最终得到一个统一、自洽的解。sopac提供原始时间序列(raw)、去粗差后的时间序列(clean)以及去除共模误差后的时间序列(filter)三种类型的数据。时间序列包含了中断、粗差以及观测网共同误差等，因此需要在时间序列分过程中首先进行预处理，探测中断，剔除粗差。然后需要从时间序列中去除常数项，线性项，同震跳跃，非地震跳跃，震后变形和局部多项式，保留可能包含周日半周日假频海潮信号的周年及半周年项。由于微小的数据缺失也会影响谱分析结果，因此为了使影响最小化，还需要使用线性插值方法对时间序列分析结果进行插值，以获得新的均匀间隔时间序列，实施例所选chur台站坐标时间序列见图3。
2.周日半周日海潮信号混频
混频就是选择一个固定频率信号乘以原有时间序列，使得原有时间序列转化成高频的频率和以及低频的频率差两个信号之和，使得后期处理可以更容易在高频或者低频上进行。
选取固定频率信号乘以原有时间序列，把原有时间序列转化成高频的(频率和)和与低频的(频率差)两个信号之和；
公式(1)中a(t)为原始高频的时间序列，其中a是幅度，ω是频率，t即采样间隔，
公式(1)乘以一个固定频率信号cosω0t后变为新的、振幅相等和频率分量ω+ω0和差频分量ω-ω0的和，即公式(2)；
公式(2)中a′(t)为混频后新的时间序列，其中a是幅度，ω是频率，ω0是固定信号频率，t即采样间隔，
3.周日半周日海潮信号滤波
与原始时间序列相比，新序列的和频分量ω+ω0变为较高频率信号，而差频分量ω-ω0变为低频信号。两个分量的振幅是原始时间序列振幅的一半。低频信号ω-ω0很容易利用低通滤波器将其和其他信号分离。然而，该方法存在一个基本问题，即镜像频率。例如，想要通过1mhz差频来检测大约20mhz的原始频率的时间序列，则有两个固定频率信号可产生1mhz差频，即为21mhz和19mhz，21mhz频率为需要的固定频率，而19mhz则为混频产生的镜像频率，可以通过选择固定信号频率的方法，降低镜像频率的影响。
4.周日半周日海潮信号提取
用公式(3)所示的周期图谱分析法中welch谱分析法处理混频后含有周日半周日海潮信号的低频信号ω-ω0；
其中，jp(ω)表示分段数为p的周期函数，xp(n)表示长度为n的信号分成p段，每段m个数据，
为了能很好的识别假频周期为1年左右周日半周日海潮信号，必须根据时间序列的长度选用合理的分段数，以便使每段的有效长度超过2.5年。
图4为利用假频谐波分析法提取chur台站坐标单日解时间序列中周日半周日海潮信号结果,从图中可以明显看出，图5为本发明方法提取chur台站坐标单日解时间序列中包含周日海潮信号结果，图6为本发明方法提取chur台站坐标单日解时间序列中包含半周日海潮信号结果。对比发现，本发明方法可以很好的检测时间序列中残留的p1、k1、s2及k2四个潮波，而假频谐波分析法则只能探测p1、k1，而且不能很好的区分开这个个潮波，k2也在假频谐波分析结果中有体现，前提是利用该台长长达23年的时间序列。根据我们的数据处理经验，如果时间序列较短，如六、七年以内，要用假频方法检测到k2是十分困难的。
5.周日半周日海潮信号振幅、相位合拟
welch谱分析法能够明确的检测出时间序列中谱峰对应的信号及其频率(差频)，进一步转换为原始频率及对应的周日半周日海潮信号。然后可以使用最小二乘在检测到的谱峰处(差频)分别拟合它们的振幅和相位。此外，由于cgps单日解是台站坐标在24小时时段的平均解，需要考虑了振幅压缩效应。周日半周日海潮信号振幅的压幅效应根据公式(6)来计算，chur拟合结果见表1。
上式中a'j为考虑压幅效应后差频信号j的振幅，fj是周期信号的频率，aj表示来利用最小二乘拟合的振幅，δt表示数据处理时间间隔，π为圆周率。
表1实施例chur振幅、相位拟合结果
表1中罗列了从台站chur的单日解时间序列中提取的周日海潮信号p1及k1，以及半周日海潮信号s2及k2，验证了混频技术从gnss单日解时间序列中提取周日半周日海潮信号的可行性。
- 还没有人留言评论。精彩留言会获得点赞！