CN116973974A - 基于北斗PPP-B2b服务的实时同震位移速度解算方法及系统 - Google Patents

基于北斗PPP-B2b服务的实时同震位移速度解算方法及系统 Download PDF

Info

Publication number
CN116973974A
CN116973974A CN202310699165.1A CN202310699165A CN116973974A CN 116973974 A CN116973974 A CN 116973974A CN 202310699165 A CN202310699165 A CN 202310699165A CN 116973974 A CN116973974 A CN 116973974A
Authority
CN
China
Prior art keywords
ppp
satellite
correction
beidou
observation
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202310699165.1A
Other languages
English (en)
Inventor
周庐艳
聂志喜
王振杰
张远帆
徐晓飞
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China University of Petroleum East China
Original Assignee
China University of Petroleum East China
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by China University of Petroleum East China filed Critical China University of Petroleum East China
Priority to CN202310699165.1A priority Critical patent/CN116973974A/zh
Publication of CN116973974A publication Critical patent/CN116973974A/zh
Pending legal-status Critical Current

Links

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/30—Analysis
    • G01V1/307—Analysis for determining seismic attributes, e.g. amplitude, instantaneous phase or frequency, reflection strength or polarity
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/13—Receivers
    • G01S19/14—Receivers specially adapted for specific applications
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/13—Receivers
    • G01S19/23—Testing, monitoring, correcting or calibrating of receiver elements
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/13—Receivers
    • G01S19/24—Acquisition or tracking or demodulation of signals transmitted by the system
    • G01S19/27—Acquisition or tracking or demodulation of signals transmitted by the system creating, predicting or correcting ephemeris or almanac data within the receiver
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/13—Receivers
    • G01S19/24—Acquisition or tracking or demodulation of signals transmitted by the system
    • G01S19/29—Acquisition or tracking or demodulation of signals transmitted by the system carrier including Doppler, related
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/13—Receivers
    • G01S19/24—Acquisition or tracking or demodulation of signals transmitted by the system
    • G01S19/30—Acquisition or tracking or demodulation of signals transmitted by the system code related
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/13—Receivers
    • G01S19/35—Constructional details or hardware or software details of the signal processing chain
    • G01S19/37—Hardware or software details of the signal processing chain
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/63—Seismic attributes, e.g. amplitude, polarity, instant phase

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Signal Processing (AREA)
  • Geology (AREA)
  • Geophysics (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • Environmental & Geological Engineering (AREA)
  • Acoustics & Sound (AREA)
  • Power Engineering (AREA)
  • Position Fixing By Use Of Radio Waves (AREA)

Abstract

本发明属于信息技术服务技术领域，尤其涉及一种基于北斗PPP‑B2b服务的实时同震位移速度解算方法及系统，利用PPP‑B2b服务播发的轨道、钟差及码偏差改正数，实时解算同震位移及速度，获取北斗三号GEO卫星播发的PPP‑B2b电文数据，结合广播星历数据，恢复基于PPP‑B2b改正数的精密钟差及轨道；获取GNSS观测数据，接收机将PPP‑B2b精密钟差与轨道应用于消电离层组合的PPP观测方程，根据卫星高度角进行定权，构建随机模型。本发明基于北斗PPP‑B2b服务，结合常加速度动力学模型，为强震近场无网络通讯的GNSS地震监测台站提供了一种稳定可靠、价格低廉的实时位移、速度及加速度计算方法，为地震监测及预警研究提供了观测数据，可广泛应用于卫星测量领域。

Description

基于北斗PPP-B2b服务的实时同震位移速度解算方法及系统
技术领域
本发明属于信息技术服务技术领域，尤其涉及一种基于北斗PPP-B2b服务的实时同震位移速度解算方法及系统。
背景技术
目前，随着GNSS(Global Navigation Satellite System)的快速发展，GNSS高精度定位技术已经广泛应用于地震监测领域。其获取的高精度的GNSS地表同震位移序列是进行震源破裂过程反演的前提条件，可以为地震应急响应提供技术支撑。
GNSS高精度同震位移提取方法主要有相对定位和绝对定位两种方式。实时动态相对定位(Real-Time Kinematic，RTK)方法需要一个或多个的参考站以去除公共误差并且可以恢复双差模糊度的特性，实现高精度的厘米级定位。然而，当参考站处于震区时，会因其发生的同震现象而干扰流动站的位移序列，并且伴随着GNSS基线产长度的增加，卫星星历及大气等误差的相关性逐渐减弱，定位精度会收到影响。
精密单点定位方法(Precise Point Positioning，PPP)仅利用单台接收机就可以获取高精度的绝对同震位移序列。但是，PPP与PPP-AR方法因模糊度固定的影响，其初始化时间较长，如果地震恰巧发生在模糊度收敛时期，其获取的同震位移序列的精度会大大影响。实时单站测速方法(Variomertic Approach Displacement Analysis Stand-aloneEngine，VADASE)采用单台接收机，基于历元间差分的载波相位观测值进行速度估计。其历元间做差的方式可以有效削弱模糊度和大气残余误差的影响，并且不需要初始化时间；但是，VADASE方法需要采用积分方式获取位移值，其位移值存在较大的漂移，并且这种漂移需要采用额外的线性滤波或者空间滤波的方式进行消除。时域点定位(Temporal PointPositioning，TPP)方法不再将历元间作差的方式局限于相邻历元，避免了积分过程，能过有效削弱历元间差分的残余误差的累积效应，并且没有初始化时间。但其需要高精度的接收机位置，当数据不连续时，需要将周跳探测与修复后的值叠加至时域点差分的GNSS观测值上。相较于TPP方法而言，精密单点定位速度估计方法(PPP Velocity Estimation，PPP-VE)及加速度估计方法(PPP-Acceleration Estimation，PPP-AE)可以直接估计测站的速度和加速度矢量，进而通过积分方法获取同震位移序列，可以有效克服上述GNSS同震位移提取方法的缺点。国际GNSS服务(International GNSS Service，IGS)从2013年开始通过互联网提供实时服务(Real-Time Service，RTS)。通过采用RTS产品，用户可以实现实时PPP应用。然而，这项服务依赖于一个固定的互联网接入。在实际地震监测过程中，通讯基站可能会因为地表的强运动而被破坏，IGS RTS精密星历产品将难以通过网络进行传输，届时将难以获取高精度的同震位移序列。
除了通过互联网通信提供的服务，卫星广播服务也已经被一些商业公司发布了高成本费用，但在大规模布设的GNSS地震监测站中配备卫星通讯设备成本较高，难以实时应用。而北斗三号系统(BDS-3)的PPP-B2b服务是一种免费开源的服务，是北斗系统的重要创新，该服务可以直接通过地球静止轨道卫星(Geosynchronous Eearth Orbit，GEO)卫星播发精密轨道及钟差改正数，可以免受网络的影响，可在强震网络中断时，实时解算高精度的位移及速度序列，应用于地震监测及预警领域。
通过上述分析，现有获取同震位移速度的技术有RTK，RTS PPP、星站差分等商业技术，存在的问题及缺陷为：RTK方法中基准站受到地震影响，导致测站解算结果误差增大，无法获取客观位移速度信息；RTS PPP方法依赖网络，强震发生时网络存在中断风险，将无法进行同震位移解算；星站差分等商业技术价格昂贵，在地震监测中无法实时应用。
发明内容
针对现有技术存在的问题，本发明提供了一种基于北斗PPP-B2b服务的实时同震位移速度解算方法及系统。
本发明是这样实现的，一种基于北斗PPP-B2b服务的实时同震位移速度解算方法，实现步骤包括：
步骤一，结合广播星历，完成基于北斗三代PPP-B2b服务的精密星历计算方法，主要包括精密轨道及钟差的计算；
步骤二，构建实时同震位移解算方法，主要包括GNSS观测数据的预处理与误差改正，PPP观测及随机模型的建立，常加速度动力学模型的构建，监测台站的三维位移、速度及加速度的估计方法；
进一步，所述步骤一中，为了保证不同信息类型所播发信息内容之间的关联性，以及改正数信息与广播星历之间的关联性，通过空间状态表达式数据期号(IOD SSR)、伪随机噪声码数据期号(IODP)、导航数据期号(IODN)和改正数据识别与数据期号(IOD Corr)四个版本号对信息进行了标识，方便匹配使用；
PPP-B2b电文的轨道改正数中轨道改正数的历元时间间隔为48s，播发的轨道改正信息包括轨道改正向量。利用轨道改正值可以计算出卫星位置改正向量。轨道改正信息的IODN与广插星历导航电文中的IODE成功匹配后，即可将卫星轨道坐标系下的精密轨道改正数表述为地心地固坐标系下的精密轨道改正数：
其中，[δOr δOa δOc]T为径向、切向和法向星固(Along-Cross-Radial，ACR)坐标系下的轨道改正数；[δOx δOy δOz]T为卫星在地心地固(Earth Center Earth Fixed，ECEF)参考框架下的轨道改正数。[er ea ec]分别对应径向、切向和法向的单位矢量，可以表示为：
式中，r和为利用广播星历计算的卫星位置和卫星速度；
对基于PPP-B2b服务的精密轨道进行计算：
其中，[X Y Z]T为卫星在地心地固(Earth Center Earth Fixed，ECEF)参考框架下的坐标。
进一步，所述步骤一中，PPP-B2b钟差改正信息中钟差改正数的历元时间间隔为6s；钟差改正信息与广播星历成功匹配后，即可利用电文中包含的钟差改正参数C0，对广播星历计算得到的钟差参数进行改正，对基于PPP-B2b服务的精密钟差进行计算：
其中，由PPP-B2b信号恢复的精密钟差，C0为利用三阶多项式计算的钟差在视线方向的改正数；当PPP-B2b电文的轨道改正数和钟差改正信息中的IOD Corr相同，且它们的IODN能够和导航电文中的IODE相匹配时，轨道和钟差改正数可匹配使用，进行误差校正。
进一步，所述步骤二中，在PPP定位模型中，常使用的GNSS原始伪距及相位观测值可以表示为：
式中，Pi和Li为原始码/相位观测值；ρ表示接收机到卫星的几何距离；c代表光速；δtr和δts分别为接收机和卫星的时钟差；T为对流层时延；fi为频率值；I1为L1的电离层延迟；λi为波长，Ni为整周模糊度；br,i和分别表示接收机和卫星相位硬件偏差；Br,i和/>分别代表接收机和卫星伪距硬件偏差；/>和/>分别表示伪距和载波相位未建模误差，包括热噪声和多径误差。此外，未指出的相对论效应、Sagnac效应、Shapiro时延、相位缠绕和潮汐影响等误差已经通过现有模型进行了修正；
在PPP-VE定位模型中，采用双频消电离层组合消除一阶电离层延迟的影响，建立的双频消电离层组合观测值方程可以表示为：
式中，
采用卫星高度角随机模型进行定权，表达式为：
式中，w为解算的观测权矩阵，pj为每颗卫星的观测权因子，elj为卫星j的卫星高度角。
进一步，所述步骤二中，PPP-VE在进行滤波计算时的动力学模型一般基于牛顿第二定律，通常采用动态位置向量、速度和加速度对载体的运动轨迹进行描述；在实际解算中，常采用有限阶微分进行近似逼近，常加速度(Constant Acceleration，CA)动力学模型可以表示为：
式中，k,k-1分别表示当前时刻和上一时刻；τ为采样间隔；w表示状态噪声矩阵；qa为加速度的功率密度，qa设置为0.01m·s-5/2，I为三维单位矩阵。
由于采用了北斗PPP-B2b的精密轨道与精密卫星钟差，卫星的轨道误差和钟误差可认为已被消除，接收机和卫星伪距硬件偏差则被接收机钟差和卫星钟差吸收；对流层延迟可分为干延迟和湿延迟，其中天顶对流层干延迟一般用萨斯塔莫宁(Saastamonien)模型进行修正，而天顶对流层湿延迟分量则需要进行估计；将每个历元线性化后的GNSS伪距、载波相位对于CA模型下的观测方程表示为：
式中，
与/>为吸收了接收机和卫星伪距硬件偏差后的模糊度；x和e分别代表测站位置增量及其相应的系数矩阵；/>为吸收了接收机和卫星伪距硬件偏差后的接收机钟差；v和a分别表示接收机的三维速度和加速度；zwd为天顶湿延迟；Mw为对流层湿分量的投影函数；
可将上式简写为卡尔曼滤波中的观测方程：
Lk＝AkXk+νk,νk～N(0,Rk)；
式中，k分别表示当前时刻，L为观测值向量，A为观测方程的系数矩阵，ν为观测噪声矩阵，R为观测噪声方差阵，可用简要矩阵表达为：
式中，和/>分别表示载波相位和伪距的观测噪声。上述观测函数模型中，可将PPPVE的待估状态参数表示为：
进一步，所述步骤二中，PPP-VE模型通过状态转移矩阵将系统为速度与加速度与其他参数联系起来，其状态转移方程如下：
Xk＝Fk,k-1Xk-1+wk,wk～N(0,Qw)；
式中，k-1表示上一时刻，w为状态噪声矩阵；F为状态转移矩阵；Qw为状态噪声的方差阵；F和Qw的具体表达式如下：
式中，；mt和mz分别表示接收机钟和天顶湿延迟的功率密度，设置为100m2·s-1，/>设置为10-9m2·s-1；
最终，利用Kalman滤波方法进行解算观测历元，测站的同震位移可以通过积分由如下式获得：
式中，s为积分获取的同震位移值，k0为初始历元。
本发明的另一目的在于提供一种应用所述的基于北斗PPP-B2b服务的实时同震位移速度解算方法的基于北斗PPP-B2b服务的实时同震位移速度解算系统，基于北斗PPP-B2b服务的实时同震位移速度解算系统包括：
精密星历计算模块，用于结合广播星历，完成基于北斗三代PPP-B2b服务的精密星历计算方法；
位移解算方法构建模块，构建实时同震位移解算方法，包括GNSS观测数据的预处理与误差改正方法，PPP观测及随机模型的建立方法，常加速度动力学模型的构建方法，监测台站的三维位移、速度及加速度的估计方法。
本发明的另一目的在于提供一种计算机设备，计算机设备包括存储器和处理器，存储器存储有计算机程序，计算机程序被处理器执行时，使得处理器执行所述的基于北斗PPP-B2b服务的实时同震位移速度解算方法的步骤。
本发明的另一目的在于提供一种计算机可读存储介质，存储有计算机程序，计算机程序被处理器执行时，使得处理器执行所述的基于北斗PPP-B2b服务的实时同震位移速度解算方法的步骤。
本发明的另一目的在于提供一种信息数据处理终端，信息数据处理终端用于实现所述的基于北斗PPP-B2b服务的实时同震位移速度解算系统。
结合上述的技术方案和解决的技术问题，本发明所要保护的技术方案所具备的优点及积极效果为：
第一，本发明公开了一种基于北斗PPP-B2b服务的实时同震位移速度解算方法，利用PPP-B2b服务播发的轨道、钟差及码偏差改正数，实时解算同震位移及速度，解决地震发生时通讯不便及同震位移及速度解算精度低的问题，技术要点在于获取北斗三号GEO卫星播发的PPP-B2b电文数据，结合广播星历数据，恢复基于PPP-B2b改正数的精密钟差及轨道；获取GNSS观测数据，完成数据预处理与误差修正后，构建消电离层组合的PPP观测方程，接收机将PPP-B2b精密钟差与轨道应用于消电离层组合的PPP观测方程，根据卫星高度角进行定权，构建随机模型；根据所述PPP观测模型和卫星高度角随机模型，采用Kalman滤波进行参数估计，主要在滤波过程建立常加速度动力学模型完成对监测台站的三维位移、速度和加速度估计。
第二，本发明基于北斗PPP-B2b服务，结合常加速度动力学模型，为强震近场无网络通讯的GNSS地震监测台站提供了一种稳定可靠、价格低廉的实时位移、速度及加速度计算方法，为地震监测及预警研究提供了观测数据，可广泛应用于卫星测量领域。
第三，作为本发明的权利要求的创造性辅助证据，还体现在以下几个重要方面：
(1)本发明的技术方案转化后的预期收益和商业价值为：本发明可提供准确可靠的实时同震位移速度提取服务，未来可应用于各省地震局、应急搜救中心等机构，可产生显著的经济收益。
(2)本发明的技术方案填补了国内外业内技术空白：PPP-B2b技术自2020年发布以来，成为大地测量领域的研究热点，但是还未有在地震监测领域的应用。本发明结合PPP-B2b定位技术和VADASE方法，克服地震发生对定位精度的影响，成功将PPP-B2b技术应用于地震监测领域，填补了PPP-B2b技术在地震监测领域的技术空白。
(3)本发明的技术方案是否解决了人们一直渴望解决、但始终未能获得成功的技术难题：我国地处世界两大地震带，是全球高发区域之一，地震监测与预警对于防震减灾至关重要。但是目前的地震监测手段无法克服强震时的网络中断风险，以及地震发生时监测精度下降的缺陷，这导致无法提供准确可靠的地震监测与预警服务，地震应急与决策将遇到阻碍。本发明采用PPP-B2b技术获取实时同震位移与速度，不依赖网络通信，并结合VADASE方法提高解算精度，解决了行业两大痛点。
附图说明
为了更清楚地说明本发明实施例的技术方案，下面将对本发明实施例中所需要使用的附图做简单的介绍，显而易见地，下面所描述的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下还可以根据这些附图获得其他的附图。
图1是本发明实施例提供的实时同震位移及速度解算方法的流程图；
图2是本发明实施例提供的实际地震事件中同震位移估计结果与PPP位移结果对比图；
图3是本发明实施例提供的际地震事件中同震速度估计结果与VADASE测速结果对比图。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
针对现有技术存在的问题，本发明提供了一种基于北斗PPP-B2b服务的实时同震位移速度解算方法及系统，下面结合附图对本发明作详细的描述。
为了使本领域技术人员充分了解本发明如何具体实现，该部分是对权利要求技术方案进行展开说明的解释说明实施例。
实施例1
本发明利用北斗三号PPP-B2b服务实现实时同震位移速度解算的方法，包括实现地震近场实时PPP-VE，其中PPP-VE为精密单点定位-测速技术。
实现地震近场实时PPP-VE的步骤：
第一步，结合广播星历，完成基于北斗三代PPP-B2b服务的精密星历计算方法，主要包括精密轨道及钟差的计算，具体步骤如下：
1)将卫星轨道坐标系下的精密轨道改正数表述为地心地固坐标系下的精密轨道改正数：
其中，[δOr δOa δOc]T为径向、切向和法向星固(Along-Cross-Radial，ACR)坐标系下的轨道改正数；[δOx δOy δOz]T为卫星在地心地固(Earth Center Earth Fixed，ECEF)参考框架下的轨道改正数。
式中，r和为利用广播星历计算的卫星位置和卫星速度。
2)对基于PPP-B2b服务的精密轨道进行计算：
其中，[X Y Z]T为卫星在地心地固(Earth Center Earth Fixed，ECEF)参考框架下的坐标。
3)结合广播星历提供的卫星钟差，对基于PPP-B2b服务的精密钟差进行计算
其中，由PPP-B2b信号恢复的精密钟差，C0为利用三阶多项式计算的钟差在视线方向的改正数。
第二步，构建实时同震位移解算方法，主要包括GNSS观测数据的预处理与误差改正，PPP观测及随机模型的建立，常加速度动力学模型的构建，监测台站的三维位移、速度及加速度的估计方法，具体步骤如下：
1)在PPP定位模型中，常使用的GNSS原始伪距及相位观测值可以表示为
式中，Pi和Li为原始码/相位观测值；ρ表示接收机到卫星的几何距离；c代表光速；δtr和δts分别为接收机和卫星的时钟差；T为对流层时延；fi为频率值；I1为L1的电离层延迟；λi为波长，Ni为整周模糊度；br,i和分别表示接收机和卫星相位硬件偏差；Br,i和/>分别代表接收机和卫星伪距硬件偏差；/>和/>分别表示伪距和载波相位未建模误差，包括热噪声和多径误差。此外，未指出的相对论效应、Sagnac效应、Shapiro时延、相位缠绕和潮汐影响等误差已经通过现有模型进行了修正。
2)在PPP-VE定位模型中，采用双频消电离层组合消除一阶电离层延迟的影响，建立的双频消电离层组合观测值方程可以表示为：
式中，
3)PPP-VE在进行滤波计算时的动力学模型一般基于牛顿第二定律，通常采用动态位置向量、速度和加速度对载体的运动轨迹进行描述。在实际解算中，常采用有限阶微分进行近似逼近，常加速度(Constant Acceleration，CA)动力学模型可以表示为：
式中，k,k-1分别表示当前时刻和上一时刻；τ为采样间隔；w表示状态噪声矩阵；qa为加速度的功率密度，qa设置为0.01m·s-5/2，I为三维单位矩阵。
4)由于采用了北斗PPP-B2b的精密轨道与精密卫星钟差，卫星的轨道误差和钟误差可认为已被消除，接收机和卫星伪距硬件偏差则被接收机钟差和卫星钟差吸收。此外，对流层延迟可分为干延迟和湿延迟，其中天顶对流层干延迟一般用萨斯塔莫宁(Saastamonien)模型进行修正，而天顶对流层湿延迟分量则需要进行估计。因此，可以将每个历元线性化后的GNSS伪距、载波相位对于CA模型下的观测方程表示为：
式中， 与为吸收了接收机和卫星伪距硬件偏差后的模糊度；x和e分别代表测站位置增量及其相应的系数矩阵；/>为吸收了接收机和卫星伪距硬件偏差后的接收机钟差；v和a分别表示接收机的三维速度和加速度；zwd为天顶湿延迟；Mw为对流层湿分量的投影函数。
可将上式简写为卡尔曼滤波中的观测方程：
Lk＝AkXk+vk,vk～N(0,Rk)
式中，k分别表示当前时刻，L为观测值向量，A为观测方程的系数矩阵，ν为观测噪声矩阵，R为观测噪声方差阵，可用简要矩阵表达为：
式中，和/>分别表示载波相位和伪距的观测噪声。上述观测函数模型中，可将PPPVE的待估状态参数表示为：
PPP-VE模型通过状态转移矩阵将系统为速度与加速度与其他参数联系起来，其状态转移方程如下：
Xk＝Fk,k-1Xk-1+wk,wk～N(0,Qw)
式中，k-1表示上一时刻，w为状态噪声矩阵；F为状态转移矩阵；Qw为状态噪声的方差阵；F和Qw的具体表达式如下：
式中，mt和mz分别表示接收机钟和天顶湿延迟的功率密度，设置为100m2·s-1，设置为10-9m2·s-1
利用Kalman滤波方法进行解算观测历元，测站的同震位移可以通过积分由如下两式获得：
式中，s为积分获取的同震位移值，k0为初始历元。
与现有技术相比，本发明实施例所提供的利用北斗三号PPP-B2b服务实现实时同震位移及速度解算方法至少具有以下有益效果:
第一，本发明实施例所提供的利用北斗三号PPP-B2b服务实现实时同震位移及速度解算方法，可在强震网络中断时，实时解算高精度的位移及速度序列，应用于地震监测及预警领域。
第二，本发明实施例所提供的利用北斗三号PPP-B2b服务实现实时同震位移及速度解算方法中的PPP-B2b服务是免费公开的，不需要通过商用卫星通讯服务播发RTS服务即可实时解算同震位移及速度，降低了GNSS台站布设的负担。
第三，本发明实施例所提供的利用北斗三号PPP-B2b服务实现实时同震位移及速度解算方法，直接采用动力学模型估计地表的速度和加速度，几乎不存在VADASE方法解算的同震位移中出现的漂移问题。
实施例2
下面对基于北斗三号PPP-B2b服务的GNSS同震位移、速度及加速度解算方法的应用场景进行举例说明：
本发明基于北斗三代的PPP-B2b服务，构建了基于PPP-VE模型的实时同震位移、速度及加速度提取方法。为了评估基于PPP-B2b服务的实时位移、速度及加速度解算方法的性能，计算了基于广播星历的VADASE测速结果进行对比，以解算的基于GFZ中心发布的多系统最终星历产品(GBM)的后处理PPP位移结果作为参考值，探讨该发明的优势与劣势(表2)。
表2验证方案
本发明的应用实施例提供了一种计算机设备，计算机设备包括存储器和处理器，存储器存储有计算机程序，计算机程序被处理器执行时，使得处理器执行所述的基于北斗PPP-B2b服务的实时同震位移速度解算方法的步骤。
本发明的应用实施例提供了一种计算机可读存储介质，存储有计算机程序，计算机程序被处理器执行时，使得处理器执行所述的基于北斗PPP-B2b服务的实时同震位移速度解算方法的步骤。
本发明的应用实施例提供了一种信息数据处理终端，信息数据处理终端用于实现所述的基于北斗PPP-B2b服务的实时同震位移速度解算系统。
作为优选，本发明的应用实施例公开了一种基于北斗PPP-B2b服务的实时同震位移速度解算方法，利用PPP-B2b服务播发的轨道、钟差及码偏差改正数，实时解算同震位移及速度，解决地震发生时通讯不便及同震位移及速度解算精度低的问题，技术要点在于获取北斗三号GEO卫星播发的PPP-B2b电文数据，结合广播星历数据，恢复基于PPP-B2b改正数的精密钟差及轨道；获取GNSS观测数据，完成数据预处理与误差修正后，构建消电离层组合的PPP观测方程，接收机将PPP-B2b精密钟差与轨道应用于消电离层组合的PPP观测方程，根据卫星高度角进行定权，构建随机模型；根据所述PPP观测模型和卫星高度角随机模型，采用Kalman滤波进行参数估计，主要在滤波过程建立常加速度动力学模型完成对监测台站的三维位移、速度和加速度估计。本发明基于北斗PPP-B2b服务，结合常加速度动力学模型，为强震近场无网络通讯的GNSS地震监测台站提供了一种稳定可靠、价格低廉的实时位移、速度及加速度计算方法，为地震监测及预警研究提供了观测数据，可广泛应用于卫星测量领域。
本发明实施例在研发或者使用过程中取得了一些积极效果，和现有技术相比的确具备很大的优势，下面内容结合试验过程的数据、图表等进行描述。
1)静态实验
静态实验选取了分布在中国周边地区的10个IGS测站，实验数据选取了2022年3月17日(年积日76)的GPS和BDS双系统的双频数据，采样率为1s。解算时，将一天的GNSS数据每4小时划分为一个解算弧段。假设测站一直处于静止状态，以0作为参考值，对VADASE和PPP-VE方法求取的速度和位移定位评估。PPP-VE的收敛判断条件设为两连续历元之间的位移变化量小于0.01m，且并且一直保持在0.01米之内。由于PPPVE存在几分钟的收敛时间，因此选取每个解算弧段开始后15分钟后的速度和位移时序进行精度评估，评估时间长度为15分钟，历元数共900个。由于地震的主震破裂时间通常小于15分钟，因此，将每四小时解算弧段内的速度按照五分钟进行划分，每15分钟后重新开始积分位移。表2给出了VADASE方法和PPP-VE方法在南北、东西、垂直及3D四个维度的速度及位移结果的RMS值。
表2静态实验精度统计表
实验结果：从表中可以看出，VADASE方法与PPP-VE方法的测速精度相当，在南北、东西和竖直方向上均优于5mm/s，而PPP-VE方法相较于VADASE方法的位移提取精度有明显改善，在南北、东西和竖直方向分别提升了60％，70％和66％。
2)实际地震事件
据中国地震台网测定，北京时间2021年05月22日02时04分11秒(2021-05-2118:04:13UTC)，位于中国青海省果洛藏族自治州玛多县(北纬34.59度，东经98.34度)发生M7.4级地震，震源深度17千米。这次地震成功被43个青海连续运行参考站(QHCORS)和16个中国大陆构造环境监测网络(Crustal Movement Observation Network of China，CMONOC)测站记录下来。本次实验选取了青海玛多(QHMD)等58个测站的发震前后3小时的高频GNSS数据，并利用PPP-VE和VADASE方法提取了全部测站的同震位移及同震速度，并以后处理PPP结算的位移结果作为参考值，进行精度评估。图2和图3分别给出了典型QHMD测站的提取的同震位移及同震速度图像。可以看出，VADASE方法与PPP-VE方法的测速序列在整体趋势上基本一致，均能捕获到明显的地震信号，二者之间的差异值小于1mm/s，而PPP-VE方法相较于VADASE方法的位移提取精度有明显改善，并且没有出现明显的漂移现象，在南北、东西和竖直方向的位移提取精度分别为2.61cm，1.52cm和2.97cm。
应当注意，本发明的实施方式可以通过硬件、软件或者软件和硬件的结合来实现。硬件部分可以利用专用逻辑来实现；软件部分可以存储在存储器中，由适当的指令执行系统，例如微处理器或者专用设计硬件来执行。本领域的普通技术人员可以理解上述的设备和方法可以使用计算机可执行指令和/或包含在处理器控制代码中来实现，例如在诸如磁盘、CD或DVD-ROM的载体介质、诸如只读存储器(固件)的可编程的存储器或者诸如光学或电子信号载体的数据载体上提供了这样的代码。本发明的设备及其模块可以由诸如超大规模集成电路或门阵列、诸如逻辑芯片、晶体管等的半导体、或者诸如现场可编程门阵列、可编程逻辑设备等的可编程硬件设备的硬件电路实现，也可以用由各种类型的处理器执行的软件实现，也可以由上述硬件电路和软件的结合例如固件来实现。
以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，都应涵盖在本发明的保护范围之内。

