幾乎周期時(shí)變MA系統(tǒng)的辨識(shí)
本文利用循環(huán)統(tǒng)計(jì)量討論幾乎周期時(shí)變滑動(dòng)平均(APTV-MA)系統(tǒng)的辨識(shí)問題.提出了基于奇數(shù)階時(shí)變累量的參數(shù)估計(jì)的法方程方法,并導(dǎo)出了基于時(shí)變累量和循環(huán)累量的系統(tǒng)模型定階方法.仿真實(shí)驗(yàn)說明了本文所提方法的性能.
關(guān)鍵詞:循環(huán)統(tǒng)計(jì)量;幾乎周期時(shí)變MA系統(tǒng);法方程;系統(tǒng)辨識(shí)
Identification of Almost Periodic Time-Varying MA Systems
LI Hong-wei
(Department of Mathermatics and Physics,China University of Geoscience,Wuhan 430074,China)
YUAN Bao-zong
(Institute of Information Science,Northern Jiaotong University,Beijng 100044,China)
Abstract:This paper addresses the problem of identification of almost periodic time-varying moving average (APTV-MA) systems using cyclic statistics.The normal equation approaches based on odd order time-varying cumulants are presented for parameters estimation.Model order detemination procedures based on time-varying and cyclic cumulants are also derived.Simulations are performed to illustrate performance of the presented methods.
Key words:cyclic statistics;almost periodic time-varying MA system;normal equation;identification of system
一、引 言
考慮q階幾乎周期時(shí)變滑動(dòng)平均(APTV-MA)信號(hào)模型
(1)
并且,假設(shè)下列假定成立.
假定1 激勵(lì)噪聲?(t)是零均值獨(dú)立同分布的,且存在一個(gè)奇數(shù)k(3),使得其k階累量γk?滿足0<|γk∈|<+∞.
假定2 對(duì)于每一個(gè)固定的n,有限沖激響應(yīng)序列h(t,n)是關(guān)于t的幾乎周期實(shí)序列(可以是非最小相位的),并且,h(0,0)=1;h(t,0)≠0,h(t,q)≠0,t.
在實(shí)際應(yīng)用中,信號(hào)x(t)是在噪聲中被觀測(cè)的:
y(t)=x(t)+v(t),t=0,1,…,T-1 (2)
其中,加性噪聲v(t)滿足
假定3 v(t)是一個(gè)與∈(t)獨(dú)立的功率譜未知的零均值高斯(有色)過程(可以是非平穩(wěn)的).
在模型(2)下,y(t)是非平穩(wěn)的,原有的討論時(shí)不變MA系統(tǒng)的辨識(shí)方法不能直接應(yīng)用.易證y(t)是循環(huán)平穩(wěn)的.因此,我們的目的是根據(jù)觀測(cè)值{y(t),t=0,1,…,T-1}利用循環(huán)統(tǒng)計(jì)量來估計(jì)
上述問題可以看作時(shí)不變MA系統(tǒng)在幾乎周期時(shí)變情形的推廣.這一類問題在通訊,水文氣象等領(lǐng)域具有廣泛的應(yīng)用[5,11,14].在高斯信號(hào)情形,文獻(xiàn)[15]提出了基于二階統(tǒng)計(jì)量的近似最大似然技術(shù)估計(jì)信號(hào)參數(shù).相關(guān)問題的討論可見文獻(xiàn)[8]和[13].最近[5],Dandawate和Giannakis對(duì)于幾乎周期時(shí)變MA信號(hào)進(jìn)行了較深入的研究,他們利用高階循環(huán)統(tǒng)計(jì)量來辨識(shí)APTV-MA系統(tǒng),提出了模型參數(shù)估計(jì)的線性和非線性算法以及一種統(tǒng)計(jì)的定階方法.
眾所周知,對(duì)于時(shí)不變MA系統(tǒng),基于高階累量的參數(shù)估計(jì)方法可分為閉式解,法方程和非線性最優(yōu)化解三大類[12].由于不涉及極值問題求解,前二種線性方法尤其令人重視[7,17].但是到目前為止,對(duì)于幾乎周期時(shí)變MA系統(tǒng)的辨識(shí),這兩種方法的研究還較小.本文討論利用循環(huán)統(tǒng)計(jì)量辨識(shí)APTV-MA系統(tǒng)的新方法,提出了參數(shù)估計(jì)的法方程方法和基于時(shí)變累量和循環(huán)累量的時(shí)域定階方法.
二、參數(shù)估計(jì)方法
在本節(jié)先假定幾乎周期時(shí)變MA系統(tǒng)的階q已知,推導(dǎo)APTV-MA信號(hào)參數(shù)的法方程.
實(shí)信號(hào)z(t)的k階時(shí)變的和循環(huán)的累量分別定義[3,5]為
其中,τ=(τ1,…,τk-1).有關(guān)循環(huán)統(tǒng)計(jì)量的詳細(xì)定義和記號(hào)參見文獻(xiàn)[3].
當(dāng)k3時(shí),由假定3和累量性質(zhì)有ckx(t;τ)=cky(t;τ).即觀測(cè)信號(hào)y(t)與無加性噪聲污染的輸出信號(hào)x(t)具有相同的時(shí)變累量,因此,在基于高階時(shí)變和循環(huán)累量的分析中不必考慮加性高斯噪聲的影響.
由累量性質(zhì)可知,滿足模型(1)和假定1~2的APTV-MA信號(hào)x(t)的k(>3)階時(shí)變累量為
(3)
它是時(shí)不變MA系統(tǒng)的BBR公式[1]在APTV-MA情形的推廣.
令τ1=τ2=…=τk-2=q,則由式(3)得到
ckx(t;q,…,q,τk-1)=γk∈h(t,0)hk-2(t+q,q).h(t+τk-1,τk-1) (4)
利用上式不難得到
ckx(t-i;q,…,q,i)=γk∈h(t-i,0)hk-2(t-i+q,q)h(t,i)
ckx(t-i;q,…,q,q)=γk∈h(t-i,0)hk-1(t-i+q,q)
ckx(t-i;q,…,q,0)=γk∈h2(t-i,0)hk-2(t-i+q,q)
由上述得到
γk∈hk(t;i)=ckkx(t-i;q,…,q,i)/[ck-2kx(t-i;q,…,q,q).ckx(t-i;q,…,q,0)] (5)
由上式(t=i=0)和假定2(h(0,0)=1)可知,
γk∈=ck-1kx(0;q,…,q,0)/ck-2kx(0;q,…,q,q) (6)
當(dāng)k為奇數(shù)時(shí),
h(t,i)=ckx(t-i;q,…,q,i)/{γ1/kkx[ck-2kx(t-i;q,…,q,q)
.ckx(t-i;q,…,q,0)]1/k} (7)
因此,當(dāng)k為奇數(shù)時(shí),由式(6)和式(7)可以獲得沖激響應(yīng)序列h(t,i).當(dāng)k為偶數(shù)時(shí),由式(6)和式(7)只能得到h(t,i)的幅值,而無法判定其符號(hào).相應(yīng)于時(shí)不變情形的結(jié)果見文獻(xiàn)[6].
特別地,當(dāng)k=3時(shí),式(6)和式(7)成為
γ3∈=c33x(0;q,0)/c3x(0;q,q)
h(t,i)=c3x(t-i;q,i)/{γ1/33∈[c3x(t-i;q,q)
.c3x(t-i;q,0)]1/3}
文獻(xiàn)[5]也得到了上式,式(6)和式(7)將文獻(xiàn)[5]的結(jié)果推廣到任意奇數(shù)階累量情形.
雖然在奇數(shù)階情形可以利用式(6)和式(7)來估計(jì)h(t,i),但是與時(shí)不變情形類似,這種方法將產(chǎn)生較大的估計(jì)誤差,考慮下面的法方程方法.
在式(3)中,令τ1=i,τ2=…=τk-1=0,得到
(8)
式(8)可以變換為
(9)
記
B1(ckx(t+q;-q,0,…,0),ckx(t+q-1;
-q+1,0,…,0),…,ckx(t-q;q,0,…,0))T.
H1(γk∈h(t,0),γk∈h(t,1),…,γk∈h(t,q-1),γk∈h(t,q))T.
在式(9)中令i=-q,-q+1,…,-1,0,1,…,q,可得
A1H1=B1 (10)
將式(6)和式(7)代入A1,求解超定線性方程組(10),可以得到γk∈和h(t,n).
方程(10)是本文得到的基于時(shí)變累量的APTV-MA系統(tǒng)參數(shù)的線性法方程.在A1中,考慮其前q+1行構(gòu)成的三角陣,由假定2可知這個(gè)三角陣的對(duì)角線元素用式(6)和式(7)中的時(shí)變累量代入均不為零,三角陣是滿秩的,故秩A1=q+1.因此,易得如下的唯一識(shí)別定理.
定理1 在模型(1)之下,假定信號(hào)的k階累量已知,則由方程(10)確定的H1的解是唯一的.
在實(shí)際的估計(jì)計(jì)算中,為求解線性方程組(10),需要將A1和B1中的元素用其估計(jì)值代替,有關(guān)循環(huán)統(tǒng)計(jì)量的估計(jì)將在第四節(jié)中討論.
由式(8)得到另一類法方程,記
A2
在式(8)中令i=-q,-q+1,…,-1,0,1,…,q,可得
A2H2=B2 (11)
將式(6)和式(7)代入A2,求解超定線性方程組(11),得到γk∈和h(t,n)的幅值,其符號(hào)可由式(6)和式(7)確定.
方程(11)是得到的基于時(shí)變累量的另一類線性法方程,類似于式(10)的討論可知,秩A2=q+1.因此,可得如下的唯一識(shí)別定理.
定理2 在模型(1)之下,假定信號(hào)的k階累量已知,則由方程(11)確定的H2的解是唯一的.
三、APTV-MA模型定價(jià)
在上節(jié)的討論中,假定APTV-MA系統(tǒng)的階q是已知的,在許多實(shí)際問題中,模型的階是未知的.因此,本節(jié)計(jì)論APTV-MA系統(tǒng)階次的確定問題.
根據(jù)式(8),有
(12)
以上兩式指出,APTV-MA的階數(shù)q是使ckx(t;i,0,…,0)≠0滿足的最大正整數(shù)imax.因此,定階問題轉(zhuǎn)化為判定時(shí)變累量是否為零的問題.與時(shí)不變MA情形一樣,在估計(jì)計(jì)算中,對(duì)于一組在噪聲中觀測(cè)到的有限長度數(shù)據(jù),由樣本估計(jì)檢驗(yàn)累量為零成立與否的閾值是難以確定的.為了解決這個(gè)問題,文獻(xiàn)[16]在討論時(shí)不變MA系統(tǒng)定價(jià)問題時(shí),提出了一種新穎的時(shí)域定階方法,該方法易于實(shí)現(xiàn),且具有較好的數(shù)值魯棒性.受該方法的啟發(fā),討論APTV-MA信號(hào)的定階.
1.基于時(shí)變累量的模型定階
利用信號(hào)x(t)的k階時(shí)變累量構(gòu)造矩陣
(13)
其中,qe(>q)是模型階次的上界(它在實(shí)際應(yīng)用中是容易取到的),由式(12)可知,秩O(t)qe=q+1.因此,模型的階確定轉(zhuǎn)化為矩陣(13)的秩的確定.
在實(shí)際應(yīng)用中,式(13)中的元素用樣本估計(jì)代替,用奇異值分解[2]求O(t)qe的有效秩,即可確定APTV-MA的階數(shù).
由式(13),矩陣O(t)qe的有效秩的確定實(shí)際上等價(jià)于驗(yàn)證
ci+1kx(t;i,0,…,0)≠0,i>0
使得上式成立的最大正整數(shù)imax即是模型的階數(shù)q.因此,基于時(shí)變累量的APTV-MA的定階方法如下,將ckx(t;i,0,…,0)用其樣本估計(jì)kx(t;i,0,…,0)代替,選擇使得
(14)
成立的最大正整數(shù)imax作為階數(shù)q的估計(jì).
2.基于循環(huán)累量的模型定價(jià)
上一小節(jié)討論的是基于t時(shí)刻時(shí)變累量的定階方法,下面討論基于循環(huán)累量的定階方法,由時(shí)變累量和循環(huán)累量的關(guān)系[3],有
t,ckx(t;τ)=0α,Ckx(α;τ)=0
由式(12),有
至少存在一個(gè)α0,使得Ckx(α0;q,0,…,q)≠0;
α,Ckx(α;i,0,…,0)=0,i>q (15)
模型定階轉(zhuǎn)化為找使上兩式成立的最大正整數(shù)imax.根據(jù)上節(jié)的討論,對(duì)于某一α0的qe>q,構(gòu)造矩陣O(c)qe(見式(16)).易知秩O(c)qe=q+1,即模型定階轉(zhuǎn)化為判定矩陣O(c)qe的秩.
在實(shí)際應(yīng)用中,盡管可以用奇異值分解來判定O(c)qe的有效秩,但必須先找α0,這在應(yīng)用中是相當(dāng)麻煩的.然而,類似上一小節(jié)的討論,對(duì)于某一α0,只須找到使Ci+1kx(α0;i,0,…,0)≠0成立的最大的正整數(shù)imax作為模型階數(shù)q的估計(jì).這等價(jià)于選擇使得[supα|Ckx(α;i,0,…,0)|]i+1≠0成立的最大的正整數(shù)imax作為模型階數(shù)q的估計(jì).
O(c)qe
(16)
因此,基于循環(huán)累量的APTV-MA的定階方法如下,將Ckx(α;i,0,…,0)用其樣本估計(jì)代替,選擇使得
(17)
成立的最大正整數(shù)imax作為階數(shù)q的估計(jì).
四、樣本估計(jì)和收斂性討論
在上兩節(jié)中,提出了APTV-MA系統(tǒng)的定階方法和參數(shù)估計(jì)的法方程方法,它們是基于時(shí)變的和循環(huán)的累量進(jìn)行的.在實(shí)際的估計(jì)計(jì)算中,這些統(tǒng)計(jì)量只能由有限長度的樣本(觀測(cè)數(shù)據(jù))進(jìn)行估計(jì)而得到.
對(duì)于滿足模型(2)和假定1-3的APTV-MA信號(hào)y(t),由循環(huán)累量的定義可知,當(dāng)k3時(shí),Cky(α;τ)=Ckx(α;τ).若要估計(jì)ckx(t;τ)和Ckx(α;τ),只須估計(jì)cky(t;τ)和Cky(α;τ).
對(duì)于一般的k階情形,cky(t;τ)和Cky(α;τ)的強(qiáng)相容估計(jì)可以按照文獻(xiàn)[9]的方法得到.
對(duì)于k=1,有c1y(t)=C1y(α)=0.
對(duì)于k=3,y(t)的三階循環(huán)矩和時(shí)變矩估計(jì)分別為
(18)
(19)
其中,,y(t)的三階時(shí)變累量和循環(huán)累量估計(jì)分別為
(20)
(T)3y(α;τ1,τ2)=(T)3y(α;τ1,τ2) (21)
在上述估計(jì)中,循環(huán)矩估計(jì)(T)3y(α;τ)是一個(gè)關(guān)鍵的估計(jì)量,時(shí)變的和循環(huán)的累量估計(jì)都是由它獲得的.另外有一個(gè)量需要估計(jì)的是Am3y.關(guān)于它的估計(jì),有兩種方法.一種是統(tǒng)計(jì)檢驗(yàn)方法[4],另一種是基于DFT的直接檢驗(yàn)方法[10].
在模型(2)的假定之下,由文獻(xiàn)[9]的定理,可知上述關(guān)于時(shí)變累量和循環(huán)累量的估計(jì)都是其真值的強(qiáng)相容估計(jì).由前兩節(jié)的參數(shù)唯一識(shí)別定理和矩陣擾動(dòng)分析理論可知,上兩節(jié)提出的參數(shù)和階數(shù)估計(jì)均是強(qiáng)收斂的.
五、仿真實(shí)驗(yàn)
為了檢驗(yàn)本文提出的基于時(shí)變和循環(huán)累量的模型定階方法和參數(shù)估計(jì)的法方程方法對(duì)APTV-MA信號(hào)的定階和參數(shù)估計(jì)的性能,進(jìn)行了如下的仿真實(shí)驗(yàn).
例 驅(qū)動(dòng)噪聲?(t)是零均值獨(dú)立的指數(shù)分布隨機(jī)數(shù),并且,σ2?=1,γ3∈=2.加性噪聲v(t)取為零均值A(chǔ)R(2)高斯過程:
v(t)-1.6v(t-1)+0.68v(t-2)=e(t) (22)
APTV沖激響應(yīng)序列h(t,n)取為關(guān)于t的周期為2的實(shí)序列:
t為偶數(shù)時(shí):h(t,0)=1,h(t,1)=0.2,t(t,2)=-0.75
t為奇數(shù)時(shí):h(t,0)=0.6,h(t,1)=-0.65,t(t,2)=1.2
即無噪聲污染的信號(hào)x(t)是一個(gè)周期為2的時(shí)變MA(2)過程.當(dāng)t為偶數(shù)時(shí),x(t)=∈(t)+0.2∈(t-1)-0.75∈(t-2).它是最小相位的,其零點(diǎn)為0.7713和-0.9713.當(dāng)t為奇數(shù)時(shí),x(t)=0.6∈(t)-0.65∈(t-1)+1.2∈(t-2).它是非最小相位的,其零點(diǎn)為0.5417±1.3064j.信噪比定義為SNR=10log10(σ2?/σ2v).
在實(shí)驗(yàn)中,信噪比取為0dB,觀測(cè)數(shù)據(jù)由上述參數(shù)和式(2)產(chǎn)生.數(shù)據(jù)長度取為8192,Monte Carlo運(yùn)行次數(shù)為100.基于三階時(shí)變累量和法方程方法(式(10))的參數(shù)估計(jì)的統(tǒng)計(jì)結(jié)果見表1.基于t=0時(shí)刻的三階時(shí)變累量的定階統(tǒng)計(jì)結(jié)果見表2.基于三階循環(huán)累量的定階統(tǒng)計(jì)結(jié)果見表3.
表1 基于三階時(shí)變累量的參數(shù)估計(jì)的統(tǒng)計(jì)結(jié)果
待估參數(shù) | γ3? | h(0,1) | h(0,2) | h(1,0) | h(1,1) | h(1,2) |
真值 | 2.00000 | 0.20000 | -0.75000 | 0.60000 | -0.65000 | 1.20000 |
估計(jì)均值 | 2.02758 | 0.26036 | -0.75511 | 0.61655 | -0.59606 | 1.16906 |
標(biāo)準(zhǔn)偏差 | 0.29703 | 0.09312 | 0.11403 | 0.12365 | 0.12205 | 0.18681 |
表2 基于三階時(shí)變累量的定階統(tǒng)計(jì)結(jié)果 |
i值 | 3y(0;i,0) | i+13y(0;i,0) | ||
均值 | 標(biāo)準(zhǔn)偏差 | 均值 | 標(biāo)準(zhǔn)偏差 | |
i=1 | -0.45716 | 0.17114 | 0.23829 | 0.16242 |
i=2* | 0.88060 | 0.28343 | 0.90109 | 0.83031 |
i=3 | -0.01843 | 0.14205 | 0.00112 | 0.00243 |
i=4 | 0.02757 | 0.20303 | 0.00085 | 0.00387 |
i=5 | 0.04614 | 0.15643 | 0.00050 | 0.00270 |
i=6 | 0.00925 | 0.18165 | 0.00006 | 0.00055 |
i=7 | 0.01157 | 0.12909 | 0.00001 | 0.00003 |
i=8 | 0.01862 | 0.21002 | 0.00004 | 0.00084 |
i=9 | 0.01646 | 0.14416 | 0.00000 | 0.00002 |
i=10 | 0.01380 | 0.17514 | 0.00000 | 0.00002 |
表3 基于三階循環(huán)累量的定階統(tǒng)計(jì)結(jié)果 |
i值 | supα|3y(α;i,0)| | supα|i+13y(α;i,0)| | ||
均值 | 標(biāo)準(zhǔn)偏差 | 均值 | 標(biāo)準(zhǔn)偏差 | |
i=1 | 0.74568 | 0.11763 | 0.56988 | 0.18373 |
i=2* | 1.17688 | 0.16014 | 1.72152 | 0.70514 |
i=3 | 0.47164 | 0.05893 | 0.05441 | 0.03061 |
i=4 | 0.44282 | 0.05312 | 0.01968 | 0.01309 |
i=5 | 0.40810 | 0.04865 | 0.00577 | 0.00514 |
i=6 | 0.39960 | 0.04397 | 0.00211 | 0.00184 |
i=7 | 0.37370 | 0.04522 | 0.00058 | 0.00074 |
i=8 | 0.38182 | 0.04109 | 0.00026 | 0.00031 |
i=9 | 0.36517 | 0.03844 | 0.00007 | 0.00009 |
i=10 | 0.38063 | 0.04190 | 0.00005 | 0.00009 |
表1說明參數(shù)估計(jì)的法方程方法具有較好的估計(jì)勻值和較小的標(biāo)準(zhǔn)偏差.特別要指出的是,例題中的系統(tǒng)不僅是周期時(shí)變的,而且是非最小相位的. 表2~3說明基于時(shí)變和循環(huán)累量的定階方法(與直接判定統(tǒng)計(jì)量相比較)具有很好的數(shù)值穩(wěn)定性,且每次運(yùn)算都能給出準(zhǔn)確的階數(shù)估計(jì)(q=imax=2). 在實(shí)驗(yàn)中,盡管信噪比較低,但是所用的數(shù)據(jù)長度較大(與傳統(tǒng)的二階統(tǒng)計(jì)量方法和累量方法比較),這是因?yàn)檠h(huán)統(tǒng)計(jì)量具有較大的估計(jì)方差. 六、結(jié)束語 |
評(píng)論
查看更多