常系數(shù)αβ濾波器參數(shù)的計(jì)算機(jī)優(yōu)化設(shè)計(jì)研究*

侯建

(北京電子工程總體研究所，北京 100854)

0 引言

αβ濾波器是一種十分常見且應(yīng)用廣泛的濾波器，在雷達(dá)數(shù)據(jù)處理中大量應(yīng)用αβ濾波器，因此針對(duì)具體的應(yīng)用需求進(jìn)行αβ濾波器的參數(shù)優(yōu)化設(shè)計(jì)是經(jīng)常遇到的問題。本文選取了目前廣泛使用的2個(gè)優(yōu)化設(shè)計(jì)準(zhǔn)則，建立了優(yōu)化設(shè)計(jì)模型，用Matlab編寫了優(yōu)化設(shè)計(jì)程序。

1 αβ濾波器簡介

為了便于后面的論述，在這里先對(duì)αβ濾波器的有關(guān)基本概念進(jìn)行一個(gè)簡要的梳理。

αβ濾波器的濾波算法如下[1]：

(1)

(2)

(3)

根據(jù)αβ濾波器的濾波方程，可以推出其傳遞函數(shù)為[1]

(4)

(5)

式中：Hf(z)，Hp(z)分別為αβ濾波器的濾波估計(jì)和預(yù)測估計(jì)的傳遞函數(shù)。

可以求出αβ濾波器濾波估計(jì)的穩(wěn)態(tài)方差壓縮系數(shù)[1]：

(6)

可以推出αβ濾波器的等效噪聲帶寬為[2]

(7)

由傳遞函數(shù)可以看出αβ濾波器是一個(gè)二階離散系統(tǒng)。

對(duì)于一個(gè)穩(wěn)定的二階離散系統(tǒng)，如果其階躍響應(yīng)為衰減振蕩或等幅振蕩(具有一對(duì)共軛極點(diǎn)或一個(gè)二重極點(diǎn))，則傳遞函數(shù)可以寫成:

(8)

系統(tǒng)傳遞函數(shù)的極點(diǎn)為

z1=e-ξω0Tejωd T,

(9)

z2=e-ξω0Te-jωd T.

(10)

對(duì)比式(4),(5)和(8)，可以得到：

1-α=e-2ξω0T,

(11)

2-α-β=2e-ξω0TcosωdT,

(12)

(13)

根據(jù)式(11)和(12)可以得到：

(14)

(15)

2 αβ濾波器參數(shù)的計(jì)算機(jī)優(yōu)化設(shè)計(jì) 對(duì)于常系數(shù)αβ濾波器，需要確定的參數(shù)有3個(gè)：α，β和T。

通常濾波周期T是由系統(tǒng)其他諸多因素確定的，因此這里僅討論在濾波周期已經(jīng)確定的情況下參數(shù)α和β的優(yōu)化設(shè)計(jì)問題。

從一般的意義上講，如果有n個(gè)參數(shù)需要確定，就必須有n個(gè)彼此獨(dú)立的方程，求解這n個(gè)方程構(gòu)成的方程組就可以得到n個(gè)參數(shù)。因此αβ濾波器的參數(shù)設(shè)計(jì)問題就轉(zhuǎn)化為如何確定2個(gè)相互獨(dú)立的方程。

在前面的論述中已經(jīng)給出了有關(guān)αβ濾波器的一些公式，如：濾波估計(jì)的穩(wěn)態(tài)方差壓縮系數(shù)公式(6)、濾波器的等效噪聲帶寬公式(7)、系統(tǒng)的阻尼系數(shù)和阻尼振蕩頻率公式(13)～(15)等?，F(xiàn)在討論如何利用這些公式并結(jié)合一些準(zhǔn)則來確定參數(shù)α和β。

αβ濾波器一般是用來對(duì)被測量進(jìn)行跟蹤濾波的，因此跟蹤濾波精度是一個(gè)十分重要的考量因素，而穩(wěn)態(tài)方差則是衡量穩(wěn)態(tài)跟蹤精度的一個(gè)重要指標(biāo)，穩(wěn)態(tài)方差壓縮系數(shù)表征了穩(wěn)態(tài)情況下濾波器輸出的方差相對(duì)輸入方差的比例系數(shù)，該數(shù)值越小，則輸出方差就越小，因此通常情況下可以將式(6)作為一個(gè)方程；當(dāng)然在某些情況下如果系統(tǒng)對(duì)等效噪聲帶寬有明確的要求也可以將式(7)作為一個(gè)方程。實(shí)際上式(6)和式(7)在很大程度上是等效的，因?yàn)槿绻麕捲綄?，則系統(tǒng)對(duì)起伏噪聲的濾波效果就越差，因此輸出方差就越大，穩(wěn)態(tài)方差壓縮系數(shù)就越大。

