CN114253962A - 一种考虑非线性因素的区域格网速度场构建方法及系统 - Google Patents

一种考虑非线性因素的区域格网速度场构建方法及系统 Download PDF

Info

Publication number
CN114253962A
CN114253962A CN202210194847.2A CN202210194847A CN114253962A CN 114253962 A CN114253962 A CN 114253962A CN 202210194847 A CN202210194847 A CN 202210194847A CN 114253962 A CN114253962 A CN 114253962A
Authority
CN
China
Prior art keywords
time sequence
fitting
adopting
parameter
gnss coordinate
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.)
Granted
Application number
CN202210194847.2A
Other languages
English (en)
Other versions
CN114253962B (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.)
Chinese Academy of Surveying and Mapping
Original Assignee
Chinese Academy of Surveying and Mapping
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 Chinese Academy of Surveying and Mapping filed Critical Chinese Academy of Surveying and Mapping
Priority to CN202210194847.2A priority Critical patent/CN114253962B/zh
Publication of CN114253962A publication Critical patent/CN114253962A/zh
Application granted granted Critical
Publication of CN114253962B publication Critical patent/CN114253962B/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
    • G06F16/00—Information retrieval; Database structures therefor; File system structures therefor
    • G06F16/20—Information retrieval; Database structures therefor; File system structures therefor of structured data, e.g. relational data
    • G06F16/21—Design, administration or maintenance of databases
    • G06F16/215—Improving data quality; Data cleansing, e.g. de-duplication, removing invalid entries or correcting typographical errors
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F16/00—Information retrieval; Database structures therefor; File system structures therefor
    • G06F16/20—Information retrieval; Database structures therefor; File system structures therefor of structured data, e.g. relational data
    • G06F16/24—Querying
    • G06F16/245—Query processing
    • G06F16/2458—Special types of queries, e.g. statistical queries, fuzzy queries or distributed queries
    • G06F16/2474—Sequence data queries, e.g. querying versioned data
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F16/00—Information retrieval; Database structures therefor; File system structures therefor
    • G06F16/20—Information retrieval; Database structures therefor; File system structures therefor of structured data, e.g. relational data
    • G06F16/29—Geographical information databases
    • 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/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Databases & Information Systems (AREA)
  • Data Mining & Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Software Systems (AREA)
  • Probability & Statistics with Applications (AREA)
  • Computational Mathematics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Computational Linguistics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Evolutionary Biology (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Operations Research (AREA)
  • Fuzzy Systems (AREA)
  • Algebra (AREA)
  • Remote Sensing (AREA)
  • Quality & Reliability (AREA)
  • Position Fixing By Use Of Radio Waves (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

本发明公开了一种考虑非线性因素的区域格网速度场构建方法及系统，涉及大地测量领域，该方法包括：对目标区域的GNSS坐标时间序列，采用稳健最小二乘法进行线性拟合；采用改进后的周期图谱法对拟合后的第一残余时间序列进行周期项提取；改进后的周期图谱算法为引入时间平移不变量，并以加权迭代的方式实现的周期图谱法；基于周期项信息和阶跃项信息，采用稳健最小二乘法对GNSS坐标时间序列进行线性拟合，得到参数解；采用随机模型法对参数解进行精度评估得到中误差；根据中误差从参数解中选取目标速度项拟合参数；根据目标速度项拟合参数，采用反距离加权法计算目标区域内各格网点的速度值，得到格网速度场。本发明能提高速度场的精度。

Description

一种考虑非线性因素的区域格网速度场构建方法及系统
技术领域
本发明涉及大地测量技术领域，特别是涉及一种考虑非线性因素的区域格网速度场构建方法及系统。
背景技术
GNSS坐标时间序列是一组按照时间顺序排列的基准站坐标组合，其包含了丰富的信息，不仅可以反映测站的线性变化运动，也可以反映出测站存在的非线性变化。线性变化主要表现为速度信号，其反映的是测站受同一区域的构造应力场控制的构造运动，而非线性变化主要表现为周期性信号，其反映的是测站受到的非潮汐海洋负载、大气负载、水文负载、冰期后回弹以及区域性共模误差等地球物理效应的作用。此外，坐标序列中还存在由地壳运动、仪器更换等因素引起的阶跃或震后变形信号。因此，GNSS坐标时间序列的分析与建模，特别是研究非线性信号的变化特征，可以更精确地分离测站速度信息，有助于合理解释板块构造运动、建立和维持动态地球参考框架，而且还能构建更高精度的区域格网速度场模型，具有重要的理论研究意义和实际应用价值。正因如此，GNSS坐标时间序列分析理论与应用研究成为当前大地测量学和地球物理学等领域的研究热点。
对于时间序列中非线性信号的估计，目前常用方法是直接考虑时间序列中的周年项和半年周项周期信号，与速度参数一起利用最小二乘理论进行估计，这一方法的明显不足之处是周期信号考虑不全，且每个测站包含的周期信号特征也有一定的差异。因此，为了建立更加准确的基准站非线性运动模型，有必要对序列中的非线性信号进行分析。而周期图谱方法，正是一种是有效适用于非均匀时间序列的周期分析方法。目前，该算法已广泛应用于天文、经济、地球物理和生物医学等学科领域的非均匀实验观测数据的频谱分析。但是，由于受到序列的非均匀性、有限长度等因素，该算法在傅里叶变换的功率谱中会产生虚假谱峰，此外，由于噪声的影响，周期信号的振幅和相位也可能存在一定的误差。可能正因如此，周期图谱法在GNSS坐标时间序列分析中仅仅作为一种辅助手段，并未得到充分和广泛的应用。基于以上考虑，GLS算法被提出，该算法在一定程度上弥补了周期图谱法存在的缺点。但是，上述提及的速度场求解方法仍然需要预先对带有缺失值的非均匀GNSS原始坐标时间序列进行插补处理才能进行后续分析，因此速度场的精度有待进一步提高。
发明内容
基于此，本发明实施例提供一种考虑非线性因素的区域格网速度场构建方法及系统，无需对带有缺失值的非均匀GNSS原始坐标时间序列进行插补处理，以提高速度场的精度。
为实现上述目的，本发明提供了如下方案：
一种考虑非线性因素的区域格网速度场构建方法，包括：
获取大地测量过程中测站测得的目标区域的GNSS坐标时间序列；
采用稳健最小二乘法对所述GNSS坐标时间序列进行线性拟合，得到第一残余时间序列；所述第一残余时间序列为去除线性速度项的GNSS坐标时间序列；
采用改进后的周期图谱法对所述第一残余时间序列进行周期项提取，得到周期项信息；所述周期项信息包括周期个数和周期长度；所述改进后的周期图谱算法为引入时间平移不变量，并以加权迭代的方式实现的周期图谱法；
基于所述周期项信息和阶跃项信息，采用稳健最小二乘法对所述GNSS坐标时间序列进行线性拟合，得到所述GNSS坐标时间序列的参数解；所述阶跃项信息为阶跃发生的历元时刻；所述参数解包括阶跃项拟合参数、速度项拟合参数和周期项拟合参数；
采用随机模型法对所述参数解进行精度评估，得到中误差；所述中误差包括第二残余时间序列的中误差值和所述参数解的中误差值；所述第二残余时间序列为去除阶跃项、速度项和周期项的GNSS坐标时间序列；
根据所述中误差从所述参数解中选取目标速度项拟合参数；所述目标速度项拟合参数为所述参数解中精度高于设定精度的速度项拟合参数；
根据所述目标速度项拟合参数，采用反距离加权法计算所述目标区域内各格网点的速度值，得到所述目标区域的格网速度场。
可选的，所述采用稳健最小二乘法对所述GNSS坐标时间序列进行线性拟合，得到第一残余时间序列，具体包括：
采用滑动中位数法对所述GNSS坐标时间序列进行粗差剔除；
采用稳健最小二乘法对粗差剔除后的GNSS坐标时间序列进行线性拟合，得到第一残余时间序列。
可选的，所述基于所述周期项信息和阶跃项信息，采用稳健最小二乘法对所述GNSS坐标时间序列进行线性拟合，得到所述GNSS坐标时间序列的参数解，具体包括：
采用滑动中位数法对所述GNSS坐标时间序列进行粗差剔除和阶跃探测，得到粗差剔除后的GNSS坐标时间序列和阶跃项信息；
基于所述周期项信息和所述阶跃项信息，采用稳健最小二乘法对所述粗差剔除后的GNSS坐标时间序列进行线性拟合，得到所述GNSS坐标时间序列的参数解。
可选的，所述采用改进后的周期图谱法对所述第一残余时间序列进行周期项提取，得到周期项信息，具体包括：
在周期图谱法的拟合模型中增加时间平移不变量，得到改变后的拟合模型；
采用所述改变后的拟合模型对所述第一残余时间序列进行分析，构建振幅因子矩阵；
在当前迭代次数下，基于当前迭代次数下的第一残余时间序列的观测误差和所述振幅因子矩阵，构建当前迭代次数下权阵因子；
根据当前迭代次数下权阵因子计算当前迭代次数下的功率谱，并由当前迭代次数下的功率谱提取当前迭代次数下的主周期和对应的振幅；
判断当前迭代次数是否达到设定最大迭代次数；所述设定最大迭代次数等于根据所述GNSS坐标时间序列设定的给定周期个数；
若是，则根据所述迭代次数下的主周期和对应的振幅确定周期项信息；
若否，则进行下一次迭代。
本发明还提供了一种考虑非线性因素的区域格网速度场构建系统，包括：
时间序列获取模块，用于获取测站测得的目标区域的GNSS坐标时间序列；
第一拟合模块，用于采用稳健最小二乘法对所述GNSS坐标时间序列进行线性拟合，得到第一残余时间序列；所述第一残余时间序列为去除线性速度项的GNSS坐标时间序列；
周期项提取模块，用于采用改进后的周期图谱法对所述第一残余时间序列进行周期项提取，得到周期项信息；所述周期项信息包括周期个数和周期长度；所述改进后的周期图谱算法为引入时间平移不变量，并以加权迭代的方式实现的周期图谱法；
第二拟合模块，用于基于所述周期项信息和阶跃项信息，采用稳健最小二乘法对所述GNSS坐标时间序列进行线性拟合，得到所述GNSS坐标时间序列的参数解；所述阶跃项信息为阶跃发生的历元时刻；所述参数解包括阶跃项拟合参数、速度项拟合参数和周期项拟合参数；
精度评估模块，用于采用随机模型法对所述参数解进行精度评估，得到中误差；所述中误差包括第二残余时间序列的中误差值和所述参数解的中误差值；所述第二残余时间序列为去除阶跃项、速度项和周期项的GNSS坐标时间序列；
速度项拟合参数选取模块，用于根据所述中误差从所述参数解中选取目标速度项拟合参数；所述目标速度项拟合参数为所述参数解中精度高于设定精度的速度项拟合参数；
速度场构建模块，用于根据所述目标速度项拟合参数，采用反距离加权法计算所述目标区域内各格网点的速度值，得到所述目标区域的格网速度场。
可选的，所述第一拟合模块，具体包括：
粗差剔除单元，用于采用滑动中位数法对所述GNSS坐标时间序列进行粗差剔除；
第一拟合单元，用于采用稳健最小二乘法对粗差剔除后的GNSS坐标时间序列进行线性拟合，得到第一残余时间序列。
可选的，所述第二拟合模块，具体包括：
预处理单元，用于采用滑动中位数法对所述GNSS坐标时间序列进行粗差剔除和阶跃探测，得到粗差剔除后的GNSS坐标时间序列和阶跃项信息；
第二拟合单元，用于基于所述周期项信息和所述阶跃项信息，采用稳健最小二乘法对所述粗差剔除后的GNSS坐标时间序列进行线性拟合，得到所述GNSS坐标时间序列的参数解。
可选的，所述周期项提取模块，具体包括：
拟合模型更新单元，用于在周期图谱法的拟合模型中增加时间平移不变量，得到改变后的拟合模型；
振幅因子矩阵构建单元，用于采用所述改变后的拟合模型对所述第一残余时间序列进行分析，构建振幅因子矩阵；
权阵因子构建单元，用于在当前迭代次数下，基于当前迭代次数下的第一残余时间序列的观测误差和所述振幅因子矩阵，构建当前迭代次数下权阵因子；
周期信息提取单元，用于根据当前迭代次数下权阵因子计算当前迭代次数下的功率谱，并由当前迭代次数下的功率谱提取当前迭代次数下的主周期和对应的振幅；
迭代更新单元，用于判断当前迭代次数是否达到设定最大迭代次数；所述设定最大迭代次数等于根据所述GNSS坐标时间序列设定的给定周期个数；若是，则根据所述迭代次数下的主周期和对应的振幅确定周期项信息；若否，则进行下一次迭代。
与现有技术相比，本发明的有益效果是：
本发明实施例提出了一种考虑非线性因素的区域格网速度场构建方法及系统，对目标区域的GNSS坐标时间序列，采用稳健最小二乘法进行线性拟合；采用改进后的周期图谱法对拟合后的第一残余时间序列进行周期项提取，得到周期项信息；改进后的周期图谱算法为引入时间平移不变量，并以加权迭代的方式实现的周期图谱法；基于周期项信息和阶跃项信息，采用稳健最小二乘法对GNSS坐标时间序列进行线性拟合，得到参数解；采用随机模型法对参数解进行精度评估，得到中误差；根据中误差从参数解中选取目标速度项拟合参数；根据目标速度项拟合参数，采用反距离加权法计算目标区域内各格网点的速度值，得到格网速度场。本发明无需预先对带有缺失值的非均匀GNSS原始坐标时间序列进行插补处理即可进行后续分析，同时考虑了序列中的非线性变化因素（同时考虑了周期项、阶跃项和速度项），进一步提高了速度解的精度，从而为区域格网速度场自动一体化构建方法提供一套完善可行的依据。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
图1为本发明实施例提供的考虑非线性因素的区域格网速度场构建方法的流程图；
图2为本发明实施例提供的考虑非线性因素的区域格网速度场构建方法的实现框图；
图3为本发明实施例提供的考虑非线性因素的区域格网速度场构建系统的结构图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。
图1为本发明实施例提供的考虑非线性因素的区域格网速度场构建方法的流程图。
参见图1，本实施例的考虑非线性因素的区域格网速度场构建方法，包括：
步骤101：获取测站测得的目标区域的GNSS坐标时间序列。
步骤102：采用稳健最小二乘法对所述GNSS坐标时间序列进行线性拟合，得到第一残余时间序列；所述第一残余时间序列为去除线性速度项的GNSS坐标时间序列。
步骤103：采用改进后的周期图谱法对所述第一残余时间序列进行周期项提取，得到周期项信息；所述周期项信息包括周期个数和周期长度；所述改进后的周期图谱算法为引入时间平移不变量，并以加权迭代的方式实现的周期图谱法。
步骤104：基于所述周期项信息和阶跃项信息，采用稳健最小二乘法对所述GNSS坐标时间序列进行线性拟合，得到所述GNSS坐标时间序列的参数解。所述阶跃项信息为阶跃发生的历元时刻；所述参数解包括阶跃项拟合参数、速度项拟合参数和周期项拟合参数。
步骤105：采用随机模型法对所述参数解进行精度评估，得到中误差；所述中误差包括第二残余时间序列的中误差值和所述参数解的中误差值；所述第二残余时间序列为去除阶跃项、速度项和周期项的GNSS坐标时间序列。
步骤106：根据所述中误差从所述参数解中选取目标速度项拟合参数；所述目标速度项拟合参数为所述参数解中精度高于设定精度的速度项拟合参数。
步骤107：根据所述目标速度项拟合参数，采用反距离加权法计算所述目标区域内各格网点的速度值，得到所述目标区域的格网速度场。
作为一种可选的实施方式，所述步骤102，具体包括：
1）采用滑动中位数法对所述GNSS坐标时间序列进行粗差剔除。
具体的，采用基于IQR准则的滑动中位数粗差探测法进行粗差剔除。该方法依据的合理假设条件为：短时间（数天或数周）内GNSS测站的坐标处于稳定状态，通常情况下GNSS静态观测站的位置是稳定不变的。正是由于中位数不受序列中离群值的影响，因此可以在设定合理滑动窗口的条件下，可以有效探测序列中的离群值。其主要实现步骤为：
S21、设定合适的滑动窗口长度w，以序列位置i为中心前后分别提取依据w/2长度的子序列，称为sub1和sub2；S22、在子序列内依据IQR准则判断i点是否为离群值，若为离群值则将其视为粗差进行剔除；S23、重复步骤S102，逐步滑动窗口至整个序列，判断并剔除所有的离群值；S24、第一次粗差探测完毕并剔除粗差后，继续进行第二次粗差探测，迭代至序列不在含有粗差为止。
2）采用稳健最小二乘法对粗差剔除后的GNSS坐标时间序列进行线性拟合，得到第一残余时间序列。具体的：
采用稳健最小二乘模型进行线性拟合并提取残余序列的。利用稳健最小二乘进行线性拟合的原理如下：
仅考虑粗差剔除后的GNSS坐标时间序列中包含的线性信号，观测方程为：
Figure 450344DEST_PATH_IMAGE001
(1)
式中，
Figure 101905DEST_PATH_IMAGE002
为
Figure 412801DEST_PATH_IMAGE003
时刻的测站分量坐标值，即粗差剔除后的GNSS坐标时间序列，
Figure 339169DEST_PATH_IMAGE004
为常量项，
Figure 950279DEST_PATH_IMAGE005
线性速率系数，
Figure 303900DEST_PATH_IMAGE003
为第
Figure 305354DEST_PATH_IMAGE006
个时刻，
Figure 533947DEST_PATH_IMAGE007
为测量误差。
改写为误差方程：
Figure 999564DEST_PATH_IMAGE008
(2)
式中，
Figure 727348DEST_PATH_IMAGE009
为残差值，
Figure 278416DEST_PATH_IMAGE010
为拟合估值，
Figure 749848DEST_PATH_IMAGE002
为观测值，
Figure 69971DEST_PATH_IMAGE011
为观测历元个数。
误差方程的矩阵形式为：
Figure 234236DEST_PATH_IMAGE012
(3)
式中，
Figure 7020DEST_PATH_IMAGE013
,
Figure 813302DEST_PATH_IMAGE014
,
Figure 987931DEST_PATH_IMAGE015
,
Figure 854256DEST_PATH_IMAGE016
。
根据公式（4）准则，采用最小二乘拟合估计，以得到第一次未知数的解
Figure 881380DEST_PATH_IMAGE017
。
Figure 694616DEST_PATH_IMAGE018
(4)
Figure 723751DEST_PATH_IMAGE019
(5)
式中，
Figure 292136DEST_PATH_IMAGE020
，
Figure 977195DEST_PATH_IMAGE021
表示第i个观测值的观测误差。
根据残差
Figure 656438DEST_PATH_IMAGE009
和定权公式（6）重新确定各观测值新的权因子，进行下一轮计算，直至前后两次估计值的变化小于限差时结束迭代。迭代完成后，将最后一次的参数结果
Figure 540081DEST_PATH_IMAGE017
代回公式（3），得到该步骤中第一残余时间序列
Figure 217050DEST_PATH_IMAGE022
。
Figure 451722DEST_PATH_IMAGE023
(6)
式中，
Figure 606760DEST_PATH_IMAGE024
为权值，
Figure 876067DEST_PATH_IMAGE025
为阈值常量，一般取二倍或三倍的单位权中误差。
作为一种可选的实施方式，步骤103是采用改进的周期图谱法提取周期项先验信息的，其主要原理如下：
周期图谱法（Lomb-Scargle periodogram，LSP）是一种基于离散傅里叶变换的周期提取方法。该算法在一定程度可以解决非均匀采样间隔对周期信号的影响，避免对非均匀采样时间序列进行插补处理，且考虑了非均匀采样对幅度和相位对信号产生的影响。该算法的基本原理是通过最小二乘法将一系列三角函数的进行线性组合来拟合时间序列，在此基础上，将信号特征从时域转换到频域上。
对于非均匀时间序列
Figure 723937DEST_PATH_IMAGE026
，即步骤102得到的第一残余时间序列结果，该第一残余时间序列的平均值为0，定义拟合方程为：
Figure 180327DEST_PATH_IMAGE027
(7)
式中，
Figure 723345DEST_PATH_IMAGE003
为离散序列的采样时间，
Figure 519262DEST_PATH_IMAGE011
为时间序列的数据统计量。离散测试频率
Figure 600351DEST_PATH_IMAGE028
，
Figure 747299DEST_PATH_IMAGE029
，
Figure 306456DEST_PATH_IMAGE030
定义为序列的极限频率，
Figure 284776DEST_PATH_IMAGE031
一般不大于尼奎斯特频率
Figure 208870DEST_PATH_IMAGE032
，测试频率的取样步长
Figure 905430DEST_PATH_IMAGE033
，测试频率数量
Figure 205962DEST_PATH_IMAGE034
。
Figure 304368DEST_PATH_IMAGE035
表示频率分量f k 的正弦变化的幅度，
Figure 399363DEST_PATH_IMAGE036
表示频率分量f k 的余弦变化的幅度。
Figure 583219DEST_PATH_IMAGE037
为第i个观测值的观测噪声。
为简便后续公式的推导，定义以下变量（仅作为临时变量辅助公式推导）：
Figure 484179DEST_PATH_IMAGE038
(8)
根据间接平差原理，构建误差方程如下：
Figure 109196DEST_PATH_IMAGE039
(9)
式中，
Figure 204453DEST_PATH_IMAGE040
为残差向量，
Figure 547710DEST_PATH_IMAGE041
为系数矩阵，
Figure 517940DEST_PATH_IMAGE042
为待估参数，
Figure 263042DEST_PATH_IMAGE043
为观测向量。
基于最小二乘原理，得到最佳估值：
Figure 27735DEST_PATH_IMAGE044
(10)
其中，
Figure 655026DEST_PATH_IMAGE045
。
则LSP的功率谱定义为：
Figure 101051DEST_PATH_IMAGE046
(11)
式中，
Figure 762976DEST_PATH_IMAGE047
为
Figure 636254DEST_PATH_IMAGE048
频率的功率谱值。
为了避免序列取样时间的平移对谱结构产生影响，引入时间平移不变量
Figure 16420DEST_PATH_IMAGE049
，其他参数变量含义同公式（7），建立新的拟合模型：
Figure 556DEST_PATH_IMAGE050
(12)
加入
Figure 782568DEST_PATH_IMAGE049
重新定义以下变量：
Figure 826747DEST_PATH_IMAGE051
(13)
选择
Figure 694209DEST_PATH_IMAGE049
值使得
Figure 482036DEST_PATH_IMAGE052
成立，推导可得：
Figure 351510DEST_PATH_IMAGE053
(14)
基于最小二乘原理，得到最佳估值
Figure 628908DEST_PATH_IMAGE017
：
Figure 655769DEST_PATH_IMAGE054
(15)
则LSP的功率谱定义为：
Figure 309605DEST_PATH_IMAGE055
(16)
使用显著性等级
Figure 738312DEST_PATH_IMAGE056
用以评价功率谱的质量，
Figure 921032DEST_PATH_IMAGE056
代表频谱的虚警概率。当
Figure 763086DEST_PATH_IMAGE056
=0.1时，频谱的置信度则为90%。
Figure 892716DEST_PATH_IMAGE057
(17)
其中，
Figure 238246DEST_PATH_IMAGE058
为估计的频谱值，
Figure 60709DEST_PATH_IMAGE011
为频谱中包含的独立周期的个数。
传统的LSP算法虽然有效，但仍然存在以下缺陷：①没有考虑观测误差对结果造成的影响；②预先假设原始序列与拟合所采用的正弦函数的均值相同；③非均匀序列往往会在真实信号两侧出现虚假谱峰；④没有考虑多重信号之间相互调制对结果造成的影响。
第一个缺陷易导致周期成分的振幅和相位偏离真实值；第二个缺陷易导致周期成分的振幅产生系统误差。针对这两个问题，Zechmeister提出了GLSP算法，采用了修正的正弦函数来对时间序列进行拟合。GLSP算法与LSP算法相似，主要区别在于GLSP算法通过加入常数的正弦三角函数对时间序列进行拟合分析，且计算时进一步考虑了观测误差的时间序列的影响。因此，从理论上说，GLSP算法的周期信号估计准确度应优于LSP算法，但GLSP算法依旧没有考虑后两个缺陷所产生的影响，为此，本发明将提出一种新的ILSP算法，主要采用加权和迭代两种方法以有效减小各频率信号之间的相互调制和观测噪声对功率谱结果产生的影响。考虑到其他频率分量和观测误差或噪声对该频率信号的影响，定义新的加权因子
Figure 124480DEST_PATH_IMAGE059
为：
Figure 57801DEST_PATH_IMAGE060
(18)
其中，
Figure 257838DEST_PATH_IMAGE061
,
Figure 549404DEST_PATH_IMAGE062
，
Figure 38154DEST_PATH_IMAGE063
表示频率分量f m 的正弦变化的幅度，
Figure 571904DEST_PATH_IMAGE064
表示频率分量f m 的余弦变化的幅度。
那么得到加权后的参数估值和功率谱：
Figure 829710DEST_PATH_IMAGE065
(19)
Figure 56292DEST_PATH_IMAGE066
(20)
式中，
Figure 563497DEST_PATH_IMAGE067
为参数估值，
Figure 838620DEST_PATH_IMAGE068
为系数矩阵，
Figure 278829DEST_PATH_IMAGE069
为加权矩阵，
Figure 348416DEST_PATH_IMAGE070
为观测向量，
Figure 874075DEST_PATH_IMAGE071
为
Figure 687310DEST_PATH_IMAGE072
频率的功率谱值。
根据频率
Figure 982026DEST_PATH_IMAGE072
与功率谱
Figure 222514DEST_PATH_IMAGE071
的一一对应关系，寻找出功率谱最大值对应的频率
Figure 969890DEST_PATH_IMAGE073
，进而根据公式(22)计算出对应的周期
Figure 147668DEST_PATH_IMAGE074
。
Figure 234573DEST_PATH_IMAGE073
=
Figure 973859DEST_PATH_IMAGE075
max[
Figure 942952DEST_PATH_IMAGE076
(21)
式中，
Figure 97990DEST_PATH_IMAGE077
n为周期个数。
Figure 367297DEST_PATH_IMAGE078
(22)
基于虚假谱峰值低于主峰值这一特性，采取迭代法消除虚假谱峰的影响，即第
Figure 949588DEST_PATH_IMAGE079
次迭代仅提取序列中主峰值对应的周期和振幅，然后从序列中减去对应的周期信号以得到残差序列重复计算，直至达到设定的迭代次数（即目标周期个数）结束迭代。每一次迭代都得到一个周期值
Figure 671557DEST_PATH_IMAGE074
（周期长度），n次迭代完成后，得到n个周期
Figure 630285DEST_PATH_IMAGE074
（周期个数），即得到最终的周期项信息。
基于上述原理，所述步骤103，具体包括：
S31：在周期图谱法的拟合模型中增加时间平移不变量，得到改变后的拟合模型。
S32：采用所述改变后的拟合模型对所述第一残余时间序列进行分析，构建振幅因子矩阵。
S33：在当前迭代次数下，基于当前迭代次数下的第一残余时间序列的观测误差和所述振幅因子矩阵，构建当前迭代次数下权阵因子。
S34：根据当前迭代次数下权阵因子计算当前迭代次数下的功率谱，并由当前迭代次数下的功率谱提取当前迭代次数下的主周期和对应的振幅。
S35：判断当前迭代次数是否达到设定最大迭代次数；所述设定最大迭代次数等于根据所述GNSS坐标时间序列设定的给定周期个数。若是，则根据所述迭代次数下的主周期和对应的振幅确定周期项信息；若否，则进行下一次迭代。每一次迭代都得到一个周期值（周期长度），多次迭代完成后，得到相应个数的周期（周期个数），这样得到最终的周期项信息。
步骤103，一个更为具体的实现过程如下：
改进周期图谱法需要经过迭代提取周期信息。结合加权与迭代法，改进周期图谱法的具体实现步骤如下：
①给定原始时间序列和待探测的周期个数N；②进行时间序列自动化预处理（主要包括粗差剔除与阶跃修复），基于稳健最小二乘模型估计并消除序列中的线性趋势项；③进行LSP分析，以得到频率对应的振幅，构造振幅因子矩阵。若N为0，进行功率谱的显著性评估，取
Figure 754099DEST_PATH_IMAGE080
=0.1，获取置信度大于99%的谱峰个数，并设为N。④进行ILSP分析，基于序列观测误差和振幅因子，构造权阵因子
Figure 772871DEST_PATH_IMAGE081
，计算功率谱，提取主周期和对应的振幅。⑤去除序列中的主周期项，以得到残差序列。⑥重复④和⑤两个步骤N次结束迭代。
作为一种可选的实施方式，所述步骤104，具体包括：
1）采用滑动中位数法对所述GNSS坐标时间序列进行粗差剔除和阶跃探测，得到粗差剔除后的GNSS坐标时间序列和阶跃项信息。
粗差剔除过程是采用基于IQR准则的滑动中位数粗差探测法实现的，其主要实现步骤与步骤102中粗差剔除的过程相同，在此不再赘述。
阶跃探测（阶跃项的定位）过程是基于K-Medoids的滑动中位数阶跃定位法进行的，该方法依据的合理假设条件为：短时间（数天或数周）内GNSS测站的坐标处于稳定状态，通常情况下GNSS静态观测站的位置是不变的。该方法同时对GNSS测站三个方向的坐标时间序列进行阶跃定位，考虑了序列间的相关性，使用滑动中位数法可以避免离群值的影响，而且使用K-Medoids对阶跃发生历元进行聚类进一步保证了结果的可靠性，因此是一种较为简便有效的阶跃定位策略。阶跃探测的主要实现步骤为：
步骤S41、设定合适的滑动窗口长度w、阶跃探测阈值e，以序列位置i为中心前后分别提取w/2长度的子序列sub1和sub2，计算sub1和sub2的中位数med1和med2的绝对值之差，若其绝对值大于探测阈值e，则判断位置i处为阶跃发生历元。
步骤S42、逐步滑动窗口至整个序列，完成GNSS坐标三个方向（NEU）的序列的阶跃探测，然后合并三个方向的阶跃发生历元并去除重复值。
步骤S43、确定阶跃发生历元的种类个数，对历元结果进行K-Medoids聚类分析，求取每一类阶跃发生历元的中位数，将其作为准确的阶跃发生历元。
步骤S44、确定阶跃发生历元后，以NEU序列阶跃发生位置为中心前后分别提取w/2长度的子序列sub1和sub2，计算sub1和sub2的中位数med1和med2的差值作为阶跃值。
在步骤S43中，是采用K-Medoids聚类算法对阶跃发生历元进行聚类分析的。K-Medoids聚类算法是一种分区方法，其类似于K-Means，两种方法的目标都是将一组测量值或观察值划分为k个子集或集群，以便子集最小化测量值与测量值集群中心之间的距离总和。在K-Means聚类算法中，子集的中心是子集中测量的平均值，而K-Medoids聚类算法中，子集的中心是子集中测量的中位数。相较于平均值，中位值不易受到集合中离群值的影响，因此K-Medoids聚类算法通常用于对异常数据、任意距离度量、平均值或中位数等没有明确定义的数据具有鲁棒性要求的领域。
2）基于所述周期项信息和所述阶跃项信息，采用稳健最小二乘法对所述粗差剔除后的GNSS坐标时间序列进行线性拟合，得到所述GNSS坐标时间序列的参数解。
采用稳健最小二乘法进行参数解拟合的，其主要原理如下：
在原始GNSS坐标时间序列中，通常需要考虑测站趋势项（线性速度）、周期变化项（主要考虑年周期和半年周期）、非地震因素（设备更换、天线高测量误差、相位中心建模误差或其他人为以及软件误差）或地震因素（地震同震断裂）引起的阶跃跳变项、震后形变项（通常呈指数或对数变化形式），以及一些未模型化的误差项。暂不考虑未模型化的误差项，且由于历元之间坐标分量的相关性很小，
Figure 716556DEST_PATH_IMAGE003
历元时刻的坐标分量可以详细模型化表示为：
Figure 275713DEST_PATH_IMAGE082
(23)
式中，
Figure 191717DEST_PATH_IMAGE083
为初始时刻
Figure 945171DEST_PATH_IMAGE084
对应的截距。
Figure 579415DEST_PATH_IMAGE085
为从
Figure 942263DEST_PATH_IMAGE084
为参考的历元时刻，单位是年。线性速率
Figure 712773DEST_PATH_IMAGE086
为震间测站的长期的构造运动。c、d、e、f为年周期和半年周期项系数，估计周期项参数则需要足够多的数据（至少12个月），年周期和半年周期的信号用
Figure 135664DEST_PATH_IMAGE087
表示，
Figure 53942DEST_PATH_IMAGE088
为幅度，
Figure 892585DEST_PATH_IMAGE089
为角速率，
Figure 845497DEST_PATH_IMAGE090
为一月一日，
Figure 376973DEST_PATH_IMAGE091
为相位。
Figure 782546DEST_PATH_IMAGE092
为发生阶跃的总次数，
Figure 690459DEST_PATH_IMAGE093
为
Figure 232299DEST_PATH_IMAGE094
历元时刻由于同震或非同震变化引起的阶跃量。
Figure 996993DEST_PATH_IMAGE095
为阶梯函数。
Figure 827545DEST_PATH_IMAGE096
为震后形变函数。
Figure 568843DEST_PATH_IMAGE097
为观测噪声。
Figure 434031DEST_PATH_IMAGE095
函数表示为：
Figure 369626DEST_PATH_IMAGE098
(24)
式中，
Figure 687475DEST_PATH_IMAGE099
为发生阶跃的时刻。
Figure 468349DEST_PATH_IMAGE096
函数通常可由指数或对数函数表示：
Figure 250360DEST_PATH_IMAGE100
(25)
式中，
Figure 294540DEST_PATH_IMAGE101
为
Figure 162001DEST_PATH_IMAGE102
地震发生时刻的阶跃系数，
Figure 949829DEST_PATH_IMAGE103
为地震衰减因子。
GNSS坐标时间序列模型中各事件的发生时刻可从地震目录、观测日志、自动检测算法或目视检查确定。由于自动检测算法可靠性差，目视检查法实施起来难度较大，本发明中主要根据地震目录和测站日志来确定测站坐标的阶跃发生时刻和震后形变发生时刻。此外，地震衰减因子通常需要通过其他方法进行单独估计。因此，剩余时间序列系数可以表示为线性模型，然后再次基于稳健最小二乘模型(公式（1）-公式（6）)可得到第二残余时间序列和参数的最佳拟合估值
Figure 586347DEST_PATH_IMAGE104
该参数估值中，
Figure 598165DEST_PATH_IMAGE105
为速度解（速度项拟合参数），[
Figure 625027DEST_PATH_IMAGE106
为周期振幅（周期项拟合参数），
Figure 278862DEST_PATH_IMAGE107
为阶跃值（周期项拟合参数），
Figure 707569DEST_PATH_IMAGE108
为震后形变系数。
其中，步骤105，其主要原理如下：
时间序列参数
Figure 657333DEST_PATH_IMAGE017
的中误差估值
Figure 437070DEST_PATH_IMAGE109
为：
Figure 629017DEST_PATH_IMAGE110
(26)
式中，
Figure 708969DEST_PATH_IMAGE111
为
Figure 531431DEST_PATH_IMAGE112
的主对角元素，
Figure 595202DEST_PATH_IMAGE113
，
Figure 262944DEST_PATH_IMAGE114
为单位权中误差估值，即第二残余时间序列的中误差，有：
Figure 728560DEST_PATH_IMAGE115
(27)
式中，
Figure 721924DEST_PATH_IMAGE116
为多余观测量，
Figure 741833DEST_PATH_IMAGE117
为观测值改正数，即步骤104得到的第二残余序列（已去除阶跃项、速度项和周期项）。参数解的中误差值的计算方式与第二残余时间序列的中误差的计算方式类似，在此不再赘述。之后，根据计算得到的中误差从参数解中选取精度高于设定精度的速度项拟合参数，即得到目标速度项拟合参数。
其中，步骤107，根据上述步骤得到符合精度要求的速度值（目标速度项拟合参数），基于反距离加权法构建区域格网速度场，其主要原理如下：
设空间待插值格网点为
Figure 541161DEST_PATH_IMAGE118
，P点邻域内有n个已知离散点
Figure 798967DEST_PATH_IMAGE119
。根据离散点的速度值(
Figure 759970DEST_PATH_IMAGE120
)，通过反距离加权法对P点的速度值
Figure 736016DEST_PATH_IMAGE121
进行插值，即：
Figure 306413DEST_PATH_IMAGE122
(28)
式中，
Figure 481042DEST_PATH_IMAGE123
，表示离散点到待定点的距离；
Figure 816208DEST_PATH_IMAGE124
一般取1~2，本文取
Figure 76289DEST_PATH_IMAGE124
=1。最终步骤得到格网速度场，即每个格网点
Figure 155103DEST_PATH_IMAGE125
的速度值
Figure 449818DEST_PATH_IMAGE126
。
上述实施例中的区域格网速度场构建方法实现框图如图2所示，由于方法中采用的周期图谱法在一定程度上解决了非均匀采样间隔产生的非周期信号问题，从而避免了对非均匀采样时间序列进行插值处理，同时考虑了非均匀采样对幅度和相位带来的影响，而且改进后的周期图谱法主要采用加权和迭代两种方法进一步考虑了序列观测噪声误差与虚假谱峰的影响，以有效减小各频率信号之间的相互调制和观测噪声对功率谱结果产生的影响。因此，此方法无需预先对带有缺失值的非均匀GNSS原始坐标时间序列进行插补处理即可进行后续分析，同时考虑了序列中的非线性变化因素，进一步提高了速度解的精度，从而为区域格网速度场自动一体化构建方法提供一套完善可行的算法依据。
本发明还提供了一种考虑非线性因素的区域格网速度场构建系统，参见图3，所述系统，包括：
时间序列获取模块301，用于获取测站测得的目标区域的GNSS坐标时间序列。
第一拟合模块302，用于采用稳健最小二乘法对所述GNSS坐标时间序列进行线性拟合，得到第一残余时间序列；所述第一残余时间序列为去除线性速度项的GNSS坐标时间序列。
周期项提取模块303，用于采用改进后的周期图谱法对所述第一残余时间序列进行周期项提取，得到周期项信息；所述周期项信息包括周期个数和周期长度；所述改进后的周期图谱算法为引入时间平移不变量，并以加权迭代的方式实现的周期图谱法。
第二拟合模块304，用于基于所述周期项信息和阶跃项信息，采用稳健最小二乘法对所述GNSS坐标时间序列进行线性拟合，得到所述GNSS坐标时间序列的参数解；所述阶跃项信息为阶跃发生的历元时刻；所述参数解包括阶跃项拟合参数、速度项拟合参数和周期项拟合参数。
精度评估模块305，用于采用随机模型法对所述参数解进行精度评估，得到中误差；所述中误差包括第二残余时间序列的中误差值和所述参数解的中误差值；所述第二残余时间序列为去除阶跃项、速度项和周期项的GNSS坐标时间序列。
速度项拟合参数选取模块306，用于根据所述中误差从所述参数解中选取目标速度项拟合参数；所述目标速度项拟合参数为所述参数解中精度高于设定精度的速度项拟合参数。
速度场构建模块307，用于根据所述目标速度项拟合参数，采用反距离加权法计算所述目标区域内各格网点的速度值，得到所述目标区域的格网速度场。
作为一种可选的实施方式，所述第一拟合模块302，具体包括：
粗差剔除单元，用于采用滑动中位数法对所述GNSS坐标时间序列进行粗差剔除。
第一拟合单元，用于采用稳健最小二乘法对粗差剔除后的GNSS坐标时间序列进行线性拟合，得到第一残余时间序列。
作为一种可选的实施方式，所述第二拟合模块304，具体包括：
预处理单元，用于采用滑动中位数法对所述GNSS坐标时间序列进行粗差剔除和阶跃探测，得到粗差剔除后的GNSS坐标时间序列和阶跃项信息。
第二拟合单元，用于基于所述周期项信息和所述阶跃项信息，采用稳健最小二乘法对所述粗差剔除后的GNSS坐标时间序列进行线性拟合，得到所述GNSS坐标时间序列的参数解。
作为一种可选的实施方式，所述周期项提取模块303，具体包括：
拟合模型更新单元，用于在周期图谱法的拟合模型中增加时间平移不变量，得到改变后的拟合模型。
振幅因子矩阵构建单元，用于采用所述改变后的拟合模型对所述第一残余时间序列进行分析，构建振幅因子矩阵。
权阵因子构建单元，用于在当前迭代次数下，基于当前迭代次数下的第一残余时间序列的观测误差和所述振幅因子矩阵，构建当前迭代次数下权阵因子。
周期信息提取单元，用于根据当前迭代次数下权阵因子计算当前迭代次数下的功率谱，并由当前迭代次数下的功率谱提取当前迭代次数下的主周期和对应的振幅。
迭代更新单元，用于判断当前迭代次数是否达到设定最大迭代次数；所述设定最大迭代次数等于根据所述GNSS坐标时间序列设定的给定周期个数；若是，则根据所述迭代次数下的主周期和对应的振幅确定周期项信息；若否，则进行下一次迭代。
本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的系统而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。
本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。

Claims (8)

1.一种考虑非线性因素的区域格网速度场构建方法，其特征在于，包括：
获取测站测得的目标区域的GNSS坐标时间序列；
采用稳健最小二乘法对所述GNSS坐标时间序列进行线性拟合，得到第一残余时间序列；所述第一残余时间序列为去除线性速度项的GNSS坐标时间序列；
采用改进后的周期图谱法对所述第一残余时间序列进行周期项提取，得到周期项信息；所述周期项信息包括周期个数和周期长度；所述改进后的周期图谱算法为引入时间平移不变量，并以加权迭代的方式实现的周期图谱法；
基于所述周期项信息和阶跃项信息，采用稳健最小二乘法对所述GNSS坐标时间序列进行线性拟合，得到所述GNSS坐标时间序列的参数解；所述阶跃项信息为阶跃发生的历元时刻；所述参数解包括阶跃项拟合参数、速度项拟合参数和周期项拟合参数；
采用随机模型法对所述参数解进行精度评估，得到中误差；所述中误差包括第二残余时间序列的中误差值和所述参数解的中误差值；所述第二残余时间序列为去除阶跃项、速度项和周期项的GNSS坐标时间序列；
根据所述中误差从所述参数解中选取目标速度项拟合参数；所述目标速度项拟合参数为所述参数解中精度高于设定精度的速度项拟合参数；
根据所述目标速度项拟合参数，采用反距离加权法计算所述目标区域内各格网点的速度值，得到所述目标区域的格网速度场。
2.根据权利要求1所述的一种考虑非线性因素的区域格网速度场构建方法，其特征在于，所述采用稳健最小二乘法对所述GNSS坐标时间序列进行线性拟合，得到第一残余时间序列，具体包括：
采用滑动中位数法对所述GNSS坐标时间序列进行粗差剔除；
采用稳健最小二乘法对粗差剔除后的GNSS坐标时间序列进行线性拟合，得到第一残余时间序列。
3.根据权利要求1所述的一种考虑非线性因素的区域格网速度场构建方法，其特征在于，所述基于所述周期项信息和阶跃项信息，采用稳健最小二乘法对所述GNSS坐标时间序列进行线性拟合，得到所述GNSS坐标时间序列的参数解，具体包括：
采用滑动中位数法对所述GNSS坐标时间序列进行粗差剔除和阶跃探测，得到粗差剔除后的GNSS坐标时间序列和阶跃项信息；
基于所述周期项信息和所述阶跃项信息，采用稳健最小二乘法对所述粗差剔除后的GNSS坐标时间序列进行线性拟合，得到所述GNSS坐标时间序列的参数解。
4.根据权利要求1所述的一种考虑非线性因素的区域格网速度场构建方法，其特征在于，所述采用改进后的周期图谱法对所述第一残余时间序列进行周期项提取，得到周期项信息，具体包括：
在周期图谱法的拟合模型中增加时间平移不变量，得到改变后的拟合模型；
采用所述改变后的拟合模型对所述第一残余时间序列进行分析，构建振幅因子矩阵；
在当前迭代次数下，基于当前迭代次数下的第一残余时间序列的观测误差和所述振幅因子矩阵，构建当前迭代次数下权阵因子；
根据当前迭代次数下权阵因子计算当前迭代次数下的功率谱，并由当前迭代次数下的功率谱提取当前迭代次数下的主周期和对应的振幅；
判断当前迭代次数是否达到设定最大迭代次数；所述设定最大迭代次数等于根据所述GNSS坐标时间序列设定的给定周期个数；
若是，则根据所述迭代次数下的主周期和对应的振幅确定周期项信息；
若否，则进行下一次迭代。
5.一种考虑非线性因素的区域格网速度场构建系统，其特征在于，包括：
时间序列获取模块，用于获取测站测得的目标区域的GNSS坐标时间序列；
第一拟合模块，用于采用稳健最小二乘法对所述GNSS坐标时间序列进行线性拟合，得到第一残余时间序列；所述第一残余时间序列为去除线性速度项的GNSS坐标时间序列；
周期项提取模块，用于采用改进后的周期图谱法对所述第一残余时间序列进行周期项提取，得到周期项信息；所述周期项信息包括周期个数和周期长度；所述改进后的周期图谱算法为引入时间平移不变量，并以加权迭代的方式实现的周期图谱法；
第二拟合模块，用于基于所述周期项信息和阶跃项信息，采用稳健最小二乘法对所述GNSS坐标时间序列进行线性拟合，得到所述GNSS坐标时间序列的参数解；所述阶跃项信息为阶跃发生的历元时刻；所述参数解包括阶跃项拟合参数、速度项拟合参数和周期项拟合参数；
精度评估模块，用于采用随机模型法对所述参数解进行精度评估，得到中误差；所述中误差包括第二残余时间序列的中误差值和所述参数解的中误差值；所述第二残余时间序列为去除阶跃项、速度项和周期项的GNSS坐标时间序列；
速度项拟合参数选取模块，用于根据所述中误差从所述参数解中选取目标速度项拟合参数；所述目标速度项拟合参数为所述参数解中精度高于设定精度的速度项拟合参数；
速度场构建模块，用于根据所述目标速度项拟合参数，采用反距离加权法计算所述目标区域内各格网点的速度值，得到所述目标区域的格网速度场。
6.根据权利要求5所述的一种考虑非线性因素的区域格网速度场构建系统，其特征在于，所述第一拟合模块，具体包括：
粗差剔除单元，用于采用滑动中位数法对所述GNSS坐标时间序列进行粗差剔除；
第一拟合单元，用于采用稳健最小二乘法对粗差剔除后的GNSS坐标时间序列进行线性拟合，得到第一残余时间序列。
7.根据权利要求5所述的一种考虑非线性因素的区域格网速度场构建系统，其特征在于，所述第二拟合模块，具体包括：
预处理单元，用于采用滑动中位数法对所述GNSS坐标时间序列进行粗差剔除和阶跃探测，得到粗差剔除后的GNSS坐标时间序列和阶跃项信息；
第二拟合单元，用于基于所述周期项信息和所述阶跃项信息，采用稳健最小二乘法对所述粗差剔除后的GNSS坐标时间序列进行线性拟合，得到所述GNSS坐标时间序列的参数解。
8.根据权利要求5所述的一种考虑非线性因素的区域格网速度场构建系统，其特征在于，所述周期项提取模块，具体包括：
拟合模型更新单元，用于在周期图谱法的拟合模型中增加时间平移不变量，得到改变后的拟合模型；
振幅因子矩阵构建单元，用于采用所述改变后的拟合模型对所述第一残余时间序列进行分析，构建振幅因子矩阵；
权阵因子构建单元，用于在当前迭代次数下，基于当前迭代次数下的第一残余时间序列的观测误差和所述振幅因子矩阵，构建当前迭代次数下权阵因子；
周期信息提取单元，用于根据当前迭代次数下权阵因子计算当前迭代次数下的功率谱，并由当前迭代次数下的功率谱提取当前迭代次数下的主周期和对应的振幅；
迭代更新单元，用于判断当前迭代次数是否达到设定最大迭代次数；所述设定最大迭代次数等于根据所述GNSS坐标时间序列设定的给定周期个数；若是，则根据所述迭代次数下的主周期和对应的振幅确定周期项信息；若否，则进行下一次迭代。
CN202210194847.2A 2022-03-02 2022-03-02 一种考虑非线性因素的区域格网速度场构建方法及系统 Active CN114253962B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210194847.2A CN114253962B (zh) 2022-03-02 2022-03-02 一种考虑非线性因素的区域格网速度场构建方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210194847.2A CN114253962B (zh) 2022-03-02 2022-03-02 一种考虑非线性因素的区域格网速度场构建方法及系统

Publications (2)

Publication Number Publication Date
CN114253962A true CN114253962A (zh) 2022-03-29
CN114253962B CN114253962B (zh) 2022-05-17

Family

ID=80797233

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210194847.2A Active CN114253962B (zh) 2022-03-02 2022-03-02 一种考虑非线性因素的区域格网速度场构建方法及系统

Country Status (1)

Country Link
CN (1) CN114253962B (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116091832A (zh) * 2023-02-16 2023-05-09 哈尔滨工业大学 基于自编码器网络的肿瘤细胞切片高光谱图像分类方法
CN116165689A (zh) * 2022-06-13 2023-05-26 兰州交通大学 一种基于gnss时间序列自适应提取构造变形变化特征的贝叶斯方法
CN117388872A (zh) * 2023-09-05 2024-01-12 武汉大学 一种北斗地基增强系统参考站坐标框架维持方法和系统

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102571652A (zh) * 2012-01-13 2012-07-11 中国科学院国家授时中心 一种gnss基带信号的评估方法
CN105572703A (zh) * 2015-12-17 2016-05-11 武汉大学 一种gps时间序列广义共模误差提取方法
CN109188466A (zh) * 2018-09-29 2019-01-11 华东交通大学 一种顾及非线性变化的gnss基准站地壳运动速度场估计方法
CN110082787A (zh) * 2019-04-11 2019-08-02 华东师范大学 一种从gnss时间序列中提取周日半周日海潮信号的方法
CN110398753A (zh) * 2019-06-28 2019-11-01 武汉大学 Gnss测站坐标时间序列周期性探测方法及系统
CN110632625A (zh) * 2019-08-19 2019-12-31 中国矿业大学 一种gnss时间序列阶跃探测与修复方法
US10852439B1 (en) * 2020-04-30 2020-12-01 Beihang University Global ionospheric total electron content prediction system
CN112799101A (zh) * 2021-01-29 2021-05-14 华东师范大学 一种构建gnss区域大地参考框架的方法

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102571652A (zh) * 2012-01-13 2012-07-11 中国科学院国家授时中心 一种gnss基带信号的评估方法
CN105572703A (zh) * 2015-12-17 2016-05-11 武汉大学 一种gps时间序列广义共模误差提取方法
CN109188466A (zh) * 2018-09-29 2019-01-11 华东交通大学 一种顾及非线性变化的gnss基准站地壳运动速度场估计方法
CN110082787A (zh) * 2019-04-11 2019-08-02 华东师范大学 一种从gnss时间序列中提取周日半周日海潮信号的方法
CN110398753A (zh) * 2019-06-28 2019-11-01 武汉大学 Gnss测站坐标时间序列周期性探测方法及系统
CN110632625A (zh) * 2019-08-19 2019-12-31 中国矿业大学 一种gnss时间序列阶跃探测与修复方法
US10852439B1 (en) * 2020-04-30 2020-12-01 Beihang University Global ionospheric total electron content prediction system
CN112799101A (zh) * 2021-01-29 2021-05-14 华东师范大学 一种构建gnss区域大地参考框架的方法

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
YINGYING REN 等: "A method based on MTLS and ILSP for GNSS coordinate time series analysis with missing data", 《ADVANCES IN SPACE RESEARCH》 *
YINGYING REN 等: "Preprocessing of GPS Coordinate Sequence Based on Singular Spectrum Analysis", 《IOP CONFERENCE SERIES: EARTH AND ENVIRONMENTAL SCIENCE》 *
任营营 等: "基于K-Means++的省内子块体划分及中国大陆水平相对运动速度场模型的建立与分析", 《地球物理学报》 *
任营营 等: "基于局部无缝Delaunay三角网反距离加权法构建中国大陆速度场", 《武汉大学学报信息科技版》 *
沈连丰等: "面向自动驾驶的车辆精确实时定位算法", 《电子与信息学报》 *
王虎 等: "大规模GNSS网数据处理一体化方案与中国大陆水平格网速度场模型构建研究", 《大地测量与地球动力学》 *

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116165689A (zh) * 2022-06-13 2023-05-26 兰州交通大学 一种基于gnss时间序列自适应提取构造变形变化特征的贝叶斯方法
CN116091832A (zh) * 2023-02-16 2023-05-09 哈尔滨工业大学 基于自编码器网络的肿瘤细胞切片高光谱图像分类方法
CN116091832B (zh) * 2023-02-16 2023-10-20 哈尔滨工业大学 基于自编码器网络的肿瘤细胞切片高光谱图像分类方法
CN117388872A (zh) * 2023-09-05 2024-01-12 武汉大学 一种北斗地基增强系统参考站坐标框架维持方法和系统
CN117388872B (zh) * 2023-09-05 2024-03-19 武汉大学 一种北斗地基增强系统参考站坐标框架维持方法和系统

Also Published As

Publication number Publication date
CN114253962B (zh) 2022-05-17

Similar Documents

Publication Publication Date Title
CN114253962B (zh) 一种考虑非线性因素的区域格网速度场构建方法及系统
Huybers Glacial variability over the last two million years: an extended depth-derived agemodel, continuous obliquity pacing, and the Pleistocene progression
Edwards et al. Determination of site amplification from regional seismicity: application to the Swiss National Seismic Networks
CN110781169B (zh) 自适应多源InSAR监测地面沉降时间序列数据拼接方法及系统
CN106814378B (zh) 一种gnss位置时间序列周期特性挖掘方法
Wang et al. An enhanced singular spectrum analysis method for constructing nonsecular model of GPS site movement
EP2304466A1 (en) Identification and analysis of persistent scatterers in series of sar images
Wang et al. An effective toolkit for the interpolation and gross error detection of GPS time series
Shao et al. What the exercise of the SPICE source inversion validation BlindTest 1 did not tell you
CN110187384B (zh) 贝叶斯时移地震差异反演方法及装置
Xu Reconstruction of gappy GPS coordinate time series using empirical orthogonal functions
CN112781616A (zh) 星敏感器在轨测量低频误差分析方法、装置和存储介质
Bao et al. Filling missing values of multi-station GNSS coordinate time series based on matrix completion
Jiang et al. Effect of removing the common mode errors on linear regression analysis of noise amplitudes in position time series of a regional GPS network & a case study of GPS stations in Southern California
Xiang et al. Characterizing the seasonal hydrological loading over the asian continent using GPS, GRACE, and hydrological model
Malkin Application of the Allan variance to time series analysis in astrometry and geodesy: A review
Aydin et al. Ability of GPS PPP in 2D deformation analysis with respect to GPS network solution
Ran et al. A truncated nuclear norm regularization model for signal extraction from GNSS coordinate time series
Zhang Temporarily coherent point SAR interferometry
Abd El-Gelil et al. Frequency-dependent atmospheric pressure admittance of superconducting gravimeter records using least squares response method
Liu et al. A constrained small baseline subsets (CSBAS) InSAR technique for multiple subsets
Hábel et al. A new tidal analysis of superconducting gravity observations in Western and Central Europe
Bobachev et al. Estimating the Error in Solving the Inverse VES Problem for Precision Investigations of Time Variations in a Geoelectric Section with a Strong Seasonal Effect
CN115598702B (zh) 一种地热资源热储空间构造分布的探测方法及装置
Li et al. Comprehensive analysis of the effects of common mode error on the position time series of a regional GPS network

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