基于改进OMP算法的稀疏目标微波关联成像方法

张 瑞, 全英汇, 朱圣棋, 李亚超, 邢孟道

(1. 西安电子科技大学雷达信号处理国家重点实验室, 陕西 西安 710071;2. 西安电子科技大学电子工程学院, 陕西 西安 710071)

0 引 言

由于成像雷达具有全天时、全天候、抗干扰能力强等显著优势,在灾害监测、对地观测等方面具有广泛应用,已逐步成为获取目标精细特性的遥感设备[1]。基于逆合成孔径或合成孔径技术的雷达设备,其成像原理均依赖于目标与雷达间的相对运动关系,通过发射大带宽信号获得距离向的高分辨率,方位向分辨率取决于合成孔径长度[2-4]。然而,在某些凝视成像或前视成像的场合中,由于在方位向上雷达与观测目标不存在相对运动关系或合成孔径阵列与视线方向相同,不能实现方位向成像,难以获得高分辨图像。

近年来,起源于光学关联成像技术[5]的雷达微波关联成像技术受到海内外越来越多学者的关注。雷达发射端通过发射随机调制的探测信号,使其在空间中叠加产生时间—空间二维随机波动的辐射场,该辐射场与目标相互作用后,信号回波在接收端与事先存储的参考辐射场模板进行关联处理,实现对波束内目标信息的解耦与提取[6-10]。

雷达微波关联成像技术不依赖于雷达与目标间的相对运动关系,与传统的高分辨雷达体制形成互补,可工作于前视模式或凝视观测模式,将其应用在静止或准静止平台上可实现对地的长时间连续监测、资源勘查等功能,具有广泛的应用前景。

目前,有关微波关联成像技术的研究主要集中在分辨力理论、空时二维随机辐射辐射源设计方法[11-14]、关联重构算法[15-17]、三维成像方法[18]、运动目标关联成像[19-21]等方面。其中,由于微波关联成像模型与压缩感知(compressed sensing, CS)存在天然一致性,对稀疏目标成像时,能够取得较好的成像效果,因此稀疏重构类成像算法成为研究热点。国防科技大学、西安电子科技大学、中国科学技术大学等单位对微波关联成像技术相继开展研究,并取得一系列研究成果[22-26]。

针对关联重构算法——正交匹配追踪(orthogonal matching pursuit, OMP)算法在每次迭代时均需求解最小二乘解而导致计算量大的问题,本文提出一种改进OMP算法的稀疏目标微波关联成像方法,该方法可获得良好的成像效果并且降低了计算量,有利于算法的实时化处理。本文首先阐述了雷达微波关联成像机理,构建了成像数学模型;然后结合共轭梯度法推导了改进OMP算法公式,给出了算法实现步骤;最后通过仿真对比实验验证了所提算法的有效性。

1 基于随机辐射场的成像机理与信号模型

1.1 成像机理描述

基于随机辐射场的雷达前视微波关联成像技术由经典热光源强度关联成像技术发展而来。经典热光源强度关联成像技术通过使用热光源在特定的成像平面上产生强度随机起伏的光场来探测目标信息,并且将接收到的探测通道和参考通道中随机起伏的光场强度进行关联得到目标物体的像。

雷达微波关联成像技术采用微波作为信号源,通过发射各种带有调制的电磁波信号(如随机相位调制信号),在感兴趣的观测空间中人为地构造空时二维随机辐射场来模拟经典热光源强度关联成像中随机起伏的光场。微波信号源发射随机调制的信号,并且该信号在空间中非相干叠加,形成分别在空间维度上和时间维度上均具有非相干特性的场。该场与观测区域的目标相互作用后,使得回波中携带有目标对不同辐射场分布的调制信息,将这些回波信号与参考随机辐射场进行关联处理,得到探测平面内目标的成像结果。微波关联成像原理示意如图1所示。

图1 微波关联成像原理示意图Fig.1 Schematic diagram of microwave correlation imaging

辐射场可用U(r)来表示,其空间相干性的定义为

(1)

式中:r1,r2表示波束覆盖区域内不同的两点;R(r1,r2,t)表示随机辐射场在时刻t的空间自相关函数。

