基于分数阶傅里叶变换的线性调频脉冲信号波达方向估计

王瑞，马艳

(西北工业大学 航海学院，陕西 西安710072)

0 引言

宽带线性调频(LFM)信号在声纳、鱼雷自导、通信和海底勘测等探测设备中有着广泛的应用，针对此类信号的波达方向(DOA)估计问题也日益受到人们的重视［1-2］。但是由于宽带LFM 类是一种典型的非平稳信号，现有的非平稳LFM 信号DOA估计算法都是建立在窄带信号的假设基础之上，因此它们不适合于宽带信号。而分数阶傅里叶变换(FRFT)作为传统傅里叶变换(FT)的推广，可以借助于快速傅里叶变换(FFT)来实现，计算更方便。其次是对LFM 信号的能量在最佳旋转角度具有很好的聚集性［3］。因此，将传统的高分辨DOA 估计方法——ESPRIT 和MUSIC 算法与FRFT 相结合在LFM 的DOA 估计领域中引起了广泛重视［4-9］。其中MUSIC 类算法需要依据信号的中心频率构造出FRFT 域的阵列流行向量和相关矩阵［7-8］。这些方法都是以观测时间等于信号时宽为前提的，此时信号的中心频率就是信号在t =0 时刻的频率，基于FRFT 直接估计出的信号频率直接对应于信号的中心频率［10］。

但在主动声纳和雷达中，常用发射脉冲信号来完成对目标距离、速度和DOA 的检测和估计。观测时间内是否包含目标回波信号、包含回波信号的长度等都是未知的，对于这种脉冲LFM 不完全包含于观测时间(包含的脉冲信号时宽不等于观测时间)的情况，由于中心频率不再对应于t=0 时刻的中心频率，因此基于FRFT 方法估计的中心频率将有一定的误差，从而给方位估计带来较大的偏差。对于这个问题，文献［11］从理论上分析了信号时宽不等于信号观测时间时对分数阶谱位置的影响，并进行了修正，但其中要求的信号确切位置在实际应用中是不知道的，因此在实际中不太适用。

文中针对雷达、声纳和鱼雷自导中这种观测时间与信号时宽不相等的LFM 脉冲信号(忽略多普勒频移)，提出了一种新的基于FRFT 的中心频率估计方法，并在此基础上对分数阶傅里叶域的阵列流行向量进行了修正，进而实现了对基于分数阶傅里叶域MUSIC 算法的宽带LFM 信号DOA 估计性能的改善。并分析了信号宽度和观测时间宽度比值对分数阶傅里叶域信噪比及估计性能的影响。该方法提高了基于FRFT 的宽带LFM 脉冲信号的参数估计准确性，并进一步扩大了用FRFT 处理宽带LFM 的适用范围。

1 分数阶傅里叶变换

FRFT 有很多种不同的定义方式，其中一种是从线性积分的角度给出了FRFT 的基本定义。信号x(t)的p 阶FRFT 定义为

式中:α = pπ/2;Fα［·］为FRFT 的算子符号;Kα(t，u)为FRFT 的变换核，

式中:Aα=

在实际工程中，信号都是经过采样后的离散时间信号，通常采用离散FRFT(DFRFT)，Pei 和Ozaktas 等在这方面做出了卓越的贡献［12-16］。Ozaktas等在文献［16］中提出了一种计算DFRFT 的快速算法，这种方法是对连续FRFT 的核函数在FRFT 和时域直接采样得到的，该算法将FRFT 分解为先乘以LFM 信号，然后和LFM 信号卷积，最后再乘以一个LFM 信号。在利用该算法对信号进行离DFRFT 之前必须对信号要进行量纲归一化处理。对离散信号直接用该算法，暗含着取信号的时域宽度为［-Td/2，Td/2］，其中Td为观测时间，而频率范围为［-fs/2，fs/2］，其中fs为采样频率，则时宽带宽积为N=Tdfs，得到的尺度因子S 和归一化宽度U［17］分别为

则u 域的范围为(-U/2，U/2)，采样间隔为Δu =1/U=(Tdfs)-1/2.

2 基于FRFT 的MUSIC 算法

鉴于FRFT 对LFM 信号良好的能量聚焦性能，对于多个宽带LFM 的方位估计提出了用分数阶傅里叶域的MUSIC 算法［7-8］，在分数阶傅里叶域构造相关矩阵和阵列流行向量，实现对多个宽带LFM 的DOA 估计。

