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

1、第七章 多元随机过程的建模与谱估计7.1 多元随机过程的表示 l维平稳随机向量过程 Y(n)由l 个平稳随机过程构成Y(n) y1(n), y2(n), ,yl(n)T(7-1)其二阶特性由均值向量Y ：Y E Y(n) y1 , y2 ,ylT(7-2)和协方差矩阵 CY (m) ：CY(m) E Y(n) YY(n m) Y TCy1y1(m)C y2 y1 ( m)Cy1y2 (m)Cy2y2(m)Cy1yl (m)Cy2yl (m)(7-3)Cyl y1(m) 决定，其中 Cyy (m)是随机过程 yi(n)和 y j ( n)的协方差，即Cyiyj(m) E yi(n)yi yj(n

2、 m)yj ，1 i l,1Cyly2 (m)Cylyl (m)jl由于Cyiyj(m) Ryi yj (m) yi yj ，1 i l,1 j l 因此，协方差矩阵 CY (m) 又可表示为(7-4)CY (m) RY (m) Y Y其中， RY(m)为l维平稳随机向量过程 Y(n)的自相关矩阵。该矩阵中的第 i 行第 j列元素是随机过程 yi(n)和 y j (n)的互相关函数 Ryy (m)，即Ry1y1(m)Ry1y2 (m)Ry1yl (m)Ry2y1(m)Ry2y2 (m)Ry2 yl (m)Ryl y1(m)Ryly2 (m)Rylyl (m)RY (m)ll(7-5)当Y(n)

3、 的均值为零时，协方差矩阵 CY (m)与互相关矩阵 RY(m)相等。一般情况下，总是将随机向量 减去其均值向量估计，构成一个零均值的、新的随机向量。然后对新的随机向量进行各种分析。举例， l 维白噪声向量 W(n) 的二阶特征量为：QW ,m 0W 0,CW (m)0, m 0其中 QW 为常数矩阵。若白噪声向量 W (n)的个分量互不相关，则其协方差矩阵QW 是对角矩阵，即2 2 2(7-6)(7-7)QW diag w21, w22, , w2l 互相关矩阵性质：)RY (m) RYT ( m)证明：因为， Ryiyj (m) E yi(n)yj(n m) E yj ( n) yi (n

4、 m) Ryjyi ( m)，所以 RY(m) Ryiyj (m) l l Ryjyi( m)l l Ryiyj( m) lT l RYT( m)】2) RY(0) 是非负定的【证明：用 l个不全为零的实数 ai，i 1,2, ,l ，作随机过程lz(n)ai yi(n) aT Y(n)i1aT a1,a2, ,al ，则有Rz(0) Ez2(n) EaTY(n)YT(n)a aTEY(n)YT(n)a aT RY(0) a 0当 Y(n) 没有零分量时， Rz (0) 正定。】7.2 向量过程的模型表示与谱向量过程的模型与功率谱用l 维 AR( p)过程模型描述的随机向量过程 Y(n)表示为

5、Y(n) A1Y(n 1)ApY(n p) W(n) (7-8)其中 Ai(i 1,2, ,p)为l l 阶参数矩阵， W(n)是l维白噪声。记 A(p,z1) I A1z1Apz p，H(p,z1) A(p,z1) 1 (7-9)则(7-8)式可改写为11 A(p,z1)Y(n) W(n) 或 Y(n) H(p,z 1)W(n)(7-10)随机向量过程 Y(n) 的功率谱密度函数矩阵为SY(ej ) H(p,ej )SW(ej )HT(p,e j ) A(p,ej )1SW(ej )A(p,e j ) T (7-11) 其中 SW(ej ) 是常数矩阵。当 W (n)的各分量互不相关时， S

