环境噪声是影响测震台站观测数据质量的主要因素之一。观测数据高频段主要受到工业活动、交通运输等人为因素的影响^[1]，微震频段主要受到海洋活动的影响^[2]，低频段主要受到温度、气流变化等因素的影响^[1]。21世纪以来，井下观测技术在地震监测中得到广泛应用，能有效降低地震观测数据低频和高频背景噪声水平^[3-5]。内蒙古地区春季沙尘暴较为严重，可能会对地震观测环境噪声水平产生影响，但目前尚未有此方面的研究。因此，本文计算二连浩特台地表和井下连续观测数据的加速度PSD值和相应的PDF值，分析沙尘暴对地震环境噪声水平的干扰特征和井下观测方式对减小沙尘暴影响的效果。

1 台站信息和数据来源

二连浩特台位于内蒙古锡林郭勒盟北部，全年干旱且风沙较大，沙尘暴主要发生在3~5月。台站同时具备地表和井下2种观测方式，配有甚宽频带地震计和24位数据采集器，仪器信息见表 1。2021年春季二连浩特地区沙尘暴天气较往年更加突出，根据二连浩特市气象局和环保局提供的信息，2021-03-14~15、2021-04-26沙尘暴较为严重(表 2)，本文以这2个时间段数据为例进行分析。

2 数据处理

首先以1 h为间隔将数据分成若干段，采样率为100 Hz，再将每小时的数据段分为42个记录段。为降低PSD方差，相邻记录段间重叠率为50%，则每个记录段长度约为160 s。为降低长周期对功率谱估计的偏差，对每个记录段进行去长周期和去均值的预处理；然后将预处理后的记录段数据进行快速傅里叶变换，得到以频率为自变量的速度PSD值，再将速度PSD值转换为加速度PSD值；最后进行平滑处理并计算与PSD值相应的PDF值。

2.1 PSD的计算

对周期时间序列y(t)的有限范围进行傅里叶变换：

$ Y\left( {f,{T_r}} \right) = \int_0^{{T_r}} y (t) \cdot {{\rm{e}}^{ - {\rm{i}}2\pi ft}}{\rm{d}}t $

(1)

式中，T_r为时间序列段长度，f为频率。

对离散频率值f_k进行傅里叶变换：

$ {Y_k} = \frac{{Y\left( {{f_k},{T_r}} \right)}}{{\Delta t}} $

(2)

式中，f_k=k/(NΔt)，k=1，2，…，N-1，Δt为采样间隔(0.01 s)，N=T_r/Δt为截取时间段的采样点数。

PSD的定义为：

$ {P_k}(f) = \frac{{2\Delta t}}{N}{\left| {{Y_k}} \right|^2} $

(3)

将速度PSD值转换为加速度PSD值：

$ {P_{a,k}}(f) = {(2{\rm{ \mathsf{ π} }}f)^2}{P_k}(f) $

(4)

消除仪器传递函数的影响：

$ {{\mathop{\rm PSD}\nolimits} _a}(f) = \frac{{{P_{a,k}}(f)}}{{|H(f){|^2}}} $

(5)

式中，PSD_a(f)为真实的地表运动加速度功率谱。

2.2 平滑处理

为了使PSD在频域对数坐标中等间隔采样，采用1/3倍频积分进行平滑处理：

$ {{\mathop{\rm PSD}\nolimits} _a}\left( {{f_c}} \right) = \frac{1}{n}\sum\limits_{{f_c}}^{{f_h}} {{{{\mathop{\rm PSD}\nolimits} }_a}} (f) $

(6)

式中，f_ι=2^－1/6f_c为低频拐角频率，f_h=2^1/6f_c为高频拐角频率，n为介于二者之间的频率f的个数。由式(6)可得中心频率f_c的PSD_a(f)平均值PSD_a(f_c)，将其作为f_c的PSD值。中心频率f_c以1/9倍频程为增加步长，即下一个中心频率为2^1/9f_c，重新计算相应的f_ι和f_h，然后将新的f_ι和f_h之间的PSD平均值作为下一个中心频率f_c的PSD取值。这样在f_c的取值范围0.01~50 Hz内，每个记录段的PSD值随频率的变化情况可由在对数坐标系中呈等间隔采样的中心频率的PSD值来表示。

2.3 PDF的计算

每个中心频率f_c的PSD的PDF为：

$ P_{\mathrm{PSD}}\left(f_{c}\right)=N_{P f_{c}} / N_{f_{c}} $

(7)

式中，N_fc为f_c频点记录段总数，N_Pfc为f_c频点PSD值落在某PSD取值范围内的记录段个数，本文PSD窗长与步长都取1 dB，变化范围为-200~-50 dB。一般以频率为横坐标、PSD为纵坐标、P_PSD(f_c)色块颜色深浅为元素绘制三维平面图，得到PDF分布，不同色块代表某频点在一定PSD窗内功率谱的概率数。

