2) 中国太原 030025 山西省太原大陆裂谷动力学国家野外科学观测研究站
2) National Continental Rift Valley Dynamics Observatory of Taiyuan, Shanxi Province, Taiyuan 030025, China
国家地震烈度速报与预警工程山西子项目现已建成由81个基准站、154个基本站及861个一般站组成的山西地震预警台网。3类站点能实时记录地表运动观测数据,是地震监测和预警的首要环节,传感器作为台站重要的专业设备之一,其性能和状态直接决定了台站观测系统的数据质量(包文超等,2023)。因此,做好台站设备运维工作,确保仪器正常运行是地震监测预警的前提。然而,预警工程新建、改建台站数量庞大,在日常工作中,通过人工浏览波形不但不易识别仪器故障,而且难以及时发现异常台站。对设备进行脉冲标定也可检验其运行是否正常,但依据预警设备标定方案:每年对地震计进行1次脉冲标定、对加速度计进行2次方波标定,可见利用该方式来不能及时检验仪器性能。因此,运用一种新的方法定期进行数据质量监控势在必行。
目前,国内外研究机构已经利用统计方法分析台站噪声功率谱密度值在指定时间段的分布,来揭示台站背景噪声的变化特征,从而发现台站的异常情况。在国外,GSN、ANSS等地震监测台网基于McNamara等(2004)提出的利用功率谱密度概率函数(PDF)方法进行台站数据质量监测;Sleeman等(2007)利用监测实时波形地噪声加速度功率谱密度(PSD)值方法,在VEBSN地震台网实现了对宽频带地震计的数据质量监测。在国内,刘旭宙等(2018)在相同观测条件下,通过对比仪器噪声概率密度函数形态,研究了不同仪器的观测性能;江苏、福建、广东、黑龙江等省地震局(王俊等,2009;徐嘉隽等,2010;陈建涛等,2019;胡宝慧等,2022)已经实现了利用PDF方法对地震监测仪器运行状态的检测,开展了波形数据质量的实时监控工作。
地震预警工程山西子项目建设台站数量众多,而台站波形数据质量监控又是地震监测预警的重要任务。因此,以3类台站记录的波形数据为基础,采用计算概率功率谱密度(PPSD)方法(Beyreuther et al,2010),计算台站记录噪声水平,生成PDF及时频分布图,同时利用PDF图对比分析部分台站波形异常与正常时的特点。该方法与PDF方法原理相同,在计算时无需去除地震、爆破、标定、仪器故障等突发事件,并将事件包含在内作为地噪声的一部分进行概率统计处理。这些信号是有用的,其概率值的大小可以用来评价台站观测系统的运行状态和数据质量,出现的频次、概率值能影响环境地噪声水平,进而影响台站观测系统的运行状态和数据质量。
1 山西地震预警台网概述山西地震预警台网新建、改建3类台站共计1 096个。根据目标任务的不同,3类台站配置了不同类型的传感器。其中,基准站配置宽频带地震计及加速度计,基本站配置加速度计,一般站配置烈度计。新建基准站和基本站观测数据采用8M专线传输方式;改造基准站数据传输使用原有20M PTN专线,汇聚节点到省预警中心数据传输采用30M专线;一般站主要采用无线网络方式进行数据传输(丁大业等,2023)。山西地震预警台网实现了测震、强震、烈度台站资源的充分融合,形成一个新型地震监测台网,重点预警区平均台站间距9.6 km、一般预警区平均台站间距42.06 km。山西地震预警台网台站分布见图 1。
对台站波形数据进行去均值、去长周期线性趋势,以消除长周期功率谱估计偏差。对经处理信号进行傅里叶变换,得
$ v_0(\omega)=\int_{-\infty}^{\infty} v_0(t) \mathrm{e}^{-i \omega t} \mathrm{~d} t $ | (1) |
对信号进行傅里叶变换平方,得到功率谱,公式如下
$ P(f)=\frac{1}{T_{\mathrm{d}}} E\left[v_0^2(f)\right] $ | (2) |
对信号做傅里叶变换得到振幅谱,根据式(2)可得功率谱,公式如下
$ P_v(f)=\frac{2}{N \Delta t}|v(f)|^2 $ | (3) |
对式(3)进行仪器响应校正,扣除仪器响应得到地噪声物理量值,公式如下
$ \operatorname{PSD}_v(f)=\frac{P_v(f)}{|H(f)|^2} $ | (4) |
将速度功率谱转换为加速度功率谱,公式如下
$ \operatorname{PSD}_a(f)=4 \pi^2 f^2 \operatorname{PSD}_v(f) $ | (5) |
在高频处,功率谱锯齿起伏严重,难以查看到其真实功率谱密度值。为得到PSD,在频率对数坐标中呈等间隔采样,需采取1/3倍频程积分做平滑处理,公式如下
$ \operatorname{PSD}_a\left(f_c\right)=\frac{1}{n} \sum_{f=f_1}^{f_{h}} \operatorname{PSD}_a(f) $ | (6) |
其中:中心频率fc以1/9倍频程为增加步长,则下一个中心频率fc = fc×21/9;fl = fc×2-1/6,为低频拐角频率;fh = fc×21/6,为高频拐角频率;n为介于两者间离散频率f的个数;PSDa(fc)为中心频率fc的加速度功率谱密度在fl和fc之间的平均值。
每个中心频率的PSD概率密度函数为
$ P_{\mathrm{PSD}}\left(f_{\mathrm{c}}\right)=N_{p_{f_{\mathrm{c}}}} / N_{f_{\mathrm{c}}} $ | (7) |
式中,Nfc为fc频点记录段总数;Npfc为fc频点的PSD值落在某PSD取值范围内的记录段个数。
2.2 技术路线采用Python语言读取连续波形文件、解析数据文件、进行PPSD计算、输出保存图文结果等。Python是开源的、解释型高级编程语言,语法简单易学,应用领域广泛,具有丰富、强大的开源库。项目使用针对地震领域研发的Obspy库,具备数据处理、格式转换、地震信号处理等常用功能。
PPSD计算程序主要包含数据处理及绘图2个部分,在运算过程中需要导入所需模块和库。其中,数据处理部分完成仪器信息文件的解析、噪声水平的计算,绘图部分用来完成PDF图及时频分布图的绘制。从Jopens6.0系统对山西预警台网1 096个台站24小时seed格式连续波形数据进行归档,并从AWS流服务器Web上导出省内所有台站的dataless文件,利用Obspy库对台站波形数据进行批处理计算,生成每个台站3个分量1天的PDF图及频谱分布图,用于反映台站地噪声水平及仪器工作状态,进而实现台站波形数据的质量监控。程序计算流程见图 2。
山西地震预警台网基准站、基本站及一般站的测震设备配置采取混搭方式,地震计型号主要是GL-CS60、JS-60,加速度计型号为JS-A2、TDA-33M,烈度计型号是MI3000及REMOS-SIT4,数据采集器型号是TDE-324CI/FI、HG-D3/D6。考虑到不同仪器型号搭配差异对噪声数据的影响,选取6个代表性台站(表 1)24小时连续波形数据进行PPSD计算。受文章篇幅所限,文中仅列出所选台站Z分向PDF图及相应时频分布图(图 3—图 5),图中00通道接地震计,20通道接加速度计,40通道接烈度计。每个PDF图由PSD概率频谱及连续数据图组成(时间轴是UTC时间,数据时间为UTC+8);右侧色标表示PSD概率对应颜色,概率值由紫色到红色逐渐升高,PSD的固有特性概率高,随机性概率低;黑色平均线越靠近低噪声模型NLNM线,表示该台址环境地噪声水平越低;绿色横条代表连续波形数据,蓝条代表PPSD计算所用数据。时频分布图反映了不同时段仪器记录噪声水平在各频段的分布。
由台站的PDF图(图 3—图 5左侧)可见,深色部分表示台站正常噪声水平所在区域,较浅的紫色线条表示台站记录的异常噪声,如地震、外界扰动等,属于小概率事件,黑色平均线代表该台站平均地噪声水平,由每个频点最大概率点连接而成,反映了台基噪声的固有特性。
对3类仪器记录的噪声水平(图 3—图 5右侧)进行分析,可知基准站地震计(00通道)记录的台基噪声水平优于加速度计(20通道)。这是因为2种仪器的结构、性能上的差异导致加速度计自噪声较大,而地震计记录的噪声水平介于高噪声模型NHNM及低噪声模型NLNM之间,其PSD值介于-170— -120 dB,噪声来源主要为10 Hz以上的高频干扰,时间段集中在8:00—20:00,多与人为活动有关;加速度计功率谱值介于-120— -100 dB,记录的背景噪声仅在大于0.1 Hz的高频段位于高噪声模型NHNM以下,可反映台基噪声水平,低于0.1 Hz的长周期噪声分布在全天各个时段,与受气压、温度等环境因素影响有关(杨千里等,2019),安装加速度计防护罩后PSD值应会有所降低;基本站与基准站加速度计记录的噪声水平基本一致,但PSD曲线在10 Hz以上的高频部分比较离散,且频谱分布呈现出时间相关性,在8:00—20:00噪声水平较高,主要因为基本站位于城镇地区,受人为活动影响较大;烈度计(40通道)的功率谱值为-80— -60 dB,记录的噪声水平基本超过高噪声模型NHNM,表明烈度计记录的是设备自噪声,无法对台基噪声进行有效记录。
3.2 异常台站分析台站PDF图不仅可以反映其地噪声水平,还可以用来检验仪器是否正常,从而保证地震预警台网的数据质量。目前,在地震预警台网日常波形巡检工作中,发现部分台站加速度计(20通道)记录的连续波形存在背景异常、高频重复周期性干扰、台阶突刺、零点偏大等异常现象,从Jopens6.0系统的截图见图 6。以上异常现象出现的原因可能是运维人员对设备操作流程不熟未能正确安装调试,加之仪器无防护罩长期运行,缺乏保温和防气流措施;另外,基本站地处城镇地区,易受外界干扰。经简单排查处理,解决了部分台站波形数据异常问题,保证了山西地震预警台网的正常运行,相应处置措施见表 2。
部分台站波形异常经处理恢复正常,文中给出HX003、BD003、ZQXLY、JKOTH台站的PDF图,见图 7(对于多分向异常,受篇幅所限,仅列出单分向异常图)。结合表 2、图 6、图 7可知:①HX003台站:NS向波形背景异常,与此相对应的噪声PSD超出NHNM高噪声模型,更换摆线后恢复正常,在大于0.1 Hz频段范围,PSD曲线位于噪声模型内;②BD003台站:三分向波形出现高频重复周期性干扰,在噪声PDF图上表现为部分PSD曲线比较离散,更换加速度计后恢复正常;③ZQXLY台站:20通道NS、Z分向波形有连续突刺现象,PDF图中表现为PSD数值变大,部分区域超出高噪声模型,更换加速度计后恢复正常;④JKOTH台站:20通道NS、Z分向零点均大于100 Gal,PSD曲线超出高噪声模型,经现场调零后恢复正常,PDF图上明显可见调零前后PSD值的变化。由此可见,不同类型的异常在PDF图上的表现也不相同。
台站波形数据质量是地震监测和预警工作的基础,而仪器设备的运行状态与数据质量密切相关。本文以山西地震预警台网台站波形数据为基础,利用PPSD方法计算台基噪声水平,绘制台站PDF图及时频分布图。结果显示,基准站台基噪声水平较低,地震计能有效记录较宽频带范围噪声,PSD数值介于NLNM与NHNM之间;加速计仅能记录大于0.1 Hz的台基噪声;一般站不能有效记录台基噪声,PSD值超过高噪声模型NHNM曲线。
在地震预警台站中,产生波形异常的设备多为加速度计,可能与仪器安装调试不规范有关,且观测设备未加防护罩,无保温和防气流措施等。利用PPSD方法能检验仪器运行状态,以及时发现异常并恢复设备的正常运行,有助于台站波形数据质量监控。
包文超, 赵铁锁, 申影, 等. 内蒙古地震预警台站异常数据检测与分析[J]. 高原地震, 2023, 35(2): 50-55. |
陈建涛, 谢剑波, 吕仲杭, 等. 基于PPSD方法的广东阳江小孔径井下型地震监测台阵环境地噪声计算与分析[J]. 华北地震科学, 2019, 37(2): 21-29. |
丁大业, 董春丽, 宫卓宏, 等. 山西地震预警新建基准站数据质量评估[J]. 山西地震, 2023(2): 15-22. |
胡宝慧, 张浩, 教智博, 等. 地震台站噪声水平和台网监测能力自动化程序的实现[J]. 防灾减灾学报, 2022, 38(3): 42-47. |
刘旭宙, 秦满忠, 郭晓, 等. 用噪声概率密度函数对比地震计观测性能[C]//2018年中国地球科学联合学术年会. 北京: 中国和平音像电子出版社, 2018: 1 072-1075.
|
王俊, 徐戈, 孙业君. 江苏省区域地表背景噪声特性的分析[J]. 地震研究, 2009, 32(2): 155-161. |
徐嘉隽, 廖诗荣, 张红才, 等. 福建测震台网观测数据质量检测软件研究[J]. 华南地震, 2010, 30(4): 97-104. |
杨千里, 郝春月, 田鑫. 新疆和田台阵PSD与PDF分析[J]. 地球物理学报, 2019, 62(7): 2 591-2 606. |
Beyreuther M, Barsch R, Krischer L, et al. ObsPy: a python toolbox for seismology[J]. Seismological Research Letters, 2010, 81(3): 530-533. |
McNamara D E, Buland R P. Ambient Noise Levels in the Continental United States[J]. Bulletin of the Seismological Society of America, 2004, 94(4): 1 517-1 527. |
Sleeman R, Vila J. Towards an automated Quality Control Manager for the Virtual European Broadband Seismograph Network (VEBSN)[J]. Orfeus Electronic Newsletter, 2007, 7(1): 5-11. |