CN113219416B - 一种基于雷达缺失采样的微动信号分离方法 - Google Patents

一种基于雷达缺失采样的微动信号分离方法 Download PDF

Info

Publication number
CN113219416B
CN113219416B CN202110475131.5A CN202110475131A CN113219416B CN 113219416 B CN113219416 B CN 113219416B CN 202110475131 A CN202110475131 A CN 202110475131A CN 113219416 B CN113219416 B CN 113219416B
Authority
CN
China
Prior art keywords
frequency
micro
time
signal
sampling
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.)
Active
Application number
CN202110475131.5A
Other languages
English (en)
Other versions
CN113219416A (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.)
Air Force Engineering University of PLA
Original Assignee
Air Force Engineering University of PLA
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 Air Force Engineering University of PLA filed Critical Air Force Engineering University of PLA
Priority to CN202110475131.5A priority Critical patent/CN113219416B/zh
Publication of CN113219416A publication Critical patent/CN113219416A/zh
Application granted granted Critical
Publication of CN113219416B publication Critical patent/CN113219416B/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/28—Details of pulse systems
    • G01S7/285—Receivers
    • G01S7/292—Extracting wanted echo-signals
    • G01S7/2923—Extracting wanted echo-signals based on data belonging to a number of consecutive radar periods
    • 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
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/02—Systems using reflection of radio waves, e.g. primary radar systems; Analogous systems
    • G01S13/50—Systems of measurement based on relative movement of target
    • G01S13/58—Velocity or trajectory determination systems; Sense-of-movement determination systems
    • G01S13/583—Velocity or trajectory determination systems; Sense-of-movement determination systems using transmission of continuous unmodulated waves, amplitude-, frequency-, or phase-modulated waves and based upon the Doppler effect resulting from movement of targets
    • 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
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/66—Radar-tracking systems; Analogous systems
    • G01S13/72—Radar-tracking systems; Analogous systems for two-dimensional tracking, e.g. combination of angle and range tracking, track-while-scan radar
    • G01S13/723—Radar-tracking systems; Analogous systems for two-dimensional tracking, e.g. combination of angle and range tracking, track-while-scan radar by using numerical data
    • 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/35—Details of non-pulse systems
    • G01S7/352—Receivers
    • G01S7/354—Extracting wanted echo-signals

Landscapes

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

Abstract

本发明涉及一种基于雷达缺失采样的微动信号分离方法，本发明首先建立缺损采样的窄带雷达微动信号回波模型，其次将信号模型转换成稀疏正则化模型并通过软阈值迭代算法(ISTA)求解得到重构的时频分布，然后提取重构时频谱的极大值点，并将极大值点关联问题定义为指派问题，通过卡尔曼滤波和匈牙利算法求解该问题得到瞬时频率轨迹，最后，采用本征线性调频分量分解(ICCD)算法分离弹道目标上每个散射点的微动回波信号。相对于传统的时频滤波组分解方法，能够较好地分解在时频域交叉重叠的信号分量，且更适用于雷达采样缺损条件下的微动信号分解。

Description

