CN112184778B - 基于gnss动态定位时间序列分割的短期位移探测方法 - Google Patents

Abstract

本发明属于位移探测技术领域，公开了基于GNSS动态定位时间序列分割的短期位移探测方法，首先通过GNSS动态定位获得原始坐标时间序列；然后利用检测窗口进行时间序列分割，获得可能变化点的集合；针对每个可能变化点，获得可能变化点的统计值，根据统计值及预设阈值判断可能变化点是否为变化点；最后针对每个变化点，获得变化点对应的历次位移估计值，基于历次位移估计值获得最终位移估计值。本发明解决了现有技术无法在短时间内有效地检测和提取位移的问题，能够在短时间内有效地检测和提取位移。

Description

基于GNSS动态定位时间序列分割的短期位移探测方法

技术领域

本发明涉及位移探测技术领域，尤其涉及基于GNSS动态定位时间序列分割的短期位移探测方法。

背景技术

位移测量对于地质灾害、矿山安全、结构健康安全监测以及工业测量领域具有积极的作用。在过去的二十年里，人们对利用全球导航卫星系统(GNSS)进行位移监测的研究激增。与传统位移监测技术相比，GNSS具有不受天气影响，监测点之间不需要通视，高精度，多时间尺度等优势。此外，随着GNSS技术的不断发展，接收机成本不断降低，而且定位精度不断提高。GNSS成为位移监测最有效的手段之一。

目前，在这些应用领域，基于GNSS的位移监测研究大多集中于长期位移分析，而对短期监测的研究较少。

近二十年来的研究为基于GNSS的位移监测提供了重要的信息，但大多只适用于离线应用。到目前为止，对短期位移检测的研究还很少。虽然利用GNSS坐标差可以直接得到位移，但由于测量误差的影响，如何在短时间内有效地检测和提取位移仍然存在挑战。目前，大多数已有的位移检测方法都是逐历元进行的，容易受到粗差和观测噪声的影响。

发明内容

本发明通过提供基于GNSS动态定位时间序列分割的短期位移探测方法，解决了现有技术无法在短时间内有效地检测和提取位移的问题。

本发明提供基于GNSS动态定位时间序列分割的短期位移探测方法，包括以下步骤：

步骤1、通过GNSS动态定位获得原始坐标时间序列；

步骤2、利用检测窗口进行时间序列分割，获得可能变化点的集合；

步骤3、针对每个可能变化点，获得可能变化点的统计值，根据统计值及预设阈值判断可能变化点是否为变化点；

步骤4、针对每个变化点，获得变化点对应的历次位移估计值，基于所述历次位移估计值获得最终位移估计值。

优选的，所述步骤1中还包括：对所述原始坐标时间序列进行去噪处理，得到去噪后的坐标时间序列；所述步骤2中，利用检测窗口对所述去噪后的坐标时间序列进行分割，获得可能变化点的集合。

优选的，所述步骤1中，所述GNSS动态定位为相对实时定位或精密单点实时定位。

优选的，采用小波分析方法对所述原始坐标时间序列进行分解和重构，得到重构的坐标时间序列，所述重构的坐标时间序列作为所述去噪后的坐标时间序列。

优选的，所述小波分析方法包括以下子步骤：

确定小波母函数；

采用离散小波变换对所述原始坐标时间序列进行分解，得到分解的近似分量参数；

基于所述近似分量参数进行小波重构，得到所述重构的坐标时间序列；

其中，对所述原始坐标时间序列进行分解表示为：

其中，x(t)表示原始坐标时间序列，m表示缩放因子，n表示平移因子，ψ(·)表示小波母函数，ψ^*(·)表示小波母函数的共轭，a₀表示指定的固定步长参数，b₀表示大于0的位置参数。

优选的，所述小波母函数采用Daubechies db1小波函数。

优选的，所述步骤2中，在时间序列分割中，根据被τ分割的两个部分的统计信息构造目标函数，所述目标函数表示为：

其中，J(τ)表示目标函数，τ表示检测窗口的分割点，x_i、x_j表示检测窗口中点的坐标值，n表示检测窗口的长度；

通过使所述目标函数最小化获得可能变化点，表示为：