辐射场U(r)的时间相干性定义为

(2)

式中：t1,t2表示两个不同时刻;R(r,t1,t2)表示随机辐射场在位置r处的时间自相关函数。

理想情况下,空时二维随机辐射场表现为理想的非相关特性,即

(3)

(4)

在波束覆盖范围内,同一时刻t时不同分辨单元r1、r2的电磁辐射场U(r1,t)、U(r2,t)具有非相干差异性分布特征,使得不同分辨单元被具有不同空间分布特征的辐射场所标度;同一分辨单元r不同时刻t1、t2的电磁辐射场U(r,t1)和U(r,t2)也具有非相干差异性分布特征,使得目标在不同时刻受到具有不同时间分布特征的辐射场所照射。

传统的成像雷达如实孔径雷达,由于发射天线发射单一形式的信号波形,其在探测空间中形成的是起伏规律的电磁辐射场,空间中不同位置处的辐射场相干性较强。一方面,分辨率受实孔径天线尺寸的限制,另一方面,位于同一天线波束范围内的不同目标被相干性较强的辐射场所标度,辐射场不能为目标分辨带来额外的有用信息,因此雷达分辨率较低。此外,发射天线通常发射的是具有周期规律的雷达信号,该信号在观测空间中形成的电磁辐射场在空间上的分布并不随时间的变化而变化,导致雷达回波信号在不同时刻也是相干的,增加照射时间导致信噪比(signal to noise ratio, SNR)提高,对雷达分辨率的提高无任何实质性帮助。而微波关联成像的信号源通过发射随机调制的信号,形成具有空间和时间二维非相干特性的随机辐射场,为雷达分辨同一波束内的不同目标提供额外信息。若能充分利用空时二维随机辐射场的非相干特性,再结合相应的成像算法,则可使得目标超分辨成像成为可能。

1.2 雷达微波关联成像数学模型

假设雷达阵列有M个发射阵元和N个接收阵元。发射阵元排列成均匀线阵,且发射独立的随机信号。在目标成像平面建立直角坐标系,并将成像平面离散划分为Q个相同大小的网格,沿X轴方向的网格尺寸为ΔX,沿Y轴方向的网格尺寸为ΔY。假设目标散射点位于成像网格中心,并以此时散射点的位置矢量r表示该网格的位置。

将目标的等效后向散射系数矢量表示为

σ=[σ1,σ2,…,σq]T

(5)

式中:σq表示第q个成像网格中心rq处等效后向散射系数,若该成像网格单元没有散射点,则σq=0。

为简化起见,现考虑具有M个发射阵元和1个接收阵元的多发单收天线结构。接收天线选定为发射天线中的参考天线(即第1个发射天线)。如此,第m(m=1,2,…,M)个发射阵元在第l(l=1,2,…,L)个脉冲内的发射信号为

(6)

(7)

信号经过调制后,每个阵元的发射信号频率在脉冲之间随机捷变。一方面,这种随机调制的信号在空间中非相干叠加,形成具有空间和时间二维非相干特性的随机辐射场,为后续微波关联成像处理提供重要基础。另一方面,这种形式信号的抗干扰能力较强,特别是对基于数字射频存储器(digital radio frequency memory, DRFM) 的转发式干扰抑制方面效果显著[28]。一般而言,对方干扰机转发的是上一个或上几个脉冲周期的信号波形,对于当前脉冲周期的信号来不及截获、存储、复制和转发。而我方雷达接收机以当前脉冲周期的信号作为参考进行脉冲压缩。由于发射信号在每个脉冲重复周期内都是随机捷变的,并且这种变化只由信号发射方掌握,干扰方只能被动转发,因此干扰信号经脉冲压缩时由于失配而被成功抑制,只有真目标才能被正确检测。

为方便公式推导,假设空间中存在1个点目标,其等效后向散射系数记为σ0,当存在多个点目标时,只需在公式上稍作修改即可,不影响公式推导。

在第l个脉冲内,接收阵元接收到的该单个点目标的回波信号为

(8)

(9)

