第七章现代谱估计电子教案

1、第七章 现代谱估计经典谱估计以傅立叶变换为基础，具有计算效率高的优点，但是由于将未观测数据认为0和数据加窗，具有频率分辨率低、旁瓣泄漏等严重的缺陷。为此，近几年来，在提高功率谱估计的分辨率方面提出了很多新的方法。以1967年burg提出的最大熵谱分析法为代表的现代谱估计法，以参数模型为基础，不认为在观察到的n个数据以外的数据全为零。因此克服了经典谱估计法的缺点，提高了谱估计的分辨率。后来发现线性预测自回归模型法(简称ar模型法)与burg的最大熵谱分析法是等价的，它们都可归结为通过yule-walker方程求解自回归模型的系数问题。目前常用的求自回归模型系数的算法有三种：为levinson递推

2、算法；为burg递推算法；为正反向线性预测最小二乘算法。除了最大熵谱分析法(包括线性预测ar模型法)外近年来出现了许多适用于不同情况的提高谱估计分辨率的新方法，如模型法中还有滑动平均(ma)模型法与自回归滑动平均(arma)模型法，另外还有pisarenko谐波分解法，prony提取极点法，prony谱线分解法以及capon最大似然法等等。本章主要讨论最大熵谱分析法(包括线性预测ar模型法)，它是目前用得最多的一种高分辨率的谱估计方法。参数模型估计法就是根据已观察到的数据，选择一个正确的模型，认为x(n)是白噪声通过此模型产生的，这样就不必认为n个以外的数据全为零了。这就有可能得到比较好的估计

3、。这种方法分以下三个步骤进行。三个处理步骤为1 确定或选择一个合适的模型依赖于对所研究随机过程进行理论分析和实验研究；2 根据观测数据估计模型参数涉及各种算法的研究；3 由模型参数计算功率谱。参数模型谱估计法的关键问题是 ：模型选择问题（ar， ma ，arma）和参数确定方法（导致产生了各种算法）7.1 自回归模型谱估计 7.1.1 建立模型 在实际中我们所遇到的随机过程，常常总是可以用一个具有有理分式的传递函数的模型来很好地表示它，因此可以用一个线性差分方程作为产生随机序列x(n)的系统的模型： （7.1）这里表示白色噪声，下图所示为离散随机信号x(n)的有理传输函数模型，输入为零均值、方

4、差为的白噪声序列。w(n)x(n)将上式(7.1)变换到z域则有该模型的传递函数为：(7.2)其中 (7.3)式中ak为自回归系数，称为ar系数；bk为滑动平均系数，称为ma系数。当输入的白噪声的功率谱密度为时，输出的功率谱密度为(7.4)将代入上式得(7.5)如果能研究各ak及bl就可求得，于是，求功率谱的实质变为确定系统参数的问题。 7.1.2 三种模型ar模型设，并不会影响式(7.1)与式(7.2)的一般性。如果除b0=1外的所有bl均为零。则式(7.1)成为：(7.6)式(7.6)的形式被称为p阶自回归模型简称ar(autoregressive)模型。将式(7.6)进行z变换，可得ar

5、模型的传递函数为 (7.7)自回归模型的h(z)只有极点，没有除原点以外的零点，如图7.1所示。因此又称为全极点模型。当我们采用自回归模型时，式(7.5)成为(7.8)此时，只要我们能求得所有ak参量，就可求得。图7.1 自回归模型ma模型如果式(7.1)中除a0=1的所有ak均为零，则式(7.1)成为(7.9)式(7.9)的形式称为q阶滑动平均模型，简称ma(moving average)模型。ma模型的传递函数为（b0=1）(7.10)ma模型的h(z)只有零点没有除原点以外的极点，因此又称为全零点模型。当我们用ma模型时，式(7.5)成为(7.11)此时，只要我们能求得与所有的bl参量，

6、就可求得。arma模型 当均不完全为零时的模型称为arma模型，即极点零点模型。式(7.1)和式(7.2)分别表示了arma模型的差分方程与传递函数。由以上的讨论可见，用模型法作功率谱估计，实际上要解决的是模型的参数估计问题。7.1.3 word分解定理wold分解定理为我们对模型的选择以及以上三种模型之间的关系提供了理论基础。 wold分解定理告诉我们：任何一个有限方差的平稳arma过程可以分为完全随机的部分和确定的部分。推论：任何有限方差的arma或ma平移过程可以用可以是无限阶的ar模型表达；同样，任何arma或ar模型可以用可能是无限阶的ma模型表示。因此，如果在这三个模型中选了一个与