一种基于雷达缺失采样的微动信号分离方法
技术领域
本发明属于雷达信号处理领域，具体涉及一种基于雷达缺失采样的微动信号分离方法，用于窄带雷达复杂电磁环境下的弹道目标微动特征提取。
背景技术
在雷达探测弹道目标场景下，弹道导弹在飞行过程中除了平动飞行以外，同时会绕自身对称轴做自旋运动并有可能绕某一定向轴作锥旋运动，这些精细的运动会对雷达回波产生除平动多普勒以外的附加频率调制，这一现象被称作微多普勒效应，为目标识别和雷达成像的研究提供了有效信息。在过去的十年中，窄带雷达弹道目标微动特征提取已得到了深入的研究。
在实际应用中，雷达回波采样的缺损和微动信号在时频域交叠已经成为阻碍微动特征有效提取的两个突出问题。在回波受到强干扰期间，采样的数据不得不被丢弃以免影响后续的数据处理，当雷达出现故障和异常值时也会影响回波信号的有效采集。文献1(X.Bai,F.Zhou,and Y.Hui,Obtaining JTF-signature of space-debris fromincomplete and phase-corrupted Data[J].IEEE Trans.Aerosp.Electron.Syst.2017,53(3):1169-1180,June 2017.doi:10.1109/TAES.2017.2667899.)利用一种非参数联合时频字典从损坏的微动信号中重构出聚焦良好的时频分布，但是没有提供针对具有多散射点微动目标信号的提取方法。多散射点微动目标的信号提取实质上可以当作多分量微动信号分解问题处理，一般情况下，基于多分量微动信号分解的参数化方法都会提前设定微动的瞬时多普勒频率多项式曲线或者正弦曲线模型，但是弹道目标微动信号有可能呈现出更为复杂的调制特性(例如滑动散射点的回波信号)，需要设置多个参数来拟合调制特性，从而增加了算法的复杂度。而大多数非参数化的方法，例如变分模态分解、结合短时傅里叶变换的同步挤压变换、稀疏短时傅里叶变换，都无法处理微动信号在时频域出现交叠的情况，然而弹道目标的多分量微动信号在时频域常常出现交叠的情况。文献2(Y.Wang,X.Huang andR.Cao,Novel Approach for ISAR Cross-Range Scaling Based on the MultidelayDiscrete Polynomial-Phase Transform Combined With Keystone Transform[J].IEEETrans.Geosci.Remote Sens.,2020,58(2):1221-1231.doi:10.1109/TGRS.2019.2944674.)提出一种岭迹重排的方法提取时频域的瞬时频率曲线，并结合本征线性调频分量分解(ICCD)算法有效分解了在时频域分解重叠的微动信号，但是在时频分布因采样缺损出现散焦和缺损的情况下，岭迹重排算法无法有效运行；短时变分模态分解(STVMD)算法可以在瞬时频率不规则调制和低信噪比情况下正确分解具有交叠状态的微动信号，但是当微动信号在时频域完全重叠时，分解效果不佳。
发明内容
要解决的技术问题
针对上述现有技术在雷达采样缺失场景下无法有效的实现微动信号分解的问题，本发明提出了一种基于雷达缺失采样的微动信号分离方法。
技术方案
一种基于雷达缺失采样的微动信号分离方法，其特征在于步骤如下：
步骤1：建立缺损采样的窄带雷达微动信号回波模型：
s＝BGf+n
其中G是HQ的伪逆矩阵，H是傅里叶变换矩阵，Q是滑动的窗函数，B是表示信号缺损的对角矩阵，n＝[η(1),...η(n),...η(N)]T是加性噪声；f是时频表征；
步骤2：将信号回波模型转换成稀疏正则化模型：
Figure BDA0003047123660000021
其中||·||1表示1范数运算符；
步骤3：通过软阈值迭代算法ISTA求解稀疏正则化模型得到重构的时频分布：
ISTA是一种近端梯度下降法，邻近算子为
Figure BDA0003047123660000031
其中
Figure BDA0003047123660000032
是函数
Figure BDA0003047123660000033
处的梯度，μ>0是步长，邻近算子可以被软阈值函数代替：
Figure BDA0003047123660000034
其中
Figure BDA0003047123660000035
λ1＝0.01，soft(z,λ1)是软阈值函数：
soft(z,λ1)＝sign(z)max{0,|z|-λ1}
通过软阈值迭代法，得到了缺失采样条件下的稀疏、去噪且锐化的时频谱分布；
步骤4：提取重构时频谱的极大值点，将重构的时频谱分布f取模并重排为由N个时刻频谱组成的时频图TF(t,f)；TF(t,f)中极大值点的位置坐标为qt,g(t)＝(t,fg(t))，其中t＝t0,...,tN-1，fg(t)代表第t时刻的第g(t)个极值点在频谱上的坐标，t时刻的频谱总共有G(t)个极值点；所述的极大值点的坐标qt,g(t)由下式得到
Figure BDA0003047123660000036
其中α代表阈值参数；为了避免交叉点附近极值点对后续算法的影响，采用下式消除交叉点的极值点
|fg(t)-fg_other(t)|≥df
其中fg_other(t)代表第t个频谱除了g(t)以外的其他极大值点在频谱上的位置，将不满足上式的极大值点剔除，得到qt,g(t)；
步骤5：通过卡尔曼滤波和匈牙利算法相结合的多目标跟踪算法将时频平面的极大值点关联成每个散射点的瞬时频率轨迹；
选取G(tm)＝P时刻的极大值点
Figure BDA0003047123660000037
作为卡尔曼滤波的初始测量值，第k个散射点瞬时频率轨迹在tm+1时刻的预测值为
Figure BDA0003047123660000041
其中
Figure BDA0003047123660000042
代表速度，dt＝1/fs，fs代表采样频率；则第k个散射点在tm+1时刻的平滑值为
Figure BDA0003047123660000043
其中
Figure BDA0003047123660000044
和
Figure BDA0003047123660000045
分别代表过程噪声和观测噪声；测量值
Figure BDA0003047123660000046
是通过选取tm+1时刻的极值点
Figure BDA0003047123660000047
得到的；本质上将G(tm+1)个测量值分配到P个瞬时频率轨迹可以被构造为一个指派问题；如果G(tm+1)＝P，该指派问题是平衡指派问题，数学模型可表示为
Figure BDA0003047123660000048
Figure BDA0003047123660000049
其中ygk＝0代表第g(tm+1)个测量值与第k测量值不匹配，ygk＝1代表g(tm+1)个测量值与第k测量值匹配；指派问题的核心问题是如何构造效率矩阵，其效率矩阵Cm+1可以通过高斯核密度函数构造：
Figure BDA00030471236600000410
其中
Figure BDA00030471236600000411
通过匈牙利算法求解平衡指派问题，因为雷达接收的信号是缺损采样且受噪声干扰，所以会出现G(tm)≠P的情况；当G(tm+1)<P时数学模型所描述的问题是一个非平衡指派问题，可以通过“虚设变量”的方法变为平衡指派问题，虚设P-G(tm+1)个测量值，并对效率矩阵Cm+1添加P-G(tm+1)行零元素，转换得到的平衡指派问题依然可以通过匈牙利算法求解，而与虚设测量值匹配的瞬时频率轨迹在该时刻的测量值则采用卡尔曼滤波的预测值，即
Figure BDA0003047123660000051
其中b代表与虚设测量值匹配的瞬时频率轨迹ID号；当G(tm+1)＝0或者G(tm+1)>P时，P个瞬时频率轨迹在该时刻的测量值采用卡尔曼滤波的预测值，即
Figure BDA0003047123660000052
第k个瞬时频率轨迹可以由平滑值表示，即fk(t)＝zt,k；
步骤6：采用本征线性调频分量分解ICCD算法分离弹道目标上每个散射点的微动回波信号。
优选地：步骤4中df＝N/20，N为信号采样的点数机离散信号的长度。
优选地：步骤5中
Figure BDA0003047123660000053
N为信号采样的点数机离散信号的长度。
优选地：步骤6具体如下：
为了将ICCD方法转换成矩阵的形式，将微动目标窄带雷达回波信号改写为以下离散的形式
Figure BDA0003047123660000054
其中
Figure BDA0003047123660000055
t＝t0,…,tN-1；因为
Figure BDA0003047123660000056
是一个平滑的函数，
Figure BDA0003047123660000057
可以用傅里叶基表示为
Figure BDA0003047123660000058
其中f0＝Fs/QN,Fs是采样频率，Q是一个正整数，时频滤波器组的带宽为Lf0，
Figure BDA0003047123660000059
是傅里叶系数；根据傅里叶基的表示，将离散回波信号转换为矩阵形式
sN×1＝FN×K(2L+1)AK(2L+1)×1+nN×1
其中A＝[A1...Ak...AK]T，矩阵A由K个列向量Ak组成，
Figure BDA00030471236600000510
F＝[F1...Fk...FK]，Fk＝CkT，其中Ck是一个N×N的对角矩阵；
Figure BDA0003047123660000061
作为一个N×K(2L+1)的矩阵，元素组成如下所示
Figure BDA0003047123660000062
为了得到未知的参数A，描述成一个优化问题：
Figure BDA0003047123660000063
因为N>K(2L+1)，其优化模型可以表示为：
Figure BDA0003047123660000064
根据吉洪诺夫正则化方法，上式的解为
Figure BDA0003047123660000065
其中I是一个尺寸为K(2L+1)的单位矩阵，得到每个散射点的信号：
Figure BDA0003047123660000066
优选地：Q＝1、λ2＝5。
一种计算机系统，其特征在于包括：一个或多个处理器，计算机可读存储介质，用于存储一个或多个程序，其中，当所述一个或多个程序被所述一个或多个处理器执行时，使得所述一个或多个处理器实现上述的方法。
一种计算机可读存储介质，其特征在于存储有计算机可执行指令，所述指令在被执行时用于实现上述的方法。
一种计算机程序，其特征在于包括计算机可执行指令，所述指令在被执行时用于实现上述的方法。
有益效果
本发明提出的一种基于雷达缺失采样的微动信号分离方法，采用软阈值迭代重构缺损采样条件下的时频分布并采用本征线性调频分量算法结合提取的瞬时频率轨迹分解出每个散射点的微动信号。将多目标跟踪算法即匈牙利算法和卡尔曼滤波算法应用于软阈值迭代算法重构的时频图中提取各散射点的瞬时频率轨迹，当重构的时频图因雷达采样缺失出现缺损时或瞬时频率轨迹交叠时采用卡尔曼滤波的预测值代替测量值，提高了该方法的鲁棒性。相对于传统的时频滤波组分解方法，能够较好地分解在时频域交叉重叠的信号分量，解决了雷达采样缺损条件下的微动信号分解的问题，填补了现有技术中这一技术空白。
附图说明
附图仅用于示出具体实施例的目的，而并不认为是对本发明的限制，在整个附图中，相同的参考符号表示相同的部件。
图1采样缺损微动信号的时频分布
图2重构的时频分布
图3时频分布的极大值
图4提取的各散射点瞬时频率轨迹
图5分解的散射点1信号
图6分解的散射点2信号
图7分解的散射点3信号
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图和实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。此外，下面描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。
本发明提出的基于雷达缺失采样的微动信号分离方法通过如下方式实现的：
首先建立缺损采样的窄带雷达微动信号回波模型；
其次将信号模型转换成稀疏正则化模型并通过软阈值迭代算法(ISTA)求解得到重构的时频分布；
然后提取重构时频谱的极大值点，并将极大值点关联问题定义为指派问题，通过卡尔曼滤波和匈牙利算法求解该问题得到瞬时频率轨迹；
最后，采用本征线性调频分量分解(ICCD)算法分离弹道目标上每个散射点的微动回波信号。
具体包括下列步骤：
步骤一：信号模型的构建，在高频电磁条件下，金属目标可以由点散射模型表示，回波信号缺损条件下的微动目标窄带雷达回波信号表示为
Figure BDA0003047123660000081
其中P表示散射点的数目，
Figure BDA0003047123660000082
Figure BDA0003047123660000083
表示回波信号的缺损情况，ak(t)>0代表瞬时幅度，
Figure BDA0003047123660000084
代表初始相位，fk(t)代表瞬时频率，瞬时幅度和瞬时频率被认为是平稳的比其相位函数变化慢得多的函数，即满足|a′k(t)|,|f′k(t)|＜＜|fk(t)|，η(t)代表加性噪声。
s(t)是一个典型的非平稳信号。时频分析是一种在时频域表征非平稳信号能量的有效方法，非平稳信号可以在时频域上表征出频域中无法表示的特征。短时傅里叶变换具备线性和没有交叉项等诸多优势，因此本发明采用短时傅里叶变换分析回波信号。
短时傅里叶变换可以表示为
TF0(t,f)＝∫s(τ)gσ(τ-t)exp(-j2πfτ)dτ      (2)
其中TF0(t,f)代表时频表征，f代表瞬时频率，gσ代表长度为σ的窗函数，用来保证加窗信号足够平稳，窗函数有多种形式，本发明使用高斯窗σ＝N/5；它可以很容易地用其他形式重写，比如矩形窗和汉明窗等。为了方便描述重构算法，结合式(1)，将式(2)转换为矩阵的形式
Figure BDA0003047123660000091
其中H＝diag{W,...,W}，WN×N是傅里叶变换矩阵，N＝Mms，其中ms＝1是滑动窗口的步数，Q＝[P1...Pm...PM]T是滑动的窗函数，其中Pm＝diag{pm(1),...pm(n),...pm(N)},n∈[1,N]是对角矩阵。B＝diag{b(1),...b(n),...b(N)},b(n)∈{0,1}是表示信号缺损的对角矩阵,s0＝[s0(1),...s0(n),...s0(N)]T代表没有缺损的理想回波信号，(·)T代表转置，根据式(3)，采样缺损时的信号如下所示
s＝BGf+n           (4)
其中G是HQ的伪逆矩阵，n＝[η(1),...η(n),...η(N)]T是加性噪声。
步骤二：为了从幅度起伏和采样缺损的回波s中得到理想的时频表征f，令s＝[s(1),...s(n),...s(N)]T，式(4)可以被描述成一个优化问题
Figure BDA0003047123660000092
其中||·||2表示2范数运算符。因为强散射点在目标上是稀疏分布的，所以
Figure BDA0003047123660000093
可以利用稀疏先验的信息被重构，式(5)表示的优化模型可以被表示L1正则化模型
Figure BDA0003047123660000094
其中||·||1表示1范数运算符。L1正则化模型可以被认为是一个套索回归问题，可以采用ISTA求解。ISTA是一种近端梯度下降法，邻近算子为
Figure BDA0003047123660000095
其中
Figure BDA0003047123660000096
是函数
Figure BDA0003047123660000097
处的梯度，μ>0是步长，邻近算子可以被软阈值函数代替
Figure BDA0003047123660000098
其中
Figure BDA0003047123660000101
λ1＝0.01，soft(z,λ1)是软阈值函数
soft(z,λ1)＝sign(z)max{0,|z|-λ1}      (9)
由图1和图2所示，通过软阈值迭代法，得到了缺失采样条件下的稀疏、去噪且锐化的时频谱分布。
步骤三：提取重构时频谱的极大值点，将重构的时频谱分布f取模并重排为由N个时刻频谱组成的时频图TF(t,f)。TF(t,f)中极大值点的位置坐标为qt,g(t)＝(t,fg(t))，其中t＝t0,...,tN-1，fg(t)代表第t时刻的第g(t)个极值点在频谱上的坐标，t时刻的频谱总共有G(t)个极值点。极大值点的坐标qt,g(t)由式(10)得到
Figure BDA0003047123660000102
其中α代表阈值参数。为了避免交叉点附近极值点对后续算法的影响，采用下式消除交叉点的极值点
|fg(t)-fg_other(t)|≥df        (11)
其中fg_other(t)代表第t个频谱除了g(t)以外的其他极大值点在频谱上的位置，将不满足式(11)的极大值点剔除，得到qt,g(t)，其在时频分布极大值的位置由图2所示，在本发明中df＝N/20。
步骤四：引入卡尔曼滤波和匈牙利算法相结合的多目标跟踪算法将时频平面的极大值点关联成每个散射点的瞬时频率轨迹。选取G(tm)＝P时刻的极大值点
Figure BDA0003047123660000103
作为卡尔曼滤波的初始测量值，瞬时频率的位置和速度由线性状态空间描述。为了简化分析，选用匀速(CV)模型，第k个散射点瞬时频率轨迹在tm+1时刻的预测值为
Figure BDA0003047123660000104
其中
Figure BDA0003047123660000105
代表速度，dt＝1/fs，fs代表采样频率。则第k个散射点在tm+1时刻的平滑值为
Figure BDA0003047123660000111
其中
Figure BDA0003047123660000112
和
Figure BDA0003047123660000113
分别代表过程噪声和观测噪声。测量值
Figure BDA0003047123660000114
是通过选取tm+1时刻的极值点
Figure BDA0003047123660000115
得到的。本质上将G(tm+1)个测量值分配到P个瞬时频率轨迹可以被构造为一个指派问题。如果G(tm+1)＝P，该指派问题是平衡指派问题，数学模型可表示为
Figure BDA0003047123660000116
Figure BDA0003047123660000117
其中ygk＝0代表第g(tm+1)个测量值与第k测量值不匹配，ygk＝1代表g(tm+1)个测量值与第k测量值匹配。指派问题的核心问题是如何构造效率矩阵，如式(14)所示的指派问题，其效率矩阵Cm+1可以通过高斯核密度函数构造。
Figure BDA0003047123660000118
其中
Figure BDA0003047123660000119
在本发明中
Figure BDA00030471236600001110
平衡指派问题可以通过匈牙利算法求解。因为雷达接收的信号是缺损采样且受噪声干扰，所以会出现G(tm)≠P的情况；当G(tm+1)<P时式(14)所描述的问题是一个非平衡指派问题，可以通过“虚设变量”的方法变为平衡指派问题，虚设P-G(tm+1)个测量值，并对效率矩阵Cm+1添加P-G(tm+1)行零元素，转换得到的平衡指派问题依然可以通过匈牙利算法求解，而与虚设测量值匹配的瞬时频率轨迹在该时刻的测量值则采用卡尔曼滤波的预测值，即
Figure BDA00030471236600001111
其中b代表与虚设测量值匹配的瞬时频率轨迹ID号；当G(tm+1)＝0或者G(tm+1)>P时，P个瞬时频率轨迹在该时刻的测量值采用卡尔曼滤波的预测值，即
Figure BDA00030471236600001112
第k个瞬时频率轨迹可以由平滑值表示，即fk(t)＝zt,k。得到的各散射点的瞬时频率轨迹由图4所示。
步骤五：采用ICCD算法结合提取的各散射点瞬时频率轨迹分解微动信号，即将步骤三提取的各散射点瞬时频率轨迹fk(t)输入到ICCD算法中以分解微动信号。ICCD算法可以被当作一种时频域的滤波器组，每个滤波器组的中心频率由fk(t)表示，滤波器组的带宽需要使用者预先设定。为了将ICCD方法转换成矩阵的形式，将式(1)改写为以下离散的形式
Figure BDA0003047123660000121
其中
Figure BDA0003047123660000122
t＝t0,…,tN-1。因为
Figure BDA0003047123660000123
是一个平滑的函数，
Figure BDA0003047123660000124
可以用傅里叶基表示为
Figure BDA0003047123660000125
其中f0＝Fs/QN,Fs是采样频率，Q是一个正整数，本发明中取Q＝1，时频滤波器组的带宽为Lf0，
Figure BDA0003047123660000126
是傅里叶系数。将式(12)代入到式(11)，可以得到式(11)的矩阵形式
sN×1＝FN×K(2L+1)AK(2L+1)×1+nN×1         (18)
其中A＝[A1...Ak...AK]T，矩阵A由K个列向量Ak组成，
Figure BDA0003047123660000127
F＝[F1...Fk...FK]，Fk＝CkT，其中Ck是一个N×N的对角矩阵
Figure BDA0003047123660000128
作为一个N×K(2L+1)的矩阵，元素组成如下所示
Figure BDA0003047123660000129
为了得到未知的参数A，式(13)可以被描述成一个优化问题
Figure BDA00030471236600001210
因为N>K(2L+1)，所以式(16)所描述的逆问题一般是一个不适定问题，其优化模型可以表示为
Figure BDA0003047123660000131
根据吉洪诺夫正则化方法，上式的解为
Figure BDA0003047123660000132
其中I是一个尺寸为K(2L+1)的单位矩阵，在本发明中λ2＝5。根据式(23)所得到的
Figure BDA0003047123660000133
每个散射点的信号可以由如下公式得到
Figure BDA0003047123660000134
至此完成了信号缺损采样条件下各散射点微动信号的分解。由图5、图6和图7所示。
以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明公开的技术范围内，可轻易想到各种等效的修改或替换，这些修改或替换都应涵盖在本发明的保护范围之内。

