朱琰虹，文光武

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

0 引言

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

1 隨機信號

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

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

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

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

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

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

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

2 功率譜

2.1 功率譜密度

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

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

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

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

2.2 功率譜估計

各態(tài)歷經(jīng)隨機信號的集合平均等于時間平均，因此從任何一個樣本即可得出隨機信號的全部信息。根據(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)的單位取樣響應(yīng)h(n)存在關(guān)系如式（7）。

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

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

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

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

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

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

圖1 單軸控制示意圖

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

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

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

3.1 激勵信號

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

圖3 M序列信號波形

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

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

圖4 Chirp信號波形

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

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

在分析掃頻信號時，輸入數(shù)據(jù)取轉(zhuǎn)矩指令信號，輸出數(shù)據(jù)取速度反饋信號。由式（10）可知，線性非時變系統(tǒng)的頻率特性可以通過自功率譜、互功率譜來確定，式（10）中的Sxx(ejω)是輸入轉(zhuǎn)矩信號的自功率譜密度函數(shù)，Syx(ejω)是速度反饋信號和輸入轉(zhuǎn)矩信號的互功率譜密度函數(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)計算轉(zhuǎn)矩指令信號的自功率譜密度以及速度反饋信號和轉(zhuǎn)矩指令信號之間的互功率譜密度，再通過式（10）計算出H(ejω)，并進一步計算出系統(tǒng)伯德圖的幅值和相位，如式（12）。

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

3.3 掃頻測試結(jié)果

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

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

4 結(jié)束語

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

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

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