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

1.本发明涉及大地测量技术领域，特别是涉及一种考虑非线性因素的区域格网速度场构建方法及系统。


背景技术：

2.gnss坐标时间序列是一组按照时间顺序排列的基准站坐标组合，其包含了丰富的信息，不仅可以反映测站的线性变化运动，也可以反映出测站存在的非线性变化。线性变化主要表现为速度信号，其反映的是测站受同一区域的构造应力场控制的构造运动，而非线性变化主要表现为周期性信号，其反映的是测站受到的非潮汐海洋负载、大气负载、水文负载、冰期后回弹以及区域性共模误差等地球物理效应的作用。此外，坐标序列中还存在由地壳运动、仪器更换等因素引起的阶跃或震后变形信号。因此，gnss坐标时间序列的分析与建模，特别是研究非线性信号的变化特征，可以更精确地分离测站速度信息，有助于合理解释板块构造运动、建立和维持动态地球参考框架，而且还能构建更高精度的区域格网速度场模型，具有重要的理论研究意义和实际应用价值。正因如此，gnss坐标时间序列分析理论与应用研究成为当前大地测量学和地球物理学等领域的研究热点。
3.对于时间序列中非线性信号的估计，目前常用方法是直接考虑时间序列中的周年项和半年周项周期信号，与速度参数一起利用最小二乘理论进行估计，这一方法的明显不足之处是周期信号考虑不全，且每个测站包含的周期信号特征也有一定的差异。因此，为了建立更加准确的基准站非线性运动模型，有必要对序列中的非线性信号进行分析。而周期图谱方法，正是一种是有效适用于非均匀时间序列的周期分析方法。目前，该算法已广泛应用于天文、经济、地球物理和生物医学等学科领域的非均匀实验观测数据的频谱分析。但是，由于受到序列的非均匀性、有限长度等因素，该算法在傅里叶变换的功率谱中会产生虚假谱峰，此外，由于噪声的影响，周期信号的振幅和相位也可能存在一定的误差。可能正因如此，周期图谱法在gnss坐标时间序列分析中仅仅作为一种辅助手段，并未得到充分和广泛的应用。基于以上考虑，gls算法被提出，该算法在一定程度上弥补了周期图谱法存在的缺点。但是，上述提及的速度场求解方法仍然需要预先对带有缺失值的非均匀gnss原始坐标时间序列进行插补处理才能进行后续分析，因此速度场的精度有待进一步提高。


技术实现要素：

