引言

地形变测量中的定点潮汐形变观测主要以地倾斜、地应变等测量指标为主，其目的在于记录地壳周期性潮汐形变，并能捕捉地震中短临前兆异常信息（孙伶俐等，2013）。但实际形变观测结果受多方面因素影响，主要包括仪器系统状态、固体潮汐、观测环境、人为因素、自然环境、地震波、地下介质变化等（蔡佩蕊等，2020）。为分析和判识观测中的各类干扰信号、潮汐波及地震前兆异常信息等问题，近年来，诸多数字信号处理技术被引入形变资料分析领域，如小波分析方法（侯跃伟等，2015）、模态分解类方法（孔祥瑞等，2018）、超限率分析方法（李娜等，2020）、奇异谱分析方法（赵佳佳等，2017）、时频分析方法（王宁等，2014）等。其中，时频分析方法以其能够表征信号频谱随时间变化规律的优势，在形变观测干扰分析（张小艳等，2019）、前兆异常判定（武善艺等，2018）、日常跟踪分析（孟庆筱等，2018）等方面得到了应用与发展。目前，制约时频分析方法进一步在台站形变数据处理中的应用主要问题为：①传统的连续小波变换（CWT）、短时傅里叶变换（STFT）、S变换（ST）、希尔伯特-黄变换（HHT）等方法时频分辨率较低，对于波形类别多、分布频段广、振幅差异大的各类形变观测信号的泛化适用能力较差；②台站缺乏相关集成处理软件程序（平建军等，2013），数据导入、函数传递、方法计算、成果绘图等功能和步骤缺乏统一设计规划，不便于操作与推广。

为解决上述问题，本研究基于MATLAB科学计算平台研发了定点形变观测数据时频分析软件，该软件通过集成数据导入、插值补全、小波分析、同步挤压时频变换、自动绘图等功能，处理过程采用命令窗参数绘图交互方式，能够灵活计算特定频段、特定时段形变观测数据的高精度时频谱，以期为台站工作人员深入挖掘形变观测资料中的各类潜在变化特征提供简便实用的处理工具。

1.   方法原理

本研究采用信号分解方法中的小波分析（DWT），时频处理方法中基于CWT和广义S变换（GST）的同步挤压小波变换（SSCWT）和同步挤压广义S变换（SSGST），以及时频辅助分析方法中的叠加幅值特征函数（SCCF）、叠加系数包络函数（SCEF）和信息熵（Renyi），全面刻画不同频段分布的复杂形变观测数据时频响应特征。

DWT是利用选定小波基和分解层数以二进塔式分解方式将原始信号分解到一系列相互独立的频带上（Mallat，1989），凸显信号中不同频段成分信息，并可通过筛选重构得到特定优势频段子信号分量的方法。

CWT是将原始信号与经过伸缩平移后的母小波作褶积运算的方法，具有多尺度表达信号频率分辨率特征，其时频谱结果频率轴尺度为指数型分布，对于中低频段信号具有较好的时频响应能力。

GST是对传统ST方法高斯窗函数引入2个可调参数（魏学强等，2016），使其窗函数随频率呈现出更多非线性变化特征而提高时频分辨率的方法，其时频谱结果频率轴尺度为线性分布，对于高频段信号具有较好的时频响应能力。

SSCWT、SSGST是分别在CWT和GST“前处理”结果的基础上，进一步利用基于微分算子的挤压“后处理”操作，沿频率轴重排时频谱瞬时频带，将特定频率区间内的时频系数叠加到理论频率点上（严海滔等，2019；潘晓等，2020），可得到时频聚集性更精确的同步挤压时频结果的方法。

SCCF、SCEF是分别描述刻画同步挤压时频谱在频率尺度上的能量分布特征和在时间尺度上的强度分布特征（Mousavi等，2016）的指标。

Renyi熵是定量检测时频谱所含信息量和时频谱能量集中度的指标（黄昱丞等，2018）。

