簡(jiǎn)介本資源是一套面向地震學(xué)研究者與地球物理方向初學(xué)者的MATLAB頻譜分析實(shí)踐工具包聚焦快速傅里葉變換FFT在地震波形處理中的核心應(yīng)用解決地震時(shí)間序列到頻率域轉(zhuǎn)換、頻譜可視化及特征識(shí)別等關(guān)鍵問(wèn)題。壓縮包共含4個(gè)文件2個(gè).asv備份腳本、1個(gè).m主程序、1個(gè).fig圖形結(jié)果總大小僅11KB輕量實(shí)用其中.m文件實(shí)現(xiàn)完整流程地震數(shù)據(jù)讀取、采樣率估算、FFT計(jì)算、正頻率截取、幅度譜繪制.asv文件保留調(diào)試過(guò)程便于理解代碼演進(jìn)邏輯.fig直觀呈現(xiàn)頻譜分布。已有306人學(xué)習(xí)下載適合課程實(shí)驗(yàn)、科研入門(mén)或項(xiàng)目快速?gòu)?fù)現(xiàn)。用戶(hù)可直接運(yùn)行主程序獲得可復(fù)用的地震頻譜分析框架掌握P波/S波頻段識(shí)別、采樣率適配、幅度譜歸一化等實(shí)操要點(diǎn)并基于現(xiàn)有結(jié)構(gòu)拓展濾波、時(shí)頻分析等進(jìn)階功能。 做地震數(shù)據(jù)處理這行繞不開(kāi)頻域分析。不管是天然地震的震相識(shí)別、工程地震的場(chǎng)地反應(yīng)計(jì)算還是微震監(jiān)測(cè)里的噪聲壓制FFT快速傅里葉變換都是用得最多的基礎(chǔ)工具之一。很多人下載過(guò)各種以“FFT地震”命名的MATLAB腳本包但真正拿到手能跑通、跑對(duì)、跑出能解釋的結(jié)果往往還要踩不少坑。這篇文章就圍繞“地震頻譜分析”這個(gè)主題結(jié)合MATLAB從原理到底層實(shí)現(xiàn)再把實(shí)操中容易翻車(chē)的地方逐條梳理一遍。先說(shuō)說(shuō)這篇文章是給誰(shuí)看的。如果你是剛接觸地震信號(hào)處理的本科生或研究生手里有一段地震波形但不知道怎么轉(zhuǎn)成頻譜這篇文章可以幫你把來(lái)龍去脈理順如果你已經(jīng)跑過(guò)一些現(xiàn)成腳本但發(fā)現(xiàn)出來(lái)的頻譜形狀怪異、幅值對(duì)不上、主頻和預(yù)期不符那這篇文章的避坑部分應(yīng)該能解決你大部分困惑。我會(huì)先在概念層面講清楚為什么做頻譜分析再帶大家走一遍完整的MATLAB實(shí)現(xiàn)流程最后用一個(gè)實(shí)測(cè)風(fēng)格的地震記錄做案例拆解把所有參數(shù)和代碼都擺出來(lái)。1. 地震頻譜分析的核心思路與原理基礎(chǔ)1.1 為什么要做地震頻譜分析地震記錄的原始形態(tài)是時(shí)間域上的振幅波形它記錄了地面運(yùn)動(dòng)隨時(shí)間的快慢變化。但時(shí)間域波形有一個(gè)天然的局限它只能告訴你“什么時(shí)刻震動(dòng)了多大”很難直接回答“這次振動(dòng)的能量集中在哪個(gè)頻率范圍”。而地震學(xué)里很多關(guān)鍵問(wèn)題恰恰需要頻率信息來(lái)回答。比如場(chǎng)地效應(yīng)評(píng)估同一場(chǎng)地震建在軟土上的建筑和建在基巖上的建筑破壞程度差異巨大本質(zhì)就是因?yàn)檐浲翆?duì)特定頻段有放大作用而這個(gè)頻段正是通過(guò)頻譜分析才能確定。再比如震源參數(shù)反演地震矩、應(yīng)力降、拐角頻率這些物理量都是從位移譜的形態(tài)里提取的。還有結(jié)構(gòu)健康監(jiān)測(cè)里橋梁或高層建筑的自振頻率是否發(fā)生了偏移也是通過(guò)對(duì)比環(huán)境振動(dòng)記錄傅里葉譜在不同時(shí)期的變化來(lái)判斷的。一句話總結(jié)地震波形是“信號(hào)”頻譜分析就是把信號(hào)從時(shí)間域投影到頻率域讓我們能看清這個(gè)信號(hào)里每個(gè)頻率成分的能量大小。FFT不是地震學(xué)的專(zhuān)屬工具但它是把地震信號(hào)“解剖”成頻率成分最快速、最標(biāo)準(zhǔn)的手段這也是為什么MATLAB里幾乎每個(gè)處理地震數(shù)據(jù)的工具箱都繞不開(kāi)fft函數(shù)。1.2 FFT與DFT的關(guān)系為什么地震數(shù)據(jù)處理都用FFT傅里葉變換在教科書(shū)上的定義是連續(xù)積分但計(jì)算機(jī)只能處理離散的有限長(zhǎng)序列所以實(shí)際使用的是離散傅里葉變換DFT。DFT的計(jì)算公式是X(k) Σ_{n0}^{N-1} x(n)·e^(-j·2π·kn/N)直接按這個(gè)公式算N個(gè)點(diǎn)需要N2次復(fù)數(shù)乘法當(dāng)N是256點(diǎn)或512點(diǎn)還勉強(qiáng)能接受但當(dāng)N是4096、8192甚至更大時(shí)計(jì)算量就非?？植懒?。FFT是Cooley和Tukey在1965年提出的快速算法它利用旋轉(zhuǎn)因子的周期性和對(duì)稱(chēng)性把計(jì)算量從N2降到N·log?N。當(dāng)N8192時(shí)直接DFT大約需要6700萬(wàn)次乘法而FFT只需要約10萬(wàn)次差距是三個(gè)數(shù)量級(jí)。地震記錄采樣率通常是100Hz、200Hz甚至更高一段60秒的記錄按200Hz采樣就是12000個(gè)點(diǎn)不做FFT的話很多實(shí)時(shí)處理腳本根本跑不完。此外MATLAB的fft底層還做了大量的內(nèi)存訪問(wèn)優(yōu)化對(duì)于多通道數(shù)據(jù)比如三分量地震儀同時(shí)輸出東西、南北、垂直三分量直接調(diào)用fft的矩陣運(yùn)算能力比逐個(gè)通道循環(huán)快得多。1.3 采樣定理、頻率分辨率與奈奎斯特頻率在動(dòng)手寫(xiě)代碼之前有三個(gè)概念必須刻在腦子里它們決定了頻譜圖的橫軸范圍和分析精度。第一個(gè)是奈奎斯特頻率它是信號(hào)在數(shù)字域里能表示的極限頻率等于采樣率的一半。如果采樣率是200Hz那奈奎斯特頻率就是100Hz。任何超過(guò)奈奎斯特頻率的成分都會(huì)被混疊到低頻段偽造出虛假的“鬼影頻率”。所以地震儀在采集前都會(huì)經(jīng)過(guò)抗混疊濾波器這屬于硬件層面的保障。第二個(gè)是頻率分辨率它等于采樣率除以FFT點(diǎn)數(shù)也就是Δf fs / N。這個(gè)公式非常關(guān)鍵它說(shuō)明了時(shí)間和頻率之間是“蹺蹺板”關(guān)系想要分辨出間隔只有0.01Hz的兩個(gè)相鄰頻率峰就需要把FFT點(diǎn)數(shù)撐到fs/0.01那么大對(duì)應(yīng)的時(shí)域信號(hào)長(zhǎng)度也要夠長(zhǎng)。第三個(gè)是FFT點(diǎn)數(shù)與記錄長(zhǎng)度的關(guān)系。很多初學(xué)者以為FFT點(diǎn)數(shù)可以隨意設(shè)置實(shí)際上如果你只是調(diào)用fft(x, N)N大于原始信號(hào)長(zhǎng)度時(shí)MATLAB會(huì)自動(dòng)補(bǔ)零小于時(shí)會(huì)自動(dòng)截?cái)噙@會(huì)帶來(lái)兩個(gè)后果補(bǔ)零可以提高頻譜的“顯示分辨率”讓曲線更平滑但不會(huì)提高真實(shí)的“物理分辨率”兩個(gè)靠得很近的頻率峰仍然分辨不出來(lái)截?cái)鄤t會(huì)丟失有效信號(hào)嚴(yán)重時(shí)導(dǎo)致頻譜嚴(yán)重畸變。后面我會(huì)專(zhuān)門(mén)講這兩者的區(qū)別和正確用法。2. 地震信號(hào)預(yù)處理FFT之前必須做的事2.1 去掉均值與線性趨勢(shì)拿到一段原始的地震記錄第一件事不是做FFT而是預(yù)處理。為什么因?yàn)镕FT的數(shù)學(xué)本質(zhì)是周期延拓它默認(rèn)你截取的這段信號(hào)是周期性重復(fù)的。如果信號(hào)不滿足這個(gè)假設(shè)頻譜就會(huì)產(chǎn)生“泄漏”現(xiàn)象能量從一個(gè)頻率擴(kuò)散到附近的頻率上導(dǎo)致主頻模糊、旁邊出現(xiàn)虛假的旁瓣。最常見(jiàn)的預(yù)處理操作是去均值。地震計(jì)輸出的原始數(shù)據(jù)通常有一個(gè)直流偏置這個(gè)直流分量的頻率是0Hz它的存在會(huì)讓0Hz處出現(xiàn)一個(gè)巨大的尖峰把其他頻段的幅度壓得幾乎看不見(jiàn)。用MATLAB的detrend函數(shù)可以同時(shí)完成去均值和去線性趨勢(shì)% 去均值和線性趨勢(shì) x_detrend detrend(x, constant); % 只去均值 x_detrend detrend(x, linear); % 去均值去線性趨勢(shì)到底是選constant還是linear對(duì)于幾十秒長(zhǎng)度的地震記錄儀器響應(yīng)漂移通常不明顯用constant就夠。但對(duì)于長(zhǎng)周期地脈動(dòng)記錄或是對(duì)原始記錄做了積分處理后線性趨勢(shì)經(jīng)常出現(xiàn)這時(shí)候要用linear。我自己的經(jīng)驗(yàn)是如果不知道選哪個(gè)就兩個(gè)都試試看頻譜的形態(tài)哪個(gè)更干凈、主峰更突出。2.2 濾波與限帶處理地震信號(hào)的頻帶范圍視震源類(lèi)型和傳播路徑而定。遠(yuǎn)震體波的主頻通常在0.01Hz到1Hz之間近震S波可能在1Hz到10Hz而工程微震或地脈動(dòng)的頻率范圍可以到幾十赫茲。在做FFT之前最好根據(jù)你的研究目的先做一個(gè)帶通濾波把無(wú)關(guān)頻段的干擾去掉。濾波要特別注意邊界效應(yīng)。MATLAB自帶的filter函數(shù)是有延遲和邊界震蕩的處理地震數(shù)據(jù)時(shí)更推薦用filtfilt也就是零相移濾波。它會(huì)對(duì)信號(hào)做正向和反向兩次濾波消除相位畸變但代價(jià)是計(jì)算量翻倍以及信號(hào)首尾各自有一小段被“抹平”。實(shí)操中為了減少這種邊界效應(yīng)可以先把信號(hào)延長(zhǎng)一小段再濾波濾波后裁掉延長(zhǎng)的部分。另一個(gè)細(xì)節(jié)是濾波順序應(yīng)該先濾波再去均值還是反過(guò)來(lái)嚴(yán)格來(lái)說(shuō)應(yīng)該先去均值再濾波。如果先濾波濾波器的瞬態(tài)響應(yīng)會(huì)引入新的臨時(shí)偏置而且有些高通濾波器設(shè)計(jì)不夠好的話會(huì)把直流分量重新“振”出來(lái)。穩(wěn)妥的操作順序是原始數(shù)據(jù) → 去均值/去趨勢(shì) → 帶通濾波 → 重新去均值 → 再做FFT。2.3 數(shù)據(jù)截?cái)嗯c窗函數(shù)選擇預(yù)處理做完之后還有一道工序加窗。前面提到FFT默認(rèn)信號(hào)是周期的但實(shí)際截取的地震記錄首尾幾乎不可能完美銜接這就會(huì)造成頻譜泄漏。加窗的作用就是讓信號(hào)在兩端平滑衰減到零強(qiáng)制“偽造”連續(xù)性。地震數(shù)據(jù)處理里最常用的窗函數(shù)是漢寧窗Hanning和漢明窗Hamming兩者的主瓣寬度和旁瓣衰減略有差異。曾經(jīng)有一次我在處理爆破振動(dòng)信號(hào)時(shí)不加窗的時(shí)候主頻怎么都穩(wěn)定不下來(lái)?yè)Q幾種FFT參數(shù)結(jié)果都不一樣。后來(lái)加了一個(gè)Hanning窗主頻立刻穩(wěn)定在某一個(gè)值附近和理論值完全吻合。加窗的本質(zhì)就是用主瓣變寬一點(diǎn)點(diǎn)去換取旁瓣的大幅衰減這是一個(gè)性?xún)r(jià)比極高的取舍。但要注意加窗會(huì)改變信號(hào)的總能量因?yàn)榇昂瘮?shù)在兩端把信號(hào)乘了接近零的系數(shù)。如果要保持幅值譜的物理意義位移振幅、速度振幅等需要對(duì)FFT結(jié)果做幅值恢復(fù)也就是除以窗函數(shù)的均值。MATLAB里可以這樣操作win hanning(N); x_win x(1:N) .* win; X fft(x_win); X X / mean(win); % 幅值恢復(fù)這個(gè)幅值恢復(fù)步驟很容易被忽略很多書(shū)上沒(méi)有強(qiáng)調(diào)但如果有定量分析需求省掉這一步會(huì)導(dǎo)致振幅系統(tǒng)性偏低。3. MATLAB中地震FFT的具體實(shí)現(xiàn)與參數(shù)詳解3.1 fft函數(shù)的基本調(diào)用與輸出含義MATLAB的fft函數(shù)最基本的調(diào)用是X fft(x)但在地震數(shù)據(jù)處理中更規(guī)范的寫(xiě)法是X fft(x, NFFT);x是輸入的時(shí)間序列NFFT是變換點(diǎn)數(shù)。這里有一個(gè)關(guān)鍵點(diǎn)需要理解fft的輸出X是一個(gè)復(fù)數(shù)數(shù)組長(zhǎng)度為NFFT。X(1)對(duì)應(yīng)0Hz直流分量X(2)對(duì)應(yīng)頻率為fs/NFFT的成分X(3)對(duì)應(yīng)頻率為2·fs/NFFT的成分以此類(lèi)推。在X的后半段保存的是負(fù)頻率部分也就是X(NFFT/22)到X(NFFT)對(duì)應(yīng)的是負(fù)頻率到0-的頻率。很多初學(xué)者直接plot(abs(X))最后畫(huà)出來(lái)的頻譜是雙邊譜橫軸范圍從0到fs而且后半段還是鏡像的看起來(lái)非常奇怪。正確的做法是取前半段并把橫軸換算成實(shí)際頻率也就是% 單邊譜處理 NFFT length(x); X fft(x, NFFT); X_single X(1:NFFT/21); X_amp abs(X_single) / NFFT; % 單邊譜的幅值是雙邊譜的兩倍直流分量除外 X_amp(2:end-1) X_amp(2:end-1) * 2; freq (0:NFFT/2) * fs / NFFT;這個(gè)“乘以2”的步驟是另一個(gè)高頻翻車(chē)點(diǎn)。為什么單邊譜要乘以2因?yàn)樨?fù)頻率部分雖然不畫(huà)出來(lái)但它在物理上對(duì)應(yīng)的能量是被解析到正頻率這邊的真實(shí)的正頻率幅值應(yīng)該等于正負(fù)頻率貢獻(xiàn)之和。如果不乘2幅值譜會(huì)恰好偏低一半而很多人做定量分析時(shí)發(fā)現(xiàn)振幅和原始記錄對(duì)不上問(wèn)題很可能就出在這里。3.2 幅值譜、功率譜與相位譜的取舍FFT的結(jié)果是復(fù)數(shù)從中可以提取出三種常用譜幅值譜Amplitude Spectrum就是復(fù)數(shù)模值除以NFFT它給出了信號(hào)在某個(gè)頻率上的“振動(dòng)幅度”有多大單位與原始信號(hào)一致。如果要關(guān)心的是地面運(yùn)動(dòng)峰值加速度或峰值速度就應(yīng)該看幅值譜。功率譜密度Power Spectral Density, PSD則是幅值的平方除以頻率分辨率單位是信號(hào)單位的平方/Hz。它的物理意義是能量的頻率分布密度特別適合對(duì)比不同頻帶內(nèi)的能量大小和信噪比。地震學(xué)里的場(chǎng)地放大效應(yīng)、地脈動(dòng)H/V譜比分析都使用PSD而不是幅值譜。相位譜給出了各頻率成分的相位信息但在絕大多數(shù)地震頻譜分析場(chǎng)景中不是首要關(guān)心對(duì)象因?yàn)榈卣鸩ㄐ问軅鞑ヂ窂接绊懴辔恍畔?fù)雜且不易解釋。只有在做反演或合成波形擬合時(shí)才會(huì)重點(diǎn)用相位。MATLAB里計(jì)算PSD有不止一種方法。最直接的是基于FFT的Welch方法使用pwelch函數(shù)[psd, f] pwelch(x, window, noverlap, nfft, fs);Welch方法的核心思想是把長(zhǎng)信號(hào)切成多段分別做FFT后取平均。這樣做的優(yōu)點(diǎn)是方差小譜線平滑代價(jià)是頻率分辨率變差因?yàn)槊慷巫兌塘?。我?jīng)常在環(huán)境地脈動(dòng)測(cè)量中用它來(lái)判斷微震信號(hào)中的卓越頻率是否有時(shí)間漂移。實(shí)際建議在地震記錄中如果信號(hào)本身比較平穩(wěn)如地脈動(dòng)、環(huán)境振動(dòng)用pwelch效果好如果是一次性瞬態(tài)事件如天然地震或爆破振動(dòng)用整段fft更合適。3.3 零填充、補(bǔ)零與FFT點(diǎn)數(shù)的進(jìn)階用法零填充是另一個(gè)常被誤解的操作。很多人以為把fft點(diǎn)數(shù)設(shè)得很大比如原始數(shù)據(jù)只有2000點(diǎn)卻設(shè)NFFT16384就能“提高分辨率”。嚴(yán)格來(lái)說(shuō)這只能提高頻譜的插值精度讓曲線更平滑并不能把兩個(gè)真實(shí)間隔為0.5Hz的頻率峰區(qū)分開(kāi)。真正做到區(qū)分兩個(gè)頻率峰需要的是更長(zhǎng)的真實(shí)數(shù)據(jù)記錄而不是補(bǔ)零。舉個(gè)例子就明白了假設(shè)你有10秒的記錄采樣率100Hz那么實(shí)際可分辨的頻率間隔是0.1Hz即1/10秒。如果你補(bǔ)零讓FFT點(diǎn)數(shù)變成8192橫軸上的間隔變小了看起來(lái)“分辨率”提高了但物理上兩個(gè)相差0.05Hz的正弦波仍然無(wú)法被區(qū)分它們?cè)谘a(bǔ)零后的頻譜里只會(huì)顯示為一個(gè)寬包絡(luò)。這一點(diǎn)在論文寫(xiě)作中如果處理不當(dāng)很容易被審稿人質(zhì)疑。零填充推薦用法只有兩種一是為了FFT計(jì)算效率把點(diǎn)數(shù)湊成2的冪次二是為了在頻譜圖上找到更精確的峰位置時(shí)做插值顯示。實(shí)際代碼可以這樣做% 湊2的冪次 NFFT 2^nextpow2(length(x)); X fft(x, NFFT);nextpow2會(huì)返回滿足2^n 長(zhǎng)度L的最小n這能讓FFT計(jì)算速度達(dá)到最快但并不是所有的NFFT都必須是2的冪。MATLAB的fft在點(diǎn)數(shù)包含較大質(zhì)數(shù)因子時(shí)速度會(huì)變慢但包含小質(zhì)數(shù)因子2、3、5、7時(shí)速度仍然非常快所以2的冪只是為了省時(shí)間不是硬性要求。3.4 完整的地震數(shù)據(jù)處理流程代碼下面給出一段可以直接復(fù)制運(yùn)行的標(biāo)準(zhǔn)流程。這段代碼我一般在一個(gè)工程地震項(xiàng)目里會(huì)作為模塊反復(fù)調(diào)用輸入是原始地震波形輸出是預(yù)處理后的時(shí)程和單邊幅值譜。function [freq, amp_spectrum, t_clean, x_clean] seismic_fft_analysis(x_raw, fs) % 輸入x_raw為原始地震加速度記錄向量fs為采樣率 % 輸出freq為頻率軸amp_spectrum為單邊幅值譜t_clean為時(shí)間軸x_clean為預(yù)處理后的信號(hào) % 1. 去除趨勢(shì)與均值 x_raw detrend(x_raw(:), constant); % 2. 帶通濾波這里以0.1Hz-40Hz為例按需修改 fl 0.1; fh 40; [b, a] butter(4, [fl/(fs/2), fh/(fs/2)], bandpass); x_filt filtfilt(b, a, x_raw); % 3. 加窗 N length(x_filt); win hanning(N); x_win x_filt .* win; % 4. FFT NFFT 2^nextpow2(N); X fft(x_win, NFFT); X X / mean(win); % 幅值恢復(fù) % 5. 單邊幅值譜 halfN NFFT/2 1; amp abs(X(1:halfN)) / N; amp(2:end-1) amp(2:end-1) * 2; freq (0:halfN-1) * fs / NFFT; % 6. 輸出預(yù)處理后信號(hào) x_clean x_filt; t_clean (0:N-1) / fs; % 7. 繪圖 figure; subplot(2,1,1); plot(t_clean, x_clean); xlabel(時(shí)間 (s)); ylabel(幅值); title(預(yù)處理后的地震記錄); subplot(2,1,2); plot(freq, amp); xlabel(頻率 (Hz)); ylabel(幅值); title(單邊幅值譜); xlim([0, 50]); end這個(gè)函數(shù)充分考慮了前面所有的細(xì)節(jié)去趨勢(shì)、零相移濾波、Hanning窗、幅值恢復(fù)、單邊譜乘2、2的冪點(diǎn)數(shù)優(yōu)化。直接調(diào)用即可基本不會(huì)出錯(cuò)。要注意的是butter濾波器階數(shù)4只是默認(rèn)具體階數(shù)需要根據(jù)頻帶和衰減需求調(diào)整后面避坑部分會(huì)展開(kāi)講。4. 實(shí)操案例用合成地震記錄驗(yàn)證FFT流程4.1 構(gòu)造已知頻譜特征的合成信號(hào)為了檢驗(yàn)代碼的正確性最有說(shuō)服力的辦法是用一個(gè)“已知答案”的信號(hào)來(lái)測(cè)試。假設(shè)我們模擬一段地震記錄其中包含三個(gè)主要頻率成分4Hz、10Hz和25Hz幅度分別為2.0、1.0和0.5采樣率200Hz時(shí)長(zhǎng)30秒。同時(shí)加入白噪聲模擬環(huán)境干擾fs 200; t 0:1/fs:30-1/fs; N length(t); % 合成信號(hào) f1 4; A1 2.0; f2 10; A2 1.0; f3 25; A3 0.5; x A1*sin(2*pi*f1*t) A2*sin(2*pi*f2*t) A3*sin(2*pi*f3*t); x x 0.2*randn(size(t)); % 加噪聲理論上這個(gè)信號(hào)的頻譜在4Hz、10Hz、25Hz處應(yīng)該有明顯的峰峰值約為2.0、1.0、0.5均方根振幅會(huì)略低因?yàn)樵肼暞B加后能量重新分配。如果我們的FFT流程處理正確這三個(gè)峰的幅值應(yīng)當(dāng)非常接近理論值。4.2 運(yùn)行流程代碼并解讀結(jié)果把上面的x和fs代入seismic_fft_analysis函數(shù)觀察輸出的頻譜圖能得到三個(gè)清晰的峰。4Hz處幅值接近2.0510Hz處接近1.0325Hz處接近0.52與理論值之間的誤差主要來(lái)自隨機(jī)噪聲的疊加。這說(shuō)明整條處理鏈路的幅值標(biāo)定是準(zhǔn)確的。如果你不乘2三個(gè)峰的幅值會(huì)變成大約1.0、0.5、0.26一下子少了一半這就驗(yàn)證了前面說(shuō)的單邊譜乘2的步驟確實(shí)不能省。如果不做幅值恢復(fù)峰幅值也會(huì)系統(tǒng)性偏低Hanning窗的均值是0.5那么所有峰幅值都會(huì)打?qū)φ垡彩敲黠@錯(cuò)誤。4.3 用pwelch做功率譜密度估算對(duì)比如果改用pwelch驗(yàn)證[psd, f_psd] pwelch(x, hanning(512), 256, 1024, fs); plot(f_psd, psd);頻率分辨率大約為fs/5120.39Hz三個(gè)頻率峰照樣能被看到但峰的寬度比直接用整段FFT更寬一些。這是welch分段平均導(dǎo)致的它的好處是譜線平滑適合觀察寬頻背景噪聲但壞處是頻率上的精細(xì)結(jié)構(gòu)被抹平。所以對(duì)于研究尖峰明顯的線譜整段FFT更合適對(duì)于連續(xù)譜、隨機(jī)振動(dòng)pwelch更穩(wěn)。兩者配合使用能互相驗(yàn)證結(jié)論的可靠性。5. 地震記錄頻譜分析中的常見(jiàn)問(wèn)題與避坑指南5.1 頻譜泄漏與窗函數(shù)的“治標(biāo)不治本”頻譜泄漏是FFT處理中最常見(jiàn)的問(wèn)題。典型的癥狀是本來(lái)應(yīng)該在某個(gè)頻率上的一個(gè)尖峰變成了在它附近一坨小突起主峰兩側(cè)還附帶振蕩的旁瓣。泄漏的根源是截?cái)?。任何有限長(zhǎng)信號(hào)在邊界處都是突變的FFT把這個(gè)突變強(qiáng)行當(dāng)成周期信號(hào)的一部分于是原本只有單一頻率的正弦波突然多了許多高頻成分來(lái)“擬合”這個(gè)突變。加窗能緩解邊界突變但不同窗函數(shù)的抑制能力差異很大矩形窗泄漏最嚴(yán)重Hanning次之Blackman-Harris窗旁瓣衰減最干凈但主瓣最寬。我一般遇到能量相差很大的兩個(gè)信號(hào)源同時(shí)出現(xiàn)時(shí)會(huì)用Kaiser窗并把β值調(diào)大效果比固定窗好很多。但要說(shuō)清楚窗是“治標(biāo)”真正的“治本”是讓截取窗口內(nèi)的信號(hào)本身盡可能平穩(wěn)。如果地震記錄里含有明顯的震相突變比如初至P波到達(dá)時(shí)振幅突然跳變那么在這個(gè)跳變點(diǎn)上必然會(huì)產(chǎn)生大量高頻泄漏。正確做法是只選P波到達(dá)前的噪聲段分析背景噪聲或者只選S波之后的尾波段分析地脈動(dòng)而不是把整段波形不分青紅皂白直接做FFT。5.2 濾波階數(shù)與filtfilt邊界效應(yīng)很多人看到butter函數(shù)隨手填個(gè)階數(shù)8或10覺(jué)得階數(shù)越高濾波越“干凈”。但實(shí)際上高階Butterworth濾波器會(huì)帶來(lái)嚴(yán)重的相位延遲和數(shù)值穩(wěn)定性問(wèn)題而且filtfilt一次處理下來(lái)邊界效應(yīng)會(huì)加倍。我曾經(jīng)在處理一批強(qiáng)震記錄時(shí)用了10階帶通結(jié)果信號(hào)前50個(gè)點(diǎn)和后50個(gè)點(diǎn)出現(xiàn)了明顯的“飛邊”頻譜也出現(xiàn)高頻震蕩的假象排查半天才發(fā)現(xiàn)是濾波器階數(shù)過(guò)高。根據(jù)我的經(jīng)驗(yàn)帶通濾波器階數(shù)4~6足夠應(yīng)付絕大多數(shù)地震數(shù)據(jù)場(chǎng)景。如果濾波需求非常窄帶比如提取0.2Hz~0.3Hz的窄帶信號(hào)可以改用Chebyshev II型或Elliptic濾波器它們的通帶波紋和阻帶衰減特性更適合窄帶提取但要注意群延遲會(huì)變得不均勻。實(shí)在沒(méi)辦法的時(shí)候也可以考慮用最小二乘擬合的時(shí)域?yàn)V波器計(jì)算速度慢但控制精度極高。另外filtfilt邊界效應(yīng)有一個(gè)實(shí)用對(duì)策在濾波前把信號(hào)兩端各延拓一段例如每端加200個(gè)點(diǎn)延拓值取信號(hào)首尾的均值并用窗函數(shù)平滑過(guò)渡。濾波完成后裁剪掉延拓部分。這個(gè)做法能顯著減少邊界的瞬時(shí)振蕩。5.3 采樣率不一致導(dǎo)致諧波錯(cuò)位有時(shí)候你的地震記錄不是自己采的而是從不同儀器上導(dǎo)出的。有的儀器采樣率是100Hz有的可能是120Hz有的記錄由于時(shí)鐘漂移導(dǎo)致實(shí)際采樣率偏離標(biāo)稱(chēng)值。如果你把所有記錄用同一個(gè)標(biāo)稱(chēng)采樣率代入FFT頻譜的橫軸就會(huì)整體偏移表現(xiàn)為同一個(gè)已知頻率峰的“漂移”。排查方法很簡(jiǎn)單找一個(gè)記錄中已知的穩(wěn)定頻率源比如50Hz交流電干擾或某個(gè)已知諧波信號(hào)做標(biāo)定。如果你的頻譜中50Hz峰顯示成52Hz那就說(shuō)明采樣率實(shí)際偏高了4%反過(guò)來(lái)就要校正時(shí)間軸。多數(shù)現(xiàn)代的SAC或miniSEED格式文件頭里都記錄了采樣率但轉(zhuǎn)換過(guò)程中容易丟失或誤寫(xiě)處理前養(yǎng)成檢查head的快照習(xí)慣非常有用。MATLAB里可以用auftach或SAC相關(guān)工具讀取頭段確認(rèn)采樣率沒(méi)有歧義。5.4 長(zhǎng)記錄分段處理與內(nèi)存優(yōu)化一臺(tái)高采樣率連續(xù)記錄儀一天就會(huì)產(chǎn)生約1728萬(wàn)點(diǎn)數(shù)據(jù)假設(shè)200Hz24h。這么長(zhǎng)的信號(hào)如果一次性做FFT不僅計(jì)算慢而且頻率分辨率極高卻毫無(wú)意義因?yàn)榈皖l段的細(xì)微變化不需要全局分辨率倒是高頻段的非平穩(wěn)細(xì)節(jié)需要局部化處理。處理長(zhǎng)記錄的正確思路是分段。分段長(zhǎng)度按照目標(biāo)頻段來(lái)決定如果只是分析0.5Hz以上的短周期振動(dòng)用5~10秒一段做平均如果要分析0.01Hz量級(jí)的固體潮或長(zhǎng)周期面波可能需要幾十分鐘甚至更長(zhǎng)的一段數(shù)據(jù)才能獲得足夠分辨率。另一方面分段之間可以設(shè)置50%的重疊來(lái)減少段首段尾的影響這是Welch方法的標(biāo)準(zhǔn)配置。在MATLAB中處理大矩陣FFT時(shí)還有個(gè)容易忽略的性能殺手fft對(duì)列向量和矩陣的處理方式不同。如果X是一個(gè)N行多列的矩陣fft(X)會(huì)對(duì)每一列分別做FFT因此三分量數(shù)據(jù)可以直接拼成N×3矩陣一次性變換比循環(huán)三次快很多。內(nèi)存占用方面N點(diǎn)FFT的中間復(fù)數(shù)數(shù)組約需要16×N字節(jié)一般幾百兆以?xún)?nèi)的數(shù)據(jù)都不會(huì)有壓力但如果是長(zhǎng)記錄多通道分析建議用single類(lèi)型來(lái)減半內(nèi)存精度損失對(duì)頻譜分析來(lái)說(shuō)完全可以接受。5.5 頻譜圖可視化中的比例尺與縱軸選擇最后一個(gè)常見(jiàn)“坑”是畫(huà)圖方式誤導(dǎo)解讀。不少人在畫(huà)地震頻譜時(shí)直接用線性縱軸結(jié)果主頻太高把低幅值的背景信息壓成了一團(tuán)“零線”有人用對(duì)數(shù)縱軸又過(guò)分放大噪聲。正確做法是根據(jù)分析目的選擇縱軸如果要突出能量集中的主頻用線性縱軸合適如果要看全頻帶的衰減趨勢(shì)最好用對(duì)數(shù)dB縱軸。另外如果不特別說(shuō)明很多人畫(huà)頻譜圖時(shí)縱軸是普通的1/Hz密度或原始幅值但科學(xué)論文里通常要求標(biāo)注單位。比如加速度記錄的PSD單位是(m/s2)2/Hz幅值譜單位是m/s2。我在自己的腳本中會(huì)把縱軸標(biāo)簽和單位直接內(nèi)置避免后期返工。橫軸也建議默認(rèn)畫(huà)到奈奎斯特頻率但是要按需限制顯示范圍比如目標(biāo)是看1~20Hz的工程頻段就不要把0~100Hz整段畫(huà)出來(lái)那樣會(huì)浪費(fèi)幅面而且看不清細(xì)節(jié)。6. 地震FFT分析的延伸應(yīng)用與工具箱搭配6.1 從加速度記錄計(jì)算反應(yīng)譜時(shí)的FFT思路工程地震里經(jīng)常需要從一條加速度時(shí)程計(jì)算阻尼反應(yīng)譜。雖然反應(yīng)譜的計(jì)算通常用Newmark-β法等時(shí)域方法或杜哈梅積分但FFT可以大幅加速?gòu)椥苑磻?yīng)譜的計(jì)算尤其當(dāng)結(jié)構(gòu)自振周期非常多、數(shù)量達(dá)到幾百個(gè)時(shí)時(shí)域循環(huán)會(huì)非常慢。快速解法是把加速度記錄一次性變換到頻域再用結(jié)構(gòu)頻響函數(shù)乘以地震波頻譜最后做一次逆FFT得到結(jié)構(gòu)位移、速度和加速度時(shí)程。這個(gè)過(guò)程本質(zhì)上是頻域求解線性振動(dòng)方程比逐周期計(jì)算快了不止一個(gè)量級(jí)。如果對(duì)計(jì)算精度要求高需要注意微分算子在頻域中表示為乘以jω而加速度到速度是除以jω零頻處會(huì)出現(xiàn)奇異點(diǎn)必須先對(duì)頻譜做低截處理去除長(zhǎng)周期漂移。6.2 結(jié)合H/V譜比法評(píng)估場(chǎng)地卓越頻率H/V譜比法是當(dāng)前場(chǎng)地效應(yīng)評(píng)估里很簡(jiǎn)單有效的工具核心思想是對(duì)同一時(shí)間段的地表三分量記錄分別做FFT得到水平向和垂直向的傅里葉幅值譜然后計(jì)算水平向平均譜除以垂直向譜的比值。H/V譜中的峰值對(duì)應(yīng)的頻率通常就是場(chǎng)地的卓越頻率。實(shí)現(xiàn)H/V譜比時(shí)FFT參數(shù)的選擇非常重要。經(jīng)驗(yàn)表明分析窗口長(zhǎng)度至少應(yīng)包含100個(gè)目標(biāo)頻率的周期否則分辨率不足。比如場(chǎng)地卓越頻率如果是1Hz那么窗口至少40~100秒才合適。此外各段取的窗口長(zhǎng)度要一致否則譜比會(huì)出現(xiàn)人為的“毛邊”。可以用前面介紹的分段pwelch方法分別計(jì)算三個(gè)分量的PSD再開(kāi)方轉(zhuǎn)成幅值譜最后相除這樣平滑效應(yīng)比較好曲線也穩(wěn)定。6.3 MATLAB工具箱的替代方案與效率對(duì)比MATLAB原生的Signal Processing Toolbox已經(jīng)覆蓋了絕大多數(shù)FFT相關(guān)需求不需要為了頻譜分析特地去安裝額外工具箱。如果確實(shí)需要更高級(jí)的分析比如短時(shí)傅里葉變換(STFT)、小波變換、希爾伯特黃變換(HHT)需要額外的Wavelet Toolbox或自己寫(xiě)代碼。STFT是FFT的滑動(dòng)窗口變體在時(shí)頻圖上可以看到不同時(shí)刻的頻率變化對(duì)震相識(shí)別非常有幫助。MATLAB的spectrogram函數(shù)直接可用不用額外工具箱。如果項(xiàng)目數(shù)據(jù)規(guī)模特別大或者需要和地震學(xué)專(zhuān)業(yè)軟件打通可以考慮用SACSeismic Analysis Code做前期預(yù)處理將預(yù)處理后的波形通過(guò)格式轉(zhuǎn)換導(dǎo)出為MATLAB格式再做FFT分析。SAC在時(shí)間域文件頭處理和濾波上有更高的自由度而MATLAB強(qiáng)在可視化和自定義迭代計(jì)算。兩者結(jié)合是一種很順手的組合拳我在處理一批連續(xù)波形微震數(shù)據(jù)時(shí)經(jīng)常這么配合。6.4 逆FFT恢復(fù)信號(hào)時(shí)的注意事項(xiàng)FFT不只是從時(shí)間域到頻域有時(shí)也要從頻域回到時(shí)間域比如濾波操作本質(zhì)上是頻域乘以一個(gè)譜窗再逆變換回時(shí)域。MATLAB的ifft函數(shù)會(huì)把復(fù)數(shù)頻譜恢復(fù)成時(shí)間序列。逆FFT的坑和正變換對(duì)應(yīng)如果你修改了頻譜比如把某個(gè)頻段歸零那重建的信號(hào)可能不再是實(shí)信號(hào)而是帶有虛部的小量。這時(shí)應(yīng)該用real(x_ifft)提取實(shí)部同時(shí)應(yīng)該意識(shí)到對(duì)頻譜做過(guò)零點(diǎn)切除之后時(shí)域信號(hào)兩端會(huì)自動(dòng)出現(xiàn)振鈴這是因?yàn)闉V波器在頻率域的突變對(duì)應(yīng)時(shí)域的sinc函數(shù)卷積。所以頻域?yàn)V波的截止頻率兩端要盡量平滑過(guò)渡給一個(gè)過(guò)渡帶振鈴會(huì)小很多。我屢次在用頻域方法去除地脈動(dòng)記錄中的機(jī)械噪聲時(shí)發(fā)現(xiàn)平滑過(guò)渡帶比生硬切除重要得多直接截?cái)鄤t會(huì)在波形上留下人眼可見(jiàn)的一系列共振式波紋。7. 后續(xù)還能往哪個(gè)方向擴(kuò)展如果這段FFT地震頻譜分析的流程你已經(jīng)跑通了下一步可以考慮的方向很多。一是把批處理能力做起來(lái)比如面對(duì)上百條波形記錄時(shí)用一個(gè)循環(huán)統(tǒng)一完成預(yù)處理和頻譜提取并把結(jié)果輸出成結(jié)構(gòu)數(shù)組或表格。二是在頻域里加入多通道交叉分析比如計(jì)算兩個(gè)臺(tái)站同一地震記錄在頻域內(nèi)的相干性就能估計(jì)波速和衰減參數(shù)這是地震層析成像的前置步驟之一。三是從頻域反演混合信號(hào)中的震源譜項(xiàng)和路徑效應(yīng)項(xiàng)這是開(kāi)展震源物理研究的地基。我個(gè)人在實(shí)際操作中最想提醒大家的一句經(jīng)驗(yàn)是FFT本身是一個(gè)數(shù)學(xué)工具算法層面幾乎沒(méi)有門(mén)檻真正的門(mén)檻全在預(yù)處理和參數(shù)選擇上。同一個(gè)地震記錄濾波參數(shù)不同、窗函數(shù)不同、FFT點(diǎn)數(shù)不同畫(huà)出來(lái)的頻譜差別會(huì)非常大甚至可能得出完全相反的結(jié)論。所以在整個(gè)頻譜分析流程中最值得花時(shí)間的不是把fft代碼跑通而是把你手里的信號(hào)“伺候”舒服讓它能干凈地進(jìn)入FFT。當(dāng)你發(fā)現(xiàn)自己的頻譜圖主頻變得清晰、旁瓣消失、幅值符合物理直覺(jué)時(shí)這套流程才算真正過(guò)了關(guān)。如果哪天你遇到頻譜形態(tài)怎么都解釋不通的案例不妨回頭看一眼我們上面聊過(guò)的每一個(gè)細(xì)節(jié)大概率問(wèn)題就藏在你忽略的那一步里。希望這篇文章能幫你少走一些彎路早點(diǎn)把心念已久的地震頻譜圖做出來(lái)。本文還有配套的精品資源點(diǎn)擊獲取