周曉華，王曉丹，賈芳芳，劉南南

(吉林大學儀器科學與電氣工程學院，長春130026)

0 引言

微動是地球表層時刻存在的微弱振動，主要由人為因素和自然因素等隨機源引起。研究表明，微動中攜帶與地殼淺部介質(zhì)密切相關(guān)的信息，可以利用微動進行場地評價、無損檢測等。由于微動本身具有一些優(yōu)良特性，人們不斷采集微動數(shù)據(jù)，并對其進行分析。但實驗獲取的數(shù)據(jù)不易進行處理，如，隨著地震觀測數(shù)字化進程的不斷發(fā)展，使地震觀測數(shù)據(jù)中包含了大量地殼活動的信息。在研究過程中需要對大量的地震記錄數(shù)據(jù)進行瀏覽分析，既費時又費力［1］。目前國內(nèi)外學者對微動信號進行了多方面研究，也提出了各種數(shù)據(jù)處理方法，但尚未有完善的微動信號波形特征分析軟件，即使有，也沒進行系統(tǒng)地整理并形成軟件。鑒于此，筆者分析微動數(shù)據(jù)的波形特征及相關(guān)處理方法并應(yīng)用Matlab編寫相應(yīng)程序，實現(xiàn)對微動信號記錄數(shù)據(jù)進行自動特征拾取、記錄和分析，并采用C#語言編制相關(guān)軟件，具有形象、直觀、可擴展性，為后續(xù)信號處理提供技術(shù)支撐［2］。該信號特征自動分析系統(tǒng)可拓展到地學、醫(yī)學、軍事和經(jīng)濟等領(lǐng)域，具有良好的應(yīng)用前景。

1 系統(tǒng)總體設(shè)計

筆者設(shè)計的微動信號波形特征自動分析系統(tǒng)，由微動信號采集模塊、歸一化處理模塊、時域分析模塊和頻域分析模塊組成，系統(tǒng)總體結(jié)構(gòu)圖如圖1所示。微動信號采集模塊用來采集實驗測得的微動信號數(shù)據(jù)，將這些數(shù)據(jù)讀取到程序中，進行數(shù)據(jù)處理［3］。歸一化數(shù)據(jù)處理模塊即對數(shù)據(jù)進行歸一化處理，使數(shù)據(jù)歸一化分布在［－1，1］區(qū)間內(nèi)。隨機過程的概率分布完全刻畫了其統(tǒng)計特性，但在實際中，難以確定隨機過程的概率分布，而確定其統(tǒng)計參數(shù)，進行時域和頻域分析要容易些［4］。時域分析模塊功能包括標準差、方差和期望的計算并列表，得到自相關(guān)性圖像和互相關(guān)性圖像，對自相關(guān)性和互相關(guān)性進行分析。頻域分析模塊包括頻譜分析和功率譜分析。

圖1 系統(tǒng)總體結(jié)構(gòu)圖Fig.1 Structure spectrum of the system

2 系統(tǒng)理論算法分析

2.1 數(shù)學期望與方差

一般地，隨機信號滿足正態(tài)分布與遍歷性，可用時間方向的特征表示空間方向特征，根據(jù)EX、DX數(shù)據(jù)可以分析微動信號特征［5］。設(shè)微動信號為X，取樣信號點數(shù)為N，第i個信號為Xi，則數(shù)學期望EX與方差DX可以表示為

2.2 自相關(guān)函數(shù)

為研究微動信息在時間方向的重復性，進行自相關(guān)計算。求得自相關(guān)圖譜，通過分析自相關(guān)譜，進一步分析微動信號時域特征。自相關(guān)函數(shù)是描述隨機信號f(t)在時刻0，T的取值之間相隔數(shù)據(jù)間的相關(guān)程度，自相關(guān)函數(shù)P(t)可以表示為

2.3 互相關(guān)函數(shù)

互相關(guān)函數(shù)表示兩個時間序列之間和同一個時間序列在任意兩個不同時刻取值之間的相關(guān)程度，即互相關(guān)函數(shù)是描述隨機信號f(t)，g(t)在任意兩個不同時刻t1，t2的取值之間的相關(guān)程度［6］。描述兩個不同的信號之間相關(guān)性的函數(shù)，t1取值為0，t2取值為T，互相關(guān)函數(shù)R(t)可表示為

2.4 頻譜

去掉信號中的直流分量后，可用快速傅里葉變換FFT(Fast Fourier Transform)求出原始信號的頻譜圖［7］。

2.5 功率譜

該系統(tǒng)主要應(yīng)用周期圖法、周期法的改進法、自相關(guān)法和Welch法求取功率譜。

周期圖法:功率譜估計的周期圖法(即直接法)，是通過計算出N個樣本信號序列X(n)的傅里葉變換，將得到的變換結(jié)果幅值的平方除以樣本序列的長度N，把最終得到的結(jié)果作為真實功率譜的估計。

周期法的改進法:周期法的改進法(Bartlett)主要分為兩個方向發(fā)展:1)可將長度為N的數(shù)據(jù)分成若干段(其區(qū)域可為重疊段，也可不重疊段)，分別求出每段的功率譜后相加再做平均，即分段平均周期圖法;2)不用矩形窗，而使用更加合適的窗函數(shù)消去由矩形窗旁瓣帶來的功率譜模糊和搶走，即為加窗平均周期圖法。

自相關(guān)法求功率譜:根據(jù)維納－辛欽定理，獲得脈動信號的自相關(guān)函數(shù)，則可求出功率譜密度。

Welch法求功率譜:Welch法也稱為加窗重疊平均周期圖法，是對Bartlett法的改進，它對信號序列重疊分段，同時采用加窗函數(shù)和FFT算法計算1個信號序列的自功率譜估計和兩個信號序列的互功率譜估計［8］。

2.6 H/V譜

波速法(近似法)是求取卓越周期的一種方法。加權(quán)平均波速法推算的場地卓越周期一般比真值小20%左右，而子層周期求和法推算的結(jié)果則偏大，約比真值大20%。一般地，加權(quán)平均波速法的計算公式為［9］

其中TG1為場地卓越周期，H為深度，Vs為剪切波速。

微動實測數(shù)據(jù)的H/V譜定義為角頻率為w的水平方向(NS和EW)的功率譜H與垂直方向(UD)的功率譜V之比，其表達式如下

其中PUD(w)表示垂直方向的Fourier功率譜，PNS(w)和PEW(w)分別表示水平方向上的兩分量的Fourier功率譜。

3 系統(tǒng)軟件設(shè)計

系統(tǒng)軟件設(shè)計流程圖如圖2所示。程序開始讀取數(shù)據(jù)，進行參數(shù)設(shè)置，然后采樣、進行歸一化。判斷是否進行時域分析，功能包括求標準差、期望和方差，并將求得的數(shù)值傳到Excel表里，畫自相關(guān)性圖像和互相關(guān)性圖像。之后判斷是否進行頻域分析，功能包括繪制經(jīng)FFT后的頻譜圖和功率譜圖，并進行頻譜分析。

圖2 系統(tǒng)軟件設(shè)計流程圖Fig.2 System software design flow

3.1 微動信號采集模塊

利用測量設(shè)備，收集在不同地理位置實測得到的微動信號數(shù)據(jù)。以一點數(shù)據(jù)處理為例(見圖3)，測量獲取1546號、1645號、1636號、1654號、1606號、1628號和1651號信號數(shù)據(jù)，通過這些數(shù)據(jù)可以分析1546號信號數(shù)據(jù)波形特征。編寫一個數(shù)據(jù)讀取函數(shù)，通過函數(shù)調(diào)用得到需要處理的信號數(shù)據(jù)。實驗測得大量的數(shù)據(jù)，為區(qū)分這些數(shù)據(jù)，每當讀取函數(shù)調(diào)用一組數(shù)據(jù)，便可顯示該組數(shù)據(jù)的采集時間、地點和方向(見圖5)。通過對讀取的數(shù)據(jù)采樣，可實現(xiàn)對實驗數(shù)據(jù)的進一步處理。

圖3 微動信號采集示意圖Fig.3 Microtremor collection spectrum

3.2 歸一化處理模塊

歸一化是將各個微動信號自身的數(shù)據(jù)特性通過數(shù)學的方法處理掉，從而可以使用一些更通用性的解決方法，將信號處理從一點轉(zhuǎn)到整體。對讀取的數(shù)據(jù)進行歸一化，以便之后的數(shù)據(jù)處理。Matlab中的歸一化處理有3種方法，可用premnmx語句進行歸一化。premnmx語句的語法格式是［Pn，Pmin，Pmax，Tn，tmin，tmax］=premnmx(P，T)，其中P，T分別為原始輸入和輸出數(shù)據(jù)，Pmin和Pmax分別為P中的最小值和最大值。tmin和tmax分別為T的最小值和最大值。premnmx函數(shù)用于對網(wǎng)絡(luò)的輸入數(shù)據(jù)或輸出數(shù)據(jù)進行歸一化，數(shù)據(jù)將分布在［－1，1］區(qū)間內(nèi)。

3.3 時域分析模塊

設(shè)定采樣間隔，對實測數(shù)據(jù)進行采樣。這里取某地的實測數(shù)據(jù)，計算實測數(shù)據(jù)的標準差、方差值和期望值，并列表進行比較。從數(shù)據(jù)計算可以判斷所獲取信號是否具有非平穩(wěn)性。

同時可以進行相關(guān)性分析，得到數(shù)據(jù)的自相關(guān)性和互相關(guān)性圖像。

3.4 頻域分析模塊

對實測信號進行頻域分析，畫出經(jīng)FFT變換后的信號頻譜圖［10］。定義采樣頻率、觀測半徑，利用freq函數(shù)求頻散曲線頻率點［11］，利用length函數(shù)求取頻散曲線的數(shù)據(jù)點數(shù)，調(diào)用自定義的數(shù)據(jù)讀取函數(shù)讀取中心測點的微動信號數(shù)據(jù)，之后存放中心測點的觀測數(shù)據(jù)，求取觀測數(shù)據(jù)長度，定義觀測數(shù)據(jù)的時間向量。對觀測半徑為200 m的數(shù)據(jù)進行處理，首先讀取觀測半徑為200 m圓周上的數(shù)據(jù)，從實測數(shù)據(jù)得到頻散曲線。利用同樣的方法再讀取和處理半徑為100 m圓周上的數(shù)據(jù)，再利用函數(shù)對數(shù)據(jù)進行歸一化處理，讀取不同長度的基本數(shù)據(jù)段，求取相應(yīng)數(shù)據(jù)段的功率譜。

4 軟件測試結(jié)果

4.1 采集信號

利用數(shù)據(jù)讀取函數(shù)導入實測信號，采用統(tǒng)計方法檢驗信號的穩(wěn)定性。該測試以某地采集的實測信號數(shù)據(jù)為測試信號數(shù)據(jù)。所采集的實測信號共有3組數(shù)據(jù)，分別為第1道東西、第2道南北和第3道垂直數(shù)據(jù)。每道數(shù)據(jù)又分為14組實測數(shù)據(jù)，共42組數(shù)據(jù)。分別檢驗信號幅度(原始記錄):觀察幅值變化情況，判斷其改變的劇烈程度［12］。

通過數(shù)據(jù)采集模塊調(diào)入野外實驗獲得的微動實測信號，其中一測點的信號記錄如圖4所示。

圖4 測點的實測信號記錄Fig.4 The measured signal of measurement points record

4.2 時域與頻域分析結(jié)果

采樣間隔取20 ms，分別對天津?qū)崪y的42組數(shù)據(jù)進行采樣。計算實測數(shù)據(jù)的標準差、方差值和期望值，并列進行比較(見表1)。從數(shù)據(jù)計算結(jié)果可以判斷，所取的3組信號具有非平穩(wěn)性。圖5是原始信號時域圖，由圖5可知，第1道東西第1組數(shù)據(jù)的時域曲線圖其幅值的非平穩(wěn)性。圖6和圖7分別為實測數(shù)據(jù)“SN060D20111104_120000.sg2”自相關(guān)性圖像和互相關(guān)圖像。頻域分析結(jié)果如圖8～圖10所示。

表1 數(shù)據(jù)特征值表Tab.1 Data eigenvalues sheet

圖5 原始信號時域圖Fig.5 The original signal time－domain spectrum

圖6 自相關(guān)性圖Fig.6 Autocorrelation spectrum

圖7 互相關(guān)性圖Fig.7 Cross－correlation spectrum

圖8 經(jīng)變換的信號頻譜圖Fig.8 The transformed signal spectrum

從圖8中可以看出，信號頻率分布在低頻約4 Hz左右，符合微動信號頻率分布特征，微動信號頻率分布特征和頻譜分析結(jié)果一致。

通過求取其峰值得到場地的卓越周期，經(jīng)測試，帶入5組數(shù)據(jù)最后得出該地層結(jié)構(gòu)的卓越周期為0.405 7 s，在0.105 s～2 s的范圍內(nèi)，符合實際地層要求。

圖9 功率譜圖Fig.9 The power spectrum

圖10 H/V圖Fig.10 H/V spectrum

4.3 軟件界面

軟件界面人性化，操作簡單，結(jié)果清晰，具有統(tǒng)計特性界面、處理前后對比界面、功率譜不同方法對比界面、卓越周期界面和數(shù)據(jù)處理界面(見圖11～圖15)。

圖11 統(tǒng)計特性界面圖Fig.11 Statistical properties of the software interface

圖12 處理前后對比界面圖Fig.12 Processing comparison of the software interface

圖13 功率譜不同方法對比界面圖Fig.13 Power with different methods comparison of the software interface

圖14 卓越周期界面圖Fig.14 Excellence periodogram of the software interface

圖15 數(shù)據(jù)處理界面圖Fig.15 Data processing of the software interface

4 結(jié) 語

筆者設(shè)計了微動信號波形特征自動分析系統(tǒng)，實現(xiàn)了應(yīng)用計算機自動對微動信號記錄數(shù)據(jù)進行異常特征拾取、記錄和分析，形成微動信號特征波形自動分析系統(tǒng)，使長時間連續(xù)觀測的數(shù)據(jù)得以及時處理，簡單方便，省時省力。通過編寫不同函數(shù)求取不同數(shù)據(jù)規(guī)律，進而分析微動信號的主要特征，為表征場地動力學特性提供技術(shù)參考，可應(yīng)用于勘探、地震監(jiān)測、房屋建設(shè)和地層結(jié)構(gòu)等研究。

［1］蘇維娜，王春澤，曾新森，等.基于m序列的可控震源編碼脈沖信號研究［J］.吉林大學學報:信息科學版，2012，30(6):569－572.SU Weina，WANG Chunze，ZENG Xinsen，et al.Design of Vibroseis Coding Pulse Signal Based on m Sequence［J］.Journal of Jilin University:Information Science Edition，2012，30(6):569－572.

［2］趙義平，王春民，吳海燕，等.核磁共振地下水探測系統(tǒng)信號源的研制［J］.吉林大學學報:信息科學版，2011，29(4):290－296.ZHAO Yiping，WANG Chunmin，WU Haiyan，et al.Development of Magnetic Resonance Sounding Signal Generator for Groundwater Detector［J］.Journal of Jilin University:Information Science Edition，2011，29(4):290－296.

［3］奧本海默.信號系統(tǒng)［M］.2版.北京:電子工業(yè)出版社，2003.OPPENHEIMER.Signals and Systems［M］.2nd ed.Beijing:The Electronic Industrial Press，2003.

［4］李杰.多道瞬態(tài)瑞雷波數(shù)字處理軟件及其應(yīng)用研究［D］.武漢:中國地質(zhì)大學資源學院，2010.LI Jie.The Study of Digital Processing Software and Application in Multi－Channel Transient Rayleigh Wave Data ［D］.Wuhan:College of Geosciences Resources，China Geological University，2010.

［5］王小平，曹立明.遺傳算法－理論、應(yīng)用與軟件實現(xiàn)［M］.西安:西安交通大學，2002.WANG Xiaoping，CAO Liming.Genetic Algorithm Theory，Application and Software Implementation［M］.Xi'an:Xi'an Jiaotong University，2002.

［6］許社教，虞柏青.計算機繪圖［M］.北京:電子工業(yè)出版社，2003.XU Shejiao，YU Baiqing.Computer Graphics［M］.Beijing:Publishing House of Electronics Industry，2003.

［7］孫祥.Matlab7.0基礎(chǔ)教程［M］.北京:清華大學出版社，2005.SUN Xiang.Matlab7.0 Guide Basis［M］.Beijing:Tsinghua University Press，2005.

［8］熊章強.淺層地震勘探［M］.北京:地震出版社，2002.XIONG Zhangqiang.Shallow Seismic Exploration［M］.Beijing:Seismological Press，2002.

［9］王家映.地球物理反演理論［M］.北京:高等教育出版社，2001.WANG Jiaying.Geophysical Inversion Theory［M］.Beijing:Higher Education Press，2001.

［10］王超凡，鄒桂高，劉金光，等.多道瞬態(tài)瑞雷波的勘探應(yīng)用研究［J］.地質(zhì)科學，2002，37(1):110－117.WANG Chaofan，ZOU Guigao，LIU Jinguang，et al.Application of Multi－Channel Transient Rayleigh Wave Exploration Study［J］.Geological Sciences，2002，37(1):110－117.

［11］杜立志.瞬態(tài)瑞雷波勘探中的數(shù)字處理技術(shù)研究［D］.長春:吉林大學建設(shè)工程學院，2005.DU Lizhi.The Study of Digital Processing Technlogy in Transient Rayleigh Wave Exploration［D］.Changchun:College of Construction Engineering，Jilin University，2005.

［12］于潤偉，朱曉慧.Matlab基礎(chǔ)及應(yīng)用［M］.2版.北京:機械工業(yè)出版社，2008.YU Runwei，ZHU Xiaohui.Foundation and Application of Matlab［M］.2nd ed.Beijing:China Machine Press，2008.