CN113376606A - 沿杂波脊快速收敛稀疏贝叶斯的杂波抑制方法 - Google Patents

沿杂波脊快速收敛稀疏贝叶斯的杂波抑制方法 Download PDF

Info

Publication number
CN113376606A
CN113376606A CN202110559117.3A CN202110559117A CN113376606A CN 113376606 A CN113376606 A CN 113376606A CN 202110559117 A CN202110559117 A CN 202110559117A CN 113376606 A CN113376606 A CN 113376606A
Authority
CN
China
Prior art keywords
clutter
matrix
vector
space
coefficient
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202110559117.3A
Other languages
English (en)
Other versions
CN113376606B (zh
Inventor
王彤
唐昊
刘程
王德根
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Xidian University
Original Assignee
Xidian University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Xidian University filed Critical Xidian University
Priority to CN202110559117.3A priority Critical patent/CN113376606B/zh
Publication of CN113376606A publication Critical patent/CN113376606A/zh
Application granted granted Critical
Publication of CN113376606B publication Critical patent/CN113376606B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/02—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
    • G01S7/41—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00 using analysis of echo signal for target characterisation; Target signature; Target cross-section
    • G01S7/414—Discriminating targets with respect to background clutter
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/02—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
    • G01S7/41—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00 using analysis of echo signal for target characterisation; Target signature; Target cross-section
    • G01S7/418—Theoretical aspects
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Engineering & Computer Science (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

本发明属于雷达技术领域，公开了一种沿杂波脊快速收敛稀疏贝叶斯的杂波抑制方法，通过先验知识和直接数据域算法联合估计杂波脊，最大程度上获得杂波脊的位置，确保不遗落杂波信息，结合杂波脊对超完备字典矩阵进行网格点选取重构字典矩阵，利用快速收敛稀疏贝叶斯算法对空时谱进行稀疏恢复，在容许杂波抑制性能下降的范围内减少了运算量，提高了运算速度，更适合在实时自适应处理中采用。

Description

沿杂波脊快速收敛稀疏贝叶斯的杂波抑制方法
技术领域
本发明属于雷达技术领域，具体涉及一种沿杂波脊快速收敛稀疏贝叶斯的杂波抑制方法。
背景技术
随着压缩感知技术的兴起，稀疏表示技术在雷达技术领域中也展现出了巨大的潜力，该技术的核心是利用信号的稀疏特性，通过在基矩阵中使用最少的取样点实现对信号的线性表示，能够在小样本的条件下实现对杂波谱的稀疏恢复。采用稀疏恢复技术去抑制杂波通常包含两个步骤：一是建立字典矩阵对训练样本数据进行稀疏表示；二是根据稀疏表示构造杂波加噪声协方差矩阵(CCM)，并利用CCM设计STAP滤波器。大多数的稀疏恢复算法需要一个或多个精确的系统参数，不精确的系统参数会严重影响稀疏恢复的准确性，进而降低杂波抑制性能。
稀疏贝叶斯收敛算法是由Tipping在2001年提出的一种稀疏恢复算法，算法的核心思想是假设稀疏系数矢量服从某一稀疏先验分布，联合先验信息和观测得到的数据，利用最大后验估计系数向量或系数矩阵。将稀疏贝叶斯收敛算法应用到STAP算法中，将空时谱等距划分网格点，网格点之间相互独立，网格点对应的稀疏系数服从复高斯先验分布。相比较传统的稀疏算法，稀疏贝叶斯收敛算法具有其他算法所不具备的优势：1.在无噪声情况下，稀疏贝叶斯收敛算法获得的真实解为最稀疏的解，不需要满足一些严格的条件，但大多数稀疏恢复算法做不到这点；2.当字典元素间的相关性很强时，稀疏贝叶斯收敛算法仍适用于这种条件。稀疏贝叶斯收敛算法本质上是迭代加权L1最小化算法，需要重复迭代取得最稀疏的解，计算成本较大。稀疏贝叶斯收敛算法利用稀疏解的结构信息不停地迭代来提高稀疏恢复的精度，直至收敛才停止迭代。正因如此，稀疏贝叶斯收敛算法受结构信息的影响比较大，先验参数的准确性直接关系到稀疏贝叶斯收敛算法的精度，先验参数不精确的直接结果会极大地降低稀疏贝叶斯收敛算法的稀疏恢复精度。
发明内容
为了解决上述稀疏贝叶斯收敛算法先验参数不精确导致稀疏恢复精度下降的问题，本发明的目的是提出一种沿杂波脊快速收敛稀疏贝叶斯杂波抑制方法，可利用有限的训练样本进行空时处理实现杂波抑制，可以显著降低训练样本数量需求，减少了运算量，提高了运算速度，更适合在实时自适应处理中采用。
本发明的技术原理为：通过先验知识和直接数据域算法联合估计杂波脊，最大程度上获得杂波脊的位置，确保不遗落杂波信息，结合杂波脊对超完备字典矩阵进行网格点选取重构字典矩阵，利用快速收敛稀疏贝叶斯算法对空时谱进行稀疏恢复，在容许杂波抑制性能下降的范围内减少了运算量，提高了运算速度，更适合在实时自适应处理中采用。
为了达到上述目的，本发明采用以下技术方案予以解决。
沿杂波脊快速收敛稀疏贝叶斯的杂波抑制方法，包括以下步骤：
步骤1，通过先验知识和直接数据域算法联合估计杂波脊位置，根据所述杂波脊位置对超完备字典矩阵进行网格点选取重构字典矩阵D；
步骤2，根据所述字典矩阵D得到正侧视阵构型下雷达接收到的杂波数据x；
步骤3，求解所述雷达接收到的杂波数据x的似然函数，通过稀疏贝叶斯迭代方法来先求解系数向量的后验概率密度函数，从而得到均值和方差，再采用期望最大化算法得到超参数向量的估计值，最终通过估计的超参数向量得到系数向量估计值；
步骤4，根据所述系数向量估计值得到杂波加噪声协方差矩阵
Figure BDA0003078248520000034
计算空时自适应滤波处理在稀疏贝叶斯迭代方法下的权向量以及在正侧视阵下雷达杂波抑制改善因子，以此来衡量正侧视阵雷达杂波谱的稀疏恢复效果。
本发明技术方案的特点和进一步的改进为：
(1)步骤1具体为：
子步骤1.1，当先验知识已知时，可以直接求得杂波脊位置，具体如下：
正侧阵时，天线锥角余弦值cosψ和多普勒频率fd间的关系为：
Figure BDA0003078248520000031
其中，v表示载机飞行速度，λ表示载波波长，ψ表示天线锥角，fr表示发射脉冲的重复频率，fd是杂波散射单元对应的多普勒频率；
斜侧阵构型时，天线锥角余弦值cosψ和多普勒频率fd间的关系为：
Figure BDA0003078248520000032
其中，
Figure BDA0003078248520000033
表示俯仰角；
通过所述多普勒频率fd与天线锥角ψ的关系可以从超完备字典矩阵中选取出对应的点，即为杂波脊对应的网格点，结合其在矩阵中相邻的点重新构建字典矩阵D；
子步骤1.2，当先验知识未知时，可以通过直接数据域算法求得杂波脊位置，具体如下：
目标信号矩阵由空域导向矢量和时域导向矢量所构成，目标信号沿行方向的相位差为
Figure BDA0003078248520000041
ψS0是目标的空间锥角；目标信号沿列方向的相位差为
Figure BDA0003078248520000042
fd0是目标的多普勒频率；目标信号沿对角线方向的相位差为
Figure BDA0003078248520000043
则空域导向矢量SS(ψS0)、时域导向矢量ST(fd0)及空时导向矢量S为：
Figure BDA0003078248520000044
Figure BDA0003078248520000045
Figure BDA0003078248520000046
其中，(i)T表示转置，j表示复数的虚部，N表示天线平面的阵元数，M表示发射脉冲数；沿空域、时域和空时域做两阵元对消，滤除目标信号并增加样本，分别得到(N-1)×M维空域低维度矩阵XS、N×(M-1)维时域低维度矩阵XT和(N-1)×(M-1)维空时域低维度矩阵XST：
Figure BDA0003078248520000047
Figure BDA0003078248520000048
Figure BDA0003078248520000049
假设空域每次的变化量为Nm，与之对应的时域每次的变化量为Mm，则所述空域低维度矩阵XS经空域滑窗取样，得到(N-Nm)×(M-Mm+1)个样本；所述时域低维度矩阵XT经时域滑窗取样，得到(N-Nm+1)×(M-Mm)个样本；所述空时域低维度矩阵XST经空时域联合滑窗取样，得到(N-Nm)×(M-Mm)个样本；因此，总共可以得到样本数为：
L0＝2((N-Nm)×(M-Mm+1)+(N-Nm+1)×(M-Mm)+(N-Nm)×(M-Mm))
训练样本采用上述L0个样本，用
Figure BDA0003078248520000051
表示，l0＝1,2,…,L0，则协方差矩阵的最大似然估计为：
Figure BDA0003078248520000052
其中，(i)H表示共轭转置，Vec(i)表示对矩阵堆叠成向量的操作；
求得所述协方差矩阵的最大似然估计
Figure BDA0003078248520000053
后，画出Capon谱即可得出杂波脊的位置，通过杂波脊的位置在超完备字典矩阵中对应的网格点结合其相邻网格点重建字典矩阵D。
(2)步骤2具体为：
雷达接收到的杂波数据x为：
x＝αi⊙Si+Di-Αi-+N0
其中，x＝[x1 x2 … xL]，L表示正侧视阵雷达接收到的杂波包含的距离门总个数；
Figure BDA0003078248520000054
αi是第i个网格点对应的信号分量的幅度，αip是第p个快拍第i个网格点对应的信号分量的幅度，1≤p≤L，
Figure BDA0003078248520000055
表示Kronecker积；⊙表示Hadamard积；Si＝[si si … si]NM×L，si表示字典矩阵D中第i个网格点对应的导向矢量；
Figure BDA0003078248520000056
表示从字典矩阵D中剔除掉si后的新的字典矩阵，Nd为归一化多普勒频率的细分格数，Ns为归一化空间频率的细分格数；
Figure BDA0003078248520000061
表示从系数矩阵Α中移除的第i个元素αi的新的系数矩阵；N0表示服从零均值复高斯分布的噪声向量。
(3)步骤3具体包含以下子步骤：
子步骤3.1，根据系数矩阵A和对角矩阵Γ，可以推导出系数矩阵A的先验概率密度函数：
Figure BDA0003078248520000062
其中，Γ＝diag(γ)为对角矩阵，
Figure BDA0003078248520000063
对应字典矩阵D中网格点系数的方差；
子步骤3.2，所述雷达接收到的杂波数据x的似然函数为：
Figure BDA0003078248520000064
其中，
Figure BDA0003078248520000065
为噪声功率，I表示单位向量，sk表示字典矩阵D中第k个网格点对应的导向矢量，γk表示对应字典矩阵D中第k个网格点系数的方差；
对系数矩阵A进行积分得到边缘概率密度函数
Figure BDA0003078248520000066
去估计超参数Γ和
Figure BDA0003078248520000067
使所述边缘概率密度函数
Figure BDA0003078248520000068
取最大，所述边缘概率密度函数
Figure BDA0003078248520000069
为：
Figure BDA00030782485200000610
直接求解上式比较复杂，因为网格点对应的系数之间相互独立，所以可以利用以下的表达式去估计超参数γi：
Figure BDA00030782485200000611
进而求得超参数Γ和
Figure BDA00030782485200000612
其中M表达式为：
Figure BDA0003078248520000071
则超参数γi的表达式为：
Figure BDA0003078248520000072
子步骤3.3，由于上式求解困难，因此引入贝叶斯定理可以得到在γi固定的条件下，系数ai的后验概率密度函数为：
Figure BDA0003078248520000073
其中，εi表示ai的方差，μi表示ai的均值，两者的表达式为：
Figure BDA0003078248520000074
Figure BDA0003078248520000075
可以得到系数ai的估计值：
Figure BDA0003078248520000076
则系数矩阵Α的估计值为：
Figure BDA0003078248520000077
子步骤3.4，采用期望最大化算法求得超参数
Figure BDA0003078248520000078
和
Figure BDA0003078248520000079
其中t表示迭代的次数；
Figure BDA00030782485200000710
表示γi第t+1次迭代估计值，
Figure BDA00030782485200000711
表示
Figure BDA00030782485200000712
第t+1次迭代估计值。
(4)子步骤3.4中，所述期望最大化算法分为求期望和最大化期望两步，在求期望中，将未知系数ai作为隐变量，在已知超参数
Figure BDA00030782485200000713
和
Figure BDA00030782485200000714
的条件下，计算
Figure BDA00030782485200000715
的后验概率密度函数
Figure BDA00030782485200000716
在最大化期望中，已知
Figure BDA00030782485200000717
的后验概率密度函数
Figure BDA00030782485200000718
时求出超参数
Figure BDA0003078248520000081
和
Figure BDA0003078248520000082
当迭代次数超过T次或||γ(t+1)-γ(t)||2/||γ(t+1)||2<η时终止迭代；其中，η表示阈值。
(5)步骤4中，所述杂波加噪声协方差矩阵
Figure BDA0003078248520000083
为：
Figure BDA0003078248520000084
其中，I表示单位矩阵，um,n表示
Figure BDA0003078248520000085
的第m行第n列的元素，1≤m≤K，1≤n≤L，K表示挑选字典元素的序列号集合；δ≥1为实常数，
Figure BDA0003078248520000086
为估计的噪声功率；
空时自适应滤波处理在稀疏贝叶斯迭代方法下的权向量为：
Figure BDA0003078248520000087
其中，s(fdt,fst)表示目标的空时导向矢量，fdt表示正侧视阵雷达接收到的杂波中包含的动目标归一化多普勒频率；fst表示正侧视阵雷达接收到的杂波中包含的动目标归一化空间频率；
空时自适应滤波处理后的非正侧视阵雷达杂波抑制改善因子为：
Figure BDA0003078248520000088
其中，tr(i)表示矩阵的迹。
与现有技术相比，本发明的有益效果为：
本发明利用先验知识和直接数据域算法联合估计杂波脊，重构字典矩阵，利用快速收敛稀疏贝叶斯算法估计雷达杂波谱，有效地提高了稀疏恢复的准确性，可以显著降低空时处理中训练样本数量需求，并且容许杂波抑制性能下降的范围内减少了运算量，提高了运算速度。与现有的稀疏贝叶斯算法相比，本发明给出的杂波抑制方法原理简单，便于工程实现，更适合实时自适应处理。
附图说明
下面结合附图和具体实施例对本发明做进一步详细说明。
图1是本发明的沿杂波脊快速收敛稀疏贝叶斯的杂波抑制方法的流程示意图；
图2(a)是利用OMP算法恢复的Capon谱；图2(b)是利用CVX算法恢复的Capon谱；图2(c)是利用SMI算法恢复的Capon谱；图2(d)是利用本发明的RM-SBL方法恢复的Capon谱；
图3是各算法在多快拍时的杂波抑制改善因子对比图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
下面结合附图对本发明的实施例及效果作进一步详细描述。
步骤1，通过先验知识和直接数据域算法联合估计杂波脊位置，根据所述杂波脊位置对超完备字典矩阵进行网格点选取重构字典矩阵D。
具体地，(a)当先验知识如载机速度、波长、天线锥角、重频等信息足够时，可以直接求得杂波脊。正侧阵时，天线锥角余弦值cosψ和多普勒频率fd间的关系为：
Figure BDA0003078248520000101
其中，v表示载机飞行速度，λ表示载波波长，ψ表示天线锥角，fr表示发射脉冲的重复频率，fd是该杂波散射单元对应的多普勒频率。
斜侧阵构型时，天线锥角余弦值cosψ和多普勒频率fd间的关系为：
Figure BDA0003078248520000102
其中，
Figure BDA0003078248520000103
表示俯仰角；
通过上述公式中多普勒频率与天线锥角的关系可以从超完备字典矩阵中选取出对应的点，即为杂波脊对应的网格点，结合其在矩阵中相邻的点重新构建字典矩阵D。
(b)当先验知识未知时，可以通过直接数据域算法求得杂波脊位置：
目标信号矩阵S由空域导向矢量和时域导向矢量所构成，其沿行方向(空域)的相位差为
Figure BDA0003078248520000104
其中ψS0是目标的空间锥角，目标沿列方向(时域)的相位差为
Figure BDA0003078248520000105
fd0是目标的多普勒频率，沿对角线方向(空时域)的相位差为
Figure BDA0003078248520000106
则空域导向矢量，时域导向矢量及空时导向矢量可写成：
Figure BDA0003078248520000107
Figure BDA0003078248520000108
Figure BDA0003078248520000109
其中，(·)T表示转置，j表示复数的虚部，N表示天线平面的阵元数，M表示发射脉冲数，沿空域、时域和空时域做两阵元(两脉冲)对消，目的是滤除目标信号并增加样本，得到三个低维度矩阵：
Figure BDA0003078248520000111
Figure BDA0003078248520000112
Figure BDA0003078248520000113
矩阵XS、XT和XST分别为(N-1)×M维、N×(M-1)维和(N-1)×(M-1)维；若假设空域每次的变化量为Nm，与之呼应的时域每次的变化量为Mm；则矩阵XS、XT和XST分别经空域滑窗取样、时域滑窗取样、空时域联合滑窗取样可以得到(N-Nm)×(M-Mm+1)、(N-Nm+1)×(M-Mm)、(N-Nm)×(M-Mm)个样本，因此总共可以得到样本数为：
L0＝2((N-Nm)×(M-Mm+1)+(N-Nm+1)×(M-Mm)+(N-Nm)×(M-Mm))
训练样本此时可以用这L0个样本，用
Figure BDA0003078248520000114
表示，l0＝1,2,…,L0，此时协方差矩阵的最大似然估计为：
Figure BDA0003078248520000115
其中，(·)H表示共轭转置，Vec(·)表示对矩阵堆叠成向量的操作，求得
Figure BDA0003078248520000116
后，画出Capon谱即可得出杂波脊的位置，通过杂波脊的位置再超完备字典矩阵中对应的网格点结合其相邻网格点重建字典矩阵D。
步骤2，根据所述字典矩阵D得到正侧视阵构型下雷达接收到的杂波数据x。
具体的，雷达接收到的杂波数据x为：
x＝αi⊙Si+Di-Αi-+N0
其中，x＝[x1 x2 … xL]，L表示正侧视阵雷达接收到的杂波包含的距离门总个数(即样本数)，xq表示第q个距离门的回波数据，1≤q≤L；Si＝[si si … si]NM×L，si表示字典矩阵D中第i个网格点对应的导向矢量，
Figure BDA0003078248520000121
αip是第p个快拍第i个网格点对应的信号分量的幅度，1≤p≤L，⊙表示Hadamard积，
Figure BDA0003078248520000122
表示Kronecker积，
Figure BDA0003078248520000123
表示从字典矩阵D中剔除掉si后的新的字典矩阵，Nd为归一化多普勒频率的细分格数，Ns为归一化空间频率的细分格数，设
Figure BDA0003078248520000124
表示从系数矩阵Α中移除的第i个元素αi的新的系数矩阵；N0表示服从零均值复高斯分布的噪声向量。
步骤3，求解所述雷达接收到的杂波数据x的似然函数，通过稀疏贝叶斯迭代方法来先求解系数向量的后验概率密度函数，从而得到均值和方差，再采用期望最大化算法得到超参数向量的估计值，最终通过估计的超参数向量得到系数向量估计值。
具体的，求解步骤2中杂波数据x的似然函数，进而得到系数向量，因为直接通过似然函数计算系数向量的计算量十分复杂，因此引入贝叶斯定理通过稀疏贝叶斯迭代方法来先求解系数向量的后验概率密度函数，得到均值和方差的表达式，再采用期望最大化算法得到超参数向量的估计值，最终通过超参数向量来估计系数向量；其中，稀疏贝叶斯迭代方法的子步骤为：
子步骤3.1，通过步骤2中的系数矩阵A联合对角矩阵Γ可以推导出系数矩阵A的先验概率密度函数：
Figure BDA0003078248520000125
其中Γ＝diag(γ)为对角矩阵，
Figure BDA0003078248520000126
对应字典矩阵D中网格点系数的方差。
子步骤3.2，因为白噪分布，杂波数据x的似然函数可以表示为：
Figure BDA0003078248520000131
其中，ai表示第i个网格点对应的幅度，
Figure BDA0003078248520000132
为噪声功率，其值未知；I表示单位向量，sk表示字典矩阵D中第k个网格点对应的导向矢量，γk表示对应的字典矩阵D第k个中网格点系数的方差。
因为系数矩阵A未知，可以把系数矩阵A当成多余参数，可以通过对系数矩阵A进行积分得到边缘概率密度函数
Figure BDA0003078248520000133
去估计超参数Γ和
Figure BDA0003078248520000134
使边缘概率密度函数
Figure BDA0003078248520000135
取最大，x的边缘概率密度函数为：
Figure BDA0003078248520000136
直接求解上式比较复杂，因为网格点对应的系数之间相互独立，所以可以利用以下的表达式去估计超参数γi：
Figure BDA0003078248520000137
其中，ai表示A的第i行，进而求得超参数Γ和
Figure BDA0003078248520000138
其中M表达式为：
Figure BDA0003078248520000139
超参数γi的表达式则可以表述成：
Figure BDA00030782485200001310
子步骤3.3，由于上式求解困难，因此引入贝叶斯定理可以得到在γi固定的条件下，系数ai的后验概率密度函数为
Figure BDA0003078248520000141
其中εi表示ai的方差，μi表示ai的均值，两者的表达式为：
Figure BDA0003078248520000142
Figure BDA0003078248520000143
可以得到ai的估计值：
Figure BDA0003078248520000144
系数矩阵Α的估计值为：
Figure BDA0003078248520000145
其中，Γ是对角元素为超参数γi(i＝1,…,NsNd)的对角矩阵，显然，要得到系数向量Α的估计值首先要估计超参数向量γ。
子步骤3.4，采用期望最大化算法求得超参数
Figure BDA0003078248520000146
和
Figure BDA0003078248520000147
其中t表示迭代的次数，该算法分为求期望(E-step)以及最大化期望(M-step)两步，在E-step中，将未知系数ai作为隐变量，在已知超参数
Figure BDA0003078248520000148
和
Figure BDA0003078248520000149
的条件下，计算
Figure BDA00030782485200001410
的后验概率密度函数
Figure BDA00030782485200001411
其中
Figure BDA00030782485200001412
是γi第t次迭代估计。M-step是在已知
Figure BDA00030782485200001413
的后验概率密度函数
Figure BDA00030782485200001414
时求出超参数
Figure BDA00030782485200001415
和
Figure BDA00030782485200001416
当迭代次数超过T次或||γ(t+1)-γ(t)||2/||γ(t+1)||2<η时终止迭代，其中η表示阈值，一般取0.0001。
步骤4，根据所述系数向量估计值得到杂波加噪声协方差矩阵
Figure BDA00030782485200001418
计算空时自适应滤波处理在稀疏贝叶斯迭代方法下的权向量以及在正侧视阵下雷达杂波抑制改善因子，以此来衡量正侧视阵雷达杂波谱的稀疏恢复效果。
具体地，杂波噪声协方差矩阵的表达式为：
Figure BDA00030782485200001417
其中I表示单位矩阵，um,n表示
Figure BDA0003078248520000151
的第m行第n列的元素，1≤m≤K，1≤n≤L，L表示样本数，K表示挑选字典元素的序列号集合。δ≥1为实常数，为对角加载量，
Figure BDA0003078248520000152
为估计的噪声功率。在上式中，采用对角加载技术来避免变形的波束形状和高旁瓣，而控制了λ加载量(在我们的实验中考虑λ＝5～10dB得到了令人满意的结果)。沿杂波脊快速收敛稀疏贝叶斯迭代方法下的STAP滤波器权向量为：
Figure BDA0003078248520000153
其中，s(fdt,fst)表示目标的空时导向矢量，fdt表示正侧视阵雷达接收到的杂波中包含的动目标归一化多普勒频率
空时自适应滤波处理后的非正侧视阵雷达杂波抑制改善因子表达式为：
Figure BDA0003078248520000154
其中tr(·)表示矩阵的迹，fst表示正侧视阵雷达接收到的杂波中包含的动目标归一化空间频率。
仿真实验
(1)杂波数据仿真及实验条件
在本实验中，发射平台坐标为(-150km,0,8km)，接收平台坐标为(0,0,9km)。天线阵面结构采用等距线阵，其中，λ表示载波波长，仿真杂波数据脉冲个数为10；
本实验采用Ward杂波模型进行杂波仿真，并添加高斯白噪声，仿真参数如表1所示：
表1仿真参数
脉冲重复频率 2500Hz
载机速度 150m/s
线阵轴向与载机速度夹角 0°
脉冲数 10
载波波长 0.1m
主波束与阵面夹角 90°
阵元个数 10
接收机带宽 2.5MHz
杂噪比 60dB
载机高度 9000m
阵元间距 0.1m
(2)仿真内容
仿真1，基于上述(1)中的杂波数据仿真及实验条件，分别采用正交匹配追踪(OMP)算法、CVX算法、采样协方差求逆(SMI)算法和本发明的沿杂波脊快速收敛稀疏贝叶斯(RM-SBL)算法对正侧视阵构型的机载雷达杂波仿真数据进行处理，通过各方法的Capon谱估计对比比较处理前后杂波恢复情况，结果如图2所示。其中，图2(a)是用OMP算法得出的Capon谱，图2(b)是用CVX算法得出的Capon谱，图2(c)是用SMI算法得出的Capon谱，图2(d)是用RM-SBL算法得出的Capon谱。
对比图2可以看出，杂波稀疏谱在空时二维平面上表现为几个不连续的杂波散射，由于稀疏恢复方法不尝试恢复杂波块的位置和幅度，而是使用尽可能少的时空导向矢量来恢复杂波子空间，因此传统的稀疏恢复算法恢复的杂波存在一些由噪声引起的虚假峰和弱峰，而本发明恢复的杂波轮廓非常清晰(峰值沿着杂波脊分布)，利用多个训练样本具有相同稀疏结构的特点，对噪声具有很强的稳健性且具有非常快的收敛速度。
图3为各方法杂波抑制改善因子对比图，可以看出本发明副瓣曲线更平稳且更高，主瓣凹口更窄，性能损失更低。
以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以所述权利要求的保护范围为准。

Claims (6)

1.沿杂波脊快速收敛稀疏贝叶斯的杂波抑制方法，其特征在于，包括以下步骤：
步骤1，通过先验知识和直接数据域算法联合估计杂波脊位置，根据所述杂波脊位置对超完备字典矩阵进行网格点选取重构字典矩阵D；
步骤2，根据所述字典矩阵D得到正侧视阵构型下雷达接收到的杂波数据x；
步骤3，求解所述雷达接收到的杂波数据x的似然函数，通过稀疏贝叶斯迭代方法来先求解系数向量的后验概率密度函数，从而得到均值和方差，再采用期望最大化算法得到超参数向量的估计值，最终通过估计的超参数向量得到系数向量估计值；
步骤4，根据所述系数向量估计值得到杂波加噪声协方差矩阵
Figure FDA0003078248510000013
计算空时自适应滤波处理在稀疏贝叶斯迭代方法下的权向量以及在正侧视阵下雷达杂波抑制改善因子，以此来衡量正侧视阵雷达杂波谱的稀疏恢复效果。
2.根据权利要求1所述的沿杂波脊快速收敛稀疏贝叶斯的杂波抑制方法，其特征在于，步骤1具体为：
子步骤1.1，当先验知识已知时，可以直接求得杂波脊位置，具体如下：
正侧阵时，天线锥角余弦值cosψ和多普勒频率fd间的关系为：
Figure FDA0003078248510000011
其中，v表示载机飞行速度，λ表示载波波长，ψ表示天线锥角，fr表示发射脉冲的重复频率，fd是杂波散射单元对应的多普勒频率；
斜侧阵构型时，天线锥角余弦值cosψ和多普勒频率fd间的关系为：
Figure FDA0003078248510000012
其中，
Figure FDA0003078248510000021
表示俯仰角；
通过所述多普勒频率fd与天线锥角ψ的关系可以从超完备字典矩阵中选取出对应的点，即为杂波脊对应的网格点，结合其在矩阵中相邻的点重新构建字典矩阵D；
子步骤1.2，当先验知识未知时，可以通过直接数据域算法求得杂波脊位置，具体如下：
目标信号矩阵由空域导向矢量和时域导向矢量所构成，目标信号沿行方向的相位差为
Figure FDA0003078248510000022
ψS0是目标的空间锥角；目标信号沿列方向的相位差为
Figure FDA0003078248510000023
fd0是目标的多普勒频率；目标信号沿对角线方向的相位差为
Figure FDA0003078248510000024
则空域导向矢量SS(ψS0)、时域导向矢量ST(fd0)及空时导向矢量S为：
Figure FDA0003078248510000025
Figure FDA0003078248510000026
Figure FDA0003078248510000027
其中，(i)T表示转置，j表示复数的虚部，N表示天线平面的阵元数，M表示发射脉冲数；沿空域、时域和空时域做两阵元对消，滤除目标信号并增加样本，分别得到(N-1)×M维空域低维度矩阵XS、N×(M-1)维时域低维度矩阵XT和(N-1)×(M-1)维空时域低维度矩阵XST：
Figure FDA0003078248510000028
Figure FDA0003078248510000029
Figure FDA0003078248510000031
假设空域每次的变化量为Nm，与之对应的时域每次的变化量为Mm，则所述空域低维度矩阵XS经空域滑窗取样，得到(N-Nm)×(M-Mm+1)个样本；所述时域低维度矩阵XT经时域滑窗取样，得到(N-Nm+1)×(M-Mm)个样本；所述空时域低维度矩阵XST经空时域联合滑窗取样，得到(N-Nm)×(M-Mm)个样本；因此，总共可以得到样本数为：
L0＝2((N-Nm)×(M-Mm+1)+(N-Nm+1)×(M-Mm)+(N-Nm)×(M-Mm))
训练样本采用上述L0个样本，用
Figure FDA0003078248510000032
表示，l0＝1,2,…,L0，则协方差矩阵的最大似然估计为：
Figure FDA0003078248510000033
其中，(i)H表示共轭转置，Vec(i)表示对矩阵堆叠成向量的操作；
求得所述协方差矩阵的最大似然估计
Figure FDA0003078248510000034
后，画出Capon谱即可得出杂波脊的位置，通过杂波脊的位置在超完备字典矩阵中对应的网格点结合其相邻网格点重建字典矩阵D。
3.根据权利要求1所述的沿杂波脊快速收敛稀疏贝叶斯的杂波抑制方法，其特征在于，步骤2具体为：
雷达接收到的杂波数据x为：
x＝αi⊙Si+Di-Αi-+N0
其中，x＝[x1 x2 … xL]，L表示正侧视阵雷达接收到的杂波包含的距离门总个数；
Figure FDA0003078248510000035
αi是第i个网格点对应的信号分量的幅度，αip是第p个快拍第i个网格点对应的信号分量的幅度，1≤p≤L，
Figure FDA0003078248510000036
表示Kronecker积；⊙表示Hadamard积；Si＝[si si … si]NM×L，si表示字典矩阵D中第i个网格点对应的导向矢量；
Figure FDA0003078248510000041
表示从字典矩阵D中剔除掉si后的新的字典矩阵，Nd为归一化多普勒频率的细分格数，Ns为归一化空间频率的细分格数；
Figure FDA0003078248510000042
表示从系数矩阵Α中移除的第i个元素αi的新的系数矩阵；N0表示服从零均值复高斯分布的噪声向量。
4.根据权利要求3所述的沿杂波脊快速收敛稀疏贝叶斯的杂波抑制方法，其特征在于，步骤3具体包含以下子步骤：
子步骤3.1，根据系数矩阵A和对角矩阵Γ，可以推导出系数矩阵A的先验概率密度函数：
Figure FDA0003078248510000043
其中，Γ＝diag(γ)为对角矩阵，
Figure FDA0003078248510000044
对应字典矩阵D中网格点系数的方差；
子步骤3.2，所述雷达接收到的杂波数据x的似然函数为：
Figure FDA0003078248510000045
其中，
Figure FDA0003078248510000046
Figure FDA0003078248510000047
为噪声功率，I表示单位向量，sk表示字典矩阵D中第k个网格点对应的导向矢量，γk表示对应字典矩阵D中第k个网格点系数的方差；
对系数矩阵A进行积分得到边缘概率密度函数
Figure FDA0003078248510000048
去估计超参数Γ和
Figure FDA0003078248510000049
使所述边缘概率密度函数
Figure FDA00030782485100000410
取最大，所述边缘概率密度函数
Figure FDA00030782485100000411
为：
Figure FDA00030782485100000412
直接求解上式比较复杂，因为网格点对应的系数之间相互独立，所以可以利用以下的表达式去估计超参数γi：
Figure FDA0003078248510000051
进而求得超参数Γ和
Figure FDA0003078248510000052
其中M表达式为：
Figure FDA0003078248510000053
则超参数γi的表达式为：
Figure FDA0003078248510000054
子步骤3.3，由于上式求解困难，因此引入贝叶斯定理可以得到在γi固定的条件下，系数ai的后验概率密度函数为：
Figure FDA0003078248510000055
其中，εi表示ai的方差，μi表示ai的均值，两者的表达式为：
Figure FDA0003078248510000056
Figure FDA0003078248510000057
可以得到系数ai的估计值：
Figure FDA0003078248510000058
则系数矩阵Α的估计值为：
Figure FDA0003078248510000059
子步骤3.4，采用期望最大化算法求得超参数
Figure FDA00030782485100000510
和
Figure FDA00030782485100000511
其中t表示迭代的次数；
Figure FDA00030782485100000512
表示γi第t+1次迭代估计值，
Figure FDA00030782485100000513
表示
Figure FDA00030782485100000514
第t+1次迭代估计值。
5.根据权利要求4所述的沿杂波脊快速收敛稀疏贝叶斯的杂波抑制方法，其特征在于，子步骤3.4中，所述期望最大化算法分为求期望和最大化期望两步，在求期望中，将未知系数ai作为隐变量，在已知超参数
Figure FDA0003078248510000061
和
Figure FDA0003078248510000062
的条件下，计算
Figure FDA0003078248510000063
的后验概率密度函数
Figure FDA0003078248510000064
在最大化期望中，已知
Figure FDA0003078248510000065
的后验概率密度函数
Figure FDA0003078248510000066
时求出超参数
Figure FDA0003078248510000067
和
Figure FDA0003078248510000068
当迭代次数超过T次或||γ(t+1)-γ(t)||2/||γ(t+1)||2<η时终止迭代；其中，η表示阈值。
6.根据权利要求3所述的沿杂波脊快速收敛稀疏贝叶斯的杂波抑制方法，其特征在于，步骤4中，所述杂波加噪声协方差矩阵
Figure FDA0003078248510000069
为：
Figure FDA00030782485100000610
其中，I表示单位矩阵，um,n表示
Figure FDA00030782485100000611
的第m行第n列的元素，1≤m≤K，1≤n≤L，K表示挑选字典元素的序列号集合；δ≥1为实常数，
Figure FDA00030782485100000612
为估计的噪声功率；
空时自适应滤波处理在稀疏贝叶斯迭代方法下的权向量为：
Figure FDA00030782485100000613
其中，s(fdt,fst)表示目标的空时导向矢量，fdt表示正侧视阵雷达接收到的杂波中包含的动目标归一化多普勒频率；fst表示正侧视阵雷达接收到的杂波中包含的动目标归一化空间频率；
空时自适应滤波处理后的非正侧视阵雷达杂波抑制改善因子为：
Figure FDA00030782485100000614
其中，tr(i)表示矩阵的迹。
CN202110559117.3A 2021-05-21 2021-05-21 沿杂波脊快速收敛稀疏贝叶斯的杂波抑制方法 Active CN113376606B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110559117.3A CN113376606B (zh) 2021-05-21 2021-05-21 沿杂波脊快速收敛稀疏贝叶斯的杂波抑制方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110559117.3A CN113376606B (zh) 2021-05-21 2021-05-21 沿杂波脊快速收敛稀疏贝叶斯的杂波抑制方法

Publications (2)

Publication Number Publication Date
CN113376606A true CN113376606A (zh) 2021-09-10
CN113376606B CN113376606B (zh) 2023-05-26

Family

ID=77571716

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110559117.3A Active CN113376606B (zh) 2021-05-21 2021-05-21 沿杂波脊快速收敛稀疏贝叶斯的杂波抑制方法

Country Status (1)

Country Link
CN (1) CN113376606B (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113777563A (zh) * 2021-09-13 2021-12-10 内蒙古工业大学 基于稀疏表示的杂波抑制方法、装置及存储介质
CN116577734A (zh) * 2023-07-13 2023-08-11 中国人民解放军空军预警学院 基于先验知识的机载雷达精细化杂波仿真方法与装置

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20110187584A1 (en) * 2010-01-29 2011-08-04 Man-On Pun Method for Suppressing Clutter in Space-Time Adaptive Processing Systems
CN104237883A (zh) * 2014-09-15 2014-12-24 西安电子科技大学 一种采用稀疏表示的机载雷达空时自适应处理方法
US20180031690A1 (en) * 2014-06-09 2018-02-01 Src, Inc. Multiplatform GMTI Radar With Adaptive Clutter Suppression
CN108387884A (zh) * 2018-05-25 2018-08-10 西安电子科技大学 基于知识辅助稀疏渐进最小方差的机载雷达杂波抑制方法
CN112415475A (zh) * 2020-11-13 2021-02-26 中国民航大学 一种基于原子范数的无网格稀疏恢复非正侧视阵stap方法
CN112415476A (zh) * 2020-11-13 2021-02-26 中国民航大学 一种基于稀疏贝叶斯学习的字典失配杂波空时谱估计方法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20110187584A1 (en) * 2010-01-29 2011-08-04 Man-On Pun Method for Suppressing Clutter in Space-Time Adaptive Processing Systems
US20180031690A1 (en) * 2014-06-09 2018-02-01 Src, Inc. Multiplatform GMTI Radar With Adaptive Clutter Suppression
CN104237883A (zh) * 2014-09-15 2014-12-24 西安电子科技大学 一种采用稀疏表示的机载雷达空时自适应处理方法
CN108387884A (zh) * 2018-05-25 2018-08-10 西安电子科技大学 基于知识辅助稀疏渐进最小方差的机载雷达杂波抑制方法
CN112415475A (zh) * 2020-11-13 2021-02-26 中国民航大学 一种基于原子范数的无网格稀疏恢复非正侧视阵stap方法
CN112415476A (zh) * 2020-11-13 2021-02-26 中国民航大学 一种基于稀疏贝叶斯学习的字典失配杂波空时谱估计方法

Non-Patent Citations (5)

* Cited by examiner, † Cited by third party
Title
FUYU TAO ET AL.: "A knowledge aided SPICE space time adaptive processing method for airborne radar with conformal array", 《SIGNAL PROCESSING》 *
HONGDA YE ET AL.: "Clutter-Ridge Matched SR-STAP Technique for Non-stationary Clutter Suppression", 《2020 IEEE RADAR CONFERENCE (RADARCONF20)》 *
同亚龙等: "基于杂波脊先验信息的非均匀杂波抑制方法", 《系统工程与电子技术》 *
段克清等: "稀疏恢复空时自适应处理技术研究综述", 《电子学报》 *
董瑞军等: "机载雷达的预滤波直接数据域算法", 《电子与信息学报》 *

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113777563A (zh) * 2021-09-13 2021-12-10 内蒙古工业大学 基于稀疏表示的杂波抑制方法、装置及存储介质
CN116577734A (zh) * 2023-07-13 2023-08-11 中国人民解放军空军预警学院 基于先验知识的机载雷达精细化杂波仿真方法与装置
CN116577734B (zh) * 2023-07-13 2023-09-22 中国人民解放军空军预警学院 基于先验知识的机载雷达精细化杂波仿真方法与装置

Also Published As

Publication number Publication date
CN113376606B (zh) 2023-05-26

Similar Documents

Publication Publication Date Title
CN110275166B (zh) 基于admm的快速稀疏孔径isar自聚焦与成像方法
CN109188344B (zh) 脉冲噪声环境下基于互循环相关music算法信源个数与来波方向角估计方法
CN104237883B (zh) 一种采用稀疏表示的机载雷达空时自适应处理方法
CN109212500A (zh) 一种基于稀疏重构的ka-stap杂噪协方差矩阵高精度估计方法
CN113376606A (zh) 沿杂波脊快速收敛稀疏贝叶斯的杂波抑制方法
CN112612006B (zh) 基于深度学习的机载雷达非均匀杂波抑制方法
CN107015214B (zh) 一种基于稀疏贝叶斯学习的空时自适应处理方法
CN105911527B (zh) 基于efa与mwf的机载雷达空时自适应处理方法
CN113933808A (zh) 机载雷达动目标检测方法、装置、设备及存储介质
CN106970358B (zh) 非正侧视阵雷达杂波谱的角度多普勒配准的优化方法
CN116699526A (zh) 一种基于稀疏与低秩模型的车载毫米波雷达干扰抑制方法
Halimi et al. Cramér-Rao bounds and estimation algorithms for delay/Doppler and conventional altimetry
CN107748364A (zh) 基于降秩多级维纳滤波器的低空风切变风场速度估计方法
CN108896963B (zh) 机载雷达空时自适应降维处理方法
CN111474527B (zh) 机载stap雷达快速去互耦的杂波协方差矩阵估计方法
CN113376607A (zh) 机载分布式雷达小样本空时自适应处理方法
CN112904298A (zh) 一种基于局部网格分裂的网格偏离空时自适应处理方法
CN111044996A (zh) 一种基于降维近似消息传递的lfmcw雷达目标检测方法
CN104035078A (zh) 一种基于阵元阶数递推的降维空时自适应权值计算方法
CN110850421A (zh) 基于混响对称谱的空时自适应处理的水下目标检测方法
CN107167782B (zh) 基于信杂噪比最大的雷达三维异构阵稀疏重构方法
CN115877380A (zh) 一种sar多运动目标成像方法、装置和存储介质
CN115421115A (zh) 一种用于联合相位校正与isar成像的重赋权交替方向乘子法
CN110609255B (zh) 一种基于特征波束的自适应波束域fsa的杂波抑制降维方法
Chen et al. A sea clutter suppression algorithm for over-the-horizon radar based on dictionary learning and subspace estimation

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant