基于数据拟合的扩展因子化空时二维自适应处理方法

摘要

本发明公开了一种基于数据拟合的扩展因子化空时二维自适应处理方法，涉及雷达技术。其具体步骤为：步骤1，得到降维回波数据和空时导向矢量。步骤2，构造待检测距离单元的数据基矩阵；步骤3，获取目标阻塞以后的数据基矩阵。步骤4，获取辅助回波数据的拟合系数。步骤5，获取待检测距离单元的数据的最小拟合误差。步骤6，对待检测距离单元的数据的最小拟合误差进行单元平均恒虚警检测。用以解决传统STAP方法运算量巨大及样本需求苛刻的问题，在实测数据处理及对STAP的性能改善中具有很好的应用前景。

说明书

技术领域

本发明属于通信技术领域，涉及雷达技术，特别涉及一种基于数据拟合的扩展因子化空时二维自适应处理方法，用于机载预警雷达信号处理。 

背景技术

机载预警雷达的主要任务是在复杂杂波背景中探测目标，并对其进行定位跟踪，而对杂波进行有效抑制是提高机载预警雷达工作性能的核心手段。空时自适应处理(space-time adaptive processing，STAP)技术充分利用空域和时域信息，在对目标信号进行相干积累的同时，通过空时自适应处理滤除地面杂波，实现机载预警雷达对目标的有效检测，如美国的E2-D机载预警雷达就采用此技术。 

在实际应用中，STAP技术主要存在以下两方面问题：一方面，在非均匀的杂波环境中，要获得足够多的用于估计协方差矩阵的独立同分布(independent and identically distributed，IID)训练样本非常困难；另一方面，即使训练样本的需求得到满足，全空时处理计算量过大的问题会导致实时性难以保证。 

为解决上述问题，推动STAP技术更加实用化，人们提出了许多改进措施或方法。清华大学申请的发明专利“非均匀杂波环境下空时自适应处理方法”(专利申请号201010129723.3，公布号CN101819269A)公开了一种在非均匀杂波环境中超分辨估计杂波空时二维谱的超完备稀疏表示方法。该方法实现了在独立同分布样本数不足的情况下，利用单帧训练样本估计杂波协方差矩阵，从而避免强非均匀杂波环境对自适应处理效果的影响。但是，该方法仍然存在以下主要不足是：一是运算量大的问题，对杂波谱进行稀疏表示的超完备基数目未定，但远大于系统自由度，而实际中系统自由度通常成千上万，这样在每一距离单元样本的协方差矩阵重构过程中所需要的运算量非常大，不利于实时处理，从而影响到在实际工程应用中的效果。二是基失配问题，该方法中使用的基实际上是一组插值离散傅里叶DFT矢量，而实际的回波数据因为存在误差等非理想因素，不能用所给的基矢量稀疏表示；另外，该方法需要对空时平面进行离散化表示，这样会使得不位于离散网格点中心的信号向所有网格点泄露，从而破坏了回波数据的稀疏性。 

发明内容

针对上述在非均匀杂波环境中超分辨估计杂波空时二维谱的超完备稀疏表示方法的 不足，本发明提出了一种基于数据拟合的扩展因子化空时二维自适应处理(ExtendedFactored Approach STAP,EFA-STAP)方法，能够实现对目标信号的检测，并解决现有技术中运算量巨大和基失配的问题，同时相对于传统扩展因子化EFA-STAP方法，能够提高动目标的检测概率，降低了检测的虚警率。 

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

一种基于数据拟合的扩展因子化空时二维自适应处理方法，其特征在于，包括以下步骤： 

步骤1，利用机载预警雷达天线接收机载预警雷达的空时二维回波数据x，对空时二维回波数据x进行扩展因子化EFA降维，得到降维回波数据z和空时导向矢量s_z-EFA； 

步骤2，构造待检测距离单元的数据基矩阵Φ_k； 

步骤3，利用空时导向矢量s_z-EFA导出阻塞矩阵B；利用阻塞矩阵B和待检测距离单元的数据z_k，获取目标阻塞以后的辅助回波数据该待检测距离单元的数据z_k依次取自步骤1中的降维回波数据z；利用阻塞矩阵B和待检测距离单元的数据基矩阵Φ_k获取目标阻塞以后的数据基矩阵

