本发明公开了一种基于协方差矩阵广对称特性的STAP方法。其步骤为：（1）接收空时二维回波数据；（2）多普勒滤波处理；（3）获取广对称训练样本；（4）估计广对称协方差矩阵；（5）计算自适应权矢量；（6）自适应滤波处理；（7）输出结果。本发明主要解决机载雷达在非均匀杂波环境中，由于满足独立同分布条件的训练样本不足而造成协方差矩阵估计不准确的问题。本发明基于协方差矩阵的广对称特性，能够有效提高回波数据的利用率，在较少独立同分布训练样本的情况下更准确地估计杂波协方差矩阵，取得更好地杂波抑制和目标检测效果，在现实中极度非均匀的杂波环境中本发明具有更好的应用前景。
1.一种基于协方差矩阵广对称特性的STAP方法，包括以下步骤：(1)接收空时二维回波数据：
利用机载雷达天线的多个阵元，在相干积累时间内接收地面反射的空时二维回波数据；
按照下式，利用扩展因子化算法EFA中的降维转换矩阵，将机载雷达接收的空时二维回波数据，从阵元-脉冲域转化到阵元-多普勒域，得到一个待处理多普勒通道的原始回波数据：其中，x表示待处理多普勒通道的原始回波数据，TEFA表示扩展因子化算法EFA中的降维转换矩阵，表示空时二维回波数据；
3a)采用空域转换矩阵，对待处理多普勒通道的原始回波数据进行空域转换，得到空域转换数据；
3b)采用时域转换矩阵，对待处理多普勒通道的原始回波数据进行时域转换，得到时域转换数据；
3c)采用空时域转换矩阵，对待处理多普勒通道的原始回波数据进行空时域转换，得到空时域转换数据；
4a)利用最大似然方法，分别估计待处理多普勒通道的原始回波数据、空域转换数据、时域转换数据以及空时域转换数据对应的原始协方差矩阵、空域协方差矩阵、时域协方差矩阵以及空时域协方差矩阵；
4b)对上述四个协方差矩阵取平均，将取平均后的结果定义为广对称协方差矩阵；
按照自适应权矢量求解公式，用广对称协方差矩阵代替传统扩展因子化算法EFA中的原始协方差矩阵，获得自适应权矢量；
利用自适应权矢量，依次对待处理多普勒通道的原始回波数据按照距离单元分别进行空时自适应滤波；
2.根据权利要求1所述的基于协方差矩阵广对称特性的STAP方法，其特征在于，步骤
其中，xs表示空域转换数据，s表示空域符号， 表示空域转换矩阵，I3表示空域转换矩阵中的3维单位矩阵，表示Kronecker积运算符号，J表示空域转换矩阵中的N维置换矩阵，N表示机载雷达天线的阵元个数，在置换矩阵中反对角线元素为1，其余元素均为零，x表示待处理多普勒通道的原始回波数据，*表示共轭运算符号。
3.根据权利要求1所述的基于协方差矩阵广对称特性的STAP方法，其特征在于，步骤
其中，xt表示时域转换数据，t表示时域符号， 表示时域转换矩阵，U表示时域转换矩阵中的3维对角矩阵，IN表示时域转换矩阵中的N维单位矩阵，N表示机载雷达天线的阵元个数，x表示待处理多普勒通道的原始回波数据。
4.根据权利要求1所述的基于协方差矩阵广对称特性的STAP方法，其特征在于，步骤
其中，xst表示空时域转换数据，st表示空时域符号， 表示空时域转换矩阵，V表示空时域转换矩阵中的3维对角矩阵，J表示空时域转换矩阵中的N维置换矩阵，N表示机载雷达天线的阵元个数，x表示待处理多普勒通道的原始回波数据。
5.根据权利要求1所述的基于协方差矩阵广对称特性的STAP方法，其特征在于，步骤(4)中所述原始协方差矩阵、空域协方差矩阵、时域协方差矩阵以及空时域协方差矩阵采用相同的最大似然方法进行估计的公式如下：其中，R表示原始协方差矩阵、空域协方差矩阵、时域协方差矩阵以及空时域协方差矩阵的其中一个，xi表示待处理多普勒通道的原始回波数据、空域转换数据、时域转换数据以及空时域转换数据的其中一个，i表示数据xi的索引号，L表示数据xi的长度，H表示共轭转置运算符号。
6.根据权利要求1所述的基于协方差矩阵广对称特性的STAP方法，其特征在于，步骤(5)中所述的自适应权矢量求解公式如下：其中，W表示自适应权矢量，μ表示由雷达系统参数设定的一个常数，RPer表示广对称协方差矩阵，per表示广对称，S表示降维后的目标导向矢量，其表达式为TEFA表示扩展因子化算法EFA中的降维转换矩阵，st表示目标时域导向矢量，t表示时域符号，ss表示目标空域导向矢量，s表示空域符号， 表示Kronecker积运算符号。
基于协方差矩阵广对称特性的STAP方法
技术领域
[0001] 本发明属于通信技术领域，更进一步涉及雷达技术领域中的一种基于协方差矩阵广对称特性的机载雷达空时自适应处理（space-time adaptive processing，STAP）方法。本发明主要用于解决机载雷达在非均匀杂波环境中，由于满足独立同分布条件的训练样本不足而造成的协方差矩阵估计不准确问题，能够有效增加训练样本，改善自适应信号处理的杂波抑制性能，提高目标的检测概率。
背景技术
[0002] 机载预警雷达的主要任务是在复杂杂波背景中探测目标并对其进行定位跟踪，而对杂波进行有效抑制是提高预警雷达工作性能的核心手段。空时自适应处理STAP技术充分利用空域和时域信息，在对目标信号进行相干积累的同时，通过空时自适应处理滤除地面杂波，实现机载雷达对目标的有效检测。如美国的E2-D机载预警雷达就采用此技术。在实际应用中，STAP处理器主要存在以下两方面问题：首先，在非均匀的杂波环境中，要获得足够多的用于估计协方差矩阵的独立同分布（independent and identically distributed，IID）训练样本非常困难；其次，即使训练样本的需求得到满足，全空时处理计算量过大的问题会导致实时性难以保证。为解决上述问题，推动STAP技术更加实用化，人们提出了许多改进措施或方法。
[0003] 清华大学申请的发明专利“非均匀杂波环境下空时自适应处理方法”（专利申请号201010129723.3，公布号CN101819269A）公开了一种在非均匀杂波环境中超分辨估计杂波空时二维谱的超完备稀疏表示方法。该方法实现了在独立同分布样本数不足的情况下，利用单帧训练样本估计杂波协方差矩阵，从而避免强非均匀杂波环境对自适应处理效果的影响。但是，该方法仍然存在的主要不足是：对杂波谱进行稀疏表示的超完备基数目未定，但远大于系统自由度，而实际中系统自由度通常成千上万，这样在每一距离单元样本的协方差矩阵重构过程中所需要的运算量非常大，不利于实时处理，从而影响到它的实际工程应用效果。
[0004] 北京理工大学申请的发明专利“一种基于协方差矩阵加权的降维空时自适应处理方法”（专利申请号201210251589.3，公布号CN102778669A）公开了一种利用协方差矩阵加权技术估计杂波协方差矩阵的方法。该方法实现了自适应地展宽杂波凹口来适应实际环境中的杂波脊，从而使得存在杂波泄露时，也能通过STAP方法有效地抑制杂波，提高了STAP在实际应用中的稳健性。但是，该方法仍然存在的不足是：杂波凹口的展宽程度是人为设定的，不能自适应地感知实际数据中的杂波脊情况，当展宽量过大时，会造成最小可检测速度变大，对低速小目标的检测非常不利。
发明内容
[0005] 本发明针对上述现有STAP技术的不足，提出了一种基于协方差矩阵广对称特性的机载雷达空时自适应处理（space-time adaptive processing，STAP）方法。本发明通过提高数据利用效率有效增加独立同分布训练样本，从而获得实际数据中的真实杂波情况，提高杂波协方差矩阵估计的准确性，改善杂波抑制性能，提高动目标的检测概率。
[0006] 为实现上述目的，本发明实现的具体步骤如下：
[0007] （1）接收空时二维回波数据：
[0008] 利用机载雷达天线的多个阵元，在相干积累时间内接收地面反射的空时二维回波数据。
[0009] （2）多普勒滤波处理：
[0010] 按照下式，利用扩展因子化算法EFA中的降维转换矩阵，将机载雷达接收的空时二维回波数据，从阵元-脉冲域转化到阵元-多普勒域，得到一个待处理多普勒通道的原始回波数据：
[0011]
[0012] 其中，x表示待处理多普勒通道的原始回波数据，TEFA表示扩展因子化算法EFA中的降维转换矩阵，表示空时二维回波数据。
[0013] （3）获取广对称训练样本：
[0014] 3a）采用空域转换矩阵，对待处理多普勒通道的原始回波数据进行空域转换，得到空域转换数据；
[0015] 3b）采用时域转换矩阵，对待处理多普勒通道的原始回波数据进行时域转换，得到时域转换数据；
[0016] 3c）采用空时域转换矩阵，对待处理多普勒通道的原始回波数据进行空时域转换，得到空时域转换数据。
[0017] （4）估计广对称协方差矩阵：
[0018] 4a）利用最大似然方法，分别估计待处理多普勒通道的原始回波数据、空域转换数据、时域转换数据以及空时域转换数据对应的原始协方差矩阵、空域协方差矩阵、时域协方差矩阵以及空时域协方差矩阵；
[0019] 4b）对上述四个协方差矩阵取平均，将取平均后的结果定义为广对称协方差矩阵。
[0020] （5）获得自适应权矢量：
[0021] 按照自适应权矢量求解公式，用广对称协方差矩阵代替传统扩展因子化算法EFA中的原始协方差矩阵，获得自适应权矢量。
[0022] （6）自适应滤波处理：
[0023] 利用自适应权矢量，依次对待处理多普勒通道的原始回波数据按照距离单元分别，进行空时自适应滤波。
[0024] （7）输出结果：
[0025] 将自适应滤波处理以后的滤波结果输出。
[0026] 与现有技术相比，本发明具有以下优点：
[0027] 第一，本发明基于协方差矩阵的广对称特性，对多普勒域原始回波数据分别进行空域变换、时域变换、空时域变换，有效地改善了现有技术STAP方法，由于可用独立同分布训练样本数不足所引起的杂波协方差矩阵估计不准确的问题，使得本发明能够更充分地挖掘并使用有限训练样本中的杂波信息，提高回波数据的利用率，在获得杂波分布统计信息的同时，又可以很大程度上缓解由杂波分布的距离依赖性带来的样本非均匀性，在杂波极度非均匀的环境中具有更好的应用前景。
[0028] 第二，本发明利用扩展因子化算法EFA中的降维转换矩阵，对机载雷达接收的空时二维回波数据进行多普勒滤波处理的同时进行降维处理；估计广对称协方差矩阵采用原始协方差矩阵、空域协方差矩阵、时域协方差矩阵、空时域协方差矩阵取平均，有效地克服了现有技术的STAP方法，由于雷达系统的大自由度而造成的计算量过大和设备成本增加的问题，使得本发明能够在有效增加训练样本的同时，也不会带来过多的运算复杂度和结构复杂度。
附图说明
[0029] 图1为本发明流程图；
[0030] 图2为现有技术PD处理和本发明处理后的距离多普勒图；
[0031] 图3为现有技术EFA和本发明处理后的改善因子对比图；
[0032] 图4为现有技术EFA和本发明处理后的检测概率对比图。
具体实施方式
[0033] 下面结合附图对本发明作进一步的描述。
[0034] 参照图1，本发明的具体实施步骤如下：
[0035] 步骤1，接收空时二维回波数据。
[0036] 利用机载雷达天线的多个阵元，在相干积累时间内接收地面反射的空时二维回波数据 。
[0037] 步骤2，多普勒滤波处理。
[0038] 按照下式，利用扩展因子化算法EFA中的降维转换矩阵，将机载雷达接收的空时二维回波数据，从阵元-脉冲域转化到阵元-多普勒域，得到一个待处理多普勒通道的原始回波数据：
[0039]
[0040] 其中，x表示待处理多普勒通道的原始回波数据，TEFA表示扩展因子化算法EFA中的降维转换矩阵，表示空时二维回波数据；
[0041] 扩展因子化算法EFA中的降维转换矩阵TEFA由下式获得：
[0042]
[0043] 其中，TEFA表示EFA降维转换矩阵，FEFA表示待处理多普勒通道的EFA傅里叶变换矩阵EFA傅里叶变换矩阵，IN表示降维转换矩阵中的N维单位矩阵，N表示机载雷达天线的阵元个数。
[0044] 待处理多普勒通道的EFA傅里叶变换矩阵FEFA由下式获得：
[0045]
[0046] 其中，FEFA表示傅里叶变换矩阵，EFA表示该矩阵是EFA算法中的傅里叶变换矩阵，f-1,f0,f1表示以待处理多普勒f0为中心的三个相邻归一化多普勒频率，M表示机载雷达发射脉冲个数。
[0047] 步骤3，获取广对称训练样本。
[0048] 采用空域转换矩阵，对待处理多普勒通道的原始回波数据进行空域转换，得到空域转换数据：
[0049]
[0050] 其中，xs表示空域转换数据，s表示空域符号， 表示空域转换矩阵，I3表示空域转换矩阵中的3维单位矩阵， 表示Kronecker积运算符号，J表示空域转换矩阵中的N维置换矩阵，N表示机载雷达天线的阵元个数，x待处理多普勒通道的原始回波数据，*表示共轭运算符号。
[0051] 采用时域转换矩阵，对待处理多普勒通道的原始回波数据进行时域转换，得到时域转换数据：
[0052]
[0053] 其中，xt表示时域转换数据，t表示时域符号， 表示时域转换矩阵，U表示时域转换矩阵中的3维对角矩阵，IN表示时域转换矩阵中的N维单位矩阵，N表示机载雷达天线的阵元个数，x表示待处理多普勒通道的原始回波数据。
[0054] 采用空时域转换矩阵，对待处理多普勒通道的原始回波数据进行空时域转换，得到空时域转换数据：
[0055]
[0056] 其中，xst表示空时域转换数据，st表示空时域符号， 表示空时域转换矩阵，V表示空时域转换矩阵中的3维对角矩阵，J表示空域转换矩阵中的N维置换矩阵，N表示机载雷达天线的阵元个数，x表示待处理多普勒通道的原始回波数据。
[0057] 步骤4，估计广对称协方差矩阵。
[0058] 利用待处理多普勒通道的原始回波数据作训练样本，通过最大似然方法估计得到原始协方差矩阵：
[0059]
[0060] 其中，R表示原始协方差矩阵，L表示待处理多普勒通道的原始回波数据的长度，xi表示待处理多普勒通道的原始回波数据，i表示数据xi的索引号，H表示共轭转置；
[0061] 利用空域转换数据作训练样本，通过最大似然方法估计得到空域协方差矩阵：
[0062]
[0063] 其中，Rs表示空域协方差矩阵，Ls表示空域转换数据的长度，xs,i表示空域转换数据，s,i分别表示空域及数据索引号；
[0064] 利用时域转换数据作训练样本，通过最大似然方法估计得到时域协方差矩阵：
[0065]
[0066] 其中，Rt表示时域协方差矩阵，Lt表示时域转换数据的长度，xt,i表示时域转换数据，t,i分别表示时域及数据索引号；
[0067] 利用空时域转换数据作训练样本，通过最大似然方法估计得到空时域协方差矩阵：
[0068]
[0069] 其中，Rst表示空时域协方差矩阵，Lst表示空时域转换数据的长度，xst,i表示空时域转换数据，st,i分别表示空时域及数据索引号；
[0070] 对上述四个协方差矩阵取算数平均，将取平均后的结果定义为广对称协方差矩阵：
[0071]
[0072] 其中，RPer表示广对称协方差矩阵，Per表示广对称。
[0073] 步骤5，获得自适应权矢量。
[0074]
[0075] 其中，W表示自适应权矢量，μ表示由雷达系统参数设定的一个常数，S表示降维后的目标导向矢量，S的表达式为： TEFA表示扩展因子化算法EFA中的降维转换矩阵，st表示目标的时域导向矢量，st的表达式为： 表示目标的时域导向矢量，t表示时域符号，ss表示目标的空域导向矢量，ss的表达式为：
s表示空域符号，fd0表示归一化多普勒频率，d0表示多普勒频率编
号，fs表示归一化空间频率，N表示机载雷达天线阵元个数，M表示雷达发射脉冲个数，T表示转置运算符号。
[0076] 步骤6，自适应滤波处理。
[0077] 利用自适应权矢量W，依次对待处理多普勒通道的各距离单元原始回波数据x按照距离单元分别，进行空时自适应滤波。
[0078] 步骤7，输出结果。
[0079] 将自适应滤波处理以后的滤波结果输出。
[0080] 下面结合附图对本发明的效果做进一步说明。
[0081] 1、仿真条件：
[0082] 本发明的仿真实验在MATLAB7.11软件下进行的。在本发明的仿真实验中，为了真实地模拟在非均匀杂波环境中训练样本不足的情景，机载雷达采用非正侧视阵天线，机载雷达的天线采用10个阵元均匀排列的线阵，阵元间距为半个波长。在本发明的仿真实验中，使用的回波数据是根据林肯实验室J.Ward提出的杂波模型仿真产生，详细的系统参数参照下表。
[0083]
[0084] 参照附图2，对采用本发明与现有技术PD处理的机载雷达数据后杂波的剩余情况。
[0085] 附图2是现有技术PD处理和本发明处理后的距离多普勒图。附图2（a）为现有技术PD处理后的距离多普勒图，附图2（b）为本发明处理后的距离多普勒图。附图2（a）和附图2（b）的横轴表示多普勒通道序号，纵轴表示距离门序号，附图2（a）和附图2（b）中的白色区域为现有技术PD处理后的杂波剩余功率分布，黑色区域为现有技术PD处理后的噪声剩余功率分布。从附图2（a）可以看出较强的剩余杂波占据了较多的距离-多普勒单元，直接影响这些区域中目标的检测，从附图2（b）可以看出白色区域已大为减小，表明回波数据中的杂波分量得到了很好抑制。综上所述，与现有技术PD处理相比本发明处理的机载雷达数据的杂波剩余明显减弱，性能明显得到改善。
[0086] 参照附图3比较现有技术EFA和本发明处理后的雷达杂波改善因子情况。
[0087] 附图3是现有技术EFA和本发明的改善因子对比图，所选择比较的数据范围是300-350号距离门，系统自由度为30。附图3的横轴表示多普勒通道序号，纵轴表示改善因子，黑体实线表示的是作为比较参考的最优处理结果，实圈线和虚圈线分别表示现有技术EFA算法使用64个和32个训练样本下的平均改善因子曲线，实线、点线和虚线分别表示本发明使用64个、32个和8个训练样本下的平均改善因子曲线。从附图3中可以看出，当样本数较多，大于两倍系统自由度时，本发明处理结果要优于EFA，更接近最优处理结果；本发明在使用一倍自由度数目样本时的改善因子曲线基本与EFA在两倍自由度数目样本情况下的曲线重合，而本发明只使用8个训练样本即可达到EFA使用32个样本所能达到的杂波抑制效果。该仿真结果可以说明，本发明能够提高样本利用效率，增加有效的IID训练样本，从而使得协方差估计更为准确，取得更好的杂波抑制效果。
[0088] 参照附图4比较现有技术EFA和本发明处理后的雷达杂波检测概率情况。
[0089] 附图4是现有技术EFA和本发明处理后的检测概率对比图，所选择比较的数据范围是161-360号距离门、40-64号多普勒通道，系统自由度为30，虚警率设置为1%。附图4的横轴表示多普勒通道序号，纵轴表示检测概率，实圈线和虚圈线分别表示现有技术EFA算法使用60个和30个训练样本下的检测概率曲线，图中实线和虚线分别表示使用60个和30个训练样本下本发明的检测概率曲线。从附图4中可以看出，当训练样本充足，即系统自由度的两倍时，相较于EFA处理，本发明方法有约2dB的性能改善；而当样本数减少到一倍系统自由度时，本发明会优于EFA约15dB。该仿真结果表明，本发明可以在较少可用训练样本情况下，比如极度非均匀环境中，取得更好的目标检测效果。
