1、1第六章 单变量随机过程的建模与谱估计6.1 概述6.1.1 引言一个平稳过程的主要统计特性可用它的协方差函数或等价的谱函数来刻画。但是协方差函数和谱函数都是用无限个数值，即一条曲线来描述，使用不方便，尤其是在预测、滤波、控制等应用场合。根据随机信号理论知，若平稳随机信号 满足表示定理条件，它可以被看作是某一线性定常系统ny对白噪声输入的响应，即)(1zH(6-1)nwzH)(1其中 是方差为 的白噪声， 是后移位算子 的有理多项式，并且 的功率谱密度函数nw2w)(1z ny可以表示为(6-2)2)(wjjyeS因此，可以用系统传递函数 的有限个参数来刻画该过程的主要统计特性。)(1z另外，
2、经典的随机信号谱估计方法是确定信号谱估计方法的推广，所获得的谱估计量是非一致的。各种改进算法，如平滑、加窗、分段平均等，但是它们也都在曲线平滑程度与峰值分辨率之间存在矛盾。用模型估计方法的基本思路是，先估计过程模型参数，再由参数估计谱，以便得到更好的估计量。6.1.2 平稳过程的线性模型表达由线性环节 的三种形式)(1zH(6-3)(1z ),(,)(1111 zpAqBzazabbzzzpqp可以导出平稳过程的三种标准的线性模型表达：1)p 阶自回归过程阿 AR(p)(6-4)nwyA),(1或(6-5)pnnaay12)q 阶滑动平均过程 MA(q)(6-6)nzqB),(1或(6-7)q
3、nnwbbwy13)p、q 阶自回归滑动平均过程 ARMA(p, q)(6-8)zqzpA),(),(1或(6-9)qnnnpnn bbyaya 111)若 是渐近稳定的，则 AR(p)过程等价于 MA()过程，即)(1zA(6-10) 001),( iinniinn whzhwzpA2)若 是 最 小 相 位 的 ， 则 MA(q)过程等 价 于 AR( )过 程 ， 即)(1zB2(6-11) 001),( iinniinn yfzfyzqBw3)若 是渐近稳定且最小相位的，则 ARMA(p, q)过程等价于 AR()过程，即,)(1pAzH(6-12) 0101 ),(),( iinni
4、inn yhzpAzfyzq4)若 是渐近稳定且最小相位的，则 ARMA(p, q)过程等价于 MA()过程，即,)(1pBz(6-13) 0101 ),(),( iinniinn wfzqBzhwzAy6.2 平 稳 过 程 的 AR(p)模 型 估 计6.2.1 AR(p)过 程 的 标 准 方 程对式(6-5)左乘 ，并取数学期望，得mny(6-14)1()npnnmayyw即(6-15)()yyyyRaR因为AR(p) 过程可以用 MA() 表示，即 0iinnh因此0()ywnminmiREw(6-16)mwiwi ni hihh 2020 )(0,2w将式(6-16)代入式(6-1
5、5)，得自相关函数的通解, (6-17)0)()1()(1 pRamRaypyy 边界条件(6-18) 0)()1()( 1011 2ypyy py wRaRap (6-18)式就是 AR(p)过程的标准方程，又称为 YuleWalker 方程（简称：方程） 。记(6-19)0()1() 10)()()()( yyy yyyy RpRpR 3, (6-20)paA1)( 0)(2wpE则式(6-18)可以改写为矩阵形式(6-21)()(ARy定义置换矩阵 P01p 是单位正交矩阵，即 。记PIPT, (6-22)1)(1apApinv 20)(winvpE则有， ，Ppinv)( Ppinv)
6、(TyyPpR)()(因此，对（621）左乘 、右乘 ，得T(6-23)EARinvivy6.2.2 AR 过程的前向线性预测所谓 阶前向线性预测就是用过程 的 时刻以前的 个观测值的线性组合 来pnYp,1pnny估计时刻的值 ，使前向预测误差ny(6-24)nfype)(,的均方差(6-25)(2,Dff最小。阶前向线性预测可以表示为(6-26) piinnyay1)(阶前向线性预测误差为(6-27)nnfn zAe),)(1, 其中ppazpazpA (111) 阶前向线性预测误差的均方差为(6-28),()(21nf yzAD由最小值条件 ，得0pamf4, (6-29)0)(, mn
7、fnypep,21即, 0)(11 mnppn yaay p,21因此，最优的阶前向线性预测方程(6-26)中的系数满足, (6-30)()( RRpmRyyy 并且，阶前向线性预测误差为(6-31)nfnzAe),()1,其中ppzapazpA )(1,( 1z 阶预测误差的方差为(6-32)()(2,eDfnfpiinfnye0. )()(将式(6-29)代入式(6-32)，得,nfnf yp(6-33)()1()0(1 pRaRaypy 显然，联立式(6-33)和式(6-30)就是方程(6-18)，只是 AR-()过程模型输入的方差 变为预2w测误差的方差 。)(pDf结论：）AR-过程
8、的最小均方意义下的阶前向线性预测参数满足阶方程。）如果 AR-过程 的真实阶次就是，则 ， 。此时，ny2w)(pDf)(paii阶前向线性预测误差 成为白噪声。因此，前向线性预测误差滤波器 又称为白化滤波)(,efn ,1zA器。5.2.3 AR-过程的后向(线性)预测所谓阶后向线性预测就是，用过程 的 n-p 时刻以后的个值 来估ny11pnnyy， 计 n-p 时刻的值 ，使后向线性预测误差pny(6-34)pnbnpe)(,的均方差(6-35)(2,eDb最小。阶后向线性预测可以表示为(6-36) pi ipnnyy1)(阶后向线性预测误差为(6-37)ninvpnbn yzBe),)
9、( 1, 其中 ppiv zpbbzB 1111 ()(zinvp),(1zpA),(1zpA nwny)(,efn5阶后向线性预测误差的均方差为(6-38),()(21ninvbyzpBD由最小值条件 ，得0)(pmb, (6-39)0,mpnnye,即, 0)()(11 mpnpp ybby p,21因此，最优的阶后向线性预测方程(6-36)中的系数满足, (6-40)( RRyyy 并且，阶后向线性预测误差为(6-41)ninvbnzpBe),()(1,其中 pppinv zpbzB 1111 )(),( )(zinv阶后向线性预测误差的均方差为(6-42)()(2,epDbnbpi i
10、pnbnye0.将式(6-39)代入式(6-42)，得,pnbnby(6-43)()1()0(1 RRypy 显然，联立式(6-43)和式(6-40)也是方程(6-18)，只是 AR-(p)过程模型输入的方差 变为2w阶后向线性预测误差的均方差 。Db由方程解的唯一性知，阶前向线性预测方程(6-26)和阶后向线性预测方程(6-36)中的系数相等，即(6-44)()(,21,piafbii因此，阶前、后向线性预测滤波器的关系为或 (6-45),(11zAzpBinvinv ),(),(11zpAzB最优的阶后向线性预测误差方程(6-41)又可改写为(6-46)nibype),()(,结论：AR-
11、过程的最小均方意义下的阶后向线性预测参数满足阶方程。6.2.4 AR-模型参数估计的 Levinson-Durbin 算法（简称-算法）算法是已知相关函数 ，估计模型参数。根据前二节的讨论知，AR-过程的阶前)(mRy向线性预测模型参数 满足阶方程，即)(pA(6-47)()(pEAy且(6-48)Rinvinvy阶前向线性预测模型参数 满足)1(p(6-49)1(pEAy6由于 )0()1()()1( )()()0()( 11)( yyyy yyytyy RpRp ppR 记(6-50) )()2()( phyyyinv 则(6-51)0()1()(yyy RphR )()1(0Rhyinv
12、invy不失一般性，记(6-52)()(pAA将式(6-52)代入(6-49)并考虑到式(6-51)、(6-47)，则有)1()(pRy 0)(pRhyy )1(0)()1(0pARhyinvinvy (6-53) )()(1( DAApEfyinv 令(6-54)1php则，式(6-53)可以改写为(6-55)1()(ApRy )1()()( AREyinv由式(6-55)的后行得(6-56) )1(0)1()ppy 由式(6-55)的第一行得(6-57)()( AhDinvff 讨论：）若 ，由式(6-56)知 。此时，为过程模型的真实阶次， 为0)1(p 0)1p )(pA过程模型的真实
13、参数。）若 ，在式(6-56)二端同除(6-58)()(pDf7有(6-59)()(0)1()1( pEDpARp invfy 比较式(6-56)与式(6-22)，有(6-60)1()(ppinv即(6-61)(AAinv将式(6-61)代入式(6-52)，得(6-62)(0)1(0)1( ppinv根据(6-23)及(6-50),显然有(6-63)()( AhAhinvinv因此,根据(6-61)和（6-63） ，并考虑到式(6-54)，得)1)1)1( ppAph invivivinv (6-64)1()()( p将上式代入式(6-57)，并考虑到式(6-58)，得(6-65)( 2DDf
14、ff 综合上述，可得如下算法：(6-66)1()()1( 00)(1)(2pDpAAphff invf初始条件：， (6-67)inv0yfR结束条件: 根据以下性质）反射系数 的绝对值小于，即 ；)1(p1)(p）前向线性预测误差的均方差是阶次的单调下降函数，即(6-68)(DDfff 可得算法递推结束条件：(6-69)1(pff或(6-70)0(ff其中， 为充分小的正整数。算法的特点：8）按阶次递推，便于阶次的选择；）计算量小；）估计出的滤波器是最小相位的（后面证明） 。6.2.5 基于观测数据的算法算法的前提是已知过程的自相关函数 ，但实际问题通常是用过程的观测数据来估计)(mRy过程
15、的模型参数与阶次。设平稳过程 的个观测数据为 ，希望将过程 表示为一 AR-过程：nY110,L ny(6-71)pi nn peyay)()(使该过程输入激励 平方和)(pen(6-72)NMnp)()(2,最小。记(6-73)Tpnnnaapyyh 21)(则式(6-71)可以改写为(6-74)()(eynTn根据 区间上的数据，构造向量与矩阵NPM,(6-75)TNMNMpNpMMTNNMeepEyyHyY)()()()()(1, 112, 1, 则式(6-72)可以改写为)(,pTNM)()(, pHYHYNMTN (6-76)2, pM式(6-76)表明，AR-模型估计问题可以化为参
16、数的最小二乘估计问题。对于给定的阶次 ，式(6-76)的最小二乘解满足正则方程p(6-77)0)()(, AYNTN其中 。最小二乘的指标函数为pA)(,1)(6-78)()(, pHMTM 将式(6-77)和(6-78)合并，得(6-79) 0)( , AYpHY NMNTNM 令)(),(11phjirjMTi9, (6-80)NMnjniypji0,则, 式(6-79)可改写为(6-81) 0)()(1),()1,()0,( ,1)()()( , parprrrrr NMNMNMNNNN 1)自相关法若数据前后加窗，即当 或 时，取 ，则利用全部数据并按(6-71)可以确定出nL0ny之
17、间的全部预测误差 。因此，式（6-80）和（6-81）中取 ，10pLn )(en 。N由于(6-82),(1,0jirpL10pLnjniyipLisjisy110Lsjisy显然， 与 p 无关，并且1,0jip(6-83),(1,ijrpL),(1,0jrpL成立。若令(6-84)(,)( 1,01,0 pDjijiRLp 则式（6-81）可以化为(6-85)0)()0()1()(01paRpR可以看出，式（6-85）与 AR-()过程的 Yule-Walker 方程完全一样，只是用估计值 代替真值)(mR。因此，可以将相关函数 的估计与 L-D 算法结合起来，形成如下递推算法：)(mR
18、m(6-86)1()(1( (00)(1)()(210pDpAAphyLinvfTpii初始条件：p=0， ， (6-87)invA102LiyR终止条件：(6-88)0(1Dp10其中， 为充分小的正整数。算法特点：1）从理论上讲，此算法具有 L-D 算法的优点；实际计算过程中，可增加模型稳定性监测，以确保反射系数的模小于 1。2）加窗是一种人为的假设，即按式（6-71）确定 时，利用了条件)(pen， (6-89)011yyp 011pLLy这将对参数估计结果产生影响。仿真实验表明，只有当样本数据长度 L 充分大时，数据窗的影响可以忽略，该算法能给出比较精确的结果。对于较短或中等长度的样本
19、数据，这一方法往往不能给出可靠的结果；也可能出现数值不稳定的情况（有限字长以及 的有偏估计可使（6-85）式中自)(mR相关矩阵实际上不具有整定性） ，特别是当被估计的模型接近于稳定域边界，即传递函数的极点靠近 z平面单位园时。2）协方差法若不加数据窗，则根据观测数据和(6-71)式无法确定出区间 和 上的前向线性1,0p,pL预测误差 。因此，(6-80)式中取 ， ，即)(pen pMLN(6-90)11,)(LpnjniLpyjir此时 与 有关，并且),(1,jirLp(6-91),(),(1,1, jijiLpLp由 ，按(6-81)式构成的自相关矩阵不具有 Toeplizs 性质，
20、也就不能采用 L-D 算法求解(6-,1,jip81)式。3）前、后加窗所谓前加窗就是，假设 ，(6-80)式中取 ， ；所谓后011yyp 0M1LN加窗就是，假设 ，(6-80)式中取 , 。1LLy p两种情况下，由 或 ，按(6-81)式构成的自相关矩阵不具有 Toeplizs 性质，),(,0jir),(1,jirp同样不能不能采用 L-D 算法求解(6-81)式。6.2.6 Burg 算法L-D 算法的关键是计算反射系数 ，Burg 算法的基本思路是，对数据前后加窗，以保留)(L-D 算法的优点；同时，通过使前、后向预测误差能量和(6-92) 12,2,1, )1()(LpnbnfLp pee最小来减少数据窗的影响，并保留预测误差的递推关系(6-93)()(1()(),1, , eebnfnbnff将上式代入指标（6-92）式，并由最小化条件(6-94)0)(,pL得(6-95)12,12, )()()(LpnbnffnpeBurg 算法:
