GNSS实时卫星钟差估计在地震监测中的应用

2023-07-08王浩浩庄文泉

导航定位与授时 2023年3期

王浩浩,郝明,庄文泉

(中国地震局第二监测中心,西安 710054)

0 引言

近年来,随着全球导航卫星系统(global navigation satellite system, GNSS)高频定位技术的不断发展,其在地壳运动监测领域发挥着越来越重要的作用,为地震监测提供了一种有效的新型技术途径[1-2]。其中,利用差分相对定位技术可以计算得到监测站相对于某一固定参考站毫米级的动态位移,但监测区域内有时因地质条件等客观因素难以布设观测环境良好、稳定的基准站,而且坐标精确已知的基准站在强震发生阶段可能会发生移动,导致动态解算的定位精度显著降低[3-4]。而采用国际GNSS服务(international GNSS service, IGS)组织发布的精密卫星轨道和钟差等产品的精密单点定位(precise point positioning, PPP)技术,具有不依赖于某一特定参考基准站、实时性强等优势,仅利用单台GNSS接收机即可获得国际地球参考框架下高精度的“绝对位置”,具备准确捕捉地震位移波形的能力,更适合长距离、大范围的地壳形变监测[3-4]。文献[5]实验结果表明,卫星钟差的采样率越高,利用PPP动态解捕捉远场站点形变信息的优势越明显。文献[6]利用武汉大学自主研发的PANDA软件对高频GNSS观测信号进行PPP后处理,能够很好地获取2016-11-13新西兰Mw7.8地震产生的动态位移特征。

地震瞬时地表动态位移的实时高可靠性监测对地震预警系统而言至关重要,能够为震中以及震级的快速确定等研究工作提供关键信息[3]。然而,实时PPP高精度动态解的实现取决于精密卫星轨道和钟差的实时估计性能。目前超快速实时预报轨道产品已经能够满足其时效性和可靠性等要求,但卫星钟差的精确预报极易受到自身时频特性以及复杂太空环境的影响,高可靠的厘米级GNSS实时卫星钟差则需利用地面站观测数据进行实时估计得到[7]。文献[8]实时估计的多模GNSS卫星钟差与武汉大学最终精密钟差互差优于0.2 ns。文献[9-10]利用均方根信息滤波实现GNSS四系统实时卫星钟差的快速估计,且与最终精密钟差产品符合性较好。文献[11]提出了一种基于序贯最小二乘的在线质量控制方法,GNSS四系统实时估计卫星钟差的标准差均值优于0.1 ns。

鉴于此,本文以长安大学北斗分析中心为平台,基于多模全球导航卫星系统实验跟踪网(multi-GNSS experiment, MGEX)监测站1 Hz的GNSS观测数据,采用模型严密、精度较高的非差估计法[12],实现了多模GNSS卫星钟差的实时估计和性能评估,并将其用于高频动态PPP实时获取2021年漾濞Mw6.4地震和玛多Mw7.4地震发生时段的地表形变波形,具体分析震时点位的运动变化情况。

1 数据处理

1.1 钟差估计

本文采用无电离层组合观测值进行GNSS卫星钟差的实时估计,非差伪距和载波相位的观测方程为

(1)

由于公式中的偏差参数具有很强甚至完全的线性相关性,将其作为未知数进行估计,不仅会引入大量待估参数,还会减弱卫星钟差的估计精度。其中,卫星端的码偏差参数具有较高的时间稳定性,可以在观测方程中将其合并到卫星的钟误差参数中,而且相位延迟偏差可以被相应的相位模糊度参数吸收。进而可得观测模型

(2)

对于全球定位系统(global positioning system, GPS)、格洛纳斯卫星导航系统(global navigation satellite system, GLONASS)、伽利略卫星导航系统(Galileo satellite navigation system, Galileo)以及我国的北斗卫星导航系统(BeiDou navigation satellite system, BDS)而言,不同卫星导航系统所采用的时空基准以及信号体制不一致,导致卫星信号在GNSS接收机内部产生系统间偏差(inter system bias, ISB)以及GLONASS所特有的频率间偏差(inter frequency bias, IFB)[11]。若以GPS时间系统为基准,则顾及ISB/IFB参数的GNSS卫星钟差实时估计模型为

(3)

