靳伯驁，高艷

（1.中國科學院噪聲與振動重點實驗室，北京 100190；2.中國科學院大學，北京 100049）

1 引言

隨著城市的快速發(fā)展，城市地下空間的開發(fā)與利用受到了越來越多的關(guān)注。城市管線定位以及地鐵軌道噪聲相關(guān)研究中使用了彈性波理論，通過實驗或仿真手段，對聲源或反射體進行定位[1～3]。其中，地表波速是一個關(guān)鍵參數(shù)，其精度直接影響定位的精度[4]。同時，工程中也可利用縱波與橫波波速確定土壤的等效體積模量以及剪切模量等相關(guān)參數(shù)。但由于城市環(huán)境和地下土況十分復雜，測試結(jié)果受周圍環(huán)境和多模式聲波影響很大。實際的波速計算中需要魯棒性較強的算法?，F(xiàn)有的計算方法有傳遞函數(shù)相位法和互相關(guān)方法，通過對兩組一維時間信號進行處理，得到兩個信號之間的時間延遲。

傳遞函數(shù)相位法和傳統(tǒng)互相關(guān)方法適用于平穩(wěn)信號的處理，對于時變或突發(fā)信號而言，需將一維時域信號映射到時頻域進行處理。關(guān)于時頻域的算法要追溯到20世紀60年代，語音編解碼技術(shù)中已開始使用短時傅里葉變換[5、6]。隨著計算技術(shù)的發(fā)展，時頻域分析技術(shù)也逐漸演進[7、8]，先進的小波技術(shù)已應用于提取和計算湍流相干結(jié)構(gòu)的研究[9]。對于互相關(guān)算法，在快速傅里葉變換之前已被大量使用，是一種十分常見的計算操作，目前已被廣泛應用于計算時間延遲相關(guān)的研究中[10]。

本文提出利用短時加窗的互相關(guān)方法求取淺表橫波波速。使用兩個傳感器進行水平橫波的實驗，分別通過傳遞函數(shù)梯度法、傳統(tǒng)互相關(guān)法以及短時互相關(guān)方法計算地震波波速。實驗發(fā)現(xiàn)地震波信號能量集中在55～250Hz內(nèi)；傳遞函數(shù)梯度法求取的波速隨頻帶選取變化；短時互相關(guān)方法比傳統(tǒng)互相關(guān)方法的結(jié)果更接近傳遞函數(shù)梯度法中線性區(qū)間段的結(jié)果。

2 理論

現(xiàn)有的頻率、時頻域或尺度變換都是將時域信號通過一系列基函數(shù)變換到所對應的域中。例如，對于短時傅里葉變換有：

式①中：i為虛數(shù)單位，ω為角頻率，t和τ分別代表時間和時間延遲，是時域信號，是時間窗函數(shù)，是短時傅里葉變換結(jié)果。由式①可發(fā)現(xiàn)，時頻域變換通過兩個算子和，可將一維的時域數(shù)據(jù)轉(zhuǎn)到二維的時頻域中。類似地，在小波變換中，可通過變化兩個參數(shù)來計算時間尺度域的結(jié)果。

短時互相關(guān)算法與短時傅里葉變換和小波變換，有很多相似的地方。假設測量的兩路信號為和，仿照短時傅里葉變換以及小波變換，短時互相關(guān)的計算表達式為：

式③中：δ代表時間延遲，τ代表時間窗起點。短時互相關(guān)通過移動時間窗函數(shù)截取信號片段與進行互相關(guān)，旨在對空間中多波速、多波數(shù)信號進行匹配，求解信號的時間延遲?；谝陨纤惴ㄟM行數(shù)值仿真。仿真中假設信號傳遞路徑是線性的，兩路信號之間相差一個時間延遲，同時信噪比較高，忽略了噪聲影響。仿真中的采樣頻率設為4096Hz，與分別使用了中心頻率為200Hz的雷克子波（又稱草帽函數(shù)），使用了32點的Hann窗作為短時窗函數(shù)，兩個信號之間的延遲為38個數(shù)據(jù)點，對應的時間延遲為9.033ms。短時互相關(guān)的結(jié)果如圖1所示。

圖1 數(shù)值仿真中兩個雷克子波的短時互相關(guān)結(jié)果

圖1的灰度圖中顯示的是短時互相關(guān)結(jié)果，橫軸是時間窗函數(shù)的起點時刻，縱軸對應的是互相關(guān)的時間延遲，灰度表示的是相應坐標下的互相關(guān)值；圖1也繪制了時域的波形圖，經(jīng)過放縮與平移后繪在灰度圖之上，其時間軸與灰度圖的時間窗起點軸一致。由圖1可知，兩個雷克子波的短時互相關(guān)的峰值在實線峰值之前，表示該點的時間窗正好包含了實線信號的大部分能量，從而導致互相關(guān)得到最大峰值。時間延遲的坐標點與生成信號時輸入的時間延遲一致。仿真結(jié)果表明，理想情況下短時互相關(guān)算法能夠準確確定兩個雷克子波之間的時間延遲。

3 地表橫波波速實驗

實驗選址在中科院聲學研究所院內(nèi)的綠化帶中，測試現(xiàn)場與周圍建筑物、柏油路面以及高大喬木的距離大于3m，地表附有稀疏草本植被。在場地中心安裝耙狀聲源底座，底座激勵的方向平行于南北經(jīng)線方向。在底座西側(cè)沿東西緯線方向布放兩個三軸加速度計（LC0161-2），加速度計、底座三點共線，間距1m。加速度計通過螺栓固定于特制金屬底座上，底座有三個長為5cm的齒狀物埋于地下，加速度計與地表緊密連接。加速度計的x軸、y軸、z軸分別平行于經(jīng)線、緯線和鉛垂線。測試時，為避免地表負載的影響，負責錘擊操作的實驗員站在底座上，數(shù)據(jù)采集設備通過信號線在5m外的柏油路面上進行，采樣頻率設為32 768Hz。由于手動錘擊操作容易發(fā)生雙擊，并且受到現(xiàn)場多因素影響，一些文獻中提出多次激勵，對統(tǒng)計平均后的數(shù)據(jù)進行處理。鑒此，先后進行了10次錘擊。

對采集到的時域信號進行初步處理，設置觸發(fā)閾值為3m/s2，對1m處的加速度計的x軸數(shù)據(jù)進行逐點掃描。當滿足觸發(fā)條件時，判斷此時刻為錘擊起始點。截取每一個錘擊起始點前3000點至起始點后12 000點的數(shù)據(jù)。通過數(shù)據(jù)的預處理剔除了存在明顯差異的2組數(shù)據(jù)。

4 數(shù)據(jù)處理與分析

大多數(shù)的地表波速測量方法都是使用已知間距的傳感器陣列，通過測量傳感器之間的地震波到達時間差進行波速估計。通過對實驗數(shù)據(jù)的預處理得到了8組錘擊實驗的數(shù)據(jù)，每組中包含兩個三軸加速度計，共計六個通道的數(shù)據(jù)，每通道都有15 000個點，持續(xù)時間接近0.5s。由離散傅里葉變換特性可知，這些數(shù)據(jù)的頻域分辨率接近2Hz。截取錘擊脈沖片段的同時，也截取了每段脈沖之后相同長度的背景噪聲片段。

首先，對每一段數(shù)據(jù)進行傅里葉變換，截取1kHz以內(nèi)的頻譜，幅值取對數(shù)后繪制于圖2。圖2中繪制了六個通道脈沖以及噪聲的頻域數(shù)據(jù)。其中（a）、（b）、（c）三個子圖分對應近端傳感器x軸、y軸、z軸三個通道，（d）、（e）、（f）對應遠端傳感器的三個通道；每個子圖中都繪制了8次錘擊實驗的結(jié)果，其中實線是脈沖實驗的數(shù)據(jù)，點圖是脈沖后靜息時的背景噪聲。對比六個通道的數(shù)據(jù)發(fā)現(xiàn)，土壤中的地震波信號在55Hz至250Hz的頻段內(nèi)信噪較高。本文使用半透明灰色矩形窗對該頻段進行了標記。

圖2 兩個三軸加速度計的脈沖與背景噪聲片段的頻譜

4.1 傳遞函數(shù)相位法的結(jié)果

由頻域分析發(fā)現(xiàn)，8組錘擊實驗的一致性較好。由于錘擊脈沖峰值幅度接近，且使用同一個判斷閾值，故假設每組錘擊實驗起始點時刻的差異可忽略。對8組錘擊實驗結(jié)果求取平均，得到平均的錘擊實驗的時域波形。計算平均波形中x軸信號的傳遞函數(shù)相位，將其繪制于圖3。

圖3 平均時域信號的x軸互相關(guān)相位圖以及三個頻率區(qū)間段內(nèi)計算得到的波速

從圖3可發(fā)現(xiàn)，傳遞函數(shù)相位在55Hz至170Hz頻段內(nèi)呈線性變化。對照圖2發(fā)現(xiàn)，該頻段與地震波峰值頻段對應。此外，在55Hz至250Hz內(nèi)，傳遞函數(shù)相位梯度存在變化。實驗中兩個三軸加速度計的間距為1m，依此計算了三個頻率區(qū)間段的等效波速，分別約為116.9m/s、82.2m/s、46.0m/s。初步分析猜測，波速變化的一種原因可能是由于地表土壤分層情況復雜，在200Hz附近存在水平橫波（SH波）模式之外的波，造成波速結(jié)果偏移；另一種可能是埋地情況復雜，傳遞路徑上存在中心頻率接近200Hz的共振結(jié)構(gòu)，導致兩個傳感器之間的相位發(fā)生偏移。

通過對錘擊信號進行頻譜分析，從圖2可發(fā)現(xiàn)，x軸的信號的能量高于其余兩軸?？紤]到加速度計安裝時存在傾斜誤差，初步判定實驗中SH波模式的能量占主導。同時，圖3的線性頻段與圖2的峰值頻段較為吻合，均在55Hz至170Hz之間。由上述分析初步判定，實驗中的水平錘擊主要激發(fā)出了SH波，能量集中在55Hz至170Hz之間，且波速約為116.9m/s。

4.2 傳統(tǒng)互相關(guān)的結(jié)果

傳統(tǒng)互相關(guān)方法通常用于描述兩個時間序列之間的相關(guān)程度。在傳遞路徑分析中，當激勵信號為寬帶脈沖或白噪聲時，互相關(guān)峰值所對應的時間延遲即為信號到達兩個傳感器的時間差。

本文計算了三個軸向的8次錘擊實驗的互相關(guān)函數(shù)，同時計算了x軸時域平均波形的互相關(guān)，繪制于圖4。由于每次錘擊的力度不同，同一個軸向的結(jié)果不同錘擊的幅度偏差較明顯，不易于觀察，因此圖4對x軸、y軸、z軸數(shù)據(jù)分別進行了能量歸一化處理。

圖4 x軸、y軸、z三軸的傳統(tǒng)互相關(guān)結(jié)果以及x軸時域平均波形的互相關(guān)結(jié)果

由圖4可發(fā)現(xiàn)，x軸互相關(guān)峰值遠高于z軸互相關(guān)峰值，且x軸與z軸互相關(guān)零點、峰值點位置十分一致。初步判定傳感器安裝時存在俯仰角方向的偏差，z軸加速度計采集到了水平方向的分量。對于壓縮波、Rayleigh波以及豎直橫波（SV波）模式而言，均存在豎直方向的運動分量，而實驗結(jié)果中豎直振動十分微弱，因此可推測實際測試中主要激發(fā)出了SH模式的地震波，與頻譜分析的結(jié)論一致。圖4實線為x軸時域平均波形的互相關(guān)，通過時間延遲以及已知的傳感器間距計算得到地震波波速約為110.0m/s，與傳遞函數(shù)相位法得到的波速相差接近6%。

4.3 短時互相關(guān)的結(jié)果

類似于短時傅里葉變換，短時互相關(guān)算法使用了時頻域局域化的時間窗函數(shù)對數(shù)據(jù)進行截斷。實際處理時，首先對x軸時域平均波形進行降采樣至4096Hz。短時截斷時，時間窗函數(shù)設為256點的哈恩窗（Hann Window）。時間窗起點從錘擊信號數(shù)據(jù)片段起點開始，步長為6，計算每個時間窗起點所對應的互相關(guān)函數(shù)，繪制于圖5。

圖5 x軸時域平均波形的短時互相關(guān)結(jié)果

圖5的灰度圖中，每一條平行于y軸的灰度數(shù)據(jù)對應一個時刻的互相關(guān)，縱軸坐標即為互相關(guān)的時間延遲。得到短時互相關(guān)的灰度圖后，又計算了每一個時間窗起點下的互相關(guān)峰值所對應的時間延遲，用黑色圓點繪制于圖5中。可發(fā)現(xiàn)，在時間窗起點為0.04s至0.1s的區(qū)間段內(nèi)，互相關(guān)峰值較為明顯，峰值坐標較為一致。這是因為在這個區(qū)間段內(nèi)，時間窗能將錘擊信號的大部分能量包含進去，使得這個區(qū)間內(nèi)的互相關(guān)峰值較高。求取所有區(qū)間段的峰值中的最大值，用黑色叉號標記于圖5中。記錄該最大值的時間延遲，得到短時互相關(guān)方法下的地震波波速約為113.8m/s。該結(jié)果介于傳遞函數(shù)相位法和傳統(tǒng)互相關(guān)方法的結(jié)果。

對比三種方法求取互相關(guān)的過程發(fā)現(xiàn)，傳遞函數(shù)相位法受到地表復雜性以及所選取頻率區(qū)間的影響，其結(jié)果存疑；傳統(tǒng)互相關(guān)方法中，8次錘擊實驗的結(jié)果一致性較差，x軸時域平均波形的結(jié)果與時域數(shù)據(jù)的結(jié)果存在不一致，有待進一步研究；短時互相關(guān)方法物理圖景較為清晰，能夠分析多模式、時變信號，但對于其窗函數(shù)選取以及算法魯棒性等問題仍有待進一步研究。

5 結(jié)論

通過水平錘擊聲源在地表激發(fā)水平橫波，使用兩個三軸加速度計對地表振動進行測量。通過數(shù)據(jù)分析初步判定，該方法所產(chǎn)生的地震波主要以SH波為主。通過傳遞函數(shù)相位方法、傳統(tǒng)互相關(guān)法以及短時互相關(guān)方法計算地震波波速，分別約為116.9m/s、110.0m/s及113.8m/s。對結(jié)果進行分析發(fā)現(xiàn)，實驗現(xiàn)場的地表水平橫波波速在110m/s至120m/s附近，但更高精度的結(jié)果仍有待進一步研究。同時，短時互相關(guān)方法擁有分析多模式、時變信號的能力，物理含義清晰，存在進一步研究的價值。