CN106646414B - 基于Bi-CGSTAB和SL0算法的MIMO雷达目标参数估计方法 - Google Patents

Abstract

本发明公开了一种基于Bi‑CGSTAB和SL0算法的MIMO雷达目标参数估计方法，属于MIMO雷达目标参数估计技术领域，通过在SL0算法中由离线计算获得的代替感知病态矩阵A的伪逆A^*(AA^*)^‑1，其中(·)^*表示矩阵的共轭转置操作，然后利用改进的SL0算法对MIMO雷达的接收信号y进行处理。避免因MIMO雷达感知病态矩阵病态而导致SL0算法失效，并改善了SL0算法的稳健性，具有较高的重构精度；离线计算并存储MIMO雷达感知病态矩阵的伪逆替代矩阵，在MIMO雷达中利用SL0算法估计目标参数时可以直接调用伪逆替代矩阵的值，节省了病态线性方程组的求解时间，加快了稀疏目标信号的重构速度，提高了MIMO雷达目标参数估计的实时性。

Description

基于Bi-CGSTAB和SL0算法的MIMO雷达目标参数估计方法

技术领域

本发明属于MIMO雷达目标参数估计技术领域，特别涉及一种基于Bi-CGSTAB和SL0算法的MIMO雷达目标参数估计方法。

背景技术

多输入多输出MIMO雷达是在相控阵雷达基础上发展起来的一种新体制雷达体制。该雷达采用了与相控阵雷达相同结构的天线阵列，但是利用波形分集技术扩展了其虚拟阵列孔径，从而提高了目标参数分辨率以及可分辨的最大目标个数，能获得更优的目标检测性能和参数估计性能。

近年来，压缩感知CS是信号处理领域和最优化领域的研究热点，它通过求解基于l_q(0≤q≤1)范数最小化的优化问题，能从少量测量值中以较高概率恢复出稀疏信号。在实际雷达探测场景中，目标个数仅占据少量的分辨单元，MIMO雷达接收到的回波信号可以稀疏表示，因此可以使用CS方法来精确恢复稀疏目标的参数信息。文献[1]利用稀疏约束的正则化迭代最小化SLIM算法从MIMO雷达的少量采样数据中恢复出目标的距离-角度-多普勒信息，然而该方法的目标多普勒分辨能力较差。为了提高目标的多普勒分辨能力，文献[2]提出了一种正则化迭代加权最小化方法RIRMA，该方法在每次迭代中给出加权l_q(0＜q＜1)范数最小化的闭式解来提高算法的运算速度，该方法能准确估计空间稀疏分布目标的角度，距离和多普勒信息。为了加快稀疏信号的重构速度，文献[3]提出了一种基于平滑l₀范数SL0的稀疏信号重构方法，该方法是利用一系列逐步逼近l₀范数的平滑连续的高斯函数，将l₀范数最小化的NP-hard问题转化为易求解的平滑函数的极值问题，能够采用较少的测量值快速重构出稀疏信号，其重构速度比基追踪算法快2～3个数量级，因此被广泛地应用于雷达的目标参数估计问题中。为了提高MIMO雷达成像的实时性，文献[4]将SL0算法应用于MIMO雷达的目标参数估计，该方法利用双曲正切函数来逼近信号的最小l₀范数，并采用修正牛顿法求解该近似l₀范数最小化问题，同时考虑到在实际应用环境中MIMO雷达的感知病态矩阵可能会呈病态，该方法还采用正则化方法来避免因感知病态矩阵导致SL0算法信号重构误差较大的问题，提高了MIMO雷达的目标参数估计速度和性能。然而，该方法的正则化参数只能根据经验来选择，而正则化参数选择不当会导致MIMO雷达目标参数估计性能严重恶化。

当MIMO雷达的感知病态矩阵相邻列之间存在近似的线性相关时，该矩阵的条件数较大，为感知病态矩阵。当利用SL0算法估计MIMO雷达目标参数时，在设定初始值和计算梯度投影中都需要求解病态线性方程组，而MIMO雷达的接收信号中不可避免存在量化造成的误差以及噪声扰动干扰，这些误差扰动会引起病态线性方程组解的剧烈波动，并与真实值相差甚大，从而导致SL0算法失效。

参考文献：

[1]Tan X,Roberts W,Li J,Stoica P.Sparse Learning via IterativeMinimization With Application to MIMO Radar Imaging[J].IEEE Transactions onSignal Processing.2011,59(3):1088-1101.

