Seismic phase recognition of coseismic data from high-rate GNSS based on S transform
-
摘要: 为研究高频全球导航卫星系统（GNSS）信号的P波和S波到时，本文利用2008年汶川MS8.0地震震时部分站点记录到的1 Hz高频GNSS数据，采用广义S变换将同震信号进行二维时频域平面分解，以此对P波和S波到时进行识别。本文通过调整广义S变换中的调节因子λa和p得出，当λa＝1.05，p＝1.05时，变换所得图像中的P波和S波到时均较为明显。结果表明，S变换在高频GNSS震相识别中效果明显，可作为GNSS地震学应用中一项有效的技术手段。Abstract: The high-rate GNSS data recorded by some stations during the 2008 Wenchuan MS8.0 earthquake was used in this paper to identify the arrive time of P-wave and S-wave in 1 Hzhigh-rate GNSS. Coseismic signals are decomposed in two-dimensional time-frequency domain using generalized S transform. The adjustment factors λa＝1.05, p＝1.05 is used to identify the arrival times of P-wave and S-wave, which could get a better result. The S transform is effec-tive and feasible in the high-rate GNSS seismic phase recognition.
-
Keywords:
- high-rate GNSS /
- S transform /
- arrival time /
- seismic phase identification
-
引言
对于地震学而言，高频全球导航卫星系统（global navigation satellite system，简写为GNSS）是能够记录到地震波传播信息的一种新手段，其在强震中可直接提供时变位移场观测信息，弥补了传统观测手段在强震观测时出现诸如震级饱和、漂移等现象的不足。但由于高频GNSS数据存在信噪比低、信号稳定性差等问题，在预警中需要与测/强震观测手段联合使用。地震波信号往往是一种非线性、非平稳信号，故研究其频率的局部特性显得至关重要，这就需要用时频分析方法将一维时域信号转换到二维的时频平面。而傅里叶变换和反变换虽然能从时间域和频率域对平稳信号进行分析，却不能同时保留时间和频率信息。故此，本文拟以2008年汶川地震时获取的高频GNSS同震信号为例，通过对其进行广义S变换，系统分析变换后同震信号的时间域和频率域，以期能够更加准确地判断P波和S波到时。同时，广义S变换的使用也为高频GNSS在地震学中的应用提供了一条从多层面提取地震波形到时的新途径。
2. S变换
美国地球物理学家Stockwell等（1996）综合了小波变换和短时傅里叶变换的特点，提出了一种分析非平稳信号的新方法—S变换。其时频随着频率的不同而变化，具有多分辨率、逆变换唯一等特点，同时获得的二维时频谱也是基于傅里叶变换的基本理论。该方法已经在地震勘探、电能信号分析、故障检测等多个领域得到了广泛的应用（齐春燕等，2010；周竹生，陈友良，2011；郑成龙，王宝善，2015）。
S变换的表示形式为
$$ {{S}}(\tau {\text{，}}f) {\text{＝}} \int _{ {\text{－}} \infty }^{ {\text{＋}} \infty }h(t)\left\{ {\frac{{|f|}}{{\sqrt {2\pi } }}\exp \left[\frac{{ {\text{－}} {f^2}{{(t {\text{－}} \tau )}^2}}}{2}\right]\exp ( {\text{－}} {\rm i}2\pi ft)} \right\}{\rm d}t{\text{，}} $$ (1) 式中：τ为窗函数的中心点，控制窗函数在时间轴上的位置；f为采样频率；h（t）为原始信号。S变换使用的高斯窗函数gf （t）和基本小波ωf （t）分别定义为
$$ {{g}_f}(t) {\text{＝}} \frac{{|f|}}{{\sqrt {2\pi } }}\exp \left({\text{－}} \frac{{{t^2}{f^2}}}{2}\right){\text{，}} $$ (2) $$ {\omega _f}(t) {\text{＝}} \frac{{|f|}}{{\sqrt {2\pi } }}\exp \left( {\text{－}} \frac{{{t^2}{f^2}}}{2} {\text{－}} {\rm i}2\pi ft\right) {\text{＝}} {g_f}(t)\exp ( {\text{－}} {\rm i}2\pi ft){\text{，}} $$ (3) 从上式可以看到，S变换在分析高频信号或低频信号时具有良好的特性：对于高频信号，设计窗口比较窄，可满足高时间分辨率的要求；对于低频信号，设计窗口比较宽，可满足高频率分辨率的要求。
相比于单分辨率的短时傅里叶变换，S变换属于多分辨率时频分析方法，可以更加有效地对复杂的地震波信号进行分析。同时由于信号x（t）可以由S变换S（τ，f ）通过逆变换进行重构，S（τ，f ）的S逆变换为
$$ {\rm{x}}(t) {\text{＝}} \int_{ {\text{－}} \infty }^{ {\text{＋}} \infty } {\left[ {\int_{ {\text{－}} \infty }^{ {\text{＋}} \infty } {S(\tau{\text{，}} f\;)\rm d\tau } } \right]} \exp ({\rm i}2\pi ft){\rm {d}}f{\text{，}} $$ (4) S变换所具有的无损逆变换的特性，使得信号在经过时频域处理后可逆变换为时间信号。
由于S变换的基本高斯窗函数固定，因此在处理复杂的地震信号时，不能根据需求对频率域和时间域的分辨率进行调节，缺乏必要的灵活性。随着后期的不断发展，多位研究人员对S变换中的高斯窗函数和基本小波进行改变，构建出了广义的S变换（Mansinha等，1997；McFadden et al，1999；高静怀等，2003；Pinnegar，Mansinha ，2003；陈学华等，2006；齐春燕等，2010；周竹生，陈友良，2011）。
本文使用的广义S变换是通过两个调节因子对高斯窗进行改造（陈学华等，2006），以便根据实际要求对高斯窗作出调整，从而更好地分析复杂频率的地震波形，其表达式为
$$ {{S}}(\tau {\text{，}} f) {\text{＝}} \int\nolimits_{ {\text{－}} \infty }^\infty {x(t)\frac{{|{\lambda _{\rm a}}||f{|^p}}}{{2\pi }}} {{\rm exp}({\frac{{ {\text{－}} {\lambda ^2_{\rm{a}}}{f^{2p}}{{(t {\text{－}} \tau )}^2}}}{2}}}){{\rm exp}({ {\text{－}} {\rm i}2\pi ft}}){\rm d}t{\text{，}} $$ (5) 式中λa和p均为调节因子。选定p值后，当λa＞1时，时窗宽度随频率呈反比变化的速度加快；当λa＜1时，时窗宽度随频率呈反比变化的速度减慢；选定λa之后，通过调节p值可以对频率域的分辨率进行调整。引入λa和p两个参数之后，可以根据实际信号的频率分布特点和时频分析的侧重点，改变高斯窗函数随频率变化的规律，对时频分辨率进行调整，但由于其运算过程与标准的S变换相似，因此并不需要额外增加运算量。当两个调节因子λa和p均设置为1时，则为标准S变换。
3. 到时分析中的S变换
本文收集了2008年汶川MS8.0地震中重庆台网和昆明台网共28个站点的高频GNSS信号（采样率为1 Hz），采用GAMIT/GLOBK软件包（Bock，2000）中的高频GNSS处理模块TRACK （Chen，1998）对附近站点的高频GNSS数据进行解算。由于所选取数据按小时文件进行解算，时间最长为1小时，故未考虑昼夜变化等周期性较长的噪音。按照GNSS站点的位置及分布情况，选取15个站点的数据使用广义S变换进行分析，并挑选其中6个站点的同震波形进行展示，如图1所示。
根据距离震中位置的远近选取了bana （距震中约346 km）、djxm （距震中381 km）、fdlh （距震中477 km）和ghxx （距震中661 km）等4个站点的同震波形，对其进行广义S变换后识别出P波、S波到时并进行展示（图1），在计算P波和S波的理论到时中使用了频率-波数法（Zhu，2002），计算到时的速度模型为一维速度模型（吴建平等，2009；殷海涛等，2013；邓文哲等，2014；王小娜等，2015；杨彦明等，2016）。
bana站点距离震中346 km，该站点记录高频GNSS信号的频率为1 Hz，根据奈奎斯特（Nyquist）采样定理，其可记录到的地震波信号的最大频率为0.5 Hz。bana站点东西向同震信号广义S变换的结果表明（图2），地震信号的能量主要集中于0—0.3 Hz，所研究各站点的P波和S波理论到时均在300 s以内。为了更加直观地显示所研究的到时区域，在后续分析中，按照S变换后信号频率为0—0.3 Hz，取地震发生后的300 s进行分析。结果显示，bana站点的P波理论到时为56.7 s，S波理论到时为96.1 s。如图1所示，经过广义S变换后的图像显示，在P波和S波理论到时（图2）的位置，均有着较明显的能量汇集。在本文分析的地震中，P波和S波周期均较长、频率均相近，故P波的识别周期为20—50 s左右的初至拾取、高于背景噪声的能量汇集，S波能量为P波能量的2.5倍以上（水平向）。
![]()
图 2 bana站点东西向同震信号的S变换结果Figure 2. S transform of coseismic signal in the EW directionin at the bana station在使用广义S变换进行频率分析时，通过对两个调节因子λa和p进行调整，以获得最优的分辨率。以bana站东西向为例进行分析，结果如图3所示。
![]()
图 3 通过微调调节因子λa和p所产生的同震信号广义S变换结果展示Figure 3. S transform results of the same coseismic signal with different adjustment factor λa and p（a） λa＝1，p＝1.1；（b） λa＝1，p＝1.2；（c） λa＝1.05，p＝1.05；（d） λa＝1.05，p＝1.05在λa不变的情况下，提高因子p值，可以增加时间域的分辨率，详见图2和图3a。但是p值过高，会导致频率域的分辨率增加，时间域的分辨率降低，从而降低了对初至波的识别能力（图3b）。当调节因子λa＝1.05，p＝1.05时，对于P波和S波的识别均较为清晰。图3d将图3c进行三维展示，更加直观地展示了震后随着时间的变化，不同周期的波形所汇聚的能量大小。
对于震中距为381 km的djxm站点，其P波理论到时为62.5 s，S波理论到时为105.9 s。将P波和S波理论到时分别标注于S变换结果上（选用的调节因子为λa＝1.05，p＝1.05），可以看出，在同震信号中，djxm站记录的到时要比理论到时延迟约4—5 s （图4a），可能与本文使用的速度结构有一定关系，而djxm站在垂向上很难识别S波到时（图4b）。在原始波形记录中，无论垂直向还是水平向均已很难判定P波和S波到时。这表明通过微调广义S变换的两个调节因子λa和p，能够得到较好的到时分辨率。鉴于目前GNSS预警中的触发方式为阈值触发（Allen，Ziv，2011；Melgar et al，2012），很难通过反演的方式进行精确的震源定位。而相比于人工或者阈值判断，通过对地震信号进行广义S变换，能够对到时作出更为精准的判断。
![]()
图 4 各高频GNSS站点同震信号S变换结果展示（a） djxm站点EW向同震信号S变换；（b） djxm站点UD向同震信号S变换；（c） fdlh站点EW向同震信号S变换；（d） fdlh站点UD向同震信号S变换；（e） ghxx站点EW向同震信号S变换；（f） ghxx站点UD向同震信号S变换Figure 4. S treansform of coseismic signals for several stations（a） S transform of EW coseismic signal in djxm；（b） S transform of UD coseismic signal in djxm；（c） S transform of the EW coseismic signal in fdlh ；（d） S transform of the UD coseismic signal in fdlh；（e） S transform of the EW coseismic signal in ghxx；（f） S transform of the UD coseismic signal in ghxx对于震中距为477 km的fdlh站，P波和S波理论到时分别为79 s和133 s。东西向同震信号的广义S变换（选用的调节因子为λa＝1.05，p＝1.05）波形识别基本与理论到时一致，相差2—3 s （图4c）。fdlh站垂直向上的分辨率也可以满足波形识别的需求（图4d）。其结果较好，可能与fdlh站的背景噪声较低有关。三维图像中，面波的识别更加清晰（图4c）。
对于震中距为661 km的ghxx站，在同震信号中很难看到P波到时和S波到时。经过计算，其P波理论到时为108.3 s，S波理论到时为183 s。在图4e中可以看到，通过广义S变换后对S波的识别较好，对P波的识别稍显模糊。而对于垂直向（图4f），通过广义S变换之后，P波已经无法识别，S波的识别时间较理论地震波形要晚将近20 s。这可能与GNSS背景噪声较高，且垂直向的噪音水平高于水平向2—3倍有关。而同震信号的分析中，P波和S波的到时基本均无法分辨。这说明通过广义S变换，即使是震中距大于600 km，也能有效地提取到地震波到时。
4. 讨论与结论
对于实时解算的高频GNSS数据，其背景噪声相对较大，使用广义S变换，可以更加准确地得到P波到时。即使远距离的站点记录，在时频域中也可以很好地指示出P波初至到时。这为高频GNSS数据在地震学中的应用和分析提供了一种新的途径。其优势在于：
1） 广义S变换集合了短时傅里叶变换和小波变换的优点，在地球物理领域及电能分析领域取得了很多成果，但依旧存在很多应用前景。本文通过对高频GNSS同震信号进行S变换，从以往的只针对单独的时间域或频率域分析高频GNSS信号，转变为从时间域和频率域同时进行分析。从频率域的角度看，P波和S波的能量均汇聚于0.025 Hz左右，面波的能量汇聚于0.04—0.25 Hz并能够明显地看出面波是分为多个频率区间；从时间域的角度看，P波和S波的理论到时均有明显的能量汇集，而面波到时则更加明显。通过对不同频率的信号进行分析，认为近震范围内P波和S波实际到时与理论到时的分辨率相差仅为1—2 s，中远震其分辨率相差3—5 s，该方法可以作为高频GNSS信号的分析及各类波形到时的提取手段。
2） 通过对高频GNSS信号使用S变换可见，有效信号与噪音信号的分离较为明显，说明S变换在降噪去噪方面具有独特的优势。鉴于信号可以进行无损S变换和S逆变换，因此可以使用该方法提高高频GNSS信号的信噪比，尤其是在压制噪音方面具有很大的优势，下一步研究中将根据各波形的到时按时间提取波形并用于专门性的分析。
感谢两位审稿专家对本文提出的意见和建议。
-
![]()
图 2 bana站点东西向同震信号的S变换结果
Figure 2. S transform of coseismic signal in the EW directionin at the bana station
![]()
图 3 通过微调调节因子λa和p所产生的同震信号广义S变换结果展示
Figure 3. S transform results of the same coseismic signal with different adjustment factor λa and p
（a） λa＝1，p＝1.1；（b） λa＝1，p＝1.2；（c） λa＝1.05，p＝1.05；（d） λa＝1.05，p＝1.05
![]()
图 4 各高频GNSS站点同震信号S变换结果展示
（a） djxm站点EW向同震信号S变换；（b） djxm站点UD向同震信号S变换；（c） fdlh站点EW向同震信号S变换；（d） fdlh站点UD向同震信号S变换；（e） ghxx站点EW向同震信号S变换；（f） ghxx站点UD向同震信号S变换
Figure 4. S treansform of coseismic signals for several stations
（a） S transform of EW coseismic signal in djxm；（b） S transform of UD coseismic signal in djxm；（c） S transform of the EW coseismic signal in fdlh ；（d） S transform of the UD coseismic signal in fdlh；（e） S transform of the EW coseismic signal in ghxx；（f） S transform of the UD coseismic signal in ghxx
-
陈学华,贺振华,黄德济. 2006. 基于广义S变换的信号提取与抑噪[J]. 成都理工大学学报:自然科学版,33(4):331–335. doi: 10.3969/j.issn.1671-9727.2006.04.001 Chen X H,He Z H,Huang D J. 2006. Signal extracting and denoising based on generalized S transform[J]. Journal of Chengdu University of Technology :Science &Technology Edition,33(4):331–335 (in Chinese). doi: 10.3969/j.issn.1671-9727.2006.04.001
邓文哲,陈九辉,郭飚,刘启元,李顺成,李昱,尹昕忠,齐少华. 2014. 龙门山断裂带精细速度结构的双差层析成像研究[J]. 地球物理学报,57(4):1101–1110. Deng W Z,Chen J H,Guo B,Liu Q Y,Li S C,Li Y,Yin X Z,Qi S H. 2014. Fine velocity structure of the Longmenshan fault zone by double-difference tomography[J]. Chinese Journal of Geophysics,57(4):1101–1110 (in Chinese).
高静怀,陈文超,李幼铭,田芳. 2003. 广义S变换与薄互层地震响应分析[J]. 地球物理学报,46(4):526–532. doi: 10.3321/j.issn:0001-5733.2003.04.015 Gao J H,Chen W C,Li Y M,Tian F. 2003. Generalized S transform and seismic response analysis of thin interbeds[J]. Chinese Journal of Geophysics,46(4):526–532 (in Chinese). doi: 10.3321/j.issn:0001-5733.2003.04.015
齐春燕,李彦鹏,彭继新,张固澜,张彦斌. 2010. 一种改进的广义S变换[J]. 石油地球物理勘探,45(2):215–218. Qi C Y,Li Y P,Peng J X,Zhang G L,Zhang Y B. 2010. An improved generalized S-transform[J]. Oil Geophysical Prospecting,45(2):215–218 (in Chinese).
王小娜,于湘伟,章文波,王思琳. 2015. 龙门山断裂带南段地壳一维P波速度结构[J]. 地震研究,38(1):16–24. doi: 10.3969/j.issn.1000-0666.2015.01.003 Wang X N,Yu X W,Zhang W B,Wang S L. 2015. 1D P wave velocity structure in the south segment of Longmenshan fault zone[J]. Journal of Seismological Research,38(1):16–24 (in Chinese).
吴建平,黄媛,张天中,明跃红,房立华. 2009. 汶川MS8.0级地震余震分布及周边区域P波三维速度结构研究[J]. 地球物理学报,52(2):320–328. Wu J P,Huang Y,Zhang T Z,Ming Y H,Fang L H. 2009. Aftershock distribution of the MS8.0 Wenchuan earthquake and three dimensional P-wave velocity structure in and around source region[J]. Chinese Journal of Geophysics,52(2):320–328 (in Chinese).
杨彦明,戴勇,张国清,王磊,赵星,王杰民. 2016. 内蒙古中西部地区地震烈度衰减关系[J]. 地震地磁观测与研究,37(1):30–37. Yang Y M,Dai Y,Zhang G Q,Wang L,Zhao X,Wang J M. 2016. Attenuation relation of seismic intensity in the middle and western regions of Inner Mongolia Autonomous Region[J]. Seismological and Geomagnetic Observation and Research,37(1):30–37 (in Chinese).
殷海涛,甘卫军,黄蓓,肖根如,李杰,朱成林. 2013. 日本M9.0级巨震对山东地区地壳活动的影响研究[J]. 地球物理学报,56(5):1497–1505. doi: 10.3969/j.issn.1000-0666.2013.03.012 Yin H T,Gan W J,Huang B,Xiao G R,Li J,Zhu C L. 2013. Study on the effects of Japan M9.0 huge earthquake on the crustal movement of Shandong area[J]. Chinese Journal of Geophysics,56(5):1497–1505 (in Chinese).
郑成龙,王宝善. 2015. S变换在地震资料处理中的应用及展望[J]. 地球物理学进展,30(4):1580–1591. Zheng C L,Wang B S. 2015. Applications of S transform in seismic data processing[J]. Progress in Geophysics,30(4):1580–1591 (in Chinese).
周竹生,陈友良. 2011. 含可变因子的广义S变换及其时频滤波[J]. 煤田地质与勘探,39(6):63–66. doi: 10.3969/j.issn.1001-1986.2011.06.015 Zhou Z S,Chen Y L. 2011. Generalized S-transform with variable-factor and its time-frequency filtering[J]. Coal Geology &Explo-ration,39(6):63–66 (in Chinese). doi: 10.3969/j.issn.1001-1986.2011.06.015
Allen R M,Ziv A. 2011. Application of real-time GPS to earthquake early warning[J]. Geophys Res Lett,38:L16310. doi: 10.1029/2011GL047947
Bock Y,Nikolaidis R M,de Jonge P J,Bevism M. 2000. Instantaneous geodetic positioning at medium distances with the Glo-bal Positioning System[J]. J Geophy Res,105:28223–28253.
Chen G. 1998. GPS Kinematics Positioning for the Airborne Laser Altimetry at Long Valley[D]. Cambridge: Massachusetts Institute of Technology: 173.
Zhu L P. 2002. A note on the dynamic and static displacements from a point source in multilayered media[J]. Geophys J Int,148:619–627. doi: 10.1046/j.1365-246X.2002.01610.x
Mansinha L,Stockwell R G,Lowe R P,Eramian M,Schincariol R A. 1997. Local S-spectrum analysis of 1-D and 2-D data[J]. Phys Earth Planet Inter,103(3/4):329–336. doi: 10.1016/S0031-9201（97）00047-2
McFadden P D,Cook J G,Forster L M. 1990. Decomposition of gear vibration signals by the generalised S transform[J]. Mech Syst Signal Process,13(5):691–707. doi: 10.1006/mssp.1999.1233
Melgar D,Bock Y,Crowell B W. 2012. Real-time centroid moment tensor determination for large earthquakes from local and regional displacement records[J]. Geophys J Int,188(2):703–718. doi: 10.1111/j.1365-246X.2011.05297.x
Pinnegar C R,Mansinha L. 2003. The S-transform with windows of arbitrary and varying shape[J]. Geophysics,68(1):381–385. doi: 10.1190/1.1543223
Stockwell R G,Mansinha L,Lowe R P. 1996. Localization of the complex spectrum:The S transform[J]. IEEE Trans Signal Process,44(4):998–1001. doi: 10.1109/78.492555
-
期刊类型引用(2)
1. 董龙飞，李艳艳，龙海立. 联合CEEMD和改进修正能量比的高频GNSS同震位移P波震相识别方法. 全球定位系统. 2025(01): 25-32+59 . 
百度学术
2. 张环宇，林吉航，陈庭. 利用高频GPS数据反演地震要素. 全球定位系统. 2020(06): 37-45 . 百度学术
其他类型引用(1)
下载:
计量
- 文章访问数: 1119
- HTML全文浏览量: 691
- PDF下载量: 41
- 被引次数: 3
