时频子空间拟合波达方向估计

时频子空间拟合波达方向估计 1 2 2 梁,殷勤业,李盈金 ()11 解放军信息工程大学 ,郑州 450002 ;21 西安交通大学信息工程研究所 ,西安 710049 摘 要 : 本文提出了一种基于信号空时特征结构的时频子空间拟合方法 ,利用双线性时频分布构造时频相关矩阵 C代替传统的阵列相关矩阵 R, 通过 C的特征分解实现了信号子空间与噪声子空间的分离. 该方法在空域和二 x x x 维时频域同时进行处理 , 能够区分具有不同时频特征的信号 , 既适用于平稳信号的场合又适用于时变 、非平稳信号的 情形 , 属于空时多维处理的范畴. 可以证明 , 基于平稳信号假设的经典子空间方法是该方法的低维特例. 由于包含了时 变滤波的过程 , 因此该方法具有信号选择性以及抗干扰和抗噪声的能力. 仿真结果证实了该方法的有效性. 关键词 :阵列处理 ; 联合时频分析 ; 波达方向 ; 非平稳信号 A () 中图分类号 :TN91117 22112 20010120071204 0372文献标识码 :文章编号 : Estimating Dire ctio n2of2Arrival Via Time2Fre que ncy Signal Sub sp a c e Fitting 1 2 2J IN Liang, YIN Qin2ye,L I Ying (11 University of Inf ormation Engineering , Zhengzhou 450002 , China ; )21 Inf ormation Engineering Institute , Xi’an Jiaotong University , Xi’an 710049 , China Ab stract : Based on the decomposition of the space2time eigenstructure of signals ,a novel time2frequency signal subspace fitting ( ) TF2SSFmethod is proposed in this paper to estimate the DOA of signals. Through the bi2linear Cohen class time2frequency distribu2 tion ,the time2frequency correlation matrix Cis constructed to replace the traditional correlation matrix R. Accordingly the signal x x subspace and noise subspace are separated with the eigen2decomposition of C. Because the observed data are processed in spatial do2 x main and 22D time2frequency domain simultaneously ,the method can separate the signals that have different time2frequency signatures and is suitable for both stationary and nonstationary signals ,while the traditional subspace methods have to assume that the signal is stationary. It is also proven herein that the traditional subspace methods is the special case of the TF2SSF method. Furthermore ,with the time2varying filtering available ,the method has the signal2selectivity and is capable of suppressing interference and noise that are difficult to handle in time or frequency domain only. The simulation results show the effectiveness of the method. () Key word s : array processing ;joint time2frequency analysis ;DOADirection2of2Arrival; nonstationary signals 时域处理 ,更充分地开发信号中的有用信息 ,成为进一步提高 1 引言 ) (DOA 波达方向估计是阵列信号处理领域研究的重要内 其性能的关键 ,也是近年来人们感兴趣的研究课题之一. 众所 容之一. 由于具有良好的超分辨力 、稳健性和相对较小的计算 周知 ,联合时频分析是对非平稳信号或谱时变信号进行处理 量 ,子空间方法已成为 DOA 估计中最常用 、最经典的方法. 然 1 ,6 的有效手段,因此我们可以将时频分析的方法与阵列信 而 ,在传统子空间方法的信号模型中并没有充分利用信号本 号处理相结合 ,通过时频分布将信号变换到时频域 ,利用时变 身蕴含的一些非空域特征 ,该类方法一般假设信号是平稳信 滤波提高空间谱估计方法的性能. 由于将一维时域信号映射 号 ,通过各态历经性假设以及对信号的随机采样 ———快拍 ,用 到二维时频域中 ,因此能够在高维空间中更精细 、更准确地刻时间平均代替集合平均从而构造出阵列协方差矩阵. 但在实 际应用中许多典型信号是非平稳的或谱时变的 ,如雷达中的 画和反映时变信号的特征和细节 ,从而可以利用时频二维掩 3 ,7 ,8 线性调频信号 、通信中的跳频信号等 ,因此传统子空间方法对 ( 模等时变滤波手段,将一些在低维空间 如单纯时域或 这些信号进行 DOA 估计时存在先天性不足. ) 频域中难以区分 、但具有不同时频特征的信号加以分离 ,同 随着阵列信号处理技术日益广泛的应用 ,在许多场合中 ( 时有效地抑制干扰 ,使得 DOA 估计方法具有信号选择性 Sig2 () 信号是配合信号 如通信中或已有一些先验知识可资利用. ) nal Selectivity以及更好的分辨力 、过载能力 、抗各种干扰和有 如何在用子空间方法进行空域处理的同时有效地引入适当的 色/ 无色噪声的能力 ,并且既适用于平稳信号又适用于时变 、 非平稳信号的 DOA 估计. 根据这一思想 ,本文提出一种基于 电子学 报 72 2001 年 ( ) 信号空时特征结构的 DOA 估计方法 ———时频子空间拟合方假设 1 已知有用信号 SOI时频分布的先验知识 , 有用 ( ) 法 , 从一个新的角度研究和探讨空间谱估计的空时多维处理 信号与其它信号具有不同的时频特征 分布;手段. ( ) 假设 2 各源信号 st之间不相关 , 信号与阵列噪声之 i 间互不相关 ; 2 信号的二次时频分布2 σ假设 3 阵元噪声为均值为零 , 方差为 的平稳随机白 ( ) 对于谱时变的确定性信号 s t, 可以通过双线性变换得 噪声 , 谱密度在时频面上均匀分布. 到它的时频分布 : 与传统模型不同 , 这里并不需要假设各源信号是平稳的. ( ) 在许多场合 如通信中 , SOI 是配合信号 , 其时频分布是已知 3 ω) (θτ) ( (τ) ) ( (τ) )( Css t ,= <,s u + / 2s u - / 2 µ ( 的 ;在其它非配合场合中 , 可以利用 SOI 的一些先验知识 如 - jωτ- jθt + jθu θτ ?e dddu()1 ) 调制方式 、频率等, 通过对阵元输出信号的时频分析来估计 [1 , 2 ] (θτ) 上式即通常所称的 Cohen 类分布, 其中 <,为核函数. 时频分布. 定义“阵列时频相关矩阵”C为 :`x 对于非平稳的随机信号 , 其时频分布可表示为数学期望 H `C> E{ Q [ xx} x 的形式 : E{ Q [ x, x} , x} E{ Q [ x, x} E{ Q [ x1 1 2 2 1 M ( ω) ( ω) (θτ) ( (τ) )Ct ,= E[ Ct ,] = E <,s u + / 2 `ss ss E{ Q [ x, x} ω > 2 1 µ ωτθθ3 - j- jt + ju ( (τ) ) θτ ?s u - / 2e ddudE{ Q [ x , x } E{ Q [ x , x } M 1 M M ωτ θ θ - j t + j u ( ω)- j( ω)( ω)C`t , `C Ct , t , ` x x x xx x ()θτ = (θτ) ( τ) 2 ddud<,Ru ,e 1 1 M 1 2 1 s µ ( ω)`Ct ,ω x x ( τ) ( ) 其中 Rt ,为随机信号 s t的时变自相关函数. 由此可见 , s 2 1 ()10 = 与平稳随机信号的情形 类 似 , 非 平 稳 信 号 的 时 变 功 率 谱 ( ω) ( Ct ,也是由时变自相关函数变换得来的 , 这里称 `Ct ,ss ss ( ω)( ω)`Ct , `Ct , x x x x M 1 M M ω) ( ) 为 s t 的“自时频分布”. 类似地 , 可以由两个随机信号 ( ω) 其中 Q [ ?]表示双线性变换. 矩阵 C`为 t ,的函数 , 其 x ( ) ( ) st、st的互相关函数得到“互时频分布”.1 2 主对角线上的元素为各阵元输出的自时频分布 , 其它元素为 τ τ 3 () ( ) ( )各阵元输出间的互时频分布. 将式 8代入 , 并考虑式 6、7 ( ω) (θτ) ( ) ( ) Ct ,= E <,su + su - `s s 1 2 1 2 µ 2 2 及假设 2 、3 得 :ωτθθ- j- jt + juθτ ?e ddudH ( ) ( ) = E{ Q [ As + nAs + n} `CxH H H τθθ- j w- jt + ju = E{ A?Q [ ss ] A+ A?Q [ s?n] θτ (θτ) ( τ) ()= ddud<,Ru ,e 3 s s 1 2 µ H H H ) ( + Q [ n?s ] A+ Q n?n} 将双线性变换看作是从一维时域空间到二维平面的映H H H = A?E{ Q [ ss ] A+ A?E{ Q [ s?n}射 , 定义为 Q [ ?, 则上式可重写为 : H H H )( } ()+ E{ Q [ n?s } ?A + E{ Q n?n4 ( ω) ( ) ( ) Ct ,= Q [ s t, s t] ss HH H ))( ( = A?E{ Q ss A + E{ Q n?n } }()( ω) ( ) ( ) 5 Ct ,= Q [ st, st] s s 1 2 1 2 H 2 ( )σI11 = A?`C?A + s () 由式 2可得下列性质 : H( ) ()其中 : C= E{ Q ss } 12 `s ( ω) ( ) ( ) Ct ,= E{ Q [ s t, s t} `ss ( ) ( ) = Q{ E[ s t, s t} ( ω)( ω)( ω) `Ct , `Ct , `Ct , s s s s s s ( τ) = Q [ Rt ,] 1 1 1 2 1 d s ()6 ( ω)ω `Ct , ( ) ( ) s s = Q{ E[ st, st} ( ω) ( ) ( ) Ct ,= E{ Q [ st, st} `1 2 s s 1 2 2 1 1 2 ()= 13 )(7 ( τ) = Q [ Rt ,] s s 1 2 ( ω)( ω)`Ct , `Ct , s s s s d 1 d d 3 时频子空间拟合的基本思想 为信号的时频相关矩阵 , 其主对角线上的元素为各信号源在 ( ω) t ,上的功率谱密度 , 其它元素为各源之间的交叉项 , 考虑 考虑采用 M 元的传感器阵列 , 假设各阵元对输入信号的 到假设 2 , 这些交叉项应为零. 则 C`为 d ×d 的对角阵 : ( ) x 响应是线性时不变的. 有 d 个窄带信号 st 分别以 DOA i ( ω) ( ω) ( ω) ( ω) ()Ct ,= diag[ Ct ,, Ct ,, , Ct ,] 14 ````θ( ) s s s s s s s i = 1 , 2 , , d入射 , 阵列输出信号矩阵为 : i 1 1 2 2 d d ( ) 假定 st为期望信号 , 选择其时频分布“峰脊”上的一点 ( ) ( )()( ) = As t+ n t 8 x tl ( ω) t,, 在该点上其它信号的时频分布相对很小 , 则上式对 其中 :l l T T () ( ) ( ) ( ) θ) θ) 角线元素中只有一项不为 0 , 这时式 11变为 : ((x t= [ xt, , xt] , A = [ a , ] , , a 1 M 1 d H 2 TT )(( ω) ( ω) (θ) (θ) σ15 Ct,= Ct,a ?a +I ``( ) l ()( ) ( ) ( ) ( ) , nt] ( ) 9 x l l s s l l l s t= [ st, , st] , n t= [ nt, M 1 d 1 l l ( ) ( ) xt , nt 分别为第 i 个阵元上的输出和噪声 , A 为 M ×d i i ( ω) 这样 , 可以用时频相关矩阵 `Ct,代替传统的阵列相关矩 x l l θ的阵列导向矩阵 , 是待估参数 的函数. 当阵列为均匀线阵 i ( ω) ( ) 阵 R, 通过对 `Ct,的特征分解得到信号 st的 DOA 估 x x l l l 时 , A 为 Vandermonde 矩阵. 计. 称这种基于信号时频分析的 DOA 估计方法为时频子空间 在以下讨论中做下列假设 :( ) 拟合法 ,简称 TF2SSF Time2Frequency Signal Subspace Fitting. 显 73 第 1 期金 梁 :时频子空间拟合波达方向估计 ω) ( 然 , 这时信号的“特征结构”已不是仅仅由信号的空域信息所1 , t ,?D M()ω) 20 ( h t ,= 决定 , 而是信号时频特征和空域结构的综合反映 , 因而是一种 (ω) , t , 0 其它 更广义的“空时特征结构”. ω) ( 对阵列输出进行时变滤波的过程 , 其中 D为{ t , | t- M l (Δ) (Δ) ω( Δω) ( Δω) ω ωt/ 2Φ t Φ t+ t/ 2, - / 2ΦΦ+ / 2} ,l l l 4 时频子空间方法与传统方法的关系ΔΔωt 、为掩模的时宽和频宽 , 取决于时频变换时离散点的个 从上述推导中不难发现 TF2SSF 方法与传统方法子空间 数. 总之 , 选择合适的时频点等价于选取适当的时变滤波器 , (θτ) 之间的联系. 实际上 , 如果选取核函数 < ,使时频分布满 足分别将混合信号中的有用信号一一提取出来 , 因而在进行 边缘条件 , 即 : DOA 估计时不必考虑其它信号的影响 , 即具有信号选择性 , 它使得 TF2SSF 算法在具有较强的抗干扰能力的同时 ,还具有 3 ()( ω) ω( ) ( ) 16 `Ct ,d= E[ xtx t]xxi j ij ?分辨率高 、过载能力强等良好性能. 则有 : 512 抗噪性能( ω) ω( ω) ω ?Ct ,d ?Ct ,d`` x x x x1 1 1 m 由于信号能量一般集中在较小的时频区域内 , 而随机噪 ( ω) ω?`Ct ,d ω x x 声在时频平面上的能量分布范围较宽且相对比较均匀 , 所以 2 1 ( ω) ω `Ct ,d>x ?( ω) 尽管在观测时段上信噪比较差 , 而在 t,的邻域 D上信 l l M ω 号能量却明显大于噪声能量 , 因此时变滤波亦能提高 TF2SSF( ω) ω( ω) ω?Ct ,d ?Ct ,d ``x x x x m 1 M M 算法的抗噪能力. 设传感器观测带宽为 B , 观测时间为 T ; 信3 3 )( )( E xx E x x 1 1 M 1 ( ) 号 s t时频分布的支集为 D , 时频面积为 s , 为讨论方便并D ()17 = 不失一般性 , 假设在 D 上信号功率谱的分布是均匀的 , 即3 3 )( )( E x x E xx 1 M M M ( ω) Ct ,| = C, 则在观测时频区域内信号能量为 E=``ss D ss s H )( ( ω) ω18 ( ) ( ) 即 :?`Ct ,d= E[ x tx t] x 2 ω σ ( ω) ω κ Ct ,dtd= C?S, 噪 声 能 量 为``E= κ dtd=n ss ss D 0 < t < TD 上式表明 , 时频相关矩阵 `C在时间轴上的投影就是传统 x ω0 << B 2 相关矩阵 R; 进而 , 由于 `C是二元函数 , 因此 R张成的空间 σ B T , 因此信噪比为 :x x x span{ R} 是 `C张成的空间 span{ `C} 在与频率正交的相应空 E`C?s x x x s ss D()SNR = = 21 2 σ E nB T 间上的投影 , 或者说 span{ C`} 是 span{ R} 的高维延拓. 所以传 x x () 经过式 21 时变滤波器后的信噪比为 :统子空间方法可看成是 TF2SSF 方法的低维特例. 显然 ,在高 TF ΔωΔ E`Ct `CS ssss B T维空间中表征信号能更精确地刻画出信号的本质与细节 ,能 ()22 SNR = = = = SNRTF TF 2 2 σ S σΔωΔEt Dn 将信号子空间与噪声子空间更好地分离. 实际上 ,这正是 TF2 一般情况下 B T µ S, 所 以 时 变 滤 波 后 信 噪 比 大 大 提 高 , D SSF 方法具有信号选择性 、抗噪性 、高分辨率等一系列优良特 B T 2 2 性的根本原因. σ() σSNRµ SNR , 同时 C= SNR µ, 则式 11可近似为 : `TF ss S DH5 TF2SSF 方法性能的简单分析 () 23 `C? AC`Ax s ( ) 511 时频点的选取与信号选择性 当各阵元上的噪声是相关噪声时 , 式 11中噪声的时频2 ( ω) 选择合适的时频点 t, 上的 `C是 TF2SSF 方法的关 σ分布矩阵不再是对角阵 I , 但由上面分析可见噪声对 C的 `l l x x ( ) ( ω) 键. 假设 st的时频分布的支集为 D, 则选择 t,?D并 l l l l l 贡献极小 , 因此 TF2SSF 方法对各种噪声均不敏感.( ω) ( ) 且t,z D, i = 1 , , d , i ?l . D为其它信号时频分布 l l i i ( ω) 的支集. 这时在 t, 上分布的时频功率谱密度基本源于 l l 6 仿真结果 假设有两个覆盖相同时段和频段的窄带线性调频信号同 ( ) ( ω) ( ω) t, | Ct,| µ | Ct,| sl s s l l s s l l l l i i 时入射于 6 元均匀线阵 , 其 Wigner 分布如图 1 所示. 两信号的 ( ω)O`Ct, ( )1 ×d - 1 s s l l l l ( )DOA 分别为 60 度和 63 度 ,信噪比为 3dB. 一个信号为 : sn ( ω) ()1 因此 : `Ct,= 19 s l l 2 2 O ( ) O d - 1×1 ( ) ( )d - 1×d - 1 ω αωα jn + j n j n + j n 1 1 2 2 ( ) π ω = e . 其中 = 2? , 另一信号为 : s n = e1 2 因此通过选择合适的时频点可以有效地筛选出所需信号 , 大ωπαα011 ,= 2?011 ,= 010005 ,= - 010005. 快拍数为 256 次 , 2 1 2 大降低其它信号的干扰 , 从而使 TF2SSF 方法有很强的信号选 择能力. 从干扰抑制的角度来看 ,许多干扰信号占据与有用信号 相同的时段和频段 ,用一般时不变滤波方法很难完全滤除 ,因 而将不可避免地影响一般 DOA 估计方法的精度与分辨率. 而 干扰信号与有用信号往往具有不同的时频分布 ,可以在时频 3 ,7 ,8 ( ) 平面上通过对信号加掩模 mask的方法进行时变滤波, 从而更有效地滤除干扰. 实际上 , 上述选取时频点的过程可看 作是用掩模 图 1 两个覆盖相同时段和频段的窄带线性调频信号的 Wigner 分布 电子学报 74 2001 年 采用 Wigner 分布. 分别选取两信号各自能量山脊上的时频点同时进行处理 , 并且包含了时变滤波的过程 , 因此具有信号选( ) () ( ) ( ) t , f = 10 , 01105和 t , f = 10 , 01223用 TF2SSF 方法对它 择性 、抗干扰能力和抗噪声的能力 , 既适合于平稳信号又适合 ( ) ( ) 们进行 DOA 估计 , 图 2 a、b给出了 50 次 Monte2Carlo 实验 于谱时变信号的 DOA 估计 ,属于空时多维处理的范畴. 可以 的结果. 作为比较 ,图 3 给出了在信噪比为 20dB 下用 MUSIC 证明适用于平稳信号 DOA 估计的传统子空间方法是 TF2SSF 方法进行估计的结果. 由于两信号具有谱时变性 ,且它们的 方法的低维特例. 仿真结果证实了上述结论. DOA 相差不大 ,所以 MUSIC 方法无法正确分辨. 而 TF2SSF 方 法可以利用信号选择性十分准确地估计出它们各自的 DOA. 需要指出的是 ,本文的方法是基于窄带假设的 ,不适合于 表 1 给出了在各种信噪比下的估计方差. 宽带信号. 对于宽带时变信号的时频 DOA 估计方法 , 作者将 在今后的工作中提出并详细讨论. 参考文献 : 2frequency distributions2a review A . Proc . of IEEE L . Cohen. Time1 () C,1989 ,77 7:941 - 981 . L . Cohen. Time2frequency analysis M . New York : Prentice 2Hall , 2 1995 . S. Qian ,D. Chen. Joint time2frequency analysis M. New York : Pren 2 3 tice Hall ,1996 . A. Belouchraini ,M. G. Amin. Blind Source separation using time2fre2 4 quency distributions :algorithm and asymptotic performance A . Proc . of IEEE International Conf . Acoustic , Speech , Signal Processing C , 1997 ,5 :3469 - 3472 . Yimin Zhang ,M. G. Amin. Spatial averaging of time2frequency distribu2 5 tions A . Proc . of IEEE International Conf . Acoustic , Speech , Signal Processing C,1999 ,5 . ( ) 图 2 TF2SSF 方法的 50 次 Monte2Carlo 仿真实验结果. a窄带信 6 A. M. Rao ,D. L . Jones. Nonstationary array signal detection using time2 ( ) ( ) 号 sn的 DOA 估计结果 ,实际 DOA 为 60 度; b窄带信 1 frequency and time2scale representations A . Proc . of IEEE Interna 2 ( ) 号 sn的 DOA 估计结果 ,实际 DOA 为 63 度 2 tional Conf . Acoustic ,Speech , Signal Processing C ,1998 ,4 : 1989 - 1992 . 7 G. Boudreaux ,T. Parks. Time2varying filtering and signal estimation us2 ing wigner distribution synthesis techniques J . IEEE Trans. ASSP , 1986 ,34 :442 - 451 . B. E. Saleh ,N. S. Subotic . Time2varying filtering of signals in the mixed 8 () time2frequency domain J . IEEE Trans. On ASSP ,1985 ,33 3:1479 - 1485 . E. P. Wigner . On the quantum correction for thermodynamic equilibri2 9 um J . Phys. Rev. ,1932 ,4 :749 - 759 . 图 3 MUSIC 方法对两信号的 50 次 Monte2Carlo 仿真实验结果 10 H. Choi , W. J . Williams. Improved time2frequency representation of 表 1 TF2SSF 方法估计方差与信噪比的关系 multicomponent signals using exponential kernels J . IEEE Trans. - 3dB 0dB 3dB 6dB 10dB SNR ASSP ,1989 ,ASSP237 . 60? 01097 01085 01053 01017 010051 63? 01086 01075 01051 01011 010048 作者简介 : 7 结论 金 梁 1969 年出生 , 教授 ,1999 年于西安本文研究了如何将联合时频分析的方法应用于阵列信号 交通大学获通信与信息系统专业博士学位 ,现工 处理的问题 , 提出了时频子空间拟合 DOA 估计方法. 这是一 作于解放军信息工程大学. 研究方向为阵列信号 种基于信号空时特征结构的多维处理方法 ,它利用二次 Cohen 处理、软件无线电、非平稳信号处理、移动通信 类时频分布构造时频相关矩阵 C, 并用 C代替传统的阵列 x x 等. 曾多次参加国家 863 、国家自然科学基金、国 相关矩阵 R, 通过 C的特征分解实现了信号子空间与噪声 x x 家重点实验室基金等项目的研究 ,已在国内外发 子空间的分离. 与传统方法不同 , 该方法在空域和二维时频域 表论文 40 多篇 ,获省部级一、二、三等奖多项.