整个方法计算流程为：（1）利用DWT方法对原始形变信号作多尺度分解；（2）根据设定的研究目标，对分解出的子信号及时间范围作选择性重构，以提取待分析的目标信号分量；（3）根据目标分量的频段分布特征，基于分频处理策略，对于短周期高频目标信号和长周期低频目标信号分别采用高频响应较好的GST方法和中低频响应较好的CWT方法进行时频计算，在此基础上利用同步挤压操作得到更高分辨率的SSGST和SSCWT时频结果；（4）进一步对计算出的同步挤压时频谱结果，采用SCCF方法、SCEF方法和Renyi熵分别从频率尺度、时间尺度、能量复杂度方面进行辅助分析。

上述计算方法基本公式和控制参数如表1所示。

表 1  软件包所用方法及基本定义

Table 1.  Summary of methods and basic definitions of the software package

类型	方法名称	基本定义	备注
信号分解方法	小波分析 （DWT）	${A}_{j}s\left(t\right)={A}_{j+1}s\left(t\right)+{D}_{j+1}s\left(t\right),j=\mathrm{0,1},\cdots ,N$	$ s\left(t\right) $为原始信号，$ N $为分解层数，$ {A}_{j}s\left(t\right) $为第$ j $阶低频趋势分量，$ {D}_{j}s\left(t\right) $为第$ j $阶高频细节分量
时频前处理 方法	连续小波变换 （CWT）	$ C W T\left(a,\tau \right)=\left\langle{s,{\psi }_{a,\tau }}\right\rangle={\displaystyle\int }_{-\infty }^{+\infty }s\left(t\right){a}^{-\frac{1}{2}}{\psi }^{*}\left(\dfrac{t-\tau }{a}\right){\rm{d}}t $	$ {a}^{-\frac{1}{2}}\psi \left(\dfrac{t-\tau }{a}\right) $为母小波族，$ a $为小波尺度因子，$ \tau $为时间平移因子
	广义S变换 （GST）	$ G S T\left(f,t\right)={\displaystyle\int }_{-\infty }^{+\infty }s\left(\tau \right)\dfrac{{\left|f\right|}^{\lambda }}{\sqrt{2{\text{π}} }\rho }{\rm{e}}^{-\dfrac{{\left(\tau -t\right)}^{2}{f}^{2\lambda }}{2{\rho }^{2}}}{\rm{e}}^{-j2{\text{π}} \tau }{\rm{d}}\tau $	$ \dfrac{\left|f\right|}{\sqrt{2{\text{π}} }}{\rm{e}}^{-\dfrac{{t}^{2}{f}^{2}}{2}} $为Gauss窗，$ \lambda $为时宽调节参数，$ \rho $为衰减趋势调节参数
时频后处理 方法	同步挤压小波变换（SSCWT）	$ S S C W T\left({\omega }_{l},\tau \right)={\Delta \omega }^{-1}\displaystyle\sum _{{a}_{k}:\left|C W T\left(a,\tau \right)-{\omega }_{l}\right|\leqslant \Delta \omega /2}C W T\left(a,\tau \right){{a}_{k}}^{-3/2}\left({\Delta a}_{k}\right) $	$ C W T\left(a,\tau \right) $为CWT时频谱，$ {a}_{k} $为第$ k $个离散化尺度，$ {\omega }_{l} $为点$ l $处离散化频率值
	同步挤压广义S变换（SSGST）	$ S S G S T\left({\tilde f}_{l},t\right)={(\Delta {\tilde f}_{l})}^{-1}\displaystyle\sum _{{f}_{k}:\left|{\tilde f}_{l}\left({f}_{k},t\right)-{\tilde f}_{l}\right|\leqslant \Delta {\tilde f}_{l}/2}\left|GST\left({f}_{k},t\right)\right|{f}_{k}\Delta {f}_{k} $	$ GST\left({f}_{k},t\right) $为GST时频谱，$ {f}_{k} $为离散频点，$ {\tilde f}_{l} $为挤压中心频率，$ \Delta {\tilde f}_{l} $为挤压带宽
