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

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

本文利用循環(huán)統(tǒng)計量討論幾乎周期時變滑動平均(APTV-MA)系統(tǒng)的辨識問題.提出了基于奇數(shù)階時變累量的參數(shù)估計的法方程方法，并導(dǎo)出了基于時變累量和循環(huán)累量的系統(tǒng)模型定階方法.仿真實驗說明了本文所提方法的性能.
　　關(guān)鍵詞:循環(huán)統(tǒng)計量；幾乎周期時變MA系統(tǒng)；法方程；系統(tǒng)辨識

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階幾乎周期時變滑動平均(APTV-MA)信號模型

　(1)

并且，假設(shè)下列假定成立.
　　假定1　激勵噪聲?(t)是零均值獨立同分布的，且存在一個奇數(shù)k(3)，使得其k階累量γk?滿足0＜|γk∈|＜+∞.
　　假定2　對于每一個固定的n,有限沖激響應(yīng)序列h(t,n)是關(guān)于t的幾乎周期實序列(可以是非最小相位的)，并且，h(0,0)=1;h(t,0)≠0,h(t,q)≠0,t.
　　在實際應(yīng)用中，信號x(t)是在噪聲中被觀測的：

y(t)=x(t)+v(t),t=0,1,…,T-1　(2)

其中，加性噪聲v(t)滿足
　　假定3　v(t)是一個與∈(t)獨立的功率譜未知的零均值高斯(有色)過程(可以是非平穩(wěn)的).
　　在模型(2)下，y(t)是非平穩(wěn)的，原有的討論時不變MA系統(tǒng)的辨識方法不能直接應(yīng)用.易證y(t)是循環(huán)平穩(wěn)的.因此，我們的目的是根據(jù)觀測值{y(t),t=0,1,…,T-1}利用循環(huán)統(tǒng)計量來估計
　　上述問題可以看作時不變MA系統(tǒng)在幾乎周期時變情形的推廣.這一類問題在通訊，水文氣象等領(lǐng)域具有廣泛的應(yīng)用［5，11，14］.在高斯信號情形，文獻［15］提出了基于二階統(tǒng)計量的近似最大似然技術(shù)估計信號參數(shù).相關(guān)問題的討論可見文獻［8］和［13］.最近［5］，Dandawate和Giannakis對于幾乎周期時變MA信號進行了較深入的研究，他們利用高階循環(huán)統(tǒng)計量來辨識APTV-MA系統(tǒng)，提出了模型參數(shù)估計的線性和非線性算法以及一種統(tǒng)計的定階方法.
　　眾所周知，對于時不變MA系統(tǒng)，基于高階累量的參數(shù)估計方法可分為閉式解，法方程和非線性最優(yōu)化解三大類［12］.由于不涉及極值問題求解，前二種線性方法尤其令人重視［7，17］.但是到目前為止，對于幾乎周期時變MA系統(tǒng)的辨識，這兩種方法的研究還較小.本文討論利用循環(huán)統(tǒng)計量辨識APTV-MA系統(tǒng)的新方法，提出了參數(shù)估計的法方程方法和基于時變累量和循環(huán)累量的時域定階方法.

二、參數(shù)估計方法
　　在本節(jié)先假定幾乎周期時變MA系統(tǒng)的階q已知，推導(dǎo)APTV－MA信號參數(shù)的法方程.
　　實信號z(t)的k階時變的和循環(huán)的累量分別定義［3，5］為

其中，τ=(τ1,…,τk-1)．有關(guān)循環(huán)統(tǒng)計量的詳細定義和記號參見文獻[3].
當k3時,由假定3和累量性質(zhì)有ckx(t;τ)=cky(t;τ).即觀測信號y(t)與無加性噪聲污染的輸出信號x(t)具有相同的時變累量，因此，在基于高階時變和循環(huán)累量的分析中不必考慮加性高斯噪聲的影響.
　　由累量性質(zhì)可知，滿足模型(1)和假定1～2的APTV-MA信號x(t)的k(＞3)階時變累量為

　(3)

它是時不變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)

當k為奇數(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)
　　因此，當k為奇數(shù)時，由式(6)和式(7)可以獲得沖激響應(yīng)序列h(t,i).當k為偶數(shù)時，由式(6)和式(7)只能得到h(t,i)的幅值，而無法判定其符號.相應(yīng)于時不變情形的結(jié)果見文獻［6］.
　　特別地，當k=3時，式(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}

文獻［5］也得到了上式，式(6)和式(7)將文獻［5］的結(jié)果推廣到任意奇數(shù)階累量情形.
　　雖然在奇數(shù)階情形可以利用式(6)和式(7)來估計h(t,i),但是與時不變情形類似，這種方法將產(chǎn)生較大的估計誤差，考慮下面的法方程方法.
　　在式(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)是本文得到的基于時變累量的APTV-MA系統(tǒng)參數(shù)的線性法方程.在A1中，考慮其前q+1行構(gòu)成的三角陣，由假定2可知這個三角陣的對角線元素用式(6)和式(7)中的時變累量代入均不為零，三角陣是滿秩的，故秩A1=q+1.因此，易得如下的唯一識別定理.
　　定理1　在模型(1)之下，假定信號的k階累量已知，則由方程(10)確定的H1的解是唯一的.
　　在實際的估計計算中，為求解線性方程組(10)，需要將A1和B1中的元素用其估計值代替，有關(guān)循環(huán)統(tǒng)計量的估計將在第四節(jié)中討論.
　　由式(8)得到另一類法方程，記
A2

在式(8)中令i=-q,-q+1,…,-1,0,1,…,q,可得

A2H2=B2　(11)

將式(6)和式(7)代入A2，求解超定線性方程組(11)，得到γk∈和h(t,n)的幅值，其符號可由式(6)和式(7)確定.
　　方程(11)是得到的基于時變累量的另一類線性法方程，類似于式(10)的討論可知，秩A2=q+1.因此，可得如下的唯一識別定理.
　　定理2　在模型(1)之下，假定信號的k階累量已知，則由方程(11)確定的H2的解是唯一的.

三、APTV-MA模型定價
　　在上節(jié)的討論中，假定APTV-MA系統(tǒng)的階q是已知的，在許多實際問題中，模型的階是未知的.因此，本節(jié)計論APTV-MA系統(tǒng)階次的確定問題.
　　根據(jù)式(8)，有

　(12)

以上兩式指出，APTV-MA的階數(shù)q是使ckx(t;i,0,…,0)≠0滿足的最大正整數(shù)imax.因此，定階問題轉(zhuǎn)化為判定時變累量是否為零的問題.與時不變MA情形一樣，在估計計算中，對于一組在噪聲中觀測到的有限長度數(shù)據(jù)，由樣本估計檢驗累量為零成立與否的閾值是難以確定的.為了解決這個問題，文獻［16］在討論時不變MA系統(tǒng)定價問題時，提出了一種新穎的時域定階方法，該方法易于實現(xiàn)，且具有較好的數(shù)值魯棒性.受該方法的啟發(fā)，討論APTV-MA信號的定階.
　　1.基于時變累量的模型定階
　　利用信號x(t)的k階時變累量構(gòu)造矩陣
　(13)

其中，qe(＞q)是模型階次的上界(它在實際應(yīng)用中是容易取到的)，由式(12)可知，秩O(t)qe=q+1.因此，模型的階確定轉(zhuǎn)化為矩陣(13)的秩的確定.
　　在實際應(yīng)用中，式(13)中的元素用樣本估計代替，用奇異值分解［2］求O(t)qe的有效秩，即可確定APTV-MA的階數(shù).
　　由式(13)，矩陣O(t)qe的有效秩的確定實際上等價于驗證

ci+1kx(t;i,0,…，0)≠0,i＞0

使得上式成立的最大正整數(shù)imax即是模型的階數(shù)q．因此，基于時變累量的APTV-MA的定階方法如下,將ckx(t;i,0,…,0)用其樣本估計kx(t;i,0,…,0)代替，選擇使得

　(14)

成立的最大正整數(shù)imax作為階數(shù)q的估計.
　　2.基于循環(huán)累量的模型定價
　　上一小節(jié)討論的是基于t時刻時變累量的定階方法，下面討論基于循環(huán)累量的定階方法，由時變累量和循環(huán)累量的關(guān)系［3］，有

t,ckx(t;τ)=0α,Ckx(α;τ)=0

由式(12)，有
　　至少存在一個α0，使得Ckx(α0;q,0,…,q)≠0;

α,Ckx(α;i,0,…,0)=0,i＞q　(15)

