CN103353596B - 基于压缩感知的波束空间域米波雷达测高方法 - Google Patents

Abstract

本发明公开了一种基于压缩感知的波束空间域米波雷达测高方法，涉及低信噪比和快拍数较少环境下的低仰角测高。其实现过程是：从雷达回波中提取目标信号并进行仰角粗测，得到目标信号仰角所在空域Θ；将空域Θ划分为P份，在空域Θ内进行波束形成得到波束变换矩阵B，并对其进行预白化得到白化波束变换矩阵T；将接收数据投影到该白化波束变换矩阵得到波束域量测信号z，并用观测矩阵φ对z进行压缩采样得到观测信号y；利用白化波束变换矩阵T和观测信号y迭代运算得到目标信号估计值；从目标信号估计值提取目标角度，进而得到目标高度。本发明减少了目标信号的采样点数和运算量，有效降低了信号功率谱和空间谱的旁瓣，提高了低信噪比下的测高精度，可用于目标定位。

Description

基于压缩感知的波束空间域米波雷达测高方法

技术领域

本发明属于雷达信号处理技术领域，特别涉及压缩感知和米波雷达测高方法，可用于目标的高分辨波达方向DOA估计和定位。

背景技术

米波雷达测高方法是阵列信号处理中的一个重要研究方向，在雷达、通信、声纳和导航等领域都得到了广泛的应用和研究。国内外雷达界普遍认为，米波雷达具有反隐身能力。但是米波雷达因受波长长、天线尺寸和架高有限等因素的限制，使得天线波束宽度宽、角分辨力低，更重要的是因地、海面反射即所谓“多径”问题使其难以在多径环境下难以测高，故米波雷达的测高问题一直是雷达界尚未很好解决的难题。

为较好地解决米波测高难题，最主要的技术途径是增大天线在高度维的孔径，以减小天线垂直面的波束宽度。而对于低空目标，即使增大天线在高度维的孔径，也无法避开“多径”问题。

目前，针对米波雷达的测高方法主要有以下四种：

1.多频平滑法。这种方法采用频率分集技术，通过改变雷达的工作频率而改变测角误差，多个频率误差的平均就可以减小测角误差。但只是靠多次平滑来减小测角误差，其测角精度仍然较差。当目标距离较近时需要的变频带宽小，效果不错，然而对于远距离目标要想获得良好的估计性能需要变频带宽大于50%，这很难满足。

2.相关高度分析法。这种方法是根据目标的高度特性和多径反射模型，在所有可能的目标高度处模拟雷达接收数据，构造关于目标测量高度和估计高度的相关函数，可以证明只有在目标真实高度处，相关函数的值达到最大，从而实现精确测高和跟踪，该方法容易出现测高模糊，很可能在多个高度上相关函数都取得最大值。

3.基于波瓣分裂的米波雷达测高方法。这种方法利用不同天线分裂波瓣的相位关系，确定目标所在仰角区间，对接收信号进行比幅处理提取归一化误差信号，最后根据归一化误差信号和仰角区间查表得到目标的高度。其在地面起伏的均方差不超过1m，信噪比达到16dB，测高精度可以达到距离的1%。陈伯孝等在2007年6月在《电子学报》上发表的论文“基于波瓣分裂的米波雷达测高方法”。这是一种在垂直维只需3根天线的米波雷达的低仰角测高方法。该方法只适合于平坦阵地，且对阵地的平坦性要求较高，测高精度也只能达到距离的1%，难以满足一些精度较高的实际使用要求。

4.阵列超分辨处理测高方法。这种方法把阵列信号处理中的超分辨技术应用于分辨直达波信号和多径信号。目前使用超分辨算法进行米波阵列雷达低仰角测高主要集中在两个方面：一是特征子空间类算法；二是最大似然及子空间拟合算法。其中：特征子空间类算法，应用于低仰角测高时必须面对多径传播所引起的直达波和多径信号的相干问题。但当信号源完全相干时，阵列接收的数据协方差矩阵的秩将为1，这会导致信号子空间的维数小于信号源数，即相干源的存在使信号子空间与噪声子空间相互渗透，因而导致某些相干源的导向矢量与噪声子空间不完全正交，使得很多经典特征子空间类算法的性能下降，甚至完全失效。

最大似然类算法思想简单，性能优越，在高信噪比和低信噪比下均有很好的性能，但是似然函数求解是一个非线性的多维最优化问题，需要进行多维网格搜索，计算量随着目标个数的增加呈指数增长，实现过程较复杂。

