一种实现Fornasini-Marchesini 模型的新方法的研究①

曹中泳 程 骅 吴徐冬子 刘 昶

(武汉科技大学信息科学与工程学院 武汉 430080)

0 引 言

在过去的几十年中，多维系统吸引了极大的关注，这主要归功于多维系统在各种工程领域的广泛使用和潜在使用，如信号和图像处理、热工艺、医疗应用、无线传感器网络等[1-3]。

多维系统理论中一个基本问题是通过确定的多维状态空间模型来实现给定的有理传递函数或传递矩阵，该模型通常是Roesser模型或Fornasini-Marchesini(F-M)模型。与传统的1维(1-D)情况不同，通常n维(n≥2)滤波器或者系统的最小状态空间实现很难获得。因此研究可以得到低阶n维状态空间的实现新方法就显得格外重要。此外，许多研究人员对Roesser状态空间模型进行了广泛和深入的研究，只有少数文献阐述了关于F-M模型实现问题。比如，Alpay和Dubi[4]给出了直接构建一个n维(n≥2)F-M实现方法。但是，这种方法求得F-M模型实现矩阵阶数却非常高[5]。

本文针对Alpay和Dubi[4]的n维F-M模型的实现矩阵求解方法，提出一种新解法，该方法获取实现矩阵的阶更低并且更易于实现。将此算法运用到多输入多输出(multiple input multiple output，MIMO)雷达系统中，提高了雷达目标检测[6-8]的效率以及反应灵敏度。

1 F-M状态空间模型

对于n维MIMO线性离散系统，F-M状态空间模型的描述如下[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分别是局部状态向量、输入向量、输出向量。A、B、C、D是实数矩阵，且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}。 若系统的传递函数h(z)=q(z)/m(z)，其中q(z)和m(z)是n维多项式，当m(0,…,0)≠0，则称这个系统为因果系统[12]。如果一个有理矩阵的每一项都是因果的，则该有理矩阵是因果有理矩阵。

对于一个给定的传递矩阵H(z)，如果式(4)成立，则A、B、C、D被称为传递矩阵H(z)的F-M 状态空间模型的实现,该实现的充分必要条件就是H(z)是因果矩阵。由Galkowski教授[13]的研究可知，对于某个因果有理传递函数或传递矩阵，存在很多F-M模型实现，但不一定具有相同的阶次。由于状态空间实现的顺序与系统的计算或硬件实现的成本密切相关，因此期望具有低阶状态空间实现。对于1维情况，已经得到很好的发展。然而，对于n维(n≥2)情况，情况变得更加复杂和困难，因此，在n维情况下，寻求低阶实现尤为重要。

本文将研究如下形式的传递矩阵的F-M 状态空间模型的实现问题。

(6)

2 Alpay和Dubi的实现方法

针对n维因果有理传递函数，Alpay和Dubi[4]提出了一个构建F-M模型实现的建设性方法。但是Alpay和Dubi的实现方法通常会产生不必要的高阶次，这就促使本文研究一个替代的实现方法，考虑到给定传递函数或传递矩阵的结构特征,阶数更低的实现方法是可行的[1]。

如果F(z)是一个矩阵值函数，本文用Μ(F(z))来表示Cm值函数的向量空间，其中值函数是F(z)所有列的线性组合。

定义1设M是Ω⊂Cn中的函数解析向量空间。如果相关的Gleason问题在Μ中是可解的，则称其为解析不变量。比如，对任意的f(z)∈Μ，ε∈Ω，存在一个函数g1,ε(z),…,gn,ε(z)∈Μ满足：

(7)

定义2设Μ是在Cn上定义的函数向量空间。若对任意f(z)∈Μ存在g1(z),…,gn(z)∈ζ满足：

(8)

则称ζ为M的后向移位空间。

定理1设H(z)：Ω⊂Cn→Cm×l(Ω包含(0,…,0))，H(z)具有满足式(4)的F-M实现矩阵(A，B，C，D)，当且仅当M(H(z))具有有限维度和解析不变后移空间ζ。

设h(z)=q(z)/m(z)是单输入单输出(SISO)的多维系统有理传递函数，q(z)和m(z)分别是分子和分母多项式。假设h(z)是因果关系，即p(0,…，0)≠0，将k={degp(z), degm(z)}称为有理函数h(z)的度。Alpay和Dubi[4]所构造空间如下：

(9)

它的维度是有限，并且还是Μ(h(z))的后向移位空间。可以看出通过Alpay和Dubi的方法得到实现矩阵的阶数r与ζ中元素的数量相同，并且可以由式(10)计算得到。

(10)

