Nautical Navigation

Size: px

Start display at page:

Download "Nautical Navigation"

遵东
6 years ago
Views:

1 统计信号分析与处理第 6 章波形估计

2 本章内容 6. 波形估计的分类 6.2 连续信号的维纳滤波 6.3 离散维纳滤波 α β 6.4 滤波 6.5 卡尔曼滤波

3 6. 波形估计的分类假设观测信号为 g( ) s( α ) x = s + n = + 0, 当前时刻为, 待估计信号为 x( ) = s( ) + n( ) s( 0 ) 0 滤波 : 估计信号 s ( 0 ) s 0 + α, α > 0 0 预测 : 估计信号 ( α ) s 0 + α > 0 s 0 + α, α < α 平滑 : 估计信号 ( α ) s 0 + α < α 0

4 波形估计目的 : 选取线性滤波器的冲激响应函数或传输函数, 使估计的均方误差达到最小拓展 :

5 6.2 连续信号的维纳滤波基本思想 : 寻找线性滤波器的最佳冲激响应或传输函数, 使滤波器的输出波形作为输入信号波形的最佳估计, 即使波形估计的均方误差达到最小主要内容 : 广义平稳随机信号的维纳滤波原理物理不可实现维纳滤波器的解物理可实现维纳滤波器的解最小均方误差

6 6.2. 广义平稳随机信号的维纳滤波原理假设信号和噪声均为零均值的广义平稳随机过程且互不相关, 则观测信号平稳随机过程, 且和 x( ) = s( ) + n( ) g( ) 联合广义平稳也是零均值广义线性滤波器的冲激响应为, 输出信号 y 为 : h + + y h τxτdτ hτx τdτ = = 估计误差为 : 均方误差为 : = ˆ = e g g g y { 2 } { 2 ()} () () E e = E g y min

7 利用线性最小均方误差估计的正交条件 { } 0, E e x τ = < τ <+ 可得非因果广义平稳随机过程的维纳 - 霍夫方程 gx, + R η = h λ R ηλ dλ < η <+ x 因果广义平稳随机过程的维纳 - 霍夫方程 + R η = h λ R ηλ dλ 0 η <+ gx 0 x

8 6.2.2 物理不可实现维纳滤波器的解由非因果维纳 - 霍夫方程, 可得 gx ( ω) = ( ω) ( ω) S H j S x 即对于 H ( jω) g() = s( + α ) = S S gx x, 有 ( ω) ( ω) H ( jω ) = S s jωα s ( ω) e ( ω) + S ( ω) S n

9 维纳滤波器的幅频相应为 H ( jω ) = S s s( ω) Sn( ω) ( ω) ( ω) + + S S SNR ω n ( ω) = ( ω)/ ( ω) 其中 SNR S S s n H ( jω) 此时估计的最小均方误差为 : { ()} + 2 min g( 0) gx E e R h λ R λ dλ = 0 SNR( ω)

10 例 6. x 考虑输入信号为高斯 - 马尔可夫信号 s 和噪声的叠加, 假定信号和噪声统计独立, 其功率谱分别为 n 3 S s = + ω ( ω) 2 和 S ω = n 求 : α 分别取 0,+ 和 - 时物理不可实现维纳滤波器的冲激响应及最小均方误差

11 6.2.3 物理可实现维纳滤波器的解对于因果的维纳 - 霍夫方程 : gx 0 + R η h λ R η λ dλ η = x <+ 其求解方法有如下两种 :. 频谱因式分解法 2. 预白化法

12 . 频谱因式分解法则有 : 引入一个新函数 : b( η ) gx 0 η 0 = 形式不作具体规定 η < 0 + R η + b η = h λ R ηλ dλ < η <+ 拉氏变换后得 : gx + = S s B s H s S s x x 可知 H(), s B() s 半平面内的极点分别只存在于 s 的左半平面和右

13 由于如果 S s 满足佩利 - 维纳条件, 则有 : x + = S s S s S s x x x * ( ω) = ( ω) = ( ω) S S S x x x ( ) = S s S s + x x, 则于是有 + ( s) H s S x + S s B s S s S s gx gx gx = + = + + x x x x x B s S s S s S s S s S s 故 H s + + αs Sgx s Ssx s e = + = + Sx ( s) Sx s Sx s Sx s

14 2. 预白化法基本思想 : 先将输入信号进行预白化处理, 然后让得到的白化信号通过其对应的维纳滤波器 x( ) z( ) H ( s ) H ( s) W y( ) 白化滤波器取为 H W ( s) = S + x ( s)

