多元随机过程的建模与谱估计

第七章 多元随机过程的建模与谱估计

7.1 多元随机过程的表示

l维平稳随机向量过程Y(n)由l个平稳随机过程构成

Y(n)[yT1(n),y2(n),,yl(n)] 其二阶性能由均值向量

YEY(n)[y,y,,T12yl] 和协方差矩阵

Cm)E[Y(n)Y(Y][Y(nm)TY]Cyiyj(m)ll 决定，其中

Cyiyj(m)是随机过程yi(n)和

yj(n)的协方差，即

Cyiyj(m)E[yi(n)yi][yj(n)yj]，1il,1jll维白噪声向量W(n)：

(7-1)

(7-2)

(7-3)

(7-4)

,m0W0,CW(m)0,m0 (7-5) 其中为常数矩阵。若白噪声向量W(n)的个分量互不相关，则其协方差矩阵是对角矩阵，即

222diag[w,,,w2wl]1 (7-6)

一般，我们总是将原随机向量减去均值向量估计，构成一个零均值的、新的随机向量。然后对新的随机向量进行各种分析。不失一般性，我们总是可以假设随机向量是零均值的。此时，协方差矩阵CY(m)与互相关矩阵RY(m)相同，即是

llRyiyj(m)CY(m)EY(n)YT(nm)Eyi(n)yj(nm)其中，

Ryiyj(m)llRY(m) (7-7)

是随机过程yi(n)和

yj(n)的互相关函数。

性质：１）

RY(m)RYT(m) (7-8)

【证明：显然，

Ryi,yj(m)Ryj,yi(m)，得

TRY(m){Ryi,yj(m)}ll{Ryj,yi(m)}ll{Ryi,yj(m)}TllRY(m)

】

2）RY(0)是非负定的

【证明：用不全为零的实数ai，作随机过程

l z(n)ai•yi(n)aT•Y(n)i1

则有

RZ(0)E{z2(n)}E{aT•Y(n)•YT(n)•a}aTE{Y(n)•YT(n)}•aaT•RY(0)•a0

当且仅当Y(n)0,RZ(0)0。

】

7.1 向量过程的模型表示与谱

①向量过程的ＡＲ模型与功率谱

用ｍ维AR(p)过程模型描述的随机向量过程表示为

Y(n)A1•Y(n1)Ap•Y(np)W(n)(7-9)

其中Ｙ(k)是ｍ维向量过程，Ａi 为ｍ×ｍ阶参数矩阵，Ｗ(k)是ｍ维白噪声。

记

A(p,z1)IA1•z1Ap•zpH(p,z1)[A(p,z1)]1, (7-10)

则(1)式可改写为

A(p,z1)•Y(n)W(n) 或 Y(n)H(p,z1)•W(n) 向量过程｛Ｙ(ｋ)｝的功率谱密度函数矩阵为

SY(ej)[H(p,ej)]•SW(ej)•HT(p,ej)[A(p,ej)]1•SjW(e)•[A(p,ej)]T 其中

SW(ej)是常数矩阵．当Ｗ(k)的各分量互不相关时，Ｓｗ (ｅｊω )是对角矩阵

Sj)diag[22W(e1,2,2l] ②向量过程的ARMA模型与功率谱

用ｍ维ARMA(p,q)过程模型描述的随机向量过程表示为

(7-11)

(7-12)

(7-13)

Y(n)A1•Y(n1)Ap•Y(np)Bq•W(nq) B0•W(n)B1•W(n1) (7-14)

其中Y(k)是ｍ维向量，Ｗ(k)是ｍ维白噪声，

Ai,Bj为ｍ×ｍ阶参数矩阵。

记

A(p,z1)IA1•z1Ap•zp11qB(q,z)B0B1•zBq•z1111H(p,q,z)[A(p,z)]•B(q,z)  (7-15)

则(7-11)式可改写为

A(p,z1)•Y(n)B(q,z1)•W(n) 或 Y(n)H(p,q,z1)•W(n) (7-16)

向量过程Y(k)的功率谱密度函数矩阵为

SY(ej)[H(p,q,ej)]•SW(ej)•H(p,q,ej)1jjjj1jA(p,e)•[B(q,e)]•S(e)•B(q,e)•A(p,e) (7-17) W

其中

SW(ej)是ｍ维白噪声的协方差矩阵。

显然，如果我们获得了过程模型参数及ｍ维白噪声的协方差矩阵Γ的估计，也就获得了过程功率谱的估计。 前面讨论的标量过程的ＡＲ、ＡＲＭＡ建模与谱估计可以推广到多变量过程。

7.2 向量ＡＲ过程的建模

１ AR(p)过程的YW方程

TY对(1)式右乘(nm)，并取期望

得

由因果性质,得

于是

展开

E{[Y(n)AT1•Y(n1)Ap•Y(np)]•YT(nm)}E{W(n)•Y(nm)}

RY(m)A1•RY(m1)Ap•RY(mp)RWY(m) R(m)QW,m0

XY0,m0 RY(m)A1•RY(m1)A•Rp)QW,m0pY(m0,m0(7-18)

(7-19)

(7-20)

RY(0)A1•RY(1)Ap•RY(p)QWRY(1)A1•RY(0)Ap•RY(p1)0R(p)A•R(p1)A•R(0)01YpYY (7-21)

对(7-21)式求转置，并考虑到相关性质

R)RTY(mY(m)，则式(7－21)可改写为

RTTTY(0)RY(1)•A1RY(p)•ATpQWR(0)•ATTp1)•ATY(1)RY1RY(p0RY(p)R(p1)•AT1R(0)•ATYYp0 (7-22)为向量AR(p)过程YW方程，令

RTY(0)RYRT(1)Y(p)IQRRTRTY,p(0)Y(1)RY(0)ATWY(p1)1p0pRY(p)RY(p1)R(0)TY，A，

0p 则(7-22)式可改写为矩阵形式

RY,p(0)•pp 1.

性质：互相关矩阵RY,p(0)是非负定的。

(7-22)

(7－23) (7－24)

【证明：

作随机过程

其中，

aTiZ(k)aT•Y(ki)ii1n

是实数向量，ai（ｉ＝１，２，，ｎ）不全为零。则 ·

n

nTEZ(k)E(ai•Y(ki))2i1nTEaj•Y(kj)•YT(ki)•aii,j12i,j1naanTj•EY(kj)•YT(ki)•ai•RY(j1)•aiTji,j1aT•RY(j1)•a0

其中，

Taa1Ta2aTnT

RY(1)RY(0)R(1)RY(0)YRY(n1)RY(n1)RY(n2)RY(n1)RY(n2)RY(0)

】

2. 向量过程ＡＲ建模的互相关矩阵法

Ｐ阶最佳线性预测

Y(n)A1•Y(n1)Ap•Y(np) (7-25)

使预测误差

ep(n)Y(n)Y(n)Y(n)A1•Y(n1)Ap•Y(np) (7-26)

的互相关矩阵

Qe,pEep(n)•eTp(n)

E{[Y(n)A1•Y(n1)Ap•Y(np)]•[Y(n)A1•Y(n1)Ap•Y(np)]T} (7-29)

最小。又记

hp(n)[YT(n)YT(n1)YT(np)]T (7-28)

则有

Tep(n)[hT(n)•]pp (7-29)

并且

因此

令

E{[hTn)•hTp((n)]}Y(n)E•[YT(n)YT(n1)YT(np)]Y(np)

RY(0)RY(1)RY(p)RY(1)RY(0)RY(p1)RY,p(0)RY(p)RY(p1)R(0)Y QTT•[hTe,pE{ep(n),eTp(n)}E{[hp(n)•p]p(n)•p]}T

TpE{[hp(n)•hTp(n)]}•pp•RY,p(0)•p (7-30)

(7-31)

A1TRY(1)QpPTApRY(p) (7-32) ,

则有：

IRY(0)TPpRY,p(0)p ，

PRY1,p(0) 于是

QTRTY(0)Ile,pIlp•P•PRY1,p(0)p 定理：随机过程的p阶最佳前向线性预测参数矩阵

p满足

RY1,p(0)•pP0 并且，最佳前向线性预测误差的方差为

QRTe,pY(0)RY(1)•AT1RTY(p)•ATp 【证明：由(7-34)得

(7-33) (7-34)

(7-35)

(7-36)

TQe,p[RY(0)p•Pp•RY,p1(0)]•IlpRY(0)T•PT•pT•RY,p1(0)•ppPpTT1TTTTRY(0)T•R•(RY,p1(0)•p)•p•RY,p1(0)•PY1,p(0)•PPPpPP1

令

CpRY(0)T•RY1,p(0)•PP1 (7-37)

则

1TQe,pCpT•[RY,p1(0)•p]•RY,p1(0)•[RY,p1(0)•pP]PpPCp[pP•RY,p1(0)]•[RY,p1(0)•pP]Cp[p•RY,p1(0)P]•RY,p1(0)•[RY,p1(0)•pP]1TT1

Cp[RY,p1(0)•pP]•RY,p1(0)•[RY,p1(0)•pP]1 (7-38)

显然，(7-37)中第一项因此，(7-35)得证。

Cp与参数

p无关，第二项是一个非负得二次型，当第二项为零时，

Qe,p达最小值。

最佳前向线性预测误差得方差

TQe,pCpRY(0)T•R(0)•pY1,p(0)•PRYPP1

将(7-32)代入上式，即可得(7-36)式。

合并(7-35)和(7-36)式，得

RY,p(0)•pp (7-39)

其中

TpQe,p00 结论：P阶最佳线性预测参数满足多变量YK方程．

依阶次递推：

按上述同样的方法可推，Ｐ＋１阶最佳线性预测参数满足P1阶YK方程：

RY,p(0)•p1P10Qe,p1RY(0)TP1•p1 记

(7-40)

(7-41)

TRY(p)pP,invTRY(1) (7-42)

则RY,p(0)又可以分块表示为

令

利用分块矩阵求逆公式可得

并且由(7-32)和(7-44)式，得

RRY,p1(0)p(0)pY,TpRY(0) GR111Y,p1(0)G112RY,p1(0)•pG1122[RY(0)Tp•RY,p1(0)•p] R1G11G12•G22•GT12G12•G22y,P(0)GT22•G12G22 GT112•pTp•RY,p1(0)•pTp•p (7-43)

(7-44)

(7-45)

(7-46)

于是，由(7-41)

记

则

G112RY,p1(0)•p

并且，预测方差

p1Tp•pGT11•pG12•G22•G12•pG12•G22•RY(p1)G22•GT12•pG22•RY(p1)pGT12•G22•p•pG12•G22•RY(p1)G•T22p•pG22•RY(p1)

pG12•G22•Tp1•pG•T22p1•p  pp10p1ll，Tp1Tp1•p G12•G22•Tp1p1G12GT•G22•Tp122•p1Ill(7-47)

(7-48)

(7-49)

Qp1RY(0)T•p1p1

RTp•pTY(0)p1•p1

QpTp1•p1 (7-50)

而根据(7-49)、(7-46)式

TTp1•pR1)•G112pI•GTY(p22•p1ll

Tp•G12RY(p1)•G22•Tp1Tp•pRY(p1)•G22•Tp1

Tp•p•G22•Tp1p1•G22•Tp1 (7-51)

将(7-51)代入(7-50)，得

QQTp1•p1QTp1ppp1•G22•p1Qp特例：当l1时，

G1T22Df,p，

G22•p1p1，G12p,inv，

p1p1•p,inv

实际应用中，互相关函数RY(m)的估计可由下式计算

(7-52)

1RY(m)N

nNm1n0Y(n)•YT(nm)，m1,2,P (7-53)

综合上述，有如下依阶次递推的算法。

第一步：初始化

按照计算式(7-53)估计互相关矩阵

nRY(0)1NN1Y(n)•YT(nm)n0，R11Y,0(0)RY(0) 第二步：计算m1

1nN2RY(1)NY(n)•YT(n1)n0，11RY(1) 1R1(0)•TTY,01，Qe,1RY(0)RY(1)•1 第三步：计算中间变量

(7-54)

(7-55)

(7-56)

GR1(0)Y,p111G12G11•pT1G[R(0)•G]Y12p221nNp2RY(p1)Y(n)•YT(np1)Nn0TTp1RY(p1)pTp1Tp1•pQp1Qpp1•G22•Tp1

第四步：判阶，即递推终止条件判断

如果

QpQp1/Q0，结束递推；否则，到第五步；

第五步：参数递推

QGQp1p12•G22•Tp10Ill

第五步：逆矩阵递推

G11G12•G22•G12TRy,P(0)TG•G2212

1G12•G22G22

PP1；

转到第三步。

特点：运算量较大；离线计算。

7.3 向量过程AR建模的RLS方法

问题：根据随机向量过程的观测数据Y(n)，建立AR模型，使前向预测误差

ep(n)Y(n)A1Y(n1)ApY(np) (7-57)

的方差

Me,pEep(n)Eep(n)•eTp(n)Eep(n)min 1) 当

