(1.北京交通大學(xué) 機(jī)械與電子控制工程學(xué)院，北京 100044； 2.北京宇航系統(tǒng)工程研究所，北京 100076；3.中國(guó)北方工業(yè)有限公司，北京 100053)

高速列車測(cè)速大多采用多普勒測(cè)速方式，其安裝在列車底部避免了傳統(tǒng)接觸式列車測(cè)速方式由于列車運(yùn)行時(shí)速過(guò)高帶來(lái)的軸承磨損、空轉(zhuǎn)等誤差[1-2]。在列車測(cè)速過(guò)程中由于運(yùn)行環(huán)境復(fù)雜，如路面變化、周邊列車影響，多普勒測(cè)速反射回波是由多個(gè)回波信號(hào)疊加而成的。一般將回波信號(hào)看作平穩(wěn)隨機(jī)信號(hào)，通過(guò)譜估計(jì)法對(duì)多普勒信號(hào)進(jìn)行頻域信號(hào)分析。譜估計(jì)分為經(jīng)典譜估計(jì)與現(xiàn)代譜估計(jì)兩種，以AR參數(shù)模型法為代表的現(xiàn)代譜估計(jì)對(duì)觀測(cè)數(shù)據(jù)進(jìn)行外推，預(yù)測(cè)求出觀測(cè)數(shù)據(jù)以外的其他數(shù)據(jù)，極大增大了頻譜的分辨率。在實(shí)際列車測(cè)速中，需要根據(jù)信號(hào)的有限采樣點(diǎn)來(lái)估計(jì)AR模型參數(shù)。常用方法有Yule-Walker法、協(xié)方差法、Burg算法和改進(jìn)協(xié)方差法[3-4]。其中改進(jìn)協(xié)方差法為Burg算法的進(jìn)一步發(fā)展。針對(duì)Burg算法的誤差來(lái)源，分析了基于窗函數(shù)的優(yōu)化算法，并針對(duì)改進(jìn)的協(xié)方差法實(shí)現(xiàn)了基于預(yù)測(cè)誤差功率最小意義的優(yōu)化算法。仿真實(shí)驗(yàn)證明，優(yōu)化算法對(duì)多普勒信號(hào)功率譜估計(jì)能力均有很大改善。

1 Burg原理及誤差分析

1.1 Burg原理

Burg 算法基于線性預(yù)測(cè)器的前、后向預(yù)測(cè)的總均方誤差之和最小的準(zhǔn)則直接從觀測(cè)數(shù)據(jù)來(lái)估計(jì)反射系數(shù)，然后通過(guò)Lenvinson-Durbin算法的遞推公式求出AR 模型參數(shù)[5]。這種算法對(duì)未知數(shù)據(jù)不需要做任何假設(shè)，估計(jì)精度較高。方法如下，前、后向預(yù)測(cè)誤差平均功率為

(1)

(2)

式中，km為反射系數(shù)即模型k階模型的第k個(gè)系數(shù)；m為模型階數(shù)，m=1,2，…，p。然后利用L-D遞推公式遞推出模型參數(shù):

(3)

式中，k=1,2，…，m-1。

am(m)=km

(4)

(5)

1.2 誤差分析

在Burg算法推導(dǎo)原理中可以看出。在L-D遞推公式中，高階系數(shù)由一階系數(shù)遞推。因此一階參數(shù)估計(jì)誤差對(duì)高階參數(shù)有一定影響。為分析一階參數(shù)影響，給定信號(hào)x(n)=sin(ωn+θ)，初相位為θ，歸一化頻率為ω。由Yule-Walker方程[6]可得一階模型參數(shù)a1(1)的真實(shí)值為a1(1)=-cosω。將信號(hào)公式代入反射系數(shù)km公式(式(2))中，求得a1(1)估計(jì)值為

(6)

式中，第二項(xiàng)為一階模型參數(shù)估計(jì)誤差，可以算出當(dāng)初相位θ為π/4的奇數(shù)倍；信號(hào)長(zhǎng)度N為1/4周期的奇數(shù)倍時(shí)誤差最大[7]。

2 兩種Burg優(yōu)化算法

2.1 基于窗函數(shù)的優(yōu)化算法

