一种基于GPS与强震仪观测自适应组合的形变监测方法

摘要

本发明提供了一种基于GPS与强震仪观测自适应组合的形变监测方法，将高分辨率的强震仪加速度观测融入到GPS‑PPP定位模型中，并将强震仪的基线漂移误差当作未知参数与其它定位参数一起进行实时估计；同时通过自适应调整动态噪声和求解时间窗口来实现参数估计的最优化，既高效解决了基线漂移的校正问题，同时又增强了GPS的求解强度，降低了高频噪声，实现了两种观测技术的优势互补，可以实时获取高精度宽频带的形变信息，具有广泛的应用前景。本发明具有算法先进，精度高、稳定性强，自动化程度高等特点，并兼具BDS、GPS、GLONASS等多卫星系统的整体封装，实现简单，可适用于多种应用场合。

说明书

技术领域

本发明涉及一种地表形变监测方法。

背景技术

GPS和强震仪观测是获取高精度地表形变(位移、速度、加速度)的两种有效手段，它们已广泛应用于自然灾害监测且各具特色。GPS易于获取高精度位移，但存在采样频率低、高频信噪比低、信号稳定性差的缺陷；同时，强震仪易于获取高分辨率加速度，但因基线漂移误差的存在，其积分后的速度和位移常存在偏差。当前的数据处理方式大多是单传感器模式，导致多传感器的观测资源没有充分利用。如何将两类观测数据进行有机组合，实现优势互补，实时提供高精度宽频带的形变信息，更好的服务于灾害等形变监测具有重要价值。

发明内容

为了克服现有技术的不足，本发明提供一种基于GPS与强震仪观测自适应组合的形变监测方法，将高分辨率的强震仪加速度观测融入到GPS-PPP定位模型中，并将强震仪的基线漂移误差当作未知参数与其它定位参数一起进行实时估计，同时通过自适应调整动态噪声和求解时间窗口来实现参数估计的最优化，既高效解决了基线漂移的校正问题，同时又增强了GPS的求解强度，降低了高频噪声，实现了两种观测技术的优势互补，可以实时获取高精度宽频带的形变信息。

本发明解决其技术问题所采用的技术方案是：通过两种观测数据的融合解算得到的基本的结果信息后，自适应的调整求解策略，实现参数估计最优化和两种观测手段的优势互补，具体包括以下步骤：

第一步，获得GPS的相位/伪距观测值、强震仪的加速度观测值、以及GPS数据处理需要的辅助产品，包括GPS卫星轨道、卫星钟差和地球自转参数；

第二步，对第一步获得的GPS的相位/伪距观测值、强震仪的加速度观测值、以及GPS数据处理需要的辅助产品进行数据检查、粗差剔除和周跳探测，得到干净的数据；对所述干净的数据进行相对论、潮汐、天线相位中心、对流层和地球自转误差的修正；

第三步，建立GPS与强震仪观测的紧组合模型，所述GPS和强震仪观测的紧组合模型如下：

${\left(\begin{array}{c}{L}_c\\ P_c\end{array}\right)}_k=\left(\begin{array}{cccc}{e}_k& m_{z,k}& 1& 1\\ e_k& m_{z,k}& 1& 0\end{array}\right)\left(\begin{array}{c}{x}_k\\ z_k\\ {\mathrm{\delta t}}_k\\ b\end{array}\right)+\left(\begin{array}{c}{ϵ}_{L_c}\\ {ϵ}_{P_c}\end{array}\right),\begin{array}{c}{ϵ}_{L_c}~N(0,{\sigma }_{L_c}^2)\\ {ϵ}_{P_c}~N(0,{\sigma }_{P_c}^2)\end{array}$

$\left(\begin{array}{cc}{a}_k={\stackrel{··}{x}}_k-u_k+{ϵ}_{a_k}& {ϵ}_{a_k}~N(0,{\sigma }_{a_k}^2)\end{array}\right)$

其中，L_c和P_c是扣除卫星和接收机间几何距离后的载波相位残差和伪距码观测值残差；e是卫星与接收机之间的单位矢量；角标k表示历元数；m_z,k是对流层湿延迟天顶方向投影函数；x、分别是接收机的位移和加速度矢量；z、δt和b分别为对流层天顶延迟、接收机钟差和相位模糊度；ε是测量噪声，ε_Lc和ε_Pc分别表示载波和伪距观测的噪声，其对应的方差为δ²；a和u分别为强震仪加速度观测值和基线漂移误差；基线漂移误差参数在设定宽度的每个估计窗口内作为常数估计，所述的天顶对流层延迟作为常数进行估计或者表现为随机游走过程，所述的接收机钟差作为高斯白噪声逐历元估计，所述的相位模糊度在连续无周跳情况下作为常数估计；