6、W(ej )是对角矩阵，即SW(ej ) diag 12, 22,l2(7-12)向量过程的 ARMA模型与功率谱 用l 维 ARMA( p,q)过程模型描述的随机向量过程表示为Y(n) A1Y(n 1)ApY(n p) B0W(n) B1W(n 1)BqW(n q) (7-13)其中Y(n)是l维向量， W(n)是l维白噪声， Ai,Bj为l l 阶参数矩阵。记A(p,z1) I A1z1Apz p1 1 qB(q,z1) B0 B1z1Bqz q (7-14)1 1 1 1H(p,q,z1) A(p,z1) 1B(q,z 1) 则(7-13)式可改写为A(p,z 1)Y(n) B(q,z

7、1)W(n) 或 Y(n) H(p,q,z 1)W(n)(7-15)向量过程 Y(n) 的功率谱密度函数矩阵为SY(ej ) H(p,q,e j )*SW(ej )H(p,q,e j )TA1(p,e j )B(q,e j )*SW(ej )A1(p,e j )B(q,e j )T A1(p,ej )B(q,ej )SW(ej )BT(q,e j )AT(p,e j ) (7-16) 其中 SW(ej )Rl l是白噪声的协方差矩阵。显然，如果我们获得了过程模型参数及 l 维白噪声的协方差矩阵 的估计， 也就获得了过程功率谱 的估计。前面讨论的标量过程的 AR 、ARMA 建模与谱估计可以推广

8、到多变量过程。7.3 向量 AR过程的建模 AR(p)过程的 Y W方程 对(7-8)式右乘 YT (n m)并取期望，得E Y(n) A1Y(n 1)Ap Y(n p) YT (n m) EW(n)YT(n m) (7-17)因此，有RY(m) A1RY(m 1)Ap RY(m p) RWYT (m)由因果性质知，当 m小于零时， Y(n m) 与W ( n)无关；同时，考虑到 W (n)的零均值特性，有EW(n)YT(n m) E W(n) E YT(n m) 0， m 0 而当m为零时，利用( 78)和( 719)有EW(n)YT(n) E W(n) W(n) A1Y(n 1)ApY(n

9、 p) TEW(n)WT (n) QW(7-18)(7-19)因此RWY(m)0Q,W ,mm 000, m 0(7-20)(7-21)RY(m) A1RY(m 1)Ap RY(m p) QW,m 00, m 0(7-22)展开RY(0) A1 RY (1)Ap RY(p) QWRY( 1) A1 RY (0)Ap RY(p 1) 0(7-23)RY( p) A1RY( p 1)Ap RY(0) 0对(7-23)式求转置，并考虑到相关性质RY(m) RYT( m)，则式 (7-22)可改写为RY(0) RYT (1) A1TRYT ( p) ApTQWRY(1) RY(0) A1TRYT(p

10、1)ApT 0(7-24)RY(p) RY(p 1)A1TRY (0) ATp 0(7-24)为向量 AR(p) 过程 Y W 方程，令RY(0),p(0)RY(1)RYT (1)RY(0)RYT ( p)RYT(p 1)(7-25)则有RY(p)RY (p 1)RY(0)l(p 1) l(p 1)IA1TpQW0ApTl( p 1) l0l(p 1) lppp(7-26)RY,p(0)(7-27)1互相关矩阵1)RY, p(0) 0；当 Y( n)中不存在零分量时， RY,p(0) 正定。RY,p(0) 的性质：RY,p(0) 是非负定的，即证明：作随机过程z(k)aiT Y(k i)i0其

11、中， aiai,1,ai,2 , ,ai,npE z2(k) E aiT Y(k i)2i0i 0,1, , p )是实数向量， ai,j ( j 1,2, , n )不全为零。则 pEaiT Y(k i)YT(k j)aji, j 0 i j ppaiT E Y(k i)YT(k j) ajaiT RY(i j)aji, j 0RY(0)RY( 1)RY ( p)a0RY (1)RY (0)RY ( p 1)a1RY (p)RY(p 1)RY(0)api,j 0RY(0)RY(1)RYT(1)RY(0)RYT ( p)a0RYT (p 1) a1其中,p(0) 0当 Y(n) 不存在零分量时

12、，z(k)必定是RY (p)RY(p 1)RY(0)ap注：以下讨论的随机向量过程2)递推性质。若记根据矩阵中的各个子块3)若已知互相关矩阵其中a0T a1T非零的，apTE z2(k)TRY,p(0)0， RY ,p (0)正定。Y(n) 均假设不含零分量。】RY(1)RY(2)RY(p)p1RY(p),pRY,p(0)RY(0)TpRY,p 1(0)RY(i) 的排列规律，显然。】RY(0)， RY ( p) ，则矩阵方程 (7-27)式有唯一解： ?pRY,1p 1(0) P， Q?e,p RY(0)Tp ?pp证明：根据 (7-31)式， (7-23)式中的参数向量A1 A2App 可

13、改写为Ilppp是， (7-29)式可改写为RY(0)RY,p T1p(0)pQWOlp l(7-28)(7-29)(7-30)(7-31)(7-32)(7-33)其中， Olp l为lp l 阶零矩阵。 (7-33)式等价于矩阵方程组RY,p 1(0) pP lp l ，RY (0)Tp p Qe, p显然，由矩阵 RY,p 1(0)的正定性质， (7-34)式有唯一解 (7-30) 式。2向量过程 AR 建模的互相关矩阵法 作 p 阶线性预测Y(n) A1Y(n 1)Ap Y(n p)则预测误差 ep(n) 为ep(n) Y(n) Y(n) Y(n) A1Y(n 1)ApY(n p)I A

14、1 ApYT(n) YT(n 1) YT(n p)T 记hpT(n) YT (n) YT(n 1)YT(n p)则(7-36)式可改写为ep(n)Tp hp (n) 或 eTp(n) hpT(n) p所谓 p阶最佳线性预测就是使预测误差 ep(n) 的互相关矩阵Qe,pEep(n) Tep(n) E Tphp n()hTp n(p)TpE hp(n)hTp n()p最小的 p 阶线性预测。由于(7-34)(7-35)(7-36)(7-37)(7-38)(7-39)E hp(n)hTp(n)Y(n)TYT (n)EY(n p)RY(0)RY(1)RY(p)RY( 1)RY(0)RY (p 1)将

15、 (7-40) 式代入 (7-39)式，得Qe,pTpE hp(n)hpT(n)考虑到互相关矩阵 RY,p(0)的性质 2)，(7-41)式可以进 IlpT RY(P0)TTYT (n 1) YT(n p)RY ( p)RY ( p 1)RY(0)RY,p(0)TTp RY , p (0) p步改写为TpRY,p 1(0)定理： l维随机向量过程的 p阶最佳前向线性预测参数矩阵p以及预测误差 ep(n) 的互相关矩阵足YW方程 (7-24)或(7-27)式，也就是RY,p(0) ? p ?pQe, pIlp其中?Tp I?pT ， ?pQ?e,p0 0,p证明：由 (7-42) 式，得Qe,p

16、 RY (0)pT PTppT RY,p 1(0) Ilp pT RY,p 1(0) p pT ppTPTpp pRY(0)RY(0)Tp RY1,p 1(0) PpT P RY,p 1(0) pTp pTp RY1,p 1(0)(7-40)(7-41)(7-42)Qe,p满(7-43)(7-44)(7-45)Cp RY(0)TpRY,p 1(0) P (7-46)则(7-45)式可以改写为Qe,p CppT RY,p 1(0) p PTp RY1,p1 (0) RY,p 1 (0) pPCppT Tp RY1,p1(0) RY,p 1 (0) p PCppT RY, p1(0) pRY,p1

17、(0)RY,p1(0)p P利用 RY,p 1(0) 的对称性，上式可以进一步改写为Qe,p Cp RY,p 1(0) pPT RY,p 1(0)RY,p 1(0) p P (7-47)显然， (7-47)中第一项 Cp 与参数 p无关，第二项是一个非负得二次型，当第二项为零时，Qe,p 达最小值。因此有RY,p 1(0) ?pP lp l (7-48)(7-49)Y(n) 的 p 1阶最佳线性RY (0)Tp ?p Q?e,p将 (7-48) 和 (7-49) 式合并，即得 (7-43) 式。】结论： p 阶最佳线性预测参数满足 p 阶多变量 YK方程 在上述推导过程中， p 是任意的正整数

18、。因此可以推断，随机向量过程预测参数满足 p 1 阶多变量0p1(7-50)P 1,invRYT (p 1)RYT(p)RYT (p 1)p,inv(7-51)Q?e,p 1 RY(0)p 1 p 1根据 (7-25)式，令RYT (1)则Y(n) 的互相关矩阵又可以表示为RY ,p (0)RY,p 1(0)p,invp,invRY(0)(7-52)I(A22 A12 A11 A12 )A12 A11A11 A12 ( A22 A12 A11 A12) T 1 1(A22 A1T2 A111A12) 10I利用分块矩阵求逆公式可得A1 1 A 1 2 ITA12 A220A111A11 A12

19、 A11I01 T 1 1 T 1A11 A12 (A22 A12 A11 A12) A12 A11T 1 1 T 1(A22 A1T2A111A12) 1 A1T2A111G11G12G221RY,p 1(0)1RY,p 1(0) p,invT1 RY (0)Tp,inv RY,p 1(0) p,inv(7-53)RY1, p(0)RY,p 1(0)1p,invp,invG11 G12G22G12G12G22p,invRY(0)G22 G12G22(7-54)并且，由 (7-53) 和 (7-48) 式，得G12p p,inv RY,p 1(0) p p,inv ?p(7-55)是，由 (7

20、-50) 式得1G11 G12G22G12p 1RY,p p 1 TG22 G12G12G22G22RY(pp 1)G11p G12 G22 G12pG12 G22 RY (p 1)将(7-48)和(7-55)式代入上式，并考虑到?p?p1G22 G12 p G22 RY (p 1)(7-51)式，得G12 G22 p,inv pG22 p,inv ?p G22RY (p 1) ?p pG12G22RY (p 1)G12 G22p 1,inv ?pG22p 1,inv ?p(7-56)不失一般性，令p1?pOl lp1(7-57)引入中间矩阵p 1 p 1,inv ? p(7-58)则由 (7

21、-56)式得【特例，当 l 1时， 又由 (7-28)和 (7-57)式，得G12 G22 p 1G12G22 p 1 I1G22 Df1(p)，G22 p 1 (p 1) 】p1Tp 1 p 1Tp RYT(p 1)因此，由 (7-50) 式中的第二子式，得Q?e,p 1 RY(0)Tp 1?p 1Q?e,pTp 1利用 (7-28)以及 (7-59)式， (7-60)式中的第二项可改写为 p 1Tp RYT (p 1) GI12p 1 I l lTp1G22 p 1(7-59)RY(0)Tp ?pTp 1p1G22p1p1TpG12 RYT(p 1) G22 p 1(7-60)(7-61)