除了穩(wěn)態(tài)跟蹤精度之外，通常也關(guān)心系統(tǒng)的暫態(tài)特性，如一般希望系統(tǒng)的過渡過程盡可能短，超調(diào)盡可能小。前文已說明αβ濾波器是一個(gè)二階系統(tǒng)，對(duì)于一個(gè)二階系統(tǒng)，系統(tǒng)的阻尼系數(shù)ξ和自然諧振頻率ω0是與系統(tǒng)暫態(tài)特性密切相關(guān)的參數(shù)；圖1給出了阻尼系數(shù)ξ取不同數(shù)值時(shí)二階系統(tǒng)的階躍響應(yīng)。

式(16～18)給出了0<ξ<1時(shí)二階系統(tǒng)階躍響應(yīng)暫態(tài)特性的計(jì)算公式[3]：

(16)

(17)

(18)

式中：Δ為二階系統(tǒng)階躍響應(yīng)的超調(diào)；Tr為二階系統(tǒng)階躍響應(yīng)的上升時(shí)間；Ts為二階系統(tǒng)階躍響應(yīng)的調(diào)整時(shí)間；ε為二階系統(tǒng)階躍響應(yīng)的調(diào)整時(shí)間對(duì)應(yīng)的歸一化誤差(階躍響應(yīng)在Ts時(shí)刻的取值與穩(wěn)態(tài)值之差相對(duì)穩(wěn)態(tài)值歸一化)。

從圖1和式(16)～(18)可以看出：

(1) 當(dāng)阻尼系數(shù)ξ=0時(shí)，系統(tǒng)的階躍響應(yīng)是一個(gè)等幅振蕩，系統(tǒng)處于臨界穩(wěn)定狀態(tài)，不能正常工作，因此要求ξ>0；

(2) 當(dāng)0<ξ<1時(shí)，系統(tǒng)的階躍響應(yīng)是一個(gè)衰減振蕩，阻尼系數(shù)ξ越大，系統(tǒng)的阻尼越大，衰減越快，超調(diào)也越??；當(dāng)ξ一定時(shí)，要使上升時(shí)間縮短，則需要提高系統(tǒng)的自然諧振頻率ω0；

(3) 當(dāng)ξ≥1時(shí)，系統(tǒng)的階躍響應(yīng)是一個(gè)單調(diào)衰減(沒有振蕩)，此時(shí)系統(tǒng)沒有超調(diào)，阻尼系數(shù)ξ越大，階躍響應(yīng)的上升越慢，調(diào)整時(shí)間也越長；雖然ξ≥1時(shí)，系統(tǒng)的階躍響應(yīng)沒有振蕩(也沒有超調(diào))，但此時(shí)系統(tǒng)的階躍響應(yīng)相對(duì)穩(wěn)態(tài)值的收斂速度要比0<ξ<1時(shí)慢。

通過上述分析可以看出如果從系統(tǒng)的暫態(tài)特性考慮，一般取0.7≤ξ≤1比較合適。因此確定參數(shù)α和β的第2個(gè)公式可以按下述原則確定：

(1) 取ξ=1，即系統(tǒng)設(shè)計(jì)工作在臨界阻尼狀態(tài)；

(2) 在0.7≤ξ≤1范圍內(nèi)選取一個(gè)合適的阻尼系數(shù)。

下面分別討論上述2種情況。

當(dāng)取ξ=1時(shí)，此時(shí)系統(tǒng)的特征方程有一對(duì)重根，這時(shí)可以推出參數(shù)α和β滿足[1]：

(19)

當(dāng)在0.7≤ξ≤1范圍內(nèi)選取一個(gè)合適的阻尼系數(shù)ξ時(shí)，可以推出參數(shù)α和β滿足：

(20)

顯然，式(19)是式(20)在ξ=1時(shí)的特例。

此外，文獻(xiàn)[1]還給出了一種優(yōu)化設(shè)計(jì)情況，即在“給定速度估值暫態(tài)平方差條件下使速度估值穩(wěn)態(tài)方差壓縮系數(shù)最小”。根據(jù)該準(zhǔn)則可以推出參數(shù)α和β滿足：

(21)