式中：c表示光速;r0表示等效后向散射系数为σ0的点目标在成像平面坐标系中的位置矢量;Rt,m和Rr,1分别表示第m个发射阵元和接收阵元在成像平面坐标系中的位置矢量。

记第l个脉冲内,M个发射阵元的频率分别为

(10)

将第l个脉冲内接收到的目标回波信号进行下变频及脉冲压缩处理后的信号可以写作:

(11)

式中:B为信号带宽。

在成像平面内,围绕目标划分出Q个大小相同的成像网格,网格的位置矢量由其中心点的坐标矢量表示。由于已知各个阵元的位置坐标(xarray(m),yarray(m),zarray(m))和各个成像网格的位置矢量(xq,yq,zq),于是可计算得到第m个阵元到第q个网格的距离为

(12)

因此,由第m个阵元发射的信号经第q个网格(σq=1)反射后,被接收阵元接收而产生的时延可表示为

(13)

这样,第l个脉冲内第q个网格处(即rq处)的参考信号可写为

(14)

因此,雷达接收端接收到来自第l个脉冲内第q个网格的回波信号可表示为

(15)

将成像模型写为矩阵形式:

y=Sσ

(16)

式中：y=[y(t1),y(t2),…,y(tL)]T为回波信号向量;σ=[σ1,σ2,…,σQ]T为待恢复的目标散射系数向量;S为采样矩阵,其具体形式为

(17)

从式(17)可以看出,观测矩阵不仅与雷达发射、接收阵列几何构型有关,还与雷达信号波形、成像网格数、脉冲数有关。此外,由式(16)还可以看出,解雷达微波关联成像方程是一个典型的线性逆问题求解过程,利用方程的解可以恢复目标、重构目标图像。

如果从另外一个角度来解释待恢复目标与空时二维随机辐射场间的关系,可以认为空时二维非相干随机辐射场的差异性分布特征相当于空间波前的随机编码,随机辐射场与目标相互作用后,目标的信息就蕴含在散射回波信号中。这种相互作用等于对感兴趣的目标进行采样或测量。从这个方面来看,式(17)描述的采样矩阵S与稀疏重构理论里的采样矩阵或者观测矩阵类似。然后,利用本地预先存储的采样矩阵和接收到的目标采样数据进行恢复和反演就可以得到原始目标图像。这种关于目标与随机辐射场间关系的解释与CS理论不谋而合。

对于微波关联成像,空时二维完全非相干随机辐射场天然地满足CS中对观测矩阵的要求:限制等距特性。此外,雷达目标回波信号可以看作是少数几个强散射点回波信号叠加的结果[29]。目标的稀疏特性是雷达成像中比较常用的先验信息。通常情况下,所观测的目标本身就具有稀疏性,此时CS中字典矩阵是一个单位阵,这样采样矩阵就与CS中的测量矩阵对应起来。因此,式(16)所示的成像模型就是CS模型。另一方面,即使当前目标不具有稀疏特性,通常利用各种稀疏变换手段(例如小波变换)便可快速得到其稀疏表示。结合本文中的应用场景,目标飞机在空域中可以看作是由K个主要散射点组成的,在构造成像模型时,需要将目标飞机所在成像平面划分成Q个成像网格,通常K≪Q。因此,飞机在空域上满足稀疏条件。稀疏度的定义可以沿用压缩感知理论中对于稀疏度的描述。如向量σ=[σ1,σ2,…,σQ]T中只有少于K(K≪Q)个非零元素,称K为信号的稀疏度。因此,可以将CS相关理论用于微波关联成像处理中重构恢复目标图像。

由于对成像平面进行网格划分时,网格数量多且尺寸小,接收信号相干性较强。因此,求解式(16)是一个欠定问题,且解通常不唯一。为了得到符合要求的唯一解,这里采用稀疏恢复算法求解。

已有文献证明[30],若σ稀疏,求解式(16)可转化为以下优化问题:

(18)

式中：η表示误差允许量；‖·‖k表示向量的lk范数。但是,求解式(18)是一个非确定性多项式难题(nondeterministic polynomially problem， NP)。最常见的一种解决方法是用范数代替l0范数:

(19)