上述方法中，方法1精度较差；方法2容易出现测高模糊；方法3只适用于平坦阵地，且精度差，无法满足实际需求；方法4运算量大，要求样本数多，在多径环境下性能下降，甚至失效，并且该方法是建立在满足空域奈奎斯特采样定理的条件下，即阵元间距必须小于或等于入射信号的半波长，否则会产生栅瓣，从而导致测角模糊。

发明内容

本发明的目的在于针对上述已有技术的不足，提出一种基于压缩感知的波束空间域米波雷达测高方法，以突破奈奎斯特采样定理的限制，降低运算量，提高低信噪比下波达方向DOA的测高精度。

为实现上述目的，本发明的技术思路是：利用目标在空域的稀疏性，将目标波达方向DOA估计问题看成是一个稀疏向量的重构问题，在波束空间域上进行压缩采样，构建出新的基于压缩感知的波达方向DOA估计模型，并通过稀疏重构算法来进行稀疏向量恢复，用该稀疏向量的非零元素及其在向量中的位置信息分别表征目标的幅度信息和角度信息，从而得到目标的高分辨波达方向DOA估计。具体实现步骤包括如下：

1）从雷达回波中提取目标信号：

x＝AS₀+n

其中，A为M×K维阵列流型，S₀＝[s₁,s₂,…,s_K]^T为K×1维的信号矢量，n为M×1维的白噪声，满足零均值、方差为σ²的复高斯分布，各阵元输出噪声统计独立；

2）根据信源是空域稀疏的特性，采用空间网格划分方法，将空间划分为N份，构造一个信号稀疏化后对应的M×N维阵列流型ψ和N×1维的空域稀疏信号S＝[s₁,s₂,…,s_N]^T，N＞＞K，那么目标信号x改写为：

x＝ψS+n；

3）使用数字波束形成法DBF对目标信号x进行仰角粗测，得到粗侧角度α，进而得到目标信号仰角所在空域O；

4）将所述空域O划分为P份，P＞＞M，M表示阵元数，得到空域矩阵Θ：

Θ＝[α_l,α_l+Δα,α_l+2Δα,…,α_r],

其中，表示Θ的左边界，表示Θ的右边界，表示半功率波束宽度，Δα为步长，Δα=0.1°；

5）在空域O内计算波束变换矩阵B；

6）对波束变换矩阵B进行预白化处理，得到白化波束变换矩阵T；

7）将目标信号x投影到白化波束变换矩阵T，得到波束域量测信号z＝Tx；

8）用F×P维观测矩阵φ对量测信号z进行压缩采样，F＜＜P，得到F×1维的观测信号y＝Φz；

9）根据观测信号y和白化波束变换矩阵T，利用贪婪追踪类算法中的正交匹配追踪方法，通过式迭代，选取一个局部最优解逐步逼近空域稀疏信号S，得到空域稀疏信号S的估计值：

$\hat{S}=[{\hat{s}}_1,{\hat{s}}_2,\·\·\·,{\hat{s}}_i,\·\·\·,{\hat{s}}_N],$

其中，|| ||₁表示求向量1-范数，s.t表示约束条件，|| ||₂表示求向量2-范数，β为噪声标准差；

10）定义目标角度范围θ＝[θ₁，θ₂，…，θ_i，…，θ_N]，根据所得估计值的元素与目标角度范围θ的元素的一一对应关系，即与θ_i一一对应，得到目标测角结果θ_d：

其中，d表示估计值中不为零的元素的下标，d∈[1,N]；

11）根据目标测角结果θ_d和已知的目标距离R，通过三角变换得到目标高度：

H＝Rsin(θ_d)。

本发明与现有技术相比具有如下优点：

（1）本发明利用目标空域稀疏的特点，将压缩感知方法引入到波达方向DOA估计问题中，突破了阵列分辨率的瑞利限，降低了信号功率谱和空间谱的旁瓣，提高了测角精度。

（2）本发明由于采用观测矩阵对量测信号进行压缩采样处理，减少了运算量，有效提高了米波雷达测高方法的性能。

（3）相比于应用最为广泛的子空间类算法，本发明在样本数较少时所得的目标信号估计结果更优。

仿真结果表明，本发明具有更高的角度分辨力。

附图说明

图1是本发明的实现流程图；

图2是用本发明与现有两种测高方法在信噪比变化时的波达方向估计性能比较图；

图3是用本发明与现有两种测高方法对目标角度估计的结果比较图。

具体实施方式

下面结合附图详细说明本发明的内容和效果。

参照图1，本发明的实现步骤如下：

步骤1：从雷达回波中提取目标信号。

