基于信杂噪比最大的雷达三维异构阵稀疏重构方法与流程

本发明属于雷达技术领域，涉及一种基于信杂噪比最大的雷达三维异构阵稀疏重构方法，适用于机载三维异构阵雷达的杂波抑制。

背景技术：

与机身共形的天线具有诸多优点，比如更好的气动特性、更好的隐身性能、更大的天线面积等，因此其应用前景广阔；但机身的多曲面特性使得天线也具有多曲面结构，进而产生三维异构阵问题；相对于二维平面阵来说，三维异构阵的杂波特性将会变得十分复杂，这种复杂性体现在单元间的互耦特性更加复杂，单元的极化改变更加复杂；三维异构阵给杂波带来的复杂性也给杂波抑制带来了很大的困难包括阵列流型更加复杂，目标和杂波的导向矢量计算困难，导向矢量间的相关性更难计算，杂波呈现距离相关性非均匀。

在杂波协方差矩阵和目标信号均确知的条件下，brennan等人在1973年，根据线性约束最小方差准则，推导出输出信杂噪比最大的全空时二维自适应处理(spacetimeadaptiveprocessing，stap)的概念和理论，全stap能够取得比较理想的杂波抑制效果，但是最优的stap处理器需要用大量的训练样本来估计杂波加噪声协方差矩阵。

在三维异构阵的情况下，由于雷达接收的空时数据维数往往很大，全维stap处理器在实际应用中面临两个问题：1)三维阵的阵面变大，阵元数变多使空时两维自由度庞大导致杂波协方差矩阵求逆的运算量十分巨大，无法满足实时处理的要求；2)缺乏用于估计杂波协方差矩阵的独立同分布(independentidenticallydistributed,iid)样本单元。根据reed,mallett,andbrennan三个人提出的准则，即rmb准则可知，为使因估计不准确而带来的输出信杂噪比损失在3db以内，要求样本数目至少为协方差矩阵阶数的两倍，但实际工程中会面临强烈的非均匀杂波环境，不同距离的杂波样本分布特性各异，这将造成独立同分布iid样本支撑困难。因此，这两个问题使得在三维阵列的情况下使用全维stap性能降低，杂波抑制效果变差。

为了摆脱在三维阵情况下雷达接接收数据很大时杂波协方差矩阵得估计对于样本的过度依赖，导致的杂波抑制效果降低的问题，有必要对三维阵列的空时自适应处理降维处理进行研究。

技术实现要素：

针对上述现有技术存在的问题，本发明的目的在于提出一种基于信杂噪比最大的雷达三维异构阵稀疏重构方法，该种基于信杂噪比最大的雷达三维异构阵稀疏重构方法能够解决三维异构阵情况下全空时自适应处理的杂波协方差矩阵估计需要大量训练样本的问题，能够降低杂波协方差矩阵估计所需的训练样本数，提高杂波抑制性能。

本发明的主要思路：利用三维异构阵杂波回波数据得到杂波的协方差矩阵和杂波加目标的协方差矩阵，推导出输出信杂噪比的行列式表示形式，用使输出信杂噪比最大的阵元选择方法得到最优的阵元选择矢量进而得到降维矩阵，得到降维后的改善因子和杂波抑制结果。

为达到上述技术目的，本发明采用如下技术方案予以实现。

一种基于信杂噪比最大的雷达三维异构阵稀疏重构方法，包括以下步骤：

步骤1，确定三维圆柱阵雷达，所述三维圆柱阵雷达的检测范围内存在目标，并获取三维异构阵机载雷达杂波空时导向矢量矩阵an×k×l；

其中，n为三维圆柱阵雷达包含的阵元个数，k为每个相干处理间隔内三维圆柱阵雷达发射的脉冲个数，l表示表示三维圆柱阵雷达接收到的杂波包含的距离门总个数，n、k、l分别为大于0的正整数；

