朱琰虹，文光武

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

0 引言

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

1 隨機(jī)信號(hào)

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

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

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

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

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

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

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

2 功率譜

2.1 功率譜密度

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

離散隨機(jī)信號(hào)的功率譜密度表示如式（5）。

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

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

2.2 功率譜估計(jì)

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

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

Bartlett法是將N點(diǎn)有限長(zhǎng)序列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法對(duì)Bartlett法進(jìn)行了2個(gè)方面的修正，一是選擇適當(dāng)?shù)拇昂瘮?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）；

隨機(jī)信號(hào)通過(guò)離散時(shí)間線性非時(shí)變系統(tǒng)

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

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

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

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

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

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

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

圖1 單軸控制示意圖

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

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

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

3.1 激勵(lì)信號(hào)

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

圖3 M序列信號(hào)波形

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

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

圖4 Chirp信號(hào)波形

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

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

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

用C語(yǔ)言實(shí)現(xiàn)功率譜密度的計(jì)算：

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

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

3.3 掃頻測(cè)試結(jié)果

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

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

4 結(jié)束語(yǔ)

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

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

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