可以证明,若σ是稀疏的且测量矩阵S满足一定条件时,式(18)和式(19)具有相同的解。

2 微波关联成像算法

从前面的分析可以看出,求解式(16)是一个病态的逆问题。导致这个逆问题病态的原因是测量数据信息量不足,实际目标场景中蕴含的信息量远远多于有限的测量数据量,仅仅采用诸如最小二乘法是无法求解或者无法获得稳定解的。因此,考虑利用目标稀疏的先验信息,结合稀疏重构类算法以获得方程稳定的解,从而恢复目标图像。

2.1 OMP算法

考虑采用OMP算法求解上述问题,其基本思想是从字典矩阵中,选择与输入信号相关性最强的原子,形成一个稀疏逼近,同时求出匹配残差,然后继续选择与残差最匹配的原子,反复迭代直到满足循环终止条件为止。需要稀疏恢复的信号可以由迭代中选择相关性最强原子的线性和再加上最后的残差值来表示。如果残差值在可以接受范围内,则信号就是这些原子的线性组合。

对于压缩观测:

y=Φx

(20)

式中:yM×1为观测所得向量;ΦM×N为测量矩阵;xN×1为原信号,其在某个变换域ΨN×N是稀疏的,即

x=Ψθ

(21)

式中:θN×1是K稀疏的,一般有K≪M≪N。

令传感矩阵A为

A=ΦΨ

(22)

则

y=Aθ

(23)

OMP重构算法的具体步骤如下。

输入:传感矩阵A=ΦΨ,观测向量y,信号的稀疏度K。

步骤 1初始化Λ0=∅,A0=∅,r0=y,i=1。

步骤 3令Λi=Λi-1∪{λi},Ai=Ai-1∪aλi。

步骤 4计算y=Aiθi的最小二乘解:

步骤 5更新残差

步骤 6迭代次数i=i+1,并判断i≤K是否满足,若满足则返回步骤2,否则终止迭代。

其中:i表示迭代数目;ri表示残差值;Λi表示第i次迭代索引构成的集合;λi表示第i次迭代中得到的索引;θi为i×1维列向量;aj表示矩阵A的第j列构成的向量;Ai为依照集合Λi选出的矩阵A的列构成的集合。

在利用OMP算法求解式(19)过程中,涉及大量向量与向量相乘、矩阵与向量相乘的运算,而且每迭代一次就要求解一次最小二乘问题,利用直接矩阵求逆法求最小二乘解时,非常耗费硬件空间、时间资源,这降低了算法的实时性,增加了算法的实现难度。因此,优化OMP算法以利于实时处理是非常有必要的。

2.2 改进OMP算法的稀疏恢复算法

观察第2.1节中OMP算法步骤4可知,每一次迭代都要利用直接矩阵求逆法来求解下式的最小二乘解:

y=Aiθi

(24)

式中:Ai为M×i维矩阵,表示按索引Λi选出的矩阵A的列集合;θi为i×1维的待重建信号。式(24)是一个超定方程,这里只能求得其最优解。定义最优解为使残差y-Aiθi的l2范数取得极小值时的解为

(25)

对式(25)关于θi求梯度并令梯度值为0可得

(26)

进一步化简可得

(27)

Gθ=b

(28)

当矩阵G维数较大时,直接矩阵求逆法的计算量会相当大,这将耗费巨大的时间和空间资源。因此,本文采用共轭梯度法来求解式(28),从而避免了矩阵求逆运算,改进了原始的OMP算法。共轭梯度法只涉及少量的矩阵与向量相乘、向量与向量相乘运算,从而大大降低计算数据量,节约计算资源和空间资源,增强整个成像恢复算法的实用性。

共轭梯度法是用来解决无约束凸二次规划问题的一种方法。其基本思想是对于一个目标函数,找到一组方向向量,依次按此方向向量中的方向对迭代点进行更新,使得函数值在该方向上取得最小值。在每一个新的更新方向对迭代点进行更新时,不会影响在之前方向上的更新结果。如果能找到这样一组方向向量,那么就可保证在多次迭代后找到函数的全局极小点。