步骤2，将三维异构阵机载雷达杂波空时导向矢量矩阵an×k×l转换为nk×l的二维导向矢量矩阵bnk×l，并计算得到目标加杂波的协方差矩阵rs的行列式；

步骤3，计算三维异构阵机载雷达的输出信杂噪比scnrout；

步骤4，计算得到三维异构阵机载雷达的输出信杂噪比scnrout的行列式表示形式；

步骤5，定义二进制挑选向量y，所述二进制挑选向量y为nk×1维列矢量，并根据三维异构阵机载雷达的输出信杂噪比scnrout的行列式表示形式，计算得到三维异构阵机载雷达的输出信杂噪比scnrout的行列式表示形式的重写式；

步骤6，根据三维异构阵机载雷达的输出信杂噪比scnrout的行列式表示形式的重写式,得到三维异构阵机载雷达的稀疏降维矩阵；

步骤7，根据三维异构阵机载雷达的稀疏降维矩阵,计算得到三维异构阵机载雷达的滤波结果；所述三维异构阵机载雷达的滤波结果为基于信杂噪比最大的雷达三维异构阵稀疏重构结果。

本发明的有益效果：在三维异构阵情况下，由于阵元数增多使得全空时自适应处理的自由度明显增大，从而需要更多的训练样本来估计杂波协方差矩阵，本发明方法通过在信杂噪比最大的情况下进行阵元挑选，挑选后自由度降低为全空时的一半，能够有效抑制杂波，并能够很好解决三维异构阵情况下杂波协方差矩阵估计的训练样本数有限的问题。

附图说明

下面结合附图和具体实施方式对本发明作进一步详细说明。

图1是本发明的一种基于信杂噪比最大的雷达三维异构阵稀疏重构方法流程图；

图2(a)为在少样本数情况下使用本发明方法与全空时处理方法分别得到的改善因子曲线对比图；

图2(b)为在多样本数情况下使用本发明方法与全空时处理方法分别得到的改善因子曲线对比图；；

图3(a)为少样本数情况下使用全空时处理方法进行杂波抑制得到的结果示意图；

图3(b)为少样本数情况下使用本发明方法进行杂波抑制得到的结果示意图。

具体实施方式

参照图1，为本发明的一种基于信杂噪比最大的雷达三维异构阵稀疏重构方法流程图；其中所述基于信杂噪比最大的雷达三维异构阵稀疏重构方法，包括以下步骤：

步骤1，确定三维圆柱阵雷达，所述三维圆柱阵雷达的检测范围内存在目标，并获取三维异构阵机载雷达杂波空时导向矢量矩阵an×k×l；本实施例是对三维圆柱阵雷达进行仿真得到三维异构阵机载雷达杂波空时导向矢量矩阵。

具体地：an×k×l为n×k×l维三维数据矩阵，n为三维圆柱阵雷达包含的阵元个数，三维圆柱阵雷达共五层，每层六个阵元，k为每个相干处理间隔cpi内三维圆柱阵雷达发射的脉冲个数，cpi表示相干处理间隔，l表示表示三维圆柱阵雷达接收到的杂波包含的距离门总个数，n、k、l分别为大于0的正整数。

步骤2，将三维异构阵机载雷达杂波空时导向矢量矩阵an×k×l转换为nk×l的二维导向矢量矩阵bnk×l，并计算得到目标加杂波的协方差矩阵rs的行列式。

具体地，将三维异构阵机载雷达杂波空时导向矢量矩阵an×k×l转换为nk×l的二维导向矢量矩阵bnk×l，计算nk×l的二维导向矢量矩阵bnk×l的秩na，na＝rank(bnk×l)，0≤na≤min{nk,l}，rank表示求秩运算，min表示求最小值操作。

然后对nk×l的二维导向矢量矩阵bnk×l进行奇异值分解，得到nk×nk维酉矩阵unk×nk，取nk×nk维酉矩阵unk×nk的前na列，记为nk×na维中间矩阵vc，

vi表示nk×nk维酉矩阵unk×nk中的第i列矢量，且vi为nk×1维列矢量。

