基于短時(shí)互相關(guān)的地表橫波波速研究

靳伯驁，高艷

（1.中國(guó)科學(xué)院噪聲與振動(dòng)重點(diǎn)實(shí)驗(yàn)室，北京 100190；2.中國(guó)科學(xué)院大學(xué)，北京 100049）

1 引言

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

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

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

2 理論

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

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

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

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

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

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

3 地表橫波波速實(shí)驗(yàn)

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

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

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

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

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

圖2 兩個(gè)三軸加速度計(jì)的脈沖與背景噪聲片段的頻譜

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

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

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

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

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

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

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

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

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

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

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

類似于短時(shí)傅里葉變換，短時(shí)互相關(guān)算法使用了時(shí)頻域局域化的時(shí)間窗函數(shù)對(duì)數(shù)據(jù)進(jìn)行截?cái)?。?shí)際處理時(shí)，首先對(duì)x軸時(shí)域平均波形進(jìn)行降采樣至4096Hz。短時(shí)截?cái)鄷r(shí)，時(shí)間窗函數(shù)設(shè)為256點(diǎn)的哈恩窗（Hann Window）。時(shí)間窗起點(diǎn)從錘擊信號(hào)數(shù)據(jù)片段起點(diǎn)開(kāi)始，步長(zhǎng)為6，計(jì)算每個(gè)時(shí)間窗起點(diǎn)所對(duì)應(yīng)的互相關(guān)函數(shù)，繪制于圖5。

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

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

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

5 結(jié)論

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