15 对于白信号输入为, 有 R τ = δ τ, 其对应的维纳滤波器满足如下方程 : 有 : z( ) + Rgz τ = h λ Rz τ λ dλ = h τ, τ 0 0 τ s z ( s) ( s) Sgx = τ τ= 0 gz = W gx = Sx H s h e d S s H s S s 则 : S ( s) gx = W = + Sx ( s) Sx s H s H s H s 最小均方误差为 : { 2 ()} () () + { } s( 0) gx E e = E e g = R h λ R λ dλ min 0 +

16 例 6.2 n x 考虑输入信号为高斯 - 马尔可夫信号 s 和噪声的叠加, 假定信号和噪声统计独立, 其功率谱分别为 3 S s = + ω ( ω) 2 和 S ω = n 求 : α 分别取 0,+ 和 - 时物理可实现维纳滤波器的冲激响应及最小均方误差

17 6.2.4 最小均方误差物理可实现维纳滤波器的最小均方误差可以表示为 : 可推得 : { ()} + 2 = s( 0) sx( + ) E e R h τ R α τ dτ min 其中 φ τ = Rsz τ = E s z τ { 2 ()} + 2 s 0 E e R φ τ d τ min α = { } 可见 : 最小均方误差随着 α 值的增大而单调地增大

18 2 φ ( τ ) 2 φ ( τ ) 2 φ ( τ ) α τ α τ α (a) α > 0 (b) α = 0 (c) α < 0 τ 预测滤波和平滑的性能示意图结论 : 平滑用到的数据最多, 对应的阴影面积最大, 其误差最小 ; 而预测用到的数据最少, 其对应的阴影面积最下, 误差最大 ; 滤波所用到的数据居中, 其误差也位于前二者之间

19 6.3 离散维纳滤波为对于随机输入序列 = ( + ) g k s k α x( k) = s( k) + n( k) 系统的输出为 + + = ( ) = ( ) y k h k j x j h j x k j j= j=, 待估计的信号序列 h( k ) 2 E{ e ( k) } = E g( k) y( k) 选取使得 y k 的估计均方误差 { 2 } 达到最小, 这时的离散线性系统称为随机序列的维纳滤波器其解满足线性最小均方误差估计的正交条件 : + E g k h j x k j x j j= ( ) ( ) 0

20 有限观测样本的维纳滤波本, 即每次对于 g k 的估计, 只使用长度为 N 的观测样则维纳滤波方程如下 : N i= 0 N = ( ) y k h j x k j j= 0 () () hir l i = R l l= 0,,, N x gx 矢量形式为 : Rh= r, 则 : h x gx = R r x gx 2 均方误差为 : N { } = g ( 0) gx E e k R h j R j min j= 0

21 x( k) = s( k) + n( k) h ( 0) z z h h( N ) h( k) s ( k ) z z z Σ g( k) = s( kα ) = ˆ = ˆ( ) y k g k s k α e( k) + +

22 下面通过举例引入 6.4 α β 滤波滤波假设常值信号中混杂了加性噪声, 现用 N 个观测样本 x = s+ n k k 来估计常值信号 s 一种处理方法是 : α β x x. 测量, 存储, 取样本均值估计量为 ŝ = x ; x2 x 测量, 存储, 取样本均值估计量为 sˆ = ( x + x ); 3. 依此类推, 测量 xn, 存储 xn, 取样本均值估计量为 sˆ N = ( x+ x2 + + xn) N 2

23 另一种处理方法是 :. 测量 x, 取样本均值估计量为 ŝ= x, 存储 ŝ, 清除 x ; 2. 测量 x2, 取样本均值估计量为 sˆ2 = ( x ˆ 2+ s), 存储 ŝ2, 2 清除 x, s ˆ ; 2 x N 3. 依此类推, 测量, 取样本均值估计量为 N sˆ ˆ N = sn + xn, 存储 s, 清除 sˆ ˆN N, xn N N

24 例 : 需要对一运动目标进行跟踪测量, 估计出目标的位移与速度假设目标在时刻的位移为 k, 采样间隔 T, 并且假设在时间区间测样本为 yk = xk + nk k [, + T] k k x 内目标作匀速直线运动观做法 : 用 x 表示对目标位移的估计值, 则 x 为该时刻 ˆk 目标的速度估计值 x k ˆk 建立预测方程 : xˆ = xˆ + Tx ˆ x ˆ k+ k k k+ = xˆ k