共轭梯度法是介于最速下降法与牛顿法之间的方法,既克服了最速下降法收敛慢的缺点,又避免了存储和计算牛顿法所需要的二阶导数信息。该算法具有收敛速度快、所需存储空间少以及二次终止等优点,因此在解决实际问题中得到了广泛的应用。

将式(28)转化为如下问题:

(29)

要应用共轭梯度法,最关键的是要确定算法中的两个参量:优化方向p和优化步长ρ,下面给出其推导过程。

(30)

在上式两边同时左乘p(k)TG,得

(31)

解得

(32)

已知迭代点θ(k)和优化方向p(k),利用一维搜索确定优化步长ρk,也就是求解下列问题:

(33)

为了求解式(31),这里记:

φ(ρ)=f(θ(k)+ρp(k))

(34)

[G(θ(k)+ρp(k))+b]Tp(k)=0

(35)

(36)

求解式(36)可获得第k次的优化步长ρk,其表达式为

(37)

利用共轭梯度法求解式(28)的详细步骤如下所示。

算法 1 共轭梯度法输入:对称正定矩阵G,等式右端项b,初始值θ0;输出:方程的解θ;计算残差:r0=b-Gθ0;初始化:p0=r0;repeat for k=0,1,… do ifpk=0then return θ0 else αk=rTkrkpTkGpk; θk+1=θk+αkpk; rk+1=rk-αkGpk; βk=rTk+1rk+1rTkrk; pk+1=rk+1+βkpk; endend

2.3 算法计算复杂度分析与对比

在OMP算法迭代过程中,相比于直接矩阵求逆,采用共轭梯度法求最小二乘解能够使得问题求解复杂度明显降低,大大节省计算时间,提高算法的有效性和可靠性。下面详细分析对比利用直接矩阵求逆法和利用共轭梯度法求最小二乘解的计算复杂度。

在第2.1节OMP重构算法的步骤4中,利用直接矩阵求逆法求得y=Aiθi的最小二乘解为

(38)

通过以上分析可以看出,由于利用直接矩阵求逆来求解最小二乘解时,不仅要频繁进行矩阵与矩阵相乘操作,还要进行矩阵求逆运算,其计算量随信号稀疏度的增加而快速增长。此外,在硬件实现上需要对矩阵、向量频繁地进行存取,耗费了大量的时间资源,不利于硬件电路性能的提升和整个成像稀疏恢复算法的实时化处理。

共轭梯度法只涉及少量的矩阵与向量相乘、向量与向量相乘、标量间运算,数据量小且数据间的一致性较强,这样就节约了计算时间资源,利于算法实现和实时性处理。

3 仿真结果及分析

根据式(14),可将常用的微波关联成像方法分为匹配滤波法、最小二乘法、稀疏重构法。本小节通过数值仿真对比这三种方法的性能。表1给出了具体的仿真参数取值。

表1 参数设置

仿真场景如图2所示,图2(a)是一个由3个散射点组成的简单目标,图2(b)是一个由8个散射点组成的飞机目标,两种目标都是稀疏的,但稀疏度不同。

图2 原始目标Fig.2 Original target

针对图2(a)的目标场景,分别利用匹配滤波法、最小二乘法、传统OMP算法和改进OMP稀疏重构算法进行仿真,其中SNR=5 dB。

由于匹配滤波法是将回波信号y与采样矩阵S进行相关处理,该方法利用了空时二维随机辐射场的非相干特性,因此回波信号中与采样矩阵相关程度较高的列会在图中产生尖峰。但实际中形成的随机辐射场不能满足理想的时空完全非相干特性,这会大大影响成像质量。如图3所示,图中红色圆圈标出的点即为目标散射点,目标勉强可见。

图3 匹配滤波法仿真结果(3个点目标)Fig.3 Simulation result of matched filtering method (three point targets)

图4 最小二乘法仿真结果(3个点目标)Fig.4 Simulation result of least square method (three point targets)

从现有的仿真结果质量上来看,传统的OMP算法和改进的OMP算法在重构结果的质量方面差异性不大，如图5和图6所示，但从算法处理用时、硬件实现难度以及算法实时性上来看,本文所提方法在运算量上具有巨大优势,有利于算法地实现和实时性处理，4种算法进行目标恢复成像的用时如表2所示。