所用阵列雷达为一垂直放置的均匀线阵，该均匀线阵由M个阵元组成，阵元间隔为d。

假设有K个远场窄带信号入射到该均匀线阵，M>K，信号入射角为α_i，i=1,2,…,K，则阵列接收到的目标信号为：

x＝AS₀+n，

其中，x为M×1维的阵元接收数据，n为M×1维的白噪声，满足零均值、方差为σ²的复高斯分布，各阵元输出噪声统计独立；S₀＝[s₁,s₂,…,s_K]^T为K×1维的信号矢量；A为M×K维阵列流型矩阵：

A＝[v(α₁),v(α₂),…,v(α_i),…,v(α_K)]，

其中，为第i个目标信号的导向矢量，上标T表示转置，λ为入射信号的波长。

步骤2：对目标信号x采用空间网格划分处理。

为了表现目标信号的空域稀疏性，对目标信号x需要采用空间网格划分处理，即将空间-180°～180°划分为ξ＝{α₁,α₂,…,α_u,…,α_N}，ξ表示目标搜索角度范围，α_u为第u个角度区间，u＝1,2,…,N，N＞＞K。该空间网格划分步骤如下：

2a）假设每个α_u都与一个目标信号s_u相对应，这样就构造一个N×1维的空域稀疏信号：S＝[s₁,s₂,…,s_u,…,s_N]^T；

2b）将目标信号x投影到S，则在S中只有实际存在目标信号的K个位置的元素不为零，其他N-K个位置的元素均为零，得到空域稀疏信号S：

S＝(x-n)ψ^-1，

其中，上标T表示转置，ψ为信号稀疏化后对应的M×N维阵列流型，也称超完备冗余字典；x和S所包含的目标信息是一致的，不同的是，x是目标信号在阵元域的表示，S是目标信号在空域的表示；

2c）根据空域稀疏信号S，将步骤1所述的目标信号x写为：

x＝ψS+n。

步骤3：使用数字波束形成法DBF对改写后的目标信号x进行仰角粗测，得到粗侧角度α，进而得到目标信号仰角所在空域O。

3a）利用导向矢量v(ξ)=[,e^-j2πsin(ξ),…,e^{-j2π(M-1)sin(ξ)}]^T对改写后的目标信号x进行加权求和，得到粗测角度α：

$\mathrm{\α}=\arg \underset{\mathrm{\ξ}}{\max }\left(\frac{1}{L}\underset{l=1}{\overset{L}{\mathrm{\Σ}}}\right|v^H\left(\mathrm{\ξ}\right)x\left(t_l\right){|}^2),$

其中，arg max表示寻找具有最大代价函数的参量，ξ表示目标搜索角度范围，-180°≤ξ≤180°,L表示快拍数，M表示阵元个数，x(t_l)表示t_l时刻采样到的目标信号，1≤l≤L，上标T表示转置，上标H表示共轭转置；

3b）计算半功率波束宽度：进而得到目标角度所在的空域Ο：

其中，λ表示雷达信号波长，d表示阵元间距。

步骤4：将目标信号仰角所在空域O划分为P份，P＞＞M，M表示阵元数，得到空域矩阵Θ：

Θ＝[α_l,α_l+Δα,α_l+2Δα,…,α_r],

其中，表示Θ的左边界，表示Θ的右边界，表示半功率波束宽度，Δα为步长，Δα=0.1°。

步骤5：在所述空域O内根据导向矢量计算波束变换矩阵B，按以下公式计算：

$B=[v\left({\mathrm{\α}}_l\right),v({\mathrm{\α}}_l+\mathrm{\Δ\α}),...,v\left({\mathrm{\α}}_j\right),...,v\left({\mathrm{\α}}_r\right){]}^H/\sqrt{M}$

其中，表示入射角为α_j时，目标信号的导向矢量，α_j∈Θ，α_l≤α_j≤α_r，Θ＝[α_l,α_l+Δα,α_l+2Δα,…,α_r]，Δα为步长，Δα=0.1°，M表示阵元个数，上标T表示转置，上标H表示共轭转置。

步骤6：对波束变换矩阵B进行预白化处理，得到白化波束变换矩阵T。

6a）对波束变换矩阵B的自相关矩阵R_B进行特征值分解，即：

R_B＝B(B^HB)^-1B^H＝QΣQ^H,

其中，Q表示正交矩阵，Q＝B，Σ表示对角矩阵,Σ＝（B^HB）^-1,上标H表示共轭转置；

