基于功率譜估計的伺服系統機械諧振分析

朱琰虹，文光武

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

0 引言

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

1 隨機信號

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

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

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

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

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

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

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

2 功率譜

2.1 功率譜密度

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

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

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

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

2.2 功率譜估計

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

用直接法的功率譜估計，當數據長度太大時，譜曲線起伏加劇，若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，Fs，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個方面的修正，一是選擇適當的窗函數，二是在分段時使各段之間有重疊。

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，Fs，range）；

[Pxx1，f]=pwelch （xn，window1，noverlap，nfft，Fs，range）；

[Pxx2，f]=pwelch （xn，window2，noverlap，nfft，Fs，range）；

plot_Pxx=10*log10（Pxx）；

plot_Pxx1=10*log10（Pxx1）；

plot_Pxx2=10*log10（Pxx2）；

隨機信號通過離散時間線性非時變系統

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

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

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

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

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

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

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

圖1 單軸控制示意圖

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

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

圖2 雙慣量旋轉系統的伯德圖

3.1 激勵信號

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

圖3 M序列信號波形

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

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

圖4 Chirp信號波形

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

3.2 采樣數據分析

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

用C語言實現功率譜密度的計算：

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程序中實現計算轉矩指令信號的自功率譜密度以及速度反饋信號和轉矩指令信號之間的互功率譜密度，再通過式（10）計算出H(ejω)，并進一步計算出系統伯德圖的幅值和相位，如式（12）。

式中：H為系統的頻率特性函數，由此可以確定系統的反諧振頻率點和諧振頻率點。

3.3 掃頻測試結果

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

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

4 結束語

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

表1 系統諧振頻率點

表2 系統反諧振頻率點