25 利用新的观测值建立更新方程 : 矢量形式 : y k + 上述递推算法称为 α ( ) xˆ = xˆ +α y xˆ k+ k+ k+ k+ ˆ ˆ β x = x + y xˆ T ( ) k+ k+ k+ k+ T T Φ =, = [ α β / T ], = [ 0] 0 k h ( y ) xˆ = Φxˆ + k hφxˆ xˆ = [ xˆ, xˆ ] T k+ k k+ k β 滤波 k k k

26 递推运算流程 y k + xˆ k + k 现刻滤波估计 h xˆ k z Φ xˆ k + 一步预测估计

27 6.5 卡尔曼滤波主要内容 : 6.5. 状态空间模型离散卡尔曼滤波

28 6.5. 状态空间模型. 状态变量状态方程和观测方程状态变量 : 描述系统状态的一组数目最少的变量状态方程 : 状态变量所满足的一阶微分方程观测方程 : 描述系统输出与状态变量之间关系的方程

29 例 :RLC 电路由基尔霍夫定律 : i () () o () = () u = R+ r i di + L + i( τ ) dτ d C u R i

30 利用 dq i () =, q () i( τ ) dτ d = 2 可得 dq () () dq L + R+ r + q = u 2 d d C dq uo = R d 令 () () () () dq q = x, = x2, y = uo d () () () i 可得 = x x2 R+ r u x 2() = x() x2() + LC L L y() = Rx2 () ()

31 写出矢量形式 () () ( + ) () x 0 x 0 = u x LC R r L + L 2 x2 x () y() = [ 0 R] x 2 () 一般的矢量形式 () () = ( ) ( ) + ( ) ( ) x F x G u () = ( ) ( ) + ( ) y H x n 状态方程观测方程如何用状态空间模型来描述随机过程?

32 2. 状态方程的解状态方程是一个一阶 P 元联立的线性微分方程组, 其齐次方程组 : x () = F() x() 有 P 个线性独立的解 2 p x x x (), (),, () 这组解构成的非奇异矩阵, 称为上式的基本解阵满足微分方程 X x x x 2 p () = [ (), (),, ()] dx() d = F() X()

33 定义 : 对于, τ 0, 矩阵 Φ(, τ ) = X() X ( τ ) 为系统的状态转移矩阵状态转移矩阵具有以下性质 : a. 对于 τ 0 Φ (, τ ) = F () Φ (, τ ) b. 对于 0 Φ(,) = I c. 对于, τ 0 Φ = X X = Φ (, τ) ( τ) () ( τ,) d. 对于 e. 对于, τ Φ(, ) Φ(, ) = Φ(, ), 2, Φ(, τ) =Φ(, τ) F( τ) τ

34 对于因果系统, 状态方程的唯一解可由状态转移矩阵表示为 : 验证初始值 : 两侧对求导 : x( ) Φ(, ) x( ) Φ(, τ ) G( τ) u( τ) dτ = (, ) (,) ()() 0 x = Φ x + Φ τ Gτ uτ dτ = x Φ(, 0) Φ(, τ) x () = x0 + Φ(, τ) G( τ) u( τ) τ = + G()() τ uτ dτ 0 = F() Φ(, ) x + G() u() + F() Φ(, τ) G( τ) u( τ) dτ 0 0 = F() x() + G() u() 0 0

35 x( ) Φ (,0) x(0) Φ (, τ ) G( τ ) u( τ ) dτ = + 0 零输入解零状态解对于线性非时变系统, 相关矩阵退化为常数矩阵 : x () = Fx() + Gu() y() = Hx() + n() 0 0 可以证明此时的解为 : x() = Φ( ) x( ) + Φ( τ ) Gu( τ) dτ 0 0 Φ(, τ) ( τ) e F ( τ ) = Φ = = 0 + k = 0 F k ( τ ) k! k

36 3. 离散状态方程和观测方程对状态方程的时域解采样 : x Φ x Φ τ Gu τ dτ k+ ( k+ ) = ( k+ k) ( k) + ( k ) + k 可得 : x = Φ x + u k+ k+, k k k y = H x + n k k k k

37 6.5.2 离散卡尔曼滤波假设在时刻 k, 基于该时刻以前所获得的全部知识, 对状态矢量 x 做一个预测估计 x 预测估计的误差为: k 预测误差的协方差矩阵为 : ˆ k e = x xˆ k k k { T } { T} E ( ˆ )( ˆ ) C = E ee = x x x x k k k k k k k 利用时刻的新观测数据 y k k 进一步改善对 xk 的估计 xˆ k, 称为更新估计 ( ) xˆ = xˆ + K y H xˆ k k k k k k

