CN103293550A

CN103293550A - 利用单频gnss接收机的实时高精度地震形变监测方法

Info

Publication number: CN103293550A
Application number: CN2013101939309A
Authority: CN
Inventors: 李敏; 施闯; 李文文; 赵齐乐; 方荣新
Original assignee: Wuhan University WHU
Current assignee: Wuhan University WHU
Priority date: 2013-05-23
Filing date: 2013-05-23
Publication date: 2013-09-11
Anticipated expiration: 2033-05-23
Also published as: CN103293550B

Abstract

本发明涉及一种利用单频GNSS接收机的实时高精度地震形变监测方法，尤其涉及一种利用单频GNSS接收机采集的高采样率数据实时确定测站速度进行地震形变监测的方法，设置一台单频GNSS接收机、一个包含有数据接口的地震形变数据处理装置以及数据存储设备。上述设备中，GNSS接收机与地震形变数据处理装置通过数据接口相连接用于实时数据处理，而地震形变数据处理装置则与数据存储设备连接用于数据以及处理结果的保存。通过地震形变数据处理装置实时处理来自接收机的单频GNSS相位观测数据，解算并输出测站的速度用于判断测站是否发生形变。本发明可以广泛应用地震、泥石流等地形变形，或者大型工程结构的变形监测领域。

Description

利用单频GNSS接收机的实时高精度地震形变监测方法

技术领域

本发明涉及一种新的实时高精度地震形变监测方法，尤其涉及一种利用单频GNSS接收机采集的高采样率数据实时确定测站速度进行地震形变监测的方法。

背景技术

近年来，全球范围内地震、滑坡、地表沉降等地质灾害频频发生对地形以及大型结构的土木建筑等造成了不同程度的形变影响，从而对人类生命财产安全造成重大威胁。因此地形形变以及大型结构的形变监测和异常预警已成为亟待解决的技术问题。其技术难点包括：测量的实时性，测量的精确性等。全球卫星导航技术（GNSS）具有布设方便、成本低、精度稳定等优点。尤其是随着GNSS精密定位技术的成熟，主要包括RTK（Real-time kinematic）和精密单点定位PPP技术（Precise Point Positioning），GNSS得以逐渐应用于上述形变监测领域。

利用PPP技术进行地震监测时，高精度的卫星轨道和钟差产品是精密定位数据处理的重要部分。目前，国际GNSS服务中心（IGS）针对不同用户的需求，分别发布了事后、快速、超快速的精密星历和钟差产品。利用IGS发布的事后精密星历和钟差产品，进行事后PPP数据处理的定位精度已经能够达到mm级，基本上能够满足绝大部分的地形形变或地震监测的需求。但是，事后精密星历和钟差产品一般需要滞后12～18天才能获取，快速或者超快速产品一般只需要滞后1～3天就可以获取，因此其并不适用于时效性要求极高的形变监测领域；另一方面，基于PPP技术的形变监测，需要利用高精度双频接收机才能达到cm～mm量级的精度。双频接收机昂贵的价格成本，并需要高精度的实时导航卫星轨道和卫星钟差严重制约了该技术在高精度实时地震形变监测中的应用，特别是对新型导航系统尚不实用。如:北斗卫星导航系统目前尚缺乏cm级的高精度卫星轨道和卫星钟差。

GNSS系统导航电文是由导航卫星向用户播发的一组反映卫星在空间的运行轨道、卫星钟差的改正参数，电离层延迟修正参数和卫星工作状态等信息的二进制代码，其与C/A码，P码同时调制在载波上。因此，接收机锁定卫星并解出C/A码后，便可以实时获取广播星历。然而，目前GPS广播星历计算卫星轨道精度为1.6m，卫星钟差的精度为7ns，利用广播星历进行PPP定位无法满足cm级甚至mm级的变形监测的要求。

为了利用广播星历进行实时的定位数据处理并得到cm级的定位精度，可以采用RTK技术。然而由于RTK技术是采用基线解算，这就要求在测区范围内先建立若干参考站，并需要流动站与参考站之间时时保持通讯连接，这使得系统建设和运行维护成本大大增加，尤其是在在监测范围较大的情况下。

发明内容

本发明通过以下技术方案实现了实时高精度的地震形变监测。

本发明的技术方案为一种利用单频GNSS接收机的实时高精度地震形变监测方法。设置一台连接有GNSS天线的单频GNSS接收机、一个包含有数据接口的地震形变数据处理装置以及数据存储设备；GNSS接收机与地震形变数据处理装置通过数据接口相连接，地震形变数据处理装置与数据存储设备连接；