[2]Gong P,Shao Z.Target estimation by iterative reweightedl_qminimization for MIMO radar[J].Signal Processing,2014,101:35-41.

[3]Mohimani H,Babaie-Zadeh M,Jutten C,A fast approach forovercomplete sparse decomposition based on smoothednorm,IEEE Transactionson Signal Processing.2009,57(1):289-301.

[4]Feng J J,Zhang G,Wen F Q.MIMO radar imaging based on smoothed l0norm[J].Mathematical Problems in Engineering,2015,2015:1-10.

发明内容

发明目的：本发明针对MIMO雷达的感知病态矩阵为感知病态矩阵的问题，提供基于稳定双共轭梯度Bi-CGSTAB和SL0算法的MIMO雷达目标参数估计方法，以解决现有技术中的问题。

技术方案：为实现上述目的，本发明采用的技术方案为：

一种基于Bi-CGSTAB和SL0算法的MIMO雷达目标参数估计方法，所述方法通过在SL0算法中由离线计算获得的矩阵代替感知病态矩阵A的伪逆A^*(AA^*)^-1，其中(·)^*表示矩阵的共轭转置操作，然后利用改进的SL0算法对MIMO雷达的接收信号y进行处理，其具体操作步骤包括：

步骤一，初始化数据：

(1a)加载伪逆替代矩阵将作为算法运行的初始值，设置初始值内循环次数L和步长μ，其中L,μ＞0；

(1b)设置衰减因子ρ，0＜ρ＜1；设置形状参数初值以及形状参数终值σ_J，其中D为在目标探测场景中所划分的距离-角度-多普勒单元总数；

步骤二，迭代求解目标参数，令σ＝σ_j，在可行解集χ＝{β|y＝Aβ}上利用最速上升法求解F_σ(β)的最大值：

(2a)令将(2b)至(2c)步循环L次；

(2b)令其中，β_i(i＝1,2,…,D)为矢量β中第i个元素；

(2c)将投影到可行解集χ＝{β|y＝Aβ}上，即

(2d)令

步骤三，验证σ，当σ＜σ_J时退出，此时得到的为目标场景向量估计值，其中J表示算法退出时j的取值；否则，j＝j+1，σ_j＝ρσ_j-1，返回步骤二；

步骤四，根据非零元素在矢量中的位置计算得出各目标的参数，所述目标场景向量估计值中的非零元素值为各目标的复散射系数值。

进一步，所述伪逆替代矩阵的计算步骤包括：

S1，任意选取满足稀疏条件的目标场景向量α_b，根据感知病态矩阵A生成叠加噪声的虚拟数据y'_b，

S2，根据虚拟数据y'_b和感知病态矩阵A构造线性方程组：A^*Aα'_b＝A^*y'_b，其中α'_b为待求得未知量；(·)^*表示矩阵的共轭转置操作；

S3，利用Bi-CGSTAB算法求解该病态方程组，获得与真实值α_b接近的解

S4，对感知病态矩阵A进行奇异值分解A＝UΣV^*；其中，和分别是获得的左和右奇异向量矩阵，表示复数集；其中，为由奇异值构成的对角矩阵，为大小为M_r(N+P-1)×[PKH-M_r(N+P-1)]的零矩阵，M_r为MIMO雷达接收阵元数，P、K和H分别为在目标探测场景中所划分的距离单元数、角度单元数和多普勒单元数；

S5，定义数据矢量y_u＝U^*y'_b和对角矩阵的估计值为其中，对角元素值α_vi和y_ui分别是数据矢量α_v和y_u中的第i个元素；

S6，计算求得感知病态矩阵A的伪逆替代矩阵并存储，所述矩阵A的伪逆替代矩阵为：其中，矩阵 由未知的元素值构成的对角矩阵。

进一步，所述伪逆替代矩阵的计算在离线的情况下计算求解。

进一步，所述步骤二中高斯函数其中σ为函数形状控制参数β_i表示向量β中的第i个元素。

进一步，所述步骤四中所述目标的参数包括：距离、角度和多普勒。

有益效果：与现有技术相比，本发明具有以下优点：

(1)本发明利用Bi-CGSTAB算法对病态线性方程组进行离线处理来获得其高精度解，根据该高精度解求解出病态感知病态矩阵A的伪逆替代矩阵并在SL0算法实现过程中将感知病态矩阵A的伪逆A^*(AA^*)^-1用矩阵来代替，避免因MIMO雷达感知病态矩阵病态而导致SL0算法失效，并改善了SL0算法的稳健性，具有较高的重构精度。

