0 引言

地倾斜、地应变等洞体形变观测能直接反映地壳介质的微动态变化，但在观测过程中不可避免地受到各种因素，如观测环境、雷电、降雨、电源电压变化、外界干扰、电子设备自身因素等的影响，观测资料中会出现各种“异常”数据，如数据曲线的毛刺、大幅突跳、台阶跃变、趋势转折等。在应用基础观测数据开展研究前，需对异常数据进行处理，因此识别各类干扰和建立各观测仪器的正常观测背景尤为重要。传统数据分析方法适用于有辅助测项和年变周期形态较好的资料，且仅可在形态上进行数据的定性分析，在日常应用中则受到诸多限制。因此，引进新的数字信号处理技术，可以有效识别干扰信息，提高对数据的认识。

功率谱密度（power spectrum density，PSD）分析是目前用于信号处理较为成熟的手段，它对信号或时间序列功率随频率的分布进行了定义，广泛用于重力仪、形变观测仪和各类地震计的干扰分析。使用功率谱密度分析方法开展的相关研究有：韩宇飞等（2015）、徐义贤（2015）研究了陆态网络太原台gPhone重力仪背景噪声水平；赵莹(2018, 2019)分析了垂直摆观测背景噪声水平；侯颉等（2019）分析了北京市测震台网数字地震台站背景噪声；胡宝慧等（2020）分析了黑龙江地震科学台阵台基背景噪声特征；李萌等（2020）分析了GPS高频信号的中强度地震频谱特性，等等。以上文献缺少利用功率谱密度对同一台站多套洞体形变仪器数据的分析，若利用该方法对同一台站洞体形变观测进行分析，将有利于干扰源频谱特性的识别。

山西省北武当地震观测站（下文简称北武当站）洞体形变观测手段齐全，配备有水平摆倾斜仪、水管倾斜仪和洞体应变仪，是该省中部地区的一个重要的地形变观测站点，然而各种干扰信号叠加在形变观测数据中，增加了日常数据的预处理难度，不利于地球物理异常的识别（王晓霞等，2019）。本文利用功率谱密度分析方法（下文简称PSD方法），对北武当站洞体形变观测数据进行处理，分析观测资料在不同频段的正常背景特征及不同干扰因素下的PSD分布规律，以进一步认识观测数据的正常背景和干扰特征，为后续异常干扰库及定量化预测指标的建立服务，并对日常震情跟踪工作起到一定推动作用。

1 数据处理

干扰因素包括自然环境（大风、降雨、雷电、气压、温度等）干扰、人为干扰、观测系统故障及震扰等。

通过查看形变观测数据曲线和观测日志，筛选未受到以上干扰的平稳时段数据和存在干扰影响时段的数据，对于存在缺数的原始时间序列，采用线性插值方法进行拟合，得到完整、连续的时间序列数据，总结各类干扰的变化形态，并按形态进行归类，利用频谱分析方法，滤除主要潮汐成分，并扣除线性趋势（在进行趋势变化分析时，仅去除潮汐成分），计算每日PSD并作图。为减少计算过程中的“频谱泄露”效应，将每日观测数据按小时划分，将每小时数据分为4个小数据段，且各段之间两两重合50%，采用快速傅里叶变换（FFT）求得每段PSD，并基于各段PSD平均值获得每日PSD。

设信号为x(j )，j = 1, 2, ..., n，则信号的PSD为

$ \mathrm{PSD}=\frac{\Delta t}{N}\left|\sum\limits_{j=1}^n x(j) \exp (2 \pi i(j-1)(k-1) / N)\right|^2 $

(1)

式中，Δt为采样间隔，N为采样点总数，i为虚数单位。

在一定频率范围内，平均功率谱密度meanPSD为

$ \text { meanPSD }=\frac{1}{M} \sum\limits_{k=k_1}^{k=k_2} \operatorname{PSD}(k) $

(2)