1.2 处理策略

无电离层组合模型是GNSS精密卫星钟差估计、精密单点定位等数据处理常用的函数模型,表1总结了GNSS卫星钟差实时估计以及实时动态PPP基于该模型的数据处理策略。对于GNSS卫星钟差的实时估计,则是在观测方程中将卫星钟误差作为白噪声进行估计以免受到钟跳的影响,将卫星轨道和测站坐标作为已知值进行改正。对于GNSS实时动态PPP,则是将测站的位置坐标作为待估参数进行估计,将精密卫星轨道和卫星钟差作为已知值进行改正。

表1 数据处理策略Tab.1 Data processing strategy

2 多模GNSS实时卫星钟差估计分析

在MGEX中选取60个均匀分布在全球的连续运行跟踪站。从IGS下载2021年5月19日—2021年5月22日连续4 d所选测站1 s采样间隔的GNSS观测数据用于卫星钟差的实时估计。在实时估计多模GNSS卫星钟差的过程中,测站始终处于静止状态且位置已知,可以将测站坐标固定到IGS周解。现阶段,超快速实时预报轨道产品与最终精密轨道产品的精度量级相当,卫星轨道可以固定为武汉大学提供的包含G/R/E/C四系统轨道的6 h超快速解,以减少待估参数的个数。同时,选择某一测站的接收机钟差作为参考基准钟进行约束[15],基于双频非差消电离层组合的载波相位和伪距观测值,采用序贯最小二乘平差和验后残差质量控制算法,最终估计得到历元间隔为5 s的多模GNSS实时卫星钟差。

2.1 估计时间分析

当利用60个测站构成的地面跟踪站网进行多模GNSS卫星钟差实时估计时,每个历元大约需要处理3 000个观测值并对将近2 000维的矩阵进行求逆。每个历元需要处理的观测值数目和法方程维数如图1所示。为加快高维矩阵的运算速度,提高钟差估计算法的计算效率,本文采用LAPACK函数库中相关的矩阵运算解算法方程,以实现多模GNSS实时卫星钟差的快速估计[16]。图2统计了LAPACK方法下,实时估计2021年年积日第141天G/R/E/C四系统卫星钟差时每个历元的计算时间。可以得出,在全球分布60个连续运行跟踪站的情况下,所有历元的钟差估计耗时均小于5 s,单历元钟差估计的平均计算时间为3.4 s,可见上述算法对高频观测数据的处理效率足以实时估计5 s历元间隔的多模GNSS卫星钟差。

图1 每个历元的观测值数目和法方程维数Fig.1 The number of observations and the dimension of the normal equation system at each epoch

图2 多模GNSS实时钟差单历元估计时间Fig.2 The computation time for the multi-GNSS real-time clock offset estimation at each epoch

2.2 钟差精度分析

采用二次差的方法计算实时估计卫星钟差与武汉大学最终精密卫星钟差的差异[17],利用该差异统计其标准差(standard deviation, STD),并与法国空间研究中心(Centre National D’études Spatiales,CNES)实时播发的卫星钟差产品进行对比,对多模GNSS实时卫星钟差估计算法的有效性进行评估。考虑到卫星钟差实时估计的收敛时间,仅对2021年5月20日—2021年5月22日(年积日第140—142天)的实时卫星钟差进行精度分析,并在图3中展示了单颗卫星连续3 d的钟误差时间序列,BDS卫星钟差的整体稳定性要弱于其他导航卫星系统。同时,GNSS实时估计钟差与CNES实时钟差产品的精度对比如表2和图4所示。其中,表2为GPS、GLONASS、Galileo、BDS各系统单天的平均STD对比,图4为GNSS每颗卫星连续3 d的平均STD对比。

图3 GPS/GLONASS/Galileo/BDS卫星钟误差的时间序列Fig.3 Time series of GPS/GLONASS/Galileo/BDS satellite clock error

图4 GPS/GLONASS/Galileo/BDS卫星钟差精度对比Fig.4 Accuracy comparison of GPS/GLONASS/Galileo/BDS satellite clock

表2 卫星钟差精度统计Tab.2 Accuracy statistics of satellites clock offset ns