步骤4，利用目标阻塞以后的数据基矩阵对辅助回波数据进行数据拟合，获取辅助回波数据的拟合系数

步骤5，利用待检测距离单元的数据基矩阵Φ_k和辅助回波数据的拟合系数对待检测距离单元的数据z_k进行数据拟合，获取待检测距离单元的数据z_k的最小拟合误差y_k；该待检测距离单元的数据z_k依次取自步骤1中的降维回波数据z； 

步骤6，对待检测距离单元的数据的最小拟合误差y_k进行单元平均恒虚警检测，并输出存在目标或不存在目标。 

上述技术方案的特点和进一步改进在于： 

(1)步骤1具体包括以下子步骤： 

1a)利用机载预警雷达天线，在相干积累时间内接收地面反射的空时二维回波数据x；空时二维回波数据x为MN维，其中M表示脉冲数，N表示机载预警雷达天线的阵元数； 

1b)根据下式得到扩展因子化EFA降维的降维转换矩阵P_EFA： 

$P_{\mathrm{EFA}}=P_t(f_{i-1},f_i,f_{i+1})\otimes I_N$

其中，P_t(f_i-1,f_i,f_i+1)表示由三个临近的多普勒滤波器组成的时域转换矩阵，f_i表示目标的归一化多普勒频率，f_i-1,f_i+1表示与目标相邻的两个归一化多普勒频率，i为小于等于M并且大于等于1的自然数，M表示脉冲数，I_N表示N×N的单位矩阵，表示克罗内克积，N表示机载预警雷达天线的阵元数，降维转换矩阵P_EFA的维数为MN×3N维； 

1c)根据降维转换矩阵P_EFA得到扩展因子化EFA降维后的降维回波数据z，表示为下式： 

$z=P_{\mathrm{EFA}}^{H}x$

其中，(·)^H表示共轭转置操作，根据降维转换矩阵P_EFA的维数为MN×3N维以及空时二维回波数据x为MN维可得降维回波数据z的维数为3N维； 

扩展因子化EFA降维后的空时导向矢量s_z-EFA，表示为下式： 

$s_{z-\mathrm{EFA}}={\left(\begin{array}{ccc}0& 1& 0\end{array}\right)}^T\otimes s_s\left(f_s\right)$

s_s(f_s)＝[1 exp(j2πf_s) exp(j2πf_s2) … exp(j2πf_s(N-1))]^T

其中，(·)^T表示转置操作，s_s(f_s)表示目标的空域导向矢量，f_s表示目标的归一化空间频率，N表示机载预警雷达天线的阵元数。 

(2)步骤2具体包括以下子步骤： 

数据基矩阵Φ_k表示为下式： 

Φ_k＝[z₁ z₂ … z_L]∈C^D×L

其中，Φ_k表示第k个待检测距离单元的数据基矩阵；k表示待检测距离单元的数据的序号，为不超过所有待检测距离单元总数的自然数；D表示降维回波数据z的维数，L表示临近距离单元的数据的个数，C^D×L表示D×L维的复矩阵空间，z_l(l＝1,2,…,L)表示临近距离单元的数据，临近距离单元表示待检测距离单元附近的距离单元，不包含两个用于保护目标的距离单元。 

(3)步骤3具体包括以下子步骤： 

3a)对空时导向矢量s_z-EFA的转置形式进行奇异值分解，得到左奇异矩阵U、右奇异矩 阵V和奇异值矩阵Λ；分解公式如下式： 

$s_{z-\mathrm{EFA}}^T=U·\Lambda ·V^H$

其中，U表示左奇异矩阵，V表示右奇异矩阵，Λ表示奇异值矩阵，(·)^T表示转置操作，(·)^H表示共轭转置操作； 

3b)阻塞矩阵B由右奇异矩阵V的第2到N列顺序排列构成，通过右奇异矩阵V和下式导出阻塞矩阵B； 

B＝[V(:,2:N)]^T

其中，V(:,2:N)表示右奇异矩阵V的所有行、第2到N列矩阵元素；(·)^T表示转置操作； 

3c)根据下式获取辅助回波数据和目标阻塞以后的数据基矩阵

${\bar{z}}_k=B{z}_k,{\bar{\Phi }}_k=B{\Phi }_k$

