第3章空时谱估计

信号处理的一项基本功能就是信号的参数测量，其中在时域参数中除到达时间外，频率
的测量也是其中一个很重要的内容，而空域参数中的角度则是重要参数。现代信号处理中谱估计通常分为经典谱估计、现代谱估计及空间谱估计，其中前两部分在本章中称为时域谱估计，其主要应用是频率参数的测量； 空间谱估计也称空域谱估计，其主要应用是角度参数的测量。本章从时域谱估计中的频率参数测量、空域谱估计中的角度参数测量入手，分析测频方法与测角方法的算法原理与结构，从而掌握方法中蕴含的空时等效性； 然后介绍空时谱中的频率和角度联合测量方法。整个介绍过程是从经典的算法开始，过渡到现代谱的高分辨算法，再介绍超分辨方法； 最后从空时等效的角度分析空时谱估计与时域谱估计和空域谱估计的关系，并对典型的测频和测角方法进行了小结。
3.1数据模型

◆

由文献［2］可知，假设有一由M个阵元组成的等距均匀线阵，入射信号有Ｎ个，则理想情况下的空域阵列模型通常由下式来表达： 


X(t)=As(θ)S(t)+N(t)(3.1)



其中，X(t)为阵列的M维快拍数据矢量； N(t)为阵列的M维噪声数据矢量； S(t)为空间信号的N维矢量； As(θ)为空间阵列的M×N流型矩阵： 



As(θ)=［as(θ1)as(θ2)…as(θN)］(3.2)


其中，空域导向矢量


as(θi)=1

exp(-j2πdsin(θi)/λ)



exp(-j2π(M-1)dsin(θi)/λ)，
i=1，2，…，N(3.3)


式中，d为等距均匀线阵的间距； λ为波长。
需要注意的是，式(3.1)的阵列模型是基于以下几个假设条件推出的。
(1) 入射的信号为窄带远场信号。这个假设保证了信号入射到各阵元上是近似的平行入射，而且入射的信号时延满足下式： 


si(t-τ)≈si(t)e-j2πfτ，i=1，2，…，N(3.4)


(2) 信号间相互独立，且噪声模型均为独立的高斯白噪声。这个假设保证了阵列接收数据的二阶统计特性满足理想情况。




(3) 阵元数M大于信号源数N。这个假设保证了模型在参数求解过程中存在唯一解。

(4) 各阵元不存在各种误差。当有误差时则需要对这个模型进行相应的修正，详细内容见有关文献，这里不作讨论。

在时域处理中，通常在窄带信号的假设下式(3.4)也是适用的，即如果对一个时间序列s(t)采样，假设采样L个点，采样间隔为τ，则得到的时间序列可以写成如下列矢量： 


s(t)=s(t)

s(t-τ)



s(t-(L-1)τ)=1

e-j2πfτ



e-j2πf(L-1)τs(t)=at(f)s(t)(3.5)


式中，at(f)定义为时域导向矢量。


如果时域采样信号包含Ｎ个信号和噪声，则也可以写成


x(t)=s1(t)+s2(t)+…+sN(t)+n(t)

s1(t-τ)+s2(t-τ)+…+sN(t-τ)+n(t-τ)



s1(t-(L-1)τ)+…+sN(t-(L-1)τ)+n(t-(L-1)τ)

=11…1

e-j2πf1τe-j2πf2τ…e-j2πfNτ



e-j2πf1(L-1)τe-j2πf2(L-1)τ…e-j2πfN(L-1)τs1(t)

s2(t)



sN(t)+n1(t)

n2(t)



nL(t)

=At(f)S(t)+N(t)(3.6)


式中，τ为时域采样间隔； At(f)为时域采样流型，由时域导向矢量构成，其表达式为


at(fi)=1

e-j2πfiτ



e-j2πfi(L-1)τ，i=1，2，…，N(3.7)


通过上述关于模型的分析可知： 

(1)  空域模型和时域模型完全等价，不同的是空域采样数据是按阵元序号排列，而时域采样数据是按采样的序号排列。
(2) 由式(3.1)和式(3.6)可知，空域的理想模型和时域采样的模型是相似的，只是空域阵列流型As(θ)由空域导向矢量as(θ)组成，而时域采样流型At(f)由时域导向矢量at(f)组成。
(3) 由式(3.3)和式(3.5)可知，空域导向矢量只与阵元间隔和角度参数θ有关，而时域导向矢量只与它的采样间隔和频率参数f有关。
(4) 对于均匀采样而言，时域采样流型At(f)是范德蒙德矩阵，时域导向矢量at(f)是等比数列。这种结构和等距均匀阵列的阵列流型和导向矢量相同。
3.2时域谱估计

◆

频率测量是信号处理领域中的一个重要的研究方向，传统的频率测量分为两大类： 一是经典谱的方法，以周期图和自相关法为代表，它们的基础都是离散傅里叶变换； 二是现代谱的方法，以自回归(auto regressive,AR)谱、最大熵法等方法为代表，它们的基础是最大信噪比准则、最小均方准则、最小方差准则及最大似然函数准则等四大准则，后续又出现了子空间类方法、高阶谱和最大似然等方法。本节主要从时域导向矢量运用的角度来介绍常用的频率测量方法，从传统的谱估计、高分辨谱估计和超分辨谱估计入手介绍典型算法。
3.2.1离散傅里叶变换

频率测量也叫频谱分析，它是数字信号处理的重要内容，也是现代谱分析的重要基础。最常用的频率测量方法就是离散傅里叶变换(DFT)，实际应用中常用快速傅里叶变换(FFT)来实现。假设对一个K点序列x(0)，x(1)，…，x(K-1)，
进行K点的DFT变换，则有


X(k)=∑K-1n=0x(n)e-j2πnk/K，k=0，1，2，…，K-1(3.8)


很显然，上式可以转化成


X(k)=1

e-j2π1·k/K



e-j2π(K-1)·k/KTx(0)

x(1)



x(K-1)=aHt(f)x(3.9)


式中，［·］T表示转置运算； ［·］H表示共轭转置运算； 矢量x就是K点采样序列构成的列矢量，即DFT就是对应fk的时域导向矢量与采样序列的内积。所以，序列
进行DFT后得到的矢量为


X(0)

X(1)



X(K-1)=11…1

e-j2π1·0Ke-j2π1·1K…e-j2π1·K-1K



e-j2π(K-1)·0Ke-j2π(K-1)·1K…e-j2π(K-1)·K-1KT
x(0)

x(1)



x(K-1)

=［at(f0)at(f1)…at(fK-1)］Hx=WHDFTx(3.10)


式中，WDFT定义为DFT的权矩阵，它由各对应频率点的时域导向矢量构成。

需要注意的是，上式中的时域导向矢量中的fi是归一化的频率(对采样频率的归一化)，即式(3.7)中


fiτ=fifts=kK(3.11)


上式中fts为采样频率。这也表明式(3.10)中各列的时域导向矢量中的频率是对归一化频率的K等分。如果K点序列做了K′点的DFT，则意味着式(3.10)中时域导向矢量构成的矩阵
维数为K×K′的，即表示此时是对归一化的频率进行了K′等分。

在雷达信号处理中，如果DFT是对相干脉冲的离散傅里叶变换，则脉冲间的采样间隔fts就是脉冲重复频率的倒数，此时的WDFT就对应动目标检测(MTD)的滤波器组，每个滤波器由时域导向矢量构成。图3.1给出了8点MTD横向滤波器幅频响应图。


图3.1MTD横向滤波器幅频响应图