(2)本发明可以离线计算并存储MIMO雷达感知病态矩阵的伪逆替代矩阵，因此在MIMO雷达中利用SL0算法估计目标参数时可以直接调用伪逆替代矩阵的值，节省了病态线性方程组的求解时间，加快了稀疏目标信号的重构速度，提高了MIMO雷达目标参数估计的实时性。

附图说明

图1是本发明离线计算A的伪逆替代矩阵的流程图；

图2是本发明的MIMO雷达在线进行目标参数估计流程图；

图3是不同方法在SNR＝0dB时的目标距离-角度估计；

图4是不同方法在SNR＝0dB时的目标距离-多普勒估计；

图5不同方法的重构信噪比SER与回波信噪比SNR的变化关系；

图6不同方法的重构均方误差MSE与回波信噪比SNR的变化关系。

具体实施方式

下面结合实施例对本发明作更进一步的说明。

一种基于Bi-CGSTAB和SL0算法的MIMO雷达目标参数估计方法，所述方法通过在SL0算法中由离线计算获得的矩阵代替感知病态矩阵A的伪逆A^*(AA^*)^-1，其中(·)^*表示矩阵的共轭转置操作，然后利用改进的SL0算法对MIMO雷达的接收信号y进行处理，其具体操作步骤包括：

步骤一，初始化数据：

(1a)加载伪逆替代矩阵将作为算法运行的初始值，设置初始值内循环次数L和步长μ，其中L,μ＞0；

(1b)设置衰减因子ρ，0＜ρ＜1；设置形状参数初值以及形状参数终值σ_J，其中D为在目标探测场景中所划分的距离-角度-多普勒单元总数；

步骤二，迭代求解目标参数，令σ＝σ_j，在可行解集χ＝{β|y＝Aβ}上利用最速上升法求解F_σ(β)的最大值：

(2a)令将(2b)至(2c)步循环L次；

(2b)令其中，β_i(i＝1,2,…,D)为矢量β中第i个元素；

(2c)将投影到可行解集χ＝{βy＝Aβ}上，即

(2d)令

步骤三，验证σ，当σ＜σ_J时退出，此时得到的为目标场景向量估计值，其中J表示算法退出时j的取值；否则，j＝j+1，σ_j＝ρσ_j-1，返回步骤二；

步骤四，根据非零元素在矢量中的位置计算得出各目标的参数，所述目标场景向量估计值中的非零元素值为各目标的复散射系数值。

进一步，所述伪逆替代矩阵的计算步骤包括：

S1，任意选取满足稀疏条件的目标场景向量α_b，根据感知病态矩阵A生成叠加噪声的虚拟数据y'_b，

S2，根据虚拟数据y'_b和感知病态矩阵A构造线性方程组：A^*Aα'_b＝A^*y'_b，其中α'_b为待求得未知量；(·)^*表示矩阵的共轭转置操作；

S3，利用Bi-CGSTAB算法求解该病态方程组，获得与真实值α_b接近的解

S4，对感知病态矩阵A进行奇异值分解A＝UΣV^*；其中，和分别是获得的左和右奇异向量矩阵，表示复数集；其中，为由奇异值构成的对角矩阵，为大小为M_r(N+P-1)×[PKH-M_r(N+P-1)]的零矩阵，M_r为MIMO雷达接收阵元数，P、K和H分别为在目标探测场景中所划分的距离单元数、角度单元数和多普勒单元数；

S5，定义数据矢量y_u＝U^*y'_b和对角矩阵的估计值为其中，对角元素值α_vi和y_ui分别是数据矢量α_v和y_u中的第i个元素；

S6，计算求得感知病态矩阵A的伪逆替代矩阵并存储，所述矩阵A的伪逆替代矩阵为：其中，矩阵 由未知的元素值构成的对角矩阵。

前述伪逆替代矩阵的计算在离线的情况下计算求解。

前述步骤二中高斯函数其中σ为函数形状控制参数β_i表示向量β中的第i个元素。

前述步骤四中所述目标的参数包括：距离、角度和多普勒。

结合图1和图2为本发明的雷达目标参数估计方法的计算流程图，其中，图1是本发明离线计算A的伪逆替代矩阵的流程图，图2是本发明的MIMO雷达在线进行目标参数估计流程图。

