基于功率譜估計的伺服系統(tǒng)機械諧振分析

朱琰虹，文光武

（1廣州數(shù)控設備有限公司，廣州 510530；2廣州數(shù)控信息科技有限公司，廣州 510603）

0 引言

機械諧振是伺服系統(tǒng)中十分常見的一種現(xiàn)象。伺服系統(tǒng)的機械傳動部分經常使用傳動軸、變速器、聯(lián)軸器等傳動裝置連接電機和負載，而實際傳動裝置并不是理想剛體，存在一定的彈性，通常會在系統(tǒng)中引發(fā)機械諧振[1]。當控制系統(tǒng)輸入的激勵力的強迫振蕩頻率接近機械系統(tǒng)的自然諧振頻率時，會導致機械系統(tǒng)在共振頻率附近產生劇烈震蕩。實際機械系統(tǒng)中柔性耦合因素的改變會改變系統(tǒng)的穩(wěn)定裕度、增加系統(tǒng)的控制難度，甚至會使系統(tǒng)不再穩(wěn)定，造成裝置劇烈震動，導致機械結構受到沖擊和變形，極易發(fā)生破壞性事故[2]。現(xiàn)代伺服技術中通常對電機轉速或轉子電流進行高速采樣，再利用傅里葉變換求得信號的幅度譜和相位譜來分析機械諧振[3]。本文考慮到機械諧振產生的隨機性，所采樣信號是確定性信號和隨機信號的疊加，而隨機信號是無始無終的并具有無限能量，不滿足絕對可積的條件[4]，因此可以通過研究其功率在頻域上的分布，即功率譜密度或功率譜，從而實現(xiàn)對伺服系統(tǒng)機械諧振的分析。

1 隨機信號

隨機信號可分為平穩(wěn)隨機信號和非平穩(wěn)隨機信號。各態(tài)歷經信號是指無限個樣本在某時刻所歷經的狀態(tài)，等同于某個樣本在無限時間里所經歷的狀態(tài)的廣義平穩(wěn)隨機信號[5]?，F(xiàn)實的隨機信號為大部分可逼近各態(tài)歷經的平穩(wěn)隨機信號，在各態(tài)歷經情況下，離散隨機信號下的集合平均等于時間平均。隨機信號不能用確定性的時間函數(shù)來描述，只能用統(tǒng)計的方法研究，其統(tǒng)計特性通常用均值、方差、相關函數(shù)與協(xié)相關函數(shù)來表征。

離散隨機信號x(n)的自相關函數(shù)如式（1）。

兩離散隨機信號 x(n)和y(n)之間的互相關函數(shù)如式（2）。

除去均值后隨機信號x(n)的相關性用自協(xié)方差函數(shù)如式（3）。

兩離散隨機信號x(n)和y(n)之間的協(xié)方差用互協(xié)方差函數(shù)來描述，如式（4）。

式（1）、式（2）、式（3）、式（4）中：n和m分別為隨機信號的不同時刻；E(x)為隨機信號的數(shù)學期望。

自相關函數(shù)或自協(xié)方差函數(shù)可用來檢測混有隨機噪聲的信號，1個系統(tǒng)可以通過系統(tǒng)輸入和輸出信號序列間的互相關函數(shù)最大值出現(xiàn)的位置來確定延遲。

2 功率譜

2.1 功率譜密度

功率譜表示隨機信號頻域的統(tǒng)計特性，有明顯的物理意義，1個信號的功率譜密度是其自相關函數(shù)的傅里葉變換[6]。

離散隨機信號的功率譜密度表示如式（5）。

兩離散隨機信號互相關函數(shù)的傅里葉變換定義為互功率譜密度函數(shù)如式（6）。

式（5）、式（6）中：m為隨機信號的時刻；ω為數(shù)字域頻率。

2.2 功率譜估計

各態(tài)歷經隨機信號的集合平均等于時間平均，因此從任何一個樣本即可得出隨機信號的全部信息。根據(jù)時間序列的一個有限觀察x(n)，(n=0，1，…，n-1)來估計功率譜密度，稱為功率譜估計。功率譜估計通常采用直接法和間接法2種方法[7]，直接法又稱周期圖法，它是把隨機序列x(n)的N個觀測數(shù)據(jù)視為1個能量有限的序列，直接計算x(n)的離散傅里葉變換x(k)，然和再取其幅值的平方，并除以N。

用直接法的功率譜估計，當數(shù)據(jù)長度太大時，譜曲線起伏加劇，若N太小，譜的分辨率不好，Bartlett法和Welch法是2種改進的直接法。

Bartlett法是將N點有限長序列x(n)分段求周期圖再平均。

Matlab代碼：

clear；

Fs=1000；

n=0:1/Fs:1；

xn=cos（2*pi*40*n） +3*cos（2*pi*100*n） +randn（size（n））；

nfft=1024；

window=boxcar（length（n））；

noverlap=0；

p=0.9；

[Pxx，Pxxc]=psd （xn，nfft，F(xiàn)s，window，noverlap，p）；

index=0:roud（nfft/2-1）；

plot_Pxx=10*log10（Pxx（index+1））；