22、将(5-55)式代入 (5 61)式，并考虑到 (7-51) 和(5-58)式，有Tp 1 p 1?pT p,inv RYT(p 1) G22 p 1RYT (p 1p,inv?Tp p 1,invG22 p 1(7-62)(7-63)(7-64)第二步：计算中间变量G11 RY,p 1(0)， G12 G11 p,inv ，G22 RY (0)p,inv G121 N p 1 R?Y(p 1) 1Y(n)YT (n p 1)NRYT(p 1) ，n0P 1,invp,invp 1 p 1,inv ?pQ?e,p 1 Q?e,pTp 1G22p1p 1G22 p 1将 (7-62) 式代入 (

23、7-60)式，得Q?e,p 1 Q?e,pp 1G22 p 1 Q?e,p实际应用中，互相关函数 RY ( p)的估计由下式计算1 N pR?Y(p) 1Y(n)YT (n p)， p 0,1,Nn0 综合上述，有如下依阶次递推的算法：第一步：初始化1 N pR?Y(p) 1Y(n)YT(n p)， p 0,1N n 0R?Y,0(0) R?Y(0) ， 1 R?Y (1)， 1,inv R?YT (1)?1R?Y,10 1 ， Q?e,1 R?Y(0)1T ?1， p 1第三步：判阶，即递推终止条件判断如果 | Tp 1G22 p 1|/ | R?Y (0) | ，结束递推，转到第四步；否则

24、，结束递推。第四步：估计参数p 1G1I2第五步：逆矩阵递推R?Y,1p(0)G11G12G22G1T2G12G22G22 G1T2G22p p 1转到第二步。 向量过程 AR 建模的 RLS 方法问题：根据随机向量过程的观测数据 Y(n) ，建立 AR模型，使前向线性预测误差ep(n) Y(n) A1Y(n 1)ApY(n p)(7-67)的方差M e, p E ep(n) E ep(n) ep(n) E ep(n)min(7-68)1)当E ep(n) 0时， Me,p E ep (n) epT (n) Qe,p2)当 ep(n) ep,1(n),ep,2 (n), ,ep,l(n) 中，

25、各分量互不相关时，即E ep,i (n) ep,j (n) E ep,i(n) E ep,j (n) 0， i j(7-69)时，有Qe,p diag E e2p,1(n) ,E e2p ,2 ( n) , ,E e2p,l (n)(7-70)使 Qe,p 最小，等价于使 ep (n)的每个分量的方差最小，即Jp,i E e2p,i (n) min ， i 1,2, ,l(7-71)等价于使lJ pep2,i(n)eTp (n)ep(n) min(7-72)i1当观测数据为有限序列 Y (0), Y(1), ,Y(N 1) ， N l (p 1) ，并且不加数据窗时，指标 函数可以改写为N1J

26、p(p, N 1)eTp(n)ep(n)npN 1 le2p,i (n)n pi 1(7-73)记p Ap,Ap 1, ,A1T(7-74)hp(n) YT (n p) YT(n p 1)YT (n 1) T(7-75)则有eTp (n) hpT (n) p YT(n)(7-76)又记ppT, Il T ， dp(n) hpT(n)YT(n)T(7-77)则(7-76)式可改写为eTp (n) dpT(n) p(7-78)Hp(p,N 1) hp(p) hp(p 1)hp(N 1) T(10)Dp(p,N 1) dp(p) dp(p 1)dp(N 1) T(11)ep(p,N 1) ep(p)

27、 ep(p 1)ep(N 1) T(12)Y(p,N 1) Y(p) Y(p 1) Y(N 1) T(13)则(14)ep(p,N 1) H p(p,N 1) p Y(p,N 1)Dp(p,N 1) p(15)其中 pRlpl， p R(p1)l l，hp(n) Rlp1，Hp(p,N 1) R(N p)lp，dp(n)l( p 1) 1 R，9Dp(p,N 1) R(N p)l(p1)，ep(p,N 1) R(N p)l，Y(p,N 1) R(N p)l。1) 正则方程求解 由(3)式得22Jp(p,N 1) ep(p,n 1)F Y(p,N 1) Hp(p,N 1) p F l2Yi (p

28、,N 1) H p(p,N 1) p,i 2 (16) i0其中Yi(p,N 1)是Y(p,N 1)的第 i列， p,i是 p的第 i列，即Y(p,N 1) Y1(p,N 1) Y2(p,N 1)Yl (p,N 1)p p,1 p,2 p,l由于 p中各列参数 p,i 是相互独立的，由指标 Jp(p,N 1) 取极小的必要条件，得HpT(p,N 1)Hp(p,N 1) p,iHpT(p,N 1)Yi(p,N 1)， i 1,2, l (17)将上式中 l 个矩阵方程合并，得HTp(p,N 1) Hp(p,N 1) pHpT(p,N 1)Y(p,N 1) (18)( 18 )式为最小二乘正则方程

