0 引言

定点形变观测通过测量钻孔孔径变化，直接或间接地测量地壳应力—应变状态的相对变化，因而可以观测到与地震的孕育、发生直接相关的应力—应变信息，使得定点形变观测方法在地震地球物理观测中具有独特优势（苏恺之等, 2003, 2020）。“十五”期间，我国“数字观测网络建设”改造和新建约100个定点钻孔应变观测点，配备YRY-4型四元件钻孔应变仪和TJ型等体积式钻孔体应变进行形变观测，采样间隔由分钟值逐步升级到秒、赫兹，观测频带为0—3 000 Hz，观测精度优于10^-9（邱泽华，2007），而且，观测台站数量逐年增多，观测数据类型不断多样化，仪器采样率不断提升，观测数据量急剧增长，形成海量的地震前兆观测数据库（刘春国等，2017）。定点形变观测可以记录到各种应变潮汐及非应变信息，为地震预测预报研究成果的取得提供了详实的基础数据，相关研究有：邱泽华等（2010）对姑咱地震台形变资料进行小波分解，发现在2008年汶川M_S 8.1地震发生前，较长周期的异常成分出现较晚，这是破裂尺度在增大的表现；池顺良等（2013）对2013年4月20日芦山M_S 7.0地震的研究表明，距震中70 km的姑咱地震台形变资料出现明显应变异常，与该地震在时间和空间上相关；侯跃伟等（2015）对体应变、四分量应变仪等数字化形变信号进行Daubechies小波分解，发现形变数据在震前均呈现出中期、短期和临震异常特征，对周边地区震前异常判定具有一定价值。此外，一些新的地震前兆数字处理软件得到开发和应用（苏淑娟等，2019；闫如玉等，2019；杨志鹏等，2022）。

在众多语言中，MATLAB语言具有用法简单、交互性高、使用方便等特点，绘图功能强大，便于数据可视化。本文利用地震监测定点形变观测数据，借助MATLAB GUI可视化交互工具，实现定点形变数据的快速处理和可视化绘图，大幅度减小地震分析人员的工作强度，提高了工作效率，保证了数据分析的客观性。

1 系统介绍

“十五”计划以来，通过“中国数字地震监测网络”“中国大陆构造环境监测网络”和“地震背景场探测”等国家重大科学工程项目建设，建成覆盖中国大陆地区的中国地震地球物理观测台网，该台网汇集多种观测手段，集数据自动汇集、数据库式管理、仪器远程监控和数据图形化处理为一体，其中的定点形变观测前兆数据分析软件使用地震前兆信息系统（EIS 2000）和中国地震前兆台网数据处理系统（集成版）进行数据分析处理，2个软件集成了流体、形变、地电等地震地球物理学科模块，具备数据下载、异常跟踪分析、日报处理等功能，为地震数据分析提供了便捷。然而，在该软件的运行中存在以下问题：①形变观测采样率由分钟值升级为秒、赫兹（Hz），高采样大数据使得内部处理软件出现卡机、闪退等现象；②扩展性不足，针对钻孔形变频谱分析使用的小波时频图、FFT分析和STFT分析等功能缺失；③绘图功能较为单一，不能绘制特定专业图形，且图形调整功能不足（表 1）。

定点形变可视化分析系统结合大量形变观测数据，通过可视化分析处理，对于理解定点形变的震前孕育过程和变化机制具有一定意义。该系统可以对形变监测大数据进行缺值、去趋势等日常预处理操作，以确保分析数据的准确性，并就定点形变数据的震前趋势变化和地震波频谱特征，从时域和频域对进行可视化分析。同时，该系统基于模块化一键可视化交互式工作界面，通过软件自动处理分析数据，可大幅减小工作强度，提高工作效率，保证数据分析的客观性。

2 系统总体设计 2.1 系统总体框架设计

定点形变数据可视化分析系统的开发环境为Matlab R2018a，通过其开发环境下的图形用户接口GUI会生成FIG文件（.fig）和M文件（.m），二者是一一对应关系。FIG文件包含序列化的图形窗口对象，用户打开GUI时，MATLAB能够通过读取FIG文件构造新的图形窗口，对象属性被设置为在图形窗口创建时保存。M文件可分为GUI初始化和回调函数2个部分，根据用户与GUI的具体交互行为调用回调函数，进行相应控件功能的修改和调试。

实现GUI操作主要包括GUI界面设计和GUI组件编程2项工作，具体实现过程如下：①设置GUIDE应用程序选项进行GUIDE组态；②使用界面设计编辑器进行GUI界面设计，并进行修改；③添加组件，实现演示内容的参数修改和函数调用；④编写用户GUI组件行为响应控制（即回调函数）代码，并逐项进行功能测试；⑤编译并封装软件、发布。

2.2 系统模块设计

定点形变可视化分析系统设置了4个模块：①数据预处理模块：具有数据输入、选择和去突跳、去趋势和插值等预处理代码编译功能；②频谱可视化显示模块：包括傅里叶变换、小波分析时频图、短时傅里叶时频图和S变换等时频图；③结果存储模块：以JPG图片格式和txt文档格式保存数据分析结果；④界面操作模块：含清除图像、系统退出、放大图像和数据提示等操作（图 1）。

2.3 系统界面设计

根据系统功能模块设计，在GUI前台界面中创建4个功能显示模块（图 2），使用触发按钮（Pushbutton）、静态文本框（Static Text）、编辑文本框（Edit Text）和坐标轴（Axes）4种控件类型进行系统界面设计，其中前3个属于UIcontrol对象，用于各功能的触发及文字显示；坐标轴（Axes）对象主要用于数据可视化处理，使数据观测起来更加直观，其中的重要部分即对各功能进行代码编译，包括Openning Fcn（初始函数）、Callback（回调函数）和Output Fcn（输出函数）的代码设计。Openning Fcn（初始函数）用于设定各个参数的初始值，可根据实际情况设定。Callback（回调函数）是核心，是对界面控件触发时的事件响应函数；Output Fcn（输出函数）用于控制输出结果，记录算法生成数据的历史信息；最后对所有功能进行统一调试、运行。

3 主要功能代码 3.1 系统主界面代码

使用Matlab中GUIDE生成的M文件，控制GUI并决定其对用户操作的响应，包含运行GUI主界面所需代码，包括主函数、Opening函数、Output函数和回调函数，设计完成定点形变数据可视化程序主界面。该界面简单、操作友好，便于研究者使用。主要代码如下

function varargout = Fixed (varargin)

gui_Singleton = 1;

gui_State = struct(‘gui_Name’, mfilename, ...

… …

else

      gui_mainfcn(gui_State, varargin{: });

end

3.2 数据载入与处理代码

定点形变数据数据格式主要包括文本格式txt、数据格式DIG和表格xls格式等。载入数据时，文件对话框由uigetfile函数创建，通过对话框获取用户输入，返回选择的路径和文件名，便于随后对该文件进行数据操作。数据载入回调函数（callback_function）代码如下

    [filename, pathname] = uigetfile({‘*.txt; *.DIG; *.TXT; *.xls; *.dat’, ’All Data Files’; ... ‘*.*’, ’All Files’ }, ’载入数据’, ... fullfile(pwd, ‘’));

    if isequal(filename, 0) || isequal(pathname, 0)

    return;

end

    data = load(fullfile(pathname, filename), ’%s’);

3.3 图像的保存代码

系统设计了hAxes1、hAxes2、hAxes3和hAxes4四个图像显示模块，使用uiputfile函数设置保存图像的文件名和路径，并使用datestr函数获取当前时间值并返回文件名，可以根据需要下载图片并自由选择路径，图像以BMP和JPG格式保存。图像保存的回调函数（callback_function）代码如下

[fName, pName, index] = uiputfile({‘*.bmp’; ‘*.jpg’}, ‘图片另存为’, datestr(now, 30));

if index

……

    delete(gca);

        end

    elseif strcmp(strEx, ‘.jpg’)

msgbox(‘保存图像成功！’, ‘信息提示框’);

3.4 短时傅里叶变换STFT代码

短时傅里叶变换STFT是与傅里叶变换相关的一种数学变换，其中傅里叶变换主要用于平稳信号的时频转换，是信号的一种整体变换，但缺乏定位能力，而短时傅里叶变换通过在时间域加窗的方法，用以确定时变信号在局部区域正弦波的频率与相位。主要代码如下

for ii=1:nstep: N

    if(ii＜L+1) %需要将前面补零以满足窗函数的要求

……

xlabel(‘时间/s’), ylabel(‘幅度/ms’);

3.5 小波变换及小波时频图代码

小波变换作为一种信号的时间—尺度分析方法，不仅具有多分辨率分析的特点，而且在时域和频域都具有表征信号局部特征的能力，适用于探测正常信号中的局部奇异性和瞬态反常现象，并可展示这些反常现象的成分。定点形变数据属于典型非平稳信号，可以使用小波变换对非平稳信号进行分析和处理。小波变换主要代码如下

[c, l]=wavedec(x, 4, ’db4’);

d1=detcoef(c, l, 1);%提取第一层的细节分量

…

% ca4=zeros(size(a4)); %将第四层的近似系数置零

小波时频图主要代码如下主要代码如下：

wavename=’cmor3-3’; %cmor是复Morlet小波，其中3-3表示Fb-Fc，Fb是带宽参数，Fc是小波中心频率。

scals=c./(1:totalscal);

…

colorbar;

xlabel(‘时间t/s’);

ylabel(‘频率f/Hz’);

title(‘小波时频图’);

4 软件功能实现

为进一步验证定点形变数据可视化系统功能，选取2018—2022年南京市地震监测台网定点形变日常观测数据及典型地震事件记录，进行软件功能分析和演示。

4.1 预处理功能可视化实现

在钻孔应变观测中，受仪器故障或者其他因素影响导致数据漏测，或因明显过失误差和突发干扰致使某些数据不可用，造成数据的连续性、完整性、可靠性降低，增加了地震监测数据的分析难度。而定点形变可视化分析系统具备去除突跳，去趋势、一阶差分和二阶差分等数据预处理功能，可以较好地解决以上问题。自该系统在南京市地震前兆台网试运行以来，较好地完成了日常观测数据的处理和分析工作。文中给出部分示例的预处理功能实现界面，见图 3。

图 3(a)示出南京市六合台2022年4—6月钻孔体应变连续3月原始观测数据，可见由于断电、数采故障等原因，原始观测数据中出现较多空值数据（999999），导致原始观测数据曲线出现一些突跳，不能反映正常的固体潮汐变化。对于离散数据，采用MATLAT插值以获得更为丰富的信息，填充缺失数据、对现有数据进行平滑处理以及预测等。从数学上来说，数据插值是一种函数逼近方法，通过插值，所选体应变观测数据固体潮清晰，见图 3(b)。钻孔应变数据存在趋势线变化，这种变化可能来自仪器自身漂移，也存在年、月等周期性变化规律，可能代表着地壳蠕变结果，这种趋势线变换容易掩盖较细微的异常信息。因此，在对应变数据的处理和分析中，需先进行去趋势，并对等间隔数据采用detrend函数进行处理；对于非等间隔数据，则需采用一元线性回归进行去趋势处理。体应变数据经去趋势预处理，结果见图 3(c)。采用一阶差分方法，对所选钻孔应变数据进行异常分析，以降低固体潮汐的低频影响，有利于提取应变异常及与地震相关的高频异常，结果见图 3(d)。

4.2 主要频谱功能可视化实现 4.2.1 FFT可视化功能实现

FFT（fast Fouriertransform）是离散傅氏变换（DFT）的快速算法，即快速傅氏变换，时域信号经傅里叶变换分解，成为不同正弦波信号的叠加，并进行频率分析，即可将信号变换到频域。快速傅里叶变换的实质是，利用权函数W_N^nk的对称性和周期性，把N点DFT进行一系列分解和组合，将计算过程变成一系列叠代运算，简化DFT运算量，为DFT及数字信号的实时处理和应用创造了良好的条件。快速傅里叶变换算法如下

$ X(n)=\sum\limits_{K=0}^{N=1} x_0(k) \mathrm{e}^{\frac{{\mathsf{π}} n k}{N}} \quad n=0, 1, 2, \ldots, N-1 $

(1)

由式（1）可知，对每一个n，计算X(n)须作N次复数乘法及(N-1)次复数加法，要完成这组变换共需N²次乘法及N(N-1)次复数加法。

以2021年9月南京市江宁地震台钻孔定点形变数据为例，采样间隔为1 min，采用分析系统对该数据进行FFT分析，结果见图 4。据固体潮理论，固体潮主要有日潮（24 h）、半日潮（12 h）和三分之一潮（8 h)，图 4所示计算结果表明，定点形变半日潮M₂波振幅最大，对应主要优势频率区间，日潮其次，三分之一潮振幅最小，与固体潮理论一致，表明定点形变仪器观测正常。