其中，z_k表示第k个待检测距离单元的数据，该待检测距离单元的数据依次取自步骤1中的降维回波数据z，k表示待检测距离单元的数据的序号，为不超过所有待检测距离单元总数的自然数；表示目标阻塞以后的辅助回波数据，表示目标阻塞以后的数据基矩阵，B为阻塞矩阵。 

(4)步骤4具体包括以下子步骤： 

4a)用目标阻塞以后的数据基矩阵对目标阻塞以后的辅助回波数据进行最小二乘拟合表示，表示为下式： 

$\underset{{\bar{\alpha }}_k}{\min }{\left|\right|{\bar{z}}_k-{\bar{\Phi }}_k{\bar{\alpha }}_k\left|\right|}_2$

其中，||·||₂表示求解矢量的2范数，表示辅助回波数据的拟合系数； 

4b)根据步骤4a)中的最小二乘拟合表示，求解辅助回波数据的拟合系数的最小二乘解表示为下式： 

${\bar{\alpha }}_k^{*}={\bar{\Phi }}_k^H{\left({\bar{\Phi }}_k{\bar{\Phi }}_k^H\right)}^{-1}{\bar{z}}_k$

其中，(·)^-1表示矩阵求逆操作，(·)^H表示共轭转置操作，表示目标阻塞以后的数据基矩阵，表示目标阻塞以后的辅助回波数据。 

(5)步骤5具体包括以下子步骤： 

5a)用待检测距离单元的数据基矩阵Φ_k，对待检测距离单元的数据z_k进行最小二乘拟合；用解辅助回波数据的拟合系数的最小二乘解计算得到优化约束边界η_k，具体通过下式实现： 

$\underset{{\alpha }_k}{\min }{\left|\right|z_k-{\Phi }_k{\alpha }_k\left|\right|}_2,\mathrm{subjectto}{\left|\right|{\alpha }_k\left|\right|}_2\le {\eta }_k$

其中，α_k表示待检测距离单元的数据的拟合系数，y_k表示最小的拟合误差矢量的2范数，η_k表示优化约束边界，

5b)通过求解上式得到待检测距离单元的数据的最优拟合系数

根据待检测距离单元的数据的最优拟合系数得到待检测距离单元的数据的最小拟合误差y_k，通过下式获取： 

$y_k={\left|\right|z_k-{\Phi }_k{\alpha }_k^{*}\left|\right|}_2$

其中，z_k表示第k个待检测距离单元的数据，该待检测距离单元的数据依次取自步骤1中的降维回波数据z，k表示待检测距离单元的数据的序号，为不超过所有待检测距离单元总数的自然数；Φ_k表示待检测距离单元的数据基矩阵；表示待检测距离单元的数据的最优拟合系数；||·||₂表示求解矢量的2范数。 

与现有技术相比，本发明具有突出的实质性特点和显著的进步。本发明与现有稀疏恢复STAP方法相比具有以下优点： 

(1)针对运算量的问题，现有稀疏恢复方法对杂波谱进行稀疏表示的超完备基数目未定，但远大于系统自由度，而实际中系统自由度通常成千上万，这样在每一距离单元样本的协方差矩阵重构过程中所需要的运算量非常大，不利于实时处理，从而影响到它的实际工程应用效果。本发明将STAP方法与EFA相结合，将原来MN维的优化问题分解为M个3N维的小优化问题来解决，大大减少了算法的计算复杂度，解决了稀疏恢复STAP方法中运算量大的问题。 

(2)针对基失配问题，现有稀疏恢复方法中使用的表示基实际上是一组插值DFT矢量，而实际数据因为存在误差等非理想因素，不能用所给基矢量稀疏表示，另外，该方法需要对空时平面进行离散化表示，这样会使得不位于离散网格点中心的信号向所有网格点泄露，从而破坏了数据的稀疏性。本发明直接使用实际数据样本作为基矩阵，而不是采用 某种确定的理想流型，从而规避了上述稀疏恢复STAP方法中基失配的问题。 

可见，现有技术中STAP通常在高维数据空间中进行机载雷达信号处理，这样就会大大增加运算复杂度和所需要的训练样本个数的问题，为了得到更好的性能，利用待检测距离单元临近距离单元的数据样本来优化拟合待检测距离单元杂波数据，从而改善自适应信号处理的杂波抑制性能，提高目标的检测概率。 

附图说明

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

图1为本发明流程图； 

图2为常规脉冲多普勒处理结果图； 

图3为检测概率随边界补偿量的变化曲线图； 

图4为检测概率随样本数的变化曲线图； 

图5为目标附近距离单元的数据的最小拟合误差的输出结果图；图5(a)为经过常规EFA方法之后一段距离单元范围的最小拟合误差的输出结果；图5(c)为经过经本发明方法处理之后一段距离单元范围的最小拟合误差的输出结果；图5(b)为经过常规EFA方法之后另一段距离单元范围的最小拟合误差的输出结果；图5(d)为经过本发明方法之后另一段距离单元范围的最小拟合误差的输出结果； 

图6为检测概率随目标功率的变化曲线图； 

图7为检测概率随虚警率的变化曲线关系图。 

具体实施方式

参照图1，说明本发明的一种基于数据拟合的扩展因子化空时二维自适应处理方法，本发明用于机载预警雷达信号处理，其具体实施步骤如下： 

步骤1，利用机载预警雷达天线接收机载预警雷达的空时二维回波数据x，对空时二维回波数据x进行扩展因子化EFA降维，得到降维回波数据z和空时导向矢量s_z-EFA； 

步骤1具体包括以下子步骤： 

1a)利用机载预警雷达天线，在相干积累时间内接收地面反射的空时二维回波数据x；空时二维回波数据x为MN维，其中M表示脉冲数，N表示机载预警雷达天线的阵元数； 

1b)根据下式得到扩展因子化EFA降维的降维转换矩阵P_EFA： 

$P_{\mathrm{EFA}}=P_t(f_{i-1},f_i,f_{i+1})\otimes I_N$

其中，P_t(f_i-1,f_i,f_i+1)表示由三个临近的多普勒滤波器组成的时域转换矩阵，f_i表示目标的 归一化多普勒频率，f_i-1,f_i+1表示与目标相邻的两个归一化多普勒频率，i为小于等于M并且大于等于1的自然数，M表示脉冲数，I_N表示N×N的单位矩阵，表示克罗内克积，N表示机载预警雷达天线的阵元数，降维转换矩阵P_EFA的维数为MN×3N维； 

1c)根据降维转换矩阵P_EFA得到扩展因子化EFA降维后的降维回波数据z，表示为下式： 

$z=P_{\mathrm{EFA}}^{H}x$

其中，(·)^H表示共轭转置操作，根据降维转换矩阵P_EFA的维数为MN×3N维以及空时二维回波数据x为MN维可得降维回波数据z的维数为3N维； 

扩展因子化EFA降维后的空时导向矢量s_z-EFA，表示为下式： 

$s_{z-\mathrm{EFA}}={\left(\begin{array}{ccc}0& 1& 0\end{array}\right)}^T\otimes s_s\left(f_s\right)$

s_s(f_s)＝[1 exp(j2πf_s) exp(j2πf_s2) … exp(j2πf_s(N-1))]^T

其中，(·)^T表示转置操作，s_s(f_s)表示目标的空域导向矢量，f_s表示目标的归一化空间频率，N表示机载预警雷达天线的阵元数。 

一般空时导向矢量s可表示为时域导向矢量与空域导向矢量的克罗内克积形式，即 f_d和f_s分别表示归一化的多普勒频率和空间频率，时域和空域导向矢量分别可表示为

步骤2，构造待检测距离单元的数据基矩阵Φ_k； 

数据基矩阵Φ_k表示为下式： 

Φ_k＝[z₁ z₂ … z_L]∈C^D×L

其中，Φ_k表示第k个待检测距离单元的数据基矩阵；k表示待检测距离单元的数据的序号，为不超过所有待检测距离单元总数的自然数；D表示降维回波数据z的维数，L表示临近距离单元的数据的个数，C^D×L表示D×L维的复矩阵空间，z_l(l＝1,2,…,L)表示临近距离单元的数据，临近距离单元表示待检测距离单元附近的距离单元，不包含两个用于保护目标的距离单元。 