7、信号不匹配的模型，利用高的阶数仍然可以得到好的逼近。由于对ar模型参数的估计，如下面将要看到的，得到的是线性方程。故ar模型比arma以及ma模型有在计算上的优点，arma或ma模型一般需要解一组非线性方程。同时，实际的物理系统往往是全极点系统。所以，ar模型得到了深入的研究和广泛的应用。研究有理分式传递函数的模型，主要研究ar模型。本节也只讨论ar模型。7.2 ar模型的yule-walker方程下面我们就集中讨论ar模型的谱估计法。已知自相关函数yule-walker方程求解： ar模型的阶数p，以及p个ar参数a(i)和激励源方差。7.2.1 yule-walker方程的推导 由上面的讨

8、论已经得到有关ar模型的如下几个方程：(7.12)（7.12a）为了得到必须求得参数a1, a2, a3, , ap及。为此，让我们来推导这些ar参数与之间的关系。 按定义将式(7.12)的关系代入上式，得 (7.13)按式(7.12)，x(n)只与相关而与无关，故式(7.13)中的第二项为代入式(7.13)得：或 （7.14）将m=1，p分别代入式(7.14)并写成矩阵形式，得再利用自相关函数的偶对称性，则有： （7.15）式（7.14）及式（7.15）称为ar模型的yulewalker方程。yulewalker方程表明：只要已知输出平稳随机信号的自相关函数，就能求出ar模型中的参数ak及，

9、并且需要的观测数据较少。n个样值x(0),x(1)x(n)ar模型谱估计流程图：自相关函数r(0),r(1).r(n)ar模型参数和a1,a2，ap激励源方差功率谱密度7.2.2 yule-walker方程的求解方法由以后的讨论可以看到，最大熵谱估计法与线性预测谱估计法都与ar谱估计法等价，它们都可归结为求解yule-walker方程中的各ar系统ak(k =1，2，p)的问题，但是直接从yule-walker方程式求解参数ak(k =1，2，n)需要作求逆矩阵的运算，当n大时，运算量很大（其运算量达到p的三次方），并且每当模型阶数增加一阶，矩阵增大一维，需要全部重新计算。levinson-d

10、urbin算法对yule-walker方程提供了一个高效率的解法。该算法是按照阶次进行递推，运算量为p的二次方。下面介绍levinsondurbin递推算法。此算法的关键就是要推导出由第k阶ar模型的参数计算第k+1阶ar模型ar(k+1)参数的迭代计算公式。首先以ar(0)和ar(1)模型参数作为初始条件，计算ar(2)模型参数，然后根据这些参数计算ar(3)模型参数，等等，一直到计算出ar(p)模型参数为止。即：依次求得。注意，附加的a的第一个下标是指ar模型的阶数，最后p阶的解即是所要求的解。那么，以一阶ar模型(求一阶参数a11及)开始，按式(7.15)一阶ar模型的yule-walk

11、er矩阵方程应为从这个矩阵方程可解a11与，分别为(7.16)(7.17)再从二阶ar模型的矩阵方程：解得a22，a21，分别为(7.18) (7.19)(7.20)以此类推得递推公式：(7.21)(7.22)(7.23)于是，当我们从式(7.16)与式(7.17)得到初始的a11与的数据以及各，就可按式(7.21)，(7.22)，(7.23)依次递推出各阶的akk，aki，。从式(7.23)有，一般来讲阶数预先是不知道的，当我们递推到第k阶，满足所允许的值，就可选阶数p=k。实际上，这里的就是误差功率（后面章节有证明）。小代表均方误差小。如果信号的正确模型是p阶的ar模型，则应有(7.24)

12、这说明已达最小均方误差值。 由式(7.23)可见，由于，对于任何k必有(7.25)以上成为反射系数。将上述估计的模型参数代入功率谱估计式中，即可计算功率谱估计值为：综上，利用l-d算法进行ar模型参数谱估计算法流程如下。给定n个观测数据xn(n),n=0,1,2,n-1，ar模型参数估计方法为：i) 由xn(n)估计出自相关函数值，m=0，1，p；ii) 利用levinson-durbin算法根据计算ar(p)模型参数的估计值。具体步骤为：令p=1，计算a11及；接着，令p=p+1，则p=2，计算app，api，；iii) 直到m=p或，满足终止规则，流程结束；iv) 计算估计值。matlab