在SL0算法的设定初始值和投影到可行解集步骤中都需要求解由感知病态矩阵A构成的病态线性方程组，由于测量数据y中不可避免存在量化造成的误差以及噪声扰动等，感知病态矩阵A往往会造成病态线性方程组解的巨大误差。Tikhonov正则化方法和截断奇异值TSVD方法等通常用于求解病态线性方程组问题。由于Tikhonov和TSVD方法中存在矩阵求逆或SVD运算，而MIMO雷达的感知病态矩阵规模较大，导致这些方法耗时较长，因此它们并不适用于改善MIMO雷达中的病态问题。迭代方法通过一系列的迭代解来快速逼近期望解，而且只存在矩阵和矢量的乘法操作，因此非常适用于求解大规模感知病态矩阵以及无特殊结构限制的病态线性方程组问题。Bi-CGSTAB是一种具有速度快、精度高和稳定性好的迭代算法，本发明利用Bi-CGSTAB算法对病态线性方程组进行离线处理来获得其高精度解，并根据该高精度解求解出病态感知病态矩阵A的伪逆替代矩阵在MIMO雷达进行目标参数估计时，可以先加载矩阵并在SL0算法的初始值和投影到可行解集步骤中，用离线计算获得的矩阵代替感知病态矩阵A的伪逆A^*(AA^*)^-1，避免因感知病态矩阵病态而导致算法失效，并改善了算法的稳健性。

本发明的具体实施方式步骤包括：

一、MIMO雷达以向量形式表示的信号模型

假设MIMO雷达发射阵列和接收阵列都为均匀线阵，其阵元数分别为M_t和M_r。发射阵列的发射信号矩阵可表示为

式中，x_m＝[x_m(1),x_m(2),...,x_m(N)]^T为第m个发射天线的发射信号，其中N为发射信号的长度，(·)^T表示转置。

将目标探测场景划分为P个距离单元、K个角度单元和H个多普勒单元，那么在目标探测区域内共有D(D＝P·K·H)个离散的距离-角度-多普勒单元{(τ_p,θ_k,ω_h)}，其中1≤p≤P,1≤k≤K,1≤h≤H，ω_h＝2πf_dhT，f_dh为多普勒频率，T为采样间隔。那么含有ω_h的第m路目标回波为

式中，⊙表示Hadamard积，将式(2)所表示的扩展成N×M_t维的信号矩阵X_d，即

因为目标场景划分成P个距离单元，因此目标回波之间最大距离单元(即第一个和最后一个距离单元反射信号之间最大可能的时延)为P-1，则发射信号矩阵表示为

式中，是补零后维度为(N+p-1)×M_t的发射信号矩阵；是维度为(P-1)×M_t的零矩阵。目标场景划分为K个角度单元，表示为θ_k,k＝1,...,K，那么接收阵列和发射阵列的导向矢量和分别为

式中，d_r和d_t分别表示接收阵列和发射天线阵列的阵元间距；λ₀为载波波长。

MIMO雷达接收回波信号可表示为

式中，(·)^*表示共轭转置；E为加性高斯白噪声矩阵；α_p,k,h(p＝1,...,P,k＝1,...,K,h＝1,...,H)表示D＝PKH个离散的距离-角度-多普勒单元内的目标复散射系数，若所划分的单元内无目标存在则其复散射系数为零；J_p是大小为(N+p-1)×(N+p-1)的矩阵，它表示用于描述从第p个距离单元反射信号时所采用的转移矩阵，即

将回波信号矩阵Y改写成向量形式，即其中vec(·)表示矩阵向量化运算。定义感知病态矩阵A和目标场景向量分别为

式中，因此MIMO雷达的接收信号向量可表示为

y＝Aα+e (11)

式中，e＝vec(E)。由于目标个数仅占据少量的分辨单元，因此目标场景向量α为稀疏信号，那么利用稀疏信号重构方法即可从式(11)中估计出目标场景向量α，根据α中非零元素的位置就能估计出目标的距离、角度和多普勒等参数信息。

二、MIMO雷达感知病态矩阵病态性分析

由式(11)可知，MIMO雷达的接收信号向量采用了基于过完备字典的稀疏表示方式，感知病态矩阵A的各列分别由D个距离-角度-多普勒划分单元的目标回波信息组成。

由式(9)可知，v_p,k,h和v_p,k+1,h是感知病态矩阵A中相邻角度划分单元所对应的列向量，令和假设z_n和z'_n分别是向量z和z'中第n(n＝1,2,…,N+P-1)个元素。那么相邻列向量v_p,k,h和v_p,k+1,h的互相关值为

式中，conj(·)表示复共轭操作，Δθ＝θ_k+1-θ_k为相邻角度划分单元的间隔。满足下式(13)时，能使得R(Δθ)≈0，即使向量v_p,k,h和v_p,k+1,h不相关。

