基于MATLAB 的定點形變觀測數(shù)據(jù)時頻分析軟件設(shè)計及應(yīng)用研究1

楊志鵬 陳秀清 張御陽 顏 歡 陳碧洪 阮 祥

（四川省地震局, 成都 610041）

引言

地形變測量中的定點潮汐形變觀測主要以地傾斜、地應(yīng)變等測量指標為主，其目的在于記錄地殼周期性潮汐形變，并能捕捉地震中短臨前兆異常信息（孫伶俐等，2013）。但實際形變觀測結(jié)果受多方面因素影響，主要包括儀器系統(tǒng)狀態(tài)、固體潮汐、觀測環(huán)境、人為因素、自然環(huán)境、地震波、地下介質(zhì)變化等（蔡佩蕊等，2020）。為分析和判識觀測中的各類干擾信號、潮汐波及地震前兆異常信息等問題，近年來，諸多數(shù)字信號處理技術(shù)被引入形變資料分析領(lǐng)域，如小波分析方法（侯躍偉等，2015）、模態(tài)分解類方法（孔祥瑞等，2018）、超限率分析方法（李娜等，2020）、奇異譜分析方法（趙佳佳等，2017）、時頻分析方法（王寧等，2014）等。其中，時頻分析方法以其能夠表征信號頻譜隨時間變化規(guī)律的優(yōu)勢，在形變觀測干擾分析（張小艷等，2019）、前兆異常判定（武善藝等，2018）、日常跟蹤分析（孟慶筱等，2018）等方面得到了應(yīng)用與發(fā)展。目前，制約時頻分析方法進一步在臺站形變數(shù)據(jù)處理中的應(yīng)用主要問題為：①傳統(tǒng)的連續(xù)小波變換（CWT）、短時傅里葉變換（STFT）、S 變換（ST）、希爾伯特-黃變換（HHT）等方法時頻分辨率較低，對于波形類別多、分布頻段廣、振幅差異大的各類形變觀測信號的泛化適用能力較差；②臺站缺乏相關(guān)集成處理軟件程序（平建軍等，2013），數(shù)據(jù)導(dǎo)入、函數(shù)傳遞、方法計算、成果繪圖等功能和步驟缺乏統(tǒng)一設(shè)計規(guī)劃，不便于操作與推廣。

為解決上述問題，本研究基于MATLAB 科學(xué)計算平臺研發(fā)了定點形變觀測數(shù)據(jù)時頻分析軟件，該軟件通過集成數(shù)據(jù)導(dǎo)入、插值補全、小波分析、同步擠壓時頻變換、自動繪圖等功能，處理過程采用命令窗參數(shù)繪圖交互方式，能夠靈活計算特定頻段、特定時段形變觀測數(shù)據(jù)的高精度時頻譜，以期為臺站工作人員深入挖掘形變觀測資料中的各類潛在變化特征提供簡便實用的處理工具。

1 方法原理

本研究采用信號分解方法中的小波分析（DWT），時頻處理方法中基于CWT 和廣義S 變換（GST）的同步擠壓小波變換（SSCWT）和同步擠壓廣義S 變換（SSGST），以及時頻輔助分析方法中的疊加幅值特征函數(shù)（SCCF）、疊加系數(shù)包絡(luò)函數(shù)（SCEF）和信息熵（Renyi），全面刻畫不同頻段分布的復(fù)雜形變觀測數(shù)據(jù)時頻響應(yīng)特征。

DWT 是利用選定小波基和分解層數(shù)以二進塔式分解方式將原始信號分解到一系列相互獨立的頻帶上（Mallat，1989），凸顯信號中不同頻段成分信息，并可通過篩選重構(gòu)得到特定優(yōu)勢頻段子信號分量的方法。

CWT 是將原始信號與經(jīng)過伸縮平移后的母小波作褶積運算的方法，具有多尺度表達信號頻率分辨率特征，其時頻譜結(jié)果頻率軸尺度為指數(shù)型分布，對于中低頻段信號具有較好的時頻響應(yīng)能力。

GST 是對傳統(tǒng)ST 方法高斯窗函數(shù)引入2 個可調(diào)參數(shù)（魏學(xué)強等，2016），使其窗函數(shù)隨頻率呈現(xiàn)出更多非線性變化特征而提高時頻分辨率的方法，其時頻譜結(jié)果頻率軸尺度為線性分布，對于高頻段信號具有較好的時頻響應(yīng)能力。

SSCWT、SSGST 是分別在CWT 和GST“前處理”結(jié)果的基礎(chǔ)上，進一步利用基于微分算子的擠壓“后處理”操作，沿頻率軸重排時頻譜瞬時頻帶，將特定頻率區(qū)間內(nèi)的時頻系數(shù)疊加到理論頻率點上（嚴海滔等，2019；潘曉等，2020），可得到時頻聚集性更精確的同步擠壓時頻結(jié)果的方法。

SCCF、SCEF 是分別描述刻畫同步擠壓時頻譜在頻率尺度上的能量分布特征和在時間尺度上的強度分布特征（Mousavi 等，2016）的指標。

Renyi 熵是定量檢測時頻譜所含信息量和時頻譜能量集中度的指標（黃昱丞等，2018）。

整個方法計算流程為：（1）利用DWT 方法對原始形變信號作多尺度分解；（2）根據(jù)設(shè)定的研究目標，對分解出的子信號及時間范圍作選擇性重構(gòu)，以提取待分析的目標信號分量；（3）根據(jù)目標分量的頻段分布特征，基于分頻處理策略，對于短周期高頻目標信號和長周期低頻目標信號分別采用高頻響應(yīng)較好的GST 方法和中低頻響應(yīng)較好的CWT 方法進行時頻計算，在此基礎(chǔ)上利用同步擠壓操作得到更高分辨率的SSGST 和SSCWT 時頻結(jié)果；（4）進一步對計算出的同步擠壓時頻譜結(jié)果，采用SCCF 方法、SCEF 方法和Renyi 熵分別從頻率尺度、時間尺度、能量復(fù)雜度方面進行輔助分析。

上述計算方法基本公式和控制參數(shù)如表1 所示。

表1 軟件包所用方法及基本定義Table 1 Summary of methods and basic definitions of the software package

2 軟件設(shè)計與功能實現(xiàn)

2.1 總體框架

軟件包總體由8 個函數(shù)（腳本）程序組成（表2），按照功能可劃分為：1 個封裝有全部功能模塊函數(shù)，執(zhí)行全計算與繪圖流程的主腳本程序；4 個分別實現(xiàn)數(shù)據(jù)導(dǎo)入、預(yù)處理、小波分解與重構(gòu)、同步擠壓時頻分析的二級功能模塊函數(shù)；3 個被相應(yīng)功能模塊函數(shù)調(diào)用，封裝有核心計算方法程序的三級子調(diào)用函數(shù)。通過在MATLAB 環(huán)境下運行主腳本程序Main_func.m，依次自動執(zhí)行各計算模塊（圖1），為增強軟件包處理不同信號的效果，對于信號分解、提取重構(gòu)、時頻變換等關(guān)鍵步驟，會在命令行窗口設(shè)置關(guān)鍵參數(shù)輸入選項和相關(guān)信息提示，用戶可根據(jù)設(shè)定研究目標不同，靈活輸入和選擇合適參數(shù)，各模塊計算完成后會立即產(chǎn)出中間繪圖圖件進行反饋，便于調(diào)參及后續(xù)操作，產(chǎn)出最終成果圖。軟件整體架構(gòu)和程序代碼簡明高效。

圖1 軟件總體處理框架流程Fig. 1 General processing framework and flow of the software package

表2 軟件包所含函數(shù)程序與功能Table 2 Summary of programs and functions of the software package

2.2 數(shù)據(jù)導(dǎo)入模塊

為方便統(tǒng)一數(shù)據(jù)接口和在臺站間進行應(yīng)用推廣，軟件包讀取數(shù)據(jù)格式統(tǒng)一為中國地震前兆臺網(wǎng)數(shù)據(jù)處理系統(tǒng)（QZProcess）軟件下載的單測項TXT 格式原始數(shù)據(jù)文件。用戶可通過QZProcess 軟件選擇和導(dǎo)出指定儀器、測項、時間范圍和采樣率下的觀測數(shù)據(jù)，并進行時頻分析。主腳本程序執(zhí)行數(shù)據(jù)導(dǎo)入模塊函數(shù)deformation_txt_read.m 后，會自動打開并顯示導(dǎo)入數(shù)據(jù)對話框，利用交互操作將數(shù)據(jù)導(dǎo)入MATLAB 中，且相應(yīng)文件信息會自動顯示在命令行窗口，包括導(dǎo)入文件位置、臺站信息、測項信息、采樣率、數(shù)據(jù)類型、采樣點長度、采樣時間范圍等。

2.3 數(shù)據(jù)預(yù)處理模塊

