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