第3章 CHAPTER 3 信号参数估计与信号 检测基础 3.1信号参数估计的基本概念 信号处理的很多问题以参数估计的形式存在,例如在雷达及声呐信号处理中,需要估计目标(飞机或船只)的距离和到达角度,这些参数可以通过对发射信号和回波信号的处理获得。又如在现代的高保真录制系统中,需要精确保持录音基准电平,这一般通过电平自动控制系统实现,具体流程是周期性地估计该基准电平,并根据估计结果做出调整。 信号处理中的参数估计一般基于观测信号与待估计参数之间的依赖关系,在不同的问题模型中,这种依赖关系可能是确定性的或随机性的。例如在通信中的信道估计问题中,发送信号x(n)经过冲激响应为h(n)的信道输出为 y(n)=h(n)x(n)(31) 如果系统采用基于训练序列的信道估计方法,x(n)为已知,这时信道输出y(n)与信道参数h(n)存在确定性依赖关系; 如果系统信道估计不采用训练序列,而是基于信道输出的统计特性,则y(n)与信道参数h(n)存在随机性依赖关系,具体而言,就是y(n)的概率分布函数与h(n)存在依赖关系。随机性依赖关系与确定性依赖关系相比,更贴近实际环境。这首先是因为实际环境中观测信号的获取总是受到干扰,如环境干扰以及观测设备自身的热噪声干扰等。考虑噪声影响后的信道接收信号一般表示为 y(n)=h(n)x(n)+v(n)(32) 其中,v(n)表示系统环境干扰或设备热噪声产生的加性干扰,一般情况下用零均值、独立同分布的高斯随机信号描述。这时y(n)也体现为具有高斯分布特性的随机信号,h(n)对y(n)的影响也体现为对其概率密度分布和统计特性的影响。 随机性依赖关系更具重要性的另一原因是实际系统中待估计参数经常具有时变性,例如雷达信号处理中目标距离就是不断变化的参数; 移动通信中,发送端或接收端的移动会造成信道环境的改变。待估计参数的变化不是毫无规律,而是服从一定的概率密度分布,例如非视距移动信道的时变性服从瑞利分布。在式(31)中,如果h(n)是具有一定概率结构的随机信号,则无论x(n)是确定性信号还是随机信号,h(n)对y(n)的影响都体现为概率密度分布和统计特性的影响。 下面以一个简单的例子说明参数估计的过程。 例31在电平自动控制系统的电平估计问题中,对某个电平值θ进行估计,由于环境干扰及测量设备的不完善,测量总会存在误差,测量误差可归结为噪声,因此,实际得到的测量值为 x=θ+w(33) 其中,w一般服从零均值高斯分布,方差为σ2。问题是,如何根据测量值x估计θ的值。 一种最简单的方法是将多个观测结果进行平均。首先定义一个时间窗口,认为实际电平值在这个时间窗口中是恒定的。在这个时间窗口中得到电平值的N个含噪测量结果记为x(0),x(1),…,x(N-1)。这些测量结果中 x(n)=θ+w(n),n=0,1,…,N-1(34) 根据x(n)得到θ的估计为 θ^=x(0)+x(1)+…+x(N-1)N(35) 式(35)的一个直观效果是通过对噪声求均值,以降低噪声的影响。一般情况下,可以认为估计结果与观测数据由一个映射关系确定: θ^=f(x(0),x(1),…,x(N-1))(36) 上面提到,由于噪声的影响,x(n)更多时候体现为随机信号的特征,因此式(36)中x(0),x(1),…,x(N-1)可看成随机变量。θ^与x(0),x(1),…,x(N-1)存在映射关系,可认为θ^也是随机变量(不管θ是常量还是随机变量)。一次估计的结果可以看成θ^的一个样本,因此,验证估计算法的性能往往需基于多次估计结果进行评估。实际应用中往往很难得到类似于式(35)的闭式表达式,参数估计可能包含复杂的步骤,但包含的要素是一致的,从上面这个例子可以看出构造一个估计问题的基本要素如下: (1) 待估计参数: 指估计的目标参数,可以是单参量,也可以是多参量。待估计的参数可能是常量,也可能是随机变量,如果是随机变量一般要求其概率密度函数已知。 (2) 信号模型: 指观测信号的产生模型,其中必须体现待估计参数与观测信号的依赖关系。观测信号可以是单参量,也可以是多参量。 (3) 误差准则: 指估计算法的性能衡量指标,很多时候误差准则可指导估计算法的设计。误差准则的确定也有助于不同估计算法间的比较,特别是针对同一个估计问题的不同估计算法。 (4) 估计算法: 对同一估计问题,估计算法不是唯一的,不同的估计算法可以得到不同的估计性能。如何得到误差小、计算量少的估计算法是信号处理的核心问题。 3.2估计算法的性能指标 依据不同的估计准则,可以得到不同的估计量,这些估计量的性能需要有评价的指标。估计量是观测的函数,而观测的是随机变量(或向量),因此,估计量也是随机变量,随着观测数据的增多,希望估计量能逐渐逼近真值,因此,对估计量的均值和方差应有一定的要求。 3.2.1性能指标 估计量的好坏可以从无偏性、有效性、一致性加以评价。 1. 无偏性 待估计量为常数θ,它的估计值为θ^,如果E(θ^)=θ,则称θ^为θ的无偏估计,否则称θ^为有偏估计。定义估计的偏差为 θ~=E(θ^)-θ(37) 通常希望估计量的均值趋于被估计量的真值或被估计量的均值,即估计应该是无偏的。如果估计θ^不是无偏估计,但随着样本数目的增加,其数学期望趋近于真的估计量,即 limN→∞[E(θ^)-θ]=0(38) 则称估计θ^是渐近无偏估计。 当被估计量θ是随机变量时,如果估计量的均值等于被估计量的均值,即 E[θ^]=E[θ](39) 则称θ^为无偏估计。对于有偏估计,如果观测数据越多,估计的性能越好,即 limN→∞[E(θ^)-E(θ)]=0(310) 则称θ^为随机变量的渐近无偏估计。 2. 有效性 估计量具有无偏性并不表明已经保证了估计的品质。当被估计量为未知常数时,不仅希望估计量的均值等于真值,而且希望估计量的取值集中在真值附近,这一品质可以通过估计的方差来描述,估计的方差为 Var(θ^)=E{[θ^-E(θ^)]2}(311) 对于无偏估计,方差越小,表明估计量的取值越集中,估计的性能越好,估计也越有效。 对于有偏估计,估计的方差小并不能说明估计是好的,因为如果估计有偏差,方差小的估计仍然可能有较大的估计误差,这时用均方误差加以描述更合适。估计的均方误差定义为 Mse(θ^)=E{|θ^-θ|2}(312) 均方误差越小,表明估计越有效。 3. 一致性 当用N个观测值估计参量时,一般来说,观测值越多,估计越趋于真值,如果 limN→∞P{θ-θ^<ε}=1(313) 则称θ^为一致估计(式(313)中,ε是任意小的正数)。 3.2.2随机信号均值及自相关函数的估计 随机信号均值和自相关函数的估计是最常用的参数估计。这一节将讨论常用均值、自相关估计算法的无偏性、有效性和一致性,这里只讨论实随机信号。 1. 均值的估计 对遍历平稳随机序列,设x(0),x(1),…,x(N-1)是观察到的N个样本,则可利用式(314)进行均值的估计 m^x=1N∑N-1n=0x(n)(314) 此估计是无偏的,且当各子样本x(0),x(1),…,x(N-1)互不相关时,此估计是一致估计,现证明这两个结论。 证明: 估计的均值为 E[m^x]=E1N∑N-1n=0x(n)=1N∑N-1n=0E[x(n)]=1N∑N-1n=0mx=mx(315) 可见估计是无偏的。 估计的均方值为 E[m^2x]=E1N∑N-1n=0x(n)1N∑N-1m=0x(m) =1N2∑N-1n=0∑N-1m=0E[x(n)x(m)] =1N2∑N-1n=0E[x2(n)]+∑N-1n=0∑N-1m=0,m≠nE[x(n)x(m)] =1N2∑N-1n=0E[x2(n)]+∑N-1n=0∑N-1m=0,m≠nE[x(n)]E[x(m)](316) 记E[x2(n)]=Dx,E[x(n)]=mx,所以 E[m^2x]=1NDx+N(N-1)N2m2x =1NDx+N-1Nm2x(317) 所以,估计的方差为 Var(m^x)=E(m^2x)-[E(m^x)]2 =1NDx+N-1Nm2x-m2x=1Nσ2x(318) 可见,当N→∞时,估计的方差趋于零,所以m^x是一致估计。 2. 自相关函数的均值 在实际工作中,对于平稳随机序列x(n),所能得到的x(n)的N个观察值为x(0),x(1),…,x(N-1),记为xN(n)。将n≥N时的x(n)的值假设为零。现在的任务是如何由这N个观察值估计出x(n)的自相关函数Rx(m)。对Rx(m)的估计方法通常有两种: 一是直接估计法,二是利用FFT求解。下面以自相关函数的直接估计法为例分析其性能。 如果观察值的点数N为有限值,则通过直接估计法 求R(m)估计值为 R^(m)=1N∑N-1n=0xN(n)xN(n-m)(319) 由于x(n)只有N个观察值,因此,对于每一个固定的迟延m,在(0,N-1)内可以利用的数据只有(N-|m|)个,所以在实际计算R^(m)时,式(319)变为 R^(m)=1N∑N-1-|m|n=0x(n)x(n-m)(320) 值得注意的是,R^(m)的长度为2N-1,对于实信号,它是以m为对称中心的偶对称函数。所以在式(320)的计算中,可先计算m≤0的部分,再根据对称性求m>0的部分。 现在讨论R^(m)估计的质量。 1) 估计的偏差 由偏差的定义,有 Bia[R^(m)]=E[R^(m)]-R(m)(321) 式中 E[R^(m)]=E1N∑N-1-|m|n=0x(n)x(n-m) =1N∑N-1-|m|n=0E[x(n)x(n-m)] =1N∑N-1-|m|n=0R(m) 即 E[R^(m)]=N-|m|NR(m)(322) 所以估计的偏差为 Bia[R^(m)]=-mNR(m)(323) 分析式(322)和式(323)可以看出,估计不是无偏估计,而且当m→N时,估计的偏差很大。对于有限的m,当N→∞时,有E[R^x(m)]=Rx(m),所以估计是渐进无偏的。由此,可以得到如下结论: (1) 对于一个固定的延迟m,当N→∞时,Bia[R^(m)]→0,因此,R^(m)是对R(m)的渐近无偏估计; (2) 对于一个固定的观察值N,只有当|m| <<N时,R^(m)的均值才接近于真值R(m),即当m越接近于N时,估计的偏差越大; (3) 由式(322)可以看出,R^(m)的均值是真值R(m)和三角窗函数w(m) w(m)=N-|m|N,0≤|m|≤N-1 0,|m|≥N(324) 图31三角窗函数w(m) 的乘积,w(m)的长度是2N-1(如图31所示)。式(324)的三角窗函数又称Bartlett窗函数,由于它对R(m)进行了加权,致使R^(m)产生了偏差。显然,由于加权窗函数的存在,产生了上述第(2)点结论。 该窗函数实际上是由于对真实数据的截短而产生的。因为xN(n)可以看作真实数据x(n)和一矩形窗函数d(n)相乘的结果,即 xN(n)=x(n)d(n)(325) 式中 d(n)=1,0≤n≤N-1 0,n<0,n≥N(326) 根据式(320)和式(325),可得 R^(m)=1N∑N-1-|m|n=0xN(n)xN(n-m) =1N∑N-1-|m|n=0x(n)d(n)x(n-m)d(n-m)(327) E{R^(m)}=1N∑N-1-|m|n=0E{x(n)x(n-m)}d(n)d(n-m) =R(m)N∑N-1-|m|n=0d(n)d(n-m) =R(m)w(m)(328) 可见,w(m)正是对矩形数据窗函数d(n)作自相关的结果。 由上面的讨论可知,当对一个信号做自然截短时,就不可避免地对该数据施加了一个矩形窗口,由此数据窗口就产生了加在自相关函数上的三角窗口。三角窗口影响了R^(m)对R(m)估计的质量。加在数据上的窗口一般称为数据窗,加在自相关函数上的窗口一般称为延迟窗,后面将看到,这些窗函数也直接影响了功率谱估计的质量。 2) 估计的方差 根据方差的定义,有 Var[R^(m)]=E[R^(m)-E[R^(m)]]2 =E[R^2(m)]-{E[R^(m)]}2(329) 由式(322)可得 {E[R^(m)]}2= N-|m|NR(m)2(330) 同时可得 E[R^2(m)]=E1N2∑N-1-|m|n=0x(n)x(n-m)∑N-1-|m|k=0x(k)x(k-m) =1N2∑n∑kE[x(n)x(k)x(n-m)x(k-m)](331) 这里要计算随机信号x(n)的四阶矩,但一般平稳信号的四阶矩很难计算,因此假定x(n)是零均值的高斯随机信号,可以证明,其四阶矩满足如下关系式: E[x(n)x(k)x(n-m)x(k-m)]=E[x(n)x(k)]E[x(n-m)x(k-m)]+ E[x(n)x(n-m)]E[x(k)x(k-m)]+ E[x(n)x(k-m)]E[x(k)x(n-m)] =R2(n-k)+R2(m)+ R(n-k+m)R(k-n+m)(332) 所以 E[R^2(m)]=1N2∑n∑k[R2(n-k)+R2(m)+R(n-k+m)R(k-n+m)] = N-|m|NR(m)2+1N2∑n∑k[R2(n-k)+ R(n-k+m)R(k-n+m)](333) 将式(330)和式(333)代入式(329),可得 Var[R^(m)]=1N2∑N-1-|m|n=0∑N-1-|m|k=0[R2(n-k)+R(n-k+m)R(k-n+m)](334) 对于任意序列g(n),式(335)成立 ∑N-1-|m|n=0∑N-1-|m|k=0g(n-k)=∑N-1-|m|i=-(N-1-|m|)(N-|m|-|i|)g(i)(335) 因此,令n-k=i,可把式(334)的双求和变成单求和,即对于一个确定的i,n-k=i的组合数为N-|m|-|i|。 Var[R^(m)]=1N∑N-1-|m|i=-(N-1-|m|)1-|m|+|i|N[R2(i)+R(i+m)R(i-m)](336) 显然,当N→∞时,Var[R^(m)]→0。又因为limN→∞Bia[R^(m)]=0,所以可得如下结论: 对于固定的延迟m,R^(m)是R(m)的一致估计。对R(m)的另一种直接估计法是对式(320)稍做修改,即 R^(m)=1N-|m|∑N-1-|m|n=0xN(n)xN(n+m)(337) 可以证明,式(337)的R^(m)是对R(m)的无偏估计,利用上述类似的推导,可分析估计子的性能,最终可得 有偏估计子的方差: Var[R^(m)]∝O1N 无偏估计子的方差: Var[R^(m)]∝O1N-|m| 可以看到无偏估计子的方差性能较差。实际应用一般采用有偏估计子,除了方差的原因外,还基于以下考虑: ①当m取值较小时,有偏估计子的偏差不明显; ②对于大多数平稳随机信号,当m值较大时,R(m)的值一般很小。式(323)显示偏差的值也和R(m)的取值有关系,而m较大时,R(m)的取值一般会变小,这时可以认为偏差也比较小。 3.3估计性能界——CRB 不同的估计方法可以得到不同的估计量,估计量的性能可以通过前面介绍的无偏性、有效性和一致性评价,但在实际中,估计量可能比较复杂,很难评价估计量的有效性和一致性。此外,在得到一个估计量后,它的性能是否已经达到最佳?是否还有更好的估计量?克拉美罗界(CramérRao Bound,CRB)揭示了无偏估计量估计方差的最小值。 3.3.1单参数实常量估计的CRB 首先考虑只有一个待估计参数的情况,令θ为待估计参数,并且为实常量(非随机量)。观测信号x(n)与θ的依赖关系用条件概率p(x|θ)表示[为描述方便,以下将x(n)简记为xn],有如下定理。 定理31令x=(x0,x1,…,xN-1)为观测样本向量,p(x|θ)是x的条件概率密度。若θ^是θ的一个无偏估计子,且f(x|θ)/θ存在,则θ^的方差满足 Var(θ^)=E[(θ^-θ)2]≥I-1(θ)(338) 式中 I(θ)=Elnp(x|θ)θ2=-E2lnp(x|θ)θ2(339) 当且仅当 lnp(x|θ)θ=K(θ)(θ^-θ)(340) 时,式(338)的等号成立。其中,K(θ)是θ的正函数,且不包含x。 证明: 由假设条件知E(θ^-θ)=0,因此 E[θ^-θ]= ∫+∞-∞…∫+∞-∞(θ^-θ)f(x|θ)dx0…dxN-1 =Δ∫+∞-∞(θ^-θ)f(x|θ)dx =0(341) 在式(341)两边对θ求偏微分,可得 θE[θ^-θ]= ∫+∞-∞θ[(θ^-θ)f(x|θ)]dx=0(342) 由式(342)得 ∫+∞-∞-f(x|θ)dx+ ∫+∞-∞(θ^-θ)θf(x|θ)dx=0(343) 因为p(x|θ)是概率密度函数,所以 ∫+∞-∞f(x|θ)dx=1(344) 同时根据组合函数导数性质可得 θf(x|θ)=θlnf(x|θ)·f(x|θ)(345) 式(343)重写为 ∫+∞-∞θlnf(x|θ)·(θ^-θ)f(x|θ)dx=1(346) 式(346)进一步重写为 ∫+∞-∞θlnf(x|θ)f(x|θ)·(θ^-θ)f(x|θ)dx=1(347) 根据泛函分析中的柯西不等式|〈z1·z2〉|2≤〈z1·z1〉·〈z2·z2〉,其中z1,z2是泛函空间的元素,〈·〉表示内积操作。对于式(347),令 z1=f1(x)=θlnf(x|θ)f(x|θ) z2=f2(x)=(θ^-θ)f(x|θ) 内积操作定义为∫+∞-∞f1(x)f2(x)dx。由式(347)得 ∫+∞-∞θlnf(x|θ)2f(x|θ)dx·∫+∞-∞(θ^-θ)2f(x|θ)dx≥1(348) 等价于 ∫+∞-∞(θ^-θ)2f(x|θ)dx≥1∫+∞-∞θlnf(x|θ)2f(x|θ)dx(349) 式(349)左边刚好就是无偏估计的方差计算结果。柯西不等式等号成立的条件是z1,z2成比例,即式(340)成立。 证毕。 特别要指出的是,式(342)积分和求导顺序可以互换,这在实际中一般可以满足。为了得到更直观的描述,将在式(344)两边对θ求导数,并将积分和求导顺序互换可得 ∫+∞-∞θf(x|θ)dx=0(350) 基于式(345)可得 ∫+∞-∞θlnf(x|θ)·f(x|θ)dx=Eθlnf(x|θ)=0(351) 式(351)称为定理31成立的规范化条件,是积分和求导顺序可以互换的直观描述。 在式(351)两边对θ求导数,则可得式(339)。 CRB给出了无偏估计量估计方差的下限。达到CRB的估计,其估计的方差是最小的,称这样的估计为有效估计。I(θ)称为数据的费希尔(Fisher)信息,CRB是费希尔信息的倒数。直观理解是,费希尔信息越多,CRB越低。 费希尔信息具有信息度量的基本性质。首先,由式(339)可以看出,费希尔信息具有非负性; 其次,对于独立观测,费希尔信息具有可加性,这是因为,对于独立的观测,有 lnp(x|θ)=∑N-1n=0lnp(xn|θ)(352) 所以 -E2lnp(x|θ)θ2=-∑N-1n=0E2lnp(xn|θ)θ2(353) 如果观测是独立同分布的,那么 I(θ)=Ni(θ)(354) 式中 i(θ)=-E2lnp(xn|θ)θ2(355) 对于非独立观测,费希尔信息要比独立观测低,例如,对于x0=x1=…=xN-1这种完全相关的情况,I(θ)=i(θ),此时,CRB并不随观测数据的增加而降低。 例32独立同分布平稳随机信号服从高斯分布,用式(314)估计该随机信号的均值,该估计子的方差是否达到CRB? 解平稳高斯信号的一维概率密度为 p(x|mx)=12πσexp-(x-mx)22σ2(356) 多变量的联合概率密度为 p(x|mx)= 12πσN∏N-1n=0exp-(xn-mx)22σ2(357) 对数联合概率密度的一阶、二阶偏导数为 lnp(x|mx)mx=∑N-1n=0xn-mxσ2(358) 2lnp(x|mx)m2x=-Nσ2(359) 可得到其费希尔信息量为 I(mx)=-E2lnp(x|mx)m2x=Nσ2(360) 对比3.2.2节中估计子的方差,可以看到该均值估计可达到CRB。注意到式(358)中的一阶导数可重写为 lnp(x|mx)mx=Nσ2(m^x-mx)(361) 即一阶导数和(m^x-mx)成正比,显示该估计子满足式(340)的条件。 例33信号观测模型为xn=Acos(2πf0n+)+w(n),其中,f0为已知参数,w(n)∝N(0,σ2),相位为未知参数。求相位估计的CRB。 解观测信号的一维概率密度为 p(x|)=12πσexp-(xn-Acos(2πf0n+))22σ2(362) 假设参数估计基于N个观测,记为x=[x0, x1, …, xN-1]T,其联合概率密度函数为 p(x|)=p(x0,x1,…,xN-1|) = 12πσN∏N-1n=0exp-(xn-Acos(2πf0n+))22σ2(363) 可得到观测信号的对数似然函数为 lnp(x|)=-N2ln(2πσ2)-12σ2∑N-1n=0(xn-Acos(2πf0n+))2(364) 可得 lnp(x|)=-1σ2∑N-1n=0[xn-Acos(2πf0n+)]Asin(2πf0n+) =-Aσ2∑N-1n=0[xnsin(2πf0n+)-A2sin(4πf0n+2)] 2lnp(x|)2=-Aσ2∑N-1n=0[xncos(2πf0n+)-Acos(4πf0n+2)]-E2lnp(x|)2 =Aσ2E∑N-1n=0[xncos(2πf0n+)-Acos(4πf0n+2)] =Aσ2∑N-1n=0E[xn]cos(2πf0n+)-Acos(4πf0n+2) =Aσ2∑N-1n=0Acos2(2πf0n+)-Acos(4πf0n+2) =A2σ2∑N-1n=012+12cos(4πf0n+2)-cos(4πf0n+2) =A22σ2∑N-1n=01-cos(4πf0n+2) =NA22σ2-A22σ2∑N-1n=0cos(4πf0n+2)(365) 式(365)中求和项为周期序列,有 limN→∞1N∑N-1n=0cos(4πf0n+2)=0(366) 当N比较大时,可认为 1N∑N-1n=0cos(4πf0n+2) 1(367) 于是可认为式(365)中第二项远小于第一项,直接忽略第二项可得 -E2lnp(x|)2=NA22σ2(368) 其倒数即为相位估计的CRB。 3.3.2多参量估计的CRB 假设有L个未知参数,用向量表示为 θ=[θ1,θ2,…,θL]T(369) 参数估计的CramerRao不等式表示为 E[(θ^-θ)(θ^-θ)T]≥I-1(θ)(370) 其中,符号“≥”用于矩阵时有不同的含义。对于矩阵A,B,A≥B表示A-B是半正定矩阵,A>B表示A-B是正定矩阵。式(370)中左边为估计子的协方差矩阵,右边的I(θ)称为费希尔信息矩阵,定义为 I(θ)=-E2lnp(x|θ)θθT(371) 式(371)中用到了向量函数的求导公式,下面给出其定义。设f(θ)是L×1向量θ的函数,函数值为标量,其导数表示为 f(θ)θ= f(θ)θ1 f(θ)θ2 ︙ f(θ)θLL×1(372) 根据该定义可得以下常用导数 aTθθ=θTaθ=a(373) θTAθθ=(AT+A)θ(374) 式(372)的定义可拓展到多个函数的组合,假设有R个关于θ的函数构成如下向量: f(θ)=f1(θ) f2(θ) ︙ fR(θ)R×1(375) 则其关于θ的导数定义为 f(θ)θT=f1(θ)θ1f1(θ)θ2…f1(θ)θL f2(θ)θ1f2(θ)θ2…f2(θ)θL ︙︙︙ fR(θ)θ1fR(θ)θ2…fR(θ)θLR×L(376) 根据上述定义,可得矩阵I(θ)中第{i,j}个元素为 [I(θ)]ij=-E2lnp(x|θ)θiθj(377) 对式(371)这里不做详细证明,下面以一个例子介绍多参量CRB的计算。 例34对独立同分布平稳高斯随机信号,计算均值估计和方差估计的CRB。 解由于待估计向量θ=[mxσ2]T,因此费希尔信息矩阵是2×2的矩阵,即 I(θ)=-E2lnp(x|θ)θθT=-E2lnp(x|θ)m2x-E2lnp(x|θ)mxσ2 -E2lnp(x|θ)σ2mx-E2lnp(x|θ)(σ2)2(378) 由式(356)得到对数似然函数为 lnp(x|θ)=-N2ln(2πσ2)-12σ2∑N-1n=0(xn-mx)2(379) 所以 lnp(x|θ)mx=1σ2∑N-1n=0(xn-mx)(380) lnp(x|θ)σ2=-N2σ2+12σ4∑N-1n=0(xn-mx)2(381) 2lnp(x|θ)m2x=-Nσ2(382) 2lnp(x|θ)(σ2)2=N2σ4-1σ6∑N-1n=0(xn-mx)2(383) 2lnp(x|θ)mxσ2=-1σ4∑N-1n=0(xn-mx)(384) 2lnp(x|θ)σ2mx=2lnp(x|θ)mxσ2(385) 费希尔信息矩阵为 I(θ)=-E-Nσ2-E-1σ4∑N-1n=0(xn-mx) -E-1σ4∑N-1n=0(xn-mx) -EN2σ4-1σ6∑N-1n=0(xn-mx)2=Nσ20 0N2σ4(386) 因此,容易得出参数估计的CRB为 Var(m^x)≥σ2N(387) Var(σ^2)≥2σ4N(388) 3.3.3参数变换的CRB 在实际中,经常遇到希望估计的参数是某个参数的函数,例如,在电平估计中,有时感兴趣的不是信号电平A的估计,而是信号的功率A2。假设估计A的CRB已知,如何得到A2估计的CRB呢?下面回答这一问题。 假定可以直接计算CRB的参数为θ,希望估计的参数为α,两者有如下函数关系: α=T(θ)。那么可以证明α^的CRB为 Var(α^)≥T(θ)θ2-E2lnp(x|θ)θ2(389) 例如,基于例34中式(387)的结果,求α=m2x的CRB,则 Var(α^)≥(2mx)2N/σ2=4m2xσ2N(390) 实际应用中经常碰到多变量的映射关系。例如,独立同分布高斯随机信号的平均功率D=m2x+σ2。前面已经得到估计mx、σ2的CRB,在此基础上可采用以下方法得到估计D的CRB。 假设P×1维向量θ包含模型中的可估计参数,有多个待估计参数,可表示为θ的函数,记为α=T(θ),其中T(θ)是R×1维的函数向量,即 α=[T1(θ)T2(θ)…TR(θ)]T(391) 则 E[(α^-α)(α^-α)T]≥T(θ)θTI-1(θ)T(θ)θTT(392) 式中,T(θ)/θT的定义见式(376)。 例35对高斯白噪声中恒定电平,求信噪比估计的CRB。 解高斯白噪声中恒定电平观测的数学模型为xn=A+wn,其中A为电平值,wn∝(0,σ2)为零均值高斯白噪声,信噪比定义为α=A2/σ2,因此可先计算估计A、σ2的CRB。注意到观测数据可以看成 非零均值高斯信号,A、σ2的估计对应均值和方差的估计,可采用例34的结果。设θ=[Aσ2]T,T(θ)=A2/σ2,由例34可知 I(θ)=Nσ20 0N2σ4(393) 又因为 T(θ)θT=T(θ)AT(θ)σ2=2Aσ2-A2σ4(394) 所以 T(θ)θTI-1(θ)T(θ)θTT=2Aσ2-1σ4σ2N0 02σ4N2Aσ2 -1σ4(395) =4A2Nσ2+2Nσ4(396) 可以证明,对向量参数的有效估计量经过线性变换后仍然是有效估计量,而向量参数的有效估计量经过非线性变换后是渐近有效估计量。 3.3.4复参数估计的CRB 当待估计参数是复数时,推导涉及复变函数的微分操作,比实参数估计复杂得多。这一节先给出复变函数的基本概念,再介绍一种复参数估计的CRB的推导方法。 令复自变量为z=x+jy,x,y∈R。定义在复数域的函数为 f(z)=u(z,z)+jv(z,z)=u(x,y)+jv(x,y)(397) 其中,u(x,y),v(x,y)为实二元函数。f(z)的导数定义为 fz=12fx-jfy fz=12fx+jfy(398) 复变函数可导的充分必要条件可参考相关复变函数文献。这里假设所涉及的复变函数均可导,以下是一些常用的性质。 假设g(z)、f(z)在区域D内可导,则其导数满足以下性质: ddz(f(z)+g(z))=df(z)dz+dg(z)dz(399) ddz(f(z)·g(z))=df(z)dzg(z)+dg(z)dzf(z)(3100) ddzf(z)g(z)=1g2(z)df(z)dzg(z)-dg(z)dzf(z)(3101) ddz(f(g(z)))=df(g)dgdg(z)dz(3102) 对多变量函数,如f(z),其中z=[z1,z2,…,zN]T为包含所有自变量的向量,其导数定义为 f(z)z= f(z)z1,f(z)z2,…,f(z)zNT(3103) 对多维函数,可定义与式(370)相同的偏导数。根据该定义可得以下对复向量的常用导数。 aHzz=a*(3104) zHAzz=ATz*(3105) 接下来考虑复参数估计的CRB,假设观测数据为Z=[z0,z1,…,zN-1],观测数据的信号产生模型已知,但其中有L个未知参数,记为θ=[θ1,θ2,…,θL]T。待估计的目标参数为θ的函数,记为 p=[p1(θ),p2(θ),…,pK(θ)]T(3106) 记估计子为r(Z),考虑无偏估计子,有 ∫Ωr(Z)f(Z|θ)dZ=p(3107) 其中,Ω表示观测空间,在式(3107)两边对θ求偏导可得 ∫Ωr(Z)f(Z|θ)θHdZ=pθH(3108) 式(3108)左边重写为 ∫Ωr(Z)f(Z|θ)θHdZ=∫Ωr(Z)lnf(Z|θ)θHf(Z|θ)dZ =Er(Z)lnf(Z|θ)θH(3109) 于是 pθH=Er(Z)lnf(Z|θ)θH(3110) 简记r(Z)为r,lnf(Z|θ)为lnf,构造如下向量 r lnfθHH(K+L)×1(3111) 可算出其自协方差矩阵,并由自协方差矩阵的半正定特性可得 Cov(r,r)pθH pθHHElnfθHHlnfθH≥0(3112) 其中已经用到式(3110)的结论。根据自协方差矩阵特性,式(3113)成立 D-pθHI-1(θ)Cov(r,r)pθH pθHHI(θ)D-pθHI-1(θ)H≥0(3113) 其中,D为单位矩阵,进一步将费希尔信息量矩阵定义为 I(θ)=ElnfθHHlnfθH(3114) 经简单计算,式(3113)矩阵可化简为 Cov(r,r)-pθHI-1(θ)pθHH≥0(3115) 可以得到估计子的误差界为 Cov(r,r)≥pθHI-1(θ)pθHH(3116) 式(3116)和式(392)实参量估计的CRB是统一的。下面以一个例子介绍CRB的具体计算过程。 例36通信中常用以下信号观测模型: xn=Asn+wn,n=0,1,…,N-1(3117) 其中,A为复值信道响应,sn为正交调制发送信号,wn是零均值复高斯白噪声序列,方差为σ2,假设发送符号有归一化功率,即sn2=1,求信道估计的CRB。 解根据式(295)给出的复高斯概率密度分布,关于A的一维条件概率密度为 p(x|A)=1πσ2exp-xn-Asn2σ2(3118) 多变量的联合概率密度为 p(x|A)= 1πσ2N∏N-1n=0exp-xn-Asn2σ2(3119) 基于式(3114)费希尔信息量矩阵的定义及式(350)所示的规范化条件,可以得到与式(339)类似的结论。 I(θ)=Elnp(x|θ)θlnp(x|θ)θ=-E2lnp(x|θ)θθ(3120) 对数联合概率密度的一、二阶偏导数为 lnp(x|A)A=∑N-1n=0 snxn-Asnσ2(3121) 2lnp(x|A)AA=-1σ2∑N-1n=0sn2=Nσ2(3122) 可得到估计信道响应的性能界为 CRBA=σ2N(3123) 3.4最大似然估计 当被估计量为未知常量时,可以采用比较简单的最大似然估计。最大似然估计可以简便地完成对复杂估计问题的求解,而且当观测数据足够多时,其性能也是非常好的。因此,最大似然估计在实际中得到了广泛采用。本节只考虑实参数的估计。 3.4.1最大似然估计的基本原理 设观测向量x=[x0x1…xN-1]T,待估计参量为θ。观测数据与待估计量的概率依赖关系以概率密度p(x|θ)描述。概率密度是以θ为参量的函数,在观测给定的条件下,如x=x0,p(x=x0|θ)反映了θ取各个值的可能性大小,称p(x=x0|θ)为θ的似然函数,p(x=x0|θ)简记为p(x|θ)。估计问题本质上就是根据观测求出未知量,也就是说,在得到某个观测值x=x0后,如何根据这个观测确定未知量θ的值。 从另外一个角度,似然函数p(x=x0|θ)反映了对不同的θ值观测到x=x0的概率。因为当前的观测结果刚好就是x0,于是可认为实际中观测到x=x0的概率应该尽可能大,即θ取值应该使得p(x=x0|θ)取得最大。使似然函数最大所对应的参数θ作为对θ的估计,称为最大似然 (Maximum Likelihood,ML) 估计,记为θ^ml,即 θ^ml=argmaxθ p(xθ)或θ^ml=argmaxθ ln p(x|θ)(3124) 下面以两个例子说明最大似然准则的应用。 例37设实随机变量满足以下分布: x∝N(0,σ2),有N次独立观测x0,…,xN-1,求σ2的最大似然估计。 解单变量似然函数为 p(x|σ2)=12πσ212exp-x22σ2(3125) 多变量联合似然函数为 p(x|σ2)=12πσ2N2exp-12σ2∑N-1n=0x2n(3126) 对数联合似然函数为 lnp(x|σ2)=-N2ln(2πσ2)-12σ2∑N-1n=0x2n(3127) 令 lnp(x|σ2)σ2=0(3128) 得到估计结果为 σ^2ml=1N∑N-1n=0x2n(3129) 很容易验证 2lnp(x|σ2)(σ2)2σ2=1N∑N-1n=0x2n<0(3130) 所以式(3129)求得的是最大值。 例38假定观测序列为 xn=Acos(2πf0n+)+wn,n=0,1,…,N-1(3131) 幅度A和频率f0是已知的,wn是零均值高斯白噪声序列,方差为σ2,求相位的最大似然估计。 解似然函数为 p(x|)=1(2πσ2)N2exp-12σ2∑N-1n=0[xn-Acos(2πf0n+)]2(3132) 对数似然函数及其导数为 lnp(x|)=-N2ln(2πσ2)-12σ2∑N-1n=0[xn-Acos(2πf0n+)]2(3133) lnp(x|)=-1σ2∑N-1n=0[xn-Acos(2πf0n+)]Asin(2πf0n+)(3134) 令式(3134)等于零,可得 ∑N-1n=0xnsin(2πf0n+^ml)=A∑N-1n=0cos(2πf0n+^ml)sin(2πf0n+^ml) =2A∑N-1n=0sin(4πf0n+2^ml)(3135) 当f0不在0或1/2附近时,式(3135)右边是周期序列求和,当N足够大时,可认为其值近似为零。因此,最大似然估计近似满足 ∑N-1n=0xnsin(2πf0n+^ml)=0(3136) 展开式(3136),得 ∑N-1n=0xnsin2πf0ncos^ml=-∑N-1n=0xncos2πf0nsin^ml(3137) ^ml=-arctan∑N-1n=0xnsin2πf0n∑N-1n=0xncos2πf0n(3138) 前面讨论的是单参量的最大似然估计问题,在实际中经常需要同时估计多个参量,例如,在雷达信号的处理中,当检测到目标后,需要同时估计目标的位置、速度等参量,这样的问题称为多参量同时估计,或称为向量估计。式(3124)所描述的最大似然估计很容易推广到向量估计的情形。 考虑实参数估计问题,设有p个参量θ0,θ1,…,θp-1需要同时估计,定义一个p维向量θ=[θ0,θ1,…,θp-1]T,那么最大似然估计定义为 θ^ml=argmaxθ p(x|θ)或者θ^ml=argmaxθ lnp(x|θ)(3139) 类似地,最大似然估计的必要条件是 p(x|θ)θ|θ=θ^ml=0或者lnp(x|θ)θ|θ=θ^ml=0(3140) 例39设有N次独立观测xn=A+wn,其中wn∝N(0,σ2),A、σ2均未知,求A、σ2的最大似然估计。 解令θ=[Aσ2]T,则 p(x|θ)=12πσ2N2exp-12σ2∑N-1n=0(xn-A)2(3141) lnp(x|θ)=-N2ln(2πσ2)-12σ2∑N-1n=0(xn-A)2(3142) lnp(x|θ)θ=Nσ21N∑N-1n=0xn-A -N2σ4σ2-1N∑N-1n=0(xn-A)2(3143) 令lnp(x|θ)θθ=θ^ml=0,可求得最大似然估计为 θ^ml=A^ml σ^2ml=x- 1N∑N-1n=0(xn- x-)2(3144) 其中,x-表示样本的直接平均,定义为 x-=1N∑N-1n=0xn(3145) 以下例子给出更一般化线性估计模型中的最大似然估计算法。 例310假定待估计向量为θ=[θ0,θ1,…,θp-1]T,观测数据x=[x0,x1,…,xN-1]T,有如下的线性模型表示 x=Hθ+w(3146) 其中,H是N×p的矩阵,且N>p,H的秩为p,w是多维高斯噪声,其概率模型为w∝N(0,C),其中C为w的协方差矩阵。求θ的最大似然估计。 解由式(291)得似然函数为 p(x|θ)=1(2π)N2[det(C)]12exp-12(x-Hθ)TC-1(x-Hθ)(3147) 要使似然函数最大,只需使式(3148)的表达式最小 J(θ)=(x-Hθ)TC-1(x-Hθ)(3148) 基于式(373)、式(374),并利用协方差矩阵的对称性,将式(3148)对θ求导 J(θ)θ=-2HTC-1(x-Hθ) (3149) 令导数等于零,可求得最大似然估计,即 HTC-1(x-Hθ^ml)=0(3150) 解式(3150)的方程式得最大似然估计为 θ^ml=(HTC-1H)-1HTC-1x(3151) 由式(3151)可以看出,线性观测模型的最大似然估计是观测结果的线性函数,计算简单,因此在噪声协方差矩阵已知情况下是常用的估计算法。 3.4.2变换参数的最大似然估计 在许多情况下,希望得到估计参量θ的一个函数,例如α=T(θ),如果利用最大似然准则求得θ^,如何进一步求解α^?下面通过几个例子说明变换参数的最大似然估计的求法。 例311设有N次独立观测xn=A+wn,n=0,1,…,N-1。其中,wn∝N(0,σ2),A为未知参数,σ2为已知,求α=eA的最大似然估计。 解α的似然函数及其导数为 p(x|α)=12πσ2N2exp-12σ2∑N-1n=0(xn-lnα)2(3152) lnp(x|α)α=1ασ2∑N-1n=0(xn-lnα)=Nασ21N∑N-1n=0xn-lnα(3153) 令式(3153)等于零,可解得 α^ml=exp1N∑N-1i=0xi=exp(x-)(3154) 由于x-刚好是A的最大似然估计,所以 α^ml=exp(A^ml)(3155) 可见,α的最大似然估计只需要用A的最大似然估计代入变换式α=eA中就可以求得。因此,如果变换α=T(θ)是一一对应的,那么变换参数后的最大似然估计可以直接由式(3156)得到。 α^ml=T(θ^ml)(3156) 这一特性称为最大似然估计的不变性。 如果变换α=T(θ)不是一一对应的,则不能简单地应用式(3156)。下面通过一个例子加以说明。 例312在例311中,假定要求的估计为α=A2,求α的最大似然估计。 解由于A=±α,所以本例的变换不是一一对应的,似然函数p(x|α)需要两个概率密度描述 p1(x|α)=12πσ2N2exp-12σ2∑N-1n=0(xn-α)2,A≥0(3157) p2(x|α)=12πσ2N2exp-12σ2∑N-1n=0(xn+α)2,A<0(3158) 那么,最大似然估计为 α^ml=argmaxα {p1(x|α),p2(x|α)}(3159) 式(3159)的计算可以分为两步进行。 (1) 给定一个α值,如α=α0,比较p1(x|α0)和p2(x|α0)的大小,如果 p1(x|α0)≥p2(x|α0)(3160) 则pT(x|α0)=p1(x|α0),对所有α的取值重复以上过程,得到pT(x|α),称pT(x|α)为修正的似然函数。 (2) 求出使pT(x|α)最大的α值作为α的最大似然估计,即 α^ml=argmaxαpT(x|α)(3161) 对于本例,要比较p1(x|α0)和p2(x|α0)的大小,只需要比较∑N-1n=0(xn-α)2和∑N-1n=0(xn+α)2的大小,令 J(x,α)=∑N-1n=0(xn-α)2-∑N-1n=0(xn+α)2(3162) 式(3162)经整理得 J(x,α)=-4∑N-1n=0xnα=-4Nαx-(3163) 可见,如果x-≥0,则J(x,α)<0,p1(x|α0)≥p2(x|α0); 如果x-<0,则J(x,α)>0,p2(x|α0)>p1(x|α0),所以修正的似然函数为 pT(x|α)=p1(x|α),x-≥0 p2(x|α),x-<0(3164) 当x-≥0时,由pT(x|α)=p1(x|α)可求得α^ml=x-2; 当x-<0时,同样由pT(x|α)=p2(x|α)可求得α^ml=x-2。综合两种情况可得 α^ml= x-2(3165) 需要注意的是,如果将A的原始估计A^ml=x-直接代入变换式α=A2中,也可以得到式(3165),但这只是一种巧合。参数变换如果不是一一对应的,那么变换后参数的最大似然估计不能简单地通过代入变换式得到。 3.5贝叶斯估计 最大似然估计没有假定被估计量的先验信息,把待估计量看成一个常量。实际环境中待估计量可能是一个随机变量,并且概率密度已知。这时最大似然方法依然可用,但未能充分利用待估计量的信息。这一节介绍贝叶斯估计,它是一类能充分利用待估计量概率信息的估计方法。本节只考虑实参数的估计。 3.5.1代价函数 在估计某个参量θ时,噪声的影响会使估计产生误差,估计误差是要付出代价的,这种代价可以用代价函数加以描述,记为c(θ,θ^)。一般而言,估计误差Δθ=θ-θ^越小,代价越小。代价函数可表示为c(θ,θ^)、c(θ-θ^)或c(Δθ),典型的代价函数有以下几种。 (1) 平方代价函数 c(θ,θ^)= (θ-θ^)2(3166) (2) 绝对值代价函数 c(θ,θ^)=|θ-θ^|(3167) (3) 均匀代价函数 c(θ,θ^)=1,|θ-θ^|≥ε2 0,|θ-θ^|<ε2(ε为一常数)(3168) (4) 二次型代价函数 如果被估计量是p维随机向量,令误差向量为Δθ=θ-θ^,则二次型代价函数定义为 c(θ,θ^)=‖Δθ‖2S=(θ-θ^)TS(θ-θ^)(3169) 式中,‖Δθ‖为误差向量的范数,S为p×p维的对称非负定的加权矩阵。二次型代价函数实际上是估计误差的加权平方和,当S为单位矩阵时,二次型代价函数是估计误差的平方和。 代价函数也可以考虑其他函数,但实际上,最佳估计对代价函数的选择并不敏感。不失一般性,以下考虑单变量,代价函数确定后,可以计算平均代价 c-=E{c(θ,θ^(x))}= ∫∞-∞∫∞-∞c(θ,θ^(x))p(θ,x)dθdx(3170) 式(3170)中用θ^(x)代替θ^,强调θ^是从x计算得到,存在函数依赖关系。所谓贝叶斯估计就是使平均代价最小的估计。式(3170)可以写成 c-= ∫∞-∞∫∞-∞c(θ,θ^(x))p(θ|x)dθp(x)dx= ∫∞-∞c-(θ|x)p(x)dx(3171) 实际应用中p(x)可能是未知的或有复杂的表达式。在一次观测中,如果观测结果为x,则可认为其对应较大的概率。同时由于p(x)是非负的,所以使平均代价最小近似等价于使条件平均代价c-(θ|x)最小,也就是使 c-(θ|x)= ∫∞-∞c(θ,θ^(x))p(θ|x)dθ(3172) 最小。贝叶斯估计与代价函数的选取有关,使用不同的代价函数将得到不同的贝叶斯估计。 3.5.2最小均方误差估计 当代价函数为平方代价函数时,平均代价刚好等于估计误差的均方值。先考虑单个参数的估计问题,有 c-= ∫∞-∞∫∞-∞[θ-θ^(x)]2p(θ,x)dθdx=E{[θ-θ^(x)]2}(3173) 使平均代价最小等价于使均方误差最小,这时的贝叶斯估计称为最小均方误差估计,记为θ^ms。 把平方代价函数代入式(3172)中,得 c-(θ|x)= ∫∞-∞[θ-θ^(x)]2p(θ|x)dθ(3174) 在式(3174)两边对θ^求导,得 c-(θ|x)θ^=-2∫∞-∞[θ-θ^(x)]p(θ|x)dθ(3175) 并令导数在θ^=θ^ms处为零,得 θ^ms= ∫∞-∞θp(θ|x)dθ=E(θ|x)(3176) 由于条件平均代价对θ^的二阶导数为正,因此θ^ms所对应的条件平均代价为极小值,即θ^ms为最小均方估计。由式(3176)可以看出,最小均方误差估计为被估计量θ的条件均值。 如果被估计量为向量θ=[θ0,θ1,…,θp-1]T,将二次型代价函数代入式(3170),得 c-= ∫∞-∞∫∞-∞[θ-θ^(x)]TS[θ-θ^(x)]p(θ,x)dθdx =E{[θ-θ^(x)]TS[θ-θ^(x)]}(3177) 将二次型代价函数代入式(3172),得 c-(θ|x)= ∫∞-∞[θ-θ^(x)]TS[θ-θ^(x)]p(θ|x)dθ(3178) 对θ^(x)求导并令导数在θ^=θ^ms处为零,得 - ∫∞-∞2S[θ-θ^(x)]p(θ|x)dθ=0(3179) 求解可得 θ^ms= ∫∞-∞θp(θ|x)dθ=E(θ|x)(3180) 由于2c-(θ|x)θ^θ^T=2S,所以式(3180)求得的是极小值,θ^ms称为最小均方误差估计。 3.5.3条件中位数估计 当代价函数为绝对值代价函数时,条件平均代价为(考虑单变量情况) c-(θ|x)= ∫∞-∞θ-θ^(x)p(θ|x)dθ = ∫θ(x)-∞[θ^(x)-θ]p(θ|x)dθ+ ∫∞θ(x)[θ-θ^(x)]p(θ|x)dθ(3181) 在式(3181)两边对θ^求导,并令导数在θ^=θ^abs处为零,得 ∫θabs-∞p(θ|x)dθ= ∫∞θabsp(θ|x)dθ(3182) 由式(3182)可以看出,采用绝对值代价函数的贝叶斯估计θ^abs刚好是条件概率密度p(θ|x)的中位数(Median),所以也称为条件中位数估计,记为θ^med。 3.5.4最大后验概率估计 当采用均匀代价函数时,条件平均代价为 c-=c(θ|x)= ∫θ-Δ2-∞p(θ|x)dθ+ ∫∞θ+Δ2p(θ|x)dθ =1- ∫θ+Δ2θ-Δ2p(θ|x)dθ(3183) 很显然,当式(3183)中后面的积分最大时,平均代价最小。要使所述积分达到最大,当p(θ|x)具有类似于高斯分布的特征时,θ^可选为p(θ|x)的最大值所对应的θ(如图32(a)所示),这时对应的贝叶斯估计就是最大后验概率(Maximum A Posteriori,MAP)估计,记为θ^map。但是,当p(θ|x)不具备上述特征时,MAP不能得到最优的估计。如图32(b)所示,p(θ|x)具有指数分布的特征,这时θ^map显然不是合理的选择。 图32θ条件概率的不同形式 高斯信号是实际应用中最常见的信号,再加上MAP估计思路直观,实现也比较简便,MAP估计成为一种重要的参数估计准则。 下面进一步讨论最大后验概率估计与最大似然估计的联系。最大后验概率估计的数学描述为 θ^map=argmaxθ p(θ|x)或θ^map=argmaxθ ln p(θ|x)(3184) 由于 p(θ|x)=p(x|θ)p(θ)p(x)(3185) 当难以得到p(x)的表达式时,对于当前样本,认为p(x)有固定值,这时后验概率密度最大和使p(x|θ)p(θ)最大近似相同,因此,最大后验概率估计为 θ^map=argmaxθ p(x|θ)p(θ)或θ^map=argmaxθ[lnp(x|θ)+lnp(θ)](3186) 式(3186)直观给出了最大后验概率估计和最大似然估计的联系,可以看到最大后验概率估计实际上是在最大似然估计的基础上利用待估计参数的先验知识进行修正。当后验概率密度是可导的函数时,最大后验概率估计的必要条件是 p(θ|x)θθ=θmap=0或lnp(θ|x)θθ=θmap=0(3187) 式(3187)称为最大后验概率方程。 式(3187)所描述的最大后验概率估计很容易推广到向量估计的情形。假定待估计量为p维向量θ=[θ0,θ1,…,θp-1]T,后验概率密度为p(θ|x),那么向量最大后验概率估计为 θ^map=argmaxθ p(θ|x)或θ^map=argmaxθ ln p(θ|x)(3188) 3.5.5贝叶斯估计举例 例313设观测信号模型为x=A+w,其中被估计量A在[-A0,A0]上均匀分布,测量噪声w∝N(0,σ2),求A的最大后验概率估计和最小均方误差估计。 解先求最大后验概率估计,因为只有一个观测值,得到条件概率为 p(x|A)=12πσexp-(x-A)22σ2 p(A)=12A0,-A0≤A≤A0 0,A>A0(3189) 根据贝叶斯公式 p(A|x)=p(x|A)p(A)p(x)(3190) 这时难以得到p(x)的表达式,用式(3191)进行最大后验概率估计 Amap=argmaxA[p(x|A)p(A)](3191) 当-A0≤x≤A0时,p(x|A)p(A)的最大值出现在A=x处,所以,A^map=x; 当x>A0时,p(x|A)p(A)的最大值出现在A=A0处,A^map=A0; 当x<-A0时,p(x|A)p(A)的最大值出现在A=-A0处,A^map=-A0,即 A^map=-A0,x<-A0 x,-A0≤x≤A0 A0,x>A0(3192) 易证最大似然估计的结果为A^ml=x。A满足在[-A0,A0]上均匀分布的假设下,显然最大后验概率估计的结果更合理。图33给出了不同噪声强度下MAP和ML估计的均方误差性能比较,可见在高噪声强度下,MAP估计性能明显优于ML估计。在实际应用中,如果能获得被估计参数的先验概率信息,可改善估计结果。 图33不同噪声强度下MAP和ML估计性能比较 下面求最小均方误差估计,由式(3176),得 A^ms= ∫∞-∞Ap(A|x)dA = ∫∞-∞Ap(x|A)p(A)p(x)dA =∫∞-∞Ap(x|A)p(A)dA∫∞-∞p(x|A)p(A)dA =∫A0-A0A2πσexp-(x-A)22σ2·12A0dA∫A0-A012πσexp-(x-A)22σ2·12A0dA =∫x+A0x-A0(x-u)·exp-u22σ2du∫x+A0x-A0exp-u22σ2du =x-2σ2∫(z+a)/2(z-a)/2u·exp(-u2)duσ∫z+az-aexp(-u2/2)du =x-σ{exp[-(z-a)2/2]-exp[-(z+a)2/2]}2π[Q(z-a)-Q(z+a)] (3193) 式中,a=A0/σ,z=x/σ表示归一化观测值,Q(·)为标准正态概率密度函数的概率右尾函数,如式(3194)所示。 Q(v)=12π∫+∞vexp(-x2)dx(3194) 例314设信号有N次独立观测模型xn=A+wn (n=0,1,…,N-1),其中wn∝N(0,σ2),A∝N(μA,σ2A),求A的最大后验概率估计和最小均方误差估计。 解先求后验概率密度 p(A|x)=p(x|A)p(A)∫∞-∞p(x|A)p(A)dA =1(2πσ2)N2exp-12σ2∑N-1n=0(xn-A)212πσ2Aexp-12σ2A(A-μA)2∫∞-∞1(2πσ2)N2exp-12σ2∑N-1n=0(xn-A)212πσ2Aexp-12σ2A(A-μA)2dA =exp-121σ2(NA2-2NAx-)+1σ2A(A-μA)2∫∞-∞exp-121σ2(NA2-2NAx-)+1σ2A(A-μA)2dA =exp-12W(A)∫∞-∞exp-12W(A)dA(3195) 式中,x-=1N∑N-1n=0xn为样本均值,且 W(A)=1σ2(NA2-2NAx-)+1σ2A(A-μA)2(3196) 注意,式(3195)的分母与A无关,W(A)是A的二次型,经过配方可以把式(3196)写成 W(A)=1σ2A|x(A-μA|x)2-μ2A|xσ2A|x+μ2Aσ2A(3197) 式中, σ2A|x=Nσ2+1σ2A-1(3198) μA|x=Nσ2x-+μAσ2Aσ2A|x(3199) 将式(3197)代入式(3195),得 p(A|x)=exp-12σ2A|x(A-μA|x)2exp-12μ2Aσ2A-μ2A|xσ2A|x∫∞-∞exp-12σ2A|xA-μA|x2exp-12μ2Aσ2A-μ2A|xσ2A|xdA =12πσ2A|xexp-12σ2A|x(A-μA|x)2(3200) 由式(3200)可以看出,后验概率密度分布满足高斯分布的。由于最小均方误差估计为被估计量的条件均值,所以 A^ms=μA|x=Nσ2x-+μAσ2Aσ2A|x=Nσ2x-+μAσ2ANσ2+1σ2A(3201) 另外,由于最大后验概率估计是使后验概率最大对应的A值,因此,由式(3183)可得 A^map=μA|x= A^ms(3202) 3.6线性最小均方误差估计 对于随机参量的估计,在3.5.2节中介绍了最小均方误差估计。最小均方误差估计是被估计量的条件均值,这个条件均值通常都是观测的非线性函数,估计器实现起来比较复杂。条件均值的计算需要用到被估计量θ的概率密度p(θ),如果并不知道概率密度p(θ),而只知道θ的一、二阶矩特性,并且希望估计器能用线性系统实现,这时可以采用线性最小均方误差估计。 3.6.1随机参量的线性最小均方误差估计 线性最小均方误差(Linear Minimum Mean Square Error, LMMSE)估计是一种使均方误差最小的线性估计。假定观测为{xn},n=0,1,…,N-1,那么估计结果定义为观测的线性组合 θ^=∑N-1n=0anxn+b(3203) 考虑实信号,估计的均方误差为 MSE(θ^)=E[(θ-θ^)2]=Eθ-∑N-1n=0anxn-b2(3204) 选择一组最佳系数an和b,使式(3204)的均方误差达到最小。注意,引入系数b是因为观测x和被估计量θ的均值可能不为零。如果观测x和被估计量θ的均值都为零,那么系数b就可以省略。 均方误差对系数b求导,并令导数等于零,得 MSE(θ^)b=-2Eθ-∑N-1n=0anxn-b=0(3205) 经整理后得 b=E(θ)-∑N-1n=0anE(xn)(3206) 将式(3206)代入式(3204),得 MSE(θ^)=E∑N-1n=0an[xn-E(xn)]-[θ-E(θ)]2(3207) 令a=[a0,a1,…,aN-1]T,x=[x0,x1,…,xN-1]T,那么,式(3203)可表示为 θ^=aTx+b(3208) 而式(3206)可改写为 b=E(θ)-aTE(x)(3209) 估计的均方误差可表示为 MSE(θ^)=E{[aT(x-E(x))-(θ-E(θ))]2} =E[aT(x-E(x))(x-E(x))Ta]-E[aT(x-E(x))(θ-E(θ))]- E[(θ-E(θ))(x-E(x))Ta]+E[(θ-E(θ))2] =aTCxa-aTcθx-cTθxa+Cθ(3210) 其中,Cx=E{[x-E(x)][x-E(x)]T}为观测的协方差矩阵,cθx=E{[θ-E(θ)][x- E(x)]}是待估计量θ与观测x的协方差向量,且Cθ是θ的方差。 将估计的均方误差对a求导,得 MSE(θ^)a=2Cxa-2cθx(3211) 令导数等于零,可解得 a=C-1xcθx(3212) 将式(3212)和式(3209)代入式(3208),得 θ^lmmse=cTθxC-1xx+E(θ)-cTθxC-1xE(x)(3213) 式(3213)经整理后得到线性最小均方误差估计为 θ^lmmse=E(θ)+cTθxC-1x[x-E(x)](3214) 如果θ和x的均值为零,则 θ^lmmse=cTθxC-1xx(3215) 将式(3215)代入式(3210)中,可以得到最小均方误差的表达式 MSE(θ^lmmse)=Cθ-cTθxC-1xcθx(3216) 线性最小均方误差能得到简洁的估计子,并且可定量分析估计误差,因此应用广泛。 例315设有N次独立观测信号 xn=A+wn,n=0,1,…,N-1(3217) 其中,A在(-A0,A0)上服从均匀分布,wn是零均值高斯白噪声,方差为σ2,且A与wn统计独立,求A的线性最小均方误差估计。 解可以把观测信号模型写成向量形式 x=AeN×1+w(3218) 其中,eN×1是N维的全1向量。根据题意,E(A)=0,σ2A=E(A2)=(2A0)2/12,E(x)=0且 Cx=E(xxT)=E[(Ae+w)(Ae+w)T]=E(A2)eeT+σ2I=σ2AeeT+σ2I(3219) cAx=E(Ax)=E[A(Ae+w)]=E(A2)e=σ2Ae(3220) 其中,I为单位矩阵,由式(3215),得 A^lmmse=cTAxC-1xx=σ2AeT[eeTσ2A+σ2I]-1x=σ2Aσ2eTeeTσ2Aσ2+I-1x(3221) 根据矩阵求逆中常用的ShermanMorrison公式 (B+cdT)-1=B-1-B-1cdTB-11+dTB-1c(3222) 式(3221)经整理后,得 A^lmmse=σ2Aσ2A+σ2/Nx-=A20/3A20/3+σ2/Nx-(3223) 3.6.2线性最小均方误差估计的几何解释 可以用向量空间的概念解释线性最小均方误差估计。在以下的讨论中假定被估计量θ和观测x都是零均值的,如果不是零均值,则总可以定义零均值的随机变量θ′=θ-E(θ),x′=x-E(x)。 考虑一个由所有随机变量构成的线性空间(见图34),θ和x0,x1,…,xN-1可以看作线性空间的点(每个点由一个向量描述)。定义两个向量的内积为(x,y)=E(xy),向量的长度为随机变量方差的均方根,即‖x‖=E(x)2,如果(x,y)=E(xy)=0,称这两个向量是正交的。由于随机变量是零均值的,因此两个向量正交和随机变量不相关等价。线性最小均方误差估计的目的是用x0,x1,…,xN-1的线性组合逼近θ。先考虑只有一个样本的情况,即用θ^=a0x0逼近θ,易见最优的估计就是θ到x0的投影(如图34(a)所示),这时误差为θ到x0的垂直距离,从几何的角度易知这时误差向量的长度最短。当使用两个样本x0,x1时,θ^=a0x0+a1x1,注意到x0,x1的线性组合张成一个平面,这时最优估计就是θ到该平面的投影(如图34(b)所示),同时误差向量ε正交于由{x0,x1}所张成的平面,易证这时误差向量的长度最短。 图34线性最小均方误差估计的几何解释 一般情况下,线性估计子用x0,x1,…,xN-1的线性组合近似θ,相当于在x0,x1,…,xN-1所张成的子空间上找一个点,使其与θ的距离最小。根据误差定义 MSE(θ^)=E [(θ-θ^)2]=Eθ-∑N-1i=0aixi2=‖θ-∑N-1i=0aixi‖2(3224) 可见最小均方误差等价于误差向量ε=θ-θ^长度的平方最小。根据上面的论述,这时误差向量将正交于x0,x1,…,xN-1所张成的子空间,即 ε⊥span{x0,x1,…,xN-1}(3225) 根据正交的定义,可以得出 E[(θ-θ^)xn]=0,n=0,1,…,N-1(3226) 式(3226)是线性最小均方误差估计非常重要的正交原理,也可由式(3224)直接对加权系数求解得到。由正交原理可知,线性最小均方误差估计可以通过使估计误差和每个观测数据正交得到。利用正交原理,加权系数很容易求出 Eθ-∑N-1n=0anxnxj=0,j=0,1,…,N-1(3227) 或者写成 ∑N-1n=0anE(xnxj)=E(θxj),j=0,1,…,N-1(3228) 用矩阵表示为 E(x20)E(x0x1)…E(x0xN-1) E(x1x0)E(x21)…E(x1xN-1) ︙︙︙ E(xN-1x0)E(xN-1x1)…E(x2N-1)a0 a1 ︙ aN-1=E(θx0) E(θx1) ︙ E(θxN-1)(3229) 或者写成 Cxa=cθx(3230) 所以 a=C-1xcθx(3231) θ的线性最小均方误差估计为 θ^lmmse=aTx=cTθxC-1xx(3232) 式(3232)与式(3215)一致。 下面考虑随机向量的线性最小均方误差估计。设待估计量为θ=[θ0,θ1,…,θp-1]T,观测向量x=[x0,x1,…,xN-1]T,线性估计可表示为 θ^=Ax+b(3233) 式中,b是p×1维的向量,A是p×N维的矩阵,所有估计的均方误差和可表示为 MSE(θ^)=E∑p-1i=0(θi-θ^i)2=E[(θ-θ^)T(θ-θ^)](3234) 将式(3233)代入式(3234),得 MSE(θ^)=E[(θ-Ax-b)T(θ-Ax-b)](3235) 分别对A和b求导,并令导数等于零,得 MSE(θ^)b=-2E[θ-Ax-b]=0(3236) MSE(θ^)A=-2E{[θ-Ax-b]xT}=0(3237) 式(3237)中对矩阵的求导与式(372)相同,定义为对矩阵每个元素求偏导数,结果为矩阵,且大小与A相同。由式(372)、式(3237),可得 b=E(θ)-AE(x)(3238) A=CθxC-1x(3239) 式中, Cθx=Cov(θ,x)=E{[θ-E(θ)][x-E(x)]T}(3240) Cx=Cov(x)=E{[x-E(x)][x-E(x)]T}(3241) 由此可得线性最小均方误差估计为 θ^lmmse=E(θ)+CθxC-1x[x-E(x)](3242) 式(3242)在表达形式上除θ是向量外,其他与式(3214)相同。线性最小均方误差估计的性质总结如下(这里不做具体证明): (1) 线性最小均方误差估计是无偏估计。 (2) 线性最小均方误差估计的均方误差阵要小于任何其他线性估计的均方误差阵,即 E[(θ-θ^lmmse)(θ-θ^lmmse)T]≤E[(θ-θ^)(θ-θ^)T](3243) 式中,θ^是任意的线性估计。式(3243)中符号≤作用于矩阵时,含义如下: 如果B≤A,则A-B是半正定矩阵。式(3243)证明从略。 (3) 对于线性变换α=Hθ+b,线性最小均方误差估计具有不变性,即 α^lmmse=Hθ^lmmse+b(3244) (4) 线性最小均方误差估计具有叠加性,即如果α=θ1+θ2,则 α^lmmse=θ^1,lmmse+θ^2,lmmse(3245) 3.7最小二乘估计 在前面介绍的几种估计方法中,最小均方误差估计、最大后验概率估计需要知道被估计量的先验概率密度,最大似然估计需要知道似然函数,线性最小均方误差估计需要知道被估计量的一、二阶矩,如果这些概率密度或矩未知,就不能采用这些方法,这时可以采用最小二乘(Least Square,LS)估计。最小二乘估计对统计特性没有做任何假定,因此,它的应用非常广泛。 在线性最小均方误差估计中,选择均方误差作为估计量性能好坏标准的度量,也就意味着估计与真实参数之间的差别平均来说达到最小。而在最小二乘估计中,则是使给定的观测数据与信号或者无噪声的观测数据之差的平方和达到最小。 采用例310的线性观测模型,待估计量为θ=[θ0,θ1,…,θp-1]T,观测信号用向量和矩阵可表示为 x=Hθ+w(3246) 式中, x=[x0,x1,…,xN-1]T,w=[w0,w1,…,wN-1]T H=h00h01…h0(p-1) h10h11…h1(p-1) ︙︙︙ h(N-1)0h(N-1)1…h(N-1)(p-1) 估计偏差的平方和可表示为 J(θ^)=[x-Hθ^]T[x-Hθ^]=∑N-1n=0xn-∑p-1j=0hnjθ^j2(3247) 在没有任何先验信息情况下,可以认为最优估计就是使得上述的平方和误差最小。当待估计参数的重要性存在差别时,可以对每个参数的误差进行加权,观测与估计偏差的加权平方和可表示为 Jw(θ^)=[x-Hθ^]TW[x-Hθ^]=∑N-1j=0∑N-1n=0xn-∑p-1k=0hnkθ^kWnjxn-∑p-1k=0hjkθ^k(3248) 其中W为加权矩阵。最小二乘估计就是使J(θ^)最小的估计,记为θ^ls; 加权最小二乘(Weighted Least Square)估计则是使Jw(θ^)最小的估计,记为θ^lsw。 接下来考虑如何求解,将式(3247)重写为 J(θ^)=xTx-xTHθ^-θ^THTx+θ^THTHθ^(3249) 根据式(373)和式(374)的向量函数求导公式,得到 xTxθ^=0(3250) xTHθ^θ^=HTx(3251) θ^THTxθ^=HTx(3252) θ^THTHθ^θ^=2HTHθ^(3253) 求J(θ^)对θ^的导数,并令导数等于零,得 J(θ^)θ^=-2HTx+2HTHθ^=0(3254) 由此可解得最小二乘估计为 θ^ls=(HTH)-1HTx(3255) 求Jw(θ^)对θ^的导数,并令导数等于零,得 Jw(θ^)θ^=-2HTW[x-Hθ^]=0(3256) 由此可解得加权最小二乘估计为 θ^lsw=(HTWH)-1HTWx(3257) 最小二乘估计具有如下特点: (1) 对于线性的观测模型,最小二乘估计和加权最小二乘估计都是线性估计,对测量噪声的统计没有做任何假定,应用十分广泛。 (2) 当测量噪声的均值为零时,即E(wi)=0时,最小二乘估计和加权最小二乘估计都是无偏估计。 (3) 对于加权最小二乘估计,如果有一些模型的知识,如E(w)=0,E[wwT]=R,当W=R-1时,估计误差的方差达到最小,这时加权最小二乘估计和式(3151)的最大似然估计是一致的。 3.8信号检测基础 信号检测的目的是在干扰、噪声背景下判断目标信号是否存在。当存在多种类型目标信号时,还需要判断信号类型。根据目标信号模型的不同,信号检测分确定性信号检测和随机信号检测进行讨论,本节只讨论实信号。 3.8.1确定性信号检测 这里以确定性信号为例介绍信号检测的基本概念和方法。假设目标信号sn的波形已知,得到了接收信号xn的N个样本,在此基础上判决其中是否包含目标信号,可定义如下两种事件(或称假设)。 H0:xn=wn,n=0,1,…,N-1 H1:xn=sn+wn,n=0,1,…,N-1(3258) 接下来考虑如何对上述事件进行检验,其过程可以看成: 将观测空间划分成不同判决区域,观测样本x=[x0,x1,…,xN-1]落在哪个区域就判决成哪个相应事件。一般可基于事件的概率来划分判决区域,定义两种事件的后验概率分别为p(H1|x),p(H0|x),后验概率给出了在得到观测样本后,对H0,H1事件发生概率的认知。一种直观的判决准则是,在观测样本的基础上,哪个事件发生的概率大就认为发生了该事件,即当p(H1|x)>p(H0|x)时,认为H1事件(假设)为真; 反之,若p(H1|x)<p(H0|x)时,则认为H0为真。若p(H1|x)=p(H0|x),则可随机认定某一事件为真。于是,可以把H1事件的判决区域定义为所有使得p(H1|x)>p(H0|x)的观测空间。具体的数学描述如式(3259)所示。 H1:p(H1|x)p(H0|x)>1(3259) 式(3259)中的后验概率一般难以直接计算。基于贝叶斯公式有 p(H1|x)p(x)=p(x|H1)p(H1)(3260) 式(3259)等价于 H1:L(x)=p(x|H1)p(x|H0)≥γ(3261) 其中,p(x|H0),p(x|H1)是3.4节最大似然估计中提到的似然函数,L(x)称为似然比(Likelihood Ratio,LR),有时为了计算便利,使用其对数形式ln[L(x)],称为对数似然比(Logarithm Likelihood Ratio,LLR)。γ=p(H0)/p(H1)表示判决门限。 计算似然函数需要知道xn的概率密度,实际中可通过假设噪声符合某种概率分布实现,这里以最常见的假设: wn是零均值、方差为σ2的高斯白噪声为例,分析如何具体应用式(3261)进行检验。根据假设可得到如下似然函数: p(x|H1)=1(2πσ2)N2exp-12σ2∑N-1n=0(xn-sn)2(3262) p(x|H0)=1(2πσ2)N2exp-12σ2∑N-1n=0x2n(3263) 利用式(3263)可计算如下似然比和对数似然比: L(x)=exp-12σ2∑N-1n=0(xn-sn)2-x2n(3264) l(x)=lnL(x)=-12σ2∑N-1n=0(xn-sn)2-x2n(3265) 采用与式(3261)等价的检验准则: lnL(x)>lnγ,经简单计算,可得式(3261)等价于对式(3266)中T(x)的检验。 H1:T(x)=∑N-1n=0xnsn=xTs>γ-(3266) 其中 γ-=σ2lnγ+12∑N-1n=0s2n(3267) 注意到T(x)可以看作匹配滤波器的输出。确定性信号的匹配滤波器定义如下: h(n)=s(N-1-n),n=0,1,…,N-1 0,其他(3268) 匹配滤波器的输出为y(n)=x(n)h(n),并且可验证T(x)=y(N-1)。 在已知噪声概率特性情况下,上述分析均可用,下面将上述推导拓展到高斯有色噪声下的信号检测问题。记x=[x0,x1,…,xN-1]T,s=[s0,s1,…,sN-1]T,w=[w0,w1,…,wN-1]T,检测问题可描述为 H0:x=w,w∝N(0,C) H1:x=s+w(3269) 判决准则重写为 H1:L(x)=p(x|H1)p(x|H0)>γ(3270) 可得到似然函数、似然比及对数似然比分别为 p(x|H1)=1(2π)N2det12(C)exp-12(x-s)TC-1(x-s)(3271) p(x|H0)=1(2π)N2det12(C)exp-12xTC-1x(3272) l(x)=lnL(x)=-12(x-s)TC-1(x-s)-xTC-1x=xTC-1s-12sTC-1s(3273) 利用lnL(x)>lnγ准则进行判决,得到 H1:T(x)=xTC-1s>γ+12sTC-1s=γ-(3274) 同理,可以由C-1s构造滤波器,而T(x)是滤波器某个时刻的输出,这种滤波器称为广义匹配滤波器。 T(x)>γ-定义了观测空间的一种划分规则,对所有满足T(x)>γ-的观测x计算积分可以得到H1成功检测的概率。根据事件的不同,可定义4种不同概率,其定义和计算方式如表31所示。 表31二元信号检测中的事件和概率计算 事件定义标号 H1为真,判H1成立H1检测概率P(H1|H1)=∫T(x)>γ-p(x|H1)dx H0为真,判H0成立H0检测概率1-P(H1|H0) H0为真,判H1成立虚警概率P(H1|H0)=∫T(x)>γ-p(x|H0)dx H1为真,判H0成立漏检概率1-P(H1|H1) 在实际应用中表31中不同的事件会有不同的代价,例如在雷达信号处理中,漏检造成的代价远大于虚警 造成的代价。考虑更一般化的情况,为表中4种不同事件定义各自的代价因子,并定义如下平均代价: C=[C00P(H0|H0)+C01P(H1|H0)]P(H0)+ [C10P(H0|H1)+C11P(H1|H1)]P(H1) (3275) 其中,C00,C01,C10,C11分别为各事件的代价因子,经简单计算可得 C=[C00(1-P(H1|H0))+C01P(H1|H0)]P(H0)+ [C10(1-P(H1|H1))+C11P(H1|H1)]P(H1) =C00P(H0)+(C01-C00)P(H0)P(H1|H0)+ C10P(H1)-(C10-C11)P(H1)P(H1|H1) =C00P(H0)+C10P(H1)+ ∫x∈Z1(C01-C00)P(H0)P(x|H0)-(C10-C11)P(H1)P(x|H1)dx (3276) 其中,Z1表示判决H1成立的判决区域。可见要使代价最小,在判决区域内的积分函数取值为负,即在判决区域内满足 (C01-C00)P(H0)P(x|H0)<(C10-C11)P(H1)P(x|H1)(3277) P(x|H1)P(x|H0)>(C01-C00)P(H0)(C10-C11)P(H1)(3278) 式(3278)称为信号检测的贝叶斯准则。从式(3278)可以看到,如果认为正确的检测不产生额外代价,即式(3275)中C00=C11=0,同时虚警和漏检有相同的代价,即式(3275)中C01=C10=1,则式(3278)的检测准则等价于式(3261)。 以上的贝叶斯准则需要确定代价因子和先验概率,实际中往往很难确定,特别是先验概率的确定。实际中往往对虚警概率和漏检概率有明确的要求,可以从虚警概率和漏检概率的约束考虑信号检测。一种常见的准则是在虚警概率恒定的约束下最小化漏检概率,称为纽曼皮儿森准则(NeymanPearson,NP)。 令PF、PM分别表示虚警概率和漏检概率,α表示允许的虚警概率上限,Z0,Z1分别表示信号检测的判决区域,问题模型可表示为 minZ0,Z1PM PF=α(3279) 即在满足虚警概率要求的条件下,寻找观测空间的划分准则,使得漏检概率尽量小。可利用拉格朗日乘子法求解上述问题,构造如下目标函数 J=PM+λ(PF-α)(3280) 其中,λ为拉格朗日乘子。目的是要确定一种观测空间的划分,使J最小,将虚警概率和漏检概率的计算表达式代入,得到 J=∫Z0p(x|H1)dx+λ∫Z1p(x|H0)dx-α =λ(1-α)+∫Z0[p(x|H1)-λp(x|H0)]dx(3281) 式(3281)中第一项非负,要使J尽可能小,则式(3281)积分结果尽可能小。一种直观的做法是将所有使p(x|H1)-λp(x|H0)<0的观测空间归入Z0,或将使p(x|H1)>λp(x|H0)成立的观测空间定义为Z1,从而得到如下H1的检测准则: H1:p(x|H1)p(x|H0)>λ(3282) 式(3282)中λ为未知量,不能直接应用。λ可根据约束条件确定,但难以得到闭式解。一般可采用迭代方式求解,即从λ的一个初始值开始,由式(3282)得到观测空间的分割Z0、Z1,再根据如下约束条件进行验证。如果满足,则停止迭代; 否则,调整λ取值重新验证。 ∫Z1p(x|H0)dx=PF(3283) 3.8.2随机信号检测 随机信号检测所采用的准则与确定性信号检测一致。但由于对随机信号只能知道其统计特征,因此具体计算有所不同。下面看一个简单例子,考虑以下两种假设。 H0:x=w,w∝N(0,σ2I) H1:x=s+w(3284) 其中,sn是平稳随机信号。先考虑一个简单情况,sn是高斯信号,均值为零,方差为σ2s,以上两种假设可描述为 H0:x∝N(0,σ2I) H1:x∝N(0,(σ2s+σ2)I)(3285) 采用似然比及对数似然比进行判断,得到如下H1的判决准则: H1:L(x)=p(x|H1)p(x|H0)>γ(3286) 得到似然函数及似然比如下: P(x|H1)=1(2π(σ2s+σ2))N2exp-12(σ2s+σ2)∑N-1n=0x2n(3287) P(x|H0)=1(2πσ2)N2exp-12σ2∑N-1n=0x2n(3288) l(x)=lnL(x)=N2lnσ2σ2s+σ2+12σ2sσ2(σ2s+σ2)∑N-1n=0x2n(3289) 采用l(x)>lnγ判决准则,对上式略加整理,可得到如下等效判决准则: H1:T(x)=∑N-1n=0x2n=xTx>lnγ-N2lnσ2σ2s+σ22σ2(σ2s+σ2)σ2s=γ-(3290) 可以看到H1的判决依赖于采样信号的能量。当能量大于门限值时,判断存在目标信号,因此上述检测方法又称为能量检测器。能量检测器可推广到更一般的情况,即检测信号具有任意的协方差矩阵Cs,考虑以下两种假设 H0:x∝N(0,σ2I) H1:x∝N(0,Cs+σ2I)(3291) 得到似然函数及似然比如下: P(x|H1)=1(2π)N2det12(Cs+σ2I)exp-12xT(Cs+σ2I)-1x(3292) P(x|H0)=1(2πσ2)N2exp-12σ2xTx(3293) l(x)=lnL(x)=12ln(σ2)Ndet(Cs+σ2I)-12xT[(Cs+σ2I)-1-1σ2I]x(3294) 可得到如下等价于l(x)>lnγ的判决准则: H1:T(x)=σ2xT1σ2I-(Cs+σ2I)-1x>γ-(3295) 利用矩阵求逆定理 (A+BCD)-1=A-1-A-1B(DA-1B+C-1)-1DA-1(3296) 在式(3296)中,令A=σ2I,B=D=I,C=Cs,有 (σ2I+Cs)-1=1σ2I-1σ41σ2I+C-1s-1(3297) 代入式(3295),得到如下判决准则: H1: T(x)=xT 1σ2C-1s+1σ2I-1x=xT(σ2C-1s+I)-1x =xT[(σ2I+Cs)C-1s]-1x=xTCs(σ2I+Cs)-1x>γ- (3298) 和式(3290)相比,式(3298)相当于在能量检测器基础上引入加权矩阵W=Cs(σ2I+Cs)-1,得到加权检测量T(x)=xTWx,称为广义能量检测器。注意到式(3298)可写成 H1:T(x)=xTs^>γ- s^=Cs(σ2I+Cs)-1x(3299) 信号检测的过程可以看成先得到目标信号的一个估计,然后再利用估计得到的信号构造一个如式(3266)所示的匹配滤波器。注意到上式中的s^实际是目标信号的线性最小均方误差(LMMSE)估计。说明如下: 根据3.6节的结论,基于观测x,目标参数θ的线性MMSE估计为 θ^=CθxC-1xx(3300) 将式(3300)用于从x中估计s,可得 Csx=E(sxT)=E(ssT+swT)=E(ssT)=Cs(3301) Cx=E(xxT)=E(ssT+swT+wsT+wwT) =E(ssT)+E(wwT) =Cs+σ2I(3302) 式(3302)中假设目标信号和噪声不相关,这在实际中一般可以满足。 以上讨论都基于假设在H0和H1条件下的似然函数是可知的,但实际情况经常没这么理想。一个最简单的例子是描述噪声特性的参数,如方差,经常是未知的。目标信号的参数也有可能是未知的,例如在无源雷达和声呐系统中,接收机对来 的波信号的频率是未知的。估计存在未知参数时,一般 分两个步骤处理,首先 是解决H0和H1条件下的未知参数估计问题,例如用最大似然准则 θ^0= maxθ0p(x|θ0,H0)θ^1= maxθ1p(x|θ1,H1)(3303) 然后将估计结果用于似然比准则进行判决 H1:L(x)=p(x|θ^1,H1)p(x|θ^0,H0)>γ(3304) 这种方法称为广义似然比准则。如果先验概率分布已知,则可以使用贝叶斯准则计算似然比 H1:L(x)=p(x|H1)p(x|H0)=∫p(x|θ1,H1)dθ1∫p(x|θ0,H0)dθ0>γ(3305) 下面以一个例子进一步说明,基于式(3290)且假设目标信号具有非零均值 H0:x∝N(0,σ2I) H1:x∝N(ue,Cs+σ2I)(3306) 其中,e是N×1全1向量,得到似然函数如下: P(x|H1)=1(2π)N2det12(Cs+σ2I)exp-12(x-ue)T(Cs+σ2I)-1(x-ue)(3307) P(x|H0)=1(2πσ2)N2exp-12σ2xTx(3308) u的最大似然估计为 u^=1NxT(Cs+σ2I)-1e(3309) 和式(3295)的推导相似,可得到以下判决准则 T(x)=-(x- u^e)T(Cs+σ2I)-1(x- u^e)+1σ2xTx>γ- xT[(Cs+σ2I)-1-1σ2I]x+Nu^2>γ- xTCs(σ2I+Cs)-1x+Nu^2>γ-(3310) 和式(3295)相比,式(3310)的不同之处在于引入了非零均值的影响。 上述讨论可推广到多种信号的检测,考虑如下模型: H0:x=w,w~N(0,C) Hi:x=si+w,i=1,2,…,K(3311) 其中,K表示有多少种不同的信号。为了实现最低的检测错误概率,可采用最大后验概率检测准则 Hi:p(Hi|x)>p(Hj|x),j=0,1,…,K; j≠i(3312) 假设各种假设具有相等概率,可采用最大似然检测准则 Hi:p(x|Hi)>p(x|Hj),j=0,1,…,K,j≠i(3313) 在实际应用中经常先计算各种假设对H0的似然比 Li(x)=p(x|Hi)p(x|H0),i=1,2,…,K(3314) 如果所有Li(x)的取值均小于1,则判决H0为真; 否则,判决具有最大似然比的假设为真。特别要指出的是,上述方法只能应用于各种假设具有相等概率的情况,如果概率不相等,则必须直接基于式(3312)计算。 本章习题 1. 假定观测x0,x1,…,xN-1是独立同分布的,且在(0,θ)上服从均匀分布,求θ的最大似然估计。 2. 证明式(339) I(θ)=Elnp(x|θ)θ2=-E2lnp(x|θ)θ2 3. 设N次观测为xi=A+wi,i=0,1,…,N-1,其中,A为未知的确定信号,噪声wi (i=0,1,…,N-1) 相互独立并服从相同的分布N(0,σ2),噪声的方差已知。 (1) 求A的最大似然估计。 (2) A的最大似然估计是否为无偏估计? (3) A的最大似然估计是否为有效估计?估计的方差等于多少? 4. 在习题3中,如果A和σ2都是未知的,试分别求其最大似然估计,并讨论它们的无偏性。 5. 设观测模型为xi=A+wi,其中wi∝N(0,σ2),如果A为未知常量,求A的最大后验概率估计和最小均方误差估计。 6. 设有N次独立同分布的观测x0,x1,…,xN-1,且单个样本在(0,θ)内服从均匀分布,证明: Elnp(x|θ)θ≠0,θ>0 7. 平坦信道的信道估计问题可以用下式描述 xi=Asi+wi,i=0,1,…,N-1 其中,wi是零均值高斯白噪声,方差为σ2,si是已知的训练序列,xi是接收信号,求估计A的CRB。证明有效估计量存在,并求它的方差。当N→∞时,估计的方差会怎样变化? 8. 设有两次观测 x0=As0+w0 x1=As1+w1 其中,w=[w0,w1]T是零均值高斯随机向量,协方差矩阵为 C=σ21ρ ρ1 求估计A的CRB,并将它与w为高斯白噪声的情况(ρ=0)进行比较,解释当ρ=±1时会怎样? 9. 考虑高斯噪声下的单频频率估计问题 xi=Aejωn+φ+wi,i=0,1,…,N-1 其中,wi是复高斯白噪声,方差为σ2。A,ω,分别为幅度、频率和相位,求: (1) A、为已知参数时频率估计的CRB; (2) A、ω、都是未知参数时频率估计的CRB。 10. 设观测的概率密度为p(x|θ),θ为未知常量,待估计量为α=g(θ),证明: 估计α的CRB为 Var(α)≥g(θ)θ2-E2lnp(x|θ)θ2 11. 假定α=g(θ)=Aθ+b,假设θ-是θ的有效估计,即Var(θ-)=I-1(θ),证明: α=Aθ+b也是有效估计量,即对于线性变换,估计量的有效性得以保持。 12. 考虑一个正弦参数的估计问题,假定观测为 xi=Acos(2πf0i+φ)+wi,i=0,1,…,N-1 其中,wi为零均值高斯白噪声,方差为σ2,A、f0、为未知参量,且A>0,0<f0<1/2,设θ=[A,f0,]T,求其费希尔信息矩阵。 13. 设有N次独立的观测为 xi=A+wi,i=0,1,…,N-1 其中,A为未知常数,wi是拉普拉斯噪声,分布为 p(wi)=12e-wi 求A的最大似然估计子,并分析当N→∞时,最大似然估计的方差是否达到CRB。 14. 设观测信号模型为 xi=A+wi,i=0,1,…,N-1 式中,随机变量A为信号幅值,E[A]=0,E[A2]=a,A是已知参数。wi是零均值白噪声,与A相互独立,求A的线性最小均方误差估计。 15. 在平稳白噪声背景中,对信号参量做线性最小均方误差估计,两次观测数据为 xi=A+wi(i=0,1),待估计参数A为信号幅度,wi为零均值,方差为σ2n的白噪声,且信号幅值与噪声不相关,求A的最佳线性估计。 16. 设观测信号模型为 xi=A+wi,i=0,1,…,N-1 未知参数A具有如下的概率密度分布: f(A)=λexp(-λA),A>0 0,A<0 其中,λ>0,wi是零均值。方差为σ2n的高斯白噪声与A相互独立,试求A的最大后验概率估计。 17. 对于无偏估计子,最小均方误差等价于最小方差。对线性估计子 θ^=∑N-1i=0aixi 引入约束条件E[xi]=siθ,则最小方差无偏估计子参数如下: a=C-1ssHC-1s 其中,C是xi的协方差矩阵,证明所得估计子是无偏估计子并求其估计误差。 18. 设x(t)=AcosΩ0t+w(t),通过采样对幅度A做线性估计。在Ω0t=0,Ω0t=π/3处采样两次,记为x0、x1,并设E[A]=0,E[A2]=2,E[w0]=E[w1]=0,E[w20]=E[w21]=1,E[w0w1]=0,E[Aw0]=E[Aw1]=0,求: (1) 线性最小均方误差估计A^lmmse; (2) 线性最小均方误差估计的误差。 19. 对于x=Hθ+w的实系数线性观测模型,假定测量噪声的均值为零,证明最小二乘估计误差的均方误差阵为 Pθls=E(θ-θls)(θ-θls)T=(HTH)HTRH(HTH)-1 其中,R=E[wwT]。 20. 观测信号模型为yi=α+βxi+wi(i=0,1,…,N-1),通过采样对幅度a做线性估计,其中,xi已知,wi为未知干扰。求α,β的最小二乘估计。 21. 设有零均值实平稳随机过程x(t),t是[0,T]内的一点。若已知x(0)和x(T),利用x(0)和x(T)求x(t)的最佳线性估计。 22. 设有如下两种假设,观测次数为N次 H0,zk=nk H1,zk=2+nk(k=1,2,…,N) 其中,nk是满足零均值,方差为σ2n的正态分布,假设p(H0)=p(H1)=0.5,求: (1) 最小错误概率准则下的判决表达式; (2) 虚警概率和检测概率(结果由误差函数表示)。 23. 设有如下两种假设,观测次数为N次 H0,zk=nk H1,zk=1+nk(k=1,2,…,N) 其中,nk是满足零均值,方差为σ2n的正态分布。试构造一个虚警概率为0.1的NP (纽曼皮儿森) 检测器,求相应的检测概率。 24. 考虑如下四元数字通信系统: H0,zk=nk H1,zk=1+nk H2,zk=2+nk H3,zk=3+nk 其中,nk是满足零均值,方差为σ2n的正态分布,各个假设的先验概率相等。试设计一个四元信号的最佳检测系统,由每一个接收符号zk检测其发送符号。