辅助分析 方法	叠加幅值特征函数（SCCF）	$ S C C F\left(t\right)=\displaystyle\sum _{t=1}^{N}\left|TFR(f,t)\right|f=1,\cdots ,M $	$ TFR(f,t)\in {\boldsymbol{C}}^{M\cdot N} $，为时频谱矩阵
	叠加系数包络函数（SCEF）	$ S C E F\left(f\right)=\displaystyle\sum _{f=1}^{M}\left|E(f,t)\right|t=1,\cdots ,N $	$ E(f,t) $为频率点$ f $处时频谱系数的包络函数
	信息熵（Renyi）	$ {H}_{a}\left(T F R\right)=\dfrac{1}{1-a}{\mathrm{log}}_{2}\displaystyle\iint T F R(f,t){\rm{d}}t{\rm{d}}f $	$ a $为Renyi熵的阶次，一般取3，$ {H}_{a} $为熵值

下载: 导出CSV 

| 显示表格

2.   软件设计与功能实现

2.1   总体框架

软件包总体由8个函数（脚本）程序组成（表2），按照功能可划分为：1个封装有全部功能模块函数，执行全计算与绘图流程的主脚本程序；4个分别实现数据导入、预处理、小波分解与重构、同步挤压时频分析的二级功能模块函数；3个被相应功能模块函数调用，封装有核心计算方法程序的三级子调用函数。通过在MATLAB环境下运行主脚本程序Main_func.m，依次自动执行各计算模块（图1），为增强软件包处理不同信号的效果，对于信号分解、提取重构、时频变换等关键步骤，会在命令行窗口设置关键参数输入选项和相关信息提示，用户可根据设定研究目标不同，灵活输入和选择合适参数，各模块计算完成后会立即产出中间绘图图件进行反馈，便于调参及后续操作，产出最终成果图。软件整体架构和程序代码简明高效。

图 1  软件总体处理框架流程

Figure 1.  General processing framework and flow of the software package

下载: 全尺寸图片 幻灯片

表 2  软件包所含函数程序与功能

Table 2.  Summary of programs and functions of the software package

函数（脚本）名称	所属功能模块	主要功能描述
Main_func.m	主脚本程序	封装全部功能模块函数，执行全计算流程
deformation_txt_read.m	数据导入模块	导入TXT文件中原始数据及测项信息
deformation_data_interp.m	数据预处理模块	利用DCT-PLS方法平滑插值缺失数据
deformation_mallat_wavelet.m	小波分解与重构模块	利用DWT分解并提取待分析目标信号
deformation_timefrequency_analysis.m	同步挤压时频分析模块	分频计算SSCWT或SSGST时频谱，并计算SCCF函数、 SCEF函数、Renyi熵等辅助分析项
dctpls_smoothn.m	模块子调用函数	封装DCT-PLS插值方法实现算法
sscwt_tfrs.m	模块子调用函数	封装SSCWT时频方法实现算法
ssgst_tfrs.m	模块子调用函数	封装SSGST时频方法实现算法

下载: 导出CSV 

| 显示表格

2.2   数据导入模块

为方便统一数据接口和在台站间进行应用推广，软件包读取数据格式统一为中国地震前兆台网数据处理系统（QZProcess）软件下载的单测项TXT格式原始数据文件。用户可通过QZProcess软件选择和导出指定仪器、测项、时间范围和采样率下的观测数据，并进行时频分析。主脚本程序执行数据导入模块函数deformation_txt_read.m后，会自动打开并显示导入数据对话框，利用交互操作将数据导入MATLAB中，且相应文件信息会自动显示在命令行窗口，包括导入文件位置、台站信息、测项信息、采样率、数据类型、采样点长度、采样时间范围等。