plot_Pxxc=10*log10（Pxx（index+1））；

Welch法對Bartlett法進行了2個方面的修正，一是選擇適當?shù)拇昂瘮?shù)，二是在分段時使各段之間有重疊。

Matlab代碼：

clear；

Fs=1000；

n=0:1/Fs:1；

xn=cos（2*pi*40*n） +3*cos（2*pi*100*n） +randn（size（n））；

nfft=1024；

window=boxcar（100）；

window1=hamming（100）；

window2=blackman（100）；

noverlap=20；

range='half'；

[Pxx，f]=pwelch （xn， window， noverlap，nfft，F(xiàn)s，range）；

[Pxx1，f]=pwelch （xn，window1，noverlap，nfft，F(xiàn)s，range）；

[Pxx2，f]=pwelch （xn，window2，noverlap，nfft，F(xiàn)s，range）；

plot_Pxx=10*log10（Pxx）；

plot_Pxx1=10*log10（Pxx1）；

plot_Pxx2=10*log10（Pxx2）；

隨機信號通過離散時間線性非時變系統(tǒng)

當廣義平穩(wěn)隨機序列x(n)輸入到離散時間線性非時變系統(tǒng)時，系統(tǒng)輸出y(n)與系統(tǒng)的單位取樣響應h(n)存在關系如式（7）。

線性非時變系統(tǒng)輸入和輸出的互相關如式（8）。

式（7）、式（8）中，n和k分別為隨機信號的不同時刻。由此可知，輸入和輸出間的互相關函數(shù)是單位取樣響應的共軛序列和輸入自相關序列的卷積，兩邊取傅里葉變換如式（9）。

顯而易見，由自功率譜和互功率譜的測量可確定系統(tǒng)的頻率特性如式（10）。

式（9）、式（10）中，ω為數(shù)字域頻率。

3 功率譜估計在機械諧振分析算法中的應用實現(xiàn)

目前，功率譜估計已經應用在廣州數(shù)控設備有限公司某伺服驅動的機械諧振分析算法中。數(shù)控機床進給傳動系統(tǒng)的單軸驅動控制示意如圖1所示。

圖1 單軸控制示意圖

圖中，輸入信號是控制系統(tǒng)產生的轉矩指令（電壓或電流模擬量信號），激勵力由永磁同步電機產生。由聯(lián)軸器、工作平臺、滾珠絲杠副組成的機械傳動環(huán)節(jié)等效為1個雙慣量旋轉系統(tǒng)，該系統(tǒng)的輸入為電機扭矩與負載扭矩，輸出為電機角速度與負載角速度。電機和負載端的轉動慣量以及軸的剛度使雙慣量旋轉系統(tǒng)產生諧振頻率點（torsional natural frequency，TNF）和反諧振頻率點（anti-resonance frequency，ARF）[8]，當系統(tǒng)激勵力的頻率接近特定頻率時就引起機械諧振。因為系統(tǒng)的頻率特性可以通過自功率譜和互功率譜來確定，本文采用對系統(tǒng)進行持續(xù)地激勵來獲取系統(tǒng)在激勵過程中的輸入和輸出數(shù)據(jù)，并對輸入和輸出數(shù)據(jù)做功率譜估計，從而得到系統(tǒng)的頻率響應。

圖2所示為雙慣量旋轉系統(tǒng)典型的伯德圖（Bode diagram）。

圖2 雙慣量旋轉系統(tǒng)的伯德圖

3.1 激勵信號

本文提出的系統(tǒng)頻率特性測試采用掃頻法進行激勵。掃頻法是以頻率在待測頻率范圍內連續(xù)平穩(wěn)變化的等幅值的正弦信號作為激勵源的測試方式，常用的掃頻信號有白噪聲信號和Chirp信號。在實際應用中，一般采用偽隨機二進制序列（Pseudo-Random Binary Sequence，PRBS）代替白噪聲信號，PRBS是1種偽隨機的二位式周期序列，它的輸出結果只有2種電平，同時，在單個周期內數(shù)值可以看作是隨機變化的。常用的PRBS有M序列（最長線性反饋移位寄存器序列）、L序列（逆重復M序列）等[9]，M序列信號如圖3所示。

圖3 M序列信號波形

Chirp信號是調頻脈沖掃頻信號是雷達和通信領域經常使用的信號[10]，在頻率特性測試中，用到的是Chirp信號的實部，如式（11）。

式中：A為掃頻幅值； β為頻率變化速率； f0為起始頻率。Chirp信號基于余弦函數(shù)的變化規(guī)律，且瞬時頻率隨時間發(fā)生線性變化。通常，掃頻測試關注的頻率段是0～1 000 Hz，以驅動單元的伺服周期為61.25μs為例，在實際采樣中采樣頻率取8 kHz絕對滿足香農定理。

圖4 Chirp信號波形

掃頻信號在驅動單元的嵌入式芯片中實現(xiàn)，TMS320F28377是一款TI高性能TMS320C28x系列32位浮點單/雙核DSP處理器，它可以輕松通過多級移位寄存器的線性反饋產生M序列，也可以快速處理余弦和平方運算產生Chirp信號。

