CN104765055A - Gps测站坐标时间序列周期性探测方法及系统 - Google Patents

Abstract

本发明提供了一种GPS测站坐标时间序列周期性探测方法及系统，包括步骤：步骤1，获取GPS测站坐标时间序列观测值；步骤2，剔除GPS测站坐标时间序列观测值粗差并修正天线相位中心偏差；步骤3，采用谐函数描述GPS测站坐标时间序列获得GPS测站坐标时间序列的谐函数模型，并设计谐函数模型矩阵；步骤4，采用最小二乘法解算谐函数模型矩阵，探测GPS测站坐标时间序列周期性。本发明可能处理含有时变周期性信号的GPS测站坐标时间序列，并保证探测准确性。

Description

GPS测站坐标时间序列周期性探测方法及系统

技术领域

本发明属于GPS数据精密处理技术领域，具体涉及一种GPS测站坐标时间序列周期性探测方法及系统。

背景技术

GPS测站坐标时间序列主要用于噪声模型和速度模型的建立。现有的GPS测站坐标时间序列预测估计模型通常包含季节性信号，并且将周期性信号如周年信号和半周年等作为常数参数广泛应用(Blewitt and Lavallée,2002)，一般采用周年周期及其谐波和常数振幅的正弦曲线建模(Nicolaidis,2002)。

GPS测站坐标时间序列中存在的季节性信号是对环境变化的响应，而隐藏在季节性信号中的环境噪声是不稳定的。研究成果表明，环境噪声的表现形式具有周期性可变特征(Davis,Wernicke et al.,2012)。

然而，数据处理分析策略与GPS相关的技术因素、固体地球质量负荷、大气压负荷、水储量变化、岩石热膨胀、降雨、多路径效应等因素的综合影响，导致GPS测站坐标时间序列中的季节性信号呈现周期性可变特征(Davis,Wernicke etal.,2012)。田云锋(2011)对青藏高原和喜马拉雅地区的GPS垂向分量分析发现存在明显的相位变化，推断是受到地表水体因子控制；而在CMONOC及周边IGS台站的位置时间序列中均发现了周期约为351/n(n＝1,…,6)天的“异常”周期项，并认为地表质量负荷并非该异常周期项的来源。

目前GPS测站坐标时间序列周期性探测方法存在不足：将周期性时变信号当作时不变信号进行描述，对速率估值及其不确定性都会产生明显影响。

发明内容

针对现有技术存在的问题，本发明提供了一种利用谐函数和最小二乘准则、顾及周期性可变信号的GPS测站坐标时间序列周期性探测方法及系统。

本发明采用谐函数对GPS测站坐标时间序列建模，通过最小二乘准则估计GPS测站周期性可变信号，从而实现GPS测站坐标时间序列周期性探测。

为解决上述技术问题，本发明采用如下的技术方案：

一、GPS测站坐标时间序列周期性探测方法，包括步骤：

步骤1，获取GPS测站坐标时间序列观测值；

步骤2，剔除GPS测站坐标时间序列观测值粗差并修正天线相位中心偏差；

步骤3，采用谐函数描述GPS测站坐标时间序列获得GPS测站坐标时间序列的谐函数模型，并设计谐函数模型矩阵其中，y表示GPS测站坐标时间序列观测值； $A=\left[\begin{array}{cc}1& t_1\\ 1& t_2\\ ...& ...\\ 1& t_m\end{array}\right],$ $A_k=\left[\begin{array}{cc}\cos {\mathrm{\ω}}_k{t}_1& \sin {\mathrm{\ω}}_k{t}_1\\ \cos {\mathrm{\ω}}_k{t}_2& \sin {\mathrm{\ω}}_k{t}_2\\ .& .\\ .& .\\ .& .\\ \cos {\mathrm{\ω}}_k{t}_m& \sin {\mathrm{\ω}}_k{t}_m\end{array}\right],$ t_i为观测历元，i＝1,2,...,m，m为观测历元数量，ω_k表示可变周期性信号的频率； $x_k=\left[\begin{array}{c}{a}_k\\ b_k\end{array}\right],$ a_k、b_k分别为正弦函数分量、余弦函数分量； $x=\left[\begin{array}{c}y\left(0\right)\\ a_0\end{array}\right],$ y(0)为GPS测站参考历元坐标，a₀为GPS测站的线性速度；k表示谐函数编号，q为谐函数数量；