根据nk×na维中间矩阵vc，计算得到三维异构阵机载雷达的杂波协方差矩阵rc，其定义式为：

分别设定目标的方位角为φ，目标的俯仰角为θ，目标的多普勒频率为fd，则将目标的空间频率记为fs，将目标的归一化多普勒频率记为n为三维圆柱阵雷达包含的阵元个数,λ表示三维圆柱阵雷达发射信号波长，fr为脉冲重复频率，cos表示求余弦操作，sin表示求正弦操作，那么目标的空域导向矢量为a，a为n×1维列矢量，n∈{0,1,…,n-1}，an表示目标的空域导向矢量a中第n个元素，n为三维圆柱阵雷达包含的阵元个数，与目标的空域导向矢量a中包含的元素个数取值相等；上标t表示转置操作，目标的时域导向矢量为b，b为k×1维列矢量，bm表示目标的时域导向矢量b中第m个元素，k为每个相干处理间隔cpi内三维圆柱阵雷达发射的脉冲个数，与目标的时域导向矢量b中包含的元素个数取值相等；cpi表示相干处理间隔；进而计算得到目标的空时导向矢量s，其计算表达式为：

其中，表示kronecker积，s为nk×1维列矢量，目标的空时导向矢量s的定义表达式为：s＝[s1,s2,...,si',...,snk]^t，上标t为转置操作，i'∈{1,2,…,nk}，si'表示目标的空时导向矢量s中第i'个元素，si'＝anbm。

定义目标加杂波的信号矩阵为vs，vs＝[s,vc]＝[s,v1,v2,...,vi,...,vna]，进而计算得到目标加杂波的协方差矩阵为rs，s^hs＝nk，上标h表示共轭转置操作；根据块矩阵求行列式的公式，得到目标加杂波的协方差矩阵rs的行列式|rs|，

步骤3，计算三维异构阵机载雷达的输出信杂噪比scnrout。

具体地：(1)将三维异构阵机载雷达的杂波加噪声协方差矩阵表示为r，且r为nk×nk的二维矩阵，表示高斯白噪声的功率，rc为三维异构阵机载雷达的杂波协方差矩阵，ink表示n×k维全1矩阵。

由步骤1和步骤2可知三维异构阵机载雷达的杂波协方差矩阵rc的计算表达式为：vi表示nk×nk维酉矩阵unk×nk中的第i列矢量，且vi为nk×1维列矢量；qi表示nk×nk维酉矩阵unk×nk中第i列矢量的杂波功率，上标h表示共轭转置操作。

定义na×1维矢量并对na×1维矢量进行对角化，即是将na×1维矢量对角化，进而得到na×na维二维矩阵q，diag表示对角化操作，得到三维异构阵机载雷达的杂波加噪声协方差矩阵的简化形式

根据矩阵求逆公式计算得到三维异构阵机载雷达的杂波加噪声协方差矩阵的简化形式的逆

假定nk×nk维酉矩阵unk×nk中每一列矢量的杂波功率都远大于噪声功率，即qi>>则得到三维异构阵机载雷达的杂波加噪声协方差矩阵的简化形式的逆的简化式进而计算得到三维异构阵机载雷达的空时自适应处理的最优权wopt，

(2)将目标的空时导向矢量s分解为sc和s⊥两个正交的子空间，s＝sc+s⊥，sc和s⊥正交，并且互为补空间，即sc∪s⊥＝s，表示空集，sc表示杂波子空间，s⊥表示杂波子空间sc的正交补空间，vc表示nk×na维中间矩阵，ink表示n×k维全1矩阵，上标h表示共轭转置，k为每个相干处理间隔cpi内三维圆柱阵雷达发射的脉冲个数，n为三维圆柱阵雷达包含的阵元个数，na表示nk×l的二维导向矢量矩阵bnk×l的秩，l表示三维圆柱阵雷达接收到的杂波包含的距离门总个数。

由步骤1可知，nk×na维中间矩阵vi表示nk×nk维酉矩阵unk×nk中的第i列矢量，且vi为nk×1维列矢量，na个nk×1维列矢量张成了杂波子空间sc，p表示正交投影矩阵，根据正交投影的概念，上标h表示共轭转置操作，上标-1表示求逆操作；由此计算得到目标的空时导向矢量s与杂波子空间sc之间夹角α的余弦值cos(α)，其中，|·|表示求绝对值操作，||s||2表示目标的空时导向矢量s的二范数，将sc的表达式代入|cos(α)|的表达式后得到|cos(α)|的展开式，再对|cos(α)|的展开式进行求平方操作，得到目标的空时导向矢量s与杂波子空间sc之间夹角α的余弦平方值|cos(α)|²，

(3)根据目标的空时导向矢量s与杂波子空间sc之间夹角α的余弦平方值|cos(α)|²，计算得到三维异构阵机载雷达的输出信杂噪比scnrout，

s表示目标的空时导向矢量，表示三维异构阵机载雷达的杂波加噪声协方差矩阵的简化形式的逆的简化式，snr表示目标的信噪比，表示目标的信号功率，代表高斯白噪声的功率，n为三维圆柱阵雷达包含的阵元个数，三维圆柱阵雷达共五层，每层六个阵元，k表示每个相干处理间隔cpi内三维圆柱阵雷达发射的脉冲个数。

步骤4：计算得到三维异构阵机载雷达的输出信杂噪比scnrout的行列式表示形式rc为三维异构阵机载雷达的杂波协方差矩阵，rs表示目标加杂波的协方差矩阵，|rs|表示目标加杂波的协方差矩阵rs的行列式，|rc|表示三维异构阵机载雷达的杂波协方差矩阵rc的行列式。

具体地，对步骤2得到的目标加杂波的协方差矩阵rs的行列式

进行变形，得到然后将代入目标的空时导向矢量s与杂波子空间sc之间夹角α的余弦平方值|cos(α)|²的表达式中，得到最后将代入到scnrout的表达式中，可将三维异构阵机载雷达的输出信杂噪比scnrout改写为用rs和rc的行列式表示的形式，进而得到三维异构阵机载雷达的输出信杂噪比scnrout的行列式表示形式

步骤5：定义一个二进制挑选向量y，所述二进制挑选向量y为nk×1维列矢量，n为三维圆柱阵雷达包含的阵元个数，k为每个相干处理间隔cpi内三维圆柱阵雷达发射的脉冲个数，nk表示一共有nk个阵元脉冲对；二进制挑选向量y中包含nk个元素，每个元素的取值分别为0或1，1代表选择了对应的阵元脉冲，0代表没有选择任何阵元脉冲，那么将三维异构阵机载雷达的杂波协方差矩阵rc和目标加杂波的协方差矩阵rs分别表示为rc(y)和rs(y)，diag表示对角化操作，rc(y)表示被挑选后的杂波协方差矩阵，rs(y)表示被挑选后的杂波加目标协方差矩阵；进而得到三维异构阵机载雷达的输出信杂噪比scnrout的行列式表示形式的重写式

因此，选择二进制挑选向量y的取值问题就转化为使三维异构阵机载雷达的输出信杂噪比scnrout的行列式表示形式的重写式最大的问题，也就是让最小的问题。

步骤6，在三维异构阵机载雷达的输出信杂噪比scnrout的行列式表示形式的重写式取值最大的情况下进行阵元挑选。

具体地：假设在总体的n个阵元与k个脉冲的乘积nk个阵元脉冲对中挑选m个阵元脉冲对，则二进制挑选向量y里面有m个1，nk-m个0；所以，在三维异构阵机载雷达的输出信杂噪比scnrout的行列式表示形式的重写式取值最大的情况下进行阵元挑选的目的是使得二进制挑选向量y里面有m个1，nk-m个0。

6.1令二进制挑选向量为y，并令将二进制挑选向量y的初始值全部为1，记为二进制挑选向量初始值y⁽⁰⁾，即y⁽⁰⁾＝[1,1,...,1]nk，二进制挑选向量初始值y⁽⁰⁾的下标集为η⁽⁰⁾，η⁽⁰⁾＝[1,2,...,nk]。

6.2定义外层循环变量k，k为第k次迭代，k的初始值为1，k∈{1,2,…,nk-m}，m<nk；定义内层循环变量g，g为第g次迭代，g随着k的变化而变化,g∈{1,2,…,nk-k+1},g的初始值为1；当k取值为1时，g取值从1循环到nk；当k取值为2时，g取值从1循环到nk-1；当k取值为m时，g取值从1循环到nk-m+1；外层循环变量每循环一次，就会将第k次迭代后的二进制挑选向量y^(k)中的一个1置为0，当外层循环变量循环m次时，就会将第m次迭代后的二进制挑选向量y^(m)中的m个1都置为0。

6.3根据第k次迭代后的二进制挑选向量y^(k),得到第k次迭代后的二进制挑选向量y^(k)的下标集η^(k),所述第k次迭代后的二进制挑选向量y^(k)的下标集η^(k)中与第k次迭代后的二进制挑选向量y^(k)中元素为0的对应位置处元素为0,y^(k)的下标集η^(k)中其余位置处的元素按照自然数顺序从1排列至nk-k+1；所述第k次迭代后的二进制挑选向量y^(k)为第k-1次迭代后的二进制挑选向量y^(k-1)。

6.4令第g次迭代后下标集η^(k)中与g取值相等的元素在η^(k)中的位置为h，即η^(k)(h)＝g，h∈{1,2,...,nk-k+1}，然后将第k次迭代后的二进制挑选向量y^(k)中位置h处的元素置为0，即y^(k)(h)＝0，并将第k次迭代后的二进制挑选向量y^(k)中位置h处的元素置为0后得到的向量,记为第g次迭代后在位置h处置零的向量且满足即向量中位置h处的元素为0。

然后计算第g次迭代后在位置h处置零的向量的杂波协方差矩阵和第g次迭代后在位置h处置零的向量的杂波加目标协方差矩阵以及分别计算杂波协方差矩阵的行列式和杂波加目标协方差矩阵的行列式进而计算杂波协方差矩阵的行列式和杂波加目标协方差矩阵的行列式的比值μ^k(g)，其表达式分别为:

其中，杂波协方差矩阵的行列式和杂波加目标协方差矩阵的行列式的维数都是nk×nk,diag表示对角化操作；然后,将杂波协方差矩阵的行列式和杂波加目标协方差矩阵的行列式的比值μ^k(g),记为第k次迭代后经过第g次迭代得到的第g个比值。

6.5令g分别取1至nk-k+1，重复执行6.3和6.4，进而分别得到第k次迭代后经过第1次迭代得到的第1个比值至第k次迭代后经过第nk-k+1次迭代得到的第nk-k+1个比值，记为第k次迭代后得到的nk-k+1个比值,比较第k次迭代后得到的nk-k+1个比值，得到比值最小值，并得到比值最小值时g的对应取值在第k次迭代后的二进制挑选向量y^(k)的下标集η^(k)中的位置，记为hmin，hmin∈{1,2,…,nk-k+1}，然后将第k次迭代后的二进制挑选向量y^(k)中位置hmin处的1置为0，即y^(k)(hmin)＝0，进而得到第k次迭代后的二进制挑选向量y^(k)，所述第k次迭代后的二进制挑选向量y^(k)中有k个0，nk-k个1。

6.6令k加1，返回6.3，直到得到第nk-m次迭代后的二进制挑选向量y^(nk-m)，迭代停止，此时得到的第nk-m次迭代后的二进制挑选向量y^(nk-m)中nk-m个位置处的值分别已置为0，其余m个位置处的值分别为1，即迭代的最终结果使得第nk-m次迭代后的二进制挑选向量y^(nk-m)中有m个1，nk-m个0。

然后计算得到对角矩阵y，y＝diag(y^(nk-m))，diag表示对角化操作，对角矩阵y为nk行nk列的二维矩阵，由于第nk-m次迭代后的二进制挑选向量y^(nk-m)中有nk-m个元素为0，因此对角矩阵y中有nk-m列全部为零，剔除对角矩阵y中全部为0的列，并将对角矩阵y中剔除全部为0的列后的矩阵，记为三维异构阵机载雷达的稀疏降维矩阵z，z是nk行m列的二维矩阵。

步骤7：由步骤3得到的空时自适应处理的最优权为使用三维异构阵机载雷达的稀疏降维矩阵z对三维异构阵机载雷达的杂波加噪声协方差矩阵r和目标的空时导向矢量s分别进行降维处理，分别计算得到稀疏后的自适应权值和稀疏后的改善因子if，进而计算得到三维异构阵机载雷达的滤波结果out，其表达式分别为：

其中，e表示三维异构阵雷达的杂波功率,上标h表示共轭转置操作，上标-1表示求逆操作；所述三维异构阵机载雷达的滤波结果out即为基于信杂噪比最大的雷达三维异构阵稀疏重构结果。

下面结合仿真实验对本发明效果作进一步验证说明：

(一)仿真参数：

三维异构阵机载雷达天线采用圆柱阵进行仿真，圆柱阵为5层，每层6个阵元，共30个阵元，波长为0.2m，阵元间距为半波长0.1m，相干积累脉冲数20，x轴向阵元间隔为0.1m，三维异构阵机载雷达的载机高度为6000m，三维异构阵机载雷达的载机速度为200m/s，载机速度方向沿x轴，三维异构阵机载雷达发射的信号带宽为2.5mhz，采样频率2.5mhz，脉冲重复频率为5000hz，信号带宽为2.5mhz，主波束指向为，方位角90°，俯仰角0°，噪声系数3db，杂噪比40db。

(二)仿真结果及分析：

为了说明本发明性能的优越性，给出了传统方法以及本发明方法的处理结果对比图。

参照图2(a)，为在少样本数情况下使用本发明方法与全空时处理方法分别得到的改善因子曲线对比图；其中，纵坐标表示改善因子的大小，单位为分贝(db)，横坐标表示多普勒通道，共20个，实线为全空时处理方法，虚线为本发明方法。

从图2(a)可以看出，由于采用了阵元挑选，选择了一半的阵元脉冲数，也就是nk/2＝300个阵元脉冲数，所以本发明方法的自由度较全空时处理的自由度nk＝600降低了一半，因此所需的训练样本数也是全空时的一半，因此在样本数为700时，对于本发明方法已经满足了两倍自由度的要求，而对于全空时处理方法并没有满足，因此在少样本情况下，本发明方法的改善因子要优于全空时处理的改善因子。

图2(b)为在多样本数情况下使用本发明方法与全空时处理方法分别得到的改善因子曲线对比图，选择的样本数为全空时情况下的五倍自由度，在多样本情况下，由于本发明方法直接去掉了一半的阵元脉冲，因此结果没有全空时的改善因子好，这也说明了，本发明方法更加适用于样本数缺乏的情况下进行杂波抑制。

图3(a)为少样本数情况下使用全空时处理方法进行杂波抑制得到的结果示意图，图3(b)为少样本数情况下使用本发明方法进行杂波抑制得到的结果示意图，图3(a)和图3(b)的横坐标均为多普勒通道数，共20个；纵坐标均为距离门数，从0到3000，共3000个距离门。

对比图3(a)和图3(b)可以看出在少样本情况下全空时处理方法的杂波剩余比较多，本发明方法在少样本情况下的滤波结果比全空时处理结果好，杂波剩余更少，因此杂波抑制性能更好。

综上所述，仿真实验验证了本发明的正确性，有效性和可靠性。

显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围；这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。