由于時頻計算的數(shù)據(jù)序列必須完整，因此執(zhí)行數(shù)據(jù)預(yù)處理模塊deformation_data_interp.m 會對導(dǎo)入的原始數(shù)據(jù)進行完整性檢查。若存在數(shù)據(jù)缺失，則利用離散余弦表示的懲罰最小二乘回歸（DCT-PLS）方法對缺失數(shù)據(jù)進行無參數(shù)自動平滑插值補全（張桂欣等，2016），處理完成后會在命令行窗口自動顯示是否缺數(shù)及缺失點位置序號，同時產(chǎn)出導(dǎo)入原始數(shù)據(jù)及插值結(jié)果繪圖圖件。

2.4 小波分解與重構(gòu)模塊

為提取觀測數(shù)據(jù)中需分析的目標信號分量，對預(yù)處理后完整數(shù)據(jù)執(zhí)行小波分解與重構(gòu)模塊deformation_mallat_wavelet.m，將信號中不同頻率成分進行分離與重組。首先，用戶需根據(jù)導(dǎo)入數(shù)據(jù)繪圖和顯示信息特征，指定小波分解參數(shù)，即在命令行窗口選擇和輸入小波分解層數(shù)和小波基函數(shù)，執(zhí)行小波分解計算，產(chǎn)出原始數(shù)據(jù)小波分解結(jié)果繪圖圖件，在定點形變數(shù)據(jù)分析中，常將分解層數(shù)設(shè)為4～10 層（張中旭等，2020），小波基函數(shù)常選用db 族小波、sym 族小波等（劉建明等，2016）。然后，根據(jù)小波分解結(jié)果中目標信號所在頻段和時間范圍分布特征，指定小波重構(gòu)參數(shù)，即以首尾序號方式輸入小波重構(gòu)（疊加）區(qū)間和提取目標時段范圍，提取目標信號，并產(chǎn)出信號提取結(jié)果繪圖圖件。

2.5 同步擠壓時頻計算模塊

最后執(zhí)行同步擠壓時頻計算模塊deformation_timefrequency_analysis.m，對小波提取的目標信號進行處理。基于分頻處理策略，軟件包內(nèi)置2 種同步擠壓時頻分析方法，用戶可根據(jù)目標信號頻段特征自主選擇采用高頻分辨率較好的SSGST 方法或低頻分辨率較好的SSCWT 方法進行時頻譜計算。選擇SSGST 方法時，需選擇和輸入決定SSGST 方法時頻成像分辨率的高斯窗參數(shù)值，該參數(shù)由時寬調(diào)節(jié)參數(shù)和衰減趨勢調(diào)節(jié)參數(shù)共同決定，一般取值范圍為0～1（楊千里等，2020），可根據(jù)實際處理信號特征進行靈活調(diào)整，達到最優(yōu)時頻成像效果。選擇SSCWT 方法時，需選擇和輸入主要影響SSCWT 方法結(jié)果的母小波類型參數(shù)和小波尺度控制參數(shù)，常用母小波包括morlet 小波、cmor 小波等（張維辰等，2019），小波尺度控制參數(shù)決定了結(jié)果在尺度方向上采樣點數(shù)和分辨率。為統(tǒng)一不同采樣率數(shù)據(jù)和不同計算方法的單位，時頻結(jié)果中頻率軸尺度采用了歸一化頻率。進一步利用3 種輔助分析方法對計算出的同步擠壓時頻譜進行統(tǒng)計分析，最終產(chǎn)出目標信號時頻分析成果繪圖圖件。

成果圖件上部為顯示導(dǎo)入原始數(shù)據(jù)及提取目標信號的時間范圍（紅色虛線內(nèi)），下部為顯示提取目標信號及其同步擠壓時頻譜、SCCF 特征譜和SCEF 包絡(luò)譜，可更全面地描述觀測數(shù)據(jù)中研究目標對象的時頻響應(yīng)特征及其在頻率、時間尺度上的變化趨勢，如圖2 所示（以西昌小廟臺2021 年3 月DSQ 型水管儀NS 向觀測數(shù)據(jù)為例）。

圖2 軟件包處理最終成果圖Fig. 2 The result plot of the software package

3 應(yīng)用實例