6b）根据正交矩阵Q和对角矩阵Σ，通过预白化公式得到白化波束变换矩阵T：

T＝Σ^1/2Q^H＝(B^HB)^-1/2B^H。

步骤7：将目标信号x投影到白化波束变换矩阵T，得到波束域量测信号z：

z＝Tx＝TψS+Tn

其中，ψ表示超完备冗余字典，n表示白噪声，S表示空域稀疏信号。

步骤8：用F×P维观测矩阵φ对量测信号z进行压缩采样，F＜＜P，P表示空域划分的份数，得到F×1维的观测信号y：

y＝Φz＝ΦTψS+ΦTn。

步骤9：根据观测信号y和白化波束变换矩阵T，利用贪婪追踪类算法中的正交匹配追踪方法，通过式迭代，选取一个局部最优解逐步逼近空域稀疏信号S，得到空域稀疏信号S的估计值：

$\hat{S}=[{\hat{s}}_1,{\hat{s}}_2,\·\·\·,{\hat{s}}_i,\·\·\·,{\hat{s}}_N],$

其中，|| ||₁表示求向量1-范数，s.t表示约束条件，|| ||₂表示求向量2-范数，β为噪声标准差。

步骤10：定义目标角度范围θ＝[θ₁，θ₂，…，θ_i，…，θ_N]，根据所得估计值的元素与目标角度范围θ的元素的一一对应关系，即与θ_i一一对应，得到目标测角结果θ_d：

其中，d表示估计值中不为零的元素的下标，d∈[1,N]。

步骤11：根据目标测角结果θ_d和已知的目标距离R，通过三角变换得到目标高度：

H＝Rsin(θ_d)。

本发明的优点和效果通过以下计算仿真进一步说明：

1.仿真条件

仿真条件为以下雷达参数：垂直布置的20个水平极化天线阵元组成的等距布阵，雷达架高20m，地面反射系数为-0.95，载波频率为300MHz，只考虑镜面反射。本发明算法中的波束数为P＝200，观测矩阵维数为20。

2.仿真内容

仿真一：选取单个目标，在目标与参考天线的距离为200km，目标直达角为2°，多径反射角为-2.01°，阵元信噪比从-10dB变化至30dB，快拍数为10的条件下，用前后向空间平滑MUSIC算法、交替投影最大似然算法和本发明分别对低仰角目标进行测角精度仿真，仿真结果如图2所示。其中横轴表示信噪比从-10分贝至20分贝变化，纵轴表示测角误差；图2中的SS-MUSIC曲线表示前后向空间平滑MUSIC算法在信噪比按照横轴变化时的测角误差；APML曲线表示交替投影最大似然算法在信噪比按照横轴变化时的测角误差；BS-CS曲线表示本发明在信噪比按照横轴变化时的测角误差。

从图2中可以得出，对低仰角目标现有的前后向空间平滑MUSIC算法、交替投影最大似然算法测角误差偏大，而本发明的测角误差最小。

仿真二：选取单个目标，在目标高度为12000m，径向从50km飞至650km，阵元间距为半波长，信噪比为10dB，快拍数为10，蒙特卡罗实验次数100次的条件下，用前后向空间平滑MUSIC算法、交替投影最大似然算法和本发明分别针对不同仰角对算法估计精度的影响进行仿真。仿真结果如图3所示。其中：

图3（a）表示前后向空间平滑MUSIC算法在目标与阵地的距离按照横轴变化时的仰角；

图3（b）表示交替投影最大似然算法在目标与阵地的距离按照横轴变化时的仰角；

图3（c）表示本发明在目标与阵地的距离按照横轴变化时的仰角。

图3中横轴表示目标与阵地的距离从0千米至650千米变化，纵轴表示仰角。从

图3中可以得出，相对于前后向空间平滑MUSIC算法、交替投影最大似然算法，

本发明角度估计误差更小。

Claims (6)

1.一种基于压缩感知的波束空间域米波雷达测高方法，包括以下步骤：

1)从雷达回波中提取目标信号：

x＝AS₀+n

其中，A为M×K维阵列流型，S₀＝[s₁,s₂,…,s_K]^T为K×1维的信号矢量，n为M×1维的白噪声，满足零均值、方差为σ²的复高斯分布，各阵元输出噪声统计独立；

2)根据信源是空域稀疏的特性，采用空间网格划分方法，将空间划分为N份，构造一个信号稀疏化后对应的M×N维阵列流型ψ和N×1维的空域稀疏信号S＝[s₁,s₂,…,s_N]^T，N＞＞K，那么目标信号x改写为：

x＝ψS+n；