13、里有专门实现l-d算法的函数可估计ar模型参数：a e=aryule(x,p),a为模型参数，e为噪声方差。例7.1 已知实数据序列的自相关为：。用levinson-durbin递推算法求ar模型的参量：。解：7.2.3 ar模型阶的确定可以证明式(7.25)正是的所有极点均在单位圆内的(即稳定性的)充分必要条件，同时它也是自相关矩阵为正定矩阵的充分必要条件。另外，由7.2.2分析可知，激励信号的均方误差随着阶数的增加而递减。所以，由于ar模型具有以上性质，故具有稳定性。关于ar模型，还需选择一个合适的阶。如果阶太低，功率谱平滑的太厉害，平滑后的谱分辨不出真实谱中的两个峰；阶太高，可以提高谱估

14、计的分辨率，但会出现许多虚假谱峰或谱的细节。图7.2为阶数太高和太低的谱估计结果。图中加粗黑色实线为真实谱，点划线为虚假谱。图7.2 阶数太高和太低所得谱估计结果图 因此，ar模型谱估计方法，既要估计模型参数，又要估计模型的阶，在这样复杂的情况下，如何评价各种谱估计的性能，目前尚无定论。下面讨论几种常见的ar模型阶的确定方法。一、确定ar模型的阶的方法一： 一般的观察方法，简单而直观u 不断增加阶数，观察预测误差功率，当它下降到最小时，对应的阶数选为模型的阶；但是，预测误差功率（即ar模型激励源的方差）是随着阶数的增加而单调下降的，很难确定降到什么时候是最小；另一方面，随阶数增加，模型参数数目

15、增加，谱估计的方差变大，出现了虚假谱峰，故一般不能依靠观察预测误差功率下降来确定模型的阶数。u 不断增加阶数，观察各阶模型预测误差序列的周期图，最接近平坦（白色谱）时对应于最佳的阶数。二、依据不同的误差准则来确定模型的阶数。准则1：fpe(最终预测误差）n为观测数据长度，为拟合残差方差，随阶增加而减小，而随着阶数增大而增大。fpe将有一个最小值。fpe的最小值对应的阶数为最后确定的阶。 准则2、akaike（aic)信息准则（适用于ar和ma过程）对于arma过程，则aic准则定义为 i为模型的阶，为模型误差，一般随着阶的增加而减小，而式中第二项随阶次增加而增加。aic试图解决减小模型误差和保

16、持较少的模型参数数目之间的矛盾。aic定义式有一个最小值，这个最小值对应的阶就是要选择的阶。此外，还有cat等准则 。通过实验发现：在将这些准则用于估计ar模型的阶，对于实际数据，所得到的谱估计结果常常无太大区别。对于短数据，以上准则都不理想。在实际应用中，应该参照实验结果对模型的阶加以适当调整。7.3 线性预测谱估计假设x(n)是一个n阶ar过程，现在时刻x(n)的值可以由过去n个时刻的取样值的加权来预测，加权系数为-ak，那么n阶线性预测器：可看作用序列x(n-n),x(n-n-1) , ,x(n-1)激励一个冲击响应为-ak的线性时不变系统的输出值。x(n-n),x(n-n-1) , ,

17、x(n-1)-ak 线性预测误差为(7.26)将上式进行z变换，得于是(7.27)由式(7.27)可见，是以x(n)作为输入信号，误差e(n)作为输出信号的滤波器的传递函数，该滤波器称为预测误差滤波器。将式(7.27)与式(7.12a)比较，当apk按最小均方误差准则求得时，有apk=ak ，故有同时，预测误差功率为，确定系数apk的一个原则是使预测误差功率最小。根据这一原则推导出的预测器系数-ak与x(n)的自相关序列rxx(m)之间的关系为（证明过程略）：且有：，其中将两个关系式写成矩阵展开式分别为：将(7.281)和(7.282)两个关系式合并为一个式子：将(7.283)写成矩阵展开形式