2.1 均匀阵列信号模型

在二维平面内M 个阵元以等间隔d 排列在一条直线上，假设信号为远场源，可看成是平行入射，入射波与阵列的法线夹角θ 为DOA，第一个阵元作为参考阵元。假设有L 个宽带LFM 入射到均匀线列阵，第j 个宽带LFM 为sj(t)= Cjexp［jπ(2fjt +Kjt2)］，j =1，2，…，L，式中:Cj是第j 个信号的幅度;fj为第j 个信号的中心频率;Kj为第j 个信号的调频率。

则第i 个阵元上收到信号为

式中:τij=(i-1)dsin(θj)/c 为第i 个阵元相对参考阵元的时延，θj为第j 个LFM 的入射角，c 为海洋中的声速，为1 500 m/s.假设信号相互之间、信号与噪声中间、各个阵元噪声之间都不相关。由(4)式可以看出，LFM 信号的时延仅仅改变了信号的中心频率和初始相位，但调频率保持不变。

2.2 基于FRFT 的MUSIC 算法

对参考阵元的信号为第j 个LFM 信号sj(t)，由FRFT 的定义式(1)式，得sj(t)关于α 角的FRFT 为

式中:T 为信号时间。由(2)式可以看出，在αj=-arccot Kj时具有最好的能量聚集特性。此时sj(αj，显然，sj(αj，u)在uj=fj/csc αj具有最大值:sj(αj，uj)=最大值点的幅值［18］为

信号分数阶傅里叶域幅度的最大值与信号的时间宽度、信号幅度呈正比，而与|sin αj|呈反比。由FRFT 的时移性质［3］可得第i 个阵元上信号的最大值［8］为

式中:

式中:Xαj= ［Xαj1(u)，Xαj2(u)，…，Xαji(u)，…，为M×N 维的阵元接收信号的αj阶FRFT;Sαj=［Sαj1(u)，Sαj2(u)，…，Sαjj(u)，…，SαjL(u)］T为L× N 维信号的αj阶FRFT;Aαj=［Aαj1，Aαj2，…，为M×L 维αj阶分数阶傅里叶域的阵列流行向量;Nαj=［Nαj1(u)，Nαj2(u)，…，NαjM(u)］T为M×N 维噪声信号的αj阶FRFT.

则Xαj的自协方差矩阵为

式中:Us是由大特征值对应的特征向量构成的信号子空间;UN是由小特征值对应的特征向量构成的噪声子空间;Es是由大特征值构成的对角矩阵;EN是由小特征值构成的对角矩阵。

根据MUSIC 算法得到第j 个信号基于FRFT 的MUSIC 空间谱为

对(12)式进行一维搜索，即可得到第j 个脉冲信号的入射角。从(8)式可以看出阵列流行向量不仅和声程差τij有关，还和信号的中心频率fj及最佳变换角度αj= - arccot Kj有关。在阵列结构确定后，声程差仅与信号的入射角有关。

由(8)式可以看出，LFM 的FRFT 方向向量与信号的中心频率fj有关，对于观测时间不等于信号时宽的脉冲信号，由于信号成分不完全包含于观测时间内，因此信号经过FRFT 后估计得到的频率不再是脉冲信号的中心频率，而是对应于0 时刻的频率，这就给传统的基于FRFT 的MUSIC 算法带来估计上的误差，为此文中提出了一种针对脉冲LFM 的FRFT 的基于MUSIC 算法，对估计的中心频率进行修正，重新构造分数阶域的方向向量。

3 基于FRFT 的脉冲信号MUSIC 算法

3.1 脉冲信号的中心频率估计

LFM 的时频分布线如图1所示，信号的时宽为［-T/2，T/2］，中心频率为f0，调频率为K.观测时间为［-Td/2，Td/2］，当T =Td，并且信号的起点和观测时间的起点相同时，信号的中心频率正好是t=0 时刻点的信号频率，可以用信号αj= -arccot Kj阶FRFT 幅度的最大值点的uj估计:fj=ujcsc αj.但是当LFM 脉冲信号沿t 轴平移τ，即信号不完全包含于观测时间［-Td/2，Td/2］内的情况，如图1所示，在观测时间内LFM 时宽为T -τ，此时信号的中心频率变为f0c，而经过FRFT 后估计出的信号频率为t=0 时刻的频率f0b(0 ＜τ ＜T/2)或f0a(τ ＞T/2)，二者已不再相等。另一方面，延迟后LFM 的调频率保持不变，也就是FRFT 幅度出现最大值的αj不变。

由图1中的几何关系得出:

当τ ＞0 时，

式中:当τ ＞T/2 时，fH为在观测时间内时延之后信号的最高频率，fL为在观测时间内时延之后信号的最低频率。如图1所示，虽然t=0 时刻没有信号成分，但是由于p 阶FRFT 相当于信号的Wigner 分布在时间-频率平面上顺时针旋转α=pπ/2，所以其在u 轴对应u0a，根据(13)式能得出t =0 时刻的频率f0a，因此(15)式仍适用。当τ ＜0 时，同理可得

3.2 修正后基于FRFT 的脉冲信号MUSIC 算法

在实际应用中，通常采用DFRFT.文中选用Ozaktas 提出的计算DFRFT 的快速算法［16］。对参考阵元信号sj(t)做DFRFT 得到sj(α，u)，用多分辨方法［19］在(α，u)平面对sj(α，u)谱峰搜索，得到峰值点对应的位置(αj，uj).由于DFRFT 算法在对信号进行DFRFT 之前必须对信号要进行量纲归一化处理，因此uj将变为进行量纲归一化后的值，结合(3)式可得做N 点DFRFT 后LFMsj(t)的中心频率为

考虑到量纲归一化的影响，(8)式变为

由图1可以看出，脉冲信号平移后，信号的中心频率不再对应于t =0 时刻。由(16)式估计出来的频率j不再是第j 个脉冲信号的中心频率fj，变为图1中所示的f0bj.而第j 个脉冲信号中心频率则对应于图1中所示的f0cj.因此需对于(17)式中的中心频率进行修正。

图1 信号的瞬时频率Fig.1 Instantaneous frequency of signal

结合(14)式得修正后的中心频率

将(18)式代入(17)式得

故修正后的方向向量:

算法具体计算步骤如下:

1)对参考阵元上的接收信号进行DFRFT 并进行二维搜索，得到第j 个脉冲信号出现，峰值坐标