3)使用数字波束形成法DBF对目标信号x进行仰角粗测，得到粗测角度α，进而得到目标信号仰角所在空域O；

4)将所述空域O划分为P份，P＞＞M，M表示阵元数，得到空域矩阵Θ：

Θ＝[α_l,α_l+Δα,α_l+2Δα,…,α_r],

其中，表示Θ的左边界，表示Θ的右边界，表示半功率波束宽度，Δα为步长，Δα＝0.1°；

5)在空域O内计算波束变换矩阵B；

6)对波束变换矩阵B进行预白化处理，得到白化波束变换矩阵T；

7)将目标信号x投影到白化波束变换矩阵T，得到波束域量测信号z＝Tx；

8)用F×P维观测矩阵φ对量测信号z进行压缩采样，F＜＜P，得到F×1维的观测信号y＝Φz；

9)根据观测信号y和白化波束变换矩阵T，利用贪婪追踪类算法中的正交匹配追踪方法，通过式迭代，选取一个局部最优解逐步逼近空域稀疏信号S，得到空域稀疏信号S的估计值

$\hat{S}=[{\hat{s}}_1,{\hat{s}}_2,...,{\hat{s}}_i,...,{\hat{s}}_N],$

其中，|| ||₁表示求向量1-范数，s.t表示约束条件，|| ||₂表示求向量2-范数，β为噪声标准差；

10)定义目标角度范围θ＝[θ₁,θ₂,…,θ_i,…,θ_N]，根据所得估计值的元素与目标角度范围θ的元素的一一对应关系，即与θ_i一一对应，得到目标测角结果θ_d：

其中，d表示估计值中不为零的元素的下标，d∈[1,N]；

11)根据目标测角结果θ_d和已知的目标距离R，通过三角变换得到目标高度：

H＝Rsin(θ_d)。

2.根据权利要求1所述的基于压缩感知的波束空间域米波雷达测高方法，其中步骤3)所述的使用数字波束形成法DBF对目标信号x进行仰角粗测，进而得到目标信号仰角所在空域Ο，按如下步骤进行：

2a)利用导向矢量v(ξ)＝[1,e^-j2πsin(ξ),…,e^{-j2π(M-1)sin(ξ)}]^T，对目标信号x进行加权求和，得到粗测角度α：

$\mathrm{\α}=\arg \underset{\mathrm{\ξ}}{\max }\left(\frac{1}{L}\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{\left|v^H\left(\mathrm{\ξ}\right)x\left(t_l\right)\right|}^2\right),$

其中，arg max表示寻找具有最大代价函数的参量，ξ表示目标搜索角度范围，-180°≤ξ≤180°,L表示快拍数，x(t_l)表示t_l时刻采样到的目标信号，1≤l≤L，上标T表示转置，上标H表示共轭转置；

2b)计算半功率波束宽度：进而得到目标角度所在的空域Ο：

其中，λ表示雷达信号波长，d表示阵元间距。

3.根据权利要求1所述的基于压缩感知的波束空间域米波雷达测高方法，其中步骤5)所述的在空域O内计算波束变换矩阵B，按如下公式进行：

$B={[v\left({\mathrm{\α}}_l\right),v({\mathrm{\α}}_l+\mathrm{\Δ\α}),...,v\left({\mathrm{\α}}_j\right),...,v\left({\mathrm{\α}}_r\right)]}^H/\sqrt{M},$

其中，表示入射角为α_j时，目标信号的导向矢量，α_j∈Θ，α_l≤α_j≤α_r，上标T表示转置，上标H表示共轭转置。

4.根据权利要求1所述的基于压缩感知的波束空间域米波雷达测高方法，其中步骤6)所述的对波束变换矩阵进行预白化处理，按如下公式进行：

T＝(B^HB)^-1/2B^H

其中，T表示白化波束变换矩阵，B表示波束变换矩阵，上标H表示共轭转置。

5.根据权利要求1所述的基于压缩感知的波束空间域米波雷达测高方法，其中步骤7)所述的将接收数据投影到白化波束变换矩阵，按如下公式进行：

z＝Tx＝TψS+Tn，

其中，z表示波束域量测信号，T表示白化波束变换矩阵，x表示接收到的目标信号。

6.根据权利要求1所述的基于压缩感知的波束空间域米波雷达测高方法，其中步骤8)所述的用F×P维观测矩阵对量测信号进行压缩采样，按如下公式进行：

y＝Φz＝ΦTψS+ΦTn

其中，y表示观测信号，Φ表示观测矩阵，z表示波束域量测信号。