因此，对于某个确定的n，Alpay和Dubi的方法所获得的实现阶数将完全由h(z)的k值确定。现在容易验证，即使对于相对较小的n和k，实现阶数通常也相当高，随着n或k的增加，r的值也快速增加。表1显示在n= 2,3,4，5和k= 2，3，…，6的情况下，通过Alpay和Dubi的方法获得的实现阶数。

表1 Alpay和Dubi的实现阶数

但是，文献[14，15]曾经提出对于n维情况下的实现阶数不仅取决于给定传递函数h(z)的k值，而且还取决于给定传递函数h(z)的结构甚至系数。特别是在n维(n≥2)情况下，具有相同k值但结构不同的传递函数可能具有不同实现阶数。这个事实意味着当h(z)缺少一些单项式时得到的实现阶数远低于Alpay和Dubi不考虑结构特征的方法获得的实现阶数[16]。

例1

其中，z=(z1,z2,z3)是独立变量，q20、q01、m10、m01为系数，由于h(z)的度为3，因为m(0,0,0)=1≠0，显然h(z)是因果的。

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

由Alpay和Dubi的方法可得F-M实现矩阵阶数r=20，与式(10)所得结果一样。利用以下矩阵来代表的实现矩阵(A,B,C,D)。

D=0不仅满足式(4)而且它的阶次为7，远低于Alpay和Dubi的方法计算出来的阶次。本文所关注的是如何构建这样一个低阶实现，这将在下一节中讨论。

3 新方法步骤

步骤1首先取得每一列中系数非0的单项式zα=zα1…zαn，然后将这些单项式组成一个如下所示的矩阵[17，18]。

(11)

(12)

其中λT为ml取反后再加上1的式子的系数。重置Bij，即将第s个位置的项置1，其他的项保持不变。重置Aij，即第i行、j列的数值置为λT。

(13)

(14)

通过式(14)，构造的矩阵(A,B,C,D)即为H(z)的F-M实现。其中A=(A1…An)，B=(B1…Bn)。

4 实 例

利用上述新方法对例1：

=q(z)/m(z)

进行实现矩阵的求解。

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。此时实现矩阵的阶数为r=7，远低于Alpay和Dubi教授所求出的20。

在雷达目标跟踪系统中，通常将被跟踪的目标建立如下状态方程[19]：

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

(15)

而雷达的观测方程可表示为

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

(16)

令k=z，上述状态方程与观测方程可以转化为

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

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

(17)

若系统矩阵v=0，则系统是严格因果系统。

基于上述公式,利用 F-M状态空间模型来表示状态方程与观测方程[20]，对MIMO雷达目标检测的研究也即是对F-M状态空间模型的研究。

图1是3个求得坐标后的目标位置信息。

图1 待测目标位置信息

通过定位算法进行坐标计算后，3个位置的坐标分别用R=(xr,yr),S=(xs,ys),W=(xw,yw)表示。则可对系统建立如下3维F-M状态空间模型：

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变换之后，可得系统传递函数为

(19)

成功实现建模后，对MIMO雷达系统的研究就转移到F-M 状态空间模型上。

例2对MIMO雷达系统进行系统建模后得到的3×2传递矩阵如下所示[20]。

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

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

通过新方法可以得到以下结果：

因为C1=NHT1,N1=NHT1Ψ1，可得：

代入式(13)可以得到：

同理，另一结果如下：

C2=NHT2,N2=NHT2Ψ2

代入式(13)得到如下结果：

由上述结果知，传递矩阵H(z)的实现矩阵如下。

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所提出方法求得的实现矩阵阶数是45，而新方法所求阶数为10。

将此算法运用到MIMO雷达系统的多维处理过程中，采用Matlab工具进行仿真实验。图2与图3分别是2种方法在x轴与y轴上速度协方差的比较。Alpay和Dubi算法的检测速度协方差在收敛前一直稍大于F-M状态空间模型的步数，并且后者的收敛速度快于前者。这表明采用新方法进行降阶后，可以提高雷达系统的定位准确度以及反应速度。

图2 x方向预测和更新速度误差协方差

图3 y方向预测和更新速度误差协方差

5 结 论

本文针对n维F-M模型实现问题提出了一种新的求解方法，与Alpay和Dubi的方法相比，它可以得到更低阶实现矩阵。具体而言，即使对于简单的有理传递函数，Alpay和Dubi给出的方法通常会产生具有不必要的高阶实现矩阵。并且，通过考虑给定传递函数或传递矩阵的结构性质，可以构造一个阶数较低的实现矩阵，并通过举例说明该方法的可实施性及有效性。将此算法运用到MIMO雷达系统的多维处理过程中，可以大幅降低系统的数据处理时间，提高雷达系统对所跟踪物体的反应速度和定位的准确度。今后更进一步的研究是尽可能减少约束条件，使之能够获取到更加低阶的系统实现矩阵。