基于DSP的快速MUSIC测角算法

2024-04-01金芳晓

现代雷达 2024年1期

金芳晓,崔伟

(1. 中国电子科技集团公司第三十八研究所, 安徽合肥 230088)(2. 合肥移瑞通信技术有限公司, 安徽合肥 230000)

0 引言

基于阵列信号处理的MUSIC算法[1],从空间分辨率上突破了瑞利限制,实现了高分辨波达方向(DOA) 估计向超分辨DOA估计的过渡。其核心原理是以信号子空间和噪声子空间的正交性为基础,划分空间来进行DOA估计[2],广泛用于雷达目标识别跟踪、无线电监测、电子对抗、参数估计等领域。但MUSIC算法在超分辨DOA估计时,存在运算复杂度高,不利于硬件实现等问题[3]。

自从MUSIC 算法提出以来,大量学者对算法的实时硬件实现、降低其运算量进行了深入研究。MUSIC算法的运算量主要集中在矩阵特征值分解和空间谱搜索两部分。针对特征值分解部分,文献[4]利用接收数据近似快速估计信号子空间和噪声子空间,简化了特征值分解问题。文献[5]中提出了一种采用多级维纳滤波求解信号特征矢量和信号空间维数,有效地降低了算法运算量。文献[6]为进一步提高DOA估计精度,提出了基于降维子空间的多重信号分类快速算法,利用最大噪声子空间对估计结果修正,提高测角精度。对于MUSIC算法中空间谱搜索部分的简化,文献[7]和文献[8]分别采用阵元降维和有限域搜索等方法研究了空间谱的快速搜索问题。文献[9]空间相移因子构建方向向量,进而利用数字傅里叶变换(DFT)实现MUSIC谱函数的快速计算。这些算法都从MUSIC算法中关键子步骤入手降低算法运算量,但是距离工程实现还有一定的差距,且算法中用于空间谱搜索的所有搜索角度的导向矢量需占用较大的内存资源进行存储,也在一定程度上影响了工程应用。

为此,本文结合TI1843 DSP处理平台,从降低空间谱搜素运算量和减少阵列导向矢量占用内存空间两方面出发,对基于数字信号处理(DSP)的快速MUSIC算法的硬件实现进行研究。在空间谱搜索方面,考虑MUSIC算法硬件工程实现中,在测角范围内存在有多目标,为保证测角角精度和角分辨率,分为粗搜和精搜两步减少谱峰搜索的次数。以角分辨率为步进进行测角范围内的粗搜;在多个目标的方向,以粗搜索的步进作为搜索范围区间,采用黄金分割算法[10]进行精搜,在保证测角精度的同时,大大降低了搜索次数。在空间谱搜索方面,本文利用阵列阵元间的关系,采用迭代的方式构建阵列导向矢量,大大缩短了导向矢量的构建时长,结合上述粗搜+精搜方式进行谱峰搜索,搜索次数大大减少,不再需要提前存储导向矢量,大大减少了内存空间。综上所述,本文算法针对MUSIC算法处理时长和内存资源消耗最大的两大部分,从工程实现上入手,大幅减少了空间谱谱峰的搜索次数,迭代生成导向矢量不再占用内存资源,具有较大的实际应用价值。

1 MUSIC算法介绍

1.1 算法模型

考虑K个远场的窄带信号s(t)=[s1(t),s2(t), …,sK(t)]T,从不同方位角θ1,θ2, …,θK入射到M个阵元组成的均匀线阵,相邻阵元间距为d,则t时刻接收到的数据x(t)=[x1(t),x2(t), …,xM(t)]T为

X(t)=AS(t)+N(t)

(1)

式中:A=[a(θ1),a(θ2), …,a(θK)]为M×K维的阵列流形,也称方向矩阵;a(θk)为M×1维的方向矢量,a(θk)=[1,ej2πdsinθk/λ, …, ej2π(M-1)dsinθk/λ]T,θk表示第k个信号的方向角,λ为波速;n(t)为M×1维的高斯噪声矢量n(t)=[n1(t),n2(t), …,nM(t)]T。

阵列数据的协方差矩阵为

R=E[XXH]=AE[SSH]AH+σ2I=ARsAH+σ2I

(2)

由于信号与噪声相互独立,数据协方差矩阵可分解为与信号、噪声相关的两部分,其中Rs是信号的协方差矩阵,ARsAH是信号部分。

对R进行特征值分解有

(3)

式中:US是由大特征值对应的特征矢量张成的子空间即信号子空间;UN是由小特征值对应的特征矢量张成的子空间即噪声子空间;ΣS为大特征值组成的对角阵;ΣN为小特征值组成的对角阵。理想条件下,数据空间中的信号子空间与噪声子空间是相互正交的,即信号子空间中的导向矢量与噪声子空间正交,如下所示。

aH(θ)UN=0

(4)

经典的MUSIC算法正是基于上述这个性质提出的,但考虑到实际接收数据是有限长的,数据协方差的最大似然估计为

(5)

(6)

MUSIC算法的谱估计公式为

(7)

1.2 算法步骤及运算量

根据上面介绍,整理MUSIC算法的步骤归纳如下。

(2) 对得到协方差矩阵进行特征值分解 (EVD):R=UΣUH。

(3) 变化搜索角度θ,按照式(7)空间谱计算谱函数并进行谱峰搜素,谱峰所对应的角度便是DOA估计值。

进一步,将算法移植到TI1843S DSP的硬件系统上,进行了100次Monte Carlo随机试验。试验条件设定如下:8阵元均匀线阵,单快拍,角度搜索次数100。表1给出了MUSIC算法各个子步骤DOA估计所需要运行时间。从表1可以看出,在三个子步骤中,耗时最多的便是空间谱构建及谱峰搜索(表1中的离散谱估计部分),是影响MUSIC算法实时性的关键因素。为此,本文提出了简化空间谱构建和搜索的方法。

表1 MUSIC算法子步骤运算量评估Tab. 1 Evaluation of sub-step computation of MUSIC algorithm

2 导向矢量搜索和构建的简化

为减少空间谱谱峰搜索次数,同时保证测角精度,本文采用黄金分割算法对空间谱搜索进行改进。该方法可以在满足算法精度要求的前提下,成倍的降低运算量,且工程实现方便。

另一方面,传统MUSIC算法在硬件实现时,一般根据系统测角范围和角精度要求,需要提前构建好空间谱矩阵,存储在系统中,用于测角时候进行各个角度的空间谱谱峰搜索。但是,这种方式需要较大内存,且测角范围越大,角精度越高,所需的内存越大。而且,采用本文提出的黄金分割算法对空间谱进行搜索,并不需要调用测角范围内的所有角度对应的空间谱,所以采用传统遍历的方式进行DOA估计很大程度上浪费了存储空间。为此,本文利用阵元间相位关系迭代构建空间谱,大大较少了空间谱的构建时间。进一步,结合黄金分割算法仅构建黄金分割点所对应角度的空间谱,从而实现在空间谱搜索过程中构建相应角度的空间谱,这样在保证系统实时性的前提和测角精度的前提下,大大地减少了内存资源的浪费。

2.1 黄金分割空间谱搜索算法

为减少搜索次数,本文采用一维黄金分割精搜算法取代传统空间谱搜索。首先简单介绍黄金分割算法原理[11],具体如下。

函数f(x):设该函数在区间(a,c)上有且仅有一个极大值,其中b点位于(a,c)区间的黄金分割点,x点位于(b,c)区间的黄金分割点,如图1所示。最大值的搜索过程如下:区间(a,c)上计算的b和x两点对应的函数值,当f(b)>f(x)时,即满足f(b)>f(a)且f(b)>f(x),则可判定极大值点位于三元点af(c)且f(x)>f(c),则可判定极大值点位于三元点组b

