实时高频（≥1 Hz）的全球导航卫星系统（global navigation satellite system，GNSS）可以克服地震仪饱和效应和强震仪基线偏差^[1-2]，在同震形变监测领域得到广泛应用^[3-7]，因此GNSS地震学这一交叉学科应运而生^[8]。由于全球定位系统（global positioning system，GPS）建成时间最早、使用最广泛、表现最成熟，GNSS地震学在一段时间内多特指GPS。随着GLONASS、北斗三号、Galileo的建设发展^[9-10]，GNSS多系统组合在同震形变监测中精度更高、稳定性更好、震源参数反演更可靠^[11-14]。2020年6月，北斗三号完成全球组网，文献[15]使用振动平台进行了一系列模拟实验，证实了北斗三号和Galileo对地震位移拾取的积极作用，但北斗三号在天然地震中的应用报道尚不多见。

随着硬件和储存技术的发展，当前GNSS接收机采样率已越来越接近地震仪和强震仪，可达到50~100 Hz^[16-18]（本文称之为超高频），中国大陆构造环境监测网络（简称陆态网络）的GNSS测站也已普遍升级，可提供50 Hz的观测记录^[19]。但是目前仍有许多国家的GNSS同震形变观测网以1 Hz、5 Hz为主流，如美国的板块边界观测网络（plate boundary observatory，PBO），仅部分站可达10 Hz^[20]。虽然提高采样率有利于削弱位移时间序列的混叠效应^[20-21]，也存在获得超高采样率的速度、加速度序列的可能性^[22-23]，使GNSS接收机实现类似强震仪、地震仪的作用，但同时也极大地增加了数据储存、管理、计算的成本。

北京时间2021-05-22T02：04：11青海玛多发生了Mw 7.4地震（https://news.ceic.ac.cn/），其同震过程被附近的陆态网络和当地的连续运行参考站（continuously operating reference stations，CORS）记录^[24-27]。本文以此次地震为例，研究多系统组合同震位移拾取的精度，分析北斗三号的地表位移拾取性能，并且针对GNSS地震监测采样率设置存在的一些争议^{[8，20，28]}，测试超高频数据获取同震位移、速度、加速度的性能，分析是否有必要将采样率提升至50 Hz。

1 数据和方法

1.1   GNSS数据与处理方法

本文选取了震中附近7个陆态网络参考站（平均距离约280 km），获取2021年玛多Mw 7.4地震的1 Hz同震位移进行分析，如图 1所示。其中G表示GPS，R表示GLONASS，E表示Galileo，C2表示北斗二号，C3表示北斗三号，C2+3表示北斗二号和北斗三号组合，其他字母组合含义类似。图 1（b）中的位移波从上到下按照测站到震中的距离排序，均可见明显同震信号。现阶段陆态网络中绝大部分参考站支持GRC2三系统，部分参考站支持GREC2四系统。2021年玛多地震发生时，陆态网络附近参考站尚不支持北斗三号，为了开展相关研究，本文选取了中国青海省CORS站中的QSHZ站观测数据，该站可观测包括北斗三号在内的四大卫星导航系统，其观测数据详细的信息见表 1。本文所涉及的GNSS数据均采用PANDA软件^[29-31]进行计算，其精密单点定位（precise point positioning，PPP）功能可以处理高达50 Hz的数据^[32]，处理策略和参数设置见表 2。

图  1  2021年玛多Mw 7.4地震的测站分布及同震位移

Figure  1.  Distribution of GNSS Stations and Coseismic Displacement of the 2021 Maduo Mw 7.4 Earthquake

下载: 全尺寸图片 幻灯片

表  1  GNSS数据信息

Table  1.  Information of GNSS Data

数据来源	测站个数	GNSS采样率/Hz	接收机型号	卫星系统
陆态网络参考站	7	1，50	Trimble NetR9	GRC2，GREC2
青海省QSHZ站	1	1	中海达VNet8	GREC2+3

下载: 导出CSV 

| 显示表格

表  2  PPP处理策略和参数设置

Table  2.  Processing Strategy and Parameters Setting of PPP

策略与参数	具体描述
组合方式	无电离层组合
截止高度角/（°）	7
观测值定权	高度角加权，伪距、载波先验误差：0.6 m、0.01周
系统定权	G∶R∶E∶C2∶C3=1.5∶1∶1.5∶1∶1.5
估计方式	最小二乘
天线改正	IGS14.atx
周跳探测	Turboedit
地球自转参数	IGS EOP
对流层估计	Saastamoinen模型+分段线性估计
精密星历/钟差	CODE：30 s钟差，15 min轨道（GREC2+3）

下载: 导出CSV 

| 显示表格

在处理陆态网络和青海省QSHZ站1 Hz数据时，考虑到卫星钟差在30 s范围内高频抖动并不剧烈，对于1 s采样率的GNSS数据，使用5 s间隔钟差与30 s间隔钟差获得的精度基本一致^[33]，故选择欧洲定轨中心（Center for Orbit Determination in Europe，CODE）提供的30 s钟差和15 min精密星历产品。

通过PPP获取位移时间序列后，对其进行一阶离散差分便可获得速度序列，再对速度差分可获得加速度序列，计算式为：

式中，s、v、a分别表示位移、速度、加速度；i表示历元；t表示时间间隔。

1.2   强震仪数据与处理方法

为了研究超高频GNSS差分获取速度和加速度的能力，本文选择QHMQ站及其附近2 km内的一个强震仪测站DAW63进行比较，DAW63与QHMQ具有高度相关性，可提供采样率为200 Hz的加速度序列。通过对原始加速度序列积分可获得速度序列，但由于基线偏差的存在，此时的速度序列包含线性漂移，以4阶巴特沃斯高通滤波的方式去除速度序列中频率在0.02 Hz以下的低频信号，再次积分即可得到地表形变。以上数据处理过程与经验基线自动校正方法^[34]有一定区别，会对静态同震位移造成一定的削弱，但考虑到本文重点关注高频同震位移波且QHMQ站静态同震位移较小，上述处理方法是合理的。

2 结果和分析

2.1   北斗二号/GNSS高频同震位移拾取

为了分析陆态网络测站1 Hz数据同震位移拾取精度，计算2021-05-21T13：00—18：00共计5 h震前时间序列与功率谱密度，出于篇幅考虑，图 2仅展示SCGZ单站的时间序列和各测站的平均功率谱密度。

图  2  SCGZ测站时间序列和各测站平均功率谱密度

Figure  2.  Time Series at SCGZ Station and Average Power Spectrum Density of All Stations

下载: 全尺寸图片 幻灯片

由图 2（a）可以直观看出，SCGZ单站的多系统组合（GRC2）相较于单系统G能够在一定程度上削弱低频噪声和跳变，使地表形变拾取更加稳健可靠。由图 2（b）可以看出，各测站的U方向功率谱密度较大，在任意频段均要高于E、N方向。随着系统数目的增加，各频段噪声普遍明显减弱，GREC2相比GRC2效果更好。单系统G在E方向噪声与GRC2差异不大，甚至部分频段噪声更低，但GREC2仍然表现最好。高频段出现个别峰值的原因在于位移时间序列存在一定的跳变，而这些跳变可能是由精密钟差的高频钟跳所导致。

计算不同系统时间序列的标准差并作站间平均，结果见表 3。单系统G在E、N方向的标准差均值约为1 cm，U方向约为4 cm，主要是因为受卫星几何构型和多路径效应等影响，U方向误差较大。GRC2三系统组合相比于单系统G，在E、N、U方向上的精度分别提升了14%、42%、56%，GREC2四系统组合相比于单系统G，在E、N、U方向上的精度分别提升了17%、47%、57%，印证了多系统组合的优越性。

表  3  不同系统时间序列的标准差均值/mm