Claims (8)

1.一种基于雷达缺失采样的微动信号分离方法，其特征在于步骤如下：
步骤1：建立缺损采样的窄带雷达微动信号回波模型：
s＝BGf+n
其中G是HQ的伪逆矩阵，H是傅里叶变换矩阵，Q是滑动的窗函数，B是表示信号缺损的对角矩阵，n＝[η(1),...η(n),...η(N)]T是加性噪声；f是时频表征；
步骤2：将信号回波模型转换成稀疏正则化模型：
Figure FDA0003047123650000011
其中||·||1表示1范数运算符；
步骤3：通过软阈值迭代算法ISTA求解稀疏正则化模型得到重构的时频分布：
ISTA是一种近端梯度下降法，邻近算子为
Figure FDA0003047123650000012
其中
Figure FDA0003047123650000013
是函数
Figure FDA0003047123650000014
在
Figure FDA0003047123650000015
处的梯度，μ>0是步长，邻近算子可以被软阈值函数代替：
Figure FDA0003047123650000016
其中
Figure FDA0003047123650000017
λ1＝0.01，soft(z,λ1)是软阈值函数：
soft(z,λ1)＝sign(z)max{0,|z|-λ1}
通过软阈值迭代法，得到了缺失采样条件下的稀疏、去噪且锐化的时频谱分布；
步骤4：提取重构时频谱的极大值点，将重构的时频谱分布f取模并重排为由N个时刻频谱组成的时频图TF(t,f)；TF(t,f)中极大值点的位置坐标为qt,g(t)＝(t,fg(t))，其中t＝t0,...,tN-1，fg(t)代表第t时刻的第g(t)个极值点在频谱上的坐标，t时刻的频谱总共有G(t)个极值点；所述的极大值点的坐标qt,g(t)由下式得到
Figure FDA0003047123650000018
其中α代表阈值参数；为了避免交叉点附近极值点对后续算法的影响，采用下式消除交叉点的极值点
|fg(t)-fg_other(t)|≥df
其中fg_other(t)代表第t个频谱除了g(t)以外的其他极大值点在频谱上的位置，将不满足上式的极大值点剔除，得到qt,g(t)；
步骤5：通过卡尔曼滤波和匈牙利算法相结合的多目标跟踪算法将时频平面的极大值点关联成每个散射点的瞬时频率轨迹；
选取G(tm)＝P时刻的极大值点
Figure FDA0003047123650000021
作为卡尔曼滤波的初始测量值，第k个散射点瞬时频率轨迹在tm+1时刻的预测值为
Figure FDA0003047123650000022
其中
Figure FDA0003047123650000023
代表速度，dt＝1/fs，fs代表采样频率；则第k个散射点在tm+1时刻的平滑值为
Figure FDA0003047123650000024
其中
Figure FDA0003047123650000025
和
Figure FDA0003047123650000026
分别代表过程噪声和观测噪声；测量值
Figure FDA0003047123650000027
是通过选取tm+1时刻的极值点
Figure FDA0003047123650000028
得到的；本质上将G(tm+1)个测量值分配到P个瞬时频率轨迹可以被构造为一个指派问题；如果G(tm+1)＝P，该指派问题是平衡指派问题，数学模型可表示为
Figure FDA0003047123650000029
Figure FDA00030471236500000210
其中ygk＝0代表第g(tm+1)个测量值与第k测量值不匹配，ygk＝1代表g(tm+1)个测量值与第k测量值匹配；指派问题的核心问题是如何构造效率矩阵，其效率矩阵Cm+1可以通过高斯核密度函数构造：
Figure FDA0003047123650000031
其中
Figure FDA0003047123650000032
通过匈牙利算法求解平衡指派问题，因为雷达接收的信号是缺损采样且受噪声干扰，所以会出现G(tm)≠P的情况；当G(tm+1)<P时数学模型所描述的问题是一个非平衡指派问题，可以通过“虚设变量”的方法变为平衡指派问题，虚设P-G(tm+1)个测量值，并对效率矩阵Cm+1添加P-G(tm+1)行零元素，转换得到的平衡指派问题依然可以通过匈牙利算法求解，而与虚设测量值匹配的瞬时频率轨迹在该时刻的测量值则采用卡尔曼滤波的预测值，即
Figure FDA0003047123650000033
其中b代表与虚设测量值匹配的瞬时频率轨迹ID号；当G(tm+1)＝0或者G(tm+1)>P时，P个瞬时频率轨迹在该时刻的测量值采用卡尔曼滤波的预测值，即
Figure FDA0003047123650000034
第k个瞬时频率轨迹可以由平滑值表示，即fk(t)＝zt,k；
步骤6：采用本征线性调频分量分解ICCD算法分离弹道目标上每个散射点的微动回波信号。
2.根据权利要求1所述的一种基于雷达缺失采样的微动信号分离方法，其特征在于步骤4中df＝N/20，N为信号采样的点数机离散信号的长度。
3.根据权利要求1所述的一种基于雷达缺失采样的微动信号分离方法，其特征在于步骤5中
Figure FDA0003047123650000035
N为信号采样的点数机离散信号的长度。
4.根据权利要求1所述的一种基于雷达缺失采样的微动信号分离方法，其特征在于步骤6具体如下：
为了将ICCD方法转换成矩阵的形式，将微动目标窄带雷达回波信号改写为以下离散的形式
Figure FDA0003047123650000041
其中
Figure FDA0003047123650000042
因为
Figure FDA0003047123650000043
是一个平滑的函数，
Figure FDA0003047123650000044
可以用傅里叶基表示为
Figure FDA0003047123650000045
其中f0＝Fs/QN,Fs是采样频率，Q是一个正整数，时频滤波器组的带宽为Lf0，
Figure FDA0003047123650000046
是傅里叶系数；根据傅里叶基的表示，将离散回波信号转换为矩阵形式
sN×1＝FN×K(2L+1)AK(2L+1)×1+nN×1
其中A＝[A1...Ak...AK]T，矩阵A由K个列向量Ak组成，
Figure FDA0003047123650000047
F＝[F1...Fk...FK]，Fk＝CkT，其中Ck是一个N×N的对角矩阵
Figure FDA0003047123650000048
作为一个N×K(2L+1)的矩阵，元素组成如下所示
Figure FDA0003047123650000049
为了得到未知的参数A，描述成一个优化问题：
Figure FDA00030471236500000410
因为N>K(2L+1)，其优化模型可以表示为：
Figure FDA00030471236500000411
根据吉洪诺夫正则化方法，上式的解为
Figure FDA00030471236500000412
其中I是一个尺寸为K(2L+1)的单位矩阵，得到每个散射点的信号：
Figure FDA00030471236500000413
Figure FDA0003047123650000051
5.根据权利要求4所述的一种基于雷达缺失采样的微动信号分离方法，其特征在于Q＝1、λ2＝5。
6.一种计算机系统，其特征在于包括：一个或多个处理器，计算机可读存储介质，用于存储一个或多个程序，其中，当所述一个或多个程序被所述一个或多个处理器执行时，使得所述一个或多个处理器实现权利要求1所述的方法。
7.一种计算机可读存储介质，其特征在于存储有计算机可执行指令，所述指令在被执行时用于实现权利要求1所述的方法。
8.一种计算机程序，其特征在于包括计算机可执行指令，所述指令在被执行时用于实现权利要求1所述的方法。
CN202110475131.5A 2021-04-29 2021-04-29 一种基于雷达缺失采样的微动信号分离方法 Active CN113219416B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110475131.5A CN113219416B (zh) 2021-04-29 2021-04-29 一种基于雷达缺失采样的微动信号分离方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110475131.5A CN113219416B (zh) 2021-04-29 2021-04-29 一种基于雷达缺失采样的微动信号分离方法