图1 黄金分割函数示意图Fig.1 Schematic diagram of golden section function

本文将上述介绍的黄金分割应用于MUSIC算法谱峰搜索中,在保证搜索精度的基础上,大大地减少了搜索次数。考虑到实际情况,在测角范围内存在有多目标,而黄金分割算法只适用于的搜索范围内有且仅有一个极大值的设定下,那么如何设定搜索的步进并采用黄金分割对多个目标进行逐一搜索,是算法工程应用的关键。为此,在已知信源个数的情况下,将空间谱搜索分为粗搜和细搜两个部分。其中,粗搜是为了检测到多个目标,精搜则是为了保证测角的角精度。具体的采用黄金分割精搜对多目标谱峰搜索的方式如下。

根据工程需要,设定测角搜索范围为[-θ,θ],所需角度分辨率为Δθ,角精度为ε,目标个数为K。阵列的波束宽度的近似计算公式为

采用上述粗搜+精搜的方式,在保证角精度和角分辨率的基础上,减少了搜索范围。传统测角算法,为保证测角的精度,应以角精度ε为搜索步进,对测角范围内全角度搜索。而采用本文算法全角度粗搜是以角分辨率Δθ为步进,由ε≈Δθ/5可知,全角度粗搜次数约为传统全角度搜索的1/5。而精搜部分只是在目标角度局部,以角精度为步进进行搜索,且结合黄金分割方法,进一步减少了局部搜索的次数,与局部遍历搜索相比搜索次数约为其1/3,粗搜和精搜总搜索次数较传统遍历的方式相比大大减少。下面以具体参数为例,对比本文搜索方式和传统遍历的方式搜索次数上减少的程度。设定测角搜索范围为[-30°,30°],所需角度分辨率为1°,角精度为0.2°,目标个数为3。采用传统遍历搜索的方式,为保证角分辨和角精度,搜索步进设为0.2°,搜索次数约为300次。而采用本文搜索方式,第一步粗搜以角分辨率1°为步进,搜索次数为60次;第二步精搜,对于第k个目标,在[θk-1°,θk+1°]的2°范围内以角精度0.2°为步进,采用黄金分割算法进行精搜,搜索次数为4次,4个目标共需12次。由此可得,本文搜索精搜+粗搜一共72次,对比传统遍历的方式,搜索次数大大减少,且根据分析可得,搜索范围越大,角精度要求越高,传统遍历方式搜索次数越多,而本文所提出的搜索方式优势越明显。

2.2 导向矢量构建

下面基于均匀线阵导向矢量的模型,推导导向矢量构建的简化方法。在DSP中复数采用实部和虚部分开表示,故将式(1)中的导向矢量a(θk)=[1,ej2πsinθkd/λ,…,ej2π(M-1)sinθkd/λ],a(θk)拆分为实部和虚部两部分,即

a(θk)=areal(θk)+j·aimag(θk)

(8)

其中,

