幾乎周期時(shí)變MA系統(tǒng)的辨識(shí)

幾乎周期時(shí)變MA系統(tǒng)的辨識(shí)

本文利用循環(huán)統(tǒng)計(jì)量討論幾乎周期時(shí)變滑動(dòng)平均(APTV-MA)系統(tǒng)的辨識(shí)問(wèn)題.提出了基于奇數(shù)階時(shí)變累量的參數(shù)估計(jì)的法方程方法，并導(dǎo)出了基于時(shí)變累量和循環(huán)累量的系統(tǒng)模型定階方法.仿真實(shí)驗(yàn)說(shuō)明了本文所提方法的性能.
　　關(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ú)立的功率譜未知的零均值高斯(有色)過(guò)程(可以是非平穩(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ì)量來(lái)估計(jì)
　　上述問(wèn)題可以看作時(shí)不變MA系統(tǒng)在幾乎周期時(shí)變情形的推廣.這一類問(wèn)題在通訊，水文氣象等領(lǐng)域具有廣泛的應(yīng)用［5，11，14］.在高斯信號(hào)情形，文獻(xiàn)［15］提出了基于二階統(tǒng)計(jì)量的近似最大似然技術(shù)估計(jì)信號(hào)參數(shù).相關(guān)問(wèn)題的討論可見(jiàn)文獻(xiàn)［8］和［13］.最近［5］，Dandawate和Giannakis對(duì)于幾乎周期時(shí)變MA信號(hào)進(jìn)行了較深入的研究，他們利用高階循環(huán)統(tǒng)計(jì)量來(lái)辨識(shí)APTV-MA系統(tǒng)，提出了模型參數(shù)估計(jì)的線性和非線性算法以及一種統(tǒng)計(jì)的定階方法.
　　眾所周知，對(duì)于時(shí)不變MA系統(tǒng)，基于高階累量的參數(shù)估計(jì)方法可分為閉式解，法方程和非線性最優(yōu)化解三大類［12］.由于不涉及極值問(wèn)題求解，前二種線性方法尤其令人重視［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)參見(jiàn)文獻(xiàn)[3].
當(dāng)k3時(shí),由假定3和累量性質(zhì)有ckx(t;τ)=cky(t;τ).即觀測(cè)信號(hào)y(t)與無(wú)加性噪聲污染的輸出信號(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)的幅值，而無(wú)法判定其符號(hào).相應(yīng)于時(shí)不變情形的結(jié)果見(jiàn)文獻(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)來(lái)估計(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í)際問(wèn)題中，模型的階是未知的.因此，本節(jié)計(jì)論APTV-MA系統(tǒng)階次的確定問(wèn)題.
　　根據(jù)式(8)，有

　(12)

以上兩式指出，APTV-MA的階數(shù)q是使ckx(t;i,0,…,0)≠0滿足的最大正整數(shù)imax.因此，定階問(wèn)題轉(zhuǎn)化為判定時(shí)變累量是否為零的問(wèn)題.與時(shí)不變MA情形一樣，在估計(jì)計(jì)算中，對(duì)于一組在噪聲中觀測(cè)到的有限長(zhǎng)度數(shù)據(jù)，由樣本估計(jì)檢驗(yàn)累量為零成立與否的閾值是難以確定的.為了解決這個(gè)問(wèn)題，文獻(xiàn)［16］在討論時(shí)不變MA系統(tǒng)定價(jià)問(wèn)題時(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(見(jiàn)式(16)).易知秩O(c)qe=q+1，即模型定階轉(zhuǎn)化為判定矩陣O(c)qe的秩.
　　在實(shí)際應(yīng)用中，盡管可以用奇異值分解來(lái)判定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ì)量只能由有限長(zhǎng)度的樣本(觀測(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)高斯過(guò)程：

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
即無(wú)噪聲污染的信號(hào)x(t)是一個(gè)周期為2的時(shí)變MA(2)過(guò)程.當(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ù)長(zhǎng)度取為8192，Monte Carlo運(yùn)行次數(shù)為100.基于三階時(shí)變累量和法方程方法(式(10))的參數(shù)估計(jì)的統(tǒng)計(jì)結(jié)果見(jiàn)表1.基于t=0時(shí)刻的三階時(shí)變累量的定階統(tǒng)計(jì)結(jié)果見(jiàn)表2.基于三階循環(huán)累量的定階統(tǒng)計(jì)結(jié)果見(jiàn)表3.

表1　基于三階時(shí)變累量的參數(shù)估計(jì)的統(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

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