步骤3，利用空时导向矢量s_z-EFA导出阻塞矩阵B；利用阻塞矩阵B和待检测距离单元的数据z_k，获取目标阻塞以后的辅助回波数据该待检测距离单元的数据z_k依次取自步骤1中的降维回波数据z；利用阻塞矩阵B和待检测距离单元的数据基矩阵Φ_k获取目标阻塞以后的数据基矩阵

步骤3具体包括以下子步骤： 

3a)对空时导向矢量s_z-EFA的转置形式进行奇异值分解，得到左奇异矩阵U、右奇异矩阵V和奇异值矩阵Λ；分解公式如下式： 

$s_{z-\mathrm{EFA}}^T=U·\Lambda ·V^H$

其中，U表示左奇异矩阵，V表示右奇异矩阵，Λ表示奇异值矩阵，(·)^T表示转置操作，(·)^H表示共轭转置操作； 

3b)阻塞矩阵B由右奇异矩阵V的第2到N列顺序排列构成，通过右奇异矩阵V和下式导出阻塞矩阵B； 

B＝[V(:,2:N)]^T

其中，V(:,2:N)表示右奇异矩阵V的所有行、第2到N列矩阵元素；(·)^T表示转置操作； 

3c)根据下式获取辅助回波数据和目标阻塞以后的数据基矩阵

${\bar{z}}_k=B{z}_k,{\bar{\Phi }}_k=B{\Phi }_k$

其中，z_k表示第k个待检测距离单元的数据，该待检测距离单元的数据依次取自步骤1中的降维回波数据z，k表示待检测距离单元的数据的序号，为不超过所有待检测距离单元总数的自然数；表示目标阻塞以后的辅助回波数据，表示目标阻塞以后的数据基矩阵，B为阻塞矩阵。 

步骤4，利用目标阻塞以后的数据基矩阵对辅助回波数据进行数据拟合，获取辅助回波数据的拟合系数

步骤4具体包括以下子步骤： 

4a)用目标阻塞以后的数据基矩阵对目标阻塞以后的辅助回波数据进行最小二乘拟合表示，表示为下式： 

$\underset{{\bar{\alpha }}_k}{\min }{\left|\right|{\bar{z}}_k-{\bar{\Phi }}_k{\bar{\alpha }}_k\left|\right|}_2$

其中，||·||₂表示求解矢量的2范数，表示辅助回波数据的拟合系数； 

4b)根据步骤4a)中的最小二乘拟合表示，求解辅助回波数据的拟合系数的最小二乘解表示为下式： 

${\bar{\alpha }}_k^{*}={\bar{\Phi }}_k^H{\left({\bar{\Phi }}_k{\bar{\Phi }}_k^H\right)}^{-1}{\bar{z}}_k$

其中，(·)^-1表示矩阵求逆操作，(·)^H表示共轭转置操作，表示目标阻塞以后的数据基矩阵，表示目标阻塞以后的辅助回波数据。 

步骤5，利用待检测距离单元的数据基矩阵Φ_k和辅助回波数据的拟合系数对待检测距离单元的数据z_k进行数据拟合，获取待检测距离单元的数据z_k的最小拟合误差y_k；该待检测距离单元的数据z_k依次取自步骤1中的降维回波数据z； 

步骤5具体包括以下子步骤： 

5a)用待检测距离单元的数据基矩阵Φ_k，对待检测距离单元的数据z_k进行最小二乘拟合；用解辅助回波数据的拟合系数的最小二乘解计算得到优化约束边界η_k，具体通过下式实现： 

$\underset{{\alpha }_k}{\min }{\left|\right|z_k-{\Phi }_k{\alpha }_k\left|\right|}_2,\mathrm{subjectto}{\left|\right|{\alpha }_k\left|\right|}_2\le {\eta }_k$

其中，α_k表示待检测距离单元的数据的拟合系数，y_k表示最小的拟合误差矢量的2范数，η_k表示优化约束边界

5b)通过求解上式得到待检测距离单元的数据的最优拟合系数

根据待检测距离单元的数据的最优拟合系数得到待检测距离单元的数据的最小拟合误差y_k，通过下式获取： 

$y_k={\left|\right|z_k-{\Phi }_k{\alpha }_k^{*}\left|\right|}_2$