至此給出了確定參數(shù)α和β可以采用的公式。為便于描述，稱式(19)～(21)為優(yōu)化準(zhǔn)則公式，稱式(6)，(7)為設(shè)計(jì)指標(biāo)公式。下面就可以利用所選用的公式，根據(jù)具體的設(shè)計(jì)要求(如要求的“穩(wěn)態(tài)方差壓縮系數(shù)/等效噪聲帶寬”、“阻尼系數(shù)——可以根據(jù)希望的系統(tǒng)暫態(tài)響應(yīng)指標(biāo)確定”)求解方程得到參數(shù)α和β。具體方法如下：

圖1 二階系統(tǒng)階躍響應(yīng)Fig.1 Step Response of second-order system

首先從式(6)或式(7)中選擇一個(gè)作為設(shè)計(jì)指標(biāo)公式并確定希望的指標(biāo)設(shè)計(jì)值(如選擇式(6)希望的指標(biāo)設(shè)計(jì)值就是方差壓縮系數(shù)，如選擇式(7)希望的指標(biāo)設(shè)計(jì)值就是等效噪聲帶寬)；然后從式(19)～(21)中選擇一個(gè)作為優(yōu)化準(zhǔn)則公式(當(dāng)選擇式(20)時(shí)，需要同時(shí)確定一個(gè)合適的阻尼系數(shù)值ξ)，這樣就構(gòu)成了一個(gè)由2個(gè)方程組成的二元方程組，求解這個(gè)方程組就得到了所要求的參數(shù)α和β。

需要說明的是：如果選擇了式(6)或式(7)再加上式(21)則可以得到一個(gè)關(guān)于α的二次方程，可以直接求解。但如果選擇了式(6)或式(7)再加上式(19)或式(20)，則得到的是一個(gè)關(guān)于α的超越方程，只能采用數(shù)值方法求解。在進(jìn)行數(shù)值解時(shí)需要首先確定一個(gè)初值，然后利用程序進(jìn)行迭代求解。在這種情況下如果初值選取的不合適則有可能導(dǎo)致迭代運(yùn)算不能收斂或收斂到錯(cuò)誤的結(jié)果。因此，如何選定初值關(guān)系到迭代運(yùn)算是否能夠正確地收斂。下面給出確定初值的方法。

首先需要確定α和β的取值范圍。

(4-2α-β)s2+2αs+β=0.

(22)

則根據(jù)系統(tǒng)穩(wěn)定性判定準(zhǔn)則，必須滿足以下條件：

(23)

圖2 αβ濾波器的參數(shù)取值范圍Fig.2 Effective parameter range of αβ filter

有了α和β的取值范圍和已經(jīng)選定的公式所確定的α和β的約束關(guān)系就可以按照下述方法確定α的計(jì)算初值，并進(jìn)一步經(jīng)迭代運(yùn)算得到最終設(shè)計(jì)值：

1) 在[0,1]范圍內(nèi)按照一定的間隔步長選取一系列的α值；

2) 對(duì)于每一個(gè)α，利用已經(jīng)選定的公式所確定的α和β的約束關(guān)系計(jì)算對(duì)應(yīng)的β；

3) 根據(jù)前面兩步獲得的每一對(duì)α和β，利用公式(6),(7),(13)～(15)計(jì)算αβ濾波器的相關(guān)性能指標(biāo)；

4) 取最接近所要求指標(biāo)所對(duì)應(yīng)的α值作為計(jì)算初值；

5) 利用所確定α計(jì)算初值啟動(dòng)迭代運(yùn)算，得到最終的α設(shè)計(jì)值；

6) 利用已經(jīng)選定的公式所確定的α和β的約束關(guān)系計(jì)算對(duì)應(yīng)的β得到β的設(shè)計(jì)值。

圖3和圖4了給出了αβ濾波器的參數(shù)計(jì)算機(jī)優(yōu)化設(shè)計(jì)流程。關(guān)于參數(shù)α的迭代算法，可以根據(jù)情況選用各種合適的數(shù)值解法，本文選擇了最簡單的牛頓切線法。

圖3 αβ濾波器的參數(shù)計(jì)算機(jī)優(yōu)化設(shè)計(jì)流程Fig.3 Optimized design flat for alpha-beita filter parameter

3 設(shè)計(jì)實(shí)例

根據(jù)本文所介紹的方法，我們按照?qǐng)D3和圖4的流程用Matlab編寫了αβ濾波器的參數(shù)計(jì)算機(jī)優(yōu)化設(shè)計(jì)程序，這里給出2個(gè)設(shè)計(jì)實(shí)例。

實(shí)例1：假定要設(shè)計(jì)的αβ濾波器的采樣周期為T=36 ms，要求的穩(wěn)態(tài)方差壓縮系數(shù)為0.36，采用臨界阻尼準(zhǔn)則。

設(shè)計(jì)結(jié)果為α=0.485 945，β=0.080 106，將參數(shù)α和β的設(shè)計(jì)結(jié)果代入式(6),(15)和(7)分別得到：

