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

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

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

關鍵詞：河道整治工程;根石;聲波;時頻分析;頻譜響應;黃河

中圖分類號：TV871.2;TV882.1 文獻標志碼：A

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

1 頻率敏感性統計

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

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

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

信號的有效信息不僅可以在原始時間域中提取，也可以轉換到其他合適的域中提取。當前能從3個方面獲取蘊含在信號中的信息：時間域、頻率域、時頻分布。根石厚度的問題類似于一定傾角的超薄地層問題，考慮河水、泥沙、石塊的反射吸收系數不同，獲取的信號資料有不同的頻率響應，將它們對應的頻率響應展現在時間維度上，獲得其頂、底界面的縱向延伸，就可以劃分出它們之間的界面?？梢哉f，信號頻率在時間上的分布顯得至關重要，區分根石的頂、底界面的反射波信號，獲取根石對應頻率分布的時間區間，就可以計算出根石的厚度。

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

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

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

短時傅里葉變換基本思想是將信號加滑動時間窗，也就是傅里葉變換乘上一個有限時間單元的窗函數，對時間窗內信號作傅里葉變換，將窗函數沿時間軸按照一定的步距移動，從而實現信號的逐段分析，得到信號的時變頻譜。短時傅里葉變換是目前研究非平穩信號有力的工具，在各個領域中應用廣泛。

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

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

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

3 時頻分布特征

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

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

4 根石厚度估算

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

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

5 結論

（1）信號的頻率響應在泥沙反射段和根石段之間的能量分布有一定差別，泥沙段優勢頻率為1300Hz，根石段優勢頻率為4500Hz，根石段的優勢頻率主要是高頻成分。

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

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

參考文獻：

[1]張賢達，保錚.非平穩信號分析與處理[M].北京：國防工業出版社，1998：35-48.

[2]李世雄，陳東方.信號瞬時參數計算方法評價[J].信號處理，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]趙淑紅，張文波.短時傅立葉變換在研究沉積旋回地質體中的應用[J].長安大學學報（地球科學版），2003，25（2）：59-62.

[5]劉麗娟.時頻分析技術及其應用[D].成都;成都理工大學，2008：20-25.

[6]何旭，李鴻光.利用經驗模式分解提高短時傅立葉變換分辨率[J].廣西師范大學學報（自然科學版），2005，23（1）：5-8.

[7]葛哲學，陳仲生.Matlab時頻分析技術及其應用[M].北京：人民郵電出版社，2006：36-40.

[8]周開明，劉賢紅，查樹貴，等.幾種時頻分析技術的性能研究及在河道砂體預測中的應用[J].天然氣工業，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]劉財，張海江.小波變換及其在薄儲層識別中的應用[J].石油物探，1995，34（3）：23-31.

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

[14]蔣龍聰.薄互層儲層地震響應時頻分析研究[D].武漢：中國地質大學，2008：55-60.