加速度積分得到速度和位移的問題小結

行路難

2018年8月17日 13:26

瀏覽：3992 收藏：4

1 前言

在振動信號采集分析中，我們常常使用加速度傳感器，是因為加速度傳感器具有可靠性高，可測頻帶寬，結構小巧，抗干擾能力強等優點。但是有些時候，必需要得到速度和位移信號，這個時候就出現了如何通過對加速度積分得到速度和位移的基礎問題，這個問題難倒了不少相關人員。在國內的各種教科書上也很少提及，也不知道是為什么。

不管是線性系統還是非線性系統，所產生的信號，在時域滿足如下微分方程

其解可以寫成

其中v₀是積分常數，積分下限t=0表示采樣起始時刻，是人為指定的時間軸原點。如果振動沒有交流分量，則加速度是常數，可以得到v(t)=v₀+at 。對速度積分得到位移的處理類似，所以有

其中d₀是積分常數，由于v₀和d₀通常是未知的，加速度積分得到速度和位移，就存在積分常數問題。

積分常數問題是非常麻煩的問題。首先可以看到，對于任意δv，得到

都對應同樣的加速度a(t)，也就是說v₀和d₀取任意值都對應于同一個測量的加速度和測量的速度，所以要想僅僅憑目前的信息通過積分獲取速度和位移，是不可能的，不管是數值積分還是模擬電路積分都不行。本文主要討論數值積分的問題，關于模擬電路積分存在的問題，推薦大家閱讀《壓電加速度計和振動前置放大電器》一書，結論就是在高頻，放大器增益會下降，不是理想放大器；在甚低頻，會受到時間常數或者回路電容限制，也不夠精確。這個回復大家可能不滿意，坊間不是有各種FFT，各種去趨勢項，各種濾波的操作嘛？有人發表在《力學學報》上的論文不也是用模擬積分電路積分了么？我的建議是，謹慎看待那些論文中的說法，你應該有自己的思考和判斷。

實際工作中，在無窮多個v₀的選擇中，可以選擇一個最靠譜的數值湊合用一下，因為在這些選項中，可確定最靠譜的那個值一定在某個有限幅值的區間[v_a，v_b]以內。所以，需要增加對被測試物體的狀態進行判斷，基于這種判斷，做出假設，基于這種假設列出約束條件，根據約束條件找到一個方程，然后求出那個看起來好像最靠譜的積分常數。如果假設成立，就得到了一個很Nice的估計；如果假設不成立，誤差很大，那就是自娛自樂，你高興就好。

2 從可觀測性角度看這個問題

下面的討論限于線性系統，是最近關注某個問題的時候順便做的筆記。

設振動方程為

令

則，對應的狀態空間模型為

其中

可以得到離散以后的狀態方程

其中

可觀性的定義：假定A，B，C和D已知，從已經觀測的加速度和載荷向量，能否計算得到位移和速度？從式的迭代解可以知道，只需要設法測試得到x(1)即可知道所有后續的狀態x(k),k=2,3,4,…。

現在，考慮如下序列

寫為矩陣形式

注意到對SDOF模型，上面的矩陣方程始終是一個p×(p+2)矩陣，屬于欠估計問題，不可能通過測試加速度同時將載荷和初始狀態x(1)求出來。現將該矩陣分割為

其中O_p就是可觀測矩陣，可以求得

也就是說要聯合測試加速度響應，知道可觀測矩陣O_p、系統矩陣E和載荷才能識別得到所有狀態變量，而且可觀測矩陣不能奇異，必需是滿秩的。2015年Computers &Structures上的一篇關于載荷識別的論文，就是依據對式(8)的簡化，假定初始狀態x(1)=0，得到

這種零初始條件的假定，在實際問題的應用中，是一廂情愿的。在數據采集器采樣了一個數據塊處理以后，不能再假定這個數據塊采集終了時候的狀態變量還是零，所以，至少在后續的計算中，需要考慮 “過去的激勵引起的未來的響應的問題”或者x(1)不為0的問題。

同樣的道理，實際情況下，我們并沒有測試激勵力?，還沒有做實驗獲取系統參數，更不能保證傳感器布置點所對應的觀測矩陣是Regular的，所以單獨想要用加速度響應，就完全構造出速度和位移是不可能實現的。就算你得到了，誤差可能在80%左右，正確的是量綱。