2)由(18)式对中心频率fj进行修正，得到修正后的f0cj.由(19)式构造第j 个脉冲信号的方向向量。

3)由(10)式构造第j 个信号的分数阶相关矩阵RαjXX.

4)对RαjXX进行特征分解，得到噪声子空间EN.由(20)式求得基于FRFT 的MUSIC 算法的空间谱P(θ)cor.

5)找出极大值对应的角度θcor就是信号的入射方向。

6)存在多个脉冲信号时，重复步骤1 ～5，分别得到各个信号的DOA 估计。

3.3 脉宽对估计结果的影响

幅度为C 的LFM 延迟τ 后，在观测时间T 中包含的信号宽度为T - τ，那么它的“最佳”阶αj=-acrcotKj阶FRFT 的幅度为一个Sinc 函数，最大值的幅度为

从(21)式可以看出信号的幅度与延迟τ 成反比，τ 越大，观测时间中包含的信号成分越少，幅度就越低，也就是信号的能量越小。而在观测时间内，噪声的宽度始终是T，对于均值为0，方差为σ2的白噪声，它的αj= - acrcotKj阶FRFT 后依然是白噪声，它的均值为0，方差为所以FRFT 域的信噪比将随着延迟的增大而降低，因而将会对估计精度产生影响。

3.4 运算量讨论

与传统的时域MUSIC 算法相比，基于FRFT 的MUSIC 算法增加了第1 步对参考阵元进行DFRFT并搜索其峰值点的位置(αj，uj)，并估计信号成分的中心频率，和对其他阵元完成αj的DFRFT 并构造分数阶相关矩阵RαjXX，这增加了一定的运算量。但是对LFM 来说，可以用多分辨搜索方法来实现第1 步的搜索［19］，并利用基于FFT 的快速DFRFT来完成对各个阵元上信号的分析［16］。长度为N 的信号，每进行一个角度的DFRFT，所需的计算量为O(Nlog2N).搜索出参考阵元的峰值点位置后，其他阵元的信号只需要做这个角度的DFRFT 就可以了。

4 仿真结果

