Design and Application of A Deformation Time-frequency Analysis Software Package Based on MATLAB
-
摘要: 为方便台站工作人员快速准确分析定点形变观测数据的时频响应特征，利用MATLAB软件研发了基于小波分析和同步挤压时频变换的交互式数据处理与成图软件包。该软件包遵循模块化设计原则，利用导入模块读取从中国地震前兆台网数据处理系统下载的原始数据，利用预处理模块对缺失数据进行插值补全，利用小波分解与重构模块从原始数据中提取待分析的目标信号分量，利用同步挤压时频分析模块对提取目标进行高精度时频分析，并在各模块关键节点中增加绘图功能，全部处理过程采用绘图-参数输入-绘图交互的方式进行，参数灵活可调，且每步计算结果直观清晰。应用该软件包对2020年1月至2021年6月西昌小庙台DSQ型水管仪和SS-Y型伸缩仪整时值采样数据进行固体潮时频计算，并与理论固体潮时频结果进行对比，结合时频辅助分析方法，从时频谱角度初步评价了2套仪器观测资料的质量情况，为台站日常数据跟踪分析提供了参考。Abstract: In order to facilitate the station staff to effectively and accurately analyze the time-frequency characteristics of the fixed point deformation data, we provide an interactive data processing and mapping software package of wavelet decomposition and synchrosqueezing time-frequency transform based on MATLAB. The package follows the principle of modularization and includes such modules as the import module reading the original data from the txt file downloaded by QZProcess, the preprocess module interpolating the missing original data, the wavelet decomposition and reconstruction module extracting the target signal component from the original data and the synchrosqueezing transform module carrying out high precision time-frequency results of the extracted target signal. Since the mapping functions are integrated into modules and the whole flow is operated interactively, the calculation parameters are adjustable and the result of each step is intuitive. We use the proposed software package to process the hourly sampling data for DSQ and SS-Y of Xiaomiao station from Jan,2020 to June,2021. The processing results of earth tide signal qualitatively evaluate the quality of observation, and also providing useful reference for precursory data tracking and analysis.
-
图 1 软件总体处理框架流程
Figure 1. General processing framework and flow of the software package
3 小庙台2020年1月至2021年6月DSQ型水管仪观测固体潮与理论固体潮时频结果对比
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
4 小庙台2020年1月至2021年6月SS-Y型伸缩仪观测固体潮与理论固体潮时频结果对比
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
表 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 软件包所含函数程序与功能
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
表 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
-
[1] 蔡佩蕊, 陈伟, 林立峰等, 2020. 基于STFT方法对沿海宽频带倾斜仪的噪声分析. 地震工程学报, 42（2）: 396—402 doi: 10.3969/j.issn.1000-0844.2020.02.396Cai P. R. , Chen W. , Lin L. F. , et al. , 2020. Noise analysis of coastal broadband tiltmeter based on STFT method. China Earthquake Engineering Journal, 42(2): 396—402. (in Chinese) doi: 10.3969/j.issn.1000-0844.2020.02.396 [2] 侯跃伟, 赵兵, 田韬, 2015. 基于Daubechies小波分析的南京数字化钻孔形变震前变化特征研究. 震灾防御技术, 10（2）: 388—396 doi: 10.11899/zzfy20150219Hou Y. W. , Zhao B. , Tian T. , 2015. Anomalies characteristics of digital data of borehole deformation before the earthquakes at the nanjing region based on daubechies wavelet analysis. Technology for Earthquake Disaster Prevention, 10(2): 388—396. (in Chinese) doi: 10.11899/zzfy20150219 [3] 黄昱丞, 郑晓东, 栾奕等, 2018. 地震信号线性与非线性时频分析方法对比. 石油地球物理勘探, 53（5）: 975—989Huang Y. C. , Zheng X. D. , Luan Y. , et al. , 2018. Comparison of linear and nonlinear time-frequency analysis on seismic signals. Oil Geophysical Prospecting, 53(5): 975—989. (in Chinese) [4] 孔祥瑞, 翟丽娜, 李梦莹等, 2018. 经验模态分解法在高频形变干扰分析中的应用. 防灾减灾学报, 34（4）: 34—37Kong X. R. , Zhai L. N. , Li M. Y. , et al. , 2018. Application of empirical mode decomposition method to analysis of high-frequency deformation interference. Journal of Disaster Prevention and Reduction, 34(4): 34—37. (in Chinese) [5] 李娜, 冯建刚, 张博等, 2020. 天水钻孔应变可靠性分析及映震效能检验. 地震工程学报, 42（5）: 1111—1116 doi: 10.3969/j.issn.1000-0844.2020.05.1111Li N. , Feng J. G. , Zhang B. , et al. , 2020. Reliability analysis of borehole strain meters and test of earthquake-reflecting efficacy in Tianshui. China Earthquake Engineering Journal, 42(5): 1111—1116. (in Chinese) doi: 10.3969/j.issn.1000-0844.2020.05.1111 [6] 刘建明, 李志海, 孙甲宁等, 2016. 基于小波分析提取地倾斜异常特征. 地震, 36（1）: 38—48 doi: 10.3969/j.issn.1000-3274.2016.01.005Liu J. M. , Li Z. H. , Sun J. N. , et al. , 2016. Extraction of ground tilt anormalies based on wavelet analysis. Earthquake, 36(1): 38—48. (in Chinese) doi: 10.3969/j.issn.1000-3274.2016.01.005 [7] 吕品姬, 赵斌, 陈志遥等, 2011. 小波分解-STFT方法在地形变观测数据中的应用. 大地测量与地球动力学, 31（5）: 136—140Lv P. J. , Zhao B. , Chen Z. Y. , et al. , 2011. Application of wavelet-decomposition and STFT method in continuous deformation observation analysis. Journal of Geodesy and Geodynamics, 31(5): 136—140. (in Chinese) [8] 孟庆筱, 吕健, 李进武等, 2018. 基于S变换的唐山四分量钻孔应力时频特征分析. 大地测量与地球动力学, 38（11）: 1202—1206Meng Q. X. , Lv J. , Li J. W. , et al. , 2018. Analysis of time-frequency characteristics of four-components borehole stress at tangshan station based on S-transform. Journal of Geodesy and Geodynamics, 38(11): 1202—1206. (in Chinese) [9] 潘晓, 曹思远, 徐彦凯等, 2020. 基于同步挤压小波变换的烃类识别技术. 地球物理学报, 63（11）: 4176—4187 doi: 10.6038/cjg2020M0663Pan X. , Cao S. Y. , Xu Y. K. , et al. , 2020. The hydrocarbon detection technology based on synchrosqueezing wavelet transform. Chinese Journal of Geophysics, 63(11): 4176—4187. (in Chinese) doi: 10.6038/cjg2020M0663 [10] 平建军, 张永仙, 单连君等, 2013. 地震前兆信息量计算软件研制及其操作说明. 震灾防御技术, 8（4）: 397—407 doi: 10.3969/j.issn.1673-5722.2013.04.007Ping J. J. , Zhang Y. X. , Shan L. J. , et al. , 2013. Software for earthquake precursor information extracting and analysis. Technology for Earthquake Disaster Prevention, 8(4): 397—407. (in Chinese) doi: 10.3969/j.issn.1673-5722.2013.04.007 [11] 孙伶俐, 李明, 蒋玲霞等, 2013. 湖北省潮汐形变观测异常及干扰识别. 大地测量与地球动力学, 33（S1）: 36—40Sun L. L. , Li M. , Jiang L. X. , et al. , 2013. Anomaly recognition of tidal deformation and disturbance factors in Hubei Province. Journal of Geodesy and Geodynamics, 33(S1): 36—40. (in Chinese) [12] 王宁, 吴云, 张燕, 2014. 时频分析方法在形变数据中的应用研究. 地震工程学报, 36（2）: 413—420 doi: 10.3969/j.issn.1000-0844.2014.02.0413Wang N. , Wu Y. , Zhang Y. , 2014. The application of time-frequency analysis methods in deformation data. China Earthquake Engineering Journal, 36(2): 413—420. (in Chinese) doi: 10.3969/j.issn.1000-0844.2014.02.0413 [13] 魏学强, 袁洪克, 秦晶晶等, 2016. 广义S变换地震信号时频分析. 震灾防御技术, 11（4）: 808—813 doi: 10.11899/zzfy20160411Wei X. Q. , Yuan H. K. , Qin J. J. , et al. , 2016. The time-frequency analysis of seismic data using generalized S transform. Technology for Earthquake Disaster Prevention, 11(4): 808—813. (in Chinese) doi: 10.11899/zzfy20160411 [14] 武善艺, 刘琦, 龚丽文等, 2018.2008年汶川MS8.0地震前定点形变高频异常特征的研究. 地震, 38（2）: 145—156Wu S. Y. , Liu Q. , Gong L. W. , et al. , 2018. High-frequency anormaly characteristics of fixed-point deformation before the 2008 Wenchuan MS8.0 earthquake. Earthquake, 38(2): 145—156. (in Chinese) [15] 严海滔, 黄饶, 周怀来等, 2019. 同步挤压广义S变换在南海油气识别中的应用. 地球物理学进展, 34（3）: 1229—1235 doi: 10.6038/pg2019CC0143Yan H. T. , Huang R. , Zhou H. L. , et al. , 2019. Application of Nanhai oil and gas identification based on synchrosqueezing generalized S transform. Progress in Geophysics, 34(3): 1229—1235. (in Chinese) doi: 10.6038/pg2019CC0143 [16] 杨千里, 王婷婷, 边银菊, 2020. 基于广义S变换的地震与爆炸识别. 地震学报, 42（5）: 613—628Yang Q. L. , Wang T. T. , Bian Y. J. , 2020. Recognition of earthquakes and explosions based on generalized S transform. Acta Seismologica Sinica, 42(5): 613—628. (in Chinese) [17] 张桂欣, 郝振纯, 祝善友等, 2016. AMSR2缺失数据重建及其土壤湿度反演精度评价. 农业工程学报, 32（20）: 137—143 doi: 10.11975/j.issn.1002-6819.2016.20.018Zhang G. X. , Hao Z. C. , Zhu S. Y. , et al. , 2016. Missing data reconstruction and evaluation of retrieval precision for AMSR2 soil moisture. Transactions of the Chinese Society of Agricultural Engineering, 32(20): 137—143. (in Chinese) doi: 10.11975/j.issn.1002-6819.2016.20.018 [18] 张维辰, 朱凯光, 池成全等, 2019. 基于小波变换的2013年芦山MS7.0地震前姑咱台钻孔应变异常时频分析. 地震学报, 41（2）: 230—238Zhang W. C. , Zhu K. G. , Chi C. Q. , et al. , 2019. Time-frequency analyses for borehole strain anomaly at guzan station before 2013 Lushan MS7.0 earthquake based on wavelet transform. Acta Seismologica Sinica, 41(2): 230—238. (in Chinese) [19] 张小艳, 熊峰, 王旭东等, 2019. 内蒙古中部地区形变主要干扰的时频响应特征分析. 中国地震, 35（4）: 718—725 doi: 10.3969/j.issn.1001-4683.2019.04.012Zhang X. Y. , Xiong F. , Wang X. D. , et al. , 2019. Analysis of time-frequency response characteristics for main deformation disturbances in central Inner Mongolia. Earthquake Research in China, 35(4): 718—725. (in Chinese) doi: 10.3969/j.issn.1001-4683.2019.04.012 [20] 张中旭, 李智蓉, 2020. 小波分析方法在定点形变日常跟踪中的应用研究. 云南大学学报（自然科学版）, 42（6）: 1121—1128Zhang Z. X. , Li Z. R. , 2020. Application of wavelet analysis method in daily tracking of fixed-point deformation. Journal of Yunnan University (Natural Sciences Edition), 42(6): 1121—1128. (in Chinese) [21] 赵佳佳, 陈志遥, 张燕等, 2017. 奇异谱分析在倾斜应变数据处理中的应用研究. 大地测量与地球动力学, 37（5）: 541—545Zhao J. J. , Chen Z. Y. , Zhang Y. , et al. , 2017. Research on application of singular spectrum analysis in data processing of earth tilt and strain. Journal of Geodesy and Geodynamics, 37(5): 541—545. (in Chinese) [22] Mallat S. G. , 1989. A theory for multiresolution signal decomposition: the wavelet representation. IEEE Transactions on Pattern Analysis and Machine Intelligence, 11(7): 674—693. doi: 10.1109/34.192463 [23] Mousavi S. M. , Langston C. A. , Horton S. P. , 2016. Automatic microseismic denoising and onset detection using the synchrosqueezed continuous wavelet transform. Geophysics, 81(4): V341—V355. doi: 10.1190/geo2015-0598.1 -
点击查看大图
图(6) / 表(3)
计量
- 文章访问数: 138
- HTML全文浏览量: 58
- PDF下载量: 30
- 被引次数: 0