Claims (10)

1.一种基于北斗PPP-B2b服务的实时同震位移速度解算方法，其特征在于，包括：
步骤一，结合广播星历，完成基于北斗三代PPP-B2b服务的精密星历计算方法，主要包括精密轨道及钟差的计算；
步骤二，构建实时同震位移解算方法，主要包括GNSS观测数据的预处理与误差改正，PPP观测及随机模型的建立，常加速度动力学模型的构建，监测台站的三维位移、速度及加速度的估计方法。
2.如权利要求1所述的基于北斗PPP-B2b服务的实时同震位移速度解算方法，其特征在于，所述步骤一中，为了保证不同信息类型所播发信息内容之间的关联性，以及改正数信息与广播星历之间的关联性，通过空间状态表达式数据期号(IOD SSR)、伪随机噪声码数据期号(IODP)、导航数据期号(IODN)和改正数据识别与数据期号(IOD Corr)四个版本号对信息进行了标识，方便匹配使用；
PPP-B2b电文的轨道改正数中轨道改正数的历元时间间隔为48s，播发的轨道改正信息包括轨道改正向量；利用轨道改正值可以计算出卫星位置改正向量；轨道改正信息的IODN与广插星历导航电文中的IODE成功匹配后，即可将卫星轨道坐标系下的精密轨道改正数表述为地心地固坐标系下的精密轨道改正数：
其中，[δOr δOa δOc]T为径向、切向和法向星固(Along-Cross-Radial，ACR)坐标系下的轨道改正数；[δOx δOy δOz]T为卫星在地心地固(Earth Center Earth Fixed，ECEF)参考框架下的轨道改正数；[er ea rc]分别对应径向、切向和法向的单位矢量，可以表示为：
式中，r和为利用广播星历计算的卫星位置和卫星速度；
对基于PPP-B2b服务的精密轨道进行计算：
其中，[X Y Z]T为卫星在地心地固(Earth Center Earth Fixed，ECEF)参考框架下的坐标。
3.如权利要求1所述的基于北斗PPP-B2b服务的实时同震位移速度解算方法，其特征在于，所述步骤一中，PPP-B2b钟差改正信息中钟差改正数的历元时间间隔为6s；钟差改正信息与广播星历成功匹配后，即可利用电文中包含的钟差改正参数C0，对广播星历计算得到的钟差参数进行改正，对基于PPP-B2b服务的精密钟差进行计算：
其中，由PPP-B2b信号恢复的精密钟差，C0为利用三阶多项式计算的钟差在视线方向的改正数；当PPP-B2b电文的轨道改正数和钟差改正信息中的IOD Corr相同，且它们的IODN能够和导航电文中的IODE相匹配时，轨道和钟差改正数可匹配使用，进行误差校正。
4.如权利要求1所述的基于北斗PPP-B2b服务的实时同震位移速度解算方法，其特征在于，所述步骤二中，在PPP定位模型中，常使用的GNSS原始伪距及相位观测值可以表示为：
式中，Pi和Li为原始码/相位观测值；ρ表示接收机到卫星的几何距离；c代表光速；δtr和δts分别为接收机和卫星的时钟差；T为对流层时延；fi为频率值；I1为L1的电离层延迟；λi为波长，Ni为整周模糊度；br,i和分别表示接收机和卫星相位硬件偏差；Br,i和/>分别代表接收机和卫星伪距硬件偏差；/>和/>分别表示伪距和载波相位未建模误差，包括热噪声和多径误差；此外，未指出的相对论效应、Sagnac效应、Shapiro时延、相位缠绕和潮汐影响等误差已经通过现有模型进行了修正；
在PPP-VE定位模型中，采用双频消电离层组合消除一阶电离层延迟的影响，建立的双频消电离层组合观测值方程可以表示为：
式中，
采用卫星高度角随机模型进行定权，表达式为：
式中，w为解算的观测权矩阵，pj为每颗卫星的观测权因子，elj为卫星j的卫星高度角。
5.如权利要求1所述的基于北斗PPP-B2b服务的实时同震位移速度解算方法，其特征在于，所述步骤二中，PPP-VE在进行滤波计算时的动力学模型一般基于牛顿第二定律，通常采用动态位置向量、速度和加速度对载体的运动轨迹进行描述；在实际解算中，常采用有限阶微分进行近似逼近，常加速度(Constant Acceleration，CA)动力学模型可以表示为：
式中，k,k-1分别表示当前时刻和上一时刻；τ为采样间隔；w表示状态噪声矩阵；qa为加速度的功率密度，qa设置为0.01m·s-5/2；
由于采用了北斗PPP-B2b的精密轨道与精密卫星钟差，卫星的轨道误差和钟误差可认为已被消除，接收机和卫星伪距硬件偏差则被接收机钟差和卫星钟差吸收；对流层延迟可分为干延迟和湿延迟，其中天顶对流层干延迟一般用萨斯塔莫宁(Saastamonien)模型进行修正，而天顶对流层湿延迟分量则需要进行估计；将每个历元线性化后的GNSS伪距、载波相位对于CA模型下的观测方程表示为：
式中， 与/>为吸收了接收机和卫星伪距硬件偏差后的模糊度；x和e分别代表测站位置增量及其相应的系数矩阵；/>为吸收了接收机和卫星伪距硬件偏差后的接收机钟差；v和a分别表示接收机的三维速度和加速度；zwd为天顶湿延迟；Mw为对流层湿分量的投影函数；
可将上式简写为卡尔曼滤波中的观测方程：
Lk＝AkXk+νk,νk～N(0,Rk)；
式中，k分别表示当前时刻，L为观测值向量，A为观测方程的系数矩阵，ν为观测噪声矩阵，R为观测噪声方差阵，可用简要矩阵表达为：
式中，和/>分别表示载波相位和伪距的观测噪声，上述观测函数模型中，可将PPPVE的待估状态参数表示为：
6.如权利要求1所述的基于北斗PPP-B2b服务的实时同震位移速度解算方法，其特征在于，所述步骤二中，PPP-VE模型通过状态转移矩阵将系统为速度与加速度与其他参数联系起来，其状态转移方程如下：
Xk＝Fk,k-1Xk-1+wk,wk～N(0,Qw)；
式中，k-1表示上一时刻，w为状态噪声矩阵；F为状态转移矩阵；Qw为状态噪声的方差阵；F和Qw的具体表达式如下：
式中，I为三维单位矩阵；mt和mz分别表示接收机钟和天顶湿延迟的功率密度，设置为100m2·s-1，/>设置为10-9m2·s-1；
最终，利用Kalman滤波方法进行解算观测历元，测站的同震位移可以通过积分由如下式获得：
式中，s为积分获取的同震位移值，k0为初始历元。
7.一种应用所述的基于北斗PPP-B2b服务的实时同震位移速度解算方法的基于北斗PPP-B2b服务的实时同震位移速度解算系统，其特征在于，基于北斗PPP-B2b服务的实时同震位移速度解算系统包括：
精密星历计算模块，用于结合广播星历，完成基于北斗三代PPP-B2b服务的精密星历计算方法；
位移解算方法构建模块，构建实时同震位移解算方法，包括GNSS观测数据的预处理与误差改正方法，PPP观测及随机模型的建立方法，常加速度动力学模型的构建方法，监测台站的三维位移、速度及加速度的估计方法。
8.一种计算机设备，计算机设备包括存储器和处理器，存储器存储有计算机程序，计算机程序被处理器执行时，使得处理器执行如权利要求1～6所述的基于北斗PPP-B2b服务的实时同震位移速度解算方法的步骤。
9.一种计算机可读存储介质，存储有计算机程序，计算机程序被处理器执行时，使得处理器执行如权利要求1～6所述的基于北斗PPP-B2b服务的实时同震位移速度解算方法的步骤。
10.一种信息数据处理终端，信息数据处理终端用于实现如权利要求7所述的基于北斗PPP-B2b服务的实时同震位移速度解算系统。
CN202310699165.1A 2023-06-14 2023-06-14 基于北斗PPP-B2b服务的实时同震位移速度解算方法及系统 Pending CN116973974A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202310699165.1A CN116973974A (zh) 2023-06-14 2023-06-14 基于北斗PPP-B2b服务的实时同震位移速度解算方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202310699165.1A CN116973974A (zh) 2023-06-14 2023-06-14 基于北斗PPP-B2b服务的实时同震位移速度解算方法及系统