18、为：可以看出：n阶线性预测器的系数ak与ar模型中的ar系数相等，即；预测误差概率最小值pmin与ar模型中的输入噪声方差相等，即。所以，线性预测谱估计与ar谱估计是等效的。7.4 最大熵谱估计（mese， maximum entropy spectral estimation）最大熵谱估计法简称为mese(maximun entropy spectral estimation)法。 7.4.1 按最大熵外推自相关函数在前面一章中我们已经讨论到经典法用已知的有限个(n个)自相关函数序列的估计求功率谱估计时，是将此有限个估值以外的自相关序列的数据认为是零，因而得不到好的分辨率。j.p.burg于

19、1967年提出的mese法与此不同，它是基于将已知的有限长度的自相关序列以外的数据用外推法求得，而不是把它们当作是零。如果假设已知问题在于按什么原则外推。在保证自相关函数的toeplitz矩阵是正定的情况下有无穷多种外推法，burg认为外推的自相关函数应使时间序列表现出最大熵，因此把burg提出的这种方法称之为最大熵谱估计法。我们知道所谓熵是代表一种不确定性的度量，最大熵为最大不确定度，其平均信息量最大，事件越不容易发生，即它的时间序列最具随机性，而它的psd应是最平滑(最白色)。按shannon对熵的定义，当x的取值为离散时，熵h定义为(7.29)这里pi为出现状态i的概率。当x的取值为连续

20、时，熵被定义为(7.30)这里p(x)为概率密度函数。当我们处理用时间序列传递信息的问题时，概率密度应由联合概率密度函数代替。设为零均值，高斯分布的随机过程。一维高斯分布为(7.31)其中(7.32)n维高斯分布为(7.33)其中 (7.34)这里(7.35)下面我们先求一维高斯分布的信号的熵，然后推广到n维。为此将式(5.78)及式(5.79)代入式(5.77)得 因为代入上式得到一维高斯分布的熵为(7.36)同理可求得n维高斯分布信号的熵为(7.37)式中代表矩阵的行列式，要使熵h最大，就要求最大。 如果已知现欲求得。由于自相关函数的矩阵必是正定的，故矩阵的行列必大于零，即(7.38)为了

21、得到最大熵，要求最大，为此用对上式微分，使，求得使最大的，满足下列方程：(7.39)上式是的一次函数，由此式可解出。于是又可以此为已知，再用类似方法求得，以此类推。这样每步都按最大熵的原则外推后一个自相关序列的值，可以外推到任意多个而不必认为它们是零。这就是最大熵谱估计法的基本思想。7.4.2 mese与ar谱估计等效可以证明这种按最大熵外推自相关函数的结果与ar模型是等价的。所以，上式(7.39)实质为yulerwalker方程。为了证明这一点，将式(7.14)，即中的m分别用1, 2, , n+1代入，并利用自相关的偶对称性，写成下列方程组：(7.40)(7.41)如果我们从式(7.40)

22、的n个线性方程中解得的n个ar参数a1, a2, an值，代入式(7.41)并将其整理成行列式的形式，即可得(7.42)注意，从ar模型得到的式(7.42)与按最大熵外推得到的式(7.39)完全相同，这就证明了当x(n)为高斯分布时最大熵谱估计法与ar模型法是等价的。事实上，当x(n)是高斯分布带限(为其最高频率)时间序列，且其功率谱密度满足时，则可以从式(5.78)出发证明其熵正比于在必须与已知的自相关函数的n+1个值符合的约束条件下，即的约束条件下，可以利用变分法证明使熵最大的功率谱为其中与an由下列矩阵方程确定(7.43)式(7.43)就是yule-walker方程，因此这就直接说明了最

23、大熵谱估计与ar模型等价的。7.5 预测误差格型滤波器及伯格(burg)递推算法 用levinson递推算法求解yule-walker方程中ar系数虽然可以简化计算，但需要知道自相关序列为阶数）。实际上自相关序列只能从时间序列x(n)的限个数据得到它的估计值。当时间序列短时，的估计误差很大，这将对ar参数ak(k =1，2，p)的计算引入很大误差，导致功率谱估计出现谱线分裂与谱峰频率偏移等现象。j.b.burg在1967年的“最大熵谱分析”一文中提出最大熵谱估计法时只说明了如何从已知的n个外推n个以外的值从而可得到高分辨率的功率谱，但该文并未涉及如何从有限时间序列来得到这个n个的问题。随后，b