步骤4，采用最小二乘法解算谐函数模型矩阵，获得可变周期性信号的频率估计值，并采用假设检验方法验证频率估计值；

步骤5，基于假设检验验证后的频率估计值，采用最小二乘方法解算谐函数模型矩阵获得正弦函数分量a_k、余弦函数分量b_k和线性速度a₀。

步骤1中通过双差定位软件工具、精密单点定位软件工具或IGS分析中心获取GPS测站坐标时间序列观测值。

步骤3中，所述的GPS测站坐标时间序列的谐函数模型为：

$y\left(t_i\right)=y\left(0\right)+a_0{t}_i+\underset{k=1}{\overset{q}{\mathrm{\Σ}}}[a_k\cos \left({\mathrm{\ω}}_k{t}_i\right)+b_k\sin \left({\mathrm{\ω}}_k{t}_i\right)]$

其中，y(t_i)为观测历元t_i对应的GPS测站坐标观测值，y(0)为GPS测站参考历元坐标；a₀为GPS测站的线性速度，t_i为观测历元；a_k、b_k分别为正弦函数分量、余弦函数分量，ω_k表示可变周期性信号的频率；k表示谐波函数编号，q为谐波函数数量。

步骤4进一步包括子步骤：

4.1解算以下最大化问题以探测频率ω_j及其对应的A_j：

${\mathrm{\ω}}_j=\arg \underset{{\mathrm{\ω}}_l}{\max }{\left|\right|P_{{\stackrel{\‾}{A}}_l}y\left|\right|}_{Q_y^{-1}}^2,{\stackrel{\‾}{A}}_l=P_{\stackrel{\‾}{A}}^{\⊥}{A}_l\∞$

$P_{\stackrel{\‾}{A}}^{\⊥}=I-\stackrel{\‾}{A}{\left({\stackrel{\‾}{A}}^T{Q}_y^{-1}\stackrel{\‾}{A}\right)}^{-1}{\stackrel{\‾}{A}}^T{Q}_y^{-1}$

$P_{{\stackrel{\‾}{A}}_l}={\stackrel{\‾}{A}}_l{\left({\stackrel{\‾}{A}}_l^T{Q}_y^{-1}{\stackrel{\‾}{A}}_l\right)}^{-1}{\stackrel{\‾}{A}}_l^T{Q}_y^{-1}$

其中，I为单位矩阵，是L2范数运算符，Q_y是GPS测站坐标时间序列观测值的方差协方差矩阵； $\stackrel{\‾}{A}=[A,A_1,...,A_{j-1}],$ $A_l=\left[\begin{array}{cc}\cos {\mathrm{\ω}}_l{t}_1& \sin {\mathrm{\ω}}_l{t}_1\\ \cos {\mathrm{\ω}}_l{t}_2& \sin {\mathrm{\ω}}_l{t}_2\\ .& .\\ .& .\\ .& .\\ \cos {\mathrm{\ω}}_l{t}_m& \sin {\mathrm{\ω}}_l{t}_m\end{array}\right],$ l＝1,2,...j-1；