38 更新估计的误差为 : e = x xˆ = IK H e K n 更新估计的均方误差为 : E k k k k k k k k { } { T } T ee k k E ( xk xk) ( xk xk) k = ˆ ˆ K min 更新估计的协方差矩阵为 : { T } E ( ˆ ) T { ( ˆ ) } C = E e e = x x x x k k k k k k k 其中 : E{ nn T } T T k k k k k k k k = IK H C I K H + K R K = Rδ k i k ik

39 令 DD = HCH + R T T k k k k k k C = C + K D A K D A AA 可得 k k ( k k )( k k ) A= C H D T T k k k T T 当此时 C k KD A= 0 k k 的主对角线元素之和达到最小 T 即 ( T) T ( T K ) k = ADk = CkHk Dk Dk = CkHk HkCkHk + Rk 此时 = ( ) C I K H C k k k k

40 xˆ, C 下面利用预测估计 xˆ k+, Ck+ 由状态方程忽略噪声作用 k k x = Φ x + u + x ˆ k k k k Φ x ˆ = k+ k k 则 k + 时刻的预测误差和协方差矩阵分别为 e = x xˆ = Φ e + u k+ k+ k+ k k k { T} C = E e e = Φ C Φ + Q T k+ k+ k+ k k k k 其中 E{ uu T } = Qδ k i k ik

41 卡尔曼滤波算法的递推过程 () 建立状态空间模型 : xk+ = Φkxk + uk y = H x + n k+ k+ k+ k+ (2) 设置初始化条件 : xˆ = E{ x }, C = Var{ x } x ˆ (3) 预测 : Φ x ˆ = k+ k k (4) 计算预测误差协方差矩阵 : C = Φ C Φ T + Q + k k k k k T (5) 计算卡尔曼增益 : ( T ) K = C H H C H + R k+ k+ k+ k+ k+ k+ k+

42 (6) 更新 : ( xˆ ˆ ˆ ) k+ = xk+ + Kk+ yk+ Hk+ xk+ (7) 计算更新误差协方差矩阵 : = ( ) C I K H C k+ k+ k+ k+ (8) 令 k = k+, 重复步骤 (3)-(8) 直到当前时刻迭代框图 : y k + + K k xˆ k 现刻滤波估计 Φ k Η k z xˆ k + 一步预测估计

43 例 6.3 飞机相对于雷达作径向直线匀加速运动目标与雷达站之间的斜距为 x k 程和加速度方程为目标相对于雷达的距离方程速度方 2 T x = x + Tx + x, x = x + Tx, x = x 2 k k k k k k k k k 式中, 为时刻目标的斜距, 为时刻目标的径向速度, 为时刻目标的径向加速度, 为采样时间间隔, T = 2s xk k k x k T k x k

44 2 已知 E{ x 0} = 0, Var{ x0} = 8( km) E{ x } = 0, Var{ x } = 0( km s) 与 k k k 都不相关当观测时刻 k = 时, 斜距的观测值 y km 为 : k 0 0 { } { } ( 2 ) 0 = 0.2, 0 = 5 E x Var x km s 且目标未受外界干扰, 动态噪声 u = 0, { } = 0, {, } = 0.5δ 2 白噪声, 即 E n Cov n n km k i j ij 是零均值的 n x, x, x, 2,,0 k 0.30,.56, 3.64, 6.44, 0.5, 4.8, 20.0, 25.2, 32.2, 40.4 k n k 2 2 求 : x, x, x k k k 的估计以及估计误差的协方差矩阵

! /. /. /> /. / Ε Χ /. 2 5 /. /. / /. 5 / Φ0 5 7 Γ Η Ε 9 5 /

! /. /. /> /. / Ε Χ /. 2 5 /. /. / /. 5 / Φ0 5 7 Γ Η Ε 9 5 / ! # %& ( %) & +, + % ) # % % ). / 0 /. /10 2 /3. /!. 4 5 /6. /. 7!8! 9 / 5 : 6 8 : 7 ; < 5 7 9 1. 5 /3 5 7 9 7! 4 5 5 /! 7 = /6 5 / 0 5 /. 7 : 6 8 : 9 5 / >? 0 /.? 0 /1> 30 /!0 7 3 Α 9 / 5 7 9 /. 7 Β Χ9