穩(wěn)態(tài)方差壓縮系數(shù)：

阻尼系數(shù)：

等效噪聲帶寬：

可見設(shè)計(jì)結(jié)果滿足設(shè)計(jì)要求(穩(wěn)態(tài)方差壓縮系數(shù)為0.36,阻尼系數(shù)為1——臨界阻尼)。

圖5～8給出了一些驗(yàn)證曲線。

圖5為按照臨界阻尼準(zhǔn)則，當(dāng)參數(shù)α在[0,1]區(qū)間變化時(shí)，參數(shù)β的變化規(guī)律曲線，圖中的“*”點(diǎn)為本次設(shè)計(jì)結(jié)果(對(duì)應(yīng)穩(wěn)態(tài)方差壓縮系數(shù)為0.36且臨界阻尼的參數(shù)α和β)；

圖6為按照臨界阻尼準(zhǔn)則，當(dāng)參數(shù)α在[0,1]區(qū)間變化時(shí)，穩(wěn)態(tài)方差壓縮系數(shù)為的變化規(guī)律曲線，圖中的“*”點(diǎn)為本次設(shè)計(jì)結(jié)果(對(duì)應(yīng)穩(wěn)態(tài)方差壓縮系數(shù)為0.36且臨界阻尼的參數(shù)α和β)；

圖4 αβ濾波器的參數(shù)計(jì)算機(jī)優(yōu)化設(shè)計(jì)流程(續(xù))Fig.4 Optimized design flat for alpha-beita filter parameter(continue)

圖7為按照臨界阻尼準(zhǔn)則，當(dāng)參數(shù)α在[0,1]區(qū)間變化時(shí)，阻尼系數(shù)的變化規(guī)律曲線(曲線為常數(shù)1表明滿足臨界阻尼準(zhǔn)則)；

圖5 按照臨界阻尼準(zhǔn)則，當(dāng)α在[0,1]變化時(shí)， β的變化曲線Fig.5 Value distribution curve of parameter beita according to the rule of critical damping while the parameter alpha varies in range of [0,1]

圖6 按照臨界阻尼準(zhǔn)則，當(dāng)α在[0,1]變化時(shí)，穩(wěn)態(tài)方差壓縮系數(shù)的變化曲線Fig.6 Value distribution curve of variance compress coefficient according to the rule of critical damping while the parameter alpha varies in range of [0,1]

圖8為求解參數(shù)α?xí)r的超越方程f(x)=0對(duì)應(yīng)的函數(shù)y=f(x)曲線。可見：由于y=f(x)曲線的一階導(dǎo)數(shù)符號(hào)不是恒定的，因此在用牛頓切線法求解方程時(shí)對(duì)迭代運(yùn)算的初值必須正確選取，否則有可能得不到正確的收斂結(jié)果。

圖7 按照臨界阻尼準(zhǔn)則，當(dāng)α在[0,1]變化時(shí)， 阻尼系數(shù)的變化曲線Fig.7 Value distribution curve of damping coefficient according to the rule of critical damping while the parameter alpha varies in range of [0,1]

圖8 求解參數(shù)α?xí)r的超越方程f(x)=0的 函數(shù)y=f(x)曲線Fig.8 Value curve of function y=f(x) where f(x) is the transcendental equation f(x)=0 to solve to determine the parameter alpha

圖9為所設(shè)計(jì)的αβ濾波器對(duì)給定輸入的濾波輸出。圖中曲線“in”，“p”，“f”分別表示輸入、預(yù)測輸出、濾波輸出。輸入信號(hào)是一個(gè)附加了噪聲的勻加速變化的信號(hào)。從圖中可以看出所設(shè)計(jì)的αβ濾波器的預(yù)測和濾波輸出能夠跟隨輸入變化，并且從濾波器的輸入噪聲和輸出起伏誤差的比較可以看出預(yù)測和濾波輸出的起伏要比輸入有一定程度的抑制，經(jīng)過對(duì)輸入和輸出的統(tǒng)計(jì)，輸入中加入的噪聲方差為1.105 547，預(yù)報(bào)輸出的起伏方差為0.604 702，濾波輸出的起伏方差為0.448 041；濾波輸出的方差壓縮倍數(shù)為0.405 266(設(shè)計(jì)值為0.36)。

實(shí)例2：假定要設(shè)計(jì)的αβ濾波器的采樣周期為T=36 ms，要求的穩(wěn)態(tài)方差壓縮系數(shù)為0.36，采用給定阻尼準(zhǔn)則，要求阻尼系數(shù)為0.78。

