曹中泳 程 驊 吳徐冬子 劉 昶
(武漢科技大學(xué)信息科學(xué)與工程學(xué)院 武漢 430080)
在過去的幾十年中,多維系統(tǒng)吸引了極大的關(guān)注,這主要歸功于多維系統(tǒng)在各種工程領(lǐng)域的廣泛使用和潛在使用,如信號和圖像處理、熱工藝、醫(yī)療應(yīng)用、無線傳感器網(wǎng)絡(luò)等[1-3]。
多維系統(tǒng)理論中一個基本問題是通過確定的多維狀態(tài)空間模型來實現(xiàn)給定的有理傳遞函數(shù)或傳遞矩陣,該模型通常是Roesser模型或Fornasini-Marchesini(F-M)模型。與傳統(tǒng)的1維(1-D)情況不同,通常n維(n≥2)濾波器或者系統(tǒng)的最小狀態(tài)空間實現(xiàn)很難獲得。因此研究可以得到低階n維狀態(tài)空間的實現(xiàn)新方法就顯得格外重要。此外,許多研究人員對Roesser狀態(tài)空間模型進行了廣泛和深入的研究,只有少數(shù)文獻闡述了關(guān)于F-M模型實現(xiàn)問題。比如,Alpay和Dubi[4]給出了直接構(gòu)建一個n維(n≥2)F-M實現(xiàn)方法。但是,這種方法求得F-M模型實現(xiàn)矩陣階數(shù)卻非常高[5]。
本文針對Alpay和Dubi[4]的n維F-M模型的實現(xiàn)矩陣求解方法,提出一種新解法,該方法獲取實現(xiàn)矩陣的階更低并且更易于實現(xiàn)。將此算法運用到多輸入多輸出(multiple input multiple output,MIMO)雷達系統(tǒng)中,提高了雷達目標(biāo)檢測[6-8]的效率以及反應(yīng)靈敏度。
對于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ù)矩陣,且A1,…,An∈Rr×r,B1,…,Bn∈Rr×l,C∈Rm×r,D∈Rm×l。
n維信號u(i1,…,in)的n維z變換的定義如下。
其中,z1,…,zn是單位延遲算子。
對式(1)進行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維多項式,當(dāng)m(0,…,0)≠0,則稱這個系統(tǒng)為因果系統(tǒng)[12]。如果一個有理矩陣的每一項都是因果的,則該有理矩陣是因果有理矩陣。
對于一個給定的傳遞矩陣H(z),如果式(4)成立,則A、B、C、D被稱為傳遞矩陣H(z)的F-M 狀態(tài)空間模型的實現(xiàn),該實現(xiàn)的充分必要條件就是H(z)是因果矩陣。由Galkowski教授[13]的研究可知,對于某個因果有理傳遞函數(shù)或傳遞矩陣,存在很多F-M模型實現(xiàn),但不一定具有相同的階次。由于狀態(tài)空間實現(xiàn)的順序與系統(tǒng)的計算或硬件實現(xiàn)的成本密切相關(guān),因此期望具有低階狀態(tài)空間實現(xiàn)。對于1維情況,已經(jīng)得到很好的發(fā)展。然而,對于n維(n≥2)情況,情況變得更加復(fù)雜和困難,因此,在n維情況下,尋求低階實現(xiàn)尤為重要。
本文將研究如下形式的傳遞矩陣的F-M 狀態(tài)空間模型的實現(xiàn)問題。
(6)
針對n維因果有理傳遞函數(shù),Alpay和Dubi[4]提出了一個構(gòu)建F-M模型實現(xiàn)的建設(shè)性方法。但是Alpay和Dubi的實現(xiàn)方法通常會產(chǎn)生不必要的高階次,這就促使本文研究一個替代的實現(xiàn)方法,考慮到給定傳遞函數(shù)或傳遞矩陣的結(jié)構(gòu)特征,階數(shù)更低的實現(xiàn)方法是可行的[1]。
如果F(z)是一個矩陣值函數(shù),本文用Μ(F(z))來表示Cm值函數(shù)的向量空間,其中值函數(shù)是F(z)所有列的線性組合。
定義1設(shè)M是Ω?Cn中的函數(shù)解析向量空間。如果相關(guān)的Gleason問題在Μ中是可解的,則稱其為解析不變量。比如,對任意的f(z)∈Μ,ε∈Ω,存在一個函數(shù)g1,ε(z),…,gn,ε(z)∈Μ滿足:
(7)
定義2設(shè)Μ是在Cn上定義的函數(shù)向量空間。若對任意f(z)∈Μ存在g1(z),…,gn(z)∈ζ滿足:
(8)
則稱ζ為M的后向移位空間。
定理1設(shè)H(z):Ω?Cn→Cm×l(Ω包含(0,…,0)),H(z)具有滿足式(4)的F-M實現(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)分別是分子和分母多項式。假設(shè)h(z)是因果關(guān)系,即p(0,…,0)≠0,將k={degp(z), degm(z)}稱為有理函數(shù)h(z)的度。Alpay和Dubi[4]所構(gòu)造空間如下:
(9)
它的維度是有限,并且還是Μ(h(z))的后向移位空間。可以看出通過Alpay和Dubi的方法得到實現(xiàn)矩陣的階數(shù)r與ζ中元素的數(shù)量相同,并且可以由式(10)計算得到。
(10)
因此,對于某個確定的n,Alpay和Dubi的方法所獲得的實現(xiàn)階數(shù)將完全由h(z)的k值確定?,F(xiàn)在容易驗證,即使對于相對較小的n和k,實現(xiàn)階數(shù)通常也相當(dāng)高,隨著n或k的增加,r的值也快速增加。表1顯示在n= 2,3,4,5和k= 2,3,…,6的情況下,通過Alpay和Dubi的方法獲得的實現(xiàn)階數(shù)。
表1 Alpay和Dubi的實現(xiàn)階數(shù)
但是,文獻[14,15]曾經(jīng)提出對于n維情況下的實現(xiàn)階數(shù)不僅取決于給定傳遞函數(shù)h(z)的k值,而且還取決于給定傳遞函數(shù)h(z)的結(jié)構(gòu)甚至系數(shù)。特別是在n維(n≥2)情況下,具有相同k值但結(jié)構(gòu)不同的傳遞函數(shù)可能具有不同實現(xiàn)階數(shù)。這個事實意味著當(dāng)h(z)缺少一些單項式時得到的實現(xiàn)階數(shù)遠低于Alpay和Dubi不考慮結(jié)構(gòu)特征的方法獲得的實現(xiàn)階數(shù)[16]。
例1
其中,z=(z1,z2,z3)是獨立變量,q20、q01、m10、m01為系數(shù),由于h(z)的度為3,因為m(0,0,0)=1≠0,顯然h(z)是因果的。
可以使用所有zα/p(z)構(gòu)造ζ,其中α=(α1,α2,α3),0≤|α|≤3。
由Alpay和Dubi的方法可得F-M實現(xiàn)矩陣階數(shù)r=20,與式(10)所得結(jié)果一樣。利用以下矩陣來代表的實現(xiàn)矩陣(A,B,C,D)。
D=0不僅滿足式(4)而且它的階次為7,遠低于Alpay和Dubi的方法計算出來的階次。本文所關(guān)注的是如何構(gòu)建這樣一個低階實現(xiàn),這將在下一節(jié)中討論。
步驟1首先取得每一列中系數(shù)非0的單項式zα=zα1…zαn,然后將這些單項式組成一個如下所示的矩陣[17,18]。
(11)
(12)
其中λT為ml取反后再加上1的式子的系數(shù)。重置Bij,即將第s個位置的項置1,其他的項保持不變。重置Aij,即第i行、j列的數(shù)值置為λT。
(13)
(14)
通過式(14),構(gòu)造的矩陣(A,B,C,D)即為H(z)的F-M實現(xiàn)。其中A=(A1…An),B=(B1…Bn)。
利用上述新方法對例1:
=q(z)/m(z)
進行實現(xiàn)矩陣的求解。
G1(z)=
h(z)-h(0)=
=G1(z)(z1B11+z2B12+z3B13)
=G1(z)(z1κ11+z2κ12+z3κ13)
進一步,令κj1=κj2=κj3=0。因為β12(z)=z1β11(z),按照步驟中的方法令κ21=1。
同理κ32(2)=κ43(2)=κ53(7)=κ61(5)=κ72(1)=1。剩余其他的項為0。Aj=A1j,Bj=B1j,j=1,2,3,G(z)=G1(z),C=G(0)=[1000000],D=h(0)=0。此時實現(xiàn)矩陣的階數(shù)為r=7,遠低于Alpay和Dubi教授所求出的20。
在雷達目標(biāo)跟蹤系統(tǒng)中,通常將被跟蹤的目標(biāo)建立如下狀態(tài)方程[19]:
X(k+1,k)=FX(k)+Qω(k-1,k-1)
(15)
而雷達的觀測方程可表示為
W(k)=h(k)+x(k)+ν
(16)
令k=z,上述狀態(tài)方程與觀測方程可以轉(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)空間模型來表示狀態(tài)方程與觀測方程[20],對MIMO雷達目標(biāo)檢測的研究也即是對F-M狀態(tài)空間模型的研究。
圖1是3個求得坐標(biāo)后的目標(biāo)位置信息。
圖1 待測目標(biāo)位置信息
通過定位算法進行坐標(biāo)計算后,3個位置的坐標(biāo)分別用R=(xr,yr),S=(xs,ys),W=(xw,yw)表示。則可對系統(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)進行z變換之后,可得系統(tǒng)傳遞函數(shù)為
(19)
成功實現(xiàn)建模后,對MIMO雷達系統(tǒng)的研究就轉(zhuǎn)移到F-M 狀態(tài)空間模型上。
例2對MIMO雷達系統(tǒng)進行系統(tǒng)建模后得到的3×2傳遞矩陣如下所示[20]。
m0(z)=1-d01z2-d02z3+d03z1z2
m2(z)=1+d11z1-d12z1z2+d13z1z2z3
通過新方法可以得到以下結(jié)果:
因為C1=NHT1,N1=NHT1Ψ1,可得:
代入式(13)可以得到:
同理,另一結(jié)果如下:
C2=NHT2,N2=NHT2Ψ2
代入式(13)得到如下結(jié)果:
由上述結(jié)果知,傳遞矩陣H(z)的實現(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所提出方法求得的實現(xiàn)矩陣階數(shù)是45,而新方法所求階數(shù)為10。
將此算法運用到MIMO雷達系統(tǒng)的多維處理過程中,采用Matlab工具進行仿真實驗。圖2與圖3分別是2種方法在x軸與y軸上速度協(xié)方差的比較。Alpay和Dubi算法的檢測速度協(xié)方差在收斂前一直稍大于F-M狀態(tài)空間模型的步數(shù),并且后者的收斂速度快于前者。這表明采用新方法進行降階后,可以提高雷達系統(tǒng)的定位準(zhǔn)確度以及反應(yīng)速度。
圖2 x方向預(yù)測和更新速度誤差協(xié)方差
圖3 y方向預(yù)測和更新速度誤差協(xié)方差
本文針對n維F-M模型實現(xiàn)問題提出了一種新的求解方法,與Alpay和Dubi的方法相比,它可以得到更低階實現(xiàn)矩陣。具體而言,即使對于簡單的有理傳遞函數(shù),Alpay和Dubi給出的方法通常會產(chǎn)生具有不必要的高階實現(xiàn)矩陣。并且,通過考慮給定傳遞函數(shù)或傳遞矩陣的結(jié)構(gòu)性質(zhì),可以構(gòu)造一個階數(shù)較低的實現(xiàn)矩陣,并通過舉例說明該方法的可實施性及有效性。將此算法運用到MIMO雷達系統(tǒng)的多維處理過程中,可以大幅降低系統(tǒng)的數(shù)據(jù)處理時間,提高雷達系統(tǒng)對所跟蹤物體的反應(yīng)速度和定位的準(zhǔn)確度。今后更進一步的研究是盡可能減少約束條件,使之能夠獲取到更加低階的系統(tǒng)實現(xiàn)矩陣。