式中，k₁和k₂分别为所选取频率范围的上限和下限，M为所选取频率范围内的采样点数。

2 洞体形变数据功率谱密度响应特征

选取北武当站水平摆倾斜仪、水管倾斜仪、洞体应变仪观测数据，由于NS、EW分量计算结果一致，故文中以NS分量为例，分析该观测站洞体形变数据功率谱密度响应特征。

2.1 日变特征分析

选取北武当站2021年6月1日—5日NS分量分钟值观测数据，计算每日PSD并作图。如图 1所示，在未受干扰情况下，3套仪器NS分量每日PSD值变化不大，且在全频段变化形态基本一致，可反映最大频率为8.3×10^-3 Hz的功率谱密度分布，该频段背景噪声源主要包括局部大气作用（＜2×10^-3 Hz）、hum信号（2×10^-3—7×10^-3 Hz）及Rayleigh波（7×10^-3—30×10^-3 Hz）（韩宇飞等，2015）。

2.2 各类干扰响应特征分析

洞体形变观测不可避免地受到各类干扰因素的影响，在观测曲线上表现为形态或趋势变化，其中：①自然环境干扰：风扰导致观测曲线抖动加粗，雷电导致观测曲线阶变和突跳；②人为干扰：观测曲线突跳、畸变；③观测系统故障：观测曲线畸变；④震扰：观测曲线上出现明显的地震波。基于观测曲线的变化形态，总结3套仪器常见典型干扰，并分析其噪声功率谱密度特征。

2.2.1 曲线加粗

北武当站水平摆倾斜仪、水管倾斜仪、洞体应变仪会受到不同程度的大风干扰，且水平摆倾斜仪所受风扰更明显。2018年4月4日—5日该观测站出现大风天气，选取水平摆倾斜仪当年4月1日—5日NS分量分钟值数据进行分析，其中绿色线代表 4月4日大风干扰日，红色线代表 4月5日大风干扰日，结果见图 2。

由图 2(a)可见，受大风干扰，北武当站当月4日—5日数据曲线加粗，毛刺现象明显，其中4日扰动幅度较小且时间较短，5日扰动明显较大，且持续时间较长。由图 2(b)可见：①水平摆倾斜仪对大风干扰较明显，风扰日观测曲线响应特征明显，扰动越大，PSD值变化越大；②干扰日能量强度明显高于平静日，且自频率10^-3 Hz起功率谱密度逐渐增大，最大幅度达18 dB/Hz；③风扰表现为高频干扰，低频部分PSD值分布和平静日较接近。

2.2.2 台阶、突跳和畸变

（1）雷电干扰导致台阶和突跳。雷电易导致洞体应变观测数据曲线出现台阶和突跳。2018年6月9日北武当站受雷电干扰，选取该观测站洞体应变仪6月7日—10日NS分量分钟值数据，计算PSD值并作图，结果见图 3。由图可见，受雷电干扰，9日观测曲线出现较大幅度的台阶和突跳，对应的PSD值较平静日全频段增大，且低频段增大幅度高于高频段，干扰日PSD最大值达61 dB/Hz。

（2）调零干扰导致数据突跳。水管倾斜仪受调零干扰较多，且调零过程中数据出现大幅突跳。2020年6月8日，北武当站水管倾斜仪进行调零，选取6月5日—8日NS分量分钟值数据进行分析，结果见图 4。由图可见：受调零干扰，6月8日观测数据曲线出现大幅突跳[图 4(a)]，对应的PSD值较平静日显著增大，且全频段增大幅度均较大，其中在(0.2—0.5)×10^-3 Hz频段内，PSD值增幅最大，变化幅度约达到40 dB/Hz[图 4(b)]。

