加速度如何在頻域上積分成位移?

01-08

加速度如何在頻域上積分成位移，自己本身學的是結構，老闆搞地震的，所以跟著老闆搞地震。。現在我想通過一條實測的地震波加速度樣本去積分出位移。。之前想過用時域的辦法發現相位差很大而且處理很不方便。我想通過頻域的辦法。。之前看過文獻發現引入傅立葉變換後求導積分真是方便。因為之前沒有具體接觸過matlab 有沒有人能夠指點下在matlab中實現我這個想法。

原理：

時域中，加速度時程積分兩次可以得到位移時程。在頻域中，將加速度的傅里葉譜稍加調整後可以得到位移的傅里葉譜。

上面兩式均省略了常數項，常數項的處理此處略去不表。

方法：
傅立葉分析是將一個函數用餘弦函數的平移和伸縮來表示。將時域中的加速度時程經過傅里葉變換（fft）後得到加速度的傅里葉譜。將加速度傅里葉譜除以譜值對應頻率的平方，再將譜值取負可以得到位移的傅里葉譜。最後利用傅里葉逆變換（ifft）得到時域中的位移時程。

附錄：
附上《MATLAB在振動信號處理中的應用》中基於頻域積分利用加速度時程求速度/位移時程的Matlab程序[1]：
clear;clc;close all;
%%%%%%%%%%%%%%%%%%%%%%%
fni=input("頻域積分-輸入數據文件名:","s");
fid=fopen(fni,"r");
sf=fscanf(fid,"%f",1); % 採樣頻率
fmin=fscanf(fid,"%f",1); % 最小截止頻率
fmax=fscanf(fid,"%d",1); % 最大截止頻率
c=fscanf(fid,"%d",1); % 單位變換係數
it=fscanf(fid,"%s",1); % 積分次數
sx=fscanf(fid,"%s",1); %橫坐標的標註
sy1=fscanf(fid,"%s",1); %縱坐標輸入單位的標註
sy2=fscanf(fid,"%s",1); %縱坐標輸出單位的標註
fno=fscanf(fid,"%s",1); %輸出數據的文件名
x=fscanf(fid,"%f",[1,inf]); % 按行讀入原始信號數據
status=fclose(fid);
% 計算輸入數據的長度
nt=length(x);
% 建立時間向量
t=0:1/sf:(n-1)/sf;
% 取大於n且與其最接近的2的整數次方為FFT長度
nf=2^nextpow2(nt);
% FFT 變換
y=fft(x,nfft);
% 計算頻率間隔
df=sf/nfft;
% 計算指定頻帶對應頻率數組的下標
ni=round(fmin/df+1);
na=round(fmax/df+1);
% 計算原頻率間隔
dw=2*pi*df;
% 建立正的離散原頻率向量
w1=0:dw:2*pi*(0.5*sf-df);
% 建立負的離散原頻率向量
w2=2*pi*(0.5*sf-df):-dw:0;
% 將正負付給一個數組
w=[w1,w2];
%以積分次數為指數，建立元頻率向量
w=w.^it;
% 進行積分的頻域變換
a=zeros(1,nfft);
a(2:nfft-1)=y(2:nfft-1)./w(2:nfft-1);
if it==2
% 進行二次積分的相位變換
y=-a;
else
% 進行一次積分的相位變換
real(y)=imag(a);
imag(y)=-real(a);
end
a=zeros(1,nfft);
% 消除指定正頻帶外的頻率成分
a(ni:na)=y(ni:na);
% 消除指定正頻帶外的頻率成分
a(nfft-na+1:nfft-ni+1)=y(nfft-na+1:nfft-ni+1);
% ifft變換
y=ifft(a,nfft);
% 積分結果
y=real(y(1:n))*c;
subplot(211)
plot(t,x);
xlabel(sx);ylabel(sy1);
grid on;
% 積分後曲線
subplot(212)
plot(t,y);
xlabel(sx);ylabel(sy2);
grid on;
% 打開文件輸出積分後的數據
fid=fopen(fno,"w");
for k=1:n
fprint(fid,"%f%f
",t(k),y(k));
end
status=fclose(fid)

參考文獻：

[1]《MATLAB在振動信號處理中的應用》，王濟，中國水利水電出版社，知道產權出版社。

（小木蟲的下載地址：《matlab在振動信號處理中的應用》（教程+程序））

Chopra 的動力學課本 Section 6.1 Earthquake excitation 結尾部分有專門的介紹，加速度積分得到速度，速度再積分一次得到位移。

用 MATLAB 的話，可以先在時域里把加速度數值積分，得到速度，比如可以用 cumtrapz；然後把速度 fft 轉換到頻域，觀察速度的頻率分布，理想狀態下的頻率分布會集中在自振周期附近，但很多時候有些比較低或者比較高的頻率會遠離自振周期，比如可能來自於實測數據里的背景噪音或者極高頻的噪音，這時候就要用 filter 把這些頻率過濾掉，比如可以用 Butterworth filter。過濾之後的速度信號，頻率應該就會集中在自振周期附近。

同樣，把速度再數值積分得到位移，然後 fft，然後 filter，最終得到過濾之後的位移。通常來說，地震工程關心的是有物理意義的最大數值，所以最終結果在時域里是更有意義的，整個時域里加速度、速度、位移的最大值都是我們關心的關鍵指標。

剛編了個程序，就是按照Chopra那本書中的辦法。

Newmark-beta法，Wilson-sita法。MATLAB程序進行迭代時，絕大多數都是二者的引伸。