要使得(13)接近于零，那么指数项中相位变化范围至少为一个圆周，即

(M_r-1)d_r|sin(θ_k+Δθ)-sin(θ_k)|/λ₀≥1 (14)

由于

sin(θ_k+Δθ)-sin(θ_k)＝sin(Δθ)cosθ_k-[1-cosΔθ]sinθ_k (15)

角度划分单元的间隔Δθ很小，上式第二项可以忽略，式(15)可简化为

sin(θ_k+Δθ)-sin(θ_k)≈Δθcosθ_k (16)

将式(16)代入式(14)，可得

式(17)给出了要使得相邻列向量v_p,k,h和v_p,k+1,h不相关，角度划分间隔Δθ应满足的条件。假设MIMO雷达的接收阵元间距d_r＝λ₀/2，接收阵元数M_r＝15，第k个角度单元θ_k＝0°时，则角度划分间隔|Δθ|≥8.2°时，感知病态矩阵A中相邻角度划分单元所对应的列向量才满足不相关的条件。然而在实际情况下，为了使得稀疏重构算法具有较高的角度分辨率，那么在构造感知病态矩阵A时，实际划分的角度间隔会远远小于8.2°，这样导致在感知病态矩阵A中不可避免地会存在近似线性相关的列向量，此时矩阵A的条件数较大，成为了感知病态矩阵。同样，当距离和多普勒划分间隔过小时也会加剧感知病态矩阵A的病态性。

三、基于Bi-CGSTAB和SL0算法的MIMO雷达目标参数估计方法

Bi-CGSTAB算法是一种基于双边Lanczos算法和基于残差正交子空间的迭代方法。该算法求解m阶病态方程问题Bx＝b的步骤如下：

(1)给定初始值x⁽⁰⁾，最大迭代次数i_max，相对容许误差ε，计算r⁽⁰⁾＝b-Bx⁽⁰⁾，令i＝1；其中B为m阶方阵；

(2)如果ρ_i-1＝0，算法失效，退出；

(3)当i＝1时，p⁽ⁱ⁾＝r^(i-1)，否则β_i-1＝(ρ_i-1/ρ_i-2)(γ_i-1/η_i-1)；p⁽ⁱ⁾＝r^(i-1)+β_i-1(p^{(i -1)}-η_i-1ψ^(i-1))；

(4)ψ⁽ⁱ⁾＝Bp⁽ⁱ⁾，s＝r^(i-1)-γ_iψ⁽ⁱ⁾；

(5)ε＝||s||₂/||b||₂，如果ε≤ε，x⁽ⁱ⁾＝x^(i-1)+γ_ip⁽ⁱ⁾，算法停止，输出x⁽ⁱ⁾；

(6)t＝Bs；η_i＝(t^Ts)/(t^Tt)；

(7)x⁽ⁱ⁾＝x^(i-1)+γ_ip⁽ⁱ⁾+η_is，r⁽ⁱ⁾＝s-η_it；

(8)如果||r⁽ⁱ⁾||₂/||b||₂≤ε或i≥i_max，输出x⁽ⁱ⁾，算法结束；否则，令i＝i+1，转向步骤(2)。

本发明基于Bi-CGSTAB和SL0算法实现MIMO雷达目标参数估计：

由于目标稀疏地分布在雷达探测场景内，则真实目标场景向量是稀疏的，因此可以利用压缩感知理论来求解如下关于稀疏向量β的重构问题：

式中，||·||₀表示l₀范数,δ为控制误差大小的参数。求解式(18)获得β值即为稀疏目标场景向量α的估计值。由于l₀范数是非连续函数，因此求解式(18)是一个难以处理问题，无法通过一般的求极值方法得到最优解，只能利用穷举非零值的所有可能排列的方式求解。SL0算法通过一类高斯函数来近似l₀范数，从而将l₀范数最小化问题转化为平滑函数的极值问题，然后可利用最速上升法和梯度投影法来求解目标函数的极值。

定义高斯函数其中σ为函数形状控制参数，则该函数满足以下特性：

定义

式中，β_i表示向量β中的第i个元素。当σ较小时可以利用函数D-F_σ(β)来近似l₀范数，即||β||₀≈D-F_σ(β)。当σ→0时，有||β||₀＝D-F_σ(β)。因此，式(18)表示稀疏信号重构问题可以近似等价为