考虑阵元个数为6 的均匀线阵。信号中心频率为10 kHz，带宽为5 kHz，时宽为0.017 s，采样频率为60 kHz，快拍数N 为1 023，入射角为20°，信噪比SNR=20 dB，信号时延为τ =0.008 s.图2给出了参考阵元上接收信号的α 角的DFRFT 仿真结果。得延迟后其中心频率理论值为f0=8.827 0 kHz，由其峰值位置可得未修正的信号中心频率估计值f0b=7.651 2 kHz，修正后的中心频率估计值为f0c=8.825 6 kHz，时延估计值为τ=0.008 s.图3是在-90°～90°范围内以0.01°为间隔得到的未修正中心频率的传统算法和本文算法的DOA 空间谱pθ.由其峰值可得传统算法入射角估计值^θ =23.35°和修正后的入射角估计值^θcor=20.04°.从图3中可以看出，时延后由于信号不完全包含于观测时间使得中心频率偏移，传统基于FRFT 的MUSIC 算法很难准确地估计出入射角。但本文算法大大提高了此类脉冲信号入射角估计的准确性。

图3 信号DOA 估计的空间谱Fig.3 Estimated spatial spectrum of DOA

其他条件同上，在不同信噪比下做100 次Monte Carlo 独立实验，结果如表1所示。从表1可以看出，在延迟较小时(τ =0.008 s)，尽管未修正的传统FRFT 估计DOA 算法，其估计结果的方差与修正后估计的方差相差无几，但是估计的均值和实际的入射角有较大差异;当延迟较大时(τ =0.015 s)，未修正的算法不仅估计的均值偏差较大，而且方差也随之增大许多，但是修正后的估计均值与理论值相差无几，方差增大了。这正如3.3 节中所述，随着观测时间中包含的信号脉宽的减少，FRFT 域的信噪比下降有关。

表1 不同SNR 和时延情况下的估计误差比较Tabl.1 Estimated errors

5 结论

针对雷达、声纳和鱼雷自导中观测时间与信号时宽不相等的LFM 脉冲信号导致对宽带LFM 信号DOA 估计误差偏大的问题，提出了一种新的基于FRFT 的中心频率估计方法，并在此基础上对分数阶傅里叶域的阵列流行向量进行了修正，改善了基于FRFT 的MUSIC 算法的宽带LFM 信号DOA 估计性能，分析了信号宽度和观测时间宽度比值对分数阶傅里叶域信噪比及估计性能的影响。仿真结果表明，修正后的算法较传统的基于FRFT 的MUSIC 算法:1)对脉冲信号中心频率估计比传统算法准确;2)DOA 估计的准确性较高;3)在脉冲信号宽度特别小时，也能够准确估计DOA，而且方差较小。该方法扩大了用FRFT 处理宽带LFM 信号的适用范围。

References)

［1］ Xiang L，Ye Z，Xu X，et al.Direction of arrival estimation for uniform circular array based on fourth-order cumulants in the presence of unknown mutual coupling［J］.IET Microwaves，Antennas＆ Propagation，2008，2(3):281 -287.

［2］ Liu Z M，Huang Z T，Zhou Y Y.Direction-of-arrival estimation of wideband signals via covariance matrix sparse representation［J］.Signal Processing，2011，59(9):4256 -4270.

［3］ Ozaktas H M，Kutay M A，Zalevsky Z.The fractional Fourier transform with applications in optics and signal processing［M］.New York:John Wiley ＆ Sons，2001:117 -183.

［4］ Cui Y，Liu K，Wang J.Direction of arrival estimation of coherent wideband LFM signals in multipath environment［C］∥2010 IEEE 10th International Conference on Signal Processing (ICSP).Beijing:IEEE，2010:58 -61.

［5］ Chong H，Xiaomin Z.DOA estimation of multi-component LFM in complex environment using ESPRIT based on FRFT［C］∥2011 IEEE International Conference on Signal Processing，Communications and Computing (ICSPCC).Xi’an:IEEE，2011:1 -4.

［6］ Jin X，Zhang T，Bai J，et al.DOA estimation of coherent wideband LFM signals based on fractionalFourier transform and virtual array［C］∥2010 3rd International Congress on Image and Signal Processing (ICISP).Yantai:IEEE，2010:4380 -4384.

