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
-
引言
高精度相对、绝对重力测量的潮汐改正精度与地球表面潮汐模型的精度密切相关。潮汐模型可利用适当的地球模型理论计算建立,也可利用实测重力时间序列解算建立。实测模型的建立主要取决于力源信号的分离和观测系统的误差改正(Habel,Meurers,2014),而大面积地建立高密度连续重力观测站面临着高成本的挑战(孙和平等,2002),因此,目前重力测量的精密潮汐改正大多采用理论潮汐模型。早在二十世纪初,Love (1909)用一个球对称、无自转的弹性地球体来描述固体潮,给出了全球潮汐因子的平均值为1.16;Molodensky (1961)建立起适合于龙格-库塔(Runge-Kutta)法求解的微分方程组后,将地球勒夫数的计算扩展到29阶,并利用均质、刚性圆球状的地球模型给出了由5个潮波组成的全球理论潮汐模型。此后,Wahr (1981a,b,c)采用无海洋和无气压影响的符合流体静力学原理的均质、刚性旋转椭球状地球模型首次提出潮汐因子可用常数项和纬度依赖项来描述,并根据不同的地球模型给出了3个全球理论潮汐模型,这已成为目前建立全球潮汐模型的主要模式。Dehant等在Wahr模式的基础上建立了流体静力弹性地球(hydrostatic elastic,缩写为He)潮汐模型和流体非静力非弹性地球(non-hydrostatic in-elastic,缩写为NHi)潮汐模型(Dehant,Ducarme,1987;Dehant et al,1999),两模型分别被命名为DDW-He和DDW-NHi模型。此后Mathews (2001)利用甚长基线干涉测量(very long baseline interferometry,缩写为VLBI)修正了初始参考地球模型(preliminary reference Earth model,缩写为PREM)的液核动力学扁率,计算了新的全球理论潮汐模型,该模型简称为M2001。该模型的建立模式与Wahr模式基本一致,只是为了描述潮汐模型的相位,常数项部分使用复数来表示。从Dehant模型开始,潮汐因子的精度从最早的±0.001提高至现今的±0.000 01。
随着全球潮汐模型精度的提高,潮汐观测精度也在提高,且两个模型参数值的相对误差逐步减小,已可利用实测资料对理论潮汐模型进行检验和评估。全球重力潮汐台网主要潮波的潮汐因子(Melchior,de Becker,1983)与Wahr (1981a,b,c)模型值的相对误差为0.060—0.015。我国的重力潮汐剖面(毛慧琴等,1989)参数与Wahr模型值的相对误差为0.030—0.006,我国1989—1993年的连续重力观测台网(魏望生,喻节林,1995)与DDW模型的潮汐因子最小相差约0.002,这得益于当时我国重力潮汐因子的观测精度从±0.005提高至±0.001左右。随着全球海潮负荷改正模型和气压改正的不断完善,超导重力仪格值的精密测定工作持续展开,早期的T型(大容量杜瓦瓶型)、CT型(紧凑型)和CD型(双球型)超导重力仪周日和半日波潮汐因子的观测精度可达到±0.000 4—±0.000 7,为建立以观测为主的全球实验重力潮汐模型(global experimental model for gravity tides,缩写为GEMGT)创造了条件(Xu et al,2004)。近一二十年里,孙和平等(2000)和韦进等(2012)利用现代观象台型超导重力仪的观测结果来检验DDW和M2001两模型,其结果显示,潮汐因子偏差仅为0.002 4—0.000 6,其中M2和O1波的偏差仅0.002 0—0.000 6。上述成果都仅利用了观测结果与模型的偏差来评价重力仪的潮汐观测能力,并未对全球潮汐模型的潮汐改正能力予以评估,且以观测为主的实验模型只用到了中国一个超导重力站(武汉)。该模型是否适合我国大面积相对重力观测数据的潮汐改正,在理论和实验模型中如何选取潮汐改正模型,这些问题依然没有解决。为研究全球理论潮汐模型对我国重力潮汐改正的精度,并指导我国大面积流动重力观测数据的处理,本文拟选取2016—2018年我国北部的10个gPhone连续重力观测站的数据,在评估其潮汐观测精度的基础上,利用观测的和理论的潮汐改正模型参数构建全球潮汐模型的评价指标,提出符合我国相对重力潮汐改正的最优模型,以期为后续建立我国潮汐改正模型提供初值,也为流动重力数据高精度潮汐改正提供改进方案。
1. 全球潮汐模型
1.1 全球潮汐模型简介
为评价全球潮汐模型在中国区域的适应性,本文首先收集了1961年以来的7个全球潮汐模型,列于表1。
表 1 全球潮汐模型的特征Table 1. Feature statistics of global body tide models来源 长周期波 周日波 半日波 1/3日波 地球模型 Molodensky (1961) - 4 1 - 圆球 Warh (1981a,b,c) 1 19 1 0 1066A,PREM-C,C2 Dehant et al (1999) 1 9 1 1 流体静力弹性地球模型 11 9 1 1 非流体非弹性地球模型 Mathews (2001) - 16 - - PREM 二十世纪中叶,Molodensky (1961)和Wahr (1981a)分别利用圆球、1066A、PREM-C、C2等地球模型反演了潮汐模型,并将潮汐因子表达为全球常数项和相对极其小的纬度依赖项之和的形式。Dehant等(1999)以及Mathews (2001)都使用VLBI技术观测地球自由核章动(free core nutation,缩写为FCN)的周期,并通过该周期改正液核动力学扁率,因而建立了更精确的重力潮汐理论模型。至此模型潮汐因子的精度可达到0.000 01。从所收集到的全球潮汐模型的潮波组成来看,Dehant模型虽然在周日波的波群数量上不及Mathews模型,但其覆盖频段最宽。
1.2 理论潮汐改正模型的正演
观测数据量决定了观测潮汐模型波群(以下简称为观测波群)的分波方式,地球模型决定了理论潮汐模型波群(以下简称为模型波群)。观测波群与模型波群之间会因频率范围的不同而导致两者进行比较时存在差异。通常一个观测波群的频段包含多个理论波群,所以以观测潮汐模型的分波方式为标准,将多个理论波群的潮汐振幅(Tamura,1987)加权平均作为潮汐改正模型的振幅。M2001模型没有长周期、半日波、1/3日波的波群,本文用DDW-NHi模型的相应波群代替。这样全球潮汐模型正演的潮汐改正模型不仅能在频域内计算重力残差矢量(Melchior,1994;Habel,Meurers,2014),还能计算正演时域内的潮汐改正值。
由于我国超导重力仪较少,弹簧重力仪数量居多,因此在数据处理和分析过程中,还大量使用全球潮汐因子平均值进行潮汐改正。本文对1961年以来的7个理论全球潮汐模型中振幅较大的5个潮波的常数项和纬度依赖项进行了对比(表2),结果表明,从Dehant模型开始,潮汐改正模型的精度提高了两个数量级,而且M2001还提供了潮汐改正模型的相位信息,这势必会为当前丰富的地面重力数据带来更精确的潮汐改正值。
表 2 7个全球潮汐模型潮汐因子的常数项和纬度依赖项Table 2. Constant term and latitude-dependent term of the gravimetric amplitude factor in seven global body tide models全球潮汐模型 O1波 P1波 K1波 Q1波 M2波 ${G_0} $ $G_ { \pm }'/10^{-5}$ ${G_0} $ $G_ { \pm }'/10^{-5}$ ${G_0} $ $G_ { \pm }'/10^{-5}$ ${G_0} $ $G_ { \pm }'/10^{-5}$ ${G_0} $ $G_ { \pm }'/10^{-5}$ Molodensky 圆球 1.16 0 1.153 0 1.137 0 1.137 0 1.16 0 1066A 1.152 −600 1.147 −600 1.132 −600 1.152 −600 1.16 −500 Warh PREM-C 1.152 −700 1.147 −700 1.132 −600 1.152 −700 1.16 −500 C2 1.151 −700 1.147 −600 1.132 −600 1.151 −700 1.159 −500 Dehant DDW-He 1.154 24 8 1.147 77 −13 1.147 77 −13 1.132 83 −60 1.160 30 7 DDW-NHi 1.154 24 8 1.149 15 −10 1.134 89 −57 1.154 03 9 1.161 72 10 Mathews M2001 1.154 07 5 1.148 91 5 1.136 1 5 1.154 1 5 - - 注: ${G_0} $ 为潮汐因子的全球常数项,$G_ { \pm }'$ 为纬度依赖项系数,其中M2001模型的${G_0} $ 为原模型常数项的模,$G_ { \pm }'$ 为原模型的周期项。测点周日波的潮汐因子可表示为
$\delta {\text{=}} {G_0} {\text{+}} G'_ { {\text{±} } }\frac{{\sqrt 3 }}{{2\sqrt 2 }}(7{\cos ^2}\theta {\text{-}} 3){\text{,}}$
(1) 半日波的潮汐因子可被表示为
$\delta {\text{=}} {G_0} {\text{+}} G_ {\text{±}} '\frac{{\sqrt 3 }}{2}(7{\cos ^2}\theta {\text{-}} 1){\text{,}}$
(2) This page contains the following errors:
error on line 1 at column 1: Start tag expected, '<' not foundBelow is a rendering of the page up to the first error.
2. 重力资料的精度评价和模型的检验
2.1 重力资料处理及其精度评价
2.1.1 数据预处理和潮汐分析
gPhone重力仪采用了零长弹簧技术,采样频率为1 Hz。在进行重力固体潮观测同时台站还提供了气压小时值数据。本研究对2016—2018年我国10台套该型重力仪的重力固体潮观测数据逐年进行降采样、格值系数校正、预处理、海潮负荷改正、联合气压潮汐分析等处理,建立了观测潮汐改正模型。本研究所用10个观测站的位置如图1所示,具体信息列于表3。
表 3 2016年重力潮汐观测站基本信息Table 3. Basic information of the gravity observatories in 2016序号 台站 省份 台站
类型观测环境
等级台站简介 中误差
/(10−8 m·s−2)1 高台 甘肃 地下室 优秀 位于黑河以北不足2 km处 ±0.637 2 格尔木 青海 山洞 优秀 位于昆仑布尔汉达山北麓山洞中 ±0.706 3 海拉尔 内蒙 山洞 优秀 位于海拉尔盆地与大兴安岭北坡海拉尔河与伊敏河交汇处 ±0.601 4 兰州 甘肃 山洞 优秀 位于兰州东盆地盐场堡以北、白塔山以东、黄河北岸Ⅲ级阶级地的后缘,
距黄河最近距离约1.5 km±0.740 5 牡丹江 黑龙江 山洞 良好 距市区约1.6 km,近年在距离台站500 m处建有成片住宅小区 ±0.718 6 沈阳 辽宁 山洞 良好 位于沈阳市东郊东陵区天柱山西南麓,南距沈抚公路北线约0.4 km ±0.747 7 泰安 山东 山洞 良好 位于泰山南麓正南,0.5 km处有小型水库 ±0.989 8 炭山 宁夏 山洞 良好 位于海原县炭山乡丘陵村,距海原县新区约20 km,台站供电系统时常故障 ±0.891 9 乌加河 内蒙 山洞 优秀 位于内蒙古自治区巴彦淖尔市乌拉特中旗乌加河镇正北1 km处 ±0.838 10 乌什 新疆 山洞 优秀 阿克苏地区乌什县西郊的七女坟旁 ±0.605 本文所选重力站大多分布在我国的中高纬度地区,观测场地与潮汐分析中误差存在一定相关性,由表3可见:观测环境等级为优秀的观测站,中误差为±0.60×10−8 m/s2;观测环境等级为良好的台站中误差虽稍大,也都小于±1×10−8 m/s2。10个观测站的中误差平均值为±0.74×10−8 m/s2,其中乌什重力站最小,约为±0.6×10−8 m/s2。
2.1.2 潮汐分析结果及比较
潮汐因子中误差是评价gPhone重力仪重力潮汐观测能力的一个重要指标。本文将2016—2018年期间的重力固体潮观测数据逐年进行调和分析,平均后统计分析振幅大于10×10−8 m/s2的5个潮波的潮汐因子δ和中误差σ,结果列于表4和表5。
表 4 连续重力观测站的潮汐因子δ和气压导纳值的统计结果Table 4. The result of gravimetric amplitude factor δ and barometric admittance of gravity observatories重力站 O1波 P1波 K1波 M2波 S2波 气压导纳/(10−8 m·s−2·hPa−1) δ σ/10−4 δ σ/10−4 δ σ/10−4 δ σ/10−4 δ σ/10−4 数值 中误差 格尔木 1.152 94 4.0 1.145 22 14.3 1.131 84 9.9 1.165 22 1.5 1.165 10 3.4 −0.269 0.015 兰州 1.157 16 4.7 1.148 77 9.0 1.136 50 3.4 1.158 77 2.4 1.160 37 5.4 −0.319 0.011 高台 1.156 37 3.8 1.150 97 6.7 1.137 49 2.5 1.161 62 1.2 1.163 21 2.7 −0.391 0.008 海拉尔 1.155 19 5.9 1.157 16 11.2 1.138 65 4.1 1.167 74 4.2 1.171 67 9.4 −0.371 0.010 牡丹江 1.153 84 4.7 1.151 92 8.4 1.137 44 3.2 1.163 77 3.4 1.167 07 7.5 −0.347 0.009 沈阳 1.153 70 7.3 1.150 00 13.1 1.137 17 4.9 1.159 45 2.8 1.159 77 6.2 −0.356 0.010 泰安 1.158 32 11.3 1.153 09 20.6 1.139 41 7.7 1.155 73 2.4 1.156 20 5.6 −0.288 0.018 炭山 1.154 74 7.8 1.145 53 16.9 1.133 36 6.4 1.161 86 2.4 1.161 45 5.9 −0.312 0.024 乌加河 1.155 54 3.1 1.150 56 5.6 1.137 52 2.1 1.157 09 1.3 1.159 06 3.1 −0.299 0.010 乌什 1.155 14 3.9 1.151 40 7.0 1.137 26 2.6 1.160 97 1.4 1.161 00 3.2 −0.430 0.010 注:δ为潮波的潮汐因子,σ为中误差。 表 5 乌什重力观测站的潮汐分析结果Table 5. The tidal analysis results of Wushi gravity observatory起始频率
/cpd结束频率
/cpd潮波名称 潮汐振幅
/(10−9 m·s−2)潮汐因子 相位滞后/° 数值 中误差 数值 中误差 0.501 37 0.911 39 Q1 58.938 1 1.154 66 0.001 80 −0.097 0.09 0.911 39 0.947 99 O1 307.830 6 1.155 14 0.000 39 0.020 0.019 0.947 99 0.981 85 M1 24.209 7 1.504 93 0.004 20 −0.427 0.160 0.981 86 0.998 63 P1 143.232 9 1.151 40 0.000 70 −0.060 0.035 0.998 63 1.001 37 S1 3.387 2 1.209 77 0.042 58 16.412 2.319 1.001 37 1.004 11 K1 432.929 1 1.137 26 0.000 26 0.007 0.014 1.004 11 1.006 85 ψ1 3.387 1 1.260 70 0.030 07 −0.923 1.368 1.006 85 1.023 62 φ1 6.164 9 1.169 38 0.016 14 −0.622 0.791 1.023 62 1.057 49 J1 24.208 8 1.153 00 0.004 61 −0.174 0.229 1.057 49 1.470 24 OO1 13.246 6 1.162 69 0.011 38 0.117 0.561 1.470 24 1.880 26 2N2 13.001 1 1.161 73 0.003 27 0.201 0.161 1.880 27 1.914 13 N2 81.411 6 1.159 33 0.000 71 −0.046 0.035 1.914 13 1.950 42 M2 425.214 3 1.160 97 0.000 14 0.013 0.007 1.950 42 1.984 28 L2 12.018 8 1.157 81 0.006 06 −0.428 0.300 1.984 28 2.002 74 S2 197.832 1 1.161 00 0.000 32 −0.289 0.020 2.002 74 2.451 94 K2 53.780 9 1.158 76 0.001 54 −0.137 0.076 2.451 94 7.000 00 M3M6 6.287 4 1.073 14 0.006 61 0.109 0.353 由表4可见:10套仪器中除泰安和海拉尔的部分波外,其它波的潮汐因子中误差均在±0.000 4—±0.000 7附近。潮汐观测精度虽不及OSG型超导重力仪的精度,但已经与二十世纪八九十年代的超导重力仪结果相当。只是gPhone重力仪测定的气压导纳平均值为(−0.319±0.028)×10−8 m·s−2·hPa−1,精度较超导重力仪结果低一个量级(表6)。此外,10个重力观测站中,乌什重力站的潮汐分析结果精度最高,其潮汐分析结果(表5)呈现潮波振幅大、精度高和振幅小、精度低的特征,其中振幅最大的M2波的中误差可达±0.000 14,振幅最小的ψ1波仅为±0.003 27。这表明,最优gPhone重力仪(乌什重力站)的潮汐观测精度虽只有OSG型超导重力仪观测精度的0.5倍,然已高于早期超导重力仪0.5个数量级(表6)。
表 6 各型重力仪潮汐观测精度和稳定度分析Table 6. Accuracy and stability analysis of tidal observatory for different gravimeters重力仪
类型数据时段 M2波潮汐因子
中误差中误差
/(10−8 m·s−2)导纳中误差
/(10−8 m·s−2·hPa−1)稳定度 来源 LRC 1973—1990 0.002 00—0.005 00 - - - Melchior和de Becher (1983) GEO 1971—1990 0.001 00—0.002 00 - - - GWR 1975—1988 0.000 08 - - - LRC 1983—1986 ≤0.005 00 - - 0.005 0 毛慧琴等(1989) GS 1989—1993 0.000 30—0.003 00 - - - 魏望生和喻节林(1995) DZW 0.000 80—0.004 00 - - - LRC 0.000 40—0.002 00 - - - GEO 0.001 00—0.005 00 - - - T,CT,CD 1989—2000 - ±0.05—±0.70 ±0.000 4—0.003 0 0.001 5 Ducarme et al (2002) ASK - - ±0.010 0 - LRC 2000—2001 0.000 51 ±0.20 ±0.022 0 - 孙和平等(2002) OSG 0.000 06 ±0.02 ±0.002 0 - T,CT,CD 1997—2001 0.000 70 - ±0.000 4—0.003 0 0.001 0 Xu et al (2004) OSG 2009—2010 0.000 04 ±1.10—±1.40 ±0.000 5 - 韦进等(2012) 2.1.3 潮汐观测的稳定性检验
重力仪的潮汐观测能力还可以通过潮汐稳定度(Ducarme et al,2002)来评估。本文将10套仪器2016—2018年各潮波的年潮汐因子与3年均值的差作为潮汐稳定度(下文简称为稳定度)来评估重力仪的潮汐观测能力。
图2统计了10个重力站4个潮波的稳定度,可见:除其中7个观测站的P1波(图2c)外,其它潮波的稳定度均小于0.001 5,这一结果已经远高于流动重力仪组网观测的重力潮汐剖面各站点的稳定度(≤0.005 0)(毛慧琴等,1989);7个重力站4个潮波的稳定度均达到了早期超导重力仪的观测水平(≤0.001 5)。
表6给出了二十世纪七十年以来重力仪潮汐观测能力的指标统计表。可见:gPhone重力仪周日和半日波的潮汐观测精度优于二十世纪八九十年代的LCR (LaCoste&Romberg)、GEO (大地型)、DZW (中国地震局地震研究所)、GS (阿斯卡尼亚公司一款重力仪)、ASK (阿斯卡尼亚公司一款重力仪)等各型弹簧重力仪约0.5—1个数量级,其中精度高的与早期的T型(大容量杜瓦瓶型)、CT型(紧凑型)和CD型(双球型)的相当,但与现代OSG型超导重力仪相比要低0.5—1个数量级。气压负荷改正参数的精度方面,gPhone可以解算出精度为±0.009×10−8—±0.028×10−8 m·s−2·hPa−1的气压导纳值,较超导重力仪结果低一个数量级。
全球潮汐模型关注的频段是gPhone重力仪的优势频段(周日和半日波),且Xu等(2004)在利用全球超导重力仪建立全球实验潮汐模型时就指出,潮汐因子中误差优于0.000 7即可进行全球潮汐模型的建立工作。本文据此选取重力仪来评估全球潮汐模型的精度。
2.2 观测和全球潮汐模型的比较
上述10个重力站的潮汐观测精度均达到了建立实验潮汐模型的标准,因此利用2016—2018年各年份的观测模型和表2中所列的7个全球潮汐模型进行比较分析,建立两者间的相对误差、单一波群的均方根(root mean square,缩写为RMS)、所有波群的和方根(root sum square,缩写为RSS)和纬度依赖关系等指标进行全球潮汐模型的精度评定。
2.2.1 基本方法
gPhone重力仪的优势频段为周日和半日波,且潮波振幅越大、精度越高。本研究采用均方根(RMS)来衡量理论模型与重力网观测模型的某潮波偏差和方根RSS来综合评估全球潮汐模型的精度(李大炜等,2012),具体表达式如下:
$ {\rm{RMS}} {\text{=}} {\left\{ {\frac{1}{K}\sum\limits_{n {\text{=}} 1}^K {[{A_0^2}{{\!\!\!\!{\text{(}}\!n\!{\text{)}}\!\!\!\!}}{{\!\!\!\!{\text{(}}\!{\delta _0}\!\!\!\!{\text{(}}\!n\!{\text{)}}\!\!\!\!\cos {\varphi _0} {\text{-}} {\delta _m}\!\!\!\!{\text{(}}\!n\!{\text{)}}\!\!\!\!\cos {\varphi _m}\!{\text{)}}\!\!\!\!}^2} {\text{+}} {A_0^2}{{\!\!\!\!{\text{(}}\!n\!{\text{)}}\!\!\!\!}}{{\!\!\!\!{\text{(}}\!{\delta _0}\!\!\!\!{\text{(}}\!n\!{\text{)}}\!\!\!\!\sin {\varphi _0} {\text{-}} {\delta _m}\!\!\!\!{\text{(}}\!n\!{\text{)}}\!\!\!\!\sin {\varphi _m}\!{\text{)}}\!\!\!\!}^2}]} } \right\}^{1/2}}{\text{,}} $
(3) $\begin{split} {\rm{RSS}} {\text{=}}& \left\{ {\frac{1}{K}\sum\limits_{i {\text{=}} 1}^I {\sum\limits_{n {\text{=}} 1}^K {\left[ {{{\!\!\!\!{\text{(}}\!{\delta _{0{\text{,}}\!\!\!\!\!i}}\!\!\!\!{\text{(}}\!n\!{\text{)}}\!\!\!\!{A_{0{\text{,}}\!\!\!\!\!i}}\!\!\!\!{\text{(}}\!n\!{\text{)}}\!\!\!\!\cos {\varphi _{0{\text{,}}\!\!\!\!\!i}} {\text{-}} {\delta _{m{\text{,}}\!\!\!\!\!i}}\!\!\!\!{\text{(}}\!n\!{\text{)}}\!\!\!\!{A_{0{\text{,}}\!\!\!\!\!i}}\!\!\!\!{\text{(}}\!n\!{\text{)}}\!\!\!\!\cos {\varphi _{m{\text{,}}\!\!\!\!\!i}}\!{\text{)}}\!\!\!\!}^2} {\text{+}} } \right.} } } \right.\\& {\left. {\left. {{{\!\!\!\!{\text{(}}\!{\delta _{0{\text{,}}\!\!\!\!\!i}}\!\!\!\!{\text{(}}\!n\!{\text{)}}\!\!\!\!{A_{0{\text{,}}\!\!\!\!\!i}}\!\!\!\!{\text{(}}\!n\!{\text{)}}\!\!\!\!\sin {\varphi _{0{\text{,}}\!\!\!\!\!i}} {\text{-}} {\delta _{m{\text{,}}\!\!\!\!\!i}}\!\!\!\!{\text{(}}\!n\!{\text{)}}\!\!\!\!{A_{0{\text{,}}\!\!\!\!\!i}}\!\!\!\!{\text{(}}\!n\!{\text{)}}\!\!\!\!\sin {\varphi _{m{\text{,}}\!\!\!\!\!i}}\!{\text{)}}\!\!\!\!}^2}} \right]} \right\}^{1/2}}{\text{,}} \end{split}$
(4) 式中:δ0,A0和φ0分别为重力站的潮汐因子、理论振幅、理论相位(通常为0°);δm和φm分别为模型在重力站点经过气压、海潮改正后的潮汐因子和相位,K为重力站的个数,n为重力站的序号,m表示观测结果,I为潮波分波数,本文为5 (M2,K1,O1,P1,Q1,共计5个)。理论潮汐振幅分别从Tamura (1987)的潮波表和郗钦文的潮波表(郗钦文,侯天航,1987)中选一。表7为利用7个理论模型计算的RMS和RSS。由表7可见:除Molodensky均质圆球模型外,考虑地球扁率的M2001和DDW-NHi模型的RSS更小,其中M2001结果与观测数据吻合得更好,这表明考虑地球扁率的模型更能反映我国潮汐的观测规律,采用郗钦文潮波表(郗钦文,侯天航,1987)计算的RSS优于采用Tamura潮汐表所得的RSS,由此可见DDW-NHi和M2001模型的精度高于其它模型。
表 7 7个全球潮汐模型计算所得重力残差的均方根RMS及和方根RSSTable 7. RMS and RSS of the residual gravity corrected by seven global body tide models潮汐模型 地球模型 RMS/(10−8 m·s−2) RSS/(10−8 m·s−2) M2波 K1波 O1波 P1波 Q1波 Tamura
潮波表郗钦文
潮波表Molodensky 均质圆球 0.143 9 0.168 7 0.173 4 0.064 8 0.047 8 0.292 8 0.291 2 Warh 1066A 0.190 9 0.179 8 0.169 3 0.068 9 0.039 9 0.322 1 0.322 0 PREM-C 0.190 9 0.179 8 0.169 4 0.069 3 0.040 0 0.322 3 0.320 6 C2 0.199 9 0.179 8 0.171 2 0.068 9 0.039 9 0.328 5 0.329 0 DDW DDW-He 0.143 6 0.227 5 0.167 1 0.066 6 0.053 7 0.328 0 0.328 1 DDW-NHi 0.144 5 0.171 3 0.167 1 0.065 2 0.040 2 0.289 8 0.289 6 Mathews M2001 0.144 5 0.168 3 0.167 2 0.065 2 0.040 2 0.288 1 0.287 9 两个全球潮汐模型与观测站潮汐因子的相对误差统计结果(表8)显示,DDW-NHi模型下各重力站的M2,O1,K1和P1波潮汐因子的相对误差均值分别为0.16%,0.17%,0.35%,0.28%,M2001模型下的各潮波潮汐因子的相对误差均值分别为0.15%,0.15%,0.37%,0.18%,这表明DDW-NHi和M2001模型下各重力站各潮波潮汐因子的相对误差为0.15%—0.37%,一致性高,其中高台重力站M2波潮汐因子的偏差在0.000 15附近,达到了早期超导重力仪的潮汐观测稳定度(Xu et al,2004;Habel,Meurers,2014)。从表8还发现,M2001模型的RMS比DDW-NHi更小,这与累计残差振幅结果(表7)一致。
表 8 两个全球潮汐模型与观测站的潮汐因子的相对误差Table 8. The relative errors of the gravimetric amplitude factors between DDW-NHi and M2001 models and the observatories台站 M2波潮汐因子相对误差 O1波潮汐因子相对误差 P1波潮汐因子相对误差 K1波潮汐因子相对误差 RMS/(10−8 m·s−2) DDW-NHi M2001 DDW-NHi M2001 DDW-NHi M2001 DDW-NHi M2001 DDW-NHi M2001 格尔木 0.251 3% 0.251 3% 0.178 2% 0.150 4% 0.253 3% 0.247 4% 0.259 5% 0.459 4% 0.120 6 0.141 0 高台 0.015 2% 0.015 2% 0.224 4% 0.253 1% 0.187 8% 0.189 1% 0.312 9% 0.094 1% 0.087 5 0.068 3 海拉尔 0.015 3% 0.495 8% 0.062 1% 0.093 9% 0.815 2% 0.800 9% 0.574 1% 0.290 9% 0.223 8 0.198 0 兰州 0.015 4% 0.211 4% 0.354 1% 0.381 8% 0.269 4% 0.275 8% 0.291 5% 0.093 8% 0.114 6 0.105 1 牡丹江 0.015 5% 0.200 8% 0.188 0% 0.157 6% 0.582 2% 0.575 4% 0.244 1% 0.008 6% 0.227 8 0.220 6 乌什 0.015 6% 0.236 3% 0.025 6% 0.003 9% 0.168 3% 0.165 9% 0.246 3% 0.011 9% 0.107 8 0.095 5 2.2.2 潮汐振幅的纬度依赖关系
在全球潮汐模型中,潮汐因子通常被表示为全球常数项与纬度依赖项之和,潮汐因子的值与重力站纬度有关。将10个重力站2016—2018年逐年的潮汐参数和3个RSS较小的全球潮汐模型的潮汐参数按纬度排列,如图3所示。
由图3a可见,潮汐振幅与纬度有明显的相关性,其中:M2波和O1波的理论振幅与观测振幅差异不大;DDW-NHi以及M2001模型的K1波一致性较好,而与Molodensky模型相差1×10−8—2×10−8 m/s2,所以观测振幅与DDW-NHi或M2001模型一致性更好。
从潮汐因子与纬度的相关性(图3b)可见,本文所选站点的纬度与潮汐因子的相关性不明显,观测值分居模型的两侧,观测值与模型值相差0.000 15—0.005 00。
综合RMS、RSS和潮汐振幅的纬度依赖关系表明,DDW-NHi和M2001模型的RSS值相对较小,处于0.28×10−8—0.29×10−8 m/s2之间(表7)。在潮汐振幅与纬度的相关性中,Molodensky模型的K1波在高纬度地区(35°N—50°N)明显偏离另外两个模型。因此本文初步认定DDW-NHi和M2001模型可被我国相对重力仪用于潮汐改正。
2.2.3 理论和观测潮汐改正模型的比较
我国流动重力观测站数目远多于连续重力观测站,流动重力数据采用的还是Molodensky均质圆形地球确定的潮汐改正模型。为评估该模型和本文筛选模型时域范围内的改正特征,依然将Molodensky模型纳入评价范围。因此本文利用3个理论的和1个观测的潮汐改正模型对观测精度最优的乌什重力站2016年的观测数据(图4a)进行时域改正,比较改正后的重力残差时间序列和功率谱密度,结果如图4所示。
图 4 潮汐模型潮汐改正后的重力残差时间序列及其功率谱密度(a) 乌什重力站固体潮预处理数据;(b) Molodensky模型潮汐改正的重力残差;(c) DDW-NHi模型潮汐改正的重力残差;(d) M2001模型潮汐改正的重力残差;(e) 观测模型潮汐改正的重力残差;(f) 重力残差的功率谱密度Figure 4. Residual gravity and power spectral density corrected by several tidal models(a) The preprocessing gravity data from Wushi observatory;(b) The residual gravity corrected by Molodensky model;(c) The residual gravity corrected by DDW-NHi model;(d) The residual gravity corrected by M2001 model;(e) The residual gravity corrected by observed model;(f) The power spectral density of the residual gravity由图4可见:观测模型的重力残差最小(图4e),其日均值中误差仅为±0.1×10−8—±0.5×10−8 m/s2,且功率谱密度最大能量不到−110 dB (1 dB=10lg (10−8 m·s−2)2/Hz),几乎被淹没在残差噪声中(图4f);3个理论模型中DDW-NHi的重力残差日均值中误差最小,仅为±0.4×10−8—±1.0×10−8 m/s2,优于Molodensky模型的中误差±0.5×10−8—±1.5×10−8 m/s2,也优于M2001模型的中误差±0.7×10−8—±1.4×10−8 m/s2。
对理论潮汐模型的潮汐改正精度进行评价时,RMS和RSS指标未考虑小振幅潮波。而利用它们在时域内进行潮汐改正,则可以分析这些小振幅潮波对潮汐改正的影响。理论模型中,DDW-NHi和M2001模型的周日波频段上O1波,P1波和K1波的残差能量相当(图5a),但与Molodensky模型的K1波存在较大差异,这与潮汐振幅与纬度的相关性结果(图3e)一致。另一方面,Molodensky模型的K1波功率谱密度大于其它两模型,应与该模型K1波的潮汐因子不准确有关。
由于M2001模型的M2波模型参数采用的是其它模型的,带入时域改正后存在−90 dB的重力残差信号(图5b),这可能是M2波模型参数相位为0°所致,因此这也是M2001模型的改正结果不及DDW-NHi模型的原因之一。此外,模型与观测的差异应该是模型中未考虑台站观测场地、观测环境等因素的影响所导致。
综上所述,DDW-NHi模型是本文统计的理论模型中最优者,而传统的Molodensky潮汐改正模型的改正精度仅为±0.5×10−8—±1.5×10−8 m/s2。
3. 讨论与结论
在全球的7个潮汐模型中,DDW-NHi模型的潮汐改正结果精度最高,能达到±0.4×10−8—±1.0×10−8 m/s2,而传统的Molodensky模型的潮汐改正结果只能达到±0.5×10−8—±1.5×10−8 m/s2。前文分析主要集中于不同地球模型的改正精度评估,然而除了地球模型之外,影响全球潮汐模型改正精度的因素还包括计算精度时所使用的海潮模型、仪器多年运行观测系统的格值系数变化等。为讨论本文结果的有效性,下面将进一步讨论上述两个因素对全球潮汐模型改正精度评估的影响。
全球潮汐模型仅含地球岩石圈及其内部圈层的影响,而实际观测数据中富含诸如海潮、气压等地球外部圈层的影响。近年来国际上推出多个高精度全球海洋潮汐模型和区域海潮模型,经过验潮站、全球导航卫星系统(global navigation satellite system,缩写为GNSS)、重力等手段的检验获知联合近海海潮模型方能获得更高精度的海潮改正。
3.1 区域海潮改正模型对利用全球模型潮汐改正的影响
本文利用重力残差矢量(Melchior,1994,Hábel,Meurers,2014)分潮波比较Nao.99b全球海潮模型、Nao.99b+osu.chinasea.2010组合以及Nao.99b+naoregional.1999组合这3种海潮模型对DDW-NHi模型精度的影响特征,结果如图6所示。由于DDW-NHi模型并不包含相位的影响,因此在进行重力固体潮频域改正时也不考虑相位。
图 6 海潮负荷改正前后观测的与模型的重力残差振幅比较图中P代表进行了气压改正;+N代表进行了Nao.99b海潮负荷改正;+NC代表进行了Nao.99b和osu.chinasea.2010海潮负荷改正;+NN代表进行了Nao.99b和naoregional.1999海潮负荷改正Figure 6. Comparison of the residual gravity amplitude before and after corrected by ocean loading models and global body tide modelP represents the gravimetric amplitude factors corrected by atmospheric pressure;+N represents those corrected by global ocean model of Nao.99b;+NC represent the those corrected by the composite ocean model of Nao. 99b and osu.chinasea.2010;+NN represent those corrected by the composite ocean model of Nao.99b and naoregional.1999未进行海潮改正时(图6中的红色线柱),各潮波重力残差振幅和各潮波累加值在0.1×10−8—2×10−8 m/s2之间,并呈现东部振幅大、西部小的特点;东部振幅最大为2×10−8 m/s2,西部最小为0.5×10−8 m/s2。海潮负荷改正后,所有潮波的重力残差振幅显著减小,分潮波的海潮负荷改正量最大为1.5×10−8 m/s2,最小只有0.3×10−8 m/s2,改正量为改正前重力残差振幅的70%—80%,这一结果与其它内陆台站的改正量(周江存等,2009)基本一致。由于海潮负荷改正量级已经远大于表7中单一潮波的重力残差振幅,因此海潮负荷改正是本文必不可少的工作。
由图6还可以看到,不同组合模型之间的重力残差振幅差异不足0.1×10−8 m/s2,该量级小于表7中单一潮波重力残差振幅,因此本文认为区域海潮模型的差异不会对全球潮汐模型的改正产生影响。
3.2 全球海潮模型对利用全球模型潮汐改正的影响讨论
为进一步讨论全球海潮模型对潮汐改正的影响,本文对海潮负荷时序改正前后的潮汐因子与DDW-NHi模型之间的差异进行对比,结果如图7所示。可见:所有东部台站的海潮负荷改正量均较大,西部台站均较小,这符合我国东部重力站受太平洋海潮影响强于西部的规律。从改正量级看,牡丹江重力站O1波的改正量最大,可达0.020 0,平均值也有0.003 0,该量级与观测系统的潮汐稳定度相当,高于潮汐因子的观测精度,说明海潮负荷改正后的潮汐因子大小均趋于DDW-NHi全球潮汐模型值,海潮负荷改正在研究非常必要。
图 7 海潮改正对潮汐因子的影响图中P代表进行了气压改正的潮汐因子;DUT10,EOT11a,FES2004,GOT4.7,HAM11a,Nao.99b和TPXO7.2分别代表进行了气压和对应全球海潮模型改正的潮汐因子;DDW-NHi代表了DDW-NHi理论潮汐因子。横坐标台站按经度由东向西排列Figure 7. Impaction of the gravimetric amplitude factors corrected by seven global ocean modelsP represents the gravimetric amplitude factors corrected by atmospheric pressure;DUT10,EOT11a,FES2004,GOT4.7,HAM11a,Nao.99b and TPXO7.2 represent those corrected by atmospheric pressure and corresponding global ocean tide model respectively;DDW-NHi represents those in DDW-NHi. The stations are arranged from east to west according to their longitude比较不同海潮负荷改正模型之间的差异可知,图7中的7个全球海潮负荷改正结果的差异较小,为0.000 5—0.002 0,小于改正前后的潮汐因子差异0.020 0—0.003 0和观测系统稳定度0.001 5,与潮汐因子的观测精度0.000 7相当。因此,在进行全球潮汐模型评估时须进行海潮负荷改正,而不同海潮负荷改正模型之间的差异对结果影响并不明显。
3.3 格值系数对模型检验的影响
相对重力仪格值系数偏差会直接影响观测潮汐改正模型的准确性(Francis et al,1998;Hábel,Meurers,2014),gPhone重力仪的格值系数的范围是0.999 9—1.019 0 (刘子维等,2011;张锐等,2011)。本研究选择青海格尔木重力站2016—2018年逐年观测数据讨论潮汐因子、残差中误差与格值系数的关系,结果如图8所示。
图 8 格尔木观测站潮汐因子和调和分析中误差随格值系数的变化(a) M2波;(b) O1波;(c) K1波;(d) P1波;(e) Q1波;(f) 中误差随格值系数的变化Figure 8. Varaition of gravimetric amplitude factor and RMS of the harmonic analysis with the scale factor for Golmud observatory(a) M2 wave;(b) O1 wave;(c) K1 wave;(d) P1 wave;(e) Q1 wave;(f) The RMS of the residual gravity以DDW-NHi模型为潮汐改正模型,利用重力残差(图8f)作为约束测定格值系数。从图8a-e可见:观测潮汐因子与理论仅相差0.001—0.002。如果不进行这样的标定(即格值系数为1),潮汐因子偏差可能大于0.005,该量级已远超当前gPhone重力仪的稳定度和潮汐观测精度。因此进行观测系统的格值系数改正是必不可少的。
3.4 结论
本文比较了7个常用的全球潮汐模型,并利用10个安装gPhone重力仪的重力站超过3年的观测资料对模型进行了检验。gPhone重力仪的周日和半日波潮汐观测能力虽然比现代OSG型超导重力仪低0.5—1.0个数量级,但是依然可达到早期超导重力仪的水平,这样的精度可以对全球潮汐模型精度进行评估。对潮汐观测和理论模型的均方根、和方根、潮汐因子相对误差、潮汐因子的纬度依赖关系和重力残差及其功率谱密度的评估结果显示DDW-NHi的全球潮汐模型能够提供最优的重力残差精度为±0.4×10−8—±1.0×10−8 m/s2,比大规模流动重力观测常用的Molodensky模型的精度±0.5×10−8—±1.5×10−8 m/s2高出1×10−8—2×10−8 m/s2。
在讨论海潮负荷和格值系数对潮汐模型精度评估影响时,仅进行全球潮汐模型的改正而未进行海潮和格值系数校正引起的潮汐因子偏差分别可达0.003 0—0.020 0和0.001—0.005,这一量级已经远大于gPhone重力仪的潮汐观测稳定度和潮汐因子的观测精度。不同海潮模型改正的影响量级(0.000 5—0.002 0)与潮汐因子的观测精度相当,因此可忽略其对模型精度评估的影响。
本文所选观测站只有10个,且多分布在我国北方中高纬度地区;无南方台站特别是沿海台站参与本次评估,难免会使得结果仅具有区域性规律。从观测资料来看,南方和沿海地区的潮湿使重力观测系统在长期连续观测过程中的故障率较高,使用现有分析方法和近期数据无法达到评估模型所需的精度。如何筛选此类观测数据、研制新的观测精度评估方法成为解决此类问题的一个重要内容。在进行潮汐改正时,气压负荷影响不容忽视,而并非所有重力站目前均具备同址同采样率的气压观测条件。为弥补站点辅助测量数据的缺失可探讨国际已经发布的全球气压模型的潮汐改正可行性。此外,本文还讨论了格值系数对模型评估的影响,认为需要研制出高效的连续重力观测数据标定方法,使得标定数据达到潮汐模型评估和模型建立所需要的精度。总之,如何科学产出符合全球潮汐模型评估精度的重力潮汐数据是本文工作特别是建立潮汐模型的根本所在。
感谢国家重力台网中心提供10套高精度gPhone重力仪超过3年的重力固体潮数据,感谢审稿专家提出的宝贵意见和建议。
-
图 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