Publications (2)

Publication Number Publication Date
CN113219416A CN113219416A (zh) 2021-08-06
CN113219416B true CN113219416B (zh) 2023-04-25

Family

ID=77090331

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110475131.5A Active CN113219416B (zh) 2021-04-29 2021-04-29 一种基于雷达缺失采样的微动信号分离方法

Country Status (1)

Country Link
CN (1) CN113219416B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114330420B (zh) * 2021-12-01 2022-08-05 南京航空航天大学 一种基于数据驱动的雷达通信混叠信号分离方法及设备

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2014048193A1 (zh) * 2012-09-28 2014-04-03 北京理工大学 一种用于舰艇编队情况下同型雷达同频干扰抑制方法
CN109917347A (zh) * 2019-04-10 2019-06-21 电子科技大学 一种基于时频域稀疏重构的雷达行人检测方法
EP3514571A1 (fr) * 2018-01-18 2019-07-24 Thales Procede de pistage d'une cible aerienne, et radar mettant en oeuvre un tel procede
CN111025256A (zh) * 2019-12-26 2020-04-17 湖南华诺星空电子技术有限公司 一种机载雷达的微弱生命体征信号的检测方法及系统
CN111551933A (zh) * 2020-04-29 2020-08-18 南京理工大学 基于稀疏表示理论的微动群目标isar成像方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2014048193A1 (zh) * 2012-09-28 2014-04-03 北京理工大学 一种用于舰艇编队情况下同型雷达同频干扰抑制方法
EP3514571A1 (fr) * 2018-01-18 2019-07-24 Thales Procede de pistage d'une cible aerienne, et radar mettant en oeuvre un tel procede
CN109917347A (zh) * 2019-04-10 2019-06-21 电子科技大学 一种基于时频域稀疏重构的雷达行人检测方法
CN111025256A (zh) * 2019-12-26 2020-04-17 湖南华诺星空电子技术有限公司 一种机载雷达的微弱生命体征信号的检测方法及系统
CN111551933A (zh) * 2020-04-29 2020-08-18 南京理工大学 基于稀疏表示理论的微动群目标isar成像方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
向虎 ; 李少东 ; 向龙 ; 陈文峰 ; 杨军 ; .欠采样下宽带自旋目标的快速高分辨成像方法.电子与信息学报.2018,(第11期),全文. *
张群 ; 屈筱钰 ; 李开明 ; 苏令华 ; .MMV模型下无源双基雷达低空目标微动特征提取.华中科技大学学报(自然科学版).2020,(第03期),全文. *