Publications (1)

Publication Number Publication Date
CN116973974A true CN116973974A (zh) 2023-10-31

Family

ID=88480538

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202310699165.1A Pending CN116973974A (zh) 2023-06-14 2023-06-14 基于北斗PPP-B2b服务的实时同震位移速度解算方法及系统

Country Status (1)

Country Link
CN (1) CN116973974A (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117647308A (zh) * 2024-01-29 2024-03-05 武汉理工大学三亚科教创新园 基于单站北斗观测值的海洋平台振动监测方法
CN117784179A (zh) * 2024-02-23 2024-03-29 中国科学院空天信息创新研究院 基于PPP-B2b的实时空间环境感知监测系统及方法
CN117991307A (zh) * 2024-04-03 2024-05-07 江苏深蓝航天有限公司 导航接收机的位移的求解方法和装置

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117647308A (zh) * 2024-01-29 2024-03-05 武汉理工大学三亚科教创新园 基于单站北斗观测值的海洋平台振动监测方法
CN117784179A (zh) * 2024-02-23 2024-03-29 中国科学院空天信息创新研究院 基于PPP-B2b的实时空间环境感知监测系统及方法
CN117784179B (zh) * 2024-02-23 2024-04-30 中国科学院空天信息创新研究院 基于PPP-B2b的实时空间环境感知监测系统及方法
CN117991307A (zh) * 2024-04-03 2024-05-07 江苏深蓝航天有限公司 导航接收机的位移的求解方法和装置

Similar Documents

Publication Publication Date Title
JP7267460B2 (ja) 高完全性衛星測位のためのシステムおよび方法
CN109709591B (zh) 一种面向智能终端的gnss高精度定位方法
US10281587B2 (en) Navigation satellite system positioning involving the generation of correction information
US9405012B2 (en) Advanced global navigation satellite systems (GNSS) positioning using precise satellite information
US5917445A (en) GPS multipath detection method and system
CN116973974A (zh) 基于北斗PPP-B2b服务的实时同震位移速度解算方法及系统
Geng Rapid integer ambiguity resolution in GPS precise point positioning
US20180252819A1 (en) Method and system for performing precise point positioning (ppp) ambiguity resolution using gnss triple frequency signals
JP2010528320A (ja) リアルタイムキネマティック（ｒｔｋ）測位における距離依存性誤差の軽減
Li et al. Review of PPP–RTK: Achievements, challenges, and opportunities
Rabbou et al. Precise point positioning using multi-constellation GNSS observations for kinematic applications
Fernández-Hernández et al. Snapshot positioning without initial information
Li Improving real-time PPP ambiguity resolution with ionospheric characteristic consideration
CN113253314A (zh) 一种低轨卫星间时间同步方法及系统
Zheng et al. Capturing coseismic displacement in real time with mixed single-and dual-frequency receivers: application to the 2018 Mw7. 9 Alaska earthquake
Aggrey Multi-GNSS precise point positioning software architecture and analysis of GLONASS pseudorange biases
Ma et al. Flight-test evaluation of integer ambiguity resolution enabled PPP
Odijk et al. Two approaches to precise kinematic GPS positioning with miniaturized L1 receivers
CN116009042A (zh) 一种单站载波历元间差分实时探测相对形变的方法及系统
Bancroft Multiple inertial measurement unit fusion for pedestrian navigation
Rabah et al. Evaluation of the IGS–Global Ionospheric Mapping model over Egypt
Li et al. Toward Wide-Area and High-Precision Positioning With LEO Constellation Augmented PPP-RTK
CN112540389A (zh) 一种利用卫星历书的时间同步方法和装置
CN103389502B (zh) 基于多个地面基站高精度确定载体加速度的方法
Preston GPS Multipath Detection and Mitigation Timing Bias Techniques

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination