首页> 中国专利> 利用单频GNSS接收机的实时高精度地震形变监测方法

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

摘要

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

著录项

  • 公开/公告号CN103293550A

    专利类型发明专利

  • 公开/公告日2013-09-11

    原文格式PDF

  • 申请/专利权人 武汉大学;

    申请/专利号CN201310193930.9

  • 申请日2013-05-23

  • 分类号G01V1/28(20060101);G01S19/14(20100101);

  • 代理机构武汉科皓知识产权代理事务所(特殊普通合伙);

  • 代理人严彦

  • 地址 430072 湖北省武汉市武昌区珞珈山武汉大学

  • 入库时间 2024-02-19 20:39:13

法律信息

  • 法律状态公告日

    法律状态信息

    法律状态

  • 2015-09-16

    授权

    授权

  • 2013-10-16

    实质审查的生效 IPC(主分类):G01V1/28 申请日:20130523

    实质审查的生效

  • 2013-09-11

    公开

    公开

说明书

技术领域

本发明涉及一种新的实时高精度地震形变监测方法,尤其涉及一种利用单频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个历元时刻tk的距离变率通过下式进行中心差分计算,

其中,λ1为单频GNSS接收机的载波波长,分别为第k-1,k+1个历元的GNSS 相位观测值,Δt为相邻历元间隔。

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

(1)获取距离观测值ρ的微分

ρ·=[(Xr-XS)(X·r-X·S)+(Yr-YS)(Y·r-Y·S)+(Zr-ZS)(Z·r-Z·S)]+c(dtr·-dt·s)+d·ρion+d·trop+ϵ

其中,c为光速,ρ为距离观测值,Xr,Yr,Zr为测站坐标,XS,YS,ZS为卫星坐标, dtr,dtS分别为接收机和卫星钟差,dion,dtrop分别为电离层和对流层延迟,ε为观测噪声, 上标点表示微分形式;

(2)利用GNSS导航星历计算卫星坐标(XS,YS,ZS),并基于伪距单点定位确定测站 的坐标(Xr,Yr,Zr);

(3)利用GNSS导航星历计算卫星速度

(4)利用下式计算卫星钟速相对论改正值改正卫星钟速,

δt·S=-4.443×10-10eacosEdE/dt

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

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

V1V2···Vi···Vk=Bx1By1Bz11Bx2By2Bz21············BxiByiBzi1············BxKByKBzK1X·rY·rZ·rcdt·r-L1L2···Li···LK

其中,

Li=ρ·-[Bxi(X·r0-X·iS)+Byi(Y·r0-Y·iS)+Bzi(Z·r0-Z·iS)+c(dt·r-dt·iS)]

Bxi=(Xr-XiS)/ρi

Byi=(Yr-YiS)/ρi

Bzi=(Zr-ZiS)/ρ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广播星历和距离变率进行单点定速求tk时刻 接收机速度(即测站速度),然后通过长期的解算测站速度结果所得速度时间序列判断是否发 生形变,并进行数据存储后返回对下一历元进行同样处理。具体实施时,也可选用单频GNSS 接收机接收到的tk历元多普勒观测值计算距离变率和精度,计算实现方式为现有技术。

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

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个历元时刻tk的距离变率,具体可以用下式表示:

λ1为L1波段的波长,分别为k-1,k+1个历元L1相位观测值(单位为周),Δt为 相邻历元间隔,为第k个历元时刻tk的距离变率。

3单点定速算法

根据基本GNSS观测方程:

ρ=(Xr-XS)2+(Yr-YS)2+(Zr-ZS)2+c(dtr-dtS)(6) +dion+dtrop+ϵ

式(6)中,c为光速,ρ为距离观测值,Xr,Yr,Zr为测站坐标,XS,YS,ZS为卫星坐 标,dtr,rtS分别为接收机和卫星钟差,dion,dtrop分别为电离层和对流层延迟,ε为观测噪 声。对式(6)两边求微分,用上标点表示微分形式,则可得到:

ρ·=[(Xr-XS)(X·r-X·S)+(Yr-YS)(Y·r-Y·S)+(Zr-ZS)(Z·r-Z·S)]/ρ+c(dt·r-dt·S)+d·ion---(7)+d·trop+ϵ

对式(7)进行变形整理,监测站速度V和监测站速度观测值残差L可表示为:

V=(Xr-XS)ρX·r+(Yr-YS)ρY·r+(Zr-ZS)ρZ·r+cdt·r-L

L=ρ·-[(Xr-XS)ρ(X·r0-X·S)+(Yr-YS)ρ(Y·r0-Y·S)  (8)  +(Zr-ZS)ρ(Z·r0-Z·S)+cdt·r0-cdt·S+d·ion+d·trop]

式(7)、(8)即为单点测速中的基本观测方程。式中,为卫地距离变率,为卫 星速度,分别为接收机和卫星钟速,分别为接收机速度和钟速的 初始值,为电离层和对流层延迟变率。

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

3.1距离观测值ρ的微分获取

可通过上述中心差分法求得,即采用式(5)计算得到。

3.2卫星轨道和接收机坐标

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

3.3卫星速度计算

利用导航星历计算卫星速度其精度优于1mm/s。

3.4卫星钟速

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

δt·S=-4.433×10-10eacosEdE/dt---(9)

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

3.5电离层对流层延迟变率

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

3.6计算测站速度

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

Bxi=(Xr-XiS)/ρi

Byi=(Yr-YiS)/ρi

Bzi=(Zr-ZiS)/ρi---(10)

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

Li=ρ·-[Bxi(X·r0-X·iS)+Byi(Y·r0-Y·iS)+Bzi(Z·r0-Z·iS)+c(dt·r-dt·iS)]

式(10)中,上/下标i表示第i颗卫星。即ρi为第i颗卫星的距离观测值,为卫地距离 变率,为第i颗卫星的卫星坐标,为第i颗卫星的卫星速度,为 第i颗卫星的卫星钟速。将当前历元所有K颗可见且无周跳的观测卫星观测方程列为矩阵形 式,并假设第i颗卫星的观测值残差为Vi,i的取值为1,2,…K,则可得观测方程:

V1V2···Vi···Vk=Bx1By1Bz11Bx2By2Bz21············BxiByiBzi1············BxKByKBzK1X·rY·rZ·rcdt·r-L1L2···Li···LK---(11)

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

利用最小二乘算法对式(11)求解即可得到当前历元时刻测站的速度与 接收机钟速通过各历元的测站速度序列即可判断是否发生形变。

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

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

去获取专利,查看全文>

相似文献

  • 专利
  • 中文文献
  • 外文文献
获取专利

客服邮箱:kefu@zhangqiaokeyan.com

京公网安备:11010802029741号 ICP备案号:京ICP备15016152号-6 六维联合信息科技 (北京) 有限公司©版权所有
  • 客服微信

  • 服务号