Also Published As

Publication number Publication date
CN113219416A (zh) 2021-08-06

Similar Documents

Publication Publication Date Title
CN108872962B (zh) 基于分数阶傅里叶变换的激光雷达微弱信号提取和分解方法
CN107894586B (zh) 一种基于同步压缩变换的激光雷达回波信号去噪方法
Wang et al. A reducing iteration orthogonal matching pursuit algorithm for compressive sensing
CN113219416B (zh) 一种基于雷达缺失采样的微动信号分离方法
CN111722188A (zh) 基于stft预分选的pri变换雷达信号分选方法
CN107390193B (zh) 基于多距离单元融合的调频连续波雷达飞机目标分类方法
CN104764714B (zh) 一种基于经验模态分解提高太赫兹频谱分辨率的方法
CN112731306B (zh) 基于CS和简化FrFT的UWB-LFM信号参数估计方法
Bai et al. Obtaining JTF-signature of space-debris from incomplete and phase-corrupted data
Liu et al. Underwater target recognition based on line spectrum and support vector machine
Li et al. SMSP jamming suppression method based on jamming reconstruction and kurtosis maximum
CN109598175A (zh) 一种基于多小波基函数和超正交前向回归的时频分析方法
CN107966687B (zh) 基于部分自相关谱的mimo雷达信号调制类型识别方法
CN113116320A (zh) 一种基于vmd的fmcw雷达生命信号检测方法
Chen et al. Study of threshold setting for rapid detection of multicomponent LFM signals based on the fourth-order origin moment of fractional spectrum
Yang et al. Time–frequency feature enhancement of moving target based on adaptive short-time sparse representation
CN111308426B (zh) 一种适用于单天线接收机的低信噪比周期调频信号检测与分离方法
CN113281714A (zh) 基于雷达微多普勒特征增强的鸟类目标探测方法
Li et al. A wideband reconnaissance receiver design based on real-time spectrum analysis technology
Qu et al. WGAN-GP-based synthetic radar spectrogram augmentation in human activity recognition
Wang et al. LFM signal perception based on Wavelet transform and Time-Frequency Technology
Cong et al. The parameter estimation of wideband LFM signal based on down-chirp and compressive sensing
Wang et al. Separation for Incomplete Micro-Doppler Signal via STVMD and Hungarian Algorithm
CN113238244B (zh) 一种fmcw激光测距差拍信号频率估计方法及系统
Wang et al. A Bayesian denoising mehod for complex radar signal with application to classification of human individuals

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