2.3   数据预处理模块

由于时频计算的数据序列必须完整，因此执行数据预处理模块deformation_data_interp.m会对导入的原始数据进行完整性检查。若存在数据缺失，则利用离散余弦表示的惩罚最小二乘回归（DCT-PLS）方法对缺失数据进行无参数自动平滑插值补全（张桂欣等，2016），处理完成后会在命令行窗口自动显示是否缺数及缺失点位置序号，同时产出导入原始数据及插值结果绘图图件。

2.4   小波分解与重构模块

为提取观测数据中需分析的目标信号分量，对预处理后完整数据执行小波分解与重构模块deformation_mallat_wavelet.m，将信号中不同频率成分进行分离与重组。首先，用户需根据导入数据绘图和显示信息特征，指定小波分解参数，即在命令行窗口选择和输入小波分解层数和小波基函数，执行小波分解计算，产出原始数据小波分解结果绘图图件，在定点形变数据分析中，常将分解层数设为4～10层（张中旭等，2020），小波基函数常选用db族小波、sym族小波等（刘建明等，2016）。然后，根据小波分解结果中目标信号所在频段和时间范围分布特征，指定小波重构参数，即以首尾序号方式输入小波重构（叠加）区间和提取目标时段范围，提取目标信号，并产出信号提取结果绘图图件。

2.5   同步挤压时频计算模块

最后执行同步挤压时频计算模块deformation_timefrequency_analysis.m，对小波提取的目标信号进行处理。基于分频处理策略，软件包内置2种同步挤压时频分析方法，用户可根据目标信号频段特征自主选择采用高频分辨率较好的SSGST方法或低频分辨率较好的SSCWT方法进行时频谱计算。选择SSGST方法时，需选择和输入决定SSGST方法时频成像分辨率的高斯窗参数值，该参数由时宽调节参数和衰减趋势调节参数共同决定，一般取值范围为0～1（杨千里等，2020），可根据实际处理信号特征进行灵活调整，达到最优时频成像效果。选择SSCWT方法时，需选择和输入主要影响SSCWT方法结果的母小波类型参数和小波尺度控制参数，常用母小波包括morlet小波、cmor小波等（张维辰等，2019），小波尺度控制参数决定了结果在尺度方向上采样点数和分辨率。为统一不同采样率数据和不同计算方法的单位，时频结果中频率轴尺度采用了归一化频率。进一步利用3种辅助分析方法对计算出的同步挤压时频谱进行统计分析，最终产出目标信号时频分析成果绘图图件。

成果图件上部为显示导入原始数据及提取目标信号的时间范围（红色虚线内），下部为显示提取目标信号及其同步挤压时频谱、SCCF特征谱和SCEF包络谱，可更全面地描述观测数据中研究目标对象的时频响应特征及其在频率、时间尺度上的变化趋势，如图2所示（以西昌小庙台2021年3月DSQ型水管仪NS向观测数据为例）。

图 2  软件包处理最终成果图

Figure 2.  The result plot of the software package

下载: 全尺寸图片 幻灯片

3.   应用实例

为验证本文所提定点形变观测数据时频分析软件包的实用性，利用该软件对2020年1月1日至2021年6月30日西昌小庙台DSQ型水管仪、SS-Y型伸缩仪4个测项整时值采样观测数据进行处理，提取原始数据中固体潮频段分量信号进行时频分析，以评估观测质量。小波分解层数设为4层，小波基函数采用db5小波。小波重构区间设为细节1～细节4，即提取周期为2～32 h的信号，并将提取目标时段范围设为整个采样时间。时频计算选择SSCWT方法，并将母小波类型设为morlet小波。通过绘制的最终成果图进行对比分析。

3.1   DSQ型水管仪固体潮分析