从图3.1中可以看出，这8个滤波器的每一个都有固定频率点，权系数对应式(3.10)中的一个a(fi)，而且正好是对归一化频率进行了8等分(-3dB交叠处)
。这也说明采用DFT方法的频率分辨率和精度都和DFT的点数密切相关，分辨率为1/K，精度为1/(2K)。这里1/K也称DFT的傅里叶限，即在一个频率分辨单元内无法分辨两个及以上的目标。另外，由图3.1可知，每个滤波器的第一旁瓣电平还是很高的(理想情况下都是-13.2dB)，此时可以通过幅度加权(如汉宁窗、海明窗等)的方法来降低旁瓣，但会导致主瓣宽度变宽且增益下降。
由于DFT每个时域的导向矢量导致了其频率指向也是固定的，但实际接收的信号频率可能和DFT划分的指向不一致，可以用改变频率指向的方法来解决，此时输出数据为


y(t)=WHx(t)=aHt(f)x(t)(3.12)


所以，整体的输出功率为


P=E［y(t)yH(t)］=aHt(f)E［x(t)xH(t)］at(f)

=aHt(f)Rxxat(f)(3.13)


式中，E［·］表示统计平均； Rxx为数据矢量的自协方差矩阵。这样就可以得到常规的波束扫描(CBF)算法


PCBF(f)=aHt(f)Rxxat(f)(3.14)


对于上述算法的说明： 
(1) 利用DFT进行频率测量实质是CBF算法的特例，式(3.14)中的变量为f，其变量范围为归一化的频率［0，1］，而DFT则固定为k/K，k的取值为 k=0，1，2，…，K-1； 
(2) CBF算法的输出是接收数据加权后的功率，而且其就是频域的匹配滤波，所以它是最大信噪比准则下的频率测量方法； 
(3) CBF算法通常是采用搜索来进行求解，且步长为自己定义，最大值对应的频率就是测量值，所以其运算量肯定大于用FFT求解的DFT算法。

图3.2给出了DFT测频的仿真分析，仿真中有两个信号，其信号频率分别为0.3MHz和1.5MHz，采样频率为2MHz，采样点数为100点，信噪比为20dB。


图3.2不同点数的DFT频率测量图



从图3.2中可见： ①DFT是可以进行频率测量的，但频率的分辨率和精度与DFT的点数密切相关； ②DFT的点数越多，则相当于对测量的不模糊范围分得越细，精度越高； ③DFT的测量范围为
［-fs/2，fs/2］，所以1.5MHz对应的频率被模糊到了-0.5MHz处，如果把图按0
~2MHz取范围，则-1~0MHz处被平移到1~2MHz处，此时就会在1.5MHz附近显示一个目标。

图3.3给出了采用式(3.14)所示的CBF算法来进行频率测量的仿真结果，其他实验条件同图3.2，只是采用频域的CBF算法，搜索步长为20kHz，搜索范围为0~2MHz。


图3.3CBF算法频率测量


对比图3.3和图3.2可知： ①由于CBF算法采用的是频率搜索，虽然其计算量大了很多，但其估计精度比较高； ②降低CBF算法频率搜索的运算量通常可以考虑采用二级搜索，第一级通过粗搜索或者直接利用DFT结果，第二级再在粗估计结果附近进行精细的搜索； ③DFT为了利用FFT算法，通常采用的变换点数为2n，而CBF则比较自由，可以是任意数； ④DFT是CBF的特例，即相当于at(f)中的频率为固定频率的搜索； ⑤由于CBF是频域的匹配滤波，
则如果把式(3.12)看成是频域的加权输出，式(3.7)中的at(f)就是最大信噪比准则下的最优权。


CBF和DFT方法属于经典谱估计范畴，由文献［20］可知其
通常采用自相关法和周期图法，但这类方法存在难以克服的固有缺点： 
①频率分辨率不高，这是因为它们的频率分辨率反比于数据记录长度，而实际应用中通常不可能获得很长的时间序列； ②经典谱估计都是以DFT为基础的，它隐含着对数据序列的加窗处理，而且加的是一个矩形窗，带来的问题是旁瓣太高，会淹没弱信号； ③可以通过幅度加权(汉宁窗、海明窗、切比雪夫窗等)的方法来改善旁瓣问题，但加权会引起主瓣的展宽和增益的下降，导致信噪比出现损失，损失量可参见文献［5］。