GNSS接收机采集GNSS观测数据，并解码生成GNSS导航星历和GNSS相位观测值，将GNSS导航星历和GNSS相位观测值输入地震形变数据处理装置用于后续数据处理；地震形变数据处理装置根据GNSS导航星历和GNSS相位观测值解算测站速度，根据测站速度的时间序列分析是否发生形变，将相关数据成果最终输入至数据存储设备中保存。

而且，根据GNSS导航星历和GNSS相位观测值解算测站速度实现方式为，根据GNSS相位观测值进行周跳探测，通过中心差分计算卫星到测站的距离变率，根据GNSS广播星历和距离变率进行单点定速解算测站速度。

而且，所述根据GNSS相位观测值进行周跳探测实现方式为，对每个卫星进行如下处理，

设取连续n个历元的GNSS相位观测值，首先计算判断因子δ，

其中，

为连续n个历元的GNSS相位观测值平均变换率，

为连续n个历元的GNSS相位观测值变换率平方和；

当δ大于等于预设阈值时，判断该卫星第n+1个历元有周跳，当δ小于预设阈值时，判断该卫星第n+1个历元没有发生周跳。

而且，所述通过中心差分计算卫星到测站的距离变率实现方式为，

设当前待求的是第k个历元时刻t_k的距离变率

通过下式进行中心差分计算，

其中，λ₁为单频GNSS接收机的载波波长，

分别为第k-1，k+1个历元的GNSS相位观测值，Δt为相邻历元间隔。

而且，所述根据GNSS广播星历和距离变率进行单点定速解算测站速度实现方式如下，

（1）获取距离观测值ρ的微分

\overset{\cdot}{ρ} = \frac{[(X_{r} - X^{S}) ({\overset{\cdot}{X}}_{r} - {\overset{\cdot}{X}}^{S}) + (Y_{r} - Y^{S}) ({\overset{\cdot}{Y}}_{r} - {\overset{\cdot}{Y}}^{S}) + (Z_{r} - Z^{S}) ({\overset{\cdot}{Z}}_{r} - {\overset{\cdot}{Z}}^{S})]}{+ c (\overset{\cdot}{{dt}_{r}} - {d \overset{\cdot}{t}}^{s}) + {\overset{ρ}{\overset{\cdot}{d}}}_{ion} + {\overset{\cdot}{d}}_{trop} + ϵ}

其中，c为光速，ρ为距离观测值，X_r，Y_r，Z_r为测站坐标，X^S，Y^S，Z^S为卫星坐标，dt_r，dt^S分别为接收机和卫星钟差，d_ion，d_trop分别为电离层和对流层延迟，ε为观测噪声，上标点表示微分形式；

（2）利用GNSS导航星历计算卫星坐标(X^S，Y^S，Z^S)，并基于伪距单点定位确定测站的坐标(X_r，Y_r，Z_r)；

（3）利用GNSS导航星历计算卫星速度

（4）利用下式计算卫星钟速相对论改正值

改正卫星钟速，

δ {\overset{\cdot}{t}}^{S} = - 4.443 \times 10^{- 10} e \sqrt{a} \cos EdE / dt

其中，e为卫星轨道偏心率，a为轨道长半径，E为偏近点角，dE/dt为偏近点角变化率；

（5）从所有可见卫星中，根据之前n个历元进行周跳探测时判断结果选择所有当前历元没有周跳的卫星，设共有K颗，K≥4；其中第i颗卫星的观测值残差为V_i，第i颗卫星计算的测站速度残差为L_i，i的取值为1,2，…K，第i颗卫星观测值在X、Y、Z三个方向的方向余弦分别用

表示，列观测方程如下，

[\begin{matrix} V_{1} \\ V_{2} \\ \cdot_{\cdot}^{\cdot} \\ V_{i} \\ \cdot_{\cdot}^{\cdot} \\ V_{k} \end{matrix}] = [\begin{matrix} B_{x}^{1} & B_{y}^{1} & B_{z}^{1} & 1 \\ B_{x}^{2} & B_{y}^{2} & B_{z}^{2} & 1 \\ \cdot_{\cdot}^{\cdot} & \cdot_{\cdot}^{\cdot} & \cdot_{\cdot}^{\cdot} & \cdot_{\cdot}^{\cdot} \\ B_{x}^{i} & B_{y}^{i} & B_{z}^{i} & 1 \\ \cdot_{\cdot}^{\cdot} & \cdot_{\cdot}^{\cdot} & \cdot_{\cdot}^{\cdot} & \cdot_{\cdot}^{\cdot} \\ B_{x}^{K} & B_{y}^{K} & B_{z}^{K} & 1 \end{matrix}] [\begin{matrix} {\overset{\cdot}{X}}_{r} \\ {\overset{\cdot}{Y}}_{r} \\ {\overset{\cdot}{Z}}_{r} \\ cd {\overset{\cdot}{t}}_{r} \end{matrix}] - [\begin{matrix} L_{1} \\ L_{2} \\ \cdot_{\cdot}^{\cdot} \\ L_{i} \\ \cdot_{\cdot}^{\cdot} \\ L_{K} \end{matrix}]

其中，

L_{i} = \overset{\cdot}{ρ} - [B_{x}^{i} ({\overset{\cdot}{X}}_{r}^{0} - {\overset{\cdot}{X}}_{i}^{S}) + B_{y}^{i} ({\overset{\cdot}{Y}}_{r}^{0} - {\overset{\cdot}{Y}}_{i}^{S}) + B_{z}^{i} ({\overset{\cdot}{Z}}_{r}^{0} - {\overset{\cdot}{Z}}_{i}^{S}) + c (d {\overset{\cdot}{t}}_{r} - d {\overset{\cdot}{t}}_{i}^{S})]

B_{x}^{i} = (X_{r} - X_{i}^{S}) / ρ_{i}

B_{y}^{i} = (Y_{r} - Y_{i}^{S}) / ρ_{i}

B_{z}^{i} = (Z_{r} - Z_{i}^{S}) / ρ_{i}

ρ_i为第i颗卫星的距离观测值，

为卫地距离变率，

为第i颗卫星的卫星坐标，

为第i颗卫星的卫星速度，

为第i颗卫星的卫星钟速，

分别为接收机速度和钟速的初始值；若当前历元在某颗卫星的GNSS相位观测值在根据之前n个历元进行周跳探测时判断结果为有周跳，列以上观测方程时不予考虑；

利用最小二乘算法对观测方程求解得到当前历元时刻测站的速度

与接收机钟速

本发明具有以下优点：

1.单频观测。本方法中直接采用单点定速的方式，利用实时获取的相位观测值和导航星历，实时解算站点的运动速度。因此本发明中不需要进行相位观测值的组合，可以直接使用单频的GNSS接收机，从而大大降低了系统成本，并扩大了应用范围。

2.单站实时监测。本发明避免了PPP技术需要高精度的星历和钟差产品无法实时处理的弊端，也避免了RTK技术需要进行组网观测的不便。利用本方法，可以实现任意单个观测站实时的形变监测。

3.方法简单，计算量小。本方法中，观测模型为单点定速模型，并进行逐历元求解。每个历元的待估参数仅仅为测站速度（3方向的分量）以及接收机钟速，采用最小二乘算法即可胜任。计算方法简单，计算量小。

4.解算精度高。本方法中虽然采用广播星历求解卫星轨道和速度，但根据研究，10m的卫星轨道误差对站星距离变率的影响最大为1.6mm/s，且目前利用广播星历计算卫星速度的精度为1mm/s。综合以上，利用本方法实时解算测站速度的精度在mm/s量级，能为后续的变形分析提供高精度的数据。

5.系统简单成本低，应用面广。利用本方法，仅仅需要一台单频GNSS接收机和天线，相应的数据处理模块以及数据存储设备，便可以解算得到实时的速度信息用于地震形变的监测分析。

附图说明

图1为本发明实施例的系统结构示意图。

图2为本发明实施例的系统运行示意图。

图3为本发明实施例的核心工作流程图。

具体实施方式

本发明提供一种利用单频GNSS接收机的实时高精度地震形变监测方法，其实现系统的组成部分仅仅包括一台单频GNSS接收机（含GNSS天线）、一个包含有数据接口的地震形变数据处理装置以及数据存储设备。上述设备中，GNSS天线连接于GNSS接收机上，GNSS接收机与地震形变数据处理装置通过数据接口相连接用于实时数据处理，而地震形变数据处理装置则与数据存储设备连接用于数据以及处理结果的保存。具体的系统结构图见图1。通过上述设备的组合，即可以实现实时高精度的地震形变监测，监测结果可送往远程的控制中心。