現在的問題是，我們得到了加速度數據，到底應該考慮哪些實際問題，有沒有辦法得到一個不是很壞的結果？答案是肯定的，只是這些要求太高了或者假設成立的可能性微乎其乎，所以得到的數據的誤差通常較大。

3 可能的假設

J.S Bendat和A.G. Piersol對信號劃分為確定性的和隨機的，一個信號是否是確定性的還是隨機的，可能存在爭議，短時間可能是確定性的，而長時間就變得隨機了，或者長時間以后就出現其他“開關條件”觸發了不可預測的因素。對確定的信號和隨機信號，又分為了如下幾種：

在確定性信號中，周期信號是稀有的；在隨機信號分類中，各態歷經的平穩信號是稀有的。最常見的信號是非周期信號和非平穩信號，是實際中遇到的高概率事件。

所以，下面討論一些“低概率”出現的假設和“高概率誤操作”。

3.1 假設1：位移是零均值的信號

被測試結構沒有發生明顯的位移現象，所以假定是零均值的。然而，我們知道，有限時間采樣以后，得到的是離散數據，這些離散數據的均值并不能嚴格等于采樣區間連續信號的真實均值。而且采樣區間真實的均值，極有可能是非零均值的，只是均值很小而已。對速度信號，進行數值積分獲取零均值位移信號的方法，同下面的對加速度數值積分獲取速度的方法。

3.2 假設2：速度是零均值的信號

設當加速度不為常數時，采樣間隔時間dt足夠小，一種離散的數值積分為梯形公式

v₀任意假定。零均值的獲取就是去掉趨勢項，得到

上面的式(10)是一種程序語言，代表將有偏估計的速度信號(組成的矩陣)，每個元素去除整體上的已知均值，這樣修正得到的v就是零均值的。去掉趨勢項的操作，可能存在各種優雅的處理，比如Matlab里面的detrend命令等等，這里不再涉及。

值得一提的是，如果去趨勢項的運算是含有非線性因素，那么可能就存在頻率混疊的困難，因為采樣頻率是一定的，這些非線性運算，比如帶入平方，開平方根，矩陣求逆，或者取符號運算等等非線性的操作，會引入帶寬以上的頻率分量，而這可能不是希望的結果，需要后續修正，最好的處理是先保證足夠的帶寬，最后做無相移的低通濾波。K. Worden教授對去趨勢項的研究結論認為，那些“優雅的”去趨勢項操作引入的誤差很難處理。有趣的是，他早年在MSSP上發表的論文認為，數值積分方案是很Nice的選擇；幾年以后在他自己的書里對這個問題又做了更細致的論證和分析，他認為數值微分方案才是更好的。

3.3 假設3：速度是周期信號

這種情況，對加速度FFT，在頻域積分，再IFFT回來就可以了。

但是我們要知道兩點潛在的問題：加速度信號，從時域圖上看，一定要是周期的或者Almost周期的，所以分辨率一定要選擇得合理才能做到整周期采樣；如果頻率分辨率選擇不合理，可能存在很嚴重的泄露，加窗號稱可以減小泄露，但是，加窗通常對積分以后的速度信號是毀滅性的，加窗的效應保留在時域速度信號里面，無法修正回來；如果有本事做到不加窗，對數據的分辨率選擇就提出了很高的要求。如果不是周期信號，怎么選擇分辨率都沒用啊。

圖1 非周期信號+加窗的FFT積分再IFFT結果，以及和精確解(藍色)對比

3.4 假設4：速度為零時，加速度最大

這個假設，是本人偶然想到的，是在研究某算法時，看到某非線性微分方程強迫響應，然后突然想到，假定速度為零時，對應于動能最小，然后慣性力最大，產生的加速度也是最大的，就這樣得到一個結果，在那個例子中，需要在眾多的加速度峰值中選擇一個，然后確定積分常數，確實產生了較好的效果，但是這個實在太巧合了，再次按照這個方法選擇其他峰值做一遍，失效了，現實中是沒有這么理想的假設的。

圖2 速度(紅色)和加速度（藍色）時域曲線對比

當時的做法是，在某一個明顯的加速度峰值附近，選取幾個離散加速度點；用多項式進行擬合，得到加速度關于時間的表達式，這一步可以表達為