Table  3.  Average Standard Deviation of Time Series of Diffferent Systems /mm

系统	E方向	N方向	U方向
G	9.4	12.4	39.1
GRC2	8.1	7.1	17.0
GREC2	7.8	6.5	16.8

下载: 导出CSV 

| 显示表格

2.2   北斗三号/GNSS高频同震位移拾取

北斗二号虽然在亚太地区卫星数与GPS相当，但其包含若干颗地球同步轨道（geosynchronous earth orbit，GEO）卫星，几何构型较差，对定位精度提升有限，单独采用北斗二号并不能取得良好的同震位移拾取效果。随着北斗三号的全面建成，这一情况得到了极大改善。

对青海省QSHZ站进行PPP解算，截取与§2.1相同的时段进行时间序列分析，并计算了不同系统组合时间序列的标准差，结果见表 4。QSHZ站震前时间序列的功率谱密度如图 3所示。

表  4  QSHZ测站不同系统时间序列标准差/mm

Table  4.  Standard Deviation of Time Series of Diffferent Systems at QSHZ Station /mm

系统	E方向	N方向	U方向
G	9.2	9.5	43.6
E	7.5	10.4	19.6
C2+3	5.3	5.4	18.1
GREC2	4.8	5.4	16.2
GREC2+3	3.4	4.7	12.4

下载: 导出CSV 

| 显示表格

图  3  QSHZ测站的时间序列和功率谱密度

Figure  3.  Time Series and Power Spectrum Density at QSHZ Station

下载: 全尺寸图片 幻灯片

由图 3（a）可以看出，附加北斗三号的四系统组合精度最高，相较于单GPS定位结果，其在E、N、U方向上的精度均有显著提高。单GPS结果受到几何构型、多路径效应等影响，难以单独提供稳定的时间序列，尤其是在N方向末段和U方向上结果较差，此部分误差和SCGZ测站较为相似，说明误差的起因更多在于卫星几何构型较差。就单系统而言，随着北斗三号的建成，北斗系统不仅可以单独实现同震位移监测，而且在单系统中效果最佳，其原因在于地震发生在亚太地区，北斗可观测卫星数更多，几何结构更好，且北斗三号搭载的原子钟更稳定，随着北斗三号系统的加入，同震形变的噪声能进一步得到削弱，提高了拾取精度和可靠性；从此时段结果来看，Galileo高频噪声较大，但低频噪声要低于GPS。由此可知，单系统受制于卫星几何构型等因素，其时间序列不如多系统稳健，采用多系统组合的方式对同震位移拾取有重要意义。

由图 3（b）可以看出，QSHZ测站震前时间序列的功率谱密度的基本特征与陆态网络平均功率谱密度一致。单系统之间相比，北斗在各频段的噪声均为最低；Galileo次之，且在0.1 Hz频段及以上噪声略大；GPS低频噪声水平最大，不够稳定，对应了图 3（a）中N方向末段和U方向的结果，四系统组合的噪声水平相比于单系统GPS在任意频段均有大幅削弱。以上结果表明，北斗不仅能单独实现高精度定位，且相比于GPS和Galileo噪声更小；在加入北斗三号后，GREC2+3的噪声相对于GREC2（代表现阶段陆态网络的普遍配置）在各个频段均有明显降低，时间序列稳定性进一步上升，说明北斗三号可以对同震位移拾取做出重要贡献。

2.3   陆态网络超高频位移/速度/加速度拾取

以QHMQ和DAW63为例，设置一系列采样率（1 Hz、2 Hz、5 Hz、10 Hz、25 Hz、50 Hz），获得相应的同震位移、速度和加速度序列。考虑到2021年青海玛多Mw 7.4地震为走滑型地震，以E、N方向的位移为主导，且GNSS在U方向的精度较差，本文重点研究E方向和N方向，GNSS与强震仪结果比较如图 4所示。为保证超高频处理结果的可靠性，将钟差间隔提升至5 s，选用武汉大学提供的5 s钟差以及15 min轨道计算，GRC2三系统参与运算，其余策略和表 2一致。

图  4  QHMQ站与DAW63站的同震位移、速度和加速度

Figure  4.  Coseismic Deformation, Velocity and Acceleration at QHMQ and DAW63 Stations

下载: 全尺寸图片 幻灯片

对图 4进行分析可得如下结论：

1）就位移拾取而言，GNSS和强震仪在E、N方向上的位移拾取保持高度一致性；混叠效应随采样率提升而削弱，在5~10 Hz处达到饱和；N方向相比于E方向振幅更大，混叠效应整体较弱，50 Hz和1 Hz相比结果无明显改善。

2）就速度拾取而言，差分噪声随着采样率的提升而增强，但同时混叠效应减弱；GNSS在N方向的2~5 Hz与强震仪保持了较好的一致性，但在E方向效果一般，原因在于N方向比E方向振幅更大，故速度序列噪声整体要更小；速度差分效果随振幅增大而改善，这表明GNSS更适合于近场以及大型地震监测。

3）就加速度拾取而言，通过二次差分获得的结果始终未取得良好效果；在10 Hz时达到了强震仪的量级，但符合程度较差；随着采样率提升，噪声水平持续增强，有效信号被掩盖。

对于其他震级的地震是否能得出类似结论有待进一步研究，值得注意的是，受限于定位精度，GNSS对中小型地震同震位移的拾取能力较弱，需要缩短测站与震中距离以实现更好的观测效果。

为了进一步研究超高频数据对时间序列存在的混叠效应的削弱，选择3个与震中距离不同且具有代表性的测站QHMD（38 km）、QHMQ（175 km）和QHGE（366 km）进行功率谱密度分析，分析时段为2021-05-21T17：00—20：00（包含同震形变），结果如图 5所示。由图 5可以看出，随着采样率提升，高频段的混叠效应逐渐被削弱，在采样率达到10 Hz时饱和；继续提高采样率并不能进一步削弱混叠效应，这一特点在E、N方向上一致，在距震源远近不同的测站也表现一致；同震形变反映在0.01~0.1 Hz处的峰值，N方向较E方向形变更大，故功率谱密度的峰值更高。

图  5  各采样率的时间序列功率谱密度

Figure  5.  Power Spectrum Density of Time Series with Different Sampling Rates

下载: 全尺寸图片 幻灯片

综上所述，50 Hz超高频数据未能显著削弱时间序列存在的混叠效应，其效果与5 Hz、10 Hz接近；在振幅达到分米级时，通过差分的方式可以获得较高精度的速度序列；二次差分获取加速度较为困难，难以取得良好的效果。

3 结语

本文以2021年青海玛多Mw 7.4地震为例，研究了GNSS多系统组合，尤其是北斗三号对地表位移拾取精度的提升，还分析了50 Hz超高频数据拾取位移、速度、加速度的能力，得到如下结论：

1）多系统组合对于削弱低频多路径噪声、提高位移拾取稳定性有显著作用。融入了北斗三号的GREC2+3四系统组合效果最优，在E、N、U方向上的精度分别达到了3.4 mm、4.7 mm、12.4 mm，相较于GREC2的组合方式分别提高了29%、13%、23%。

2）对2021年玛多地震设置1~50 Hz内多个采样率，位移序列混叠效应的削弱在5~10 Hz时达到饱和；采样率在5 Hz时差分速度序列效果最佳；在1~50 Hz，通过二次差分获取加速度的效果均不理想；超高频GNSS在同震位移监测中的作用需要进一步挖掘。

3）将采样率提升到50 Hz，增大了数据运算、储存、管理的成本（数据量约20 GB/站/天），既未能显著削弱混叠效应，也难以提供有效的超高频速度、加速度序列。5~10 Hz已经能满足大部分场景的需求，加密基准站或升级原有接收机使其可观测北斗三号卫星，将会取得更好的同震形变监测效果。