为了改进经典谱的相关问题，广大学者经过不懈努力，提出了很多性能优异的现代谱估计方法，下面就从信号处理准则的角度介绍几种现代谱估计方法。
3.2.2最小均方误差算法

最小均方误差算法是自适应滤波的一种基本算法，它的结构如图3.4所示。由图可知，这种结构就是一个AR滤波器。



图3.4最小均方误差自适应原理图


对于一个输入序列u=［u(n)，u(n-1)，…，u(n-M+1)，u(n-M)］T，假设u(n)是要得到的期望数据，其他是输入数据，很明显这些输入数据的加权求和就是u(n)的估计值： 



u^(n)=∑Mk=1w*ku(n-k)=WHfuf(3.15)


式中，数据矢量uf=［u(n-1)，u(n-2)，…，u(n-M)］T； 自适应权矢量Wf=［w1，w2，…，wM］T。

显然估计值与真实值之间存在一个误差


e(n)=u(n)-u^(n)=u(n)-WHfuf(3.16)



则最小均方误差准则： 


ε(n)=E［e2(n)］=min(3.17)


将式(3.15)和式(3.16)代入上式可得


ε(n)=E［e2(n)］=E［u2(n)］+WHfE［ufuHf］Wf-2E［ufuH(n)］Wf(3.18)



定义协方差矩阵和互相关矢量分别为


Rf=E［ufuHf］，
rf=E［ufuH(n)］(3.19)


则式(3.18)可以简化成


ε(n)=E［u2(n)］+WHfRfWf-2rfWf(3.20)


将上式对Wf求导，并令其等于０，即得


ε(n)Wf=2RfWf-2rf=0(3.21)


式中，协方差矩阵如果可以直接求逆(通常Rf是正定的且对角线绝对占优的
矩阵，可以直接求逆)，则得到最优的最小均方误差权矢量


Wf=R-1frf(3.22)


上述的推导过程说明对于自适应线性组合器整体的权矢量为


W=1

-Wf=1

-R-1frf(3.23)


此时，得到的AR谱就是


PAR(f)=1‖aHt(f)W‖2=1aHt(f)WWHat(f)(3.24)


经典的最大熵算法则可以由下式的最优化问题求解： 



minWWHRW



uT0W=1
(3.25)


式中，u0=［1，0，…，0］T； R=E［uuH］。利用拉格朗日常数法，就可以求出上式定义的目标函数最优权


W=μR-1u0(3.26)


式中，μ=1/（R-1)11，即为协方差逆矩阵的第1行第1列元素的倒数。
因为常数不影响谱的形状，忽略常数就得到了Burg提出的最大熵算法(MEM)



PMEM(f)=1‖aHt(f)R-1u0‖2(3.27)



对比一下式(3.25)和式(3.22)，可以发现两者存在一定的差别，因为u和uf的维数差一个数u(n)，所以R和Rf之间具有如下关系： 


R=E［uuH］=
u(n)u*(n)rHf

rfRf(3.28)


很显然，如果上式中的R是Toeplitz矩阵(一般情况都满足)，则其逆矩阵可以写成如下形式： 


R-1=d11dH

dD(3.29)


式中，d11表示R-1的第1行第1列元素； d为逆矩阵第1列去除第1个元素的矢量，此时下式成立： 


dd11=-R-1frf=-Wf(3.30)


所以，最大熵的权为


W=μR-1u0=
1

dd11=1

-Wf(3.31)



上述的分析说明： 在一定的情况下，AR算法与MEM算法是完全等价的。另外，从权的形式来看，AR算法的权矢量第1个元素为1，而MEM算法的约束条件就是权矢量的第1个元素为1，这也恰好说明了AR算法和MEM算法之间的等价性。