3 数据分析 3.1 1~20 Hz频段PSD值分析

为研究沙尘暴对1~20 Hz频段环境噪声的影响特征，选取二连浩特台无沙尘暴天气(2021-03-01~02)和沙尘暴天气(2021-03-14~15、2021-04-26)的连续观测数据，分别计算观测数据三分向PSD值，结果见图 1。

对比图 1(a)~1(c)和图 1(d)~1(i)可见，有大风、沙尘暴时间段的背景噪声水平明显高于无沙尘暴时间段的背景噪声水平，这种差异主要体现在4 Hz和9 Hz附近及10~20 Hz频段。由图 1(d)~1(f)可见，2021-03-14 19:00~2021-03-15 01:00 10~20 Hz频段PSD值有一个突然增加的过程。参考表 1可知，这个时间段与气象局统计的强沙尘暴时间段吻合。2021-03-14 13：00~2021-03-15 01：00 4 Hz、9 Hz附近的PSD值一直保持较高状态，这个时间段基本与气象局统计的大风时间段吻合。由图 1(g)~1(i)可见，2021-04-26 06:00~19:00 4 Hz、9 Hz附近及10~20 Hz频段的PSD值一直保持较高状态。参考表 1可知，这个时间段与气象局和环境监测站统计的扬沙时间段重合。综上所述，沙尘暴天气对于地震观测数据背景噪声高频的干扰主要体现在4 Hz、9 Hz附近及10~20 Hz频段，其中沙尘主要影响10~20 Hz频段的背景噪声，大风主要影响4 Hz、9 Hz附近的背景噪声。

3.2 沙尘暴干扰特征分析与压制

分别选取沙尘暴最突出(2021-03-14 13:00~2021-03-15 01:00、2021-04-26 08:00~20:00)和无沙尘暴(2021-04-04 00:00~2021-04-05 00:00)各24 h的地表和井下连续波形数据，分别计算最高PDF值对应的PSD平均值(图 2)，研究沙尘暴对环境噪声水平的干扰特征。

3.2.1 有无沙尘暴天气的地表PSD值对比

由图 2(a)、2(b)可见，在4~10 Hz频段，无沙尘暴天气观测数据三分向PSD值变化较为平缓，有沙尘暴天气观测数据PSD值高于无沙尘暴天气，尤其是在4~5 Hz和7~10 Hz频段内，沙尘暴天气观测数据三分向PSD值均有突变升高，且变化趋势较为一致。说明沙尘暴在4~10 Hz频段对环境噪声水平有较大的影响，其中，4~5 Hz和7~10 Hz频段对环境噪声水平的影响更明显，在大于10 Hz频段，有沙尘暴天气观测数据三分向PSD值相比无沙尘暴天气有一定的增加，其中垂直向最显著，东西向次之，南北向变化不大。说明在大于10 Hz频段内，沙尘暴对观测数据的环境噪声有一定影响，且对三分向影响程度不同。在10 s~4 Hz频段，有无沙尘暴天气的观测数据三分向PSD值基本一致，9.5 s和2.5 Hz处沙尘暴天气观测数据的PSD峰值甚至低于无沙尘暴天气的PSD峰值。说明在10 s~4 Hz频段内沙尘暴对观测数据没有影响。在低于10 s频段，沙尘暴天气观测数据东西向PSD值高于无沙尘暴天气15 dB以上，垂直向和南北向PSD值变化不大。参考表 1可知，本次沙尘暴风向以偏西为主，所以导致低于10 s频段的东西向背景噪声受干扰最大。

3.2.2 井下观测对减小沙尘暴影响的效果

由图 2(c)、2(d)可见，整个频段内有无沙尘暴天气的观测数据PSD值基本一致，只有在低于60 s频段内沙尘暴天气观测数据的PSD值才有一定的升高，说明井下观测可明显压制沙尘暴对观测数据环境噪声的干扰。由图 2(a)、2(c)可见，在10 s~1 Hz频段，2种观测方式的PSD值一致；在低于10 s频段，井下观测数据三分向PSD值的一致性更高；在高于1 Hz频段，井下观测数据三分向PSD值低于地表观测数据PSD值10~20 dB左右。再次说明井下观测能够有效压制沙尘暴对环境噪声的干扰，且大于1 Hz的频段压制效果最好。

4 结语

1) 沙尘暴对测震台站环境噪声水平有较大干扰，主要体现在4 Hz、9 Hz附近及10~20 Hz和低于10 s频段。其中，沙尘主要影响10~20 Hz频段的背景噪声；大风主要影响4 Hz、9 Hz附近的背景噪声；低于10 s频段内三分向背景噪声主要受风向的影响。10 s~4 Hz频段内沙尘暴对环境噪声水平没有影响。

2) 井下观测能够有效压制沙尘暴对环境噪声的干扰。