由式(3)～式(5)可知，模型參數(shù)直接由反射系數(shù)求出。為減小誤差可直接對(duì)反射系數(shù)km進(jìn)行修正。引人補(bǔ)償系數(shù)ωn(m)，其相當(dāng)于對(duì)信號(hào)進(jìn)行加窗處理，進(jìn)一步提高其方差能力，從而降低頻譜偏移程度、改善頻譜分辨率。Kaveh M等人基于平均頻率誤差方差最小原則得出最優(yōu)窗函數(shù)為[8]

(7)

(8)

2.2 基于預(yù)測(cè)誤差功率最小意義的優(yōu)化算法

由誤差分析可知，高階參數(shù)由低階參數(shù)遞推得到。低階參數(shù)誤差對(duì)高階參數(shù)推導(dǎo)有一定影響，導(dǎo)致譜估計(jì)質(zhì)量下降。針對(duì)改進(jìn)的協(xié)方差法原理中?ρm/?a(m,i)=0對(duì)原算法做如下改進(jìn)[10-11]：在預(yù)測(cè)誤差功率最小的意義下直接求解高階模型系數(shù)，為保證不增加過(guò)多計(jì)算量，可先求二階預(yù)測(cè)誤差功率。由先計(jì)算一階模型參數(shù)a1(1)，改為先計(jì)算二階參數(shù)a2(1)、a2(2)。如式(9)所示，先求二階誤差平均功率:

(9)

接下來(lái)計(jì)算二階模型系數(shù)a2(1)、a2(2)。

(10)

式中，各參數(shù)如下：

(11)

(12)

由式(12)可以看出，改進(jìn)算法通過(guò)二階參數(shù)的推導(dǎo)間接修正了一階參數(shù)估計(jì)誤差，避免了誤差項(xiàng)的影響，因此譜估計(jì)精度有了極大提高。

3 仿真分析與實(shí)驗(yàn)驗(yàn)證

3.1 算法性能對(duì)比

輸入仿真信號(hào)為

x(t)=sin(2π·50·t+π/4)+0.5randn(size(n))

(13)

給定信號(hào)頻率為50 Hz，初相位為π/4，采樣點(diǎn)數(shù)N為55，其相當(dāng)于1/4信號(hào)周期的11倍。由1.2節(jié)誤差分析可知，此時(shí)的數(shù)據(jù)長(zhǎng)度和初相位為最不利情況。信號(hào)混有均值為0，方差為1的隨機(jī)噪聲干擾。

取采樣頻率為1 kHz，模型階數(shù)m=7。對(duì)原Burg算法及兩種優(yōu)化算法進(jìn)行功率譜估計(jì)對(duì)比，結(jié)果如圖1所示。從圖1中可以看出，原算法頻譜分辨率最差與給定信號(hào)可達(dá)到5.66 Hz偏差，兩種優(yōu)化算法的頻譜分辨率均有一定改善，可有效抑制頻譜偏移。其中基于預(yù)測(cè)誤差功率最小意義的Burg優(yōu)化算法2的譜估計(jì)能力最好，與原信號(hào)只有0.29 Hz的頻率偏差。

圖1 不同算法的譜估計(jì)性能對(duì)比

3.2 算法復(fù)雜度分析

對(duì)原算法及兩種優(yōu)化算法進(jìn)行復(fù)雜度分析，當(dāng)模型階數(shù)m=7時(shí)，結(jié)果如表1所示。

表1 算法復(fù)雜度分析表

從表1中可以看出，優(yōu)化算法1中運(yùn)算次數(shù)較多，而優(yōu)化算法2只增加了2m-13次乘法，m-6次加法，對(duì)算法復(fù)雜度影響不大。3種算法運(yùn)行時(shí)間大致相同，說(shuō)明優(yōu)化算法在不增加運(yùn)行時(shí)間的前提下提高了頻譜的估計(jì)能力，對(duì)原算法起到了明顯優(yōu)化作用。

3.3 模型參數(shù)對(duì)算法性能影響

給定信號(hào)：

x(t)=5sin(2π·300·t)+3sin(2π·305·t)+0.5randn(size(n))

(14)

信號(hào)包含300 Hz和305 Hz頻率分量。為滿足采樣定理，取采樣頻率為1000 Hz。首先固定FFT點(diǎn)數(shù)為1024，給定不同采樣點(diǎn)數(shù)N。