（3）观测系统故障、场地干扰较多导致数据畸变。选取水平摆倾斜仪2020年2月23日—27日NS分量分钟值数据进行分析，其中红色线代表 2月23日畸变日，绿色线代表 2月24日畸变日，结果见图 5，可见23日00:00—23:59、24日00:00—05:32观测数据畸变[图 5(a)]，对应PSD值全频段呈增大现象，且与观测数据畸变时长及幅度呈一定相关性，表现为畸变时间越长，幅度越大，PSD值增幅越大[图 5(b)]。

综上可知：观测数据原始曲线台阶、突跳和畸变对应的PSD变化特征较相似，均表现为全频段增大，且其增大幅度与干扰的持续时间及幅度成正比。

2.2.3 趋势改变

趋势改变是固体潮曲线发生趋势性上升、下降或者转折的现象。选取水平摆倾斜仪2016年1月26日—29日NS分量分钟值数据进行分析，结果见图 6，可见26日观测曲线加速南倾[图 6(a)]，对应PSD值增大幅度较小[图 6(b)]，表明趋势转折对功率谱密度值的变化影响不大。

2.2.4 地震波

相对于正常背景而言，地震波可以看作是一种对观测数据产生影响的干扰信号。一般情况下，洞体形变观测可记录到清晰地震波形，其固体潮曲线形态表现为震荡。2013年7月22日07:45:55甘肃省定西市岷县、漳县交界发生M_S 6.6地震（34.52°N，104.23°E），震源深度20 km，距震中734 km的北武当站洞体形变同震响应明显。选取水平摆倾斜仪和水管倾斜仪2013年7月22日—26日NS分量分钟值数据进行PSD计算，并对2套仪器记录的同时段地震波形数据的频谱响应特征进行对比分析，结果见图 7、图 8。

由图 7(a)、图 8(a)对比可见，水平摆倾斜仪同震效应较水管倾斜仪明显，其同震响应曲线变化幅度较大、持续时间较长。

由图 7(b)、图 8(b)对比可见：①7月22日地震发生当天，洞体形变观测数据PSD值较平静日明显增大；②水平摆倾斜仪PSD值全频段增大，但高频部分增大幅度较大，水管仪仅高频部分出现增大现象，且所在频段范围较小，这是因为，对于同一地震，水平摆倾斜仪记录的地震波最大振幅更大、持续时间更长；③水平摆倾斜仪和水管倾斜仪PSD值变化趋势基本一致，尤其在高频部分（频段200—600 s），地震波形呈“V”字型变化，该现象与地震波自身特性相关。总之，不难发现，地震波的变化幅度越大、持续时间越长，PSD值的频段范围和变化幅度就会越大。

3 结果与讨论

利用功率谱密度分析方法（PSD），对北武当站洞体形变观测的正常背景及存在干扰的PSD值特征进行分析，结果表明：

在正常背景下，形变观测数据PSD值变化不大，且3套数据变化形态基本一致，基于分钟值数据的分析结果，可反映最大频率为8.3×10^-3 Hz的PSD分布。

风扰日PSD响应明显高于平静日，且扰动越大，PSD变化越大；风扰日能量强度明显高于平静日，且频率大于10^-3 Hz时逐渐增大，最大幅度达18 dB/Hz；风扰表现为高频干扰，低频段PSD值变化不大，与平静日相当。

人为干扰、观测系统故障、雷电等引起的突跳、台阶、畸变，对应PSD值变化较为相似，均表现为PSD值全频段增大，干扰日数值增大显著，与观测数据干扰的幅度和时长相关，干扰数据变化幅度越大，干扰时间越长，PSD值增幅就越大。

观测数据曲线的趋势转折变化对PSD值影响不大，在计算频段内，干扰日和平静日的PSD值较为接近，可能与观测数据曲线形态基本一致有关。

地震发生当日，噪声信号PSD值较平静日明显增大，且高频段增大幅度显著；水平摆倾斜仪和水管倾斜仪观测数据PSD值变化趋势基本一致，尤其在高频段，均呈“V”字型变化；随着同震效应增大，地震波形数据PSD值的频段范围和变化幅度均会增大。