图5 传统OMP稀疏重构结果(3个点目标)Fig.5 Sparse reconstruction result of traditional OMP method (three point targets)

图6 改进OMP稀疏重构结果(3个点目标)Fig.6 Sparse reconstruction result of improved OMP method (three point targets)

表2 4种算法运行时间对比

由于利用传统OMP算法求解成像方程时,需要涉及大量的向量与向量相乘、矩阵与向量相乘运算,而且每迭代一次就要求解一次最小二乘问题,非常耗费硬件空间、时间资源,降低了算法的实时性,增加了算法的实现难度。此外,由于改进的OMP算法增加了稀疏先验信息的约束,所以削弱了逆问题求解过程中的病态性,降低了采样矩阵S病态性造成的成像不稳定,提高了雷达成像质量。

针对图2(b)的飞机目标,分别利用匹配滤波法、最小二乘法、传统OMP算法和改进OMP稀疏重构算法进行仿真。其中SNR=5 dB。

如图7所示,由仿真结果很难分辨出飞机目标。这是由于匹配滤波法是将回波信号y与采样矩阵S进行相关处理,其本质还是利用随机辐射场的时空非相干特性。若采样矩阵S中各列的相关性很强,则经相关处理后,图上会出现很多幅度值较大的点,导致原始目标被其淹没,难以分辨。

图7 匹配滤波法仿真结果(飞机目标)Fig.7 Simulation result of matched filtering method (aircraft target)

最小二乘法的仿真结果如图8所示。由于最小二乘法对噪声和误差极其敏感,且辐射场的时空相干性较强,所以该算法不能正确得到飞机目标图像,成像失败。

图8 最小二乘法仿真结果(飞机目标)Fig.8 Simulation result of least square method (aircraft target)

图9和图10分别给出了利用传统OMP算法和所提的改进OMP稀疏重构算法的仿真结果。两种算法在有噪声的情况下能够准确恢复目标,得到正确的目标图像。但所提算法优化了传统OMP算法以利于实时处理,改进后的OMP算法数据量小且数据间的一致性较强,节约了计算时间资源,有利于算法硬件实现和实时性处理。

图9 传统OMP稀疏重构结果(飞机目标)Fig.9 Sparse reconstruction result of traditional OMP method (aircraft target)

图10 改进OMP稀疏重构结果(飞机目标)Fig.10 Sparse reconstruction result of improved OMP method (aircraft target)

由于实现匹配滤波法只需进行矩阵的相关运算,因此耗时最短。但如前面所述,该方法的性能受众多因素影响,其成像效果差。改进OMP算法耗时与匹配滤波法相当,但不论是在有噪声的情况下还是当采样矩阵S的非相干性不理想时,其成像结果均令人满意,为3种算法中最优异的一种。最小二乘法耗时最长,原因在于该方法计算过程中需要进行大型矩阵与矩阵相乘操作以及矩阵求逆运算,计算量随着网格划分数目、脉冲数、信号稀疏度的增加而迅速增加,因此非常耗时。此外,由于受到噪声和采样矩阵S非相干性不理想等因素的影响,该方法的成像效果非常差。传统OMP算法步骤涉及大量的向量与向量相乘、矩阵与向量相乘运算,而且每迭代一次就要求解一次最小二乘问题,利用直接矩阵求逆法求最小二乘解时,非常耗费硬件空间与时间资源。

4 结 论

针对传统OMP算法在每次迭代时,均需利用直接矩阵求逆法求解目标函数的最小二乘解而导致计算量大的问题,本文提出一种改进OMP算法的稀疏目标微波关联成像方法。首先,阐明了微波关联成像原理并构建了微波关联成像数学模型;然后,利用共轭梯度法对OMP算法中的直接矩阵求逆步骤进行了改进,并分析了改进后算法的计算量;最后,通过仿真实验验证了所提算法的有效性。实验结果表明,与匹配滤波法、最小二乘法、传统OMP算法相比,所提算法能够获得良好的成像效果,降低了算法计算量,有利于算法的实时化处理。