接收机的运动和基线漂移误差的状态方程如下：

$\left(\begin{array}{c}{x}_k\\ {\stackrel{·}{x}}_k\\ {\stackrel{··}{x}}_k\\ u_k\end{array}\right)=\left(\begin{array}{cccc}I& 0.5I\tau & 0.5{\mathrm{I\tau }}^2& 0\\ 0& I& I\tau & 0\\ 0& 0& I& 0\\ 0& 0& 0& I\end{array}\right)\left(\begin{array}{c}{x}_{k-1}\\ {\stackrel{·}{x}}_{k-1}\\ a_k-{\hat{u}}_{k-1}\\ u_{k-1}\end{array}\right)+w_k$

其中，为接收机的速度矢量，I是3x3的单位矩阵，τ是GPS和强震仪的最小采样率，w_k是动态噪声，其数学期望E(w_k)＝0，方差q_u为基线漂移的功率密度，k-1表示上一历元；采用基线漂移改正后加速度来代替加速度预测值a_k-1；给定观测值权比和动态噪声；采用卡尔曼滤波进行参数估计，同时进行残差编辑和精度统计，得到接收机的位置、速度、加速度、基线漂移、对流层延迟和接收机钟差信息及精度指标；

第四步，根据第三步的求解信息计算基线漂移的时间序列，并构建四个学习统计量，确定统计量的标准；四个统计量分别为：

基线漂移的位移影响

原始加速度的位移影响

改正后的加速度的位移影响

信号强度

其中，t₁、t₂是定义的估计窗口的起始时间和终止时间，分别为起始时间和终止时间对应的历元数，是a的均值，每个学习统计量的标准A₀，B₀，C₀，D₀是静态观测期间每个学习统计量在设定时间序列中的三倍标准差；

第五步，根据第四步确定的标准将基线漂移误差参数分成四个阶段，并计算新的动态噪声和时间估计窗口；四个阶段的确定如下：

初始化阶段：满足(A₂&B₁)或(B₁&C₂)，

静态阶段：满足(A₁)，

瞬时基线漂移阶段：满足(A₂&D₂)，

永久基线漂移阶段：满足(A₂&D₁)&(B₂&C₁)；

在瞬时基线漂移阶段，动态噪声调整如下：

$q_u^{\prime }=\sqrt{\frac{|u_{t_2}-u_{t_1}|}{(t_2-t_1)}}$

其中，q′_u为基线漂移的新的功率密度；基线漂移估计窗口的长度改变为每秒一次；对于其它三个阶段，动态噪声q_u＝0.0001-0.001m^1/2/s^3/2，窗口长度保持为5秒；

第六步，根据第五步中确定的新的动态噪声和新的时间估计窗口，采用卡尔曼滤波进行自适应组合解算，获得最优组合结果。

本发明的有益效果是：

第一，实现两种观测手段的优势互补，提供高精度宽频带的形变信息。

本发明将高精度GPS和高分辨率强震仪观测进行紧组合处理，弥补了GPS技术信噪比差和强震仪技术基线漂移的局限性，实现了他们的优势互补，其组合结果可以实时输出高精度宽频带的位移、速度、加速度信息，可直接应用于灾害等形变监测。

第二，有效识别并校正了强震仪的基线漂移，为基线漂移的校正提供了新途径。

强震仪的基线漂移误差的识别是校正是地震领域棘手的问题，通常采用经验的方法进行完成，但经验的方法不能准确描述基线漂移误差的特性，校正后的结果存在不同程度的系统偏差，且不能实时进行。本发明将强震仪的基线漂移误差进行模型化参数估计，既可以有效的识别并校正强震仪的基线漂移，同时为基线漂移的校正提供了新途径，对旋转地震学的研究具有重要意义。

第三，有效的降低了GPS的观测噪声，加快了GPS定位的收敛速度。

高分辨率的加速度观测融合到GPS定位方程中，可以约束GPS求解强度，加快GPS定位的收敛速度，同时可以有效的降低GPS的噪声，提高GPS定位精度和稳定性。

附图说明

图1是紧组合处理中的四个阶段示意图；

图2是GPS与强震仪观测组合框图；

图3是GPS与强震仪自适应组合框图。

具体实施方式

下面结合附图和实施例对本发明进一步说明，本发明包括但不仅限于下述实施例。

本发明将高分辨率的强震仪加速度观测融入到GPS-PPP定位模型中，并将强震仪的基线漂移误差当作未知参数与其它定位参数一起进行实时估计，同时通过自适应调整动态噪声和求解时间窗口来实现参数估计的最优化，既高效解决了基线漂移的校正问题，同时又增强了GPS的求解强度，降低了高频噪声，实现了两种观测技术的优势互补，可以实时获取高精度宽频带的形变信息。

本发明的技术方案主要包括两大核心算法：

(1)GPS与强震仪观测的紧组合处理

GPS和强震仪观测的紧组合采用下列观测方程：

${\left(\begin{array}{c}{L}_c\\ P_c\end{array}\right)}_k=\left(\begin{array}{cccc}{e}_k& m_{z,k}& 1& 1\\ e_k& m_{z,k}& 1& 0\end{array}\right)\left(\begin{array}{c}{x}_k\\ z_k\\ {\mathrm{\delta t}}_k\\ b\end{array}\right)+\left(\begin{array}{c}{ϵ}_{L_c}\\ {ϵ}_{P_c}\end{array}\right),\begin{array}{c}{ϵ}_{L_c}~N(0,{\sigma }_{L_c}^2)\\ {ϵ}_{P_c}~N(0,{\sigma }_{P_c}^2)\end{array}---\left(1\right)$

$\begin{array}{cc}{a}_k={\stackrel{··}{x}}_k-u_k+{ϵ}_{a_k}& {ϵ}_{a_k}~N(0,{\sigma }_{a_k}^2)\end{array}---\left(2\right)$

其中，L_c和P_c是扣除卫星和接收机间几何距离后的载波相位和伪距码观测值残差，e是卫星与接收机之间的单位矢量，k表示历元。m是对流程湿延迟投影函数，x，分别是接收机位移、加速度矢量，z，δt和b分别为对流层天顶延迟、接收机钟差和相位模糊度，ε是测量噪声，其对应的方差为δ²。a和u分别为强震仪加速度观测值和基线漂移误差，k表示历元数。基线漂移参数在每个估计窗口中作为常数估计，其它的参数，如天顶对流层延迟，在短时间内(一个小时或两个小时)可以作为常数进行估计，或者表现为随机游走过程，而接收机钟差当作高斯白噪声逐历元估计，模糊度参数在连续无周跳情况下作为常数估计。

测站运动和基线漂移的状态方程表示如下：

$\left(\begin{array}{c}{x}_k\\ {\stackrel{·}{x}}_k\\ {\stackrel{··}{x}}_k\\ u_k\end{array}\right)=\left(\begin{array}{cccc}I& 0.5I\tau & 0.5{\mathrm{I\tau }}^2& 0\\ 0& I& I\tau & 0\\ 0& 0& I& 0\\ 0& 0& 0& I\end{array}\right)\left(\begin{array}{c}{x}_{k-1}\\ {\stackrel{·}{x}}_{k-1}\\ a_k-{\hat{u}}_{k-1}\\ u_{k-1}\end{array}\right)+w_k,---\left(3\right)$

其中为速度矢量，I是3x3的单位矩阵，τ是GPS和强震仪的最小采样率，w_k是动态噪声，其数学期望E(w_k)＝0，方差q_u为基线漂移的功率密度。k-1表示上一历元，这里采用基线漂移改正后加速度来代替加速度预测值a_k-1，避免加速度和基线漂移具有大的动态噪声。应用公式(1)，(2)及(3)，给定观测值权比和动态噪声(这里经验权重由测量噪声来确定，通常的取值范围分别为：动态噪声给定为常数，q_u＝0.003～0.01m^1/2/s^3/2，窗口时间长度为1-10秒)，可用卡尔曼滤波进行参数估计既可以得到位置，速度，加速度、基线漂移、对流层延迟和接收机钟差。另外，上述随机模型可自适应地调整使得参数估计最优化。

(2)GPS与强震仪观测自适应组合处理

对于GPS和强震仪观测的紧组合的稳健估计处理，需要解决两个关键问题。第一个是基线漂移是否发生？如果发生，则需要估算基线漂移的持续时间。第二个是基线漂移的变化情况，于是本发明可以确定状态向量的最佳动态噪声。

为了解决这些问题，常采用自适应滤波处理。首先定义一些学习的统计量来检测基线漂移是否发生，并确定它是否显著，然后调整动态噪声来实现参数求解的自适应估计。在卫星导航领域中，这些统计量通常是基于状态差异、方差分量和预测残差来定义。但它们不能很好的应用于GPS和强震紧组合处理。一方面，对于GPS和强震仪紧组合处理，引入了基线漂移和相位模糊参数，单历元的有限观测量估计的几何状态参数偏差较大甚至不能估计。在这种情况下，状态差异或方差分量的统计特性会很差。另一方面，强震记录的基线漂移对积分后速度位移非常敏感。虽然基线漂移的值相比峰值加速度的值通常是非常小的，但它们将显著影响积分后的速度和位移。此外，基线漂移误差主要是由地面倾斜或旋转造成，并有一个重要特征，在没有运动时为零或恒定值，因此，选择基线漂移作为学习统计量更为合适。再者，高分辨率的加速度信号也可以用于检测运动状态。这里先假设静态状态下基线漂移和加速度为零，然后定义了四个学习统计量，来识别是否存在基线漂移以及基线漂移变化情况。由于强震仪采样率比GPS高得多，选择GPS采样间隔持续的时间作为基线漂移估计窗口的长度。同时，为避免随机误差的影响，在整个基线漂移估计窗口将基线漂移和加速度积分成位移进行分析。

