基于時(shí)頻分析法研究黃河河道壩垛根石厚度

張宏杰 郭乙霏 李廣超 侯超普 陳萌

摘要：黃河中下游河道整治工程壩垛水下的散拋石塊體大小不一，散亂堆砌成有坡度的堆石薄層，采用非接觸式探All方法探測(cè)時(shí)繞射波比較強(qiáng)烈。根石探測(cè)數(shù)據(jù)主要在時(shí)間域處理，僅能人工拾取得到根石頂界面，很難識(shí)別根石底界面，無法在時(shí)間域獲取根石厚度。通過對(duì)疊后時(shí)間偏移的探測(cè)數(shù)據(jù)作頻譜分析與s變換，獲取河床泥沙段與根石段的頻譜特征;根據(jù)河床泥沙與根石段的時(shí)頻特征差異，設(shè)計(jì)識(shí)別根石的算法;利用時(shí)頻分布圖中頻率在時(shí)間維度的延續(xù)計(jì)算根石的厚度，實(shí)現(xiàn)自動(dòng)化識(shí)別，避免了人工識(shí)別波形造成的誤差。

關(guān)鍵詞：河道整治工程;根石;聲波;時(shí)頻分析;頻譜響應(yīng);黃河

中圖分類號(hào)：TV871.2;TV882.1 文獻(xiàn)標(biāo)志碼：A

黃河中下游河道整治工程以壩垛調(diào)控河流走勢(shì)，受河水沖刷影響，經(jīng)常發(fā)生根石走失現(xiàn)象。探測(cè)清楚根石的缺失狀況，在汛后的加固處理中就會(huì)有針對(duì)性，避免治理工作的盲目性。非接觸式根石探測(cè)方法主要是采用淺地層剖面儀獲取根石的聲波反射資料，并加以分析處理。泥沙、根石、河水之間存在較大的波阻抗差異，存在明顯的反射界面，在資料處理過程中，人工能夠較清晰地識(shí)別根石頂界面，但根石厚度較小，堆砌有一定的坡度，且塊石大小不均勻、堆砌散亂，在淺地層剖面儀獲取的記錄上形成較強(qiáng)的繞射，根石界面的反射波和繞射波及各種干擾波混疊在一起，沒有辦法在時(shí)間域獲取根石厚度和壩前的最大沖刷深度。通過對(duì)疊后時(shí)間偏移的淺地層剖面數(shù)據(jù)作頻譜分析與S變換，可獲取河床泥沙與根石段的頻譜特征。筆者根據(jù)河床泥沙與根石段的時(shí)頻特征差異，設(shè)計(jì)識(shí)別根石的算法;利用時(shí)頻分布圖中頻率在時(shí)間維度的延續(xù)計(jì)算根石的厚度，以實(shí)現(xiàn)自動(dòng)化識(shí)別。

1 頻率敏感性統(tǒng)計(jì)

原始記錄經(jīng)過疊后偏移處理后可以看出，1～80道的記錄為泥沙段反射波，平滑連續(xù)，能量集中，并且在泥沙界面之下沒有其他反射波的存在，只有一些多次波的存在;81～150道記錄為根石段反射波，有較強(qiáng)的繞射波，有一定的延續(xù)度，見圖1。

從時(shí)間記錄很難識(shí)別出根石的厚度，需要將聲波道集進(jìn)行時(shí)頻分析，并將時(shí)頻分析后的數(shù)據(jù)按照頻率進(jìn)行振幅統(tǒng)計(jì)，并做歸一化處理。圖2為所有150道數(shù)據(jù)的時(shí)頻分析頻率統(tǒng)計(jì)，可以看出統(tǒng)計(jì)圖上有兩個(gè)峰值，其中的低頻段在1300Hz附近，整體能量較強(qiáng)，均高于高頻段的峰值;高頻段的峰值在4000Hz附近，整體能量分布不太集中。

根據(jù)波形圖中反映的根石和泥沙界面段，拾取相對(duì)應(yīng)的兩個(gè)時(shí)間窗口，將這兩個(gè)時(shí)窗段內(nèi)的頻率譜分別進(jìn)行統(tǒng)計(jì)。圖3為泥沙段的頻率統(tǒng)計(jì)，可以看出和圖2之間有著很明顯的差異。圖3中低頻能量過強(qiáng)，為了便于比較做了限幅，在高頻段有非常微弱的峰值跡象。圖4為根石段的頻率統(tǒng)計(jì)，有兩個(gè)比較明顯的峰值，兩個(gè)峰值對(duì)應(yīng)頻率分別為1300Hz和4500Hz。這和圖2是一致的，因此可以從頻率域來分析泥沙和根石的界面，根據(jù)高頻成分的延續(xù)來獲取根石的厚度。2頻率域分析處理