［7］ 陶然，周云松.基于分数阶傅里叶变换的宽带LFM 信号波达方向估计新算法［J］.北京理工大学学报，2005，25(10):895-899.TAO Ran，ZHOU Yun-song.A novel method for the direction of arrival estimation of wideband linear frequency modulated sources based on fractional Fourier transform［J］.Transactions of Beijing Institute of Technology，2005，25(10):895 -899.(in Chinese)

［8］ 刘小河，王建英，杨美英.基于分数阶傅里叶变换的宽带LFM相干信号的DOA 估计［J］.数据采样与处理，2008，23(5):547 -550.LIU Xiao-he，WANG Jian-ying，YANG Mei-ying.DOA estimation of coherent wideband LFM sources based on fractional Fourier transform［J］.Journal of Data Acquisition ＆ Processing，2008，23(5):547 -550.(in Chinese)

［9］ 张艳红，齐林，穆晓敏，等.基于分数阶傅立叶变换的WLFM信号DOA 估计［J］.信号处理，2005，21(4):57 -60.ZHANG Yan-hong，QI Lin，MU Xiao-min，et al.DOA estimation of WLFM signal based on fractional Fourier transforms［J］.Signal Processing，2005，21(4):57 -60.(in Chinese)

［10］ 马艳，罗美玲.基于分数阶傅里叶变换水下目标距离及速度的联合估计［J］.兵工学报，2011，32(8):1030 -1053.MA Yan，LUO Mei-ling.FRFT-based joint range and radial velocity estimation of underwater target［J］.Acta Armamentarii，2011，32(8):1030 -1053.(in Chinese)

［11］ 吴晓涛.基于分数阶傅里叶变换的信道估计算法研究［D］，哈尔滨:哈尔滨工业大学，2010.WU Xiao-tao.Research on channel estimation algorithm based on fractional Fourier transform［D］.Harbin:Harbin Institute of Technology，2010.(in Chinese)

［12］ Hsue W L，Pei S C.Rational-ordered discrete fractional Fourier transform［C］∥2012 Proceedings of the 20th European Signal Processing Conference (EUSIPCO).Bucharest:IEEE，2012:2124 -2127.

［13］ Pei S C，Liu C L，Lai Y C.The generalized fractional Fourier transform［C］∥2012 IEEE International Conference on Acoustics，Speech and Signal Processing (ICASSP).Kyoto:IEEE，2012:3705 -3708.

［14］ Koç A，Ozaktas H M，Candan C，et al.Digital computation of linear canonical transforms［J］.Signal Processing，2008，54(6):2383 -2394.

［15］ Oktem F S，Ozaktas H M.Linear algebraic analysis of fractional Fourier domain interpolation［C］∥17th Signal Processing and Communications Applications Conference (SIU)，Antalya:IEEE，2009:873 -875.

［16］ Ozaktas H M，Arikan O，Kutay M A，et al.Digital computation of the fractional Fourier transform［J］.IEEE Transactions on Signal Processing，1996，44(9):2141 -2150.

［17］ 赵兴浩，邓兵，陶然.分数阶傅立叶变换数值计算中的量纲归一化［J］.北京理工大学学报，2005，25(4):360 -364.ZHAO Xing-hao，DENG Bing，TAO Ran.Dimensional normalization in the digital computation of the fractional Fourier transform［J］.Transactions of Beijing Institute of Technology，2005，25(4):360 -364.(in Chinese)

［18］ 刘峰，徐会法，陶然，等.分数阶Fourier 域多分量LFM 信号间的分辨研究，中国科学:信息科学，2012，42(2):136 -148.LIU Feng，XU Hui-fa，TAO Ran，et al.Research on resolution among multi-component LFM signals in the fractional Fourier domain［J］.Scientia Sinica Informationis，2012，42(2):136 -148.(in Chinese)

［19］ 马艳，王瑞，陈航.一种FRFT 变换幅度峰值的快速搜索［C］∥2012 鱼雷自导引信技术研讨会论文集.北京:中国造船工程学会水中兵器委员会，2012，9:95 -98.MA Yan，WANG Rui，CHEN hang.A method of transform amplitude quick search in FRFT［C］∥2012 Torpedo Homing Letter Technology Conference Proceedings.Beijing:China Shipbuilding Engineering Society Underwater Weapons Committee，2012，9:95-98.(in Chinese)