4.基于此，本发明实施例提供一种考虑非线性因素的区域格网速度场构建方法及系统，无需对带有缺失值的非均匀gnss原始坐标时间序列进行插补处理，以提高速度场的精度。
5.为实现上述目的，本发明提供了如下方案：一种考虑非线性因素的区域格网速度场构建方法，包括：获取大地测量过程中测站测得的目标区域的gnss坐标时间序列；采用稳健最小二乘法对所述gnss坐标时间序列进行线性拟合，得到第一残余时间
序列；所述第一残余时间序列为去除线性速度项的gnss坐标时间序列；采用改进后的周期图谱法对所述第一残余时间序列进行周期项提取，得到周期项信息；所述周期项信息包括周期个数和周期长度；所述改进后的周期图谱算法为引入时间平移不变量，并以加权迭代的方式实现的周期图谱法；基于所述周期项信息和阶跃项信息，采用稳健最小二乘法对所述gnss坐标时间序列进行线性拟合，得到所述gnss坐标时间序列的参数解；所述阶跃项信息为阶跃发生的历元时刻；所述参数解包括阶跃项拟合参数、速度项拟合参数和周期项拟合参数；采用随机模型法对所述参数解进行精度评估，得到中误差；所述中误差包括第二残余时间序列的中误差值和所述参数解的中误差值；所述第二残余时间序列为去除阶跃项、速度项和周期项的gnss坐标时间序列；根据所述中误差从所述参数解中选取目标速度项拟合参数；所述目标速度项拟合参数为所述参数解中精度高于设定精度的速度项拟合参数；根据所述目标速度项拟合参数，采用反距离加权法计算所述目标区域内各格网点的速度值，得到所述目标区域的格网速度场。
6.可选的，所述采用稳健最小二乘法对所述gnss坐标时间序列进行线性拟合，得到第一残余时间序列，具体包括：采用滑动中位数法对所述gnss坐标时间序列进行粗差剔除；采用稳健最小二乘法对粗差剔除后的gnss坐标时间序列进行线性拟合，得到第一残余时间序列。
7.可选的，所述基于所述周期项信息和阶跃项信息，采用稳健最小二乘法对所述gnss坐标时间序列进行线性拟合，得到所述gnss坐标时间序列的参数解，具体包括：采用滑动中位数法对所述gnss坐标时间序列进行粗差剔除和阶跃探测，得到粗差剔除后的gnss坐标时间序列和阶跃项信息；基于所述周期项信息和所述阶跃项信息，采用稳健最小二乘法对所述粗差剔除后的gnss坐标时间序列进行线性拟合，得到所述gnss坐标时间序列的参数解。
8.可选的，所述采用改进后的周期图谱法对所述第一残余时间序列进行周期项提取，得到周期项信息，具体包括：在周期图谱法的拟合模型中增加时间平移不变量，得到改变后的拟合模型；采用所述改变后的拟合模型对所述第一残余时间序列进行分析，构建振幅因子矩阵；在当前迭代次数下，基于当前迭代次数下的第一残余时间序列的观测误差和所述振幅因子矩阵，构建当前迭代次数下权阵因子；根据当前迭代次数下权阵因子计算当前迭代次数下的功率谱，并由当前迭代次数下的功率谱提取当前迭代次数下的主周期和对应的振幅；判断当前迭代次数是否达到设定最大迭代次数；所述设定最大迭代次数等于根据所述gnss坐标时间序列设定的给定周期个数；若是，则根据所述迭代次数下的主周期和对应的振幅确定周期项信息；若否，则进行下一次迭代。
9.本发明还提供了一种考虑非线性因素的区域格网速度场构建系统，包括：
时间序列获取模块，用于获取测站测得的目标区域的gnss坐标时间序列；第一拟合模块，用于采用稳健最小二乘法对所述gnss坐标时间序列进行线性拟合，得到第一残余时间序列；所述第一残余时间序列为去除线性速度项的gnss坐标时间序列；周期项提取模块，用于采用改进后的周期图谱法对所述第一残余时间序列进行周期项提取，得到周期项信息；所述周期项信息包括周期个数和周期长度；所述改进后的周期图谱算法为引入时间平移不变量，并以加权迭代的方式实现的周期图谱法；第二拟合模块，用于基于所述周期项信息和阶跃项信息，采用稳健最小二乘法对所述gnss坐标时间序列进行线性拟合，得到所述gnss坐标时间序列的参数解；所述阶跃项信息为阶跃发生的历元时刻；所述参数解包括阶跃项拟合参数、速度项拟合参数和周期项拟合参数；精度评估模块，用于采用随机模型法对所述参数解进行精度评估，得到中误差；所述中误差包括第二残余时间序列的中误差值和所述参数解的中误差值；所述第二残余时间序列为去除阶跃项、速度项和周期项的gnss坐标时间序列；速度项拟合参数选取模块，用于根据所述中误差从所述参数解中选取目标速度项拟合参数；所述目标速度项拟合参数为所述参数解中精度高于设定精度的速度项拟合参数；速度场构建模块，用于根据所述目标速度项拟合参数，采用反距离加权法计算所述目标区域内各格网点的速度值，得到所述目标区域的格网速度场。
10.可选的，所述第一拟合模块，具体包括：粗差剔除单元，用于采用滑动中位数法对所述gnss坐标时间序列进行粗差剔除；第一拟合单元，用于采用稳健最小二乘法对粗差剔除后的gnss坐标时间序列进行线性拟合，得到第一残余时间序列。
11.可选的，所述第二拟合模块，具体包括：预处理单元，用于采用滑动中位数法对所述gnss坐标时间序列进行粗差剔除和阶跃探测，得到粗差剔除后的gnss坐标时间序列和阶跃项信息；第二拟合单元，用于基于所述周期项信息和所述阶跃项信息，采用稳健最小二乘法对所述粗差剔除后的gnss坐标时间序列进行线性拟合，得到所述gnss坐标时间序列的参数解。
12.可选的，所述周期项提取模块，具体包括：拟合模型更新单元，用于在周期图谱法的拟合模型中增加时间平移不变量，得到改变后的拟合模型；振幅因子矩阵构建单元，用于采用所述改变后的拟合模型对所述第一残余时间序列进行分析，构建振幅因子矩阵；权阵因子构建单元，用于在当前迭代次数下，基于当前迭代次数下的第一残余时间序列的观测误差和所述振幅因子矩阵，构建当前迭代次数下权阵因子；周期信息提取单元，用于根据当前迭代次数下权阵因子计算当前迭代次数下的功率谱，并由当前迭代次数下的功率谱提取当前迭代次数下的主周期和对应的振幅；迭代更新单元，用于判断当前迭代次数是否达到设定最大迭代次数；所述设定最
大迭代次数等于根据所述gnss坐标时间序列设定的给定周期个数；若是，则根据所述迭代次数下的主周期和对应的振幅确定周期项信息；若否，则进行下一次迭代。
13.与现有技术相比，本发明的有益效果是：本发明实施例提出了一种考虑非线性因素的区域格网速度场构建方法及系统，对目标区域的gnss坐标时间序列，采用稳健最小二乘法进行线性拟合；采用改进后的周期图谱法对拟合后的第一残余时间序列进行周期项提取，得到周期项信息；改进后的周期图谱算法为引入时间平移不变量，并以加权迭代的方式实现的周期图谱法；基于周期项信息和阶跃项信息，采用稳健最小二乘法对gnss坐标时间序列进行线性拟合，得到参数解；采用随机模型法对参数解进行精度评估，得到中误差；根据中误差从参数解中选取目标速度项拟合参数；根据目标速度项拟合参数，采用反距离加权法计算目标区域内各格网点的速度值，得到格网速度场。本发明无需预先对带有缺失值的非均匀gnss原始坐标时间序列进行插补处理即可进行后续分析，同时考虑了序列中的非线性变化因素（同时考虑了周期项、阶跃项和速度项），进一步提高了速度解的精度，从而为区域格网速度场自动一体化构建方法提供一套完善可行的依据。
附图说明
14.为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
15.图1为本发明实施例提供的考虑非线性因素的区域格网速度场构建方法的流程图；图2为本发明实施例提供的考虑非线性因素的区域格网速度场构建方法的实现框图；图3为本发明实施例提供的考虑非线性因素的区域格网速度场构建系统的结构图。
具体实施方式
16.下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
17.为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。
18.图1为本发明实施例提供的考虑非线性因素的区域格网速度场构建方法的流程图。
19.参见图1，本实施例的考虑非线性因素的区域格网速度场构建方法，包括：步骤101：获取测站测得的目标区域的gnss坐标时间序列。
20.步骤102：采用稳健最小二乘法对所述gnss坐标时间序列进行线性拟合，得到第一
残余时间序列；所述第一残余时间序列为去除线性速度项的gnss坐标时间序列。
21.步骤103：采用改进后的周期图谱法对所述第一残余时间序列进行周期项提取，得到周期项信息；所述周期项信息包括周期个数和周期长度；所述改进后的周期图谱算法为引入时间平移不变量，并以加权迭代的方式实现的周期图谱法。
22.步骤104：基于所述周期项信息和阶跃项信息，采用稳健最小二乘法对所述gnss坐标时间序列进行线性拟合，得到所述gnss坐标时间序列的参数解。所述阶跃项信息为阶跃发生的历元时刻；所述参数解包括阶跃项拟合参数、速度项拟合参数和周期项拟合参数。
23.步骤105：采用随机模型法对所述参数解进行精度评估，得到中误差；所述中误差包括第二残余时间序列的中误差值和所述参数解的中误差值；所述第二残余时间序列为去除阶跃项、速度项和周期项的gnss坐标时间序列。
24.步骤106：根据所述中误差从所述参数解中选取目标速度项拟合参数；所述目标速度项拟合参数为所述参数解中精度高于设定精度的速度项拟合参数。
25.步骤107：根据所述目标速度项拟合参数，采用反距离加权法计算所述目标区域内各格网点的速度值，得到所述目标区域的格网速度场。
26.作为一种可选的实施方式，所述步骤102，具体包括：1）采用滑动中位数法对所述gnss坐标时间序列进行粗差剔除。
27.具体的，采用基于iqr准则的滑动中位数粗差探测法进行粗差剔除。该方法依据的合理假设条件为：短时间（数天或数周）内gnss测站的坐标处于稳定状态，通常情况下gnss静态观测站的位置是稳定不变的。正是由于中位数不受序列中离群值的影响，因此可以在设定合理滑动窗口的条件下，可以有效探测序列中的离群值。其主要实现步骤为：s21、设定合适的滑动窗口长度w，以序列位置i为中心前后分别提取依据w/2长度的子序列，称为sub1和sub2；s22、在子序列内依据iqr准则判断i点是否为离群值，若为离群值则将其视为粗差进行剔除；s23、重复步骤s102，逐步滑动窗口至整个序列，判断并剔除所有的离群值；s24、第一次粗差探测完毕并剔除粗差后，继续进行第二次粗差探测，迭代至序列不在含有粗差为止。
28.2）采用稳健最小二乘法对粗差剔除后的gnss坐标时间序列进行线性拟合，得到第一残余时间序列。具体的：采用稳健最小二乘模型进行线性拟合并提取残余序列的。利用稳健最小二乘进行线性拟合的原理如下：仅考虑粗差剔除后的gnss坐标时间序列中包含的线性信号，观测方程为：(1)式中，为时刻的测站分量坐标值，即粗差剔除后的gnss坐标时间序列，为常量项，线性速率系数，为第个时刻，为测量误差。
29.改写为误差方程：(2)式中，为残差值，为拟合估值，为观测值，为观测历元个数。
30.误差方程的矩阵形式为：
(3)式中，,,,。
31.根据公式（4）准则，采用最小二乘拟合估计，以得到第一次未知数的解。
32.(4)(5)式中，，表示第i个观测值的观测误差。
33.根据残差和定权公式（6）重新确定各观测值新的权因子，进行下一轮计算，直至前后两次估计值的变化小于限差时结束迭代。迭代完成后，将最后一次的参数结果代回公式（3），得到该步骤中第一残余时间序列。
34.(6)式中，为权值，为阈值常量，一般取二倍或三倍的单位权中误差。
35.作为一种可选的实施方式，步骤103是采用改进的周期图谱法提取周期项先验信息的，其主要原理如下：周期图谱法（lomb-scargle periodogram，lsp）是一种基于离散傅里叶变换的周期提取方法。该算法在一定程度可以解决非均匀采样间隔对周期信号的影响，避免对非均匀采样时间序列进行插补处理，且考虑了非均匀采样对幅度和相位对信号产生的影响。该算法的基本原理是通过最小二乘法将一系列三角函数的进行线性组合来拟合时间序列，在此基础上，将信号特征从时域转换到频域上。
36.对于非均匀时间序列，即步骤102得到的第一残余时间序列结果，该第一残余时间序列的平均值为0，定义拟合方程为：(7)式中，为离散序列的采样时间，为时间序列的数据统计量。离散测试频率，，定义为序列的极限频率，一般不大于尼奎斯特频率，测试频率的取样步长，测试频率数量。表示频率分量f k
的正弦变化的幅度，表示频率分量f k
的余弦变化的幅度。为第i个观测值的观测噪声。
37.为简便后续公式的推导，定义以下变量（仅作为临时变量辅助公式推导）：
(8)根据间接平差原理，构建误差方程如下：(9)式中，为残差向量，为系数矩阵，为待估参数，为观测向量。
38.基于最小二乘原理，得到最佳估值：(10)其中，。
39.则lsp的功率谱定义为：(11)式中，为频率的功率谱值。
40.为了避免序列取样时间的平移对谱结构产生影响，引入时间平移不变量，其他参数变量含义同公式（7），建立新的拟合模型：(12)加入重新定义以下变量：
(13)选择值使得成立，推导可得：(14)基于最小二乘原理，得到最佳估值：(15)则lsp的功率谱定义为：(16)使用显著性等级用以评价功率谱的质量，代表频谱的虚警概率。当=0.1时，频谱的置信度则为90%。
41.(17)
其中，为估计的频谱值，为频谱中包含的独立周期的个数。
42.传统的lsp算法虽然有效，但仍然存在以下缺陷：
①
没有考虑观测误差对结果造成的影响；
②
预先假设原始序列与拟合所采用的正弦函数的均值相同；
③
非均匀序列往往会在真实信号两侧出现虚假谱峰；
④
没有考虑多重信号之间相互调制对结果造成的影响。
43.第一个缺陷易导致周期成分的振幅和相位偏离真实值；第二个缺陷易导致周期成分的振幅产生系统误差。针对这两个问题，zechmeister提出了glsp算法，采用了修正的正弦函数来对时间序列进行拟合。glsp算法与lsp算法相似，主要区别在于glsp算法通过加入常数的正弦三角函数对时间序列进行拟合分析，且计算时进一步考虑了观测误差的时间序列的影响。因此，从理论上说，glsp算法的周期信号估计准确度应优于lsp算法，但glsp算法依旧没有考虑后两个缺陷所产生的影响，为此，本发明将提出一种新的ilsp算法，主要采用加权和迭代两种方法以有效减小各频率信号之间的相互调制和观测噪声对功率谱结果产生的影响。考虑到其他频率分量和观测误差或噪声对该频率信号的影响，定义新的加权因子为： (18)其中，,，表示频率分量f m
的正弦变化的幅度，表示频率分量f m
的余弦变化的幅度。
44.那么得到加权后的参数估值和功率谱： (19) (20)式中，为参数估值，为系数矩阵，为加权矩阵，为观测向量，为频率的功率谱值。
45.根据频率与功率谱的一一对应关系，寻找出功率谱最大值对应的频率，进而根据公式(22)计算出对应的周期。
46.=max[ (21)式中，n为周期个数。
[0047] (22)基于虚假谱峰值低于主峰值这一特性，采取迭代法消除虚假谱峰的影响，即第次迭代仅提取序列中主峰值对应的周期和振幅，然后从序列中减去对应的周期信号以得到残差序列重复计算，直至达到设定的迭代次数（即目标周期个数）结束迭代。每一次迭代都得
到一个周期值（周期长度），n次迭代完成后，得到n个周期（周期个数），即得到最终的周期项信息。
[0048]
基于上述原理，所述步骤103，具体包括：s31：在周期图谱法的拟合模型中增加时间平移不变量，得到改变后的拟合模型。
[0049]
s32：采用所述改变后的拟合模型对所述第一残余时间序列进行分析，构建振幅因子矩阵。
[0050]
s33：在当前迭代次数下，基于当前迭代次数下的第一残余时间序列的观测误差和所述振幅因子矩阵，构建当前迭代次数下权阵因子。
[0051]
s34：根据当前迭代次数下权阵因子计算当前迭代次数下的功率谱，并由当前迭代次数下的功率谱提取当前迭代次数下的主周期和对应的振幅。
[0052]
s35：判断当前迭代次数是否达到设定最大迭代次数；所述设定最大迭代次数等于根据所述gnss坐标时间序列设定的给定周期个数。若是，则根据所述迭代次数下的主周期和对应的振幅确定周期项信息；若否，则进行下一次迭代。每一次迭代都得到一个周期值（周期长度），多次迭代完成后，得到相应个数的周期（周期个数），这样得到最终的周期项信息。
[0053]
步骤103，一个更为具体的实现过程如下：改进周期图谱法需要经过迭代提取周期信息。结合加权与迭代法，改进周期图谱法的具体实现步骤如下：
①
给定原始时间序列和待探测的周期个数n；
②
进行时间序列自动化预处理（主要包括粗差剔除与阶跃修复），基于稳健最小二乘模型估计并消除序列中的线性趋势项；
③
进行lsp分析，以得到频率对应的振幅，构造振幅因子矩阵。若n为0，进行功率谱的显著性评估，取=0.1，获取置信度大于99%的谱峰个数，并设为n。
④
进行ilsp分析，基于序列观测误差和振幅因子，构造权阵因子，计算功率谱，提取主周期和对应的振幅。
⑤
去除序列中的主周期项，以得到残差序列。
⑥
重复
④
和
⑤
两个步骤n次结束迭代。
[0054]
作为一种可选的实施方式，所述步骤104，具体包括：1）采用滑动中位数法对所述gnss坐标时间序列进行粗差剔除和阶跃探测，得到粗差剔除后的gnss坐标时间序列和阶跃项信息。
[0055]
粗差剔除过程是采用基于iqr准则的滑动中位数粗差探测法实现的，其主要实现步骤与步骤102中粗差剔除的过程相同，在此不再赘述。
[0056]
阶跃探测（阶跃项的定位）过程是基于k-medoids的滑动中位数阶跃定位法进行的，该方法依据的合理假设条件为：短时间（数天或数周）内gnss测站的坐标处于稳定状态，通常情况下gnss静态观测站的位置是不变的。该方法同时对gnss测站三个方向的坐标时间序列进行阶跃定位，考虑了序列间的相关性，使用滑动中位数法可以避免离群值的影响，而且使用k-medoids对阶跃发生历元进行聚类进一步保证了结果的可靠性，因此是一种较为简便有效的阶跃定位策略。阶跃探测的主要实现步骤为：步骤s41、设定合适的滑动窗口长度w、阶跃探测阈值e，以序列位置i为中心前后分别提取w/2长度的子序列sub1和sub2，计算sub1和sub2的中位数med1和med2的绝对值之差，若其绝对值大于探测阈值e，则判断位置i处为阶跃发生历元。
[0057]
步骤s42、逐步滑动窗口至整个序列，完成gnss坐标三个方向（neu）的序列的阶跃
探测，然后合并三个方向的阶跃发生历元并去除重复值。
[0058]
步骤s43、确定阶跃发生历元的种类个数，对历元结果进行k-medoids聚类分析，求取每一类阶跃发生历元的中位数，将其作为准确的阶跃发生历元。
[0059]
步骤s44、确定阶跃发生历元后，以neu序列阶跃发生位置为中心前后分别提取w/2长度的子序列sub1和sub2，计算sub1和sub2的中位数med1和med2的差值作为阶跃值。
[0060]
在步骤s43中，是采用k-medoids聚类算法对阶跃发生历元进行聚类分析的。k-medoids聚类算法是一种分区方法，其类似于k-means，两种方法的目标都是将一组测量值或观察值划分为k个子集或集群，以便子集最小化测量值与测量值集群中心之间的距离总和。在k-means聚类算法中，子集的中心是子集中测量的平均值，而k-medoids聚类算法中，子集的中心是子集中测量的中位数。相较于平均值，中位值不易受到集合中离群值的影响，因此k-medoids聚类算法通常用于对异常数据、任意距离度量、平均值或中位数等没有明确定义的数据具有鲁棒性要求的领域。
[0061]
2）基于所述周期项信息和所述阶跃项信息，采用稳健最小二乘法对所述粗差剔除后的gnss坐标时间序列进行线性拟合，得到所述gnss坐标时间序列的参数解。
[0062]
采用稳健最小二乘法进行参数解拟合的，其主要原理如下：在原始gnss坐标时间序列中，通常需要考虑测站趋势项（线性速度）、周期变化项（主要考虑年周期和半年周期）、非地震因素（设备更换、天线高测量误差、相位中心建模误差或其他人为以及软件误差）或地震因素（地震同震断裂）引起的阶跃跳变项、震后形变项（通常呈指数或对数变化形式），以及一些未模型化的误差项。暂不考虑未模型化的误差项，且由于历元之间坐标分量的相关性很小，历元时刻的坐标分量可以详细模型化表示为： (23)式中，为初始时刻对应的截距。为从为参考的历元时刻，单位是年。线性速率为震间测站的长期的构造运动。c、d、e、f为年周期和半年周期项系数，估计周期项参数则需要足够多的数据（至少12个月），年周期和半年周期的信号用表示，为幅度，为角速率，为一月一日，为相位。为发生阶跃的总次数，为历元时刻由于同震或非同震变化引起的阶跃量。为阶梯函数。为震后形变函数。为观测噪声。
[0063]
函数表示为： (24)式中，为发生阶跃的时刻。
[0064]
函数通常可由指数或对数函数表示：
ꢀ
(25)式中，为地震发生时刻的阶跃系数，为地震衰减因子。
[0065]
gnss坐标时间序列模型中各事件的发生时刻可从地震目录、观测日志、自动检测算法或目视检查确定。由于自动检测算法可靠性差，目视检查法实施起来难度较大，本发明中主要根据地震目录和测站日志来确定测站坐标的阶跃发生时刻和震后形变发生时刻。此外，地震衰减因子通常需要通过其他方法进行单独估计。因此，剩余时间序列系数可以表示为线性模型，然后再次基于稳健最小二乘模型(公式（1）-公式（6）)可得到第二残余时间序列和参数的最佳拟合估值该参数估值中，为速度解（速度项拟合参数），[为周期振幅（周期项拟合参数），为阶跃值（周期项拟合参数），为震后形变系数。
[0066]
其中，步骤105，其主要原理如下：时间序列参数的中误差估值为： (26)式中，为的主对角元素，，为单位权中误差估值，即第二残余时间序列的中误差，有： (27)式中，为多余观测量，为观测值改正数，即步骤104得到的第二残余序列（已去除阶跃项、速度项和周期项）。参数解的中误差值的计算方式与第二残余时间序列的中误差的计算方式类似，在此不再赘述。之后，根据计算得到的中误差从参数解中选取精度高于设定精度的速度项拟合参数，即得到目标速度项拟合参数。
[0067]
其中，步骤107，根据上述步骤得到符合精度要求的速度值（目标速度项拟合参数），基于反距离加权法构建区域格网速度场，其主要原理如下：设空间待插值格网点为，p点邻域内有n个已知离散点。根据离散点的速度值()，通过反距离加权法对p点的速度值进行插值，即：(28)
式中，，表示离散点到待定点的距离；一般取1~2，本文取=1。最终步骤得到格网速度场，即每个格网点的速度值。
[0068]
上述实施例中的区域格网速度场构建方法实现框图如图2所示，由于方法中采用的周期图谱法在一定程度上解决了非均匀采样间隔产生的非周期信号问题，从而避免了对非均匀采样时间序列进行插值处理，同时考虑了非均匀采样对幅度和相位带来的影响，而且改进后的周期图谱法主要采用加权和迭代两种方法进一步考虑了序列观测噪声误差与虚假谱峰的影响，以有效减小各频率信号之间的相互调制和观测噪声对功率谱结果产生的影响。因此，此方法无需预先对带有缺失值的非均匀gnss原始坐标时间序列进行插补处理即可进行后续分析，同时考虑了序列中的非线性变化因素，进一步提高了速度解的精度，从而为区域格网速度场自动一体化构建方法提供一套完善可行的算法依据。
[0069]
本发明还提供了一种考虑非线性因素的区域格网速度场构建系统，参见图3，所述系统，包括：时间序列获取模块301，用于获取测站测得的目标区域的gnss坐标时间序列。
[0070]
第一拟合模块302，用于采用稳健最小二乘法对所述gnss坐标时间序列进行线性拟合，得到第一残余时间序列；所述第一残余时间序列为去除线性速度项的gnss坐标时间序列。
[0071]
周期项提取模块303，用于采用改进后的周期图谱法对所述第一残余时间序列进行周期项提取，得到周期项信息；所述周期项信息包括周期个数和周期长度；所述改进后的周期图谱算法为引入时间平移不变量，并以加权迭代的方式实现的周期图谱法。
[0072]
第二拟合模块304，用于基于所述周期项信息和阶跃项信息，采用稳健最小二乘法对所述gnss坐标时间序列进行线性拟合，得到所述gnss坐标时间序列的参数解；所述阶跃项信息为阶跃发生的历元时刻；所述参数解包括阶跃项拟合参数、速度项拟合参数和周期项拟合参数。
[0073]
精度评估模块305，用于采用随机模型法对所述参数解进行精度评估，得到中误差；所述中误差包括第二残余时间序列的中误差值和所述参数解的中误差值；所述第二残余时间序列为去除阶跃项、速度项和周期项的gnss坐标时间序列。
[0074]
速度项拟合参数选取模块306，用于根据所述中误差从所述参数解中选取目标速度项拟合参数；所述目标速度项拟合参数为所述参数解中精度高于设定精度的速度项拟合参数。
[0075]
速度场构建模块307，用于根据所述目标速度项拟合参数，采用反距离加权法计算所述目标区域内各格网点的速度值，得到所述目标区域的格网速度场。
[0076]
作为一种可选的实施方式，所述第一拟合模块302，具体包括：粗差剔除单元，用于采用滑动中位数法对所述gnss坐标时间序列进行粗差剔除。
[0077]
第一拟合单元，用于采用稳健最小二乘法对粗差剔除后的gnss坐标时间序列进行线性拟合，得到第一残余时间序列。
[0078]
作为一种可选的实施方式，所述第二拟合模块304，具体包括：预处理单元，用于采用滑动中位数法对所述gnss坐标时间序列进行粗差剔除和阶
跃探测，得到粗差剔除后的gnss坐标时间序列和阶跃项信息。
[0079]
第二拟合单元，用于基于所述周期项信息和所述阶跃项信息，采用稳健最小二乘法对所述粗差剔除后的gnss坐标时间序列进行线性拟合，得到所述gnss坐标时间序列的参数解。
[0080]
作为一种可选的实施方式，所述周期项提取模块303，具体包括：拟合模型更新单元，用于在周期图谱法的拟合模型中增加时间平移不变量，得到改变后的拟合模型。
[0081]
振幅因子矩阵构建单元，用于采用所述改变后的拟合模型对所述第一残余时间序列进行分析，构建振幅因子矩阵。
[0082]
权阵因子构建单元，用于在当前迭代次数下，基于当前迭代次数下的第一残余时间序列的观测误差和所述振幅因子矩阵，构建当前迭代次数下权阵因子。
[0083]
周期信息提取单元，用于根据当前迭代次数下权阵因子计算当前迭代次数下的功率谱，并由当前迭代次数下的功率谱提取当前迭代次数下的主周期和对应的振幅。
[0084]
迭代更新单元，用于判断当前迭代次数是否达到设定最大迭代次数；所述设定最大迭代次数等于根据所述gnss坐标时间序列设定的给定周期个数；若是，则根据所述迭代次数下的主周期和对应的振幅确定周期项信息；若否，则进行下一次迭代。
[0085]
本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的系统而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。
[0086]
本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。