其中下標k=1,2,3...s只是代表選擇的鄰近的有限個數據點。

假定系統的阻尼很小，在最大慣性力和最大彈性恢復力平衡時粘性阻尼力為零，這個假設在某些時候是成立的，那么加速度最大的時候，速度可以為零。令上面的表達式為零，選擇離散數據中加速度最大的時刻t_max，解方程找到加速度真正的峰值時刻t_p，t_p非常靠近t_max，即

對加速度的表達式積分得到速度的表達式，積分常數v₀待求

t_p時刻的速度為零，從而解出積分常數v₀

在這個區間的速度表達式就已知了，回到t_max時刻，可以直接得到這個時刻下的速度，應該是一個非常小的量

這樣，我們終于找到了“估計的速度小量”，可以根據上面的離散公式，繼續完成剩下的時間積分。

因為這個方法并不總是成立，失效以后的某個結果如下所示。

圖3 梯形數值積分結果和參考結果(紅色)對比，存在一個整體偏移量δv

4 一個數值積分案例

在網上蕩到了一個96年Steven D. Glaser等人在(Preliminary Processing of the Lotung LSST Data報告中，對一個測試到的加速度數據進行數值積分的案例，附錄中有源代碼，這個數據是地震加速度信號，頻率特別低，把源代碼貼在這里，給感興趣的人：

function[a,v,d]=vd(f,Lcut,dT,pre,name)

% function[a,v,d]=vd(f,Lcut,dT,pre,name)

% detrends andfilters the input acc and integrates twice to give velo and disp

% f is the inputacc vector

% Lcut is thelow-end cutoff frequence in Hz

% dT is the timestep

% n is the orderof the Butterworth filter

% pre is thepre-event segment length to be zeroed

zip=1:pre;

zip=zeros(size(zip));

[b,c]=butter(n,Lcut*dT*2,'high');

a=filtfilt(b,c,f);

a=dtrend(a(:,1),1,pre);%%mxl:這個函數已經過時了.

a=detrend(a);

a=a-a(pre+1);

a(1:pre)=zip;

v=inttrap(a,dT);%%mxl:這個是數值積分，這個函數過時了

v=filtfilt(b,c,v);

v=dtrend(v(1,:)',1,pre);%%mxl:這個函數已經過時了

v=detrend(v);%去趨勢項

v=v-v(pre+1);%整體平移量

v(1:pre)=zip;%前pre個數值設置為0

d=inttrap(v,dT);%%mxl:這個是數值積分，這個函數過時了

[b,c]=butter(n,Lcut*dT*2,'high');

d=filtfilt(b,c,d);

d=dtrend(d(:,1),1,pre);

d=detrend(d); %去趨勢項

d=d-d(pre+1); %整體平移量

d(1:pre)=zip; %前pre個數值設置為0

可以看到，他們首先使用了低通濾波，做了去趨勢項，然后整體平移了一個量，再將前pre個數值人為置為零。

5 終極建議

如果同時測試了加速度和位移，那么就特別好辦了。所以，如果有條件，盡可能測試加速度和位移。

設在某個區間 [t₁,t₂] 里面，加速度和位移都可以擬合一個表達式(expression)

在兩次積分的情況下可以得到時間常數的表達式，從而解出未知系數

可見速度也同時得到了。

測試位移+測試加速度，可以通過數值微分+低通濾波得到速度；當然也可以直接通過Kalman濾波得到速度，方便、快捷、高效又穩健；或者帶有微分算法的低通濾波得到速度，等等。方法總是有很多種。但并不建議測試位移，兩次數值微分得到加速度，因為數值微分算法是高通濾波器，會引入高頻噪聲，采樣頻率越高，高頻噪聲問題越嚴重，兩次數值微分得到的加速度誤差較大，但也并不是完全不能使用，加上低通濾波可以降低這個誤差，總比從加速度積分得到位移要靠譜得多。使用低通濾波算法要注意不要引入相移，平滑處理（Smooth）也是一種低通濾波，效果還行。下面是把不同點數差分算法看成是高通濾波器，然后計算其頻率響應得到的曲線，注意是對采樣頻率f_s做了歸一化以后的結果。結論是：用5點差分算法是可行的，向前差分或者Matlab的diff函數，最糟糕。