2. 中国成都 610041 四川省地震局
2. Sichuan Earthquake Agency, Chengdu 610041, China
对于地球背景噪声的产生机制,国内外目前主流观点是,其来源主要是海洋,其次是大气(特别是风与地表的相互作用),再次是地震和火山,最后是人类活动。目前,地震计不仅能够记录到各种类型的地震波在地下介质的传播,也能记录到外界条件(风、海浪、人类活动等)所产生的微弱震动(地脉动)。随着地震观测技术的不断发展,观测数据质量逐渐提升,微弱的地脉动也能被地震仪器所记录,台站背景噪声源的识别能力得以提高。因此,整体分析噪声来源和谱异常,可以有效评估台站监测能力,具有重要的实际意义。
Peterson(1993)收集了全球8个台网中75个地震台站的噪声记录,给出全球新高噪声模型(new high noise model,简写为NHNM)和全球新低噪声模型(new low noise model,简写为NLNM)。杨千里等(2019)利用Welch平均周期法,对和田地震台阵记录的噪声进行功率谱分析,认为噪声功率谱密度随季节变化显著,具有周期性,且受采石场影响,谱密度曲线在频段4—8 Hz出现高频尖刺。这为定量计算测震台站台基背景噪声提供了可行途径。何思源等(2019)利用成都地震台测震资料,得出该台2—4 Hz频段的背景噪声为固有干扰,且背景噪声水平于2011—2012年陡增并达到峰值,认为背景噪声水平变化主要由地质构造运动和台站周边城市扩建等引起,为分析台站形成的噪声异常提供参考。
若羌测震台(下文简称若羌台)属国家无人值守测震台,建台时观测场地环境地噪声水平属于Ⅰ级。近年来受人类活动影响,干扰持续增多,台基噪声水平发生了相应变化。为此,笔者收集整理2009—2019年若羌台地震波形数据,采用噪声功率谱密度(power spectral density,简称PSD),分析气压、人类活动、远震和近震等引起的谱异常,结合概率密度函数(probability density function,简称PDF),总结该台背景噪声水平的变化规律,以期对本地区背景噪声异常分析提供借鉴,为台站观测数据质量评估、本地区地震分析、台址勘选等提供参考。
1 台站概况及资料选取 1.1 台站概况若羌国家地震台位于新疆若羌县县城以东53 km处,海拔高度为1 390 m(图 1),场地气候条件属于大陆性荒漠干旱气候,全年降水量极少,年温差变化大,冬季最低温度可达零下25℃以下。台站位于米兰河西岸近EW向支沟内,地处阿尔金山造山带与塔里木地块交汇部位,构造上属青藏亚板块阿尔金强烈隆起区北缘。台站以北为塔里木凹陷沉积区,以南为元古代地层及岩浆岩组成的中高山区。台站附近发育江尕勒萨依断裂,属阿尔金断裂带次级断裂。该台站台基岩性为变质泥岩、砂岩,观测山洞进深22 m,2009年以来共架设2套宽频带地震计,具体参数见表 1。
选取若羌台2009—2020年共11年的波形数据,在每天24 h波形数据中截取无震、无明显干扰、长度1 h的记录,得到4 984个小时文件。对符合的波形数据进行分道处理,结合地震计灵敏度和数采转换因子,计算得到地面运动的速度记录,进行窗函数平滑处理及傅里叶变换,获得各段功率谱估计。将速度功率谱转换为加速度功率谱,利用地震计传递函数获得实际地面运动,在进一步转化为以dB为单位的台基加速度功率谱。然后将不同频段的加速度功率谱进行概率统计,获得台基的概率密度函数。
2 方法原理 2.1 台基噪声功率谱密度(PSD)计算采用Welch平均周期法(Welch,1967),进行噪声功率谱估计。设一离散信号x(n)的数据长度为N,分为L段(允许数据重叠),每一段为M,取x(n)进行傅里叶变换,得到xN(e-j2πfn),求取每段数据功率谱进行平均。每段数据均选用汉宁窗,以改善矩形窗边瓣较大产生的谱失真。其中x(n)的功率谱密度PPERi(f)可以表示为
$ P_{{\rm{PER}}}^i\left(f \right) = \frac{1}{{MU\Delta f}}{\left| {\sum\limits_{n = 0}^{M - 1} {x_N^i\left(n \right){\rm{d}}\left(n \right){{\rm{e}}^{ - j2{\rm{ \mathsf{ π} }}fn}}} } \right|^2} $ | (1) |
式中,PPERi(f)表示为第i段的加速度功率谱密度,U为归一化因子,d(n)为汉宁窗口,f为尼奎斯特频率,n为有限随机信号的个数。几段数据平均功率谱PPER(f)可表示为
$ {P_{{\rm{PER}}}}\left(f \right) = \frac{1}{{MUL\Delta f}}\sum\limits_{i = 1}^L {{{\left| {\sum\limits_{n = 0}^{M - 1} {x_N^i\left(n \right){\rm{d}}\left(n \right){{\rm{e}}^{ - j2{\rm{ \mathsf{ π} }}fn}}} } \right|}^2}} $ | (2) |
目前,功率谱概率密度函数(PDF)是对台站噪声水平整体评估的有效途径之一,无需进行特殊事件(地震、仪器标定等)筛选,是对数据全时段的分析,能够全面反映仪器的记录状态和能力。为了提高运算速率,将计算所得PSD结果进行1/8倍频程间隔采样,将各中心频率功率值求平均。功率谱密度值在上端频率fu与下端频率fl之间进行平均,获得中心频率fc,可表示为
$ {f_{\rm{c}}} = \sqrt {2{f_{\rm{u}}}{f_1}} $ | (3) |
fu以1/8倍频程的功率窗口fu= fu×21/8增加,计算下一间隔平均功率,对于每一个中心频率fc的概率密度函数,可以表示为
$ P\left({{f_{\rm{c}}}} \right){\rm{ = }}\frac{{{N_{P{f_{\rm{c}}}}}}}{{{N_{{f_{\rm{c}}}}}}} $ | (4) |
式中,NPfc表示频率为fc时,功率落在某个1 dB窗口中的个数,Nfc表示每一中心频率功率值的总和。
3 噪声特征分析通过浏览若羌台大量PSD曲线和PDF图,分析认为,谱异常主要受以下因素影响:①自然条件:气压变化;②人类活动:米兰河水库施工建设。结合显著地震事件波形数据的PSD曲线,总结若羌台背景噪声功率谱特征。
3.1 自然条件引起的谱异常若羌县属大陆性荒漠干旱气候,昼夜温差较大,测震观测长期受气压干扰影响。为了解若羌台不同季节、不同频段的功率谱变化情况,选取2009年1月11日17时(冬季)和2009年7月7日5时(夏季)记录的测震数据,计算PSD,获得功率谱密度曲线,结果见图 2,其中(a)图表示冬季时段的PSD和波形数据,(b)图表示夏季时段的PSD和波形数据。
由冬、夏季节2个时段PSD曲线对比可知,气压变化对EW向和NS向数据影响明显,对UD向影响不明显。这是因为:①在气压的微弱作用下,观测洞室被覆山体发生变形;②观测洞室位于山体北坡,受夏季昼夜温差大的影响,山体南坡形变大于北坡,垂直向耦合到水平向,使得水平向噪声得到放大。由图 2(b)可见,在夏季,在小于0.1 Hz长周期频段内,噪声水平接近高噪声模型(NHNM),在0.008—0.015 Hz频段,噪声水平甚至高于NHNM。由图 2中原始波形数据曲线可知,夏季水平向(EW向、NS向)噪声水平明显高于冬季。
将2009年无断记数据PSD曲线进行PDF统计,获得三分向功率谱概率密度图(左侧纵坐标代表功率谱密度噪声值,右侧纵坐标代表某一单位为1 dB的频率出现的概率,横坐标为具体频率),结果见图 3。由图 3可见,水平向(NS、EW向)长周期频段噪声水平在-150— -125 dB,垂直向(UD向)噪声水平相对收敛,集中在8 dB左右;在高频拐点后,三分向出现高频异常干扰,应为仪器自身影响所致。笔者认为,气压、气温变化导致水平向噪声具有相同的变化幅度。
选取若羌台记录的4个远震事件(地震参数见表 2)进行功率谱计算,PSD曲线见图 4,(a)、(b)、(c)、(d)图所示地震对应表 2中a、b、c、d事件。4个远震事件波形功率谱变化均出现在长周期频段,噪声水平高于NHNM,拐点集中出现在-90— -80 dB,但PSD曲线变化趋势各异。
事件a、b长周期频段集中在0.04—0.06 Hz,应为长周期面波频段。二者震中距相近,震源深度不同,分属浅源地震和深源地震。在微震频段0.2—2 Hz,事件a的背景噪声水平受地震的影响低于事件b,噪声水平更接近NLNM,这是因为,深源地震事件b体波发育,能量衰减缓慢,地震计可较为详细地记录到地震波形。事件a、d均为浅源地震,但随着震中距的增加,事件d在频率0.01 Hz附近出现拐点,这是因为,该事件属极远震,面波发育,周期大,持续时间长,在0.2—2 Hz频段噪声变化大,而事件a则变化平稳。在0.1—1 Hz频段,事件d受震中距影响,体波衰减快,噪声水平明显低于事件a,接近NLNM。事件b、c震源深度相近,震中距不同,由PSD曲线可见,震中距增大,长周期拐点出现“左移”现象。在频段0.1—1 Hz,事件c的噪声水平明显高于事件d,这是因为,前者属于深源地震,体波衰减慢,台站记录到该频段波形数据噪声水平较大。事件c、d均属于远震,在频段0.1—1 Hz内不产生面波。结合事件c、d的背景噪声水平,可知在微震频段0.1—1 Hz,相较震中距而言,震源深度对背景噪声水平变化的影响较大。
3.3 近震记录功率谱特征选取4个中小近震事件(地震参数见表 3)进行功率谱计算,结果见图 5,(a)、(b)、(c)、(d)图对应表 3中a、b、c、d事件。由波形图可见,近震事件a、b的面波基本不发育,长周期频段噪声低,在频段0.3—6 Hz,噪声水平均高于NHNM,形态一致,表明此为近震引起噪声水平变化的主要频段。由事件a、c的波形可见,事件a较短周期体波发育,事件c较长周期面波发育,在PSD图上,起伏波动不一致。从整体看,不同近震记录的功率谱特征相差不大,主要是因为,震中距较近,不同周期的震相发育不完全。
在2016年底,距若羌台约1 km的米兰河水库建设完成并投入运营。选取该台2009年和2018年2个平静时段测震数据,利用快速傅里叶变换FFT,计算2个时段UD向振幅谱,见图 6。由图 6明显可见:①受持续增加的人类活动影响,在频段3 Hz附近,与2009年相比,2018年背景噪声水平明显增高;②2018年的频谱数据在10 Hz频段附近出现高频异常,属新增噪声干扰。笔者对若羌台BBVS-60地震计2015年7月—2016年3月架设期间的记录数据进行频谱分析,发现无10 Hz频段,表明此为后期某段时间形成的固有频段,与仪器设备无关。
笔者将2009—2019年无断记波形数据进行功率谱概率密度计算,获得UD向PDF图,见图 7。由每年的PDF图可知,在1—20 Hz高频段,2011—2016年的噪声水平明显高于其他年份,2017年之后趋于稳定。据调查,2011—2016年为米兰河水库施工时段,人类活动影响较大。在2012年PDF图上,在高频段可见曲线高于NHNM噪声模型并快速下降到高频部分,结合同年噪声PSD曲线,发现该异常为水库建设时施工爆破所致。由2017年以后PDF图可见,在频率10 Hz附近出现“方形”拐点,而施工期间的大量PSD曲线中并不存在该频段。据悉,2017年初米兰河水库投入运营并网发电,对若羌台测震观测造成固定干扰,由此形成此处谱异常,采用滤波方法剔除相应频段数据,即可获得高质量波形资料。
笔者收集整理若羌测震台2009—2019年无断记波形数据,计算并绘制PSD和PDF图,分析该台近年噪声水平,得到以下认识:①若羌台噪声水平受气压影响显著,其中水平向受影响显著,应为噪声水平从垂直向耦合到水平向的结果;②若羌台受持续增加的人类活动影响,背景噪声水平不断升高;③结合PDF图,10 Hz频段干扰是米兰河水库投入运营形成的固定干扰,与仪器设备无关联性。
若羌台地处阿尔金山和塔里木盆地交汇处,地广人稀,台基环境较好,是测震监测台网的重要组成部分。近年来受米兰河水库建设的影响,台站环境地噪声水平持续升高,形成固定干扰,使得该台环境噪声水平持续下降,应合理改善观测条件,保护地震台基监测环境。
文中计算PDF采用的Matlab程序由四川省地震局谢江涛工程师提供,图 1绘制使用了GMT软件,在此一并表示感谢。
高伟亮, 吕睿, 梁艳, 等. 山西数字测震台站观测动态范围和台基背景噪声分析[J]. 山西地震, 2015(4): 17-24, 42. DOI:10.3969/j.issn.1000-6265.2015.04.005 |
葛洪魁, 陈海潮, 欧阳飚, 等. 流动地震观测背景噪声的台基响应[J]. 地球物理学报, 2013(3): 857-868. |
何思源, 李贵元. 关于成都地震台测震资料中2-4 Hz干扰频率范围的研究[J]. 内陆地震, 2019, 33(1): 80-90. |
刘旭宙, 沈旭章, 李秋生, 等. 青藏高原东北缘宽频带地震台阵远震记录波形及背景噪声分析[J]. 地球学报, 2014, 35(6): 759-768. |
刘洋君, 杨毅, 王燕, 等. 浙江测震台网地震监测能力分析[J]. 地震地磁观测与研究, 2016, 37(4): 57-61. DOI:10.3969/j.issn.1003-3246.2016.04.011 |
王芳, 王伟涛, 龙剑锋, 等. 中国大陆地区宽频带地震台网台基噪声特征[J]. 地震学报, 2019(5): 569-584. |
谢江涛, 林丽萍, 谌亮, 等. 地震台站台基噪声功率谱概率密度函数Matlab实现[J]. 地震地磁观测与研究, 2018, 39(2): 84-89. DOI:10.3969/j.issn.1003-3246.2018.02.012 |
许可, 刘瑞瑞, 孔繁旭. 天津地区台基背景噪声特征分析[J]. 内陆地震, 2015(2): 170-175. |
杨千里, 郝春月, 田鑫. 新疆和田台阵PSD与PDF分析[J]. 地球物理学报, 2019(7): 2591-2606. |
张荣杉, 冯武, 彭澎, 等. 宿迁地震台地震监测能力分析[J]. 地震地磁观测与研究, 2017, 38(4): 199-202. DOI:10.3969/j.issn.1003-3246.2017.04.033 |
赵铁锁, 张晖, 高昌志, 等. 钻井对呼和浩特地震台地震观测影响[J]. 地震地磁观测与研究, 2016, 37(6): 90-94. DOI:10.3969/j.issn.1003-3246.2016.06.016 |
中国地震局. 测震学原理与方法[M]. 北京: 地震出版社, 2015.
|
Welch P. The Use of Fast Fourier Transform for the Estimation of Power Spectra: A Method Based on Time Averaging over Short, Modified Periodograms[J]. IEEE Transactions on Audio and Electroacoustics, 1967, 15(2): 70-73. DOI:10.1109/TAU.1967.1161901 |