统计量1:

统计量2:

统计量3:

统计量4:

其中t₁，t₂是定义的估计时间窗口(一般定义为采样间隔的倍数)，为历元数，是a的均值，A代表基线漂移的位移影响，B代表原始加速度的位移影响，C代表改正后的加速度的位移影响，D代表信号强度。每个学习统计量根据定义的标准A₀，B₀，C₀，D₀都有两种情况。这些标准是在静态观测期间，通过每个学习统计量(即A，B，C，D)的时间序列(3-5分钟)的三倍标准差计算得到。

基于所定义的统计量和标准，可以分析基线漂移估值并自适应识别。它们可以分成四个阶段，如图1所示。

初始化阶段(1)：满足(A₂&B₁)或(B₁&C₂)，图1.b中的0-47秒。数据解算的开始。静止状态，原始加速度几乎为零(如图1.a)，但解算的基线漂移不是零(如图1.b)。解算基线漂移需要时间来收敛，完成初始化。基线漂移的动态噪声应该是非常小的值。

静态阶段(2)：满足(A₁)，图1.b中的47-95秒。在此期间，解算的基线漂移几乎为零，基线漂移的动态噪声也非常小。

瞬时基线漂移阶段(3)：满足(A₂&D₂)，图1.b中的95-150秒。在此期间，检测到基线漂移，并且信号显著；基线漂移的动态噪声应调整为较大的值。

永久基线漂移阶段(4)：满足(A₂&D₁)&(B₂&C₁)，图1.b中的150秒至最后。在此期间，也检测到基线漂移，但没有显著信号，基线漂移的动态噪声应该是较小的值。

在检测到瞬时基线漂移的时间间隔内，动态噪声调整如下：

$q_u^{\prime }=\sqrt{\frac{|u_{t_2}-u_{t_1}|}{(t_2-t_1)}}---\left(8\right)$

其中q′_u为基线漂移的新的功率密度，反映了加速度基线漂移的实际变化率。另外，基线漂移估计窗口的长度可以根据信号存在情况自适应的调整：如果不存在信号，如情况D₁时，基线漂移主要是由环境变化引起，并且变化缓慢，所以本发明可以在几秒或更长的时间间隔来估计基线漂移；否则，如情况D₂，基线漂移主要是由地面运动的倾斜和/或旋转引起，且变化非常迅速，因此估计的频率应该为一秒或更少的时间间隔，来发现基线漂移的变化并调整动态噪声，本发明采用的是一秒钟估计一次。

本发明的实施例包括以下步骤：

第一步，GPS与强震仪观测的紧组合处理，如图2所示。分为数据输入，数据处理和结果输出三个核心模块。在数据输入模块，需要提供GPS的相位/伪距观测值，强震仪的加速度观测值，以及GPS数据处理需要的辅助产品如轨道、钟差、地球自转参数产品。在数据解算模块，首先是数据的预处理，进行数据检查，粗差剔除和周跳探测得到干净的数据；随后，进行相对论、潮汐、天线相位中心、对流层、地球自转误差的修正工作；接着进行紧组合模型的建立，将强震仪的基线漂移误差当作未知参数和其它定位参数一起估计；最后进行残差编辑和精度统计得到组合解(即基本组合解)。

第二步，GPS与强震仪自适应组合处理，如图3所示。自适应组合处理采用两个卡尔曼滤波器来实现。第一个滤波器用于基本组合解算，即是第一步的工作，用来判断运动状态和确定最佳的动态噪声以及基线漂移估计窗口的长度。第二个滤波器用于自适应组合解算，以获得稳健的组合结果。由经验权重和动态噪声(这里经验权重由测量噪声来确定，通常的取值范围分别为：比如动态噪声给定为常数，比如q_u＝0.003-0.01m^1/2/s^3/2，窗口时间长度为1-10秒)，首先采取第一个滤波器来估计基线漂移，然后计算和确定其时间序列以及四种学习统计数据标准。通过分析，若判定为瞬时基线漂移阶段，动态噪声根据瞬时基线漂移的变化进行改变，基线漂移估计窗口的长度改变为每秒一次。对于其它状态期间，动态噪声限制为一个较小的值(取值范围为q_u＝0.0001-0.001m^1/2/s^3/2)，窗口长度保持为5秒。然后，根据确定的最佳动态噪声和时间估计窗口，使用第二个滤波器进行自适应组合解算，以获得最优组合结果。