29、，最小二乘解为?pH pT(p,N 1) Hp(p,N 1) 1HpT(p,N 1)Y(p,N 1) (19)利用最小二乘递推算法 (RLS) 可以避免矩阵求逆。2) 正交变换求解 根据 (16) 式，有Jp(p,N 1) ep(p,n 1)F Y(p,N 1) Hp(p,N 1) p FTrace Y(p, N 1) Hp(p,N 1) pT Y(p,N 1) Hp(p,n 1) p T r a ceTp DpT(p,N 1)Dp(p,N 1) p (20) 对于信息矩阵 Dp(p,N 1) ，存在正交变换矩阵 T ，使其上三角化，即Rp(p,N 1)T Dp (p,N 1) p 0 其中，

30、 Rp(p,N 1)是(p 1)l (p 1)l 阶上三角阵。于是有Jp(p,N 1) Tr aceTpDTp(p,N 1) TTT Dp(p,N 1) pTrace Tp RpT (p,N 1) Rp(p,N 1) pRp(p,N 1) p F (21)对上三角阵 Rp (p,N 1) 进行块划分Rp(p,N 1)其中， Rp' (p,N 1) Rpl pl ，G1(p,N 1) 到F-范数的性质， (21) 式可改写为R'p(p,N 1) G1(p,N 1)0 G2(p,N 1)Rpl l ，G2(p,N 1) Rl l 。将(22)代入(21)式，并考虑(22)1) p

31、G1(p,N 1) F G2(p,N 1)J p(p,N 1) Rp(p,N显然， (23)式中的第二项与参数p 无关；当第一项为零时，指标函数 Jp(p,N 1) 到达最小值。解为(23)因此，?pR'p1(p,N 1) G1(p,N 1)J?p ( p, N 1) G2(p,N 1) 正交变换求解也可以利用递推算法。'1(24)107.4 多元 ARMA 过程的建模1)因果性和可逆性A(p,z 1) Y(n) B(q,z 1) W(n) (25) 稳定性(因果性)判据若对所有 | z| 1， detA(p,z 1) 0，则 (25)式具有唯一的平稳解Y(n) H(z 1)W

32、(n)HiW(n i) (26)i0其中，矩阵 Hi 由下式唯一确定H (z 1) Hi z i A 1(p,z 1) B(q,z 1)(27)i0用 A 1(p,z 1 )左乘(25)式，并比较等式两边 z i的系数，得H 0 IiH i BiAjHi j，1 i q (28)j1pHiAj Hi j， i qj1最小相位(可逆性)判据若对所有 | z| 1， det B(q, z 1) 0，且 Y(n)是(25)式的平稳解，则W(n) F(z 1)Y(n)Fi Y(n i) (29)i0其中矩阵 Fi 由下式唯一确定F(z 1)Fi z i B 1(q,z 1)A(p,z 1)(30)i0

33、用 B 1(q,z 1)左乘 (30)式，并比较等式两边 z i 的系数，得F0 Imin( i ,q)Fi AiBj Fi j，1 i p (31)j1qFiBj Fi j， i pj12)向量 ARMA 过程模型估计的 LS两步法根据随机向量过程的观测数据 Y(0), Y(1), 建立ARMA 模型是一非线性问题。然而与单变量 AMAR 过程相同的是，如果向量 ARMA 过程是最小相位 (可逆 )的，则可表示为等价的无限阶的向量 AR 过程。因此，我们可以直接将单变量过程 ARMA 建模的两步法推广到向量过程， 具体做法如下所述。 【注：在向量情况下， F(z 1 )可能是有限阶次的，这是因为1 1 1 1 1 1 1F(z 1) B 1(q,z 1)A(p,z 1)1 B (z 1)A(p,z 1)detB(q,z 1)其中， B (z 1)是 B(q,z 1 )的伴随矩阵， detB(q,z 1) 为 B(q,z 1)的行列式。由于有限阶次矩阵多项 式的伴随矩阵也是有限阶次的。 因此