Eep(n)0时，

Me,pEep(n)•eTp(n)Qe,p

2) 当

ep(n)ep,1,ep,2,,ep,l中，各分量互不相关时，即

Eep,i(n)•ep,j(n)Eep,i(n)•Eep,j(n)0,ij 时，有

Q2e,pdiagEe2p,1,Eep,2,,Ee2p,l 使

Qe,p最小，等价于使

ep(n)的每个分量的方差最小，即

(7-58)

(7-59)

(7-60)

Jp,ie2p,i(n)min,ijnpN1 (7-61)

等价于使

TJpe2p,iep(n)•ep(n)minnpN1 (7-62)

当观测数据为有限序列Y(0),Y(1),,Y(N1)，Nl(p1)，并且不加数据窗时，指标函数可以改写为

Jp(p,N1)e(n)•ep(n)e2p,i(n)Tpnpnpi1N1N1l (7-63)

其中

ep(n)[ep,1,ep,2,,ep,l]T (7-)

由于

TTTeTp(n)Y(n)Y(n1)•A1YT(np)•ATp (7-65)

记

p[Ap,Ap1,,A1]T (7-66)

pTp,Il (7-67)

Thp(n)YT(np)YT(np1)YT(n1) (7-68)

TTdp(n)hTp(n)Y(n) (7-69)

T则 HTp(p,N1)hp(p)hp(p1)hp(N1) DTp(p,N1)dp(p)dp(p1)dp(N1) e(p,N1)e1)Tpp(p)ep(p1)ep(N Y(p,N1)Y(p)Y(p1)Y(N1)T eT(n)hTTTpp(n)•pY(n)dp(n)•p ep(p,N1)Hp(p,N1)•pY(p,N1)

Dp(p,N1)•p (7-70)

(7-71) (7-72)

(7-73)

(7-74)

(7-75)

lpl(p1)lllp1(Np)lpl(p1)1(Np)l(p1)其中

pR，

pR，

hp(n)R，

Hp(p,N1)R，

dp(n)R，

Dp(p,N1)R(p,N1)R(Np)l，Y(p,N1)R(Np)l。

1) 正则方程求解

由(7-63)式得

J,N1)e1)2Y(p,N1)H2p(pp(p,nFp(p,N1)•pF

lY1)H2i(p,Np(p,N1)•p,ii02 其中

Yi(p,N1)是Y(p,N1)的第i列，p,i是p的第i列，即

Y(p,N1)Y1(p,N1)Y2(p,N1)Yl(p,N1)

pp,1p,2p,l

由指标

Jp(p,N1)取极小的必要条件，得

HTp(p,N1)Hp(p,N1)Tp,iHp(p,N1)Yi(p,N1)，i1,2,l 将上式中l个矩阵方程合并，得

，

(7-76)

(7-77)

ep

THTp(p,N1)•Hp(p,N1)pHp(p,N1)•Y(p,N1) (7-78)

（7-78）式为最小二乘正则方程，最小二乘解为

1TˆpHp(p,N1)•Hp(p,N1)HTp(p,N1)•Y(p,N1) (7-79)

利用最小二乘递推算法(RLS)可以避免矩阵求逆。

2) 正交变换求解

根据(7-76)式，有

J1)e2p,N1)•2