其中，z_k表示第k个待检测距离单元的数据，该待检测距离单元的数据依次取自步骤1中的降维回波数据z，k表示待检测距离单元的数据的序号，为不超过所有待检测距离单 元总数的自然数；Φ_k表示待检测距离单元的数据基矩阵；表示待检测距离单元的数据的最优拟合系数；||·||₂表示求解矢量的2范数。 

步骤6，对待检测距离单元的数据的最小拟合误差y_k进行单元平均恒虚警检测，并输出存在目标或不存在目标。 

对待检测距离单元的数据的最小拟合误差y_k进行单元平均恒虚警检测(Cell-averaging constant false alarm rate，CA-CFAR)，即将最小拟合误差y_k与待检测距离单元的周围的距离单元的数据的最小拟合误差的平均值比较，根据最小拟合误差y_k与平均值比值大小来判断是否有存在目标，最后将存在目标或不存在目标输出。 

下面结合仿真实验对本发明的效果做进一步说明。 

(1)实验条件： 

本发明的实验在MATLAB7.11软件下进行的。在本发明的实验设计中，为了检验本发明方法在实际应用中的有效性，处理数据来自国外录取的实测数据MCARM数据集，其雷达工作在正侧阵模式下，阵面由2行11列的可用于自适应信号处理的空域阵元组成，其部分系统参数参照表1。为了试验方便，本发明抽取MCARM数据中的第一行11个阵元、前32个脉冲的回波数据进行验证。 

表1 

(2)实验结果比较 

附图2是常规脉冲多普勒处理结果图，也叫距离多普勒图。横轴表示多普勒通道序号，纵轴表示距离单元序号。在本发明中距离单元也称距离门。从图2中可以看出虽然该雷达工作在正侧阵情况下，但空时二维回波数据中的强杂波分量，受到地形反射率不同的调制及实际环境中各种虚假目标的存在，具有一定的非均匀性。 

附图3为检测概率随边界补偿量的变化曲线图。横轴表示边界补偿，纵轴表示检测概 率。边界补偿是指给求解步骤5中优化约束边界时所得的优化上界附加一个补偿量。从图3中可以看出，基本在0补偿量附近，检测概率可以取得最大值，这也就说明了本发明中步骤5计算优化约束边界是有效且最优的。 

附图4为检测概率随样本数的变化曲线图。在本发明中样本数指临近距离单元的数据的个数。横轴表示样本数，纵轴表示检测概率。本发明中经过数据EFA降维处理后，数据维数(即自由度)降为3N，其中N为阵元数，在本实验中3N＝33。从图中可以看出，检测概率在样本数超过100时基本达到最大值，这就说明，本发明中数据基矩阵可由待检测距离单元附近3倍自由度数目的临近距离单元的数据来构成，得到最优的检测性能。 

图5是目标附近距离单元的数据的最小拟合误差的输出结果图。在本发明中，目标附近距离单元的数据包括多组待检测距离单元的数据。横轴表示距离单元的序号，纵轴表示功率，用dB表示。图5中曲线表示目标多普勒通道内一段距离单元范围的最小拟合误差的输出结果，水平直线表示该段距离单元范围的最小拟合误差的平均输出值，圆圈表示目标距离单元的最小拟合误差输出值，图5中“差值”表示目标距离单元的最小拟合误差输出值与该段距离单元范围的最小拟合误差的平均输出值之间的差值。从图5中可以看出，经本发明处理后，目标输出更为明显，图5(a)与5(c)相比较在距离单元为150时功率提高超过15dB，图5(b)与5(d)相比较在距离单元为200时功率提高超过25dB，从而可以更有利于目标信号的检测。 

图6是检测概率随目标功率的变化曲线图。横轴表示人为设置的目标功率，纵轴表示该目标的检测概率。图6中比较了直接数据域稀疏恢复方法(D3SR)、传统扩展因子化方法(EFA)及本发明提出的基于数据拟合的扩展因子化方法(DF-EFA)。从图6中可以看出，本发明方法的检测性能远好于其他两种方法。 

图7是检测概率随虚警率的变化曲线关系图。横轴表示虚警率，纵轴表示检测概率。图7中比较了传统扩展因子化方法(EFA)和本发明提出的基于数据拟合的扩展因子化方法(DF-EFA)。从图7中可以看出，本发明方法相对于传统EFA方法，在一定虚警率下具有更高的检测概率，而在相同的检测概率下具有更低的虚警率。