模型定階轉(zhuǎn)化為找使上兩式成立的最大正整數(shù)imax.根據(jù)上節(jié)的討論，對于某一α0的qe＞q，構(gòu)造矩陣O(c)qe(見式(16)).易知秩O(c)qe=q+1，即模型定階轉(zhuǎn)化為判定矩陣O(c)qe的秩.
　　在實際應(yīng)用中，盡管可以用奇異值分解來判定O(c)qe的有效秩，但必須先找α0，這在應(yīng)用中是相當麻煩的.然而，類似上一小節(jié)的討論，對于某一α0，只須找到使Ci+1kx(α0;i,0,…,0)≠0成立的最大的正整數(shù)imax作為模型階數(shù)q的估計.這等價于選擇使得［supα|Ckx(α;i,0,…,0)|］i+1≠0成立的最大的正整數(shù)imax作為模型階數(shù)q的估計.
O(c)qe
　(16)
　　因此，基于循環(huán)累量的APTV-MA的定階方法如下，將Ckx(α;i,0,…,0)用其樣本估計代替，選擇使得

　(17)

成立的最大正整數(shù)imax作為階數(shù)q的估計.

四、樣本估計和收斂性討論
　　在上兩節(jié)中，提出了APTV-MA系統(tǒng)的定階方法和參數(shù)估計的法方程方法，它們是基于時變的和循環(huán)的累量進行的.在實際的估計計算中，這些統(tǒng)計量只能由有限長度的樣本(觀測數(shù)據(jù))進行估計而得到.
　　對于滿足模型(2)和假定1-3的APTV-MA信號y(t)，由循環(huán)累量的定義可知，當k3時，Cky(α;τ)=Ckx(α;τ).若要估計ckx(t;τ)和Ckx(α;τ)，只須估計cky(t;τ)和Cky(α;τ).
　　對于一般的k階情形，cky(t;τ)和Cky(α;τ)的強相容估計可以按照文獻［9］的方法得到.
　　對于k=1,有c1y(t)=C1y(α)=0.
　　對于k=3,y(t)的三階循環(huán)矩和時變矩估計分別為

　(18)
　(19)

其中，，y(t)的三階時變累量和循環(huán)累量估計分別為

　(20)
(T)3y(α;τ1,τ2)=(T)3y(α;τ1,τ2)　(21)

　　在上述估計中，循環(huán)矩估計(T)3y(α;τ)是一個關(guān)鍵的估計量，時變的和循環(huán)的累量估計都是由它獲得的.另外有一個量需要估計的是Am3y.關(guān)于它的估計，有兩種方法.一種是統(tǒng)計檢驗方法［4］，另一種是基于DFT的直接檢驗方法［10］.
　　在模型(2)的假定之下，由文獻［9］的定理，可知上述關(guān)于時變累量和循環(huán)累量的估計都是其真值的強相容估計.由前兩節(jié)的參數(shù)唯一識別定理和矩陣擾動分析理論可知，上兩節(jié)提出的參數(shù)和階數(shù)估計均是強收斂的.

五、仿真實驗
　　為了檢驗本文提出的基于時變和循環(huán)累量的模型定階方法和參數(shù)估計的法方程方法對APTV-MA信號的定階和參數(shù)估計的性能，進行了如下的仿真實驗.
　　例　驅(qū)動噪聲?(t)是零均值獨立的指數(shù)分布隨機數(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的實序列：
t為偶數(shù)時：h(t,0)=1,h(t,1)=0.2,t(t,2)=-0.75
t為奇數(shù)時：h(t,0)=0.6,h(t,1)=-0.65,t(t,2)=1.2
即無噪聲污染的信號x(t)是一個周期為2的時變MA(2)過程.當t為偶數(shù)時，x(t)=∈(t)+0.2∈(t-1)-0.75∈(t-2).它是最小相位的，其零點為0.7713和-0.9713.當t為奇數(shù)時，x(t)=0.6∈(t)-0.65∈(t-1)+1.2∈(t-2).它是非最小相位的，其零點為0.5417±1.3064j.信噪比定義為SNR=10log10(σ2?/σ2v).
　　在實驗中，信噪比取為0dB，觀測數(shù)據(jù)由上述參數(shù)和式(2)產(chǎn)生.數(shù)據(jù)長度取為8192，Monte Carlo運行次數(shù)為100.基于三階時變累量和法方程方法(式(10))的參數(shù)估計的統(tǒng)計結(jié)果見表1.基于t=0時刻的三階時變累量的定階統(tǒng)計結(jié)果見表2.基于三階循環(huán)累量的定階統(tǒng)計結(jié)果見表3.

表1　基于三階時變累量的參數(shù)估計的統(tǒng)計結(jié)果

i值	3y(0;i,0)		i+13y(0;i,0)
	均值	標準偏差	均值	標準偏差
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)|
	均值	標準偏差	均值	標準偏差
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