一種實(shí)現(xiàn)Fornasini-Marchesini 模型的新方法的研究①

曹中泳 程 驊 吳徐冬子 劉 昶

(武漢科技大學(xué)信息科學(xué)與工程學(xué)院 武漢 430080)

0 引 言

在過(guò)去的幾十年中，多維系統(tǒng)吸引了極大的關(guān)注，這主要?dú)w功于多維系統(tǒng)在各種工程領(lǐng)域的廣泛使用和潛在使用，如信號(hào)和圖像處理、熱工藝、醫(yī)療應(yīng)用、無(wú)線傳感器網(wǎng)絡(luò)等[1-3]。

多維系統(tǒng)理論中一個(gè)基本問(wèn)題是通過(guò)確定的多維狀態(tài)空間模型來(lái)實(shí)現(xiàn)給定的有理傳遞函數(shù)或傳遞矩陣，該模型通常是Roesser模型或Fornasini-Marchesini(F-M)模型。與傳統(tǒng)的1維(1-D)情況不同，通常n維(n≥2)濾波器或者系統(tǒng)的最小狀態(tài)空間實(shí)現(xiàn)很難獲得。因此研究可以得到低階n維狀態(tài)空間的實(shí)現(xiàn)新方法就顯得格外重要。此外，許多研究人員對(duì)Roesser狀態(tài)空間模型進(jìn)行了廣泛和深入的研究，只有少數(shù)文獻(xiàn)闡述了關(guān)于F-M模型實(shí)現(xiàn)問(wèn)題。比如，Alpay和Dubi[4]給出了直接構(gòu)建一個(gè)n維(n≥2)F-M實(shí)現(xiàn)方法。但是，這種方法求得F-M模型實(shí)現(xiàn)矩陣階數(shù)卻非常高[5]。

本文針對(duì)Alpay和Dubi[4]的n維F-M模型的實(shí)現(xiàn)矩陣求解方法，提出一種新解法，該方法獲取實(shí)現(xiàn)矩陣的階更低并且更易于實(shí)現(xiàn)。將此算法運(yùn)用到多輸入多輸出(multiple input multiple output，MIMO)雷達(dá)系統(tǒng)中，提高了雷達(dá)目標(biāo)檢測(cè)[6-8]的效率以及反應(yīng)靈敏度。

1 F-M狀態(tài)空間模型

對(duì)于n維MIMO線性離散系統(tǒng)，F(xiàn)-M狀態(tài)空間模型的描述如下[9-11]：

x(i1+1,i2+1,…,in+1)

=A1x(i1,i2+1,…,in+1)+…

+Anx(i1+1,…,in-1+1,in)

+B1u(i1,i2+1,…,in+1)+…

+Bnu(i1+1,…,in-1+1,in)

(1)

y(i1,…,in)=Cx(i1,…,in)+Du(i1,…,in)

(2)

其中，x(i1,…,in)∈Rr、u(i1,…,in)∈Rl、y(i1,…,in)∈Rm分別是局部狀態(tài)向量、輸入向量、輸出向量。A、B、C、D是實(shí)數(shù)矩陣，且A1,…,An∈Rr×r,B1,…,Bn∈Rr×l,C∈Rm×r,D∈Rm×l。

n維信號(hào)u(i1,…,in)的n維z變換的定義如下。

其中，z1,…,zn是單位延遲算子。

對(duì)式(1)進(jìn)行n維z變換后可得：

將z1z2…zn整理可得：

X(z1,…,zn)=(A1z1+…+Anzn)X(z1,…,zn)

+(B1z1+…+Bnzn)U(z1,…,zn)

若將上述等式按X(z1,…,zn)整理可得：

X(z1,…,zn)=(I-A1z1-…-Anzn)

=(B1z1+…+Bnzn)U(z1,…,zn)

(3)

式(2)的n維z變換如下。

Y(z1,…,zn)=CX(z1,…,zn)+DU(z1,…,zn)

代入式(3)后可得如下所示傳遞矩陣：

(4)

(5)

其中，η=degm(z)=max{|γ|?γs.t.mγ≠0}。 若系統(tǒng)的傳遞函數(shù)h(z)=q(z)/m(z)，其中q(z)和m(z)是n維多項(xiàng)式，當(dāng)m(0,…,0)≠0，則稱(chēng)這個(gè)系統(tǒng)為因果系統(tǒng)[12]。如果一個(gè)有理矩陣的每一項(xiàng)都是因果的，則該有理矩陣是因果有理矩陣。

對(duì)于一個(gè)給定的傳遞矩陣H(z)，如果式(4)成立，則A、B、C、D被稱(chēng)為傳遞矩陣H(z)的F-M 狀態(tài)空間模型的實(shí)現(xiàn),該實(shí)現(xiàn)的充分必要條件就是H(z)是因果矩陣。由Galkowski教授[13]的研究可知，對(duì)于某個(gè)因果有理傳遞函數(shù)或傳遞矩陣，存在很多F-M模型實(shí)現(xiàn)，但不一定具有相同的階次。由于狀態(tài)空間實(shí)現(xiàn)的順序與系統(tǒng)的計(jì)算或硬件實(shí)現(xiàn)的成本密切相關(guān)，因此期望具有低階狀態(tài)空間實(shí)現(xiàn)。對(duì)于1維情況，已經(jīng)得到很好的發(fā)展。然而，對(duì)于n維(n≥2)情況，情況變得更加復(fù)雜和困難，因此，在n維情況下，尋求低階實(shí)現(xiàn)尤為重要。

本文將研究如下形式的傳遞矩陣的F-M 狀態(tài)空間模型的實(shí)現(xiàn)問(wèn)題。

(6)

2 Alpay和Dubi的實(shí)現(xiàn)方法

針對(duì)n維因果有理傳遞函數(shù)，Alpay和Dubi[4]提出了一個(gè)構(gòu)建F-M模型實(shí)現(xiàn)的建設(shè)性方法。但是Alpay和Dubi的實(shí)現(xiàn)方法通常會(huì)產(chǎn)生不必要的高階次，這就促使本文研究一個(gè)替代的實(shí)現(xiàn)方法，考慮到給定傳遞函數(shù)或傳遞矩陣的結(jié)構(gòu)特征,階數(shù)更低的實(shí)現(xiàn)方法是可行的[1]。

如果F(z)是一個(gè)矩陣值函數(shù)，本文用Μ(F(z))來(lái)表示Cm值函數(shù)的向量空間，其中值函數(shù)是F(z)所有列的線性組合。

定義1設(shè)M是Ω?Cn中的函數(shù)解析向量空間。如果相關(guān)的Gleason問(wèn)題在Μ中是可解的，則稱(chēng)其為解析不變量。比如，對(duì)任意的f(z)∈Μ，ε∈Ω，存在一個(gè)函數(shù)g1,ε(z),…,gn,ε(z)∈Μ滿足：

(7)

定義2設(shè)Μ是在Cn上定義的函數(shù)向量空間。若對(duì)任意f(z)∈Μ存在g1(z),…,gn(z)∈ζ滿足：

(8)

則稱(chēng)ζ為M的后向移位空間。

定理1設(shè)H(z)：Ω?Cn→Cm×l(Ω包含(0,…,0))，H(z)具有滿足式(4)的F-M實(shí)現(xiàn)矩陣(A，B，C，D)，當(dāng)且僅當(dāng)M(H(z))具有有限維度和解析不變后移空間ζ。

設(shè)h(z)=q(z)/m(z)是單輸入單輸出(SISO)的多維系統(tǒng)有理傳遞函數(shù)，q(z)和m(z)分別是分子和分母多項(xiàng)式。假設(shè)h(z)是因果關(guān)系，即p(0,…，0)≠0，將k={degp(z), degm(z)}稱(chēng)為有理函數(shù)h(z)的度。Alpay和Dubi[4]所構(gòu)造空間如下：

(9)

它的維度是有限，并且還是Μ(h(z))的后向移位空間?？梢钥闯鐾ㄟ^(guò)Alpay和Dubi的方法得到實(shí)現(xiàn)矩陣的階數(shù)r與ζ中元素的數(shù)量相同，并且可以由式(10)計(jì)算得到。

(10)

因此，對(duì)于某個(gè)確定的n，Alpay和Dubi的方法所獲得的實(shí)現(xiàn)階數(shù)將完全由h(z)的k值確定?，F(xiàn)在容易驗(yàn)證，即使對(duì)于相對(duì)較小的n和k，實(shí)現(xiàn)階數(shù)通常也相當(dāng)高，隨著n或k的增加，r的值也快速增加。表1顯示在n= 2,3,4，5和k= 2，3，…，6的情況下，通過(guò)Alpay和Dubi的方法獲得的實(shí)現(xiàn)階數(shù)。

表1 Alpay和Dubi的實(shí)現(xiàn)階數(shù)

但是，文獻(xiàn)[14，15]曾經(jīng)提出對(duì)于n維情況下的實(shí)現(xiàn)階數(shù)不僅取決于給定傳遞函數(shù)h(z)的k值，而且還取決于給定傳遞函數(shù)h(z)的結(jié)構(gòu)甚至系數(shù)。特別是在n維(n≥2)情況下，具有相同k值但結(jié)構(gòu)不同的傳遞函數(shù)可能具有不同實(shí)現(xiàn)階數(shù)。這個(gè)事實(shí)意味著當(dāng)h(z)缺少一些單項(xiàng)式時(shí)得到的實(shí)現(xiàn)階數(shù)遠(yuǎn)低于Alpay和Dubi不考慮結(jié)構(gòu)特征的方法獲得的實(shí)現(xiàn)階數(shù)[16]。

例1

其中，z=(z1,z2,z3)是獨(dú)立變量，q20、q01、m10、m01為系數(shù)，由于h(z)的度為3，因?yàn)閙(0,0,0)=1≠0，顯然h(z)是因果的。

可以使用所有zα/p(z)構(gòu)造ζ，其中α=(α1,α2,α3)，0≤|α|≤3。

由Alpay和Dubi的方法可得F-M實(shí)現(xiàn)矩陣階數(shù)r=20，與式(10)所得結(jié)果一樣。利用以下矩陣來(lái)代表的實(shí)現(xiàn)矩陣(A,B,C,D)。

D=0不僅滿足式(4)而且它的階次為7，遠(yuǎn)低于Alpay和Dubi的方法計(jì)算出來(lái)的階次。本文所關(guān)注的是如何構(gòu)建這樣一個(gè)低階實(shí)現(xiàn)，這將在下一節(jié)中討論。

3 新方法步驟

步驟1首先取得每一列中系數(shù)非0的單項(xiàng)式zα=zα1…zαn，然后將這些單項(xiàng)式組成一個(gè)如下所示的矩陣[17，18]。

(11)

(12)

其中λT為ml取反后再加上1的式子的系數(shù)。重置Bij，即將第s個(gè)位置的項(xiàng)置1，其他的項(xiàng)保持不變。重置Aij，即第i行、j列的數(shù)值置為λT。

(13)

(14)

通過(guò)式(14)，構(gòu)造的矩陣(A,B,C,D)即為H(z)的F-M實(shí)現(xiàn)。其中A=(A1…An)，B=(B1…Bn)。

4 實(shí) 例

利用上述新方法對(duì)例1：

=q(z)/m(z)

進(jìn)行實(shí)現(xiàn)矩陣的求解。

G1(z)=

h(z)-h(0)=

=G1(z)(z1B11+z2B12+z3B13)

=G1(z)(z1κ11+z2κ12+z3κ13)

進(jìn)一步，令κj1=κj2=κj3=0。因?yàn)棣?2(z)=z1β11(z)，按照步驟中的方法令κ21=1。

同理κ32(2)=κ43(2)=κ53(7)=κ61(5)=κ72(1)=1。剩余其他的項(xiàng)為0。Aj=A1j,Bj=B1j,j=1，2，3，G(z)=G1(z)，C=G(0)=[1000000]，D=h(0)=0。此時(shí)實(shí)現(xiàn)矩陣的階數(shù)為r=7，遠(yuǎn)低于Alpay和Dubi教授所求出的20。

在雷達(dá)目標(biāo)跟蹤系統(tǒng)中，通常將被跟蹤的目標(biāo)建立如下?tīng)顟B(tài)方程[19]：

X(k+1,k)=FX(k)+Qω(k-1,k-1)

(15)

而雷達(dá)的觀測(cè)方程可表示為

W(k)=h(k)+x(k)+ν

(16)

令k=z，上述狀態(tài)方程與觀測(cè)方程可以轉(zhuǎn)化為

X(z+1,z)=FX(z)+Qω(z,z)

W(z)=h(z)x(z)+v

(17)

若系統(tǒng)矩陣v=0，則系統(tǒng)是嚴(yán)格因果系統(tǒng)。

基于上述公式,利用 F-M狀態(tài)空間模型來(lái)表示狀態(tài)方程與觀測(cè)方程[20]，對(duì)MIMO雷達(dá)目標(biāo)檢測(cè)的研究也即是對(duì)F-M狀態(tài)空間模型的研究。

圖1是3個(gè)求得坐標(biāo)后的目標(biāo)位置信息。

圖1 待測(cè)目標(biāo)位置信息

通過(guò)定位算法進(jìn)行坐標(biāo)計(jì)算后，3個(gè)位置的坐標(biāo)分別用R=(xr,yr),S=(xs,ys),W=(xw,yw)表示。則可對(duì)系統(tǒng)建立如下3維F-M狀態(tài)空間模型：

x(xr,yr,t+1)

=A1x(xs,ys,t+1)+A2x(xw,yw,t+1)

+A3x(xr,yr,t)+B1u(xs,ys,t+1)

+B2u(xw,yw,t+1)+B3u(xr,yr,t+1)

y(xr,yr,t)=Cx(xr,yr,t)+Du(xr,yr,t)

(18)

將式(15)進(jìn)行z變換之后，可得系統(tǒng)傳遞函數(shù)為

(19)

成功實(shí)現(xiàn)建模后，對(duì)MIMO雷達(dá)系統(tǒng)的研究就轉(zhuǎn)移到F-M 狀態(tài)空間模型上。

例2對(duì)MIMO雷達(dá)系統(tǒng)進(jìn)行系統(tǒng)建模后得到的3×2傳遞矩陣如下所示[20]。

m0(z)=1-d01z2-d02z3+d03z1z2

m2(z)=1+d11z1-d12z1z2+d13z1z2z3

通過(guò)新方法可以得到以下結(jié)果：

因?yàn)镃1=NHT1,N1=NHT1Ψ1，可得：

代入式(13)可以得到：

同理，另一結(jié)果如下：

C2=NHT2,N2=NHT2Ψ2

代入式(13)得到如下結(jié)果：

由上述結(jié)果知，傳遞矩陣H(z)的實(shí)現(xiàn)矩陣如下。

A1=diag{A11,A21}A2=diag{A12,A22}

A3=diag{A13,A23}

B1=diag{B11,B21}B2=diag{B12,B22}

B3=diag{B13,B23}

C=[C1C2]D=0

若使用Alpay和Dubi所提出方法求得的實(shí)現(xiàn)矩陣階數(shù)是45，而新方法所求階數(shù)為10。

將此算法運(yùn)用到MIMO雷達(dá)系統(tǒng)的多維處理過(guò)程中，采用Matlab工具進(jìn)行仿真實(shí)驗(yàn)。圖2與圖3分別是2種方法在x軸與y軸上速度協(xié)方差的比較。Alpay和Dubi算法的檢測(cè)速度協(xié)方差在收斂前一直稍大于F-M狀態(tài)空間模型的步數(shù)，并且后者的收斂速度快于前者。這表明采用新方法進(jìn)行降階后，可以提高雷達(dá)系統(tǒng)的定位準(zhǔn)確度以及反應(yīng)速度。

圖2 x方向預(yù)測(cè)和更新速度誤差協(xié)方差

圖3 y方向預(yù)測(cè)和更新速度誤差協(xié)方差

5 結(jié) 論

本文針對(duì)n維F-M模型實(shí)現(xiàn)問(wèn)題提出了一種新的求解方法，與Alpay和Dubi的方法相比，它可以得到更低階實(shí)現(xiàn)矩陣。具體而言，即使對(duì)于簡(jiǎn)單的有理傳遞函數(shù)，Alpay和Dubi給出的方法通常會(huì)產(chǎn)生具有不必要的高階實(shí)現(xiàn)矩陣。并且，通過(guò)考慮給定傳遞函數(shù)或傳遞矩陣的結(jié)構(gòu)性質(zhì)，可以構(gòu)造一個(gè)階數(shù)較低的實(shí)現(xiàn)矩陣，并通過(guò)舉例說(shuō)明該方法的可實(shí)施性及有效性。將此算法運(yùn)用到MIMO雷達(dá)系統(tǒng)的多維處理過(guò)程中，可以大幅降低系統(tǒng)的數(shù)據(jù)處理時(shí)間，提高雷達(dá)系統(tǒng)對(duì)所跟蹤物體的反應(yīng)速度和定位的準(zhǔn)確度。今后更進(jìn)一步的研究是盡可能減少約束條件，使之能夠獲取到更加低階的系統(tǒng)實(shí)現(xiàn)矩陣。