现代谱估计法及应用效果

2009 年 11 月
第 44 卷
增刊 1
# 处理技术 #
现代谱估计法及应用效果
刘志刚*
¹
李录明 º
赵冬梅 »
( ¹ 东方地球物理公司研究院 , 河北涿州 072751; º成都理工大学信息学院 , 四川成都 610059; » 东方地球物理公司物探技术研究中心 , 河北涿州 072751)
刘志刚 , 李录明 , 赵冬梅 . 现代谱估计法及应用 效果 . 石油地球物理勘探 , 2009 , 44 ( 增刊 1) : 5~ 9 摘要 本文针对 Burg 谱估计法中存在的问题 , 讨论了改进 Burg 谱估计法和改进协方差谱 估计法 , 以理论 信号为 测试对象 , 对不同谱估计法的应用效果进行了对比 , 结果表明 : ¹ Burg 谱估计法分辨率明显高于 Welch 谱估计法 , 但 Burg 谱估计法存在明显的峰值偏移 , 改进 Bur g 谱估计法几乎没有峰值偏移 ; º改进协方差谱估计法和 Burg 谱 估计法都具有较高的分辨率 , 而前者的波峰较后者更明显 、 尖锐 , 对于短数据 、 信号频率差异较小的信号 , 前者具有 更好的分辩效果, 还能抑制谱线分裂和出现假谱峰等问题 ; » Itakur a 算法求得的反射系数大于或等于 Bur g 谱估计 法求出的反射系数 , 使得接收到的信号更接近于实际输入信号 , 因此可用 Itakura 算法 替代 Burg 谱估计法 。 实际 地震剖面去噪结果表明 , 以高分辨率谱估计方法为基础的信噪分离方法具有较好的去噪效果 。 关键词 现代谱估计法 L evinson 递推算法 改进 Burg 谱估计法 改进协方差谱估计法 分辨率 反射系数
f 0 b
典( 非参数 ) 谱估计方法, 其以相关和傅里叶变换为 基础 , 对于长数据记录较适用 , 但无法根本解决频率 分辨率低和谱估计稳定性的问题 , 特别是在数据记 录很短的情况下, 这一问题尤其突出。以随机过程 的参数模型为基础的现代参数法功率谱估计具有更 高的频率分辨率和更好的适应性 , 可实现信号检测 或信噪分离, 对语音、 声纳雷达、 电磁波及地震波等 信号处理具有重要意义, 并广泛应用于通信、 自动控 制、 地球物理等领域。在现代参数法功率谱估计方 法中 , 比较有效且实用的是 Burg 谱估计法, 该法避 免了计算相关, 对短数据具有更强的适应性, 从而弥 补了经典谱估计法的不足, 但其也有一些自身的缺 陷。本文针对 Bur g 谱 估计法中存在 的问题, 讨论 了改进 Burg 谱 估计法和改进协方 差谱估计法 , 理 论模型及对实际地震剖面的处理结果表明, 现代谱 估计法的分辨率更高、 适应性更强[ 1~ 5] 。
N- 1
Q ^ m = *
n= m N- 1 n= m
E ^e
f m - 1
* ( n) eb m - 1 ( n - 1) 2 N- 1 2
E
^e m- 1 ( n)
f
*
n= m
E
e m- 1 ( n - 1) ( 3g)
b
第 44 卷
增刊 1
刘志刚等 : 现代谱估计法及应用效果
LS LS E p
7
实验结果表明 , Burg 谱估计法求出的 ^ Q m 小 于或等
* 于^ Q m 时( 即| ^ Q m| [ | Q ^* m | ) , 即由式 ( 3g ) 求得的反射
Px ( f ) = 1+
系数大于 Burg 谱估计法求出的反 射系数, 从 而使 得接收到的信号更接近于实际输入信号 , 减少失真。 实际上 , 在谱估计性能方面 Burg 谱估计法与 It aku r a 算法差别很小 , 在语音信号处理中常用式 ( 3g) 估 计反射系数。
数给功率谱 带来的影响 , Burg 提出了一 种新的算
1
引言
周期图、 自相关法及其改进方法 ( Welch) 为经
法[ 6 ] , 其基本思想是 , 利用由线性预测器直接获得的 观测数据的前向和后向预测的总均方误差之和的最 小值来估计预测系数, 进而通过 L ev inso n 算法的递 推公式求出 AR 模型优化参 数。 Burg 谱估计法步 骤如下 : ( 1) 计算初始值 e ( n) = e0 ( n) = x ( n)
k+ 1
^ p ( k) e- j2Pf k Ea
p
2
( 2e)
( 3b)
综上所述, 谱估计法的优点在于: ¹ 频率分辨 率高 ; º所得的 AR 模型稳定 ; » 计算 效率高。然 而, Burg 谱估计法也有公认的不足。首先 , 在高信 噪比时, 其谱线呈现出分裂 ; 其次 , 对于高阶模型来 讲, 该方法也可能引入假谱峰 ; 第三 , 对于噪声中的 正弦信号 , Burg 谱估计法对正弦信号的初始相位呈 现出敏感性, 特别是对短的数据记录更是如此。
k= 1
^ p ( k) e- j2Pkf Ea
p
( 4e)
考虑到以 r x ( l , k ) 为元素的相关矩阵不是 T o eplit z 阵, 因此不能利用 L ev inson 算法求解。但相 关矩阵的结构足以设计复杂性与 p 成正比的有效 计算方法。 Marple 设计了这样的算法 , 该算法具 有格型结构, 利用了 L evinso n 类型的阶递推和另外 的时间递推。这种形式的改进协方差 ( 无约束最小 二乘法 ) 方法也一直被称为数据未加窗最小二乘方 法。改进协方差方法基本上克服了谱线分裂、 谱峰 偏移和出现伪峰等缺点 , 但不能保证 AR 模型稳定, 并且所需运算量也偏大 , 然而在谱估计中这不是一 个大问题。
E
2 [ efm- 1 ( n) ] 2 + [ eb m- 1 ( n) ]
( 2a)
m = 1 , 2, ,, p
* 河北省涿州市东方地球物理公司研究院海外业务部 , 072751 本文于 2009 年 3 月 10 日收到。
6
石油地球物理勘探
2009 年
a p ( k) = ap- 1 ( k ) + Q p a p- 1 ( p - k) k = 1, 2, ,, p - 1 ap ( p ) = Q p
2
e ( n) = e e ( n) = e
b p
f p
f p- 1 b p- 1
( n) + Q pe
[ x ( n - 2) + a2 ( 1 ) x ( n - 1) + a 2 ( 2) x ( n) ] ( 2d) 同时采用两个使得 E 2 最小的必要条件, 即 5E 2 5E2 = 0, = 0 5 a2 ( 1 ) 5 a2 ( 2) 经推导可得关于 a2 ( 1) , a2 ( 2 ) 的方程组 c( 1 ) + b( 1) a 2 ( 1) + c( 1 ) a2 ( 2) = 0 c( 2 ) + c( 1 ) a2 ( 1) + b( 2 ) a2 ( 2) = 0 其中
n = 0, 1 , ,, N - 1 ( 1a)
此时均方差等于预测器平均输出功率 , 故有
2 1 R w 0 = R ( 0) = N N- 1 n= 0
Ex
2
( n) ( 1b)
R = E [ e ( 0) ]
2 w0 2
其中: e0f ( n) 、 e0b ( n) 分别为前、 后向 预测误差 初值; x ( n) 为观测数据 ; R w 0 为白噪声方差; E[ e ( 0) ] 为最
p , 而其余的预测系数 中, 一旦求出 Q p , 则 ap ( p ) = Q 就由低一阶模型的系数 ap - 1 ( k) 来确定 , 不能灵活
5
现代谱估计法应用结果
5. 1 不同谱估计方法的效果对比 5. 1. 1 不同谱估计方法的分辨率对比 图 1 为 Bur g 谱估计法与 Welch 谱估计法分辨 率对比图 , 由图中可见, Bur g 谱估计法能分辨出间 隔为 0. 02H z 的频率 , 而 Welch 谱估计法不能分辨 出 , 因此前者的分辨率明显高于后者。 图 2 为 Bur g 谱估计法与改进 Burg 谱估计法分 辨率对比图, 由图中可见, Burg 谱估计法的峰值有 明显偏 移, 而 改进 Burg 谱 估计 法 几乎 没 有峰 值 偏移。
[ 2] 2
4
改进协方差谱估计法原理
确定 AR 模型参数的 Burg 谱估计法的实质是 预测系数满足 L evinson 递推的、 附加约束的最小二 乘格型算法 , 由于这一约束的结果, 增加 AR 模型的 阶仅需要在各级最优化一个参数。与 Burg 谱估计 法相比, 可以利用无约束改进协方差谱估计法 ( 最小 二乘算法) 确定 AR 参数, 该算法的主要思路是摆脱 因采用递推运算对确定预测系数的约束 ( 在式 ( 2a)
3
改进Burg 谱估计法
众所周知, Burg 谱估计法首先从低阶开始 , 根
2 a2 ( 2) = - b( 1 ) c( 2 ) - c2 ( 1) b( 1) b( 2 ) - c ( 1) 利用 L evison 递推公式得
( 3d)
据数据估计反射系数 , 然后利用 L evinson 递推法计 算预测滤波器 ( PEF) 的系数, 从而避免了直接计算 相关值和矩阵求逆。而 Burg 谱估计法中一阶反射 系数 a 1 ( 1) 在计算功率谱密度( PSD) 中起着关键作 用, 这是由于在运用递推公式时, 各阶 PEF 系数均 与其有关。然而 , 在 Burg 谱估计法中 a1 ( 1 ) 的计算 本身就存在误差 [ 4] , 结果使谱估计的峰值产生偏差 , 特别是当信号长度为 1/ 4 信号周期的奇数倍且初相 位为 P/ 4( 或其奇数倍) 时误差更大。改进算法不直 接计算一阶反射系数 , 而是通过二阶 PEF 系数再求 一阶反射系数, 结果与初始相位无关, 从而使谱峰偏 差减小, 改善了谱估计的质量。 由于一阶反射系数的误差是引起谱估计偏差的 主要原因之一 , 因此从改进 a1 ( 1) 的计算入手, 将能改 善谱 估 计 的性 能。改 进 算法 的 基 本 思想 是 不 用 式( 2a) 直接计算 a1 ( 1) , 而是从二阶误差总功率着手, 先计算二阶 PEF 系数 a2 ( 1) , a2 ( 2) , 再计算 a1 ( 1) 。
N- 1
b p- 1
( n- 1)
f
( n - 1) + Q p e p- 1 ( n)
( 3a)
( 4) 仿照式 ( 2b) 的递推关系 , 重复上述过程, 直 到阶数 p 等于所需 AR 模型阶数, 求出 AR 模型参 数估值, 再按下式进行功率谱密度估计 P ^ BU ( f ) = 1+ E ^p
c( i ) =
n= 2
E [ x(n+ E [x
2
i - 2 ) x ( n- 2) + ( 3c)
x ( n - i) x ( n) ]
N- 1
b( i) = 由此解得
( n + i - 2) + x 2 ( n - i) ]
n= 2
a2 ( 1) = -
c( 1) [ b( 2 ) - c( 2 ) ] b( 1) b( 2 ) - c2 ( 1)
2 Ep = ( 1 - Q p ) E p- 1
二阶 PEF 输出误差功率为 ( 2b)
N- 1
E2 = =
n= 2 N- 1
E
2 [ ef2 ( n) ] 2 + [ eb 2 ( n) ]
E = R ( 3) 递推高一阶前、 后向预测误差, 即
p 2 wp
( 2c)
n= 2
E
[ x( n) + a2 ( 1) x( n- 1Leabharlann Baidu + a2 ( 2) x( n- 2) ] 2 +
2 2
小均方预测误差。 ( 2 ) 估计 p 阶预测系数 ( AR 模型参数) 及最小 预测误差的均方值 , 即 Q m = - 2 E e fm- 1 ( n) eb m- 1 ( n - 1)
n= m N- 1 n= m N- 1
2
Burg 谱估计法
为了克服 Levinson 递推算法中因估计 相关函
a 2 ( 1) = a 1 ( 1) + a2 ( 2) a1 ( 1) ( 3e) 由式( 3e) 可求出 a1 ( 1 ) , 再利用式( 3d) 得 c( 1 ) a2 ( 1) = ( 3f) a 1 ( 1) = b( 1) 1 + a2 ( 2) 综上所述 , 改进 Burg 谱估计法计 算步骤可归 纳为: ¹ 由式 ( 3c) ~ ( 3f ) 计算 a2 ( 2) , a2 ( 1 ) , a1 ( 1) ; º 当阶数大于 2 时, 利用式 ( 3a) ~ ( 3e) 递推计算反 射系数和 PEF 系数, 直至所需的阶数。 另外 , 还有一种改进方法称为 It akura 算法, 此 算法的关键 是改 变反射 系数 的形式 , 基 本思 路与 Bur g 谱估计法一致, 不同之处在于计算各阶矩时, 是用时间均值 来代替集合均值 , 推出 It akura 算法 反射系数, 即