p(p,Np(p,n1)FY(p,N1)Hp(pF

TraceTY(p,N1)Hp(p,N1)•pY(p,N1)Hp(p,n1)•p

TraceTTp•Dp(p,N1)•Dp(p,N1)•p 对于信息矩阵

Dp(p,N1)，存在正交变换矩阵T，使其上三角化，即

T•Dp,N1)Rp(p,N1)p(0

(7-80)

其中，

Rp(p,N1)是(p1)l(p1)l阶上三角阵。于是有

TTTJp(p,N1)TraceD(p,N1)•TT•Dp(p,N1)ppp

TTTrace•Rpp(p,N1)•Rp(p,N1)•p

Rp(p,N1)•p2F (7-81)

对上三角阵

Rp(p,N1)进行块划分

'Rp(p,N1)G1(p,N1)Rp(p,N1)0G2(p,N1) (7-82)

其中，

'Rp(p,N1)RplplpllllG(p,N1)RG(p,N1)R12，，。将(7-82)代入(7-81)式，并考虑到F-范

数的性质，(7-81)式可改写为

Jp(p,N1)R(p,N1)pG1(p,N1)'p2FG2(p,N1)2F (7-83)

显然，(23)式中的第二项与参数Ｓ解为

p无关；当第一项为零时，指标函数

Jp(p,N1)到达最小值。因此，Ｌ

ˆR'1(p,N1)•G(p,N1)p1p2ˆJ(p,N1)G(p,N1)2pF (7-84)

正交变换求解也可以利用递推算法。

7.5 多元ARMA过程的建模

1) 因果性和可逆性

A(p,z1)•Y(n)B(q,z1)•W(n) (7-85)

稳定性（因果性）判据

11detA(p,z)0A(p,z)的根都在z-平面单位元内，|z|1若对所有，，即ARMA过程得极点多项式矩阵

则(7-85)式具有唯一的平稳解

Y(n)H(z)•W(n)Hi•W(ni)1i0 (7-86)

其中，矩阵Hi由下式唯一确定

H(z)HiziA1(p,z1)•B(q,z1)1i0 (7-87)

用A1(p,z1)左乘(25)式，并比较等式两边zi的系数，得

H0IiHiBi1iqAjHij，j1pHiAjHij，iqj1 多元ARMA(p,q)过程可以用AR()过程表示。

最小相位（可逆性）判据

若对所有|z|1，

detB(q,z1)0，且Y(n)是(7-85)式的平稳解，则 W(n)F(z1)•Y(n)Fi•Y(ni)i0 其中矩阵Ｆi 由下式唯一确定

F(z1)F11iziB(q,z1)•A(p,z)i0用

B1(q,z1)左乘(7-90)式，并比较等式两边zi的系数，得 (7-88)

(7-)

(7-90)

F0Imin(i,q)1ipFiAiBjFij，j1qFiBjFij，ipj1 (7-91)

1F(z)可能是有限阶次的，这是因为 【注：在向量情况下，

11F(z)B(q,z)A(p,z1)1111B(z)A(p,z)1detB(q,z)

111detB(q,z)为行列式。由于有限阶次矩阵多项式的伴随矩阵也是B(z)B(q,z)其中，是的伴随矩阵，1111detB(q,z)B(q,z)F(z)也有限阶数的。因此，当行列式为一常数矩阵时，矩阵多项式也是有限阶数的；

是有限阶数的。

】

2) 向量ARMA过程模型估计的LS两步法