4.2采用假设 $\left\{\begin{array}{c}{H}_0:y=\mathrm{Ax}+\underset{k=1}{\overset{j-1}{\mathrm{\Σ}}}\left(A_k{x}_k\right)\\ H_a:y=\mathrm{Ax}+\underset{k=1}{\overset{j}{\mathrm{\Σ}}}\left(A_k{x}_k\right)\end{array}\right.$ 判断是否向谐函数模型增加新频率，即根据步骤4.1获得的探测频率ω_j及其对应的A_j判定是否接受零假设，若接受，则增加新频率，然后执行步骤4.1；否则，执行步骤4.3；H₀为零假设，H_a为备选假设；

4.3采用Q_y＝σ²I和统计量评估假设检验，其中，为零假设下的最小二乘残差；为备选假设下的验后方差，为备选假设下的最小二乘残差，统计量T₂≈F(2,m-n-2s)，s为假设检验过程中确定的备选谐函数数量。

步骤5具体为：

采用最小二乘法，通过解算以下矩阵获得正弦函数分量a_k和余弦函数分量b_k：

y＝Bx+v

$\hat{x}={\left(B^T{Q}_y^{-1}B\right)}^{-1}{B}^T{Q}_y^{-1}y$

$Q_{\hat{x}}={\left(B^T{Q}_Y^{-1}B\right)}^{-1}$

其中，矩阵B＝[A A_k]，k＝1,...,j；为x的最小二乘估值；Q_y是GPS测站坐标时间序列观测值的方差协方差矩阵；y是GPS测站坐标时间序列观测值，v是随机误差矢量，采用最小二乘法解算获得。

二、GPS测站坐标时间序列周期性探测系统，包括：

观测值获得模块，用来获取GPS测站坐标时间序列观测值；

修正模块，用来剔除GPS测站坐标时间序列观测值粗差并修正天线相位中心偏差；

谐函数模型构建模块，用来采用谐函数描述GPS测站坐标时间序列获得GPS测站坐标时间序列的谐函数模型，并设计谐函数模型矩阵其中，y表示GPS测站坐标时间序列观测值； $A=\left[\begin{array}{cc}1& t_1\\ 1& t_2\\ ...& ...\\ 1& t_m\end{array}\right],$ $A_k=\left[\begin{array}{cc}\cos {\mathrm{\ω}}_k{t}_1& \sin {\mathrm{\ω}}_k{t}_1\\ \cos {\mathrm{\ω}}_k{t}_2& \sin {\mathrm{\ω}}_k{t}_2\\ .& .\\ .& .\\ .& .\\ \cos {\mathrm{\ω}}_k{t}_m& \sin {\mathrm{\ω}}_k{t}_m\end{array}\right],$ t_i为观测历元，i＝1,2,...,m，m为观测历元数量，ω_k表示可变周期性信号的频率； $x_k=\left[\begin{array}{c}{a}_k\\ b_k\end{array}\right],$ a_k、b_k分别为正弦函数分量、余弦函数分量； $x=\left[\begin{array}{c}y\left(0\right)\\ a_0\end{array}\right],$ y(0)为GPS测站参考历元坐标，a₀为GPS测站的线性速度；k表示谐函数编号，q为谐函数数量；

频率估计模块，用来采用最小二乘法解算谐函数模型矩阵，获得可变周期性信号的频率估计值，并采用假设检验方法验证频率估计值；

周期性探测模块，用来基于假设检验验证后的频率估计值，采用最小二乘方法解算谐函数模型矩阵获得正弦函数分量a_k、余弦函数分量b_k和线性速度a₀。

与现有技术相比，本发明具有特点：

本发明采用谐函数建立GPS测站坐标时间序列预估模型，突破了现有技术中仅采用周期性时不变信号描述坐标时间序列的现状；建立顾及周期性可变特征信号模型，一方面有助于进一步精化已有的GPS测站坐标时间序列的速度模型，提高速度估值的可靠性；另一方面，顾及周期性可变特征信号模型有助于真实反映与当地环境效应相关的(如降雨、温度、地表负荷及含水层抽水等)季节性地球物理信号，以进一步开展地球物理过程等方面的建模及解释。

现有GPS测站坐标时间序列周期性探测中，将周期性信号当作时不变信号进行估计，这样导致GPS测站坐标时间序列速度模型的不完善、速率估值及其不确定性的较大偏差，本发明可解决上述问题，可处理含有时变周期性信号的GPS测站坐标时间序列，并保证周期性探测的准确性。

附图说明

图1为本发明方法的具体流程示意图。

具体实施方式

为了使本发明目的、技术方案及有益效果更加清楚明白，下面将结合附图及具体实施方式，进一步说明本发明。应当理解，以下描述的具体实施方式仅用以解释本发明，并不用于限定本发明。

本发明具体步骤如下：

步骤1，获取GPS测站坐标时间序列观测值。

本步骤属于本技术领域的常规技术，具体可通过双差定位软件(如成熟的GAMIT/GLOBK数据处理软件)、精密单点定位软件(如成熟的GIPSY/OASIS数据处理软件)、或IGS分析中心获取GPS测站坐标时间序列观测值。

步骤2，剔除GPS测站坐标时间序列观测值中粗差，修正与更换硬件相关的天线相位中心偏差。

本步骤属于本技术领域内常规的数据预处理方式，在此不作赘述。

步骤3，采用谐函数描述GPS测站坐标时间序列，获得GPS测站坐标时间序列谐函数模型。

GPS测站坐标时间序列谐函数模型如下：

$y\left(t_i\right)=y\left(0\right)+a_0{t}_i+\underset{k=1}{\overset{q}{\mathrm{\Σ}}}[a_k\cos \left({\mathrm{\ω}}_k{t}_i\right)+b_k\sin \left({\mathrm{\ω}}_k{t}_i\right)]---\left(1\right)$

式(1)中：

y(t_i)为观测历元t_i对应的GPS测站坐标观测值；

y(0)为GPS测站参考历元坐标；

a₀为GPS测站的线性速度；

系数a_k、b_k分别为正弦函数分量、余弦函数分量，用来描述GPS测站周期性运动的振幅，对应于频率ω_k和观测历元t_i，其中，i＝1,2,...,m，m为观测历元数量；

q为谐函数数量，k为谐函数编号，这里的谐函数指a_kcos(ω_kt_i)+b_ksin(ω_kt_i)。

采用矩阵表达公式(1)，获得GPS测站坐标时间序列谐函数模型矩阵：

$y=\mathrm{Ax}+\underset{k=1}{\overset{q}{\mathrm{\Σ}}}\left(A_k{x}_k\right)---\left(2\right)$

式(2)中：

y表示GPS测站坐标时间序列观测值；

A为GPS测站坐标时间序列谐函数模型的线性回归系数构成的2列矩阵，即 $A=\left[\begin{array}{cc}1& t_1\\ 1& t_2\\ ...& ...\\ 1& t_m\end{array}\right];$

A_k为cos(ω_kt)和sin(ω_kt)构成的2列矩阵，即 $A_k=\left[\begin{array}{cc}\cos {\mathrm{\ω}}_k{t}_1& \sin {\mathrm{\ω}}_k{t}_1\\ \cos {\mathrm{\ω}}_k{t}_2& \sin {\mathrm{\ω}}_k{t}_2\\ .& .\\ .& .\\ .& .\\ \cos {\mathrm{\ω}}_k{t}_m& \sin {\mathrm{\ω}}_k{t}_m\end{array}\right];$

$x_k=\left[\begin{array}{c}{a}_k\\ b_k\end{array}\right];$

$x=\left[\begin{array}{c}y\left(0\right)\\ a_0\end{array}\right],$ y(0)为GPS测站参考历元坐标。

GPS测站坐标时间序列谐函数模型中待估参数包括频率ω_k及系数a_k、b_k。

本具体实施方式中，计算GPS测站坐标时间序列时不仅考虑测站的线性速度和周年/半周年等周期性信号影响，更考虑了可变周期性信号。

与以往GPS测站坐标时间序列建模方法不同，本发明充分考虑了不同因素对GPS测站坐标时间序列造成的可变周期性，同时利用利用三角函数系数及其对应的频率将这些可变周期性纳入计算体系，并通过假设检验方法对待估频率ω_k(即可变周期性信号的频率)进行数值计算、验证。

步骤4，解算GPS测站坐标时间序列谐函数模型矩阵获得待估频率ω_k，并采用统计检验方法验证待估频率ω_k。

采用最小二乘法确定待估频率ω_k，具体如下：

$\left\{\begin{array}{c}{H}_0:y=\mathrm{Ax}+\underset{k=1}{\overset{j-1}{\mathrm{\Σ}}}\left(A_k{x}_k\right)\\ H_a:y=\mathrm{Ax}+\underset{k=1}{\overset{j}{\mathrm{\Σ}}}\left(A_k{x}_k\right)\end{array}\right.---\left(3\right)$

式(3)中：

H₀表示零假设，H_a表示备选假设，j＝1、2、…、q。

该过程迭代完成，j从1到q取值；每次运行，如果接受零假设，即向函数模型增加一个新的频率。

4.1解算以下最大化问题以探测待估频率ω_j(Teunissen,2000)及其对应的A_j：

${\mathrm{\ω}}_j=\arg \underset{{\mathrm{\ω}}_l}{\max }{\left|\right|P_{{\stackrel{\‾}{A}}_l}y\left|\right|}_{Q_y^{-1}}^2,{\stackrel{\‾}{A}}_l=P_{\stackrel{\‾}{A}}^{\⊥}{A}_l---\left(4\right)$

$P_{\stackrel{\‾}{A}}^{\⊥}=I-\stackrel{\‾}{A}{\left({\stackrel{\‾}{A}}^T{Q}_y^{-1}\stackrel{\‾}{A}\right)}^{-1}{\stackrel{\‾}{A}}^T{Q}_y^{-1}---\left(5\right)$

$P_{{\stackrel{\‾}{A}}_l}={\stackrel{\‾}{A}}_l{\left({\stackrel{\‾}{A}}_l^T{Q}_y^{-1}{\stackrel{\‾}{A}}_l\right)}^{-1}{\stackrel{\‾}{A}}_l^T{Q}_y^{-1}---\left(6\right)$

公式(4)～(6)中：

I为单位矩阵；

运算符表示L2范数；

Q_y是GPS测站坐标时间序列观测值y的方差协方差矩阵；

A_l为的子矩阵，l＝1,2,...j-1，A_l结构同A_k相同，即 $A_l=\left[\begin{array}{cc}\cos {\mathrm{\ω}}_l{t}_1& \sin {\mathrm{\ω}}_l{t}_1\\ \cos {\mathrm{\ω}}_l{t}_2& \sin {\mathrm{\ω}}_l{t}_2\\ .& .\\ .& .\\ .& .\\ \cos {\mathrm{\ω}}_l{t}_m& \sin {\mathrm{\ω}}_l{t}_m\end{array}\right],$ 用于假设检验待估频率ω_l对应的系数a_l和b_l。满足公式(4)-(6)中取最大值时的ω_l即待估频率ω_j，ω_l对应的A_l被设置为系数矩阵A_j。

采用公式(4)～(6)所示的最大化问题解算待估频率ω_j，采用公式(3)的假设检验判定是否增加新的待估频率。

采用数值方法解算上述问题，GPS测站坐标时间序列的功率谱通过应用不同频率ω_j的谱值具有最大谱值的连续频率用来构建矩阵A_j。

4.2采用Q_y＝σ²I评估假设检验，先验单位权方差σ²未知。采用下述统计量进行评估：

${\stackrel{\‾}{A}}_i=P_{\stackrel{\‾}{A}}^{\⊥}{A}_i---\left(8\right)$

式(7)～(8)中：

为零假设下的最小二乘残差；

为备选假设下的验后方差，采用计算，其中，df为自由度，为备选假设下的最小二乘残差；

公式(7)中统计量T₂的分布为自由度为2和m-n-2s的F-分布，即T₂≈F(2,m-n-2s)，m为观测历元数量，n为未知参数数量，本发明中n为3，m-n表示模型冗余度；s为假设检验过程中确定的备选谐函数数量。

本步骤为假设检验统计标准过程，为便于理解，特说明本步骤过程：T₂为统计量，根据统计量T₂对应的T₂≈F(2,m-n-2s)分布可以得到一个具体分布值；另根据公式(7)也可计算T₂的分布值，将两个分布值比较，判断假设检验中是接受零假设、还是接受备选假设。

步骤5，基于假设检验验证后的频率ω_k，采用最小二乘方法解算系数a_k、b_k。

采用最小二乘法，解算以下矩阵：

y＝Bx+v        (9)

$\hat{x}={\left(B^T{Q}_y^{-1}B\right)}^{-1}{B}^T{Q}_y^{-1}y---\left(10\right)$

$Q_{\hat{x}}={\left(B^T{Q}_Y^{-1}B\right)}^{-1}---\left(11\right)$

式(9)～(11)中：

设计矩阵B包括矩阵A和A_k,k＝1,...,j，即B＝[A A_k]；

为x的最小二乘估值；

v是随机误差矢量，采用最小二乘法解算获得。

采用上述方法估计任一GPS测站可变周期性信号频率及其对应振幅、线性速度等参数，从而实现更完善的速度模型的建立，提高速度估值可靠性。

本文中所描述的具体实施仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施做各种的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (6)

1.GPS测站坐标时间序列周期性探测方法，其特征是，包括步骤：

步骤1，获取GPS测站坐标时间序列观测值；

步骤2，剔除GPS测站坐标时间序列观测值粗差并修正天线相位中心偏差；

步骤3，采用谐函数描述GPS测站坐标时间序列获得GPS测站坐标时间序列的谐函数模型，并设计谐函数模型矩阵其中，y表示GPS测站坐标时间序列观测值； $A=\left[\begin{array}{cc}1& t_1\\ 1& t_2\\ ...& ...\\ 1& t_m\end{array}\right],$ $A_k=\left[\begin{array}{cc}\cos {\mathrm{\ω}}_k{t}_1& \sin {\mathrm{\ω}}_k{t}_1\\ \cos {\mathrm{\ω}}_k{t}_2& \sin {\mathrm{\ω}}_k{t}_2\\ .& .\\ .& .\\ .& .\\ \cos {\mathrm{\ω}}_k{t}_m& \sin {\mathrm{\ω}}_k{t}_m\end{array}\right],$ t_i为观测历元，i＝1,2,...,m，m为观测历元数量，ω_k表示可变周期性信号的频率； $x_k=\left[\begin{array}{c}{a}_k\\ b_k\end{array}\right],$ a_k、b_k分别为正弦函数分量、余弦函数分量； $x=\left[\begin{array}{c}y\left(0\right)\\ a_0\end{array}\right],$ y(0)为GPS测站参考历元坐标，a₀为GPS测站的线性速度；k表示谐函数编号，q为谐函数数量；

步骤4，采用最小二乘法解算谐函数模型矩阵，获得可变周期性信号的频率估计值，并采用假设检验方法验证频率估计值；

步骤5，基于假设检验验证后的频率估计值，采用最小二乘方法解算谐函数模型矩阵获得正弦函数分量a_k、余弦函数分量b_k和线性速度a₀。

2.如权利要求1所述的GPS测站坐标时间序列周期性探测方法，其特征是：

步骤1中通过双差定位软件工具、精密单点定位软件工具或IGS分析中心获取GPS测站坐标时间序列观测值。

3.如权利要求1所述的GPS测站坐标时间序列周期性探测方法，其特征是：

所述的GPS测站坐标时间序列的谐函数模型为：

$y\left(t_i\right)=y\left(0\right)+a_0{t}_i+\underset{k=1}{\overset{q}{\mathrm{\Σ}}}[a_k\cos \left({\mathrm{\ω}}_k{t}_i\right)+b_k\sin \left({\mathrm{\ω}}_k{t}_i\right)]$

其中，y(t_i)为观测历元t_i对应的GPS测站坐标观测值，y(0)为GPS测站参考历元坐标；a₀为GPS测站的线性速度，t_i为观测历元；a_k、b_k分别为正弦函数分量、余弦函数分量，ω_k表示可变周期性信号的频率；k表示谐波函数编号，q为谐波函数数量。

4.如权利要求1所述的GPS测站坐标时间序列周期性探测方法，其特征是：

步骤4进一步包括子步骤：

4.1解算以下最大化问题以探测频率ω_j及其对应的A_j：

${\mathrm{\ω}}_j=\arg \underset{{\mathrm{\ω}}_l}{\max }{\left|\right|P_{{\stackrel{\‾}{A}}_l}y\left|\right|}_{Q_y^{-1}}^2,{\stackrel{\‾}{A}}_l=P_{\stackrel{\‾}{A}}^{\⊥}{A}_l\∞$

$P_{\stackrel{\‾}{A}}^{\⊥}=I-\stackrel{\‾}{A}{\left({\stackrel{\‾}{A}}^T{Q}_y^{-1}\stackrel{\‾}{A}\right)}^{-1}{\stackrel{\‾}{A}}^T{Q}_y^{-1}$

$P_{{\stackrel{\‾}{A}}_l}={\stackrel{\‾}{A}}_l{\left({\stackrel{\‾}{A}}_l^T{Q}_y^{-1}{\stackrel{\‾}{A}}_l\right)}^{-1}{\stackrel{\‾}{A}}_l^T{Q}_y^{-1}$

其中，I为单位矩阵，是L2范数运算符，Q_y是GPS测站坐标时间序列观测值的方差协方差矩阵； $\stackrel{\‾}{A}=[A,A_1,...,A_{j-1}],$ $A_l=\left[\begin{array}{cc}\cos {\mathrm{\ω}}_l{t}_1& \sin {\mathrm{\ω}}_l{t}_1\\ \cos {\mathrm{\ω}}_l{t}_2& \sin {\mathrm{\ω}}_l{t}_2\\ .& .\\ .& .\\ .& .\\ \cos {\mathrm{\ω}}_l{t}_m& \sin {\mathrm{\ω}}_l{t}_m\end{array}\right],$ l＝1,2,...j-1；

4.2采用假设 $\left\{\begin{array}{c}{H}_0:y=\mathrm{Ax}+\underset{k=1}{\overset{j-1}{\mathrm{\Σ}}}\left(A_k{x}_k\right)\\ H_a:y=\mathrm{Ax}+\underset{k=1}{\overset{j}{\mathrm{\Σ}}}\left(A_k{x}_k\right)\end{array}\right.$ 判断是否向谐函数模型增加新频率，即根据步骤4.1获得的探测频率ω_j及其对应的A_j判定是否接受零假设，若接受，则增加新频率，然后执行步骤4.1；否则，执行步骤4.3；H₀为零假设，H_a为备选假设；

4.3采用Q_y＝σ²I和统计量评估假设检验，其中， 为零假设下的最小二乘残差；为备选假设下的验后方差，为备选假设下的最小二乘残差，统计量T₂≈F(2,m-n-2s)，s为假设检验过程中确定的备选谐函数数量。

5.如权利要求1所述的GPS测站坐标时间序列周期性探测方法，其特征是：

步骤5具体为：

采用最小二乘法，通过解算以下矩阵获得正弦函数分量a_k和余弦函数分量b_k：

y＝Bx+v

$\hat{x}={\left(B^T{Q}_y^{-1}B\right)}^{-1}{B}^T{Q}_y^{-1}y$

$Q_{\hat{x}}={\left(B^T{Q}_Y^{-1}B\right)}^{-1}$

其中，矩阵B＝[A A_k]，k＝1,...,j；为x的最小二乘估值；Q_y是GPS测站坐标时间序列观测值的方差协方差矩阵；y是GPS测站坐标时间序列观测值，v是随机误差矢量，采用最小二乘法解算获得。

6.GPS测站坐标时间序列周期性探测系统，其特征是，包括：

观测值获得模块，用来获取GPS测站坐标时间序列观测值；

修正模块，用来剔除GPS测站坐标时间序列观测值粗差并修正天线相位中心偏差；

谐函数模型构建模块，用来采用谐函数描述GPS测站坐标时间序列获得GPS测站坐标时间序列的谐函数模型，并设计谐函数模型矩阵其中，y表示GPS测站坐标时间序列观测值； $A=\left[\begin{array}{cc}1& t_1\\ 1& t_2\\ ...& ...\\ 1& t_m\end{array}\right],$ $A_k=\left[\begin{array}{cc}\cos {\mathrm{\ω}}_k{t}_1& \sin {\mathrm{\ω}}_k{t}_1\\ \cos {\mathrm{\ω}}_k{t}_2& \sin {\mathrm{\ω}}_k{t}_2\\ .& .\\ .& .\\ .& .\\ \cos {\mathrm{\ω}}_k{t}_m& \sin {\mathrm{\ω}}_k{t}_m\end{array}\right],$ t_i为观测历元，i＝1,2,...,m，m为观测历元数量，ω_k表示可变周期性信号的频率； $x_k=\left[\begin{array}{c}{a}_k\\ b_k\end{array}\right],$ a_k、b_k分别为正弦函数分量、余弦函数分量； $x=\left[\begin{array}{c}y\left(0\right)\\ a_0\end{array}\right],$ y(0)为GPS测站参考历元坐标，a₀为GPS测站的线性速度；k表示谐函数编号，q为谐函数数量；

频率估计模块，用来采用最小二乘法解算谐函数模型矩阵，获得可变周期性信号的频率估计值，并采用假设检验方法验证频率估计值；

周期性探测模块，用来基于假设检验验证后的频率估计值，采用最小二乘方法解算谐函数模型矩阵获得正弦函数分量a_k、余弦函数分量b_k和线性速度a₀。