3.2 采樣數(shù)據(jù)分析

在分析掃頻信號時，輸入數(shù)據(jù)取轉矩指令信號，輸出數(shù)據(jù)取速度反饋信號。由式（10）可知，線性非時變系統(tǒng)的頻率特性可以通過自功率譜、互功率譜來確定，式（10）中的Sxx(ejω)是輸入轉矩信號的自功率譜密度函數(shù)，Syx(ejω)是速度反饋信號和輸入轉矩信號的互功率譜密度函數(shù)。

用C語言實現(xiàn)功率譜密度的計算：

intnfft；

if（winlength%2==1） winlength-=-1；

nfft=winlength；

intistart=0，step，iter，i，j；

double scale=0；

CComplexX_tmp；

CComplexY_tmp；

double*window= （double*）malloc（sizeof（dou?ble）*winlength）；

double*X_re= （double*） malloc（sizeof（dou?ble）*winlength）；

double*X_im= （double*） malloc（sizeof（dou?ble）*winlength）；

double*Y_re= （double*） malloc（sizeof（dou?ble）*winlength）；

double*Y_im= （double*） malloc（sizeof（dou?ble）*winlength）；

CComplex*cpsd= （CComplex*） malloc（sizeof（structCComplex_t） * （nfft/2+1））；

step=winlength-noverlap；

iter=1+（datalength-winlength）/step；

for （i=0；i＜winlength；i++）

{

window[i]=0.54-0.46*cos（2*PI*i/（win?length-1））；

scale+=window[i]*window[i]；

}

scale=scale*iter*fs；

for （i=0；i＜=nfft/2；i++）

{

cpsd[i].re=0.0；

cpsd[i].im=0.0；

}

for （i=0；i＜iter；i++）

{

f

or （j=0；j＜winlength；j++）

{

X_re[j]=x[j+istart]*window[j]；

Y_re[j]=y[j+istart]*window[j]；

X_im[j]=0.0；//虛部為0

Y_im[j]=0.0；//虛部為0

}

if （!Fft_transform （X_re， X_im， win?length）） exit（1）；

if （!Fft_transform （Y_re， Y_im， win?length）） exit（1）；

for （j=0；j＜=nfft/2；j++）

{

X_tmp.re=X_re[j]；

X_tmp.im=X_im[j]；

Y_tmp.re=Y_re[j]；

Y_tmp.im=Y_im[j]；

cpsd[j] =add （cpsd[j]， multiply（X_tmp，myconj（Y_tmp）））；

}

istart=istart+step；

}

for （i=0；i＜=nfft/2；i++）

{

cpsd[i]=multiplybynum（cpsd[i]，2.0/scale）；

}

cpsd[0]=multiplybynum（cpsd[0]，0.5）；

cpsd[nfft/2]=multiplybynum （cpsd[nfft/2]， 0.5）；

free（window）；

free（X_re）；

free（X_im）；

free（Y_re）；

free（Y_im）；

在C程序中實現(xiàn)計算轉矩指令信號的自功率譜密度以及速度反饋信號和轉矩指令信號之間的互功率譜密度，再通過式（10）計算出H(ejω)，并進一步計算出系統(tǒng)伯德圖的幅值和相位，如式（12）。

式中：H為系統(tǒng)的頻率特性函數(shù)，由此可以確定系統(tǒng)的反諧振頻率點和諧振頻率點。

3.3 掃頻測試結果

系統(tǒng)激勵采用Chirp信號，掃頻時間設置為10 s，初始頻率為1 Hz，終止頻率為1 kHz，采樣頻率為8 kHz。轉速由反饋的電機轉角位置計算而得，負載慣量盤的數(shù)量依次為：0、1、2、3、4、5。激勵信號的幅值依次為：500、712、935、1 024（該數(shù)值是模數(shù)轉化器電流輸入的數(shù)字量，其中1 024對應的電流值為電機額定轉矩6 A）。對每一種工況均重復4次掃頻測試，分析不同的負載慣量和掃頻幅值下系統(tǒng)的頻率特性差異。表1是分析每次的測試數(shù)據(jù)得到的諧振點位置信息，表2是分析每次的測試數(shù)據(jù)得到的反諧振點位置信息，每種工況重復測試4次。

觀察表1和表2的數(shù)據(jù)可以看出：系統(tǒng)諧振頻率隨負載的增加而變小，并且在相同工況下重復4次測試得到的系統(tǒng)諧振頻率點的一致性很好。

4 結束語

本文提出把數(shù)字信號處理中的功率譜估計運用在伺服系統(tǒng)的機械諧振分析中，依據(jù)各態(tài)歷經情況下離散隨機信號的集合平均等于時間平均的特性，通過功率譜估計算法分析系統(tǒng)的頻率響應更具有實驗價值。實驗結果表明：計算出來的反諧振頻率和諧振頻率反映了系統(tǒng)的頻率特性，為伺服單元速度環(huán)和位置環(huán)閉環(huán)帶寬的配置提供了依據(jù)，具有較好的工程應用價值和可行性。

表1 系統(tǒng)諧振頻率點

表2 系統(tǒng)反諧振頻率點