由圖2中看出，隨著采樣點(diǎn)數(shù)N的增加，頻譜分辨率越來(lái)越好。但由于采樣點(diǎn)數(shù)增加會(huì)造成計(jì)算量的增大，不能滿足列車測(cè)速的實(shí)時(shí)性要求。綜合考慮頻譜計(jì)算精度與實(shí)時(shí)性要求，取采樣點(diǎn)數(shù)N為1024點(diǎn)。FFT計(jì)算點(diǎn)數(shù)參數(shù)選擇過(guò)程同理，取FFT點(diǎn)數(shù)為1024點(diǎn)。

圖2 不同采樣點(diǎn)數(shù)N對(duì)算法性能影響

3.4 列車測(cè)速范圍內(nèi)算法的譜估計(jì)驗(yàn)證

列車速度與多普勒頻率的計(jì)算公式為[8]

(15)

式中，v為列車速度(km/h)；fd為多普勒頻率；f0為雷達(dá)發(fā)射波頻率，為24.125 GHz；σ為雷達(dá)發(fā)射波方向與地面的夾角，假設(shè)為0°；c為光速3×108m/s。

由于高速列車被測(cè)速度范圍為0.2～600 km/h，由式(15)可知，多普勒頻率fd的范圍在8 Hz～27 kHz之間。

為驗(yàn)證優(yōu)化算法在列車測(cè)速過(guò)程中的適用性能，分別給定待測(cè)量頻率范圍內(nèi)不同頻段信號(hào)頻率組，三組分別為51 Hz、306 Hz、2006 Hz；4050 Hz、9000 Hz、14 kHz；22.2 kHz、24.9 kHz、26.3 kHz。綜合考慮頻率分辨率與系統(tǒng)實(shí)時(shí)性處理要求，對(duì)第一組頻率分量采用8 kHz采樣頻率，后兩組采用80 kHz采樣頻率。使用Matlab對(duì)算法譜估計(jì)能力進(jìn)行仿真，結(jié)果如圖3所示，圖中可以看出優(yōu)化算法在各個(gè)頻段內(nèi)適應(yīng)性良好，可準(zhǔn)確估計(jì)出各個(gè)頻段的信號(hào)頻率。

圖3 測(cè)速范圍內(nèi)不同頻段Burg算法譜估計(jì)

3.5 基于實(shí)測(cè)數(shù)據(jù)的的譜估計(jì)驗(yàn)證

為進(jìn)一步驗(yàn)證算法實(shí)用性，實(shí)物如圖4所示，通過(guò)信號(hào)發(fā)生器給定上述三組頻率，通過(guò)DSP進(jìn)行信號(hào)處理測(cè)試算法譜估計(jì)能力。測(cè)試結(jié)果如表2所示，表中可以看出，算法可以準(zhǔn)確識(shí)別各個(gè)頻段頻率且誤差小于1%。

圖4 DSP信號(hào)處理實(shí)物圖

理論頻率/Hz實(shí)測(cè)頻率/Hz理論速度/km·h-1實(shí)測(cè)速度/km·h-1頻率誤差/%5150.78 1.14 1.13 0.43 306308.59 6.90 6.84 0.85 20062007.81 44.91 44.87 0.09 40504062.50 90.86 90.58 0.31 90008984.38 200.95 201.30 0.17 1400013984.38 312.78 313.13 0.112220022187.50 496.25 496.53 0.062490024921.88 557.41 556.92 0.092630026328.13 588.86 588.24 0.11

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

針對(duì)多普勒信號(hào)處理的Burg算法進(jìn)行了研究，分析其推導(dǎo)原理及產(chǎn)生誤差原因。對(duì)于模型低階參數(shù)估計(jì)誤差來(lái)源進(jìn)行了兩種修正優(yōu)化算法分析。一種是基于窗函數(shù)的優(yōu)化，另一種基于預(yù)測(cè)誤差功率最小意義的優(yōu)化。通過(guò)Matlab仿真分析對(duì)比了原算法與兩優(yōu)化算法的功率譜估計(jì)能力，選擇出模型最優(yōu)參數(shù),并對(duì)列車測(cè)速范圍內(nèi)進(jìn)行了優(yōu)化算法的譜估計(jì)驗(yàn)證。實(shí)驗(yàn)結(jié)果表明，兩種優(yōu)化算法在不增加運(yùn)行時(shí)間的基礎(chǔ)上對(duì)頻譜估計(jì)性能均有極大改善，可以識(shí)別出列車測(cè)速范圍內(nèi)各個(gè)頻段頻率，譜估計(jì)頻率誤差小于1%。