根据随机向量过程的观测数据Y(0),Y(1),建立ARMA模型是一非线性问题。然而与单变量AMAR过程相同的是，如果向量ARMA过程是最小相位(可逆)的，则可表示为等价的无限阶的向量AR过程。因此，我们可以直接将单变量过程ARMA建模的ＬＳ两步法推广到向量过程，具体做法如下所述。

第一步，利用观测数据Y(0),Y(1),建立高阶的AR模型

Y(n)F1•Y(n1)FM•Y(nM)W(n) (7-92)

第二步，利用高阶AR模型白化观测数据Y(0),Y(1),，获得W(n)的估计值

ˆ W(n)Y(n)F•1Y(n1)FM•Y(nM) (7-93)

该估计实际上就是最优预测误差。

W(0),W(1),拟合方程 第三步，利用Y(0),Y(1),和Y(n)A1Y(n1)ApY(np)W(n)B1W(n1)BqW(nq)

使拟合误差

ep,q(n)Y(n)A1•Y(n1)Ap•Y(np)W(n)B1•W(n1)Bq•W(nq) (7-94)

的(方差)二次型指标

Jp,qEep,q(n)22EeTp,q(n)•ep,q(n) (7-95)

最小。

当观测数据为有限序列Y(0),Y(1),,Y(N1)，N(p1)，并且不加数据窗时，指标函数可以改写为

Jp,q(p,N1)eTp,q(n)•ep,q(n)npN1 (7-96)

利用最小二乘方法，即可求解(7-96)式。

3) 向量ARMA过程模型估计的交叉相乘法

第一步，利用观测数据Y(0),Y(1),，按(7-92)式建立高阶的AR模型。其中，阶次满足Mpq。

11F(M,z)F(z)的泰勒级数展开式的前M项的估计，第二步，由于是可逆性矩阵方程(30)式中矩阵多项式

根据可逆性方程，近似有

F(M,z1)B1(q,z1)A(p,z1) (7-97)

11A(p,z)，具体过程如下。 B(q,z)因此，可根据(7-97)式确定过程的零．极点多项式矩阵和

由（7-97）式得

111A(p,z)B(q,z)•F(M,z) (7-98)

即

IA1z1Apzp

(IB1z1Bqzq)•(IF1z1FMzM) (7-99)

比较(7-99)式两边zi的系数，得

min(i,q)AiBj•Fij,0ipj0 qBj•Fij,piMj0 根据(7-101)式，有

FTpFTp1FTTp1qBT1FFTTpFpFTp1TTp2qB2Fp2•FTM1FTM2FTBTTMqqFM 当Mpq时，(7-102)式有唯一解；否则，有最小二乘解。又根据(40)式，有

AT1FT1I0BT1ATqFTFTI•BT2q1qATpFTpFTp1FTBTqqp (7-100)

(7-101)

(7-102) (7-103)