為驗證本文所提定點形變觀測數(shù)據(jù)時頻分析軟件包的實用性，利用該軟件對2020 年1 月1 日至2021年6 月30 日西昌小廟臺DSQ 型水管儀、SS-Y 型伸縮儀4 個測項整時值采樣觀測數(shù)據(jù)進行處理，提取原始數(shù)據(jù)中固體潮頻段分量信號進行時頻分析，以評估觀測質(zhì)量。小波分解層數(shù)設(shè)為4 層，小波基函數(shù)采用db5小波。小波重構(gòu)區(qū)間設(shè)為細節(jié)1～細節(jié)4，即提取周期為2～32 h 的信號，并將提取目標時段范圍設(shè)為整個采樣時間。時頻計算選擇SSCWT 方法，并將母小波類型設(shè)為morlet 小波。通過繪制的最終成果圖進行對比分析。

3.1 DSQ 型水管儀固體潮分析

從時頻譜角度考察觀測數(shù)據(jù)精度，主要依據(jù)為觀測條件好、儀器運行穩(wěn)定的測項可獲得信噪比強的固體潮信號，其反映在時頻譜上為能量主要集中在優(yōu)勢波群頻段（呂品姬等，2011）。一般通過時頻譜視覺上的定性對比確定時頻能量集中性的優(yōu)劣，本研究為定量評估時頻譜所含信息量和復(fù)雜度，引入了Renyi 熵作為參考指標，熵值越小，表明信號在時頻域的能量集中度越高，反之越低。

由處理結(jié)果可知，研究時段內(nèi)西昌小廟臺DSQ 型水管儀NS 分量觀測固體潮數(shù)據(jù)（圖3（a））與EW分量觀測固體潮數(shù)據(jù)（圖3（b））在時域變化形態(tài)上均較規(guī)律，潮汐波形態(tài)清晰，時頻譜中半日波、全日波等信號頻率能量分布集中，除部分由地震波、降雨造成的突跳信號外，基本不存在非固體潮頻段的雜波信號能量，將觀測固體潮時頻結(jié)果與理論固體潮時頻結(jié)果（圖3（c）、圖3（d））進行對比，發(fā)現(xiàn)西昌小廟臺DSQ 型水管儀觀測數(shù)據(jù)在時域信號形態(tài)、時頻譜能量分布、SCCF 特征譜頻率分布、SCEF 包絡(luò)譜強度變化等方面一致性較好，觀測固體潮時頻譜Renyi 熵值較小，且接近理論固體潮熵值（表3），總體上，西昌小廟臺DSQ 型水管儀觀測精度較好，觀測質(zhì)量較高。

圖3 小廟臺2020 年1 月至2021 年6 月DSQ 型水管儀觀測固體潮與理論固體潮時頻結(jié)果對比Fig. 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

3.2 SS-Y 型伸縮儀固體潮分析

西昌小廟臺SS-Y 型伸縮儀NS 分量固體潮（圖4（a））與EW 分量固體潮（圖4（b））在研究時段內(nèi)信號存在較大噪聲干擾，時頻譜中主要波群信號能量較分散，且被強雜波能量覆蓋。由SCCF 特征譜可知，相較于理論固體潮（圖4（c）、圖4（d））的頻率分布，觀測數(shù)據(jù)頻率能量明顯分散在固體潮中心頻率以外的其他頻帶上，且時頻譜Renyi 熵值偏大，說明西昌小廟臺SS-Y 型伸縮儀觀測精度較差，數(shù)據(jù)質(zhì)量有待提高。這是因為近年來西昌小廟臺應(yīng)變儀器相比傾斜儀器更易受到臺站周邊施工、交通、機井抽水等觀測環(huán)境干擾，且傳感器老化嚴重，導(dǎo)致日常觀測數(shù)據(jù)常出現(xiàn)高頻噪聲和波形畸變。

圖4 小廟臺2020 年1 月至2021 年6 月SS-Y 型伸縮儀觀測固體潮與理論固體潮時頻結(jié)果對比Fig. 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 結(jié)論

本研究所提基于MATLAB 的定點形變觀測數(shù)據(jù)時頻分析軟件包采用了小波分析、同步擠壓時頻分析和時頻輔助分析等多種方法，程序基于模塊化集成封裝原則，以繪圖-參數(shù)輸入-繪圖交互運行方式，實現(xiàn)了提取原始形變觀測數(shù)據(jù)中特定頻段、特定時段的目標信號分量，并對其進行高精度時頻計算與自動繪圖。將該軟件包應(yīng)用于形變觀測資料的固體潮數(shù)據(jù)分析，通過對比理論固體潮，可從時頻譜角度評估觀測數(shù)據(jù)質(zhì)量，具有一定推廣應(yīng)用價值。