上面推导了AR算法、MEM算法及两者之间的关系，仔细观察式(3.15)还可以发现，
u^(n)的估计值也可以看成是利用以前的数据序列u(n-M)，…，u(n-2)，u(n-1)来预测当前的数据u(n)，也就是说，可以进一步外推来预测将来的数据。所以，按最小均方误差推导出来的u^(n)也是采用前向线性预测算法(FLP)，其权矢量就是式(3.23)，其谱为


PFLP(f)=1‖aHt(f)W‖2=1aHt(f)WWHat(f)(3.32)



这里只是从最小均方误差的角度讨论频率估计问题，由上述的讨论可以看出： AR谱算法、MEM算法和FLP算法是完全等效的，它们均可以看成是最小均方误差意义下的谱估计算法。
图3.5对比了CBF算法和MEM算法的频率估计谱，仿真中设置了两个信号，采样频率为2MHz，采样点数
为16点，信噪比为20dB，搜索步长为1kHz。其中图3.6(a)中两个信号的频率分别为0.5MHz和0.8MHz，图
3.6(b)中两个信号的频率分别为0.5MHz和0.55MHz。


图3.5CBF和MEM算法频谱图


(a) 0.5MHz和0.8MHz； (b) 0.5MHz和0.55MHz




图3.6AR、FLP和MEM算法频谱图



从图3.5(a)中可知，当两个信号的频率差为0.3MHz时，超过了频率分辨单元0.125MHz(采样频率的十六分之一)，所以两种方法均实现了两个信号的估计，很明显MEM算法的分辨力要强于CBF算法； 由于采样点数只有16点，所以CBF算法的分辨宽度要远远大于MEM算法。从图3.5(b)中可以看出，当两个信号的频率差只有0.05MHz时，小于频率分辨单元，CBF算法(DFT类的算法)无法实现一个频率分辨单元内多信号的估计，而MEM算法还能高精度地分辨两个信号。

图3.6的仿真条件同图3.5，信号为两个，频率分别为0.5MHz和0.55MHz，图中对比了AR、FLP和MEM算法的频谱。虽然，这三种方法分别对应不同的公式，但从图中可以看出，这三种方法的谱线几乎一样，均能分辨在一个频率分辨单元内的两个目标，这再一次验证了本节中的理论分析。

通过本小节的分析可以得出如下结论： ①AR算法、MEM算法和FLP算法是不同专家在不同时期提出的方法，但从最小均方准则的角度看它们是完全等价的，也就是说这三种方法均可得出最小均方意义下的最优权，其由式(3.22)表示； ②AR算法、MEM算法和FLP算法突破了频率分辨单元的限制，实现了频率分辨单元内的多个信号同时高分辨估计； ③CBF算法只能实现频率分辨单元之外信号的分辨，所以其分辨能力是比较低的。对比图3.5和图3.3可知： CBF类的算法提升分辨能力的方法需要在提高采样点数的同时提高DFT的点数。
3.2.3最小方差算法
最小均方误差算法是一种自适应算法，它是从估计值与期望值均方误差最小的角度推出来的。Capon则从自适应滤波的另一个角度提出了最小方差算法，自适应滤波除了图3.4所示的结构外，还可以采用如图3.7所示的结构。


图3.7全自适应结构原理图



仔细对比这两种结构可以发现： ①最小均方误差结构的自适应属于部分自适应，即第一个数据作为期望信号其权值恒等于1，它没有参与自适应权值计算，其他的数据均参与了自适应权值计算，所以这种结构也称为旁瓣对消结构； ②图3.7所示的自适应结构，所有的数据均参与了自适应权值计算，所以这种结构也称为全自适应结构； ③如果将由图3.7计算出来的权值同时除以第1个数据权值w1，则得到的结构和图3.4完全一样，这也说明这两种结构其实是等价的。

由图3.7可知，自适应的输出为


y(t)=∑M-1k=0w*ku(n-k)=WHuf(3.33)