由于卫星钟差与轨道误差的耦合性,在实时估计多模GNSS卫星钟差时,90%的轨道径向误差会被钟差所吸收,从而导致实时估计钟差和最终精密钟差两者之间并不完全吻合[18]。结合表2和图4可以得出,实时估计GPS卫星钟差的STD数值区间为0.044～0.094 ns,每颗GPS卫星的钟差精度均优于0.1 ns,整体稳定性较好,所有卫星的STD均值为0.064 ns。除G18卫星的钟差STD值为0.215 ns外,CNES实时播发的GPS卫星钟差产品与自主估计的GPS实时钟差精度基本一致,这可能与G18卫星的型号为Block Ⅲ,发射时间较晚,CNES的轨道模型还未对其完全精化有关。对于Galileo系统而言,实时估计卫星钟差与CNES实时卫星钟差二次差结果序列相差不大,所有Galileo卫星的STD均值分别为0.058 ns和0.080 ns,且每颗卫星的钟差精度基本上均能优于0.1 ns,能够达到与GPS卫星钟差相当的精度水平,这可能与Galileo卫星所搭载的高精度氢原子钟有关。CNES实时播发的GLONASS钟差产品中,R09卫星连续3 d的钟差精度均将达到1.5 ns,这是因为该卫星的钟差序列出现多次中断,导致精度异常。对于其他GLONASS卫星,二次差结果序列的STD数值区间为0.087～0.598 ns,STD均值为0.232 ns,稍差于自主估计的GLONASS实时卫星钟差精度,所有卫星的STD均值为0.163 ns。对于CNES实时BDS卫星钟差产品,地球静止轨道(geostationary Earth orbit,GEO)卫星的钟差精度明显差于中圆地球轨道(medium Earth orbit,MEO)卫星和倾斜地球同步轨道(inclined geosynchronous orbit,IGSO)卫星,导致BDS卫星钟差整体精度较差。除GEO卫星外,其余MEO/IGSO卫星的STD均值为0.294 ns,大约是自主估计BDS MEO/IGSO实时卫星钟差STD均值的2倍。从整体上看,自主估计的GPS、GLONASS、Galileo、BDS GEO、BDS IGSO和BDS MEO实时卫星钟差STD均值分别为0.064、0.163、0.058、0.232、0.154和0.179 ns,要稍优于CNES实时播发的卫星钟差产品,两者G/R/E/C四系统的STD均值分别为0.142 ns和0.347 ns,这是由于在实时估计GNSS卫星钟差时,二者选取测站的数据质量和解算策略存在差异所致。

3 震时地表形变分析

3.1 2021年云南漾濞Mw 6.4地震

为进一步评估实时估计多模GNSS卫星钟差的精度和稳定性,利用中国大陆构造环境监测网络(简称“陆态网络”)基准站2021年5月21日21:45—21:59时段内的高频GNSS观测数据,基准站分布情况如图5所示,采用表1中实时动态PPP的数据处理策略,基于5 s历元间隔的GNSS实时卫星钟差,在实时处理的定位模式下,通过序贯最小二乘估计方法,历元之间不继承坐标信息,分别进行GPS单系统、GNSS多系统的高频动态PPP,并与武汉大学最终精密卫星钟差的GNSS多系统组合动态PPP结果进行对比(表3),从而得到不同实验下2021年云南漾濞Mw6.4地震发震时段内单历元高频动态解的时间序列。上述3种高频动态PPP实验分别对应实验1、实验2和实验3,在不同定位实验下,云南下关站(XIAG)、云南云龙站(YNYL)和云南丽江站(YNLJ)在东(East,E)、北(North,N)和高程(Up,U)方向上的动态位移时间序列如图6所示,图中垂直于横轴的红色实线为漾濞地震的发震时刻。

图5 漾濞地震震中周边GNSS基准站分布Fig.5 Distribution of GNSS reference station around the epicenter of Yangbi earthquake

图6 漾濞地震发生时段的动态位移时间序列Fig.6 The time series of dynamic displacement during the Yangbi earthquake

表3 不同测站高频动态PPP的标准差Tab.3 Standard deviation of high-rate kinematic PPP at different stations