24、urg又于另一篇文章中提出一种直接由时间序列计算ar模型参数的方法，被人们称为burg算法，这种算法与预测误差滤波器有密切关系。它是在levinson关系式(7.22)的约束下，用使前向与后向预测误差能量之和为最小的方法来求得各ak(k =1，2，p)的值。 7.5.1 预测误差格型滤波器按线性预测理论，已知n个观测数据x(1),x(2),x(n-1)，利用p阶线性预测滤波器估计x(n)得 x (n)的估计值，可用x(n)的各过去值的加权之和表示，即(7.44)误差e(n)，在这里用ep(n)表示(p为阶数)(7.45)ep(n)称为线性预测器的前向误差，因为是由x(n)以前的各数据：x(n-

25、1)，x(n-2)，x(n-p)加权之和得到的。由levinson关系式(7.22)可得(7.46)又令，这里的kp称为部分相关系数(parcor)或反射系数。将这些关系代入式(7.45)得 (7.47)这里 所以(7.48)这里 (7.49)是由x(n-p)以后的各数据：加权之和得到的，故bp(n)称为后向预测误差。比较前向预测误差方程(7.45)与后向预测误差方程式(7.48)可见，它们具有相同的系数。 将式(7.46)代入式(7.48)，用证明式(7.47)：类似的方法，可以证明(7.50)按式(7.45)及式(7.48)，当p=0时有(7.51)式(7.47)与式(7.50)为前向及后

26、向预测误差的递推公式。由式(7.49)、(7.50)及(7.51)可得格型预测误差滤波器，如图5.10所示。图5.10 格型预测误差滤波器如果将式(7.45)进行z变换，得于是有(7.52)它是预测误差滤波器的传递函数。因此线性预测误差滤波器可以用格型结构的fir滤波器实现(如图5.10所示)。将式(7.45)与式(7.48)分别写成展开形式，有(7.45)(7.48)比较上二式可见，后向预测误差bp(n)的系数正九是前向预测误差ep(n)的系数在次序上的逆转。因此，如果我们将横向结构的预测误差滤波器的输入不是从第一级(最左端)，而是从末级(最右端)送入，如图5.11所示，则输出将不是前向误差

27、，而是后向误差，有图5.11 后向预测误差(7.49a) (7.49b) (7.49c)这里(7.50)如果我们令，于是，式(7.49c)成为(7.49d)而其中(7.51)将式(7.51)与式(7.44)比较可更清楚的看到，后向预测不是像前向预测那样，用它以前的p个值的线性组合来预测它，而是用它以后的p个值的线性组合来预测它(即式(7.51)右边k前的符号与式(7.44)右边k前的符号相反)。前向预测误差滤波器的传递函数为 (7.52)这里，它与he(z)是一对z变换对，因此，如果输入序列，则输出序列。 由式(7.49)可得后向预测误差滤波器的传递函数(用hb(n)表示)(7.53)这里，它

28、与hb(z)是一对z变换对。将式(7.53)与式(7.52)比较显然有(7.54)由式(7.54)可见，当z=z1是hb(z)的一个零点(或极点)时，则它也将是he(z-1)的一个零点(或极点)。因此1/z1将是he(z)的一个零点(或极点)。于是，如果he(z)是一个最小相移网络，其零极点全部在单位圆内，则hb(z)将是一个最大相移网络，其零极点将全部在单位圆外。当时，he(z)将是稳定的和最小相移的，而hb(z)将是稳定的和最大相移的。全零点的格型预测误差滤波器是itakura和saita于1971年首先提出的。burg虽然在1967年就提出了有关的方法，但他并没有具体地把格型网络的性质和

