说明书多阶段分解空时自适应信号处理方法
技术领域
本发明属于雷达信号处理技术领域，特别是一种为雷达在杂波环境下的运动目标检测方法。
背景技术
空时自适应处理STAP技术是阵列自适应技术的推广和延伸，它将脉冲和阵元采样的二维数据排列等效成一个长的一维阵列数据，在时空联合域进行自适应滤波，全维空时自适应处理技术能大大提高雷达在杂波下的动目标检测性能，理论上可以实现最优处理。但是由于雷达接收的空时数据维数往往很大，进行全空时域的自适应处理主要存在两个问题，一是运算量太大，硬件难以实现；二是估计杂波协方差矩阵所需的独立同分布参考距离单元太多，难以满足实际要求。
由于上述问题的存在，人们进行了大量的降维STAP研究，通过降低自适应自由度有效的降低相应的计算量和参考距离单元样本数，以提高自适应处理对非均匀环境杂波抑制的性能。已有的降维方法多种多样，如互谱法CSM、多级维纳滤波法MWF和辅助通道法ACP等。这些降维方法可以直接或等效成这样一个过程：通过降维算子有效地降低全维空时数据向量的长度，然后根据线性约束最小均方误差准则构建二次优化方程，从而解得最优权向量。通过降维处理，有效降低了运算量和估计协方差矩阵的独立同分布IID参考距离单元数据长度，从而提高自适应处理对非均匀环境杂波抑制的性能。对于线性约束二次方程的优化问题，自适应自由度和杂波协方差矩阵维数相等，保证自适应损益小于3dB的最少训练样本数为杂波协方差矩阵维数的2倍，与降维后的维数呈线性增长关系，同时，计算复杂度为O(M3)，与降维后的维数呈3阶关系，这里M为降维后的维数。为了改善杂波抑制性能，需要增加自适应自由度，即增大降维后的维数，但训练样本数和计算量也会相应增加。
发明内容
本发明的目的在于克服上述已有技术的缺陷，提出一种多阶段分解空时自适应信号处理方法，以提高自适应自由度，降低训练样本和计算复杂度，实现对雷达信号的自适应处理。
实现本发明目的的技术方案包括如下步骤：
1)对回波信号进行第一次滤波；
2)对回波信号同时进行空域和时域降维；
3)对降维后的回波信号再进行滤波，并计算滤波后回波信号的剩余杂波功率；
4)将计算出的剩余杂波功率与第一次滤波后的杂波功率进行比较，如果该剩余杂波功率远小于第一次滤波后的杂波功率，则进行步骤6)，反之进行步骤5)；
5)重复步骤2)～步骤4)，直到剩余杂波功率远小于第一次滤波后的杂波功率为止；
6)对降维后的滤波信号进行目标检测。
本发明与现有的技术相比较有如下特点：
1)本发明由于在降维处理过程中采用两个降维矩阵，因而能够分别在空域和时域同时降低数据矩阵的维数，同时提高了计算效率；
2)本发明由于分阶段的降维处理，保证了自适应自由度没有降低；
3)本发明由于各级滤波器分量之间存在正交关系，降低了对训练样本的需求；
4)本发明由于采用有限的滤波次数和迭代次数，可获得较优处理性能，同时具有很好的格型结构，实际处理中可以视情况进行简单的截断处理，降低了计算复杂度。
附图说明：
图1是本发明多阶段分解空时自适应信号处理流程图；
图2是本发明计算滤波器分量时采用的双迭代方法流程图；
图3是本发明仿真主杂波区的其中一个多普勒通道在各阶段迭代处理中改善因子IF随迭代次数的变化曲线图；
图4是本发明与已有方法对改善信杂噪比的改善因子对比图；
图5是本发明与已有方法的空时二维响应对比图；
图6是本发明采用的双迭代方法在各阶段处理中输出信杂噪比SCNR随迭代次数的变化曲线图；
图7是基于实测数据情况下，本发明与已有方法的滤波结果对比图。
具体实施方式
参照图1，本发明的具体步骤如下：
步骤1，对回波信号进行第一次滤波。
设雷达接收到的回波信号用X1表示，为了最小化对X1滤波后的剩余杂波功率，构建如下代价函数：
minu‾1,v‾1J(u~1,v~1)=E{||y1||2}=E{||u~1HX1v~1||2}---(1)]]>
式中，和为第1次滤波的两个滤波器分量，且满足u~1Hs1(ωs)s1H(ωt)v~1=1,]]>s1(ωs)和s1(ωt)分别为空域和时域导向矢量，E{·}表示期望值。
采用如图2所示的双迭代方法对代价函数(1)求解，得到第一次滤波的滤波器分量和具体实现步骤如下：
1a)选取初值
1b)计算v~1(k)=Rv,1-1(k)s1(ωt)/((u~1H(k-1)s1(ωs))s1H(ωt)Rv,1-1(k)s1(ωt)),]]>其中Rv,1(k)=E{d~1(k)d~1H(k)},]]>d~1(k)=X1Hu~1(k-1),]]>k表示迭代次数，k＝1，2，…；
1c)计算u~1(k)=Ru,1-1(k)s1(ωs)/||Ru,1-1(k)s1(ωs)||,]]>其中
Ru,1=E{c~1(k)c~1H(k)},]]>c~1(k)=X1v~1(k);]]>
1d)更新迭代次数k，重复步骤1b)和步骤1c)，直到||u~1(k)-u~1(k-1)||<ϵ(0<ϵ<1)]]>为止，且令u~1=u~1(k),]]>v~1=v~1(k),]]>其中ε表示阈值。
利用上述求得的滤波器分量和对回波信号进行滤波，得到滤波后的信号：y1=u~1HX1v~1.]]>
步骤2，对回波信号同时进行空域和时域降维。
对第p(p＝1，2，…)次滤波的滤波器分量和进行归一化处理：
分别定义归一化处理后的滤波器分量和的Householder矩阵：
式中，N为天线阵元个数，M为一个相干处理时间内的脉冲数，eN-p+1＝[1，0，…，0]T为(N-p+1)×1维的单位列向量，eM-p+1＝[1，0，…，0]T为(M-p+1)×1维的单位列向量，IN-p+1为(N-p+1)×(N-p+1)的单位矩阵，IM-p+1为(M-p+1)×(M-p+1)的单位矩阵，和分别为和的第一个元素，C表示复数域。
由Householder矩阵的性质可知，GpH的最后N-p列是的正交补空间，HpH的最后M-p列是的正交补空间，于是构建出关于和的阻塞矩阵：
Gp=[G‾p](2:N-p+1,1:N-p+1)]]>
Hp=[H‾p](2:M-p+1,1:M-p+1)]]>
式中，0N-p＝[0，0，…，0]T为(N-p)×1维的零向量，0M-p＝[0，0，…，0]T为(M-p)×1维的零向量，IN-p为(N-p)×(N-p)的单位矩阵，IM-p为(M-p)×(M-p)的单位矩阵，eN-p＝[1，0，…，0]T为(N-p)×1维的单位列向量，eM-p＝[1，0，…，0]T为(M-p)×1维的单位列向量。
Gp和Hp即为第p次降维的降维矩阵，进而对回波信号进行降维处理，得到第p次降维后的回波信号：Xp+1=GpXpHpH,]]>式中Xp为第p-1次降维后的回波信号。
步骤3，对降维后的回波信号再进行滤波，并计算滤波后回波信号的剩余杂波功率。
对于第p次降维后的回波信号Xp+1，为了最小化p+1次滤波后的剩余杂波功率，再构建如下代价函数：
minu‾p+1,v‾p+1J(u~p+1,v~p+1)=E{||Σq=1pyq+yp+1||2}=E{||Σq=1pu~qHXqv~q+u~P+1HXp+1v~p+1||2}---(2)]]>
式中，和为第p+1次滤波的两个滤波器分量，并满足u~p+1Hsp+1(ωs)sp+1H(ωt)v~P+1=0,]]>u~p+1Hu~q=0]]>和v~p+1Hv~q=0,(q=1,2,...,p),]]>sn+1(ωs)和sn+1(ωt)分别为第p次降维后的空域和时域导向矢量。
采用如图2所示的双迭代方法对代价函数(2)求解，得到第p+1次滤波的滤波器分量和具体实现步骤如下：
3a)选取初值
3b)计算v~p+1(k)=-Rv,p+1-1(k)[bv,p+1(k)+μp+1(k)sp+1(ωt)],]]>其中
Rv,p+1(k)=E{d~p+1(k)d~p+1H(k)},]]>bv,p+1(k)=E{d~p+1(k)Σq=1pyq},]]>
d~p+1(k)=Xp+1Hu~p+1(k-1),]]>μp+1(k)=-(sp+1H(ωt)Rv,p+1-1(k)bv,p+1(k))/(sp+1H(ωt)Rv,p+1-1(k)sp+1(ωt));]]>
3c)计算u~p+1(k)=-Ru,p+1-1(k)bu,p+1(k),]]>其中：
Ru,p+1(k)=E{c~p+1(k)c~p+1H(k)},]]>bu,p+1(k)=E{c~p+1(k)(Σq=1pyq*)},]]>c~p+1(k)=Xp+1v~p+1(k);]]>
3d)更新迭代次数k，重复步骤3b)和步骤3c)，直到||u~p+1(k)-u~p+1(k-1)||/||u~p+1(k)||<ϵ(0<ϵ<1)]]>为止，且令u~p+1=u~p+1(k),]]>v~p+1=v~p+1(k).]]>
利用上述求得的滤波器分量和对回波信号进行滤波处理，得到第p+1次滤波后的信号：yp+1=u~p+1HXp+1v~p+1.]]>
根据第p+1次滤波后的信号yp+1，计算第p+1次滤波后的剩余杂波功率：P＝E{||yp+1||2}。
步骤4，将计算出的剩余杂波功率与第一次滤波后的杂波功率进行比较，如果该剩余杂波功率远小于第一次滤波后的杂波功率，则进行步骤6)，反之进行步骤5)；
步骤5，重复步骤2)～步骤4)，直到剩余杂波功率远小于第一次滤波后的杂波功率为止；
步骤6，对滤波后的信号y＝y1+y2+…+yr进行目标检测，式中r为满足要求时滤波的次数。
本发明的效果可以通过以下实验进一步说明：
(一)仿真条件：
1.基于计算机仿真数据的仿真条件。
雷达系统采用16×16的正侧视面阵阵列结构，发射时俯仰向和方位向均加30dB的Chebychev权，雷达工作波长λ＝0.23m，阵元间距d＝0.115m，接收数据首先微波合成为16个子天线的均匀线阵，载机的飞行速度为V＝115m/s，飞行高度为6km，脉冲重复频率为fr＝2000Hz，一个相干处理时间内的脉冲数为32，输入杂噪比为60dB，波束指向偏离阵面法向90°，阵元随机幅相误差为3％，距离单元数为496。
2.基于实测数据的仿真条件。
采用基本反映机载雷达工作真实环境的一批L波段相控阵机载雷达实测数据研究本发明对杂波抑制的性能，主要参数为：雷达工作波长λ＝0.2419m，载机速度v＝100.2m/s，偏航角阵元间距d＝0.1092m，脉冲重复频率fr＝1984Hz，仰角θ＝4°。
实验中主要处理320-410号距离门中的11个子阵的前32个脉冲数据，为了便于性能分析，在第375号距离单元注入一信杂噪比为-43dB的动目标信号，目标信号方位角为90°，归一化多普勒频率为0.125，位于主杂波区。
(二)仿真结果：
1.基于计算机仿真数据的仿真结果。
如图3给出了主杂波区的其中一个多普勒通道在各阶段迭代处理中改善因子随迭代次数的变化曲线。其中图3(a)为第一次滤波时的改善因子随迭代次数的变化曲线，6次迭代后实现收敛；图3(b)为第二次滤波时的改善因子随迭代次数的变化曲线，6次迭代后实现收敛，此时与第一次滤波相比有3.12dB的性能增益；图3(c)为第三次滤波时的改善因子随迭代次数的变化曲线，6次迭代后实现收敛，此时与第二次滤波相比有0.53dB的性能增益；图3(d)为第四次滤波时的改善因子随迭代次数的变化曲线，6次迭代后实现收敛，此时与第三次滤波相比仅有0.08dB的性能增益，表明在经过3次滤波处理后，剩余杂波已经很小，再进行滤波处理性能改善有限。从3图中可以看出，在每一阶段的迭代处理中，均有快速的收敛性。
图4给出了本发明在不同阶段处理的改善因子对比，以及与已有的因子化方法FA和扩展的因子化方法EFA的改善因子的对比，从图4中可以看出，本发明随着r的增大性能不断提高，但r增到3时，该方法可以得到很好的性能，有助于改善和提高主杂波区的性能，再过多的增加r也将毫无意义的增加计算复杂度。同时，通过与现有FA和EFA方法的改善因子曲线的比较可知，本发明方法的性能优于FA和EFA方法，特别是本发明方法能够有效提高主瓣杂波区的性能，有利于对慢速运动目标的检测。
本发明在r＝1，r＝3时的空时二维响应，以及全维STAP处理方法的空时二维自适应响应如图5所示。其中：图5(a)为全维STAP处理方法的自适应二维响应图，可以看到，滤波器在目标信号处形成无损益的主瓣，并且在杂波所对应的斜带上均形成了很好的阻带，能有效地抑制杂波，此时改善因子为82.45dB。图5(b)和图5(c)分别为本发明在r＝1和r＝3时的空时二维响应图，可以看出，随着r的增大，滤波器在杂波所对应的斜带上形成了很好的阻带，表明r＝3时能够实现对杂波的良好抑制，此时改善因子为82.37dB，与全维STAP处理方法相比仅有0.08dB的性能损失。
2.基于实测数据的仿真结果。
本发明分别在第一次、第二次和第三次滤波时输出SCNR随迭代次数的变化曲线如图6所示，其中图6(a)为第一次滤波时输出SCNR随迭代次数的变化曲线，经过10次迭代可以实现收敛，此时输出SCNR为30.36dB；图6(b)为在第一次滤波的基础上进行第二次滤波时输出SCNR随迭代次数的变化曲线，经过9次迭代可以实现收敛，此时输出SCNR为35.11dB，与第一次相比有4.75dB的性能改善；图6(c)为在第二次滤波的基础上进行第三次滤波时输出SCNR随迭代次数的变化曲线，经过10次迭代可以实现收敛，此时输出SCNR为36.20dB，与第二次相比仅有1.09dB的性能改善，此时剩余杂波已经很小，再进行滤波处理性能改善有限。以上结果表明，本发明在处理实测数据时同样具有很好的收敛性能。
图7给出了本发明在r＝3时的滤波结果，以及FA方法和EFA方法的滤波结果，几种方法都能有效滤除杂波，将动目标从残余杂波中提取出来。图7(a)为本发明在r＝3时的滤波结果，目标信号高出剩余杂波噪声峰值功率值为28.77dB；图7(b)为FA方法的滤波结果，目标信号高出剩余杂波噪声峰值功率值为11.52dB；图7(c)为EFA方法的滤波结果，目标信号高出剩余杂波噪声峰值功率值为17.66dB。以上结果表明，本发明的杂波抑制性能要优于FA方法和EFA方法。
从图6和图7可以得到，本发明在提高输出SCNR的同时，进一步提高了目标信号高出剩余噪声杂波峰值功率值，避免了大的剩余杂波对检测性能的影响，使得该方法在虚警概率一定的情况下可以获得较高的检测概率。