从时频谱角度考察观测数据精度，主要依据为观测条件好、仪器运行稳定的测项可获得信噪比强的固体潮信号，其反映在时频谱上为能量主要集中在优势波群频段（吕品姬等，2011）。一般通过时频谱视觉上的定性对比确定时频能量集中性的优劣，本研究为定量评估时频谱所含信息量和复杂度，引入了Renyi熵作为参考指标，熵值越小，表明信号在时频域的能量集中度越高，反之越低。

由处理结果可知，研究时段内西昌小庙台DSQ型水管仪NS分量观测固体潮数据（图3（a））与EW分量观测固体潮数据（图3（b））在时域变化形态上均较规律，潮汐波形态清晰，时频谱中半日波、全日波等信号频率能量分布集中，除部分由地震波、降雨造成的突跳信号外，基本不存在非固体潮频段的杂波信号能量，将观测固体潮时频结果与理论固体潮时频结果（图3（c）、图3（d））进行对比，发现西昌小庙台DSQ型水管仪观测数据在时域信号形态、时频谱能量分布、SCCF特征谱频率分布、SCEF包络谱强度变化等方面一致性较好，观测固体潮时频谱Renyi熵值较小，且接近理论固体潮熵值（表3），总体上，西昌小庙台DSQ型水管仪观测精度较好，观测质量较高。

下载: 全尺寸图片 幻灯片

图 3  小庙台2020年1月至2021年6月DSQ型水管仪观测固体潮与理论固体潮时频结果对比

Figure 3.  Comparison of time-frequency results between observed and theoretical earth tide of DSQ data of Xiaomiao seismic station from Jan, 2020 to June,2021

下载: 全尺寸图片 幻灯片

表 3  观测固体潮与理论固体潮时频谱Renyi熵值

Table 3.  The Renyi entropy of observed earth tide and theoretical earth tide

项目	DSQ（NS向）	DSQ（EW向）	SS-Y（NS向）	SS-Y（EW向）
观测固体潮时频谱Renyi熵值	6.26	5.56	8.44	8.17
理论固体潮时频谱Renyi熵值	5.50	5.50	5.13	5.13
Renyi熵误差	0.76	0.06	3.31	3.04

下载: 导出CSV 

| 显示表格

3.2   SS-Y型伸缩仪固体潮分析

西昌小庙台SS-Y型伸缩仪NS分量固体潮（图4（a））与EW分量固体潮（图4（b））在研究时段内信号存在较大噪声干扰，时频谱中主要波群信号能量较分散，且被强杂波能量覆盖。由SCCF特征谱可知，相较于理论固体潮（图4（c）、图4（d））的频率分布，观测数据频率能量明显分散在固体潮中心频率以外的其他频带上，且时频谱Renyi熵值偏大，说明西昌小庙台SS-Y型伸缩仪观测精度较差，数据质量有待提高。这是因为近年来西昌小庙台应变仪器相比倾斜仪器更易受到台站周边施工、交通、机井抽水等观测环境干扰，且传感器老化严重，导致日常观测数据常出现高频噪声和波形畸变。

下载: 全尺寸图片 幻灯片

图 4  小庙台2020年1月至2021年6月SS-Y型伸缩仪观测固体潮与理论固体潮时频结果对比

Figure 4.  Comparison of time-frequency results between observed and theoretical earth tide of SS-Y data of Xiaomiao seismic station from Jan,2020 to June,2021

下载: 全尺寸图片 幻灯片

4.   结论

本研究所提基于MATLAB的定点形变观测数据时频分析软件包采用了小波分析、同步挤压时频分析和时频辅助分析等多种方法，程序基于模块化集成封装原则，以绘图-参数输入-绘图交互运行方式，实现了提取原始形变观测数据中特定频段、特定时段的目标信号分量，并对其进行高精度时频计算与自动绘图。将该软件包应用于形变观测资料的固体潮数据分析，通过对比理论固体潮，可从时频谱角度评估观测数据质量，具有一定推广应用价值。