当参数σ较小时，函数f_σ因高度非光滑而导致许多局部极大值出现，不易进行优化求解；而当参数σ较大时，虽然函数f_σ较为光滑，但是稀疏信号重构误差较大。因此，SL0算法采用2个嵌套迭代运算从式(21)中求解出最稀疏表示解：在外循环中，选择作为算法运行的初始值，通过逐步减小σ的方式来避免优化过程中F_σ(β)时陷入局部极大值；在内循环中，针对每个σ值，在可行解集χ＝{β|y＝Aβ}上寻找使得F_σ(β)达到最大值的β值。

因此，SL0算法在初始值和投影到可行解集步骤中需要求解以下两个线性方程组：

式中，Δβ为投影到可行解集χ＝{β|y＝Aβ}上的修正值。由式(23)可知，那么将投影到可行解集χ上的操作为由于MIMO雷达的感知病态矩阵A通常是感知病态矩阵，因此(22)和(23)都是病态线性方程组，常数向量y往往会受到量化误差和噪声扰动的影响，当y中存在甚至是微小的误差扰动时，就会导致方程(22)和(23)的最小二乘解和剧烈波动并与真值相差较大，致使SL0算法的初始值和投影到可行解集的计算精度较低，进而导致SL0算法失效。为了提高SL0算法的稳健性，寻求病态感知病态矩阵A的伪逆A^*(AA^*)^-1的合理替代矩阵，以提高初始值和投影到可行解集的计算精度。

由于目标探测场景划分情况、阵列结构以及发射信号在目标探测前能预先确定，MIMO雷达的感知病态矩阵A是已知的，因此A的伪逆替代矩阵可以离线计算的。任意选取满足稀疏条件的目标场景向量α_b，α_b为大小为PKH×1的列矢量。利用矩阵A和目标场景向量α_b的乘积来构造信号数据y_b＝Aα_b，并在所构造的数据y_b上叠加服从零均值、方差为δ²的复高斯噪声矢量即一般情况下，信号y'_b的信噪比可以选择得比较大，其中信噪比定义为

构造以下线性方程组：

Aα'_b＝y'_b (25)

式中，α'_b为未知量。则式(25)即为病态线性方程组，而且y'_b中存在扰动分量e，方程组(25)的最小二乘解α'_b＝A^*(AA^*)^-1y'_b不稳定，并与真值α_b相差较大，因此采用Bi-CGSTAB算法求解式(25)。由于Bi-CGSTAB算法要求线性方程组的系数矩阵为方阵，因此在式(25)两端乘上A^*，即

A^*Aα'_b＝A^*y'_b (26)

利用Bi-CGSTAB算法求解病态方程组(26)，即可获得与真实值α_b接近的解本发明利用Bi-CGSTAB算法的估计值和信号矢量y'_b来获得病态感知病态矩阵A的伪逆A^*(AA^*)^-1的替代矩阵即满足以下关系式：

对矩阵A进行奇异值分解(Singular Value Decomposition，SVD)，可得

A＝UΣV^* (28)

式中，和分别是获得的左和右奇异向量矩阵，其中为由奇异值构成的对角矩阵，为大小为M_r(N+P-1)×[PKH-M_r(N+P-1)]的零矩阵。那么矩阵A的伪逆可表示为

式中，为了能求解矩阵A的伪逆替代矩阵令

式中，其中由元素值构成的对角矩阵，它是一个待求解的未知矩阵。将式(30)代入式(27)，即可构造一个等式方程组：

由于V是酉矩阵，则V^*V＝I，其中I为单位矩阵。将式(31)两边同乘V^*，可得

令数据矢量y_u＝U^*y'_b和代入式(32)可得

则对角矩阵的对角元素值d_i可以通过下式估计

式中，α_vi和y_ui分别是矢量α_v和y_u中的第i个元素。那么对角矩阵的估计值为其中diag{·}表示对角化操作,因此对角矩阵将对角矩阵代入式(30)可获得矩阵A伪逆的替代矩阵

在SL0算法的初始值和投影到可行解集步骤中，用求解得到的代替A^*(AA^*)^-1，避免因感知病态矩阵病态而导致算法失效，并改善了算法的稳健性。

本发明的技术效果可以通过以下仿真结果进一步说明。为了验证本发明方法在改善MIMO雷达病态问题方面的优势，进行了几组分别利用RIRMA方法、SL0方法、SL0_Tikhonov方法、SL0_TSVD方法和本文方法进行MIMO雷达目标参数估计的对比实验，其中SL0_Tikhonov方法和SL0_TSVD方法分别采用Tikhonov正则化方法和TSVD方法解决SL0算法中的病态问题，这两种方法都采用L-曲线法来确定其正则化参数。