其中，

表示满足目标函数最小化的最佳估值。

优选的，所述步骤3中，对于一个历元，在一个大小为n的检测窗口中，该历元参与窗口分割的次数为n次；在窗口分割过程中，将第i个历元被划分为可能变化点的次数记为count_i，count_i作为该可能变化点的统计值；若该可能变化点的统计值大于预设阈值，则将该历元视为变化点。

优选的，所述步骤4中，对于count_i次被划分为可能变化点的第i个历元，存在count_i次的位移估值，表示为：

Δh_ij＝h_ij0-h_ij1,1≤j≤count_i

其中，Δh_ij表示第i个历元对应的第j次分割的位移估值，count_i表示第i个历元被划分为可能变化点的次数，h_ij0表示第i个历元对应的第j次分割的前半窗口中的坐标均值，h_ij1表示第i个历元对应的第j次分割的后半窗口中的坐标均值；

基于count_i次的位移估值划分位移区域，采用如下公式：

bin_ik＝div(min(Δh_ij),max(Δh_ij),s),1≤k≤s

其中，bin_ik表示第i个历元对应的第k个位移区间，div(.)操作符表示根据count_i次的位移估值中的最大值和最小值均匀划分为指定的份数，s表示份数；

基于落入位移估值最多的位移区间，获得最终位移估计值，表示为：

其中，Δh_i表示最终位移估计值，count(bin_ik)表示Δh_ij落入bin_ik的数量,f为从位移区间到位移的映射函数。

本发明中提供的一个或多个技术方案，至少具有如下技术效果或优点：

在发明中，首先通过GNSS动态定位获得原始坐标时间序列；然后利用检测窗口进行时间序列分割，获得可能变化点的集合；针对每个可能变化点，获得可能变化点的统计值，根据统计值及预设阈值判断可能变化点是否为变化点；最后针对每个变化点，获得变化点对应的历次位移估计值，基于历次位移估计值获得最终位移估计值。本发明能够在短时间内有效地检测和提取位移，且位提取的位移精度较高。

附图说明

图1为本发明实施例提供的基于GNSS动态定位时间序列分割的短期位移探测方法的流程示意图；

图2为本发明实施例提供的基于GNSS动态定位时间序列分割的短期位移探测方法中基于时间序列分割的位移检测原理图；

图3为本发明实施例提供的基于GNSS动态定位时间序列分割的短期位移探测方法中GNSS动态定位以及小波重构的结果示意图；

图4为本发明实施例提供的基于GNSS动态定位时间序列分割的短期位移探测方法中窗口分割可能变化点的示意图；

图5为本发明实施例提供的基于GNSS动态定位时间序列分割的短期位移探测方法中可能变化点统计计数图；

图6为本发明实施例提供的基于GNSS动态定位时间序列分割的短期位移探测方法中可能变化点位移统计结果图。

具体实施方式

本发明旨在探索一种基于GNSS动态定位时间序列分割的短期位移探测方法，以GNSS动态定位获得的当前历元坐标和前一个历元坐标形成的坐标时间序列作为分割的基本处理单元。将从检测窗口中得到的分割点视为潜在的变化点，并对潜在变化点进行统计检验以确定分割点是否为变化点。

为了更好的理解上述技术方案，下面将结合说明书附图以及具体的实施方式对上述技术方案进行详细的说明。

本实施例提供了基于GNSS动态定位时间序列分割的短期位移探测方法，主要包括以下步骤：

步骤1、通过GNSS动态定位获得原始坐标时间序列。

优选的方案中还包括对所述原始坐标时间序列进行去噪处理，得到去噪后的坐标时间序列。例如，采用小波分析方法对所述原始坐标时间序列进行分解和重构，得到重构的坐标时间序列，所述重构的坐标时间序列作为所述去噪后的坐标时间序列。

步骤2、利用检测窗口进行时间序列分割，获得可能变化点的集合。

具体的，在时间序列分割中，根据被τ分割的两个部分的统计信息构造目标函数，所述目标函数表示为：

其中，J(τ)表示目标函数，τ表示检测窗口的分割点，x_i、x_j表示检测窗口中点的坐标值，n表示检测窗口的长度。