areal(θk)=[1,cos(2πsinθkd/λ,…,cos(2π(M-1)

sinθkd/λ)]

aimag(θk)=[0,sin(2πsinθkd/λ),…,sin(2π(M-1)

sinθkd/λ)]

便可利用areal(θk)和aimag(θk)来构建导向矢量a(θk)。此方式虽然简单直接,但是其运算量却较为复杂。以M个阵元为例,构建θk方向的导向矢量a(θk)需进行4(M-1)次三角函数运算。为此,常规方法一般为提前存储好所需的导向矢量,方便角度搜索直接调用,且角精度要求越高,所需的存储空间越大。

为简化构建,减少存储空间,基于导向矢量a(θk)的数学表达式,本文利用三角函数和差化积和推导阵元间导向元素的递进关系,大大减少运算量。具体推导如下:

取areal(θk)和aimag(θk)中的任意元素,令

qreal(m,θk)=cos(2π(m-1)sinθkd/λ),

m=2,3,…,M-1

(9)

qimag(m,θk)=sin(2π(m-1)sinθkd/λ),

m=2,3,…,M-1

(10)

利用和差化积,可得

qreal(m,θk)=cos(2π(m-2)sinθkd/λ+2πsinθkd/λ)=

cos(2π(m-2)sinθkd/λ)cos(2πsinθkd/λ)-

sin(2π(m-2)sinθkd/λ)sin(2πsinθkd/λ)=

qreal(m-1,θk)qreal(1,θk)-qimag(m-1,θk)·

qimag(1,θk)

(11)

进一步推导,可得

qimag(m,θk)=sin(2π(m-2)sinθkd/λ+2πsinθkd/λ=

sin(2π(m-2)sinθkd/λ)cos(2πsinθkd/λ)+

cos(2π(m-2)sinθkd/λ)sin(2πsinθkd/λ)=

qimag(m-1,θk)qreal(1,θk)+qreal(m-1,θk)·

qimag(1,θk)

基于此可以看出,只需求解qreal(2,θk)和qimag(2,θk),后面m-2个导向矢量元素qreal(m,θk)和qimag(m,θk),便可利用式(11)和式(12)迭代得到。采用此方法,仅需要4次三角函数运算来构建qreal(2,θk)和qimag(2,θk),大大地降低了导向矢量构建的运算复杂度。

采用此方式迭代生成对应角度的导向矢量,计算量小。同时,配合前文提出的粗搜+精搜的搜索方式,空间谱无需按角精度为步进遍历搜索,搜索次数减少。为此,采用本文方法空间谱搜索的次数少,且每次搜索导向矢量构建的运算量少,故每次搜索可根据所需的角度生成相应的导向矢量即可,无需提前存储导向矢量,大大减少了存储空间,具有较强的实际应用价值。

3 试验结果

为了验证本文提出的简化的导向矢量构建和空间谱搜索算法,在实际DSP中运用的有效性,对传统方式和本文方式从空间谱搜索、导向矢量构建时长和所占内存大小三个方面进行工程化对比分析。试验的硬件平台采用TI1843 DSP环境,其主频为600 MHz进行了100次Monte Carlo随机试验。试验测试条件采用8阵元均匀线阵,目标个数为2,为满足角分辨率要求,设定两个目标的角度分别为5°和10°。针对导向矢量构建时长和内存占用情况以及空间谱搜索时长的验证分析,设定如下试验:一、导向矢量构建时长和内存占用情况试验,设定100个搜索角度,分别采用传统和本文两种方式构建100个角度的导向矢量,并统计总的构建时长和内存占用情况进行对比;二、分析对比多目标下空间谱搜索时长,对角精度(也就是搜索步进)分别为1°和0.1°两种情况,采用传统遍历搜索和本文粗搜+精搜两种方式进行空间谱搜索的时长对比。

试验对比结果见表2,由表2可以看出:一、采用本文所提出的递进的导向矢量构建方式,在构造100个方向角度导向矢量的相同条件下,与传统构建方式相比构建时间大大缩短。进一步分析内存资源的占用情况,由于本文采用阵元间迭代构建导向矢量的方式,构建时间大大缩短,且黄金分割谱峰搜索也不需遍历所有角度,故可采用搜索时生成对应角度的空间谱的方式进行角度估计,无需再设置内存用来存储100个搜索角度8阵元的导向矢量,大大地减少了内存空间的占用。二、空间谱搜索时长方面,传统的方式是以角精度为步进遍历所有角度进行搜索,由表2可以看出,采用传统遍历方式,以0.1°为角精度与1°为角精度相比时长约为其10倍,而采用本文粗搜+精搜的方式,时长上大幅缩短,且随着角精度的提高,黄金分割谱峰搜索的优势越明显。