由图6和表3可见,随着其他系统的加入,基于实时估计卫星钟差的GNSS多系统融合精密单点定位的动态位移波形更为稳定,其高频动态解在E、N、U方向上的平均标准差分别为0.4、0.3、1.0 cm,相对于GPS单系统动态定位结果的平均标准差,在三个方向上分别提升了33%、25%和50%,印证了多系统融合定位的优越性[19]。同时,自主估计的GNSS实时卫星钟差与武汉大学最终精密卫星钟差的定位性能相当,其GNSS多系统组合PPP动态解的平均标准差能够达到水平方向0.5 cm左右,垂直方向1.0 cm左右的精度水平。对于2021年5月21日云南漾濞地震,GNSS高频动态解可以观测到距震中39.7 km的XIAG站在水平方向上具有相对明显的波动,这与文献[20]漾濞地震的水平位移主要发生在震中距50 km范围内相吻合。可见,基于上述数学模型估计的多模GNSS卫星钟差能够用于实时监测地震发生阶段地表的动态形变特征。

3.2 2021年青海玛多Mw 7.4地震

据中国地震台网正式测定:北京时间2021年5月22日2时4分,青海省果洛藏族自治州玛多县(震中34.59°N,98.34°E)发生7.4级地震,震源深度17 km[21]。此次地震是中国大陆地区内继2008年汶川Mw8.0地震后发生的震级最大的地震,其剧烈的地壳形变信号能够被区域内连续运行的GNSS基准站可靠监测[22]。为具体分析玛多Mw7.4地震发生阶段对周围地表形变的影响,利用中国地震台网中心提供的震源周边地区10个陆态网络基准站的高频GNSS观测数据,这10个基准站分别为青海玛多站(QHMD)、青海玛沁站(QHMQ)、青海都兰站(QHDL)、青海班玛站(QHBM)、青海德令哈站(DLHA)、甘肃玛曲站(GSMA)、四川甘孜站(SCGZ)、青海格尔木站(QHGE)、青海西宁站(XNIN)和西藏昌都站(XZCD),基准站分布情况如图7所示。以地震发生前3 d的精密坐标的平均值为参考值,基于自主估计的2021年年积日第141天实时卫星钟差进行GNSS多系统组合实时动态PPP,获取玛多地震震后5 min时段内不同基准站高频动态解在E、N、U方向上随时间变化的位移序列,如图8所示。同时,对上述10个基准站震前3 d以及震后1 d的低频GNSS观测数据进行事后静态PPP处理得到精密坐标的单天解,然后分别计算3个坐标分量上地震前后基准站位置的坐标平均值,通过差分获得玛多地震高可靠性的同震形变[23],并与GNSS多系统组合实时动态PPP获得的同震形变进行对比,具体形变结果如表4所示。

图7 玛多地震震中周边GNSS基准站分布Fig.7 Distribution of GNSS reference stations around the epicenter of Madoi earthquake

图8 不同测站动态位移时间序列Fig.8 The time series of dynamic displacement at different stations

表4 基于高频和低频GNSS观测的地表形变Tab.4 The surface deformation from 1 s and 30 s sampling GNSS observation data

由图8可见,当玛多地震发生后,10个陆态网络基准站在不同时刻均发生不同程度的波动,且开始波动时刻随基准站震中距的增大而延迟,影响区域可达距震中约400 km。其中,距震中约39 km的陆态网络GNSS基准站QHMD最先感知到地震波,受玛多地震的影响也最大,其位置在东西向和南北向均发生不可逆的永久性位移,东西向为23.3 cm,南北向为8.5 cm。不同于东西方向上的动态位移特征,QHMQ和GSMA相对于其他基准站在南北向受到玛多地震的影响最为显著,其震时最大动态位移分别约为23.5 cm和15.6 cm,说明地震对震源周围站点的影响不仅和震中距有关,与其方位角也有一定的相关性[24]。然而,玛多地震在垂直方向上对震源周边基准站的影响相对较弱,10个基准站均未出现显著的位移形变。从表4可以得出,受益于GNSS多系统融合具有更为稳健的定位精度,其实时PPP高频动态解能够有效地获取地震发生阶段各测站在水平方向上较为准确的地表形变量,与低频GNSS观测的形变结果基本一致,而两者在垂直方向上的差值最大为1.1 cm,可靠性相对较差。综上所述,基于多模GNSS实时卫星钟差的PPP高频动态解可以有效地揭示玛多地震的震时地表形变特征。