通过使所述目标函数最小化获得可能变化点，表示为：

其中，

表示满足目标函数最小化的最佳估值。

若步骤1中进行了去噪处理，则利用检测窗口对所述去噪后的坐标时间序列进行分割，获得可能变化点的集合。

步骤3、针对每个可能变化点，获得可能变化点的统计值，根据统计值及预设阈值判断可能变化点是否为变化点。

具体的，对于一个历元，在一个大小为n的检测窗口中，该历元参与窗口分割的次数为n次；在窗口分割过程中，将第i个历元被划分为可能变化点的次数记为 coun_it，count_i作为该可能变化点的统计值；若该可能变化点的统计值大于预设阈值，则将该历元视为变化点。

步骤4、针对每个变化点，获得变化点对应的历次位移估计值，基于所述历次位移估计值获得最终位移估计值。

具体的，对于count_i次被划分为可能变化点的第i个历元，存在count_i次的位移估值，表示为：

Δh_ij＝h_ij0-h_ij1,1≤j≤count_i

其中，Δh_ij表示第i个历元对应的第j次分割的位移估值，count_i表示第i个历元被划分为可能变化点的次数，h_ij0表示第i个历元对应的第j次分割的前半窗口中的坐标均值，h_ij1表示第i个历元对应的第j次分割的后半窗口中的坐标均值。

基于count_i次的位移估值划分位移区域，采用如下公式：

bin_ik＝div(min(Δh_ij),max(Δh_ij),s),1≤k≤s

其中，bin_ik表示第i个历元对应的第k个位移区间，div(.)操作符表示根据count_i次的位移估值中的最大值和最小值均匀划分为指定的份数，s表示份数。

基于落入位移估值最多的位移区间，获得最终位移估计值，表示为：

其中，Δh_i表示最终位移估计值，count(bin_ik)表示Δh_ij落入bin_ik的数量,f为从位移区间到位移的映射函数。

下面对本发明做进一步的说明。

参见图1，本实施例提供的基于GNSS动态定位时间序列分割的短期位移探测方法主要包括以下几个步骤：动态定位获得坐标时间序列、去噪处理、窗口分割探测可能变化点、判定变化点，提取变化点的位移。

下面对本实施例的几个主要部分进行解释说明。

1.1通过GNSS动态定位获得坐标时间序列。

坐标时间序列GNSS动态定位获得坐标组成，预测更新根据前一历元的状态估值来预测当前历元的状态估值，测量更新根据新的测量值对状态估值进行更新，在GNSS动态定位中将状态估值中的坐标分量的预测更新权重设置很小以确保定位结果的动态特性。动态定位方式包括相对实时定位，精密单点实时定位等实时定位模式。

1.2去噪处理。

对于以上通过GNSS动态定位构造的坐标时间序列处理单元，如果突变信号不明显，选择一定的方法对处理单元进行去噪，目的是削弱噪声影响，突出突变信号。这里以小波分析为例，选择合适的小波函数,进行小波变换提取信号中的突变成分，这里选择Daubechies小波函数。其中，db1小波的正则系数与信号在突变处的Lipschitz 指数具有最大的相似性，可以最有效地刻画突变。确定小波母函数之后，采用离散小波变换对原始信号进行分解。

其中，x(t)是需要分析的坐标时间序列；m和n分别是缩放和平移因子。ψ(·)是所选择的小波母函数，ψ^*(·)是小波母函数的共轭。a₀是指定的固定步长参数，通常大于1；b₀是大于0的位置参数。对于选择的db1小波母函数，一般将分解层级设置为 5-7层，将最终的分解的近似分量参数进行小波重构。

1.3窗口分割探测可能变化点。

(1)基于时间序列分割的位移检测。

时间序列表示为x₁,x₂,…x_n。对于经典的时间序列分割问题，假设时间序列中存在m个变化点，则将发生时间记录为τ＝[τ₁,τ₂,…τ_m]。整个时间序列通过这些m个变化点被分成m+1段。各段的统计平均值表示如下：

其中，τ₀和τ_m+1分别定义为τ₀＝1，τ_m+1＝n+1。

对于使用GNSS进行的滑坡短期位移检测，更多的是关注变化发生的时间和程度，而不是对整个时间序列进行分段。此外，多变化点检测通常是计算密集型的，不适合在线实现。因此，如图2所示的检测窗口被定义用于短期位移检测。被τ分割的两个部分的平均值和方差分别表示为h₀，h₁和ρ² ₀，ρ² ₁。变化点前后的位移估计为：

Δh＝h₀-h₁

在定义了这样一个检测窗口之后，需要找到最可能的变化点(PCP)。因此，下面介绍检测指标。

(2)检测窗口分割指标。

对于所定义的检测窗口的分割指标不局限于本发明提到的方法。在时间序列分割中，根据各分割部分的统计信息构造目标函数，最常用的统计值是方差：

其中，J(τ)为目标函数。

对于位移探测问题，目标函数简化为：

其中，τ为检测窗口的分割点，x_i、x_j表示检测窗口中点的坐标值，n表示检测窗口的长度。

可能变化点的探测问题变化为：对于任何检测窗口，都可以找到一个τ使下面的公式最小化。

其中，

表示满足目标函数最小化的最佳估值。

下面需要进一步根据可能变化点的统计值来进一步确定变化点。

1.4判定变化点。

对于一个历元，在一个大小为n的检测窗口中，该历元数据参与窗口分割的次数是从进入检测窗口开始到完全离开检测窗口结束，共n次。在窗口分割过程中，如果一个历元的数据被划分为可能变化点的次数越多，就越有可能成为变化点，因此，将统计值count_i定义为历元i被划分为可能变化点的次数。该统计值大于某一阈值的历元被视为变化点(CP)。

1.5提取变化点的位移。

对于count_i次被划分为变化点的第i个历元，存在count_i次的位移估值，表示为：

Δh_ij＝h_ij0-h_ij1,1≤j≤count_i

其中，Δh_ij是第i个历元对应的第j次分割的位移估值，h_ij0表示第i个历元对应的第j次分割的前半窗口中的坐标均值，h_ij1表示第i个历元对应的第j次分割的后半窗口中的坐标均值。

由于可能分割点对应于多个分割窗口，因此不可能使用全局指标来评估这些位移的可靠性。此处认为位移估值越集中，越接近真实值。在此基础上，根据以上count_i次位移估计范围划分位移区间，采用如下公式：

bin_ik＝div(min(Δh_ij),max(Δh_ij),s),1≤k≤s

其中，bin_ik表示第i个历元对应的第k个位移区间，div(.)操作符表示根据count_i个估值中的最大值和最小值均匀划分为指定的份数，s表示份数。例如，可设计为100 份，即s取100。

然后，基于落入位移估值最多的位移区间获得最终位移估计值，表示为：

其中，Δh_i表示最终位移估计值，count(bin_ik)表示Δh_ij落入bin_ik的数量，f为从位移区间到位移的映射函数，此实施例中采用位移区间的中值。

为了验证本发明的效果，下面结合实验结果进行说明。

根据短期滑坡位移变化特征，设计了场地实验，接收机天线在一个小时内发生三次阶跃变动。根据上述1.1中计算的动态定位结果，以及根据上述1.2进行小波重构结果如图3所示。可以看出经过小波重构之后，噪声得到了有效滤除，时间序列中阶跃成分特征明显。根据上述1.3计算的可能变化点如图4所示，在突变点附近比较集中。小波重构对应的可能变化点相对集中。对于窗口大小设置为360的检测窗口，根据上述1.4统计的所有的可能变化点的统计计数如图5所示，三次阶跃的统计计数分别为339,353,290，明显大于附近可能变化点的统计计数。对于图4中的可能变化点，根据上述1.5中介绍的方法提取位移如图6所示，提取的位移精度可以达到亚厘米级别。

最后所应说明的是，以上具体实施方式仅用以说明本发明的技术方案而非限制，尽管参照实例对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明技术方案的精神和范围，其均应涵盖在本发明的权利要求范围当中。

Claims (8)

1.基于GNSS动态定位时间序列分割的短期位移探测方法，其特征在于，包括以下步骤：

步骤1、通过GNSS动态定位获得原始坐标时间序列；

步骤2、利用检测窗口进行时间序列分割，获得可能变化点的集合；

在时间序列分割中，根据被τ分割的两个部分的统计信息构造目标函数，所述目标函数表示为：

其中，J(τ)表示目标函数，τ表示检测窗口的分割点，x_i、x_j表示检测窗口中点的坐标值，n表示检测窗口的长度；

通过使所述目标函数最小化获得可能变化点，表示为：

其中，

表示满足目标函数最小化的最佳估值；

步骤3、针对每个可能变化点，获得可能变化点的统计值，根据统计值及预设阈值判断可能变化点是否为变化点；

步骤4、针对每个变化点，获得变化点对应的历次位移估计值，基于所述历次位移估计值获得最终位移估计值。

2.根据权利要求1所述的基于GNSS动态定位时间序列分割的短期位移探测方法，其特征在于，所述步骤1中还包括：对所述原始坐标时间序列进行去噪处理，得到去噪后的坐标时间序列；所述步骤2中，利用检测窗口对所述去噪后的坐标时间序列进行分割，获得可能变化点的集合。

3.根据权利要求1所述的基于GNSS动态定位时间序列分割的短期位移探测方法，其特征在于，所述步骤1中，所述GNSS动态定位为相对实时定位或精密单点实时定位。

4.根据权利要求2所述的基于GNSS动态定位时间序列分割的短期位移探测方法，其特征在于，采用小波分析方法对所述原始坐标时间序列进行分解和重构，得到重构的坐标时间序列，所述重构的坐标时间序列作为所述去噪后的坐标时间序列。

5.根据权利要求4所述的基于GNSS动态定位时间序列分割的短期位移探测方法，其特征在于，所述小波分析方法包括以下子步骤：

确定小波母函数；

采用离散小波变换对所述原始坐标时间序列进行分解，得到分解的近似分量参数；

基于所述近似分量参数进行小波重构，得到所述重构的坐标时间序列；

其中，对所述原始坐标时间序列进行分解表示为：

其中，x(t)表示原始坐标时间序列，m表示缩放因子，n表示平移因子，ψ(·)表示小波母函数，ψ^*(·)表示小波母函数的共轭，a₀表示指定的固定步长参数，b₀表示大于0的位置参数。

6.根据权利要求5所述的基于GNSS动态定位时间序列分割的短期位移探测方法，其特征在于，所述小波母函数采用Daubechies db1小波函数。

7.根据权利要求1所述的基于GNSS动态定位时间序列分割的短期位移探测方法，其特征在于，所述步骤3中，对于一个历元，在一个大小为n的检测窗口中，该历元参与窗口分割的次数为n次；在窗口分割过程中，将第i个历元被划分为可能变化点的次数记为count_i，count_i作为该可能变化点的统计值；若该可能变化点的统计值大于预设阈值，则将该历元视为变化点。

8.根据权利要求1所述的基于GNSS动态定位时间序列分割的短期位移探测方法，其特征在于，所述步骤4中，对于count_i次被划分为可能变化点的第i个历元，存在count_i次的位移估值，表示为：

Δh_ij＝h_ij0-h_ij1,1≤j≤count_i

其中，Δh_ij表示第i个历元对应的第j次分割的位移估值，count_i表示第i个历元被划分为可能变化点的次数，h_ij0表示第i个历元对应的第j次分割的前半窗口中的坐标均值，h_ij1表示第i个历元对应的第j次分割的后半窗口中的坐标均值；

基于count_i次的位移估值划分位移区域，采用如下公式：

bin_ik＝div(min(Δh_ij),max(Δh_ij),s),1≤k≤s

其中，bin_ik表示第i个历元对应的第k个位移区间，div(.)操作符表示根据count_i次的位移估值中的最大值和最小值均匀划分为指定的份数，s表示份数；

基于落入位移估值最多的位移区间，获得最终位移估计值，表示为：

其中，Δh_i表示最终位移估计值，count(bin_ik)表示Δh_ij落入bin_ik的数量,f为从位移区间到位移的映射函数。

Images