本发明中涉及的GNSS接收机及GNSS天线，用于采集和输出高频（大于1Hz）GNSS观测数据，其类型包含但不限于可接收美国的GPS，俄罗斯的GLONASS，欧盟Galileo，中国北斗（BDS），以及日本QZSS，印度IRNSS等卫星导航系统信号的接收机。本发明中的地震形变数据处理装置具体实施时可采用电脑实现，或集成于接收机芯片。地震形变数据处理装置是整个系统的核心，其与GNSS接收机和数据存储装置连接，并根据本发明的核心形变数据处理算法实时处理GNSS接收机输出的GNSS观测数据，快速判断是否发生形变，实时数据处理的结果最终存储于数据存储设备中。具体的处理流程图见图2。

本发明中的数据存储设备通过数据接口与地震形变数据处理装置连接，用于存储GNSS观测数据以及地震形变数据处理装置实时数据处理的结果，以便于事后的处理分析等。具体实施时可采用硬盘等作为数据存储设备。

如图2所示，本发明由地震形变数据处理装置执行的核心形变数据处理算法思想为由于在一般状态下，测站以及大型结构可以认为是处于静止状态，其速度为0。如果测站的速度发生突变，则可以认为是测站或大型结构发生了变形。单频GNSS接收机采集GNSS观测数据，并解码数据生成GNSS导航星历和高频的GNSS相位观测值后，地震形变数据处理装置利用实时获取的GNSS广播星历（即GNSS导航星历），以及通过高采样率的GNSS相位观测值中心差分计算所得卫星到测站的距离变率，可以通过单点定速算法逐历元实时求解测站速度。通过对测站速度的时间序列进行分析，最终可以判断是否发生形变，并可以准确定位发生形变的时间，为后续的分析提供参考依据。具体的算法流程见图3：通过GNSS相位观测值中心差分计算卫星到测站的距离变率，包括根据L1相位观测值周跳探测，采用无周跳相位数据利用中心差分法求距离变率。根据GNSS广播星历和距离变率进行单点定速求t_k时刻接收机速度（即测站速度），然后通过长期的解算测站速度结果所得速度时间序列判断是否发生形变，并进行数据存储后返回对下一历元进行同样处理。具体实施时，也可选用单频GNSS接收机接收到的t_k历元多普勒观测值计算距离变率和精度，计算实现方式为现有技术。

具体各个主要部分的算法包括：单频相位观测值周跳探测、中心差分求t_k时刻的距离变率和单点定速算法。其中单频相位观测值周跳探测是单频GNSS数据处理的必要前提，通过探测周跳，以保证距离变化率计算的正确性；中心差分求t_k时刻的距离变率是获取地震形变观测数据，是计算测站速度的基础；单点定速算法是整个地震形变数据处理的关键，通过计算各历元时刻速度反应是否发生地震形变。

1单频相位观测值周跳探测

由载波相位测量的特点可知：发生周跳前的载波相位观测值是随时间连续而平滑的函数，周跳后该函数值将发生跳变，但其变化率则是连续函数，且为载波相位的严格一阶导数。因此，可以利用载波相位变化率对单频载波相位观测值的进行周跳探测。

1.1求相邻载波相位变化率

Δt(e)为第f+1个历元和第f个历元之间的时间间隔，为第f+1个历元的L1相位观测值，

为第f个历元的L1相位观测值。

1.2求载波相位观测值的平均变化率

其中，设第f个参与计算的历元的载波相位变化率记为

n为参与计算的连续

的个数，即取连续n个历元的GNSS相位观测值数据，可根据经验进行设置。

1.3求载波相位观测值的变化率平方和

1.4计算判断因子δ

目前GNSS接收机量测的精度大约是一个码元长度的百分之一，而单频GPS接收机载波波长L1为0.119m，则载波相位测量精度以长度为单位来衡量为1.19mm，并考虑观测噪声的影响，因此将δ限差取为3mm。在实际操作中，假设前n个历元的载波相位观测值没有周跳，按上述公式(1)到(5)依次计算第n+1历元中每颗卫星的δ值，根据δ值与预设阈值的比较判断周跳。具体实施时，本领域技术人员可根据情况自行设定阈值，实施例的阈值设为3。当δ≥3时，说明该卫星第n+1个历元可能有周跳，应在最后式（11）所提供观测方程进行计算时舍弃第n+1历元该颗卫星的相位观测值。当δ<3时，说明此时所对应历元的载波相位观测值没有发生周跳，无需进行舍弃。在对第n+1历元所有卫星处理完毕之后，选取下一历元的观测数据后可以重复上述过程进行周跳探测，例如根据第2～n+1历元的载波相位观测值计算第n+2历元中每颗卫星的δ值…