29、他的最大熵技术等同起来，现在已经知道itakura和burg的方法是更一般的格型概念的特例。格型结构的fir滤波器与横向结构的比较有很多优点。在数字滤波器中有二个重要问题：有限字长效应的影响和参数值扰动对滤波器性能的敏感性。在这二点上格型结构的滤波器都优于横向结构的滤波器。由图5.10可见，一个p阶的格型预测误差滤波器是由p级组成，而前面各级的输出，正好依次是各低于p阶的预测误差滤波器的输出。当滤波器各阶输出信号都满足最小均方误差时，并在输入信号是平稳的情况下，可以证明，格型预测误差滤波器前后各级的输出误差之间的正交的。这种正交性导致前后各级之间的无耦性能(decoupleing proper

30、ty)，这种性能十分有用，它使全局最优可以用各级局部最优来实现。当用作自适应滤波器时，各级可选择不同的自适应步长，使其收敛速度提高。7.5.1 burg递推算法kp的确定 在实际应用中，根据信号的有限个取样值估计ar模型参数的方法，通常有自相关法、协方差法和burg递推法。自相关法和协方差法都是直接估计ar参数，而burg法是先估计反射系数，然后利用levinson-durbin递推算法由反射系数求得ar参数。burg递推算法的优点是不需要估计自相关函数，可以直接从已知的x(n)序列求得参数kp0另外，这种算法可保证满足稳定性的充要条件：。算法准则是前向均方误差和后向均方误差之和最小。 如果k

31、p按前向均方误差最小的准则确定并用kep表示，则按式(7.47)令即所以(7.55) 如果kp按后向均方误差最小的准则确定并用kbp表示，则按式(7.50)即所以(7.56)burg算法是以前向均方误差与后向均方误差之和最小为准则求得kp0令即所以(7.57)对于平稳随机过程，集合平均可用时间平均代替，因此上式可写成：(7.58)由式(7.52)、(7.53)、(7.54)不难看出(7.59)如果已知x(n)为有限长序列x0，x1，xn-1，(即x(0)，x(1)，x(n-1)，当p=1时，则按式(7.58)可得待添加的隐藏文字内容2(7.60)而按式(7.47)及(7.50)可得(7.61)

32、 (7.62)将e1(n)及b1(n)代入式(7.58)，又可得(7.63)再代入式(7.47)及(7.50)又可得e2(n)及b2(n)： 这里(见式(7.46) 再将e2(n)与b2(n)代入式(7.58)又可求得k3，将k3与e2(n)及b2(n)代入式(7.47)与(7.50)又可求得e3(n)与b3(n)以此类推，可直接从时间序列x(n)求得各阶的kp以及前向与后向误差ep(n)与bp(n)。将各kp(=app)代入式(7.46)又可求得各apk。于是将这些求得的ar系数代入式(7.8)(其中)即可求得功率谱的估计值。注意，为了使式(7.58)求和的x(n)值不超出已知的n个x(n)

33、(x0，x1，xn-1)的范围，在p=1时求和的上下限应为n=1到n-1(见式(7.60)，在p=2时求和的上下限应从n=2开始到n-1(见式(7.63)，p阶时的n应取(5.128)burg递推算法，对于短的时间序列x(n)仍能得到较正确的估计，因此得到普遍应用。但burg算法受levinson关系式(7.46)的约束，(仍应用了自相关矩阵的toeplitz性质，实际上只有开始无终的平稳随机序列才有这种性质)，因此不能完全克服levinson算法中的缺点，仍存在某些谱线分裂与频率偏移现象。综上，burg法估计ar(p)模型参数的具体步骤为：1、确定初始条件：2、按照公式 ， 计算kp。3、按

34、照公式和，计算ep(n)和bp(n)。4、计算均方误差：。5、p=p+1。6、重复第25步，直至满足条件为止。例7-2、设n=5的数据记录为x(0)=1,x(1)=2,x(2)=3,x(3)=4,x(4)=5,ar模型的阶次p=3,试用相关函数法确定ar参量及预测值。 .解：先由数据求自相关函数式： 用levinson-durbin递推算法求ar模型的参量分别是：根据所得的ar(3)参量，预测值：若使用的是二阶线性预测器，有例51所得的结果，则可分别由前向与后向预测得到如下： 例7-3、设仍利用例72中的记录数据，试用伯格法求ar(2)的参量。解：用上述递推公式，i=1时：e1(n)和b1(n) p=2时： 若使用此二阶线性预测，可得： 算法比较 levinson-durbin burg算法 真实值：前向误差功率 25.9090 3.0650：后向误差功率 25.5403 0.17415 1.2700 0.8549 1.0000 2.89834.5825 5.0000从以上比较看出：显然，伯格算法要比莱文森德宾算法优越得多。短数据！例子：现代谱估计和经典谱估计方法的比较。比较welch方法和bur