仿真参数设置：MIMO雷达系统的发射阵列阵元个数M_t＝15，接收阵列阵元个数M_r＝15，它们按均匀线阵布置，其中发射天线间距为d_t＝2.5λ₀，接收天线间距为d_r＝0.5λ₀；发射阵列各阵元发射相互正交的Hadamard编码信号，发射波形的采样个数N＝32。目标场景的距离单元数P＝12；雷达扫描的角度范围为[-30°,30°]，角度划分间隔为1°，则划分后的角度单元数K＝61；目标多普勒频率单位采用度，即Φ_h＝ω_hN(180°/π)，感兴趣的多普勒范围为[-25°,25°]，多普勒角度划分间隔为5°，则划分后的多普勒单元数H＝11。在SL0方法、SL0_Tikhonov方法、SL0_TSVD方法和本文方法中都采用了SL0算法，其中SL0算法的运行参数设置为σ_J＝0.03，ρ＝0.8，内循环次数L＝30，步长因子μ＝2。在RIRMA方法中，选取q＝0.3，迭代次数为6。在本发明方法中，选择相对容许误差为ε＝10^-1，同时设定接收回波虚拟数据的信噪比SNR＝10dB和向量α_b的非零元素个数K_s＝7来计算感知病态矩阵A的伪逆替代矩阵

仿真内容1：MIMO雷达距离-角度-多普勒目标参数估计

图3为MIMO雷达在多普勒5°处的距离-角度像。其中，图3(a)为真实目标的距离-角度分辨单元分布，图3(b)、图3(c)、图3(d)、图3(e)和图3(f)分别为SL0方法、RIRMA方法、SL0_Tikhonov方法、SL0_TSVD方法和本文方法估计获得的目标的距离-角度成像图。图4为MIMO雷达在多普勒单元-10°处的距离-多普勒成像，其中，图4(a)为真实目标的距离-多普勒分辨单元分布，图4(b)、图4(c)、图4(d)、图4(e)和图4(f)分别为SL0方法、RIRMA方法、SL0_Tikhonov方法、SL0_TSVD方法和本文方法估计获得的目标距离-多普勒成像图，其中回波信噪比为0dB。由图3和图4可知，由于感知病态矩阵的病态性导致SL0算法失效，无法估计目标参数，而RIRMA方法的距离-角度和距离-多普勒成像的旁瓣电平较高，不利于目标检测；SL0_Tikhonov方法和SL0_TSVD方法虽然采用Tikhonov正则化方法和TSVD方法解决SL0算法中的病态问题的方法，但是它们的距离-角度和距离-多普勒成像的旁瓣电平要高于本发明方法；本发明方法在SL0算法实现过程中将感知病态矩阵A的伪逆A^*(AA^*)^-1用离线计算出的来代替，避免因MIMO雷达感知病态矩阵病态而导致SL0算法失效的问题，实现目标距离-角度-多普勒的精确估计。

仿真内容2：稀疏目标信号重构性能与回波信噪比的变化关系

图5和图6分别是不同方法的重构信噪比SER和重构均方误差MSE与回波信噪比SNR的变化关系，其中信号重构均方误差MSE定义为其中为真实目标场景向量α的估计值。由于SL0算法因存在病态感知病态矩阵导致算法失效，因此在图5和图6中该算法不参与比较。由图5和图6可知，本发明方法对目标信号的重构性能始终要优于RIRMA方法；虽然SL0_Tikhonov方法和SL0_TSVD方法也采用常用的Tikhonov正则化方法和TSVD方法来求解MIMO雷达的病态问题，但是这些方法的目标信号重构性能要劣于本发明方法。

仿真内容4：不同算法的运行时间对比

虽然CPU运行时间对算法复杂度无法准确测量评估，但是可以粗略评价RIRMA方法、SL0方法、SL0_Tikhonov方法、SL0_TSVD方法和本文方法的运算复杂度。本实验在MATLABR2013a中完成，计算机配置为：Intel(R)Core(TM)i5-M560处理器、主频为2.67GHz、内存为4GB。由于本发明方法可以通过离线计算感知病态矩阵的伪逆替代矩阵因此在目标参数估计之前可以预先计算并存储在从MIMO雷达接收回波信号中估计目标参数时能直接调用该矩阵的值，因此离线计算可不计入该方法的运行时间内。表1为不同算法的运行时间、重构信噪比和重构均方误差。由表1可知，RIRMA方法在每次迭代中存在大规模更新矩阵的求逆运算，导致耗时较长，而本发明方法的运行时间比RIRMA方法约降低了96.7％；在利用SL0算法估计MIMO雷达目标参数时因病态感知病态矩阵导致SL0算法在运行过程中发散，增加了算法运行时间；虽然SL0_Tikhonov方法和SL0_TSVD方法分别利用两种正则化方法即Tikhonov方法和TSVD方法来求解MIMO雷达病态问题，但是这两种方法在计算正则解时利用L-曲线准则来确定与原始回波数据的误差水平相匹配的正则化参数，因此无法离线计算正则化参数，从而这两种方法的运行时间要高于本发明方法；相比于RIRMA方法、SL0方法、SL0_Tikhonov方法和SL0_TSVD方法，本发明方法不仅所需计算时间最少，而且重构性能是最好的。

表1不同算法的运行时间、重构信噪比和重构均方误差

以上所述仅是本发明的优选实施方式，应当指出：对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (4)

1.一种基于Bi-CGSTAB和SL0算法的MIMO雷达目标参数估计方法，其特征在于：所述方法通过在SL0算法中由离线计算获得的伪逆替代矩阵代替感知病态矩阵A的伪逆A^*(AA^*)^-1，其中(·)^*表示矩阵的共轭转置操作，然后利用改进的SL0算法对MIMO雷达的接收信号y进行处理，其具体操作步骤包括：

步骤一，初始化数据：

(1a)加载伪逆替代矩阵将作为算法运行的初始值，设置初始值内循环次数L和步长μ，其中L,μ＞0；

(1b)设置衰减因子ρ，0＜ρ＜1；设置形状参数初值以及形状参数终值σ_J，其中D为在目标探测场景中所划分的距离-角度-多普勒单元总数；

步骤二，迭代求解目标参数，令σ＝σ_j，j＝1,2,3…，在可行解集χ＝{β|y＝Aβ}上利用最速上升法求解F_σ(β)的最大值：

(2a)令将(2b)至(2c)步循环L次；

(2b)令其中，β_i(i＝1,2,…,D)为矢量β中第i个元素；

(2c)将投影到可行解集χ＝{β|y＝Aβ}上，即

(2d)令

步骤三，验证σ，当σ＜σ_J时进行步骤四，此时得到的为目标场景向量估计值，其中J表示算法退出时j的取值；否则，j＝j+1，σ_j＝ρσ_j-1，返回步骤二；

步骤四，根据非零元素在矢量中的位置计算得出各目标的参数，所述目标场景向量估计值中的非零元素值为各目标的复散射系数值；

所述伪逆替代矩阵的计算步骤包括：

S1，任意选取满足稀疏条件的目标场景向量α_b，根据感知病态矩阵A生成叠加噪声的虚拟数据y'_b，

S2，根据虚拟数据y′_b和感知病态矩阵A构造线性方程组：A^*Aα'_b＝A^*y'_b，其中α'_b为待求得未知量；(·)^*表示矩阵的共轭转置操作；

S3，利用Bi-CGSTAB算法求解该病态方程组，获得与真实值α_b接近的解

S4，对感知病态矩阵A进行奇异值分解A＝UΣV^*；其中，和分别是获得的左奇异向量矩阵和右奇异向量矩阵,表示复数集；其中，为由奇异值构成的对角矩阵，为大小为M_r(N+P-1)×[PKH-M_r(N+P-1)]的零矩阵，M_r为MIMO雷达接收阵元数，P、K和H分别为在目标探测场景中所划分的距离单元数、角度单元数和多普勒单元数；

S5，定义数据矢量y_u＝U^*y'_b和对角矩阵的估计值为其中，对角元素值i＝1,2,…M_r(N+P-1)；α_vi和y_ui分别是数据矢量α_v和y_u中的第i个元素；

S6，计算求得感知病态矩阵A的伪逆替代矩阵并存储，所述矩阵A的伪逆替代矩阵为：其中，矩阵 由未知的元素值构成的对角矩阵。

2.根据权利要求1所述的基于Bi-CGSTAB和SL0算法的MIMO雷达目标参数估计方法，其特征在于：所述伪逆替代矩阵的计算在离线的情况下计算求解。

3.根据权利要求1所述的基于Bi-CGSTAB和SL0算法的MIMO雷达目标参数估计方法，其特征在于：所述步骤二中高斯函数其中σ为函数形状控制参数，β_i表示向量β中的第i个元素。

4.根据权利要求1所述的基于Bi-CGSTAB和SL0算法的MIMO雷达目标参数估计方法，其特征在于：所述步骤四中所述目标的参数包括：距离、角度和多普勒。