2通过中心差分求取距离变率

距离变化率为计算观测站速度的观测值，而一般GNSS接收机仅输出伪距与相位观测值，本发明基于GNSS接收机输出的单频相位观测值通过中心差分计算距离变化率，设当前待求的是第k个历元时刻t_k的距离变率，具体可以用下式表示：

λ₁为L1波段的波长，

分别为k-1，k+1个历元L1相位观测值（单位为周），Δt为相邻历元间隔，

为第k个历元时刻t_k的距离变率。

3单点定速算法

根据基本GNSS观测方程：

ρ = \sqrt{{(X_{r} - X^{S})}^{2} + {(Y_{r} - Y^{S})}^{2} + {(Z_{r} - Z^{S})}^{2}} + c ({dt}_{r} - {dt}^{S})

(6)

+ d_{ion} + d_{trop} + ϵ

式（6）中，c为光速，ρ为距离观测值，X_r,Y_r，Z_r为测站坐标，X^S，Y^S，Z^S为卫星坐标，dt_r，rt^S分别为接收机和卫星钟差，d_ion，d_trop分别为电离层和对流层延迟，ε为观测噪声。对式（6）两边求微分，用上标点表示微分形式，则可得到：

\overset{\cdot}{ρ} = [(X_{r} - X^{S}) ({\overset{\cdot}{X}}_{r} - {\overset{\cdot}{X}}^{S}) + (Y_{r} - Y^{S}) ({\overset{\cdot}{Y}}_{r} - {\overset{\cdot}{Y}}^{S})

+ (Z_{r} - Z^{S}) ({\overset{\cdot}{Z}}_{r} - {\overset{\cdot}{Z}}^{S})] / ρ + c (d {\overset{\cdot}{t}}_{r} - d {\overset{\cdot}{t}}^{S}) + {\overset{\cdot}{d}}_{ion} - - - (7)

+ {\overset{\cdot}{d}}_{trop} + ϵ

对式（7）进行变形整理，监测站速度V和监测站速度观测值残差L可表示为：

V = \frac{(X_{r} - X^{S})}{ρ} {\overset{\cdot}{X}}_{r} + \frac{(Y_{r} - Y^{S})}{ρ} {\overset{\cdot}{Y}}_{r} + \frac{(Z_{r} - Z^{S})}{ρ} {\overset{\cdot}{Z}}_{r} + cd {\overset{\cdot}{t}}_{r} - L

L = \overset{\cdot}{ρ} - [\frac{(X_{r} - X^{S})}{ρ} ({\overset{\cdot}{X}}_{r}^{0} - {\overset{\cdot}{X}}^{S}) + \frac{(Y_{r} - Y^{S})}{ρ} ({\overset{\cdot}{Y}}_{r}^{0} - {\overset{\cdot}{Y}}^{S})

(8)

+ \frac{(Z_{r} - Z^{S})}{ρ} ({\overset{\cdot}{Z}}_{r}^{0} - {\overset{\cdot}{Z}}^{S}) + cd {\overset{\cdot}{t}}_{r}^{0} - cd {\overset{\cdot}{t}}^{S} + {\overset{\cdot}{d}}_{ion}

+ {\overset{\cdot}{d}}_{trop}]

式（7）、（8）即为单点测速中的基本观测方程。式中，

为卫地距离变率，

为卫星速度，

分别为接收机和卫星钟速，

分别为接收机速度和钟速的初始值，

为电离层和对流层延迟变率。

实施例进行单定定速处理时，具体实现如下：

3.1距离观测值ρ的微分

获取

可通过上述中心差分法求得，即采用式（5）计算得到。

3.2卫星轨道和接收机坐标

由式（7）中可见，卫星轨道和接收机坐标误差会导致站星方向余弦的计算误差。由于站星距数值大，算法对卫星和测站的位置小误差并不敏感，目前导航星历计算卫星位置的精度约1.6m，伪距单点定位精度也在5m左右，二者联合对测速精度的影响小于1mm/s。因此可利用广播星历计算卫星坐标(X^S，Y^S，Z^S)，并基于伪距单点定位确定测站的坐标(X_r，Y_r，Z_r)。

3.3卫星速度计算

利用导航星历计算卫星速度

其精度优于1mm/s。

3.4卫星钟速