4.2.2 短时傅里叶频谱特征可视化实现

使用短时傅里叶频谱变换进行分析时，由于信号长度有限，会导致频谱出现泄漏现象，即信号在频域上呈现明显的畸变、初值和谐波等非理想特征。为了减少这种泄漏，采用汉明窗口（Hamming）来对原始信号进行加窗处理，其主要作用是在时域上对信号进行平滑加权，抑制信号的端点泄漏，以减小分析过程中信号的泄漏现象，从而得到更精确的频域信息。以2018年9月28日印度尼西亚M_S 7.4地震的记录为例，采用短时傅里叶变换，对南京市高淳地震台钻孔分量应变数据进行频谱分析。为了克服数据的关联性，选择汉明（Hamming）窗，相邻的窗部分重叠，窗长为2 s，重叠1 s，窗口长度61，滑动步长10，对采样数据进行预处理，并进行db4小波分解，分解层数为6层，提取细节项，分析高频信号的频谱特征及能量在不同频段的分布，结果见图 5。由图 5可见，当地震波到达时，对应频点0.2 Hz和0.3 Hz处的信号能量较强，频带响应为0—0.4 Hz，之后2个强信号逐渐减弱，并恢复到正常背景。

4.2.3 小波变换与小波时频谱可视化实现

以2018年9月28日印度尼西亚M_S 7.4地震记录为例，采用小波变换与小波时频变换，分析南京市高淳地震台定点钻孔形变YRY-4钻孔分量应变数据频谱特征，结果见图 6。此次M_S 7.4地震波形正常背景曲线光滑，固体潮清晰，同震波响应明显[图 6(a)]。采用db4小波进行分解，分解层数为6层，去除原始观测曲线的趋势项，细节部分显示，此次M_S 7.4地震同震波形清晰，振幅为40×10^-3″[图 6(b)]。对细节部分进行小波时频分析，与短时傅里叶变换不同的是，小波时频图（图 7）显示，地震波到达时，在低频部分会产生虚假频率信息，频带响应为0—0.5 Hz，且在频率0.1 Hz、0.2—0.4 Hz处，对应高频信号的能量更加集中。

5 结论与展望

基于MATLAB GUI可视化交互式工具，设计开发一套基于MATLAB的定点形变可视化分析处理系统，软件设置数据导入模块、预处理模块、数据输出模块及图像保存模块，并提供常见时频分析处理方法，由Command Window和Delpoyment tool将Fixpd.m文件打包封装为可在Windows操作系统独立运行的“Fixpd.exe”程序。

该定点形变可视化分析处理系统取得一项计算机软件著作权，2021年7月，在江苏省南京市地震前兆台网试运行。定点形变数据经程序自动分析，输出预处理、小波分析等可视化频谱图像，在地震会商处理、震后资料分析中具有明显的实际意义，具有一定推广价值，目前在江苏省徐州地震台网得到推广和应用。