CN110398753B - Gnss测站坐标时间序列周期性探测方法及系统 - Google Patents
Gnss测站坐标时间序列周期性探测方法及系统 Download PDFInfo
- Publication number
- CN110398753B CN110398753B CN201910579540.2A CN201910579540A CN110398753B CN 110398753 B CN110398753 B CN 110398753B CN 201910579540 A CN201910579540 A CN 201910579540A CN 110398753 B CN110398753 B CN 110398753B
- Authority
- CN
- China
- Prior art keywords
- frequency
- matrix
- alternative
- time sequence
- coordinate time
- 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
Links
- 230000000737 periodic effect Effects 0.000 title claims abstract description 106
- 238000001514 detection method Methods 0.000 title claims abstract description 38
- 239000011159 matrix material Substances 0.000 claims abstract description 102
- 238000000034 method Methods 0.000 claims abstract description 38
- 238000012360 testing method Methods 0.000 claims abstract description 26
- 238000012163 sequencing technique Methods 0.000 claims abstract description 7
- 238000010183 spectrum analysis Methods 0.000 claims abstract description 7
- 238000012795 verification Methods 0.000 claims abstract description 7
- 238000013461 design Methods 0.000 claims description 25
- 238000004458 analytical method Methods 0.000 claims description 8
- 238000012937 correction Methods 0.000 claims description 7
- 241001061260 Emmelichthys struhsakeri Species 0.000 claims description 5
- 238000010276 construction Methods 0.000 claims description 5
- 238000004364 calculation method Methods 0.000 claims description 4
- 230000008569 process Effects 0.000 claims description 3
- 230000001932 seasonal effect Effects 0.000 description 4
- 238000013178 mathematical model Methods 0.000 description 3
- 230000002159 abnormal effect Effects 0.000 description 2
- 230000008859 change Effects 0.000 description 2
- 230000007547 defect Effects 0.000 description 2
- 230000007613 environmental effect Effects 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 230000004044 response Effects 0.000 description 2
- 201000004569 Blindness Diseases 0.000 description 1
- 241000764238 Isis Species 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 230000002301 combined effect Effects 0.000 description 1
- 230000007812 deficiency Effects 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 230000001747 exhibiting effect Effects 0.000 description 1
- 230000002349 favourable effect Effects 0.000 description 1
- 230000007774 longterm Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000005086 pumping Methods 0.000 description 1
- 238000007670 refining Methods 0.000 description 1
- 239000011435 rock Substances 0.000 description 1
- 239000000523 sample Substances 0.000 description 1
- 239000007787 solid Substances 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 239000002352 surface water Substances 0.000 description 1
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
- G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
- G01S19/03—Cooperating elements; Interaction or communication between different cooperating elements or between cooperating elements and receivers
- G01S19/08—Cooperating elements; Interaction or communication between different cooperating elements or between cooperating elements and receivers providing integrity information, e.g. health of satellites or quality of ephemeris data
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
- G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
- G01S19/13—Receivers
- G01S19/20—Integrity monitoring, fault detection or fault isolation of space segment
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
- G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
- G01S19/13—Receivers
- G01S19/24—Acquisition or tracking or demodulation of signals transmitted by the system
- G01S19/243—Demodulation of navigation message
Landscapes
- Engineering & Computer Science (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Computer Networks & Wireless Communication (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Computer Security & Cryptography (AREA)
- Position Fixing By Use Of Radio Waves (AREA)
Abstract
本发明提供一种GNSS测站坐标时间序列周期性探测方法及系统，包括获取GPS测站坐标时间序列观测值，剔除粗差并修正天线相位中心偏差；采用周期图法对GPS测站坐标时间序列进行初步频谱分析，获取对应测站若干个主要周期频率，按照主频的振幅大小进行排序；采用谐函数描述GPS测站坐标时间序列，获得GPS测站坐标时间序列的谐函数模型，并构建谐函数模型矩阵；采用主频作为先验约束，获得多个周期性信号的备选频率；基于最小二乘准则解算谐函数模型，获取优选的频率并采用假设检验方法验证备选频率，构建包含多个优选频率的谐函数模型；根据假设检验验证后的备选频率，基于最小二乘准则解算谐函数模型矩阵，获得对任一GNSS测站周期性信号的探测结果。
Description
技术领域
本发明属于GNSS数据精密处理技术领域，具体涉及一种GNSS测站坐标时间序列周期性探测方法及系统。
背景技术
GNSS测站坐标时间序列反映了测站位置随时间变化的规律特征。现有的GNSS测站坐标时间序列预测估计模型通常包含长期趋势项和季节性信号。目前最广泛的做法为，将周期性信号简单地作为周年信号和半周年信号及其组合作为常数参数广泛应用(BlewittandLavallée,2002)，一般采用周年周期及其谐波和常数振幅的正弦曲线建模(Nicolaidis,2002)。
GNSS测站坐标时间序列中存在的季节性信号是对环境变化的响应。然而，数据处理分析策略与GPS相关的技术因素、固体地球质量负荷、大气压负荷、水储量变化、岩石热膨胀、降雨、多路径效应等因素的综合影响，导致GNSS测站坐标时间序列中的季节性信号呈现周期性可变特征(Davis,Wernicke et al.,2012)。田云锋(2011)对青藏高原和喜马拉雅地区的GPS垂向分量分析发现存在明显的相位变化，推断是受到地表水体因子控制；而在CMONOC及周边IGS台站的位置时间序列中均发现了周期约为351/n(n＝1,…,6)天的“异常”周期项，并认为地表质量负荷并非该异常周期项的来源。
目前GNSS测站坐标时间序列周期性探测方法存在不足：将周期性信号简单地当作周年信号和半周年信号及其组合进行描述，对速率估值及其不确定性都会产生明显影响。针对此项不足，Amiri-Simkooei et al.,(2007)提出了利用最小二乘谐波(Least SquareHarmonic Estimation，简称LSHE)估计方法对周期性信号进行探测。然而，采用LSHE方法存在两大问题：(1)理论上而言，不采用任何先验约束的情况下，周期信号及其谐波信号有无限多个组合的可能性，这将造成估计结果的可靠性不强；(2)任意谐波信号的组合，周期性信号频率相近时，易出现奇异情况，可能造成相近频率周期性信号的判定错误，最终导致周期性信号建模的错误。
发明内容
针对上述现有技术存在的两大问题，本发明提供了一种GNSS测站坐标时间序列周期性探测方法及系统，以周期图法结果作为约束信息，减少组合周期性信号的备选项；顾及线性模型系数矩阵的结构，解决临近频率的奇异问题；并利用谐函数和最小二乘准则，实现GPS测站坐标时间序列周期性探测方法。
本发明所采用的技术方案提供一种GNSS测站坐标时间序列周期性探测方法，包括以下步骤，
步骤1，获取GPS测站坐标时间序列观测值；
步骤2，剔除GPS测站坐标时间序列观测值粗差并修正天线相位中心偏差，得到修正后的GPS测站坐标时间序列；
步骤3，采用周期图法对步骤2中获取的GPS测站坐标时间序列进行初步频谱分析，获取对应测站若干个主要周期频率，按照主频的振幅大小进行排序；
步骤4，采用谐函数描述GPS测站坐标时间序列，获得GPS测站坐标时间序列的谐函数模型，并构建谐函数模型矩阵：
步骤5，采用步骤3中主频作为先验约束，获得多个周期性信号的备选频率；基于最小二乘准则解算谐函数模型，获取优选的频率并采用假设检验方法验证备选频率，构建包含多个优选频率的谐函数模型；
步骤6，根据假设检验验证后的备选频率，基于最小二乘准则解算谐函数模型矩阵获得正弦函数分量ak、余弦函数分量bk和线性速度a0，获得对任一GNSS测站周期性信号的探测结果。
而且，步骤1中，通过双差定位软件工具、精密单点定位软件工具或IGS分析中心获取积累的GPS测站坐标时间序列观测值。
而且，步骤4中，所述的GPS测站坐标时间序列的谐函数模型为，
其中，y(ti)为观测历元ti对应的GPS测站坐标观测值，y(0)+a0ti为线性运动部分，y(0)为GPS测站参考历元坐标，a0为GPS测站的线性速度；
为周期性信号部分，ak、bk分别为正弦函数分量、余弦函数分量，ωk表示周期性信号的频率；k表示谐波函数编号，q为谐波函数数量。
所述谐函数模型矩阵为，
其中，y表示GNSS测站坐标时间序列观测值；
表示线性运动部分的设计矩阵；
表示周期性信号部分的设计矩阵，ti为观测历元，i＝1,2,...,m，m为观测历元数量，ωk表示周期性信号的频率；
表示周期性信号部分的待定未知参数，ak、bk分别为正弦函数分量、余弦函数分量；
表示线性运动部分的待定未知参数，y(0)为GPS测站参考历元坐标，a0为GPS测站的线性速度；k表示谐函数编号，q为谐函数数量；
而且，步骤5包括以下子步骤，
步骤5.1，确定周期性信号的备选频率如下，
基于假设
计算各备择假设Ha下模型设计矩阵的条件数：
cond([M|Aj])＝||[M|Aj]||·||[M|Aj]-1||
其中，H0为零假设，Ha为备选假设；j-1为零假设下模型中已有的谐函数数量；Al、xl表示第l个谐函数的设计矩阵与待定未知参数；cond()表示求取括号内矩阵的条件数，M＝[A A1 … Aj-1]，表示零假设下模型的设计矩阵，Ajxj表示备选频率在备选假设模型中新增的谐函数。
若
小于预设的阈值，则认为对应备选频率的引入可能导致矩阵奇异，将该备选频率删去；
步骤5.2，基于不同的备选频率，改变备选假设模型，并通过解算最大化问题以探测优选的频率ωj及对应的Aj；
步骤5.3，根据步骤5.1获得的探测频率ωj及其对应的Aj，计算统计量T2；
设
采用统计量T2进行假设检验，判断是否向谐函数模型增加新频率；T2～F(2,m-2-2s)，s为备选假设下谐波函数的数量，若计算得到的T2值大于的F分布显著水平，判定接受备选假设，增加频率ωj并返回步骤5.1；否则，执行步骤6。
而且，步骤6中，采用最小二乘法，通过解算以下矩阵获得正弦函数分量ak和余弦函数分量bk，
y＝Bx+v
其中，矩阵B＝[A A1 … Aj]；
为x的最小二乘估值；Qy是GPS测站坐标时间序列观测值的方差协方差矩阵，
为协因数矩阵；y是GPS测站坐标时间序列观测值，v是随机误差矢量，采用最小二乘准则解算获得。
本发明还提供一种GNSS测站坐标时间序列周期性探测系统，包括以下模块，
观测值获得模块，用于获取GPS测站坐标时间序列观测值；
修正模块，用于剔除GPS测站坐标时间序列观测值粗差并修正天线相位中心偏差，得到修正后的GPS测站坐标时间序列；
主频分析模块，用于采用周期图法对修正模块中获取的GPS测站坐标时间序列进行初步频谱分析，获取对应测站若干个主要周期频率，按照主频的振幅大小进行排序；
谐函数模型构建模块，用于采用谐函数描述GPS测站坐标时间序列，获得GPS测站坐标时间序列的谐函数模型，并构建谐函数模型矩阵：
频率估计模块，用于采用主频分析模块所得主频作为先验约束，获得多个周期性信号的备选频率；基于最小二乘准则解算谐函数模型，获取优选的频率并采用假设检验方法验证备选频率，构建包含多个优选频率的谐函数模型；
周期性探测模块，用根据假设检验验证后的备选频率，基于最小二乘准则解算谐函数模型矩阵获得正弦函数分量ak、余弦函数分量bk和线性速度a0，获得对任一GNSS测站周期性信号的探测结果。
而且，观测值获得模块中，通过双差定位软件工具、精密单点定位软件工具或IGS分析中心获取积累的GPS测站坐标时间序列观测值。
而且，谐函数模型构建模块中，所述的GPS测站坐标时间序列的谐函数模型为，
其中，y(ti)为观测历元ti对应的GPS测站坐标观测值，y(0)+a0ti为线性运动部分，y(0)为GPS测站参考历元坐标，a0为GPS测站的线性速度；
为周期性信号部分，ak、bk分别为正弦函数分量、余弦函数分量，ωk表示周期性信号的频率；k表示谐波函数编号，q为谐波函数数量。
所述谐函数模型矩阵为，
其中，y表示GNSS测站坐标时间序列观测值；
表示线性运动部分的设计矩阵；
表示周期性信号部分的设计矩阵，ti为观测历元，i＝1,2,...,m，m为观测历元数量，ωk表示周期性信号的频率；
表示周期性信号部分的待定未知参数，ak、bk分别为正弦函数分量、余弦函数分量；
表示线性运动部分的待定未知参数，y(0)为GPS测站参考历元坐标，a0为GPS测站的线性速度；k表示谐函数编号，q为谐函数数量；
而且，频率估计模块执行以下子步骤，
步骤5.1，确定周期性信号的备选频率如下，
基于假设
计算各备择假设Ha下模型设计矩阵的条件数：
cond([M|Aj])＝||[M|Aj]||·||[M|Aj]-1||
其中，H0为零假设，Ha为备选假设；j-1为零假设下模型中已有的谐函数数量；Al、xl表示第l个谐函数的设计矩阵与待定未知参数；cond()表示求取括号内矩阵的条件数，M＝[A A1 … Aj-1]，表示零假设下模型的设计矩阵，Ajxj表示备选频率在备选假设模型中新增的谐函数。
若
小于预设的阈值，则认为对应备选频率的引入可能导致矩阵奇异，将该备选频率删去；
步骤5.2，基于不同的备选频率，改变备选假设模型，并通过解算最大化问题以探测优选的频率ωj及对应的Aj；
步骤5.3，根据步骤5.1获得的探测频率ωj及其对应的Aj，计算统计量T2；
设
采用统计量T2进行假设检验，判断是否向谐函数模型增加新频率；T2～F(2,m-2-2s)，s为备选假设下谐波函数的数量，若计算得到的T2值大于的F分布显著水平，判定接受备选假设，增加频率ωj并返回步骤5.1；否则，命令周期性探测模块工作。
而且，周期性探测模块中，采用最小二乘法，通过解算以下矩阵获得正弦函数分量ak和余弦函数分量bk，
y＝Bx+v
其中，矩阵B＝[A A1 … Aj]；
为x的最小二乘估值；Qy是GPS测站坐标时间序列观测值的方差协方差矩阵，
为协因数矩阵；y是GPS测站坐标时间序列观测值，v是随机误差矢量，采用最小二乘准则解算获得。
与现有技术相比，本发明具有以下区别特点：
本发明采用周期图法初步获取备选周期信号的主频信息，采用谐函数建立GPS测站坐标时间序列预估模型，有效避免了单一采用谐函数进行周期信号频率估计的盲目性，而且也突破了现有技术中仅采用周年加半周年信号组合的周期性信号描述坐标时间序列的现状；建立顾及真实情况的周期性信号模型，一方面有助于进一步精化已有的GPS测站坐标时间序列的速度模型，提高速度估值的可靠性；另一方面，包含完整周期性信号模型有助于真实反映与当地环境效应相关的(如降雨、温度、地表负荷及含水层抽水等)季节性地球物理信号，以进一步开展地球物理过程等方面的建模及解释。
现有GPS测站坐标时间序列周期性探测中，将周期性信号当作周年与半周年信号进行估计，这样导致GPS测站坐标时间序列速度模型的不完善、速率估值及其不确定性的较大偏差，本发明可解决上述问题，可处理含有时变周期性信号的GPS测站坐标时间序列，并保证周期性探测的准确性。
附图说明
图1为本发明实施例的流程图。
具体实施方式
下面结合附图和实施例对本发明做进一步说明。
本发明利用周期图法结果作为先验信息，采用谐函数对GNSS测站坐标时间序列建模，通过最小二乘准则估计GNSS测站多个周期性信号，从而实现GNSS测站坐标时间序列周期性探测。
参见图1，本发明实施例提供的一种GPS测站坐标时间序列周期性探测方法，包括步骤：
步骤1，获取GPS测站坐标时间序列观测值；
具体实施时，可以通过双差定位软件工具、精密单点定位软件工具或IGS分析中心获取积累的GPS测站坐标时间序列观测值。
步骤2，剔除GPS测站坐标时间序列观测值粗差并修正天线相位中心偏差，从而获取干净的GPS测站坐标时间序列，即修正后的GPS测站坐标时间序列；
具体实施时，可以基于拉依达准则、IQR准则等方法对GPS测站坐标时间序列观测值粗差进行探测与剔除；天线相位中心偏差即硬件相关的天线相位中心改正，可以通过IGS发布的天线相位中心参数，使用模型进行修正更换。
步骤3，采用周期图法对步骤2中获取的GPS测站坐标时间序列进行初步频谱分析，获取对应测站的前N个主要周期频率，按照主频的振幅大小进行排序；
周期图法是成熟方法，可用于非均匀采样数据的频率获取。然而该方法获得的低频信号会影响高频信号的获取，因此本专利提出，将其结果作为先验信息，与最小二乘谐波估计方法相结合，同时对高频与低频信号进行建模，得到完整的GNSS坐标时间序列函数模型。
具体实施时，N可采用预设的取值。实施例中，N优选取值20。
步骤4，采用谐函数描述GPS测站坐标时间序列，获得GPS测站坐标时间序列的谐函数模型，并构建谐函数模型矩阵：
所述的GPS测站坐标时间序列的谐函数模型为：
其中，y(ti)为观测历元ti对应的GPS测站坐标观测值，y(0)+a0ti为线性运动部分，y(0)为GPS测站参考历元坐标，a0为GPS测站的线性速度；
为周期性信号部分，ak、bk分别为正弦函数分量、余弦函数分量，ωk表示周期性信号的频率；k表示谐波函数编号，q为谐波函数数量。
所述谐函数模型矩阵为：
其中，y表示GNSS测站坐标时间序列观测值；
表示线性运动部分的设计矩阵；
表示周期性信号部分的设计矩阵，ti为观测历元，i＝1,2,...,m，m为观测历元数量，ωk表示周期性信号的频率；
表示周期性信号部分的待定未知参数，ak、bk分别为正弦函数分量、余弦函数分量；
表示线性运动部分的待定未知参数，y(0)为GPS测站参考历元坐标，a0为GPS测站的线性速度；k表示谐函数编号，q为谐函数数量；
本步骤实现构建GNSS测站坐标时间序列的数学模型，包括参考历元的坐标、线性速度和周期性信号，确立数学模型的涉及矩阵及对应的待定未知参数。
步骤5，采用步骤3中主频作为先验约束，获得多个周期性信号的备选频率；基于最小二乘准则解算谐函数模型，获取优选的频率并采用假设检验方法验证备选频率，构建包含多个优选频率的谐函数模型；
实施例中，步骤5进一步包括子步骤：
5.1确定周期性信号的备选频率：
基于假设
计算各备择假设Ha下模型设计矩阵的条件数：
cond([M|Aj])＝||[M|Aj]||·||[M|Aj]-1||
其中，H0为零假设，Ha为备选假设；j-1为零假设下模型中已有的谐函数数量，即其最初取值为1，并随谐函数增加递增；Al、xl表示第l个谐函数的设计矩阵与待定未知参数；cond()表示求取括号内矩阵的条件数，M＝[A A1 … Aj-1]，表示零假设下模型的设计矩阵，Ajxj表示备选频率在备选假设模型中新增的谐函数。
若
小于设定的阈值(具体实施时预先设置，可优选设为10-6)，则认为对应备选频率的引入可能导致矩阵奇异，将该备选频率删去；
5.2基于最小二乘准则，获取优选频率：基于不同的备选频率，改变备选假设模型，并通过解算以下最大化问题以探测优选的频率ωj及其对应的Aj：
其中，
表示优选的ωj值应使
最大化；
表示备选假设模型中待定未知参数xj的估值，
为其转置矩阵，
为矩阵Aj的转置矩阵；N0为平差因子阵，I表示单位矩阵，
为备选假设下的验后方差，Qy是GPS测站坐标时间序列观测值的方差协方差矩阵，
为其逆矩阵。
5.3基于假设检验，判断是否新增周期信号的频率：基于假设
采用统计量T2进行假设检验，判断是否向谐函数模型增加新频率，即根据步骤5.1获得的探测频率ωj及其对应的Aj，计算统计量T2；由于T2～F(2,m-2-2s)，s为备选假设下谐波函数的数量，即该统计量服从F分布，
若计算得到的T2值大于的F分布显著水平，判定接受备选假设，增加频率ωj并返回步骤5.1；
否则，执行步骤6。
T2的计算公式如下：
其中，
与前文相同，
为备选假设下验后方差的估值，
为备选假设下的最小二乘残差。
步骤6，计算并评估任一测站的所有未知参数，实现对任一GNSS测站周期性信号频率及其对应振幅、线性速度等参数的估计：
实施例中，基于假设检验验证后的备选频率，基于最小二乘准则解算谐函数模型矩阵获得正弦函数分量ak、余弦函数分量bk和线性速度a0。据此可以获得对任一GNSS测站周期性信号的探测结果。
实施例中步骤6中，进行以下操作：
采用最小二乘法，通过解算以下矩阵获得正弦函数分量ak和余弦函数分量bk：
y＝Bx+v
其中，矩阵B＝[A A1 … Aj]；
为x的最小二乘估值；Qy是GPS测站坐标时间序列观测值的方差协方差矩阵，
为协因数矩阵；y是GPS测站坐标时间序列观测值，v是随机误差矢量，采用最小二乘准则解算获得。
具体实施时，本发明所提供方法可基于软件技术实现自动运行流程，也可采用模块化方式实现相应系统。
实施例提供的GPS测站坐标时间序列周期性探测系统，包括以下模块，
观测值获得模块，用于获取GPS测站坐标时间序列观测值；
修正模块，用于剔除GPS测站坐标时间序列观测值粗差并修正天线相位中心偏差，得到修正后的GPS测站坐标时间序列；
主频分析模块，用于采用周期图法对修正模块中获取的GPS测站坐标时间序列进行初步频谱分析，获取对应测站若干个主要周期频率，按照主频的振幅大小进行排序；
谐函数模型构建模块，用于采用谐函数描述GPS测站坐标时间序列，获得GPS测站坐标时间序列的谐函数模型，并构建谐函数模型矩阵：
频率估计模块，用于采用主频分析模块所得主频作为先验约束，获得多个周期性信号的备选频率；基于最小二乘准则解算谐函数模型，获取优选的频率并采用假设检验方法验证备选频率，构建包含多个优选频率的谐函数模型；
周期性探测模块，用根据假设检验验证后的备选频率，基于最小二乘准则解算谐函数模型矩阵获得正弦函数分量ak、余弦函数分量bk和线性速度a0，获得对任一GNSS测站周期性信号的探测结果。
特别需要说明的是，当候选频率非常接近时，会出现矩阵病态的情况，在周期性探测模块需要首先选择振幅较大的那一个候选频率进行计算，如果满足假设检验条件，作为新增周期性信号，并且带入原有的时间序列进行建模，并获取残差时间序列，然后，将振幅较小的那一个候选频率再次重复假设检验，判断是否作为新增周期性信号。这种处理，主要是为了避免当候选频率接近时出现的奇异情况。
应当理解的是，本说明书未详细阐述的部分均属于现有技术。
应当理解的是，上述针对较佳实施例的描述较为详细，并不能因此而认为是对本发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下，在不脱离本发明权利要求所保护的范围情况下，还可以做出替换或变形，均落入本发明的保护范围之内，本发明的请求保护范围应以所附权利要求为准。
Claims (8)
1.一种GNSS测站坐标时间序列周期性探测方法，其特征在于：采用周期图法初步获取备选周期信号的主频信息，包括以下步骤，
步骤1，获取GPS测站坐标时间序列观测值；
步骤2，剔除GPS测站坐标时间序列观测值粗差并修正天线相位中心偏差，得到修正后的GPS测站坐标时间序列；
步骤3，采用周期图法对步骤2中获取的GPS测站坐标时间序列进行初步频谱分析，获取对应测站若干个主要周期频率，按照主频的振幅大小进行排序，所得结果用于作为先验信息，与最小二乘谐波估计方法相结合，同时对高频与低频信号进行建模，得到完整的GNSS坐标时间序列函数模型；
步骤4，采用谐函数描述GPS测站坐标时间序列，获得GPS测站坐标时间序列的谐函数模型，并构建谐函数模型矩阵：
步骤5，采用步骤3中主频作为先验约束，获得多个周期性信号的备选频率；基于最小二乘准则解算谐函数模型，获取优选的频率并采用假设检验方法验证备选频率，构建包含多个优选频率的谐函数模型；实现过程如下，
5.1)确定周期性信号的备选频率：
基于假设
计算各备择假设Ha下模型设计矩阵的条件数：
cond([M|Aj])＝||[M|Aj]||·||[M|Aj]-1||
其中，H0为零假设，Ha为备选假设；j-1为零假设下模型中已有的谐函数数量；Al、xl表示第l个谐函数的设计矩阵与待定未知参数；cond()表示求取括号内矩阵的条件数，M＝[AA1...j-1]，表示零假设下模型的设计矩阵，Ajxj表示备选频率在备选假设模型中新增的谐函数；
若
小于设定的阈值，则认为对应备选频率的引入可能导致矩阵奇异，将该备选频率删去；
5.2)基于最小二乘准则，获取优选频率，包括基于不同的备选频率，改变备选假设模型，并通过解算以下最大化问题以探测优选的频率ωj及对应的Aj：
其中，
表示优选的ωj值应使
最大化；
表示备选假设模型中待定未知参数xj的估值，
为其转置矩阵，
为矩阵Aj的转置矩阵；N0为平差因子阵，I表示单位矩阵，
为备选假设下的验后方差，Qy是GPS测站坐标时间序列观测值的方差协方差矩阵，
为其逆矩阵；
5.3)基于假设检验，判断是否新增周期信号的频率：
基于假设
采用统计量T2进行假设检验，判断是否向谐函数模型增加新频率，即根据步骤5.1获得的探测频率ωj及其对应的Aj，计算统计量T2；T2～F(2,m-2-2s)，s为备选假设下谐波函数的数量，该统计量服从F分布，
若计算得到的T2值大于的F分布显著水平，判定接受备选假设，增加频率ωj并返回步骤5.1；
否则，执行步骤6；
T2的计算公式如下：
其中，
为备选假设下验后方差的估值，
为备选假设下的最小二乘残差；
步骤6，根据假设检验验证后的备选频率，基于最小二乘准则解算谐函数模型矩阵获得正弦函数分量ak、余弦函数分量bk和线性速度a0，获得对任一GNSS测站周期性信号的探测结果。
2.根据权利要求1所述GNSS测站坐标时间序列周期性探测方法，其特征在于：步骤1中，通过双差定位软件工具、精密单点定位软件工具或IGS分析中心获取积累的GPS测站坐标时间序列观测值。
3.根据权利要求1或2所述GNSS测站坐标时间序列周期性探测方法，其特征在于：步骤4中，所述的GPS测站坐标时间序列的谐函数模型为，
其中，y(ti)为观测历元ti对应的GPS测站坐标观测值，y(0)+a0ti为线性运动部分，y(0)为GPS测站参考历元坐标，a0为GPS测站的线性速度；
为周期性信号部分，ak、bk分别为正弦函数分量、余弦函数分量，ωk表示周期性信号的频率；k表示谐波函数编号，q为谐波函数数量；
所述谐函数模型矩阵为，
其中，y表示GNSS测站坐标时间序列观测值；
表示线性运动部分的设计矩阵；
表示周期性信号部分的设计矩阵，ti为观测历元，i＝1,2,...,m，m为观测历元数量，ωk表示周期性信号的频率；
表示周期性信号部分的待定未知参数，ak、bk分别为正弦函数分量、余弦函数分量；
表示线性运动部分的待定未知参数，y(0)为GPS测站参考历元坐标，a0为GPS测站的线性速度；k表示谐函数编号，q为谐函数数量。
4.根据权利要求3所述GNSS测站坐标时间序列周期性探测方法，其特征在于：步骤6中，采用最小二乘法，通过解算以下矩阵获得正弦函数分量ak和余弦函数分量bk，
y＝Bx+v
其中，矩阵B＝[A A1…Aj]；
为x的最小二乘估值；Qy是GPS测站坐标时间序列观测值的方差协方差矩阵，
为协因数矩阵；y是GPS测站坐标时间序列观测值，v是随机误差矢量，采用最小二乘准则解算获得。
5.一种GNSS测站坐标时间序列周期性探测系统，其特征在于：采用周期图法初步获取备选周期信号的主频信息，包括以下模块，
观测值获得模块，用于获取GPS测站坐标时间序列观测值；
修正模块，用于剔除GPS测站坐标时间序列观测值粗差并修正天线相位中心偏差，得到修正后的GPS测站坐标时间序列，所得结果用于作为先验信息，与最小二乘谐波估计方法相结合，同时对高频与低频信号进行建模，得到完整的GNSS坐标时间序列函数模型；
主频分析模块，用于采用周期图法对修正模块中获取的GPS测站坐标时间序列进行初步频谱分析，获取对应测站若干个主要周期频率，按照主频的振幅大小进行排序；
谐函数模型构建模块，用于采用谐函数描述GPS测站坐标时间序列，获得GPS测站坐标时间序列的谐函数模型，并构建谐函数模型矩阵：
频率估计模块，用于采用主频分析模块所得主频作为先验约束，获得多个周期性信号的备选频率；基于最小二乘准则解算谐函数模型，获取优选的频率并采用假设检验方法验证备选频率，构建包含多个优选频率的谐函数模型；实现过程如下，
5.1)确定周期性信号的备选频率：
基于假设
计算各备择假设Ha下模型设计矩阵的条件数：
cond([M|Aj])＝||[M|Aj]||·||[M|Aj]-1||
其中，H0为零假设，Ha为备选假设；j-1为零假设下模型中已有的谐函数数量；Al、xl表示第l个谐函数的设计矩阵与待定未知参数；cond()表示求取括号内矩阵的条件数，M＝[AA1…Aj-1]，表示零假设下模型的设计矩阵，Ajxj表示备选频率在备选假设模型中新增的谐函数；
若
小于设定的阈值，则认为对应备选频率的引入可能导致矩阵奇异，将该备选频率删去；
5.2)基于最小二乘准则，获取优选频率，包括基于不同的备选频率，改变备选假设模型，并通过解算以下最大化问题以探测优选的频率ωj及对应的Aj：
其中，
表示优选的ωj值应使
最大化；
表示备选假设模型中待定未知参数xj的估值，
为其转置矩阵，
为矩阵Aj的转置矩阵；N0为平差因子阵，I表示单位矩阵，
为备选假设下的验后方差，Qy是GPS测站坐标时间序列观测值的方差协方差矩阵，
为其逆矩阵；
5.3)基于假设检验，判断是否新增周期信号的频率：
基于假设
采用统计量T2进行假设检验，判断是否向谐函数模型增加新频率，即根据步骤5.1获得的探测频率ωj及其对应的Aj，计算统计量T2；T2～F(2,m-2-2s)，s为备选假设下谐波函数的数量，该统计量服从F分布，
若计算得到的T2值大于的F分布显著水平，判定接受备选假设，增加频率ωj并返回步骤5.1；
否则，执行步骤6；
T2的计算公式如下：
其中，
为备选假设下验后方差的估值，
为备选假设下的最小二乘残差；
周期性探测模块，用根据假设检验验证后的备选频率，基于最小二乘准则解算谐函数模型矩阵获得正弦函数分量ak、余弦函数分量bk和线性速度a0，获得对任一GNSS测站周期性信号的探测结果。
6.根据权利要求5所述GNSS测站坐标时间序列周期性探测系统，其特征在于：观测值获得模块中，通过双差定位软件工具、精密单点定位软件工具或IGS分析中心获取积累的GPS测站坐标时间序列观测值。
7.根据权利要求5或6所述GNSS测站坐标时间序列周期性探测系统，其特征在于：谐函数模型构建模块中，所述的GPS测站坐标时间序列的谐函数模型为，
其中，y(ti)为观测历元ti对应的GPS测站坐标观测值，y(0)+a0ti为线性运动部分，y(0)为GPS测站参考历元坐标，a0为GPS测站的线性速度；
为周期性信号部分，ak、bk分别为正弦函数分量、余弦函数分量，ωk表示周期性信号的频率；k表示谐波函数编号，q为谐波函数数量；
所述谐函数模型矩阵为，
其中，y表示GNSS测站坐标时间序列观测值；
表示线性运动部分的设计矩阵；
表示周期性信号部分的设计矩阵，ti为观测历元，i＝1,2,...,m，m为观测历元数量，ωk表示周期性信号的频率；
表示周期性信号部分的待定未知参数，ak、bk分别为正弦函数分量、余弦函数分量；
表示线性运动部分的待定未知参数，y(0)为GPS测站参考历元坐标，a0为GPS测站的线性速度；k表示谐函数编号，q为谐函数数量。
8.根据权利要求7所述GNSS测站坐标时间序列周期性探测系统，其特征在于：周期性探测模块中，采用最小二乘法，通过解算以下矩阵获得正弦函数分量ak和余弦函数分量bk，
y＝Bx+v
其中，矩阵B＝[A A1…Aj]；
为x的最小二乘估值；Qy是GPS测站坐标时间序列观测值的方差协方差矩阵，
为协因数矩阵；y是GPS测站坐标时间序列观测值，v是随机误差矢量，采用最小二乘准则解算获得。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910579540.2A CN110398753B (zh) | 2019-06-28 | 2019-06-28 | Gnss测站坐标时间序列周期性探测方法及系统 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910579540.2A CN110398753B (zh) | 2019-06-28 | 2019-06-28 | Gnss测站坐标时间序列周期性探测方法及系统 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN110398753A CN110398753A (zh) | 2019-11-01 |
| CN110398753B true CN110398753B (zh) | 2023-06-06 |
Family
ID=68322602
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910579540.2A Active CN110398753B (zh) | 2019-06-28 | 2019-06-28 | Gnss测站坐标时间序列周期性探测方法及系统 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN110398753B (zh) |
Families Citing this family (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111241473B (zh) * | 2019-12-27 | 2023-09-29 | 中国空间技术研究院 | 一种提高区域地下水储量估计精度的方法 |
| CN111339483B (zh) * | 2020-01-19 | 2022-03-11 | 武汉大学 | 一种基于去趋势互相关分析的gnss影像生成方法 |
| CN111443366B (zh) * | 2020-04-28 | 2022-04-29 | 武汉大学 | 一种gnss区域网中异常点探测方法及系统 |
| CN111965669B (zh) * | 2020-08-14 | 2021-09-03 | 长江空间信息技术工程有限公司(武汉) | 一种gnss时间序列中观测墩热膨胀信号的分离方法 |
| CN111965670B (zh) * | 2020-08-14 | 2023-05-12 | 长江空间信息技术工程有限公司(武汉) | 量化gnss时间序列高频观测墩热膨胀信号混叠的方法 |
| CN112556563B (zh) * | 2020-11-30 | 2022-03-29 | 深圳大学 | 一种北斗定位长期监测数据的处理方法及系统 |
| CN112612822B (zh) * | 2020-12-11 | 2023-04-28 | 中铁第四勘察设计院集团有限公司 | 一种北斗坐标时间序列的预测方法、装置、设备和存储介质 |
| CN113341439B (zh) * | 2021-06-22 | 2022-04-15 | 武汉大学 | 一种顾及周期信号的gnss测站速度稳健估测方法 |
| CN114253962B (zh) * | 2022-03-02 | 2022-05-17 | 中国测绘科学研究院 | 一种考虑非线性因素的区域格网速度场构建方法及系统 |
| CN114922646B (zh) * | 2022-04-02 | 2024-05-14 | 中铁隧道局集团有限公司 | 超小半径缓和曲线段盾构割线始发施工方法 |
| CN116204756B (zh) * | 2023-04-28 | 2023-07-07 | 武汉大学 | 一种多分析中心精密站坐标产品综合方法及系统 |
Family Cites Families (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104765055A (zh) * | 2015-04-14 | 2015-07-08 | 武汉大学 | Gps测站坐标时间序列周期性探测方法及系统 |
| CN106772498B (zh) * | 2016-11-21 | 2019-05-10 | 华东交通大学 | 一种gps位置时间序列噪声模型建立方法 |
| CN106814378B (zh) * | 2017-01-17 | 2019-05-10 | 华东交通大学 | 一种gnss位置时间序列周期特性挖掘方法 |
-
2019
- 2019-06-28 CN CN201910579540.2A patent/CN110398753B/zh active Active
Also Published As
| Publication number | Publication date |
|---|---|
| CN110398753A (zh) | 2019-11-01 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110398753B (zh) | Gnss测站坐标时间序列周期性探测方法及系统 | |
| CN107132558B (zh) | 惯性辅助的多频多模gnss周跳修复方法及系统 | |
| CN109738917B (zh) | 一种北斗变形监测中的多路径误差削弱方法及装置 | |
| Wing et al. | Consumer-grade global positioning system (GPS) accuracy and reliability | |
| RU2591953C2 (ru) | Навигационная система и способ разрешения целочисленных неоднозначностей с использованием ограничения неоднозначности двойной разности | |
| CN109932735A (zh) | 北斗短基线单频单历元解算的定位方法 | |
| RU2615172C2 (ru) | Адаптивный способ для оценки электронного содержания ионосферы | |
| JP2011517770A (ja) | 衛星ナビゲーションシステムのインテグリティのリアルタイム監視用の装置および方法 | |
| CN105849589A (zh) | 全球导航卫星系统、定位终端、定位方法以及记录介质 | |
| CN102124365A (zh) | 具有质量测量的换算的gnss信号处理方法和装置 | |
| CN101950024B (zh) | 用于局域增强系统的码载一致性检测方法 | |
| CN106093991A (zh) | 一种用于gnss定位的模糊度快速恢复方法及系统 | |
| Gunning et al. | Multi-GNSS constellation anomaly detection and performance monitoring | |
| CN105425248B (zh) | 单频gnss相位稳定性监测的高频逐历元相位差方法 | |
| CN113758469A (zh) | 一种基于多模多频gnss接收机的潮位监测方法及系统 | |
| Varbla et al. | Assessment of marine geoid models by ship-borne GNSS profiles | |
| CN105652298B (zh) | 一种bds三频伪距相位组合的周跳探测修复方法及装置 | |
| CN114384557A (zh) | 星基增强系统的服务性能评估方法及装置 | |
| CN105511481A (zh) | 一种星载定轨优化方法 | |
| Fangchao et al. | A STEP CYCLE SLIP DETECTION AND REPAIR METHOD BASED ON DOUBLECONSTRAINT OF EPHEMERIS AND SMOOTHED PSEUDORANGE. | |
| CN112444825A (zh) | 一种基于北斗geo的电离层tec同化建模方法 | |
| CN109143286A (zh) | 一种顾及非模型化误差的卫星导航定位方法 | |
| CN114690210A (zh) | 一种基于多普勒观测值的北斗卫星机动探测方法 | |
| JP7329980B2 (ja) | 測位アルゴリズムの設定パラメータ決定方法 | |
| Yoon et al. | Multi-dimensional verification methodology of ionospheric gradient observation during plasma bubble events in the Brazilian region |
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 |