GPS卫星配置的原子钟的稳定度约为10^-10～10^-12,其对站星距离变化率的影响量级为0.001ns/s。可以通过广播星历钟差参数进行改正，由于相对论效应对钟速影响可以达到0.01ns/s量级，因此需要利用下式进行改正，计算卫星钟速相对论改正值

δ {\overset{\cdot}{t}}^{S} = - 4.433 \times 10^{- 10} e \sqrt{a} \cos EdE / dt - - - (9)

式（9）中e为卫星轨道偏心率，a为轨道长半径，E为偏近点角，dE/dt为偏近点角变化率。

3.5电离层对流层延迟变率

由于电离层和对流层延迟变化缓慢，在高采样率时，可以忽略其影响。

3.6计算测站速度

综合以上，利用卫星广播星历计算卫星位置和速度以及卫星钟速，通过单点定位确定接收机坐标，并忽略电离层和对流层延迟率的影响，从所有可见卫星中，根据之前n个历元进行周跳探测时判断结果选择所有当前历元没有周跳的卫星，设共有K颗，K≥4。利用式（8）逐一列观测方程，假设第i颗卫星观测值在X、Y、Z三个方向的方向余弦分别用

表示，即

B_{x}^{i} = (X_{r} - X_{i}^{S}) / ρ_{i}

B_{y}^{i} = (Y_{r} - Y_{i}^{S}) / ρ_{i}

B_{z}^{i} = (Z_{r} - Z_{i}^{S}) / ρ_{i} - - - (10)

则第i颗卫星计算的测站速度残差L_i可表示为：

L_{i} = \overset{\cdot}{ρ} - [B_{x}^{i} ({\overset{\cdot}{X}}_{r}^{0} - {\overset{\cdot}{X}}_{i}^{S}) + B_{y}^{i} ({\overset{\cdot}{Y}}_{r}^{0} - {\overset{\cdot}{Y}}_{i}^{S}) + B_{z}^{i} ({\overset{\cdot}{Z}}_{r}^{0} - {\overset{\cdot}{Z}}_{i}^{S}) + c (d {\overset{\cdot}{t}}_{r} - d {\overset{\cdot}{t}}_{i}^{S})]

式（10）中，上/下标i表示第i颗卫星。即ρ_i为第i颗卫星的距离观测值，

为卫地距离变率，

为第i颗卫星的卫星坐标，为第i颗卫星的卫星速度，为第i颗卫星的卫星钟速。将当前历元所有K颗可见且无周跳的观测卫星观测方程列为矩阵形式，并假设第i颗卫星的观测值残差为V_i，i的取值为1,2，…K，则可得观测方程：

[\begin{matrix} V_{1} \\ V_{2} \\ \cdot_{\cdot}^{\cdot} \\ V_{i} \\ \cdot_{\cdot}^{\cdot} \\ V_{k} \end{matrix}] = [\begin{matrix} B_{x}^{1} & B_{y}^{1} & B_{z}^{1} & 1 \\ B_{x}^{2} & B_{y}^{2} & B_{z}^{2} & 1 \\ \cdot_{\cdot}^{\cdot} & \cdot_{\cdot}^{\cdot} & \cdot_{\cdot}^{\cdot} & \cdot_{\cdot}^{\cdot} \\ B_{x}^{i} & B_{y}^{i} & B_{z}^{i} & 1 \\ \cdot_{\cdot}^{\cdot} & \cdot_{\cdot}^{\cdot} & \cdot_{\cdot}^{\cdot} & \cdot_{\cdot}^{\cdot} \\ B_{x}^{K} & B_{y}^{K} & B_{z}^{K} & 1 \end{matrix}] [\begin{matrix} {\overset{\cdot}{X}}_{r} \\ {\overset{\cdot}{Y}}_{r} \\ {\overset{\cdot}{Z}}_{r} \\ cd {\overset{\cdot}{t}}_{r} \end{matrix}] - [\begin{matrix} L_{1} \\ L_{2} \\ \cdot_{\cdot}^{\cdot} \\ L_{i} \\ \cdot_{\cdot}^{\cdot} \\ L_{K} \end{matrix}] - - - (11)

这样当前历元在某颗卫星的载波相位观测值在根据之前n个历元进行周跳探测时判断结果为可能周跳，列以上观测方程时舍弃了该历元该颗卫星的相位观测值，可以避免影响结果精度。

利用最小二乘算法对式（11）求解即可得到当前历元时刻测站的速度

与接收机钟速

通过各历元的测站速度序列即可判断是否发生形变。

持续进行以上步骤，即可长期监测地震。

本文中所描述仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。