圖9 所設(shè)計(jì)濾波器的輸入和輸出Fig.9 Input and output of the filter designed

設(shè)計(jì)結(jié)果為α=0.457 137，β=0.113 303，將參數(shù)α和β的設(shè)計(jì)結(jié)果代入式(6),(15)和(7)分別得到：

穩(wěn)態(tài)方差壓縮系數(shù)：

阻尼系數(shù)：

等效噪聲帶寬：

可見設(shè)計(jì)結(jié)果滿足設(shè)計(jì)要求(穩(wěn)態(tài)方差壓縮系數(shù)為0.36,阻尼系數(shù)為0.78)。

圖10～13給出了一些驗(yàn)證曲線。

圖10為按照給定阻尼系數(shù)為0.78作為設(shè)計(jì)準(zhǔn)則，當(dāng)參數(shù)α在[0,1]區(qū)間變化時(shí)，參數(shù)β的變化規(guī)律曲線，圖中的“*”點(diǎn)為本次設(shè)計(jì)結(jié)果(對(duì)應(yīng)穩(wěn)態(tài)方差壓縮系數(shù)為0.36且阻尼系數(shù)為0.78的參數(shù)α和β)；

圖11為按照給定阻尼系數(shù)為0.78作為設(shè)計(jì)準(zhǔn)則，當(dāng)參數(shù)α在[0,1]區(qū)間變化時(shí)，穩(wěn)態(tài)方差壓縮系數(shù)為的變化規(guī)律曲線，圖中的“*”點(diǎn)為本次設(shè)計(jì)結(jié)果(對(duì)應(yīng)穩(wěn)態(tài)方差壓縮系數(shù)為0.36且阻尼系數(shù)為0.78的參數(shù)α和β)；

圖12為按照給定阻尼系數(shù)為0.78作為設(shè)計(jì)準(zhǔn)則，當(dāng)參數(shù)α在[0,1]區(qū)間變化時(shí)，阻尼系數(shù)的變化規(guī)律曲線(曲線為常數(shù)0.78表明滿足給定阻尼準(zhǔn)則)；

圖13為求解參數(shù)α?xí)r的超越方程f(x)=0對(duì)應(yīng)的函數(shù)y=f(x)曲線?？梢姡河捎趛=f(x)曲線的一階導(dǎo)數(shù)符號(hào)不是恒定的，因此在用牛頓切線法求解方程時(shí)對(duì)迭代運(yùn)算的初值必須正確選取，否則有可能得不到正確的收斂結(jié)果。

圖10 按照給定阻尼0.78準(zhǔn)則，當(dāng)α在[0,1]變化時(shí)， β的變化曲線Fig.10 Value distribution curve of parameter beita according to the rule of given damping with value of 0.78 while the parameter alpha varies in range of [0,1]

圖11 按照給定阻尼0.78準(zhǔn)則，α在[0,1]區(qū)間變化時(shí)，穩(wěn)態(tài)方差壓縮系數(shù)的變化曲線Fig.11 Value distribution curve of variance compress coefficient according to the rule of given damping with value of 0.78 while the parameter alpha varies in range of [0,1]

圖12 按照給定阻尼0.78準(zhǔn)則，當(dāng)α在[0,1] 區(qū)間變化時(shí)，阻尼系數(shù)的變化曲線Fig.12 Value distribution curve of damping coefficient according to the rule of given damping with value of 0.78 while the parameter alpha varies in range of [0,1]

圖13 求解參數(shù)α?xí)r的超越方程f(x)=0的函數(shù)y=f(x)曲線Fig.13 Value curve of function y=f(x) where f(x) is the transcendental equation f(x)=0 to solve to determine the parameter alpha

4 結(jié)論

本文給出了3個(gè)優(yōu)化設(shè)計(jì)準(zhǔn)則:

準(zhǔn)則1:臨界阻尼(ξ=1);

準(zhǔn)則2:給定阻尼(0<ξ<1);

準(zhǔn)則3:給定速度估值暫態(tài)平方差條件下使速度估值穩(wěn)態(tài)方差壓縮系數(shù)最小。

準(zhǔn)則1和準(zhǔn)則2都已經(jīng)過實(shí)例驗(yàn)證，得到了預(yù)期效果，表明了設(shè)計(jì)模型的正確性。由于采用準(zhǔn)則3時(shí)得到一個(gè)關(guān)于α的二次方程，可以直接求解，比較簡單，所以本文沒有給出設(shè)計(jì)算例。

利用本文所建立的模型可以實(shí)現(xiàn)對(duì)αβ濾波器的參數(shù)優(yōu)化設(shè)計(jì)。