信號(hào)的有效信息不僅可以在原始時(shí)間域中提取，也可以轉(zhuǎn)換到其他合適的域中提取。當(dāng)前能從3個(gè)方面獲取蘊(yùn)含在信號(hào)中的信息：時(shí)間域、頻率域、時(shí)頻分布。根石厚度的問題類似于一定傾角的超薄地層問題，考慮河水、泥沙、石塊的反射吸收系數(shù)不同，獲取的信號(hào)資料有不同的頻率響應(yīng)，將它們對(duì)應(yīng)的頻率響應(yīng)展現(xiàn)在時(shí)間維度上，獲得其頂、底界面的縱向延伸，就可以劃分出它們之間的界面。可以說，信號(hào)頻率在時(shí)間上的分布顯得至關(guān)重要，區(qū)分根石的頂、底界面的反射波信號(hào)，獲取根石對(duì)應(yīng)頻率分布的時(shí)間區(qū)間，就可以計(jì)算出根石的厚度。

傅里葉變換是信號(hào)處理中最常用的方法，該變換是在時(shí)間域和頻率域的整個(gè)空間展現(xiàn)信號(hào)的，是全局的整體變換，可得到信號(hào)的頻率分布情況，但其時(shí)間分辨率為零[1]。還有常用的Hilbert變換方法Teager -Kaiser方法、Shekel方法，屬于非時(shí)頻方法的瞬時(shí)頻率估計(jì)[2]，得到的是單值函數(shù)，能處理任何時(shí)刻單一頻率的信號(hào)。淺地層剖面儀發(fā)射頻帶范圍為0.5～7.0kHz，經(jīng)過水下介質(zhì)濾波作用后被探頭接收的信號(hào)包含多種頻率成分，用這些非時(shí)頻方法得到的瞬時(shí)頻率和實(shí)際信號(hào)的自身頻率完全不能對(duì)應(yīng)，也就失去了其大部分的應(yīng)用價(jià)值。

時(shí)頻分析技術(shù)能夠表征信號(hào)的頻率隨時(shí)間的變化，其核心思想是構(gòu)建一個(gè)時(shí)間和頻率的密度函數(shù)，將每一道的一維時(shí)間信號(hào)映射到二維的時(shí)間一頻率平面，展示信號(hào)頻率是如何隨時(shí)間變化的，不僅解決了傅里葉變換沒有時(shí)間分辨率的問題，也克服了Hilbert變換等方法只能處理單一頻率的不足。利用時(shí)頻分析結(jié)果能夠做出時(shí)頻分布圖形（二維或三維），指示每一時(shí)刻的信號(hào)在瞬時(shí)頻率附近的能量聚集情況，等效的時(shí)域信號(hào)可以通過時(shí)頻表示的求逆過程得到[3]。趙淑紅等[4]利用短時(shí)傅里葉變換的時(shí)頻分析研究沉積旋回，取得了很好的效果。

時(shí)頻分析的線性時(shí)頻表示主要有短時(shí)傅里葉變換、小波變換、S變換等，二次型時(shí)頻表示反映信號(hào)能量分布的有Wigner-Ville分布等。

短時(shí)傅里葉變換基本思想是將信號(hào)加滑動(dòng)時(shí)間窗，也就是傅里葉變換乘上一個(gè)有限時(shí)間單元的窗函數(shù)，對(duì)時(shí)間窗內(nèi)信號(hào)作傅里葉變換，將窗函數(shù)沿時(shí)間軸按照一定的步距移動(dòng)，從而實(shí)現(xiàn)信號(hào)的逐段分析，得到信號(hào)的時(shí)變頻譜。短時(shí)傅里葉變換是目前研究非平穩(wěn)信號(hào)有力的工具，在各個(gè)領(lǐng)域中應(yīng)用廣泛。

短時(shí)傅里葉變換（STET）定義為[5]式中：f為頻率;t為時(shí)間;h（τ）為信號(hào);w（τ-t）為窗函數(shù)。

根據(jù)Heisenberg測(cè)不準(zhǔn)原理[6]，一旦窗函數(shù)選定，時(shí)頻分辨率便確定下來，在窗口形狀不變的情況下，窗口面積越小，時(shí)頻表征局部化的能力就越強(qiáng)。要提高時(shí)間分辨率，就要選擇的窗函數(shù)盡可能短;要提高頻率分辨率，則要求窗函數(shù)的時(shí)間寬度盡可能長(zhǎng)。短時(shí)傅里葉變換雖然能描述時(shí)頻的變化情況，但它在時(shí)域、頻域的分辨率Δt、Δf不隨時(shí)間t和頻率f的變化而變化，不能敏感地反映信號(hào)的突變，適用于對(duì)緩變信號(hào)的分析，對(duì)突變信號(hào)和非平穩(wěn)信號(hào)的分析存在局限性。對(duì)于采集的時(shí)間序列信號(hào)，很難找到一個(gè)合適的時(shí)間窗函數(shù)來適應(yīng)不同的時(shí)段[7]。

S變換是介于短時(shí)傅里葉變換（STFT）和小波變換（WT）之間的一種非平穩(wěn)信號(hào)分析和處理的方法，綜合了短時(shí)傅里葉變換和小波變換的優(yōu)點(diǎn)[8]。它的時(shí)間一頻率譜分辨率與頻率相關(guān)，不但有多尺度聚焦性，而且直接與Fourier譜聯(lián)系，保持頻率的絕對(duì)相位，基本變換函數(shù)不必滿足容許性條件[9-13]。信號(hào)h（t）的S變換表示為[14]式中：w（t）為窗函數(shù);σ為頻率的函數(shù)。

3 時(shí)頻分布特征

分析由淺地層剖面儀獲取的資料，抽取其中包含河床泥沙段和根石段的150道記錄進(jìn)行分析，單道時(shí)間域采樣點(diǎn)數(shù)1446，采樣頻率21739Hz，采樣時(shí)間0.046ms（見圖1）。利用幾種時(shí)頻分析方法，分別對(duì)每一道記錄進(jìn)行時(shí)頻分析，獲取單道數(shù)據(jù)的時(shí)頻分布，對(duì)比總結(jié)出反射波在泥沙和根石界面的時(shí)頻分布特征：河床泥沙段的時(shí)頻分布特征是頻率能量具有較小的時(shí)間延續(xù)度，并且能量較強(qiáng)，強(qiáng)度超過180;根石段的時(shí)頻分布特征是在頻率1000～6500Hz范圍內(nèi)有較大的時(shí)間延續(xù)度，能量比較分散，能量強(qiáng)度小于100。圖5為河床泥沙段（CDP25、CDP45）和根石段（CDP100、CDP130）單道記錄的S變換時(shí)頻分布。

圖6為S變換的時(shí)頻分布，分析的時(shí)間點(diǎn)數(shù)為480。可以看出1～85道泥沙段的頻譜具有明顯的特征：能量強(qiáng)，連續(xù)圓滑，一致性好。86～140道的根石段頻譜能量特點(diǎn)是能量比較強(qiáng)，能夠看出根石的輪廓，有比較圓滑清晰的下邊界，也就是根石的底部反射，利用這些信息可以確定根石的厚度。

4 根石厚度估算

根據(jù)泥沙段和根石段的時(shí)頻分布特征，反射記錄經(jīng)過S變換，按照?qǐng)D1記錄道集方式展示時(shí)頻分布圖，利用式（5）計(jì)算時(shí)頻分布屬性，并設(shè)定：任意點(diǎn)數(shù)據(jù)的反射能量強(qiáng)度超過120后，將其能量強(qiáng)度降低到20，其目的是突出根石反射的能量強(qiáng)度，便于利用時(shí)頻分布圖自動(dòng)識(shí)別根石（見圖7）。式中：A（i，j）為頻率分布的能量強(qiáng)度;tf（i，j，k）為第i道記錄第j個(gè)采樣點(diǎn)、第K個(gè)頻點(diǎn)的時(shí)頻分布。

圖7中較深色塊表示根石反射范圍，其邊界是明顯的，可以方便地估算水的深度和根石的厚度。剖面中任何位置的根石厚度和水的深度都可以根據(jù)圖7中深顏色區(qū)域的頻率分布區(qū)間對(duì)應(yīng)的時(shí)間值估算。例如，在位置CDP120，考慮到傳播時(shí)間聲反射雙程旅行時(shí)間，該處水深1.875m或雙程時(shí)2.5ms，根石厚度5.5m或雙程時(shí)5ms，根石中的平均傳播速度設(shè)定為2200m/s。計(jì)算每一處的根石厚度都需要聲波在該處根石中傳播的速度，這個(gè)速度受到根石塊大小、堆放情況、充填物情況的影響，沒有辦法獲取精確的速度，需要采用平均速度，平均速度的獲取需要一定的經(jīng)驗(yàn)積累和驗(yàn)證。

5 結(jié)論

（1）信號(hào)的頻率響應(yīng)在泥沙反射段和根石段之間的能量分布有一定差別，泥沙段優(yōu)勢(shì)頻率為1300Hz，根石段優(yōu)勢(shì)頻率為4500Hz，根石段的優(yōu)勢(shì)頻率主要是高頻成分。

（2）河床泥沙段的時(shí)頻分布能量比較集中，光滑連續(xù)，在時(shí)間維度的延續(xù)比較短;根石段的時(shí)頻分布能量相對(duì)集中，上下界面清晰但不規(guī)則、不光滑，時(shí)頻分布在時(shí)間維的延伸較大。

（3）S變換比較適用于傾斜薄層信號(hào)處理，獲取了根石段信號(hào)在時(shí)間維度的延續(xù)，就能夠求取散拋的薄層根石厚度。

參考文獻(xiàn)：

[1]張賢達(dá)，保錚.非平穩(wěn)信號(hào)分析與處理[M].北京：國防工業(yè)出版社，1998：35-48.

[2]李世雄，陳東方.信號(hào)瞬時(shí)參數(shù)計(jì)算方法評(píng)價(jià)[J].信號(hào)處理，2003，19（1）：59-63.

[3]ROESSGEN M，BOASHASH B.Time-Frequency Peak Filte-ring Applied to FSK Signals[C]//Ieee-Sp International Sym-posium on Time-Frequency and Time-Scale Analysis.IEEE，1994：516-519.

[4]趙淑紅，張文波.短時(shí)傅立葉變換在研究沉積旋回地質(zhì)體中的應(yīng)用[J].長(zhǎng)安大學(xué)學(xué)報(bào)（地球科學(xué)版），2003，25（2）：59-62.

[5]劉麗娟.時(shí)頻分析技術(shù)及其應(yīng)用[D].成都;成都理工大學(xué)，2008：20-25.

[6]何旭，李鴻光.利用經(jīng)驗(yàn)?zāi)Ｊ椒纸馓岣叨虝r(shí)傅立葉變換分辨率[J].廣西師范大學(xué)學(xué)報(bào)（自然科學(xué)版），2005，23（1）：5-8.

[7]葛哲學(xué)，陳仲生.Matlab時(shí)頻分析技術(shù)及其應(yīng)用[M].北京：人民郵電出版社，2006：36-40.

[8]周開明，劉賢紅，查樹貴，等.幾種時(shí)頻分析技術(shù)的性能研究及在河道砂體預(yù)測(cè)中的應(yīng)用[J].天然氣工業(yè)，2007（s1）：451-453.

[9]COHEN L.Time-Frequency Analysis：Theory and Application[M].Englewood Cliffs，NJ;Prentice Hal，Inc，1995：20-25.

[10]GARBOR D.Theory of Communication[J].The Journal of theInstitute of Electrical Engineers，1946，93（3）：429-456.

[11]劉喜武，年靜波，黃文松.利用廣義S變換提取地震旋回的方法[J].石油物探，2006，45（2）：129-134.

[12]劉財(cái)，張海江.小波變換及其在薄儲(chǔ)層識(shí)別中的應(yīng)用[J].石油物探，1995，34（3）：23-31.

[13]孫魯平，鄭曉東，首皓，等.薄層地震峰值頻率與厚度關(guān)系研究[J].石油地球物理勘探，2010，45（2）;254-259.

[14]蔣龍聰.薄互層儲(chǔ)層地震響應(yīng)時(shí)頻分析研究[D].武漢：中國地質(zhì)大學(xué)，2008：55-60.