CN114814830A - 一种基于鲁棒主成分分析降噪的米波雷达低仰角测高方法 - Google Patents

一种基于鲁棒主成分分析降噪的米波雷达低仰角测高方法 Download PDF

Info

Publication number
CN114814830A
CN114814830A CN202210325727.1A CN202210325727A CN114814830A CN 114814830 A CN114814830 A CN 114814830A CN 202210325727 A CN202210325727 A CN 202210325727A CN 114814830 A CN114814830 A CN 114814830A
Authority
CN
China
Prior art keywords
matrix
low
signal
noise
rank
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.)
Pending
Application number
CN202210325727.1A
Other languages
English (en)
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 CN202210325727.1A priority Critical patent/CN114814830A/zh
Publication of CN114814830A publication Critical patent/CN114814830A/zh
Pending legal-status Critical Current

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
    • 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/88—Radar or analogous systems specially adapted for specific applications
    • G01S13/882—Radar or analogous systems specially adapted for specific applications for altimeters
    • 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
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
    • G06F18/213—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods
    • G06F18/2135—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods based on approximation criteria, e.g. principal component analysis
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
    • G06F2218/02—Preprocessing
    • G06F2218/04—Denoising

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Remote Sensing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Data Mining & Analysis (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Theoretical Computer Science (AREA)
  • Artificial Intelligence (AREA)
  • General Engineering & Computer Science (AREA)
  • Evolutionary Computation (AREA)
  • Evolutionary Biology (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Electromagnetism (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

本发明涉及一种基于鲁棒主成分分析降噪的米波雷达低仰角测高方法，包括步骤：将米波雷达的阵列接收信号建模为二维接收信号矩阵，得到阵列接收信号模型；采用基于鲁棒主成分分析降噪方法建立低秩信号矩阵和稀疏噪声矩阵的凸优化模型；使用改进的非精确增广拉格朗日乘子法求解凸优化模型，得到无噪声污染的最优低秩信号矩阵；利用最优低秩信号矩阵建立协方差矩阵，并结合协方差矩阵和合成导向矢量构建直达波仰角波达方向的最大似然估计，得到目标仰角估计值；利用目标仰角估计值计算目标高度。该方法克服现有技术在非均匀噪声、低信噪比、少快拍数情形下测角性能下降、角度分辨力不足等问题，具有较高的测角稳定性以及良好的角度分辨力。

Description

一种基于鲁棒主成分分析降噪的米波雷达低仰角测高方法
技术领域
本发明属于雷达信号处理技术领域，具体涉及一种基于鲁棒主成分分析降噪的米波雷达低仰角测高方法。
背景技术
米波雷达一度因其测量精度低、体型庞大，难以满足现代战场高精度、信息化与高机动的要求等缺陷被西方雷达界摒弃。反隐身的需求给米波雷达发展带来了新的契机，各国开始纷纷研制先进的米波雷达并将其做为骨干反隐身装备，米波雷达广泛的应用于防空预警任务中，起到中远程警戒、目标跟踪及指示的作用。然而，米波雷达受阵列孔径及其波长限制，波束宽度相对较宽，难以分辨位于同一个波束宽度内的直达波和反射波信号，严重多径效应的影响导致测角出现误差。因此，如何提高米波雷达低仰角的波达方向DOA估计是其测高技术的关键。
阵列超分辨技术因突破瑞利限的限制，所以被广泛应用在米波雷达低仰角测高中。其中，以多重信号分类算法为代表的空间平滑算法SSMUSIC在特定条件下均有很高的分辨力。空间平滑法是一种典型的降维方法，它将原阵列分为多个重叠的子阵，得到各子阵的协方差矩阵后对其求平均以便恢复协方差矩阵的秩。然而该算法是以降低阵列有效孔径为代价，低仰角范围内测角稳定性较低，难以满足低信噪比、少快拍数的要求，并且难以适应非均匀高斯噪声环境。以最大似然类算法为代表的合成导向矢量最大似然SVML估计算法是另一种性能优良的参数化DOA估计算法，该算法均不受信号间相关性影响，对阵列形状没有特殊要求，其对多目标进行DOA估计时都需要多维搜索，多维搜索的快速实现需要利用天线架高、实际阵地反射面的反射系数等先验信息将一个复杂的多维搜索转化为多个简单的一维搜索来减少运算量。但是此类算法在非均匀高斯噪声的影响下，测角性能下降，难以满足低信噪比和少快拍条件。
综上，现有非均匀高斯噪声背景下米波雷达低仰角测高算法存在角度分辨力不高、测角估计精度降低，难以在低信噪比和少快拍数条件下保持良好的DOA估计性能的问题。
发明内容
为了解决现有技术中存在的上述问题，本发明提供了一种基于鲁棒主成分分析降噪的米波雷达低仰角测高方法。本发明要解决的技术问题通过以下技术方案实现：
本发明实施例提供了一种基于鲁棒主成分分析降噪的米波雷达低仰角测高方法，包括步骤：
S1、将米波雷达的阵列接收信号建模为二维接收信号矩阵，得到阵列接收信号模型；
S2、采用基于鲁棒主成分分析降噪方法对所述阵列接收信号模型分解为低秩信号矩阵和稀疏噪声矩阵，并建立所述低秩信号矩阵和所述稀疏噪声矩阵的凸优化模型；
S3、使用改进的非精确增广拉格朗日乘子法求解所述凸优化模型，交替迭代更新所述低秩信号矩阵和所述稀疏噪声矩阵，得到无噪声污染的最优低秩信号矩阵；
S4、利用所述最优低秩信号矩阵建立协方差矩阵，并结合所述协方差矩阵和合成导向矢量构建直达波仰角波达方向的最大似然估计，得到目标仰角估计值；
S5、利用所述目标仰角估计值计算目标高度。
在本发明的一个实施例中，所述阵列接收信号模型为：
X＝A(θ)S+N
其中，X＝[x(1),…x(L)]为M×L维的接收信号矩阵，N＝[n(1),…,n(L)]为M×L维的非均匀高斯噪声矩阵，S＝[s1(t),…,sK(t)]T为K×L维的信号复包络矩阵，A(θ)＝[a(θ1),…a(θK)]为M×K维的阵列流型矩阵，
Figure BDA0003573345110000031
为导向矢量矩阵，θk为入射角范围，θk∈[-90°,90°]，M为阵元个数，K＜M＜L，L表示采样快拍数，d为阵元间距，λ0为入射信号波长。
在本发明的一个实施例中，步骤S2包括：
S21、采用基于鲁棒主成分分析降噪方法将所述阵列接收信号模型分解为所述低秩信号矩阵和所述稀疏噪声矩阵，并建立所述低秩信号矩阵和所述系稀疏噪声矩阵的低秩优化模型：
Figure BDA0003573345110000032
其中，P为低秩信号矩阵，
Figure BDA0003573345110000033
S22、引入折中因子，将所述低秩优化模型由双目标优化问题转化为单目标优化问题，得到转化后的优化模型：
Figure BDA0003573345110000034
其中，λ为折中因子，λ＞0，||·||0为矩阵的l0范数；
S23、将所述转化后的优化模型松弛为所述凸优化模型：
Figure BDA0003573345110000041
其中，||·||*为矩阵核范数，||·||1为矩阵的l1范数。
在本发明的一个实施例中，所述折中因子的取值范围为：
Figure BDA0003573345110000042
其中，M为阵元个数，L为采样快拍数。
在本发明的一个实施例中，步骤S3包括：
S31、获取阵列接收信号、阵元数、快拍数和折中因子；
S32、初始化设置拉格朗日乘子、噪声矩阵、收敛因子、放大因子、收敛判断因子和迭代次数；
S33、对目标矩阵进行奇异值分解：
(U,∑,V)＝svd(X-Nk+Yk/μk)
其中，svd表示矩阵X-Nk+1+Yk/μk的奇异值分解，U为矩阵X-Nk+1+Yk/μk的左奇异矩阵，∑为非零对角矩阵，其对角线上的元素为矩阵X-Nk+1+Yk/μk的奇异值，V为矩阵X-Nk+1+Yk/μk对应的右奇异矩阵，X为接收信号矩阵，Nk为稀疏噪声矩阵，Yk为拉格朗日乘子，μk为惩罚参数；
S34、根据所述奇异值分解结果更新所述低秩信号矩阵：
Figure BDA0003573345110000043
其中，Pk+1为更新后的低秩信号矩阵，S*(·)为奇异值收缩因子或软阈值算子，
Figure BDA0003573345110000044
σij为矩阵∑上第i行第j列的奇异值，μk为惩罚参数；
S35、利用更新后的低秩信号矩阵更新所述稀疏噪声矩阵：
Figure BDA0003573345110000045
其中，Nk+1为更新后的稀疏噪声矩阵，
Figure BDA0003573345110000051
λ为入射信号波长，μk为惩罚参数，X为接收信号矩阵，Pk+1为更新后的低秩信号矩阵，Yk为拉格朗日乘子；
S36、利用更新后的稀疏噪声矩阵更新拉格朗日乘子：
Figure BDA0003573345110000052
其中，Yk+1为更新后的拉格朗日乘子，Yk为拉格朗日乘子，μk为惩罚参数，X为接收信号矩阵，
Figure BDA0003573345110000053
为低秩信号矩阵的最优解，
Figure BDA0003573345110000054
为稀疏噪声矩阵的最优解；
S37、利用更新后的拉格朗日乘子更新惩罚参数：
μk+1＝min(ρμk,μk.max)
其中，μk+1为更新后的惩罚参数，ρ＞1为常数，μk为惩罚参数，μk.max为k次迭代中最大的惩罚参数；
S38、判断是否满足迭代终止条件，若否，则令迭代次数加1次，返回步骤S33继续迭代，若是，则循环结束，输出所述最优低秩信号矩阵。
在本发明的一个实施例中，所述迭代终止条件为：
Figure BDA0003573345110000055
其中，X为接收信号矩阵，Pk+1为更新后的低秩信号矩阵，Nk+1为更新后的稀疏噪声矩阵。
在本发明的一个实施例中，所述迭代终止条件为：迭代次数达到预设迭代次数。
在本发明的一个实施例中，步骤S4包括：
S41、利用所述最优低秩信号矩阵建立协方差矩阵：
Figure BDA0003573345110000061
其中，RP为协方差矩阵，N为阵元个数，E为，P*为最优低秩信号矩阵；
S42、确定角度搜索范围，利用先验信息将直达波强相干信号与反射波强相干信号合成为一个信号，构建所述合成导向矢量：
asum(θ)＝a(θ)-exp(-j4πhasinθ/λ)a(-θ)
其中，反射系数设置为-1，asum(θ)为合成导向矢量，a(θ)为直达波导向矢量，a(-θ)为对应的反射波导向矢量，ha为雷达架高，θ为入射角度，λ为入射信号波长；
S43、利用所述合成导向矢量构建最大似然投影矩阵：
Figure BDA0003573345110000062
其中，Psum(θ)为最大似然投影矩阵，asum(θ)为合成导向矢量；
S44、利用所述协方差矩阵和所述最大似然投影矩阵构建所述直达波仰角波达方向的最大似然估计，得到所述目标仰角估计值：
Figure BDA0003573345110000063
其中，θd为目标仰角估计值，Psum(θ)为最大似然投影矩阵，RP为协方差矩阵。
在本发明的一个实施例中，当不考虑地球曲率半径时，所述目标高度为：
ht＝Rdsinθd+ha
其中，ht为目标高度，Rd为已知的目标距离，θd为目标仰角估计值，ha为雷达架高。
在本发明的一个实施例中，当考虑地球曲率半径时，所述目标高度为：
Figure BDA0003573345110000071
其中，ht为目标高度，Re为地球的曲率半径，Re＝4R0/3，R0为地球平均半径，R0＝6370km，Rd为已知的目标距离，θd为目标仰角估计值。
与现有技术相比，本发明的有益效果：
本发明的米波雷达低仰角测高方法，采用鲁棒主成分分析降噪凸优化模型不断地优化阵列接收信号矩阵，降低噪声分量的同时又使得信号矩阵最优，从而逐渐消除阵列信号处理过程中噪声的影响，特别是非均匀高斯白噪声的影响；然后利用无噪声干扰的最优低秩信号协方差矩阵结合合成导向矢量最大似然进行DOA估计，在提高测角稳定性的同时又实现了超分辨，克服现有技术在非均匀噪声、低信噪比、少快拍数情形下测角性能下降、角度分辨力不足等问题，具有较高的测角稳定性以及良好的角度分辨力。
附图说明
图1为本发明实施例提供的一种基于鲁棒主成分分析降噪的米波雷达低仰角测高方法的流程示意图；
图2为本发明实施例提供的米波雷达接收信号模型示意图；
图3为本发明实施例提供的一种改进的非精确增广拉格朗日乘子算法的流程示意图；
图4a-图4b为本发明实施例提供的几种算法空间谱在高斯噪声和非均匀高斯噪声背景下DOA估计示意图；
图5为本发明实施例提供的几种算法测角性能在非均匀高斯噪声背景下随信噪比的变化示意图；
图6为本发明实施例提供的几种算法测角性能在非均匀高斯噪声背景下随快拍数的变化示意图。
具体实施方式
下面结合具体实施例对本发明做进一步详细的描述，但本发明的实施方式不限于此。
实施例一
请参加图1，图1为本发明实施例提供的一种基于鲁棒主成分分析降噪的米波雷达低仰角测高方法的流程示意图，该方法包括步骤：
S1、将米波雷达的阵列接收信号建模为二维接收信号矩阵，得到阵列接收信号模型。
请参见图2，图2为本发明实施例提供的米波雷达接收信号模型示意图。假设有K个远场窄带非相干源信号从不同的方向入射到垂直放置的M个各向同性阵元组成的均匀线阵，阵元间距均为d，入射信号波长为λ，d≤λ/2，入射角范围为θk∈[-90°,90°]，k∈[1,2,…,K]，则t时刻阵列接收信号模型可表示为：
Figure BDA0003573345110000081
其中，x(t)为接收信号矢量，
Figure BDA0003573345110000082
为第k个信源的阵列导向矢量，sk(t)为信号复包络矩阵，n(t)＝[n1(t),…,nM(t)]为互不相关的非均匀高斯白噪声。
进一步的，将t时刻阵列接收信号模型改写为：
x(t)＝A(θ)s(t)+n(t) (2)
其中，A(θ)＝[a(θ1),…a(θK)]为M×K维的阵列流型矩阵，一般情况下阵元数大于信源数，a(θ)＝[a1(θ),…,aM(θ)]T。
对于L次采样快拍，阵列接收信号模型可表示为：
X＝A(θ)S+N (3)
其中，X＝[x(1),…x(L)]为M×L维的接收信号矩阵，N＝[n(1),…,n(L)]为M×L维的非均匀高斯噪声矩阵，S＝[s1(t),…,sK(t)]T为K×L维的信号复包络矩阵，A(θ)＝[a(θ1),…a(θK)]为M×K维的阵列流型矩阵，
Figure BDA0003573345110000091
为导向矢量矩阵，θk为入射角范围，θk∈[-90°,90°]，M为阵元个数，K＜M＜L，L表示采样快拍数，d为阵元间距，λ0为入射信号波长。
本实施例的目的是从阵列接收的混合信号中准确估计出A(θ)S待处理的目标信号，采用所得无噪声污染的最优低秩目标信号进行DOA估计，从而实现相干信号和非均匀噪声条件下DOA估计性能差和分辨率低的问题。针对上述目的，本实施例进行以下步骤S2～S6的处理。
S2、采用基于鲁棒主成分分析降噪方法将所述阵列接收信号模型分解为低秩信号矩阵和稀疏噪声矩阵，并建立所述低秩信号矩阵和所述系稀疏噪声矩阵的凸优化模型。具体包括：
S21、采用基于鲁棒主成分分析降噪方法将所述阵列接收信号模型分解为所述低秩信号矩阵和所述稀疏噪声矩阵，并建立所述低秩信号矩阵和所述系稀疏噪声矩阵的低秩优化模型。
具体的，为了消除阵列接收信号中影响波达角DOA估计精度的非均匀高斯噪声干扰，本实施例对阵列接收信号进行鲁棒主成分分析降噪分解，通过求解鲁棒主成分分析凸优化模型得到无噪声干扰的最优低秩信号矩阵。
由于阵列流型矩阵A(θ)是列满秩矩阵并且很容易证明rank(A(θ)S)＝r≤K，因此A(θ)S是去除噪声后的低秩信号矩阵。此外非均匀高斯白噪声背景下，噪声协方差矩阵Rn＝σn 2I，σn 2＝[σ1 2,…,σM 2]，rank(Rn)＝M是满秩矩阵，且除对角线元素外其余元素均为零，满足稀疏特性。因此，令
Figure BDA0003573345110000101
为了得到最优化的低秩矩阵P而免除噪声的干扰，建立鲁棒主成分分析降噪的凸优化模型，从而实现信号空间与噪声空间的分离。
进一步的，采用鲁棒主成分分析降噪方法建立所述低秩信号矩阵和所述系稀疏噪声矩阵的低秩优化模型：
Figure BDA0003573345110000102
其中，P为低秩信号矩阵，
Figure BDA0003573345110000103
N为稀疏噪声矩阵，X为接收信号矩阵。
S22、引入折中因子，将所述低秩优化模型由双目标优化问题转化为单目标优化问题，得到转化后的优化模型。
具体的，对于式(4)，经典的PCA算法求取最优低秩矩阵P，对数据矩阵N进行SVD分解即可；但是当N为稀疏大噪声矩阵时，PCA方法将不再适用。此时求取最优低秩矩阵P是双目标优化问题，希望低秩矩阵P的秩最小的同时噪声N也最小，因此，本实施例通过引入折中因子λ(＞0)进而控制稀疏噪声的权重，将低秩优化模型的双目标优化问题转化为单目标优化问题从而简化求解，转化后的优化模型为：
Figure BDA0003573345110000104
其中，λ为折中因子，λ＞0，||·||0为矩阵的l0范数(矩阵中非零元素个数)。
S23、将所述转化后的优化模型松弛为所述凸优化模型。
具体的，式(5)的优化问题属于NP难问题，可以将上述NP问题转化为凸优化问题求解。由于矩阵的核范数是矩阵秩的包络，矩阵的l0范数在一定条件下可以和l1范数等价，所以转化后的优化模型可以松弛为基于鲁棒主成分分析的凸优化模型：
Figure BDA0003573345110000111
其中，||·||*为矩阵核范数(矩阵奇异值的和)，||·||1为矩阵的l1范数(列向量绝对值之和的最大值)。
在一个具体实施例中，折中因子λ的取值范围为：
Figure BDA0003573345110000112
其中，M为阵元个数，L为采样快拍数。
本实施例中，将最优低秩信号矩阵的求解转化为鲁棒主成分分析的凸优化问题，在提高信噪比的同时也降低噪声对波达方向DOA估计精度的影响。
S3、使用改进的非精确增广拉格朗日乘子法求解所述凸优化模型，交替迭代更新所述低秩信号矩阵和所述稀疏噪声矩阵，得到无噪声污染的最优低秩信号矩阵。
具体的，首先构造增广拉格朗日函数为：
Figure BDA0003573345110000113
其中，Y为拉格朗日乘子；μ为惩罚参数也称收敛因子，μ0＞0且为比较小的正数，例如μ0＜10-6；X为接收信号矩阵；P为低秩信号矩阵；N为稀疏噪声矩阵；λ为折中因子；<·>表示矩阵的欧式内积；||·||F是矩阵的Frobenius范数，
Figure BDA0003573345110000121
如果
Figure BDA0003573345110000122
σi为其矩阵的奇异值。
根据增广拉格朗日函数表达式，将鲁棒主成分分析的凸优化模型的解描述为：
Figure BDA0003573345110000123
其中，(P*,N*)为矩阵P,N的最优解。
由于精确拉格朗日乘子法在计算过程中，每次迭代都包含计算量较大的奇异值分解，计算复杂度高、耗时长；而非精确拉格朗日乘子法在这方面有所改善，它不需要求
Figure BDA0003573345110000124
的精确解，而是求出近似解。因此，本实施例采用改进的非精确增广拉格朗日乘子法求解所述凸优化模型，以使得计算复杂度降低。
进一步的，假设(Pk,Nk,Yk,μk)代表(P,N,Y,μ)在第k次迭代中的值，那么可通过迭代求解下面2个子问题，进而求得最优解(P*,N*)。这2个子问题表述为：
Figure BDA0003573345110000125
其中，U为矩阵X-Nk+1+Yk/μk的左奇异矩阵，∑为非零对角矩阵，其对角线上的元素为矩阵X-Nk+1+Yk/μk的奇异值，V为矩阵X-Nk+1+Yk/μk对应的右奇异矩阵，
Figure BDA0003573345110000131
是一个矩阵算子，具体而言是矩阵X-Nk+1+Yk/μk的SVD分解；Sε(∑)为奇异值收缩因子或软阈值算子，
Sε(Σ)＝sgn(σij)max(|σij|-ε,0) (11)
Sε(∑)具体而言为：
Figure BDA0003573345110000132
其中，σij为矩阵∑上第i行第j列的奇异值。
如果Pk+1，Nk+1分别收敛于
Figure BDA0003573345110000133
则拉格朗日乘子的更新公式为：
Figure BDA0003573345110000134
其中，Nk+1为更新后的稀疏噪声矩阵，λ为入射信号波长，μk为惩罚参数，X为接收信号矩阵，Pk+1为更新后的低秩信号矩阵，Yk为拉格朗日乘子。
由上述过程可见，该算法的计算量主要在于SVD分解。同时，为保证迭代求解的(Pk,Nk)收敛到RPCA问题的最优解(P*,N*)，序列{μk}为非增序列且
Figure BDA0003573345110000135
因此选择合适的μk时应尽可能保证RPCA收敛且SVD分解总次数最小。
本实施例中，惩罚参数μ的更新公式为：
μk+1＝min(ρμk,μk.max) (14)
其中，μk+1为更新后的惩罚参数，ρ＞1为常数，μk为惩罚参数，μk.max为k次迭代中最大的惩罚参数。
本实施例中，惩罚参数μk对改进的非精确增广拉格朗日乘子算法的收敛速度和稳态误差起着关键性作用。在更新迭代低秩信号矩阵与稀疏噪声矩阵的过程中，选择合适的惩罚参数μk可以提高***在收敛速度和稳态误差方面的平衡性。
最后，判断迭代是否满足终止条件；当判断迭代满足迭代终止条件，则迭代终止；当判断迭代不满足迭代终止条件，则令k＝k+1继续进行迭代。
综上，请参见图3，图3为本发明实施例提供的一种改进的非精确增广拉格朗日乘子算法的流程示意图。该算法具体包括步骤：
S31、获取阵列接收信号、阵元数、快拍数和折中因子。
具体的，获取如下参数：阵列接收信号
Figure BDA0003573345110000141
阵元数M、快拍数L和折中因子λ。
S32、初始化设置拉格朗日乘子、噪声矩阵、收敛因子、放大因子、收敛判断因子和迭代次数。
具体的，设置拉格朗日乘子和噪声矩阵Y0＝N0＝0；收敛因子μ0＞0且较小，例如μ0＜10-6，μ0＝1.5/||X||2，μk.max＝104；放大因子ρ＝1.1；收敛判断因子ε＝10-6；迭代次数k＝0。
在初始化设置完成后，更新矩阵参数，具体包括：
S33、对目标矩阵进行奇异值分解。
具体的，对目标矩阵X-Nk+Yk/μk进行奇异值分解：
(U,∑,V)＝svd(X-Nk+Yk/μk) (15)
其中，svd表示矩阵X-Nk+1+Yk/μk的奇异值分解，U为矩阵X-Nk+1+Yk/μk的左奇异矩阵，Σ为非零对角矩阵，其对角线上的元素为矩阵X-Nk+1+Yk/μk的奇异值，V为矩阵X-Nk+1+Yk/μk对应的右奇异矩阵，X为接收信号矩阵，Nk为稀疏噪声矩阵，Yk为拉格朗日乘子，μk为惩罚参数。
S34、根据所述奇异值分解结果更新所述低秩信号矩阵：
Figure BDA0003573345110000151
其中，Pk+1为更新后的低秩信号矩阵，S*(·)为奇异值收缩因子或软阈值算子，
Figure BDA0003573345110000152
σij为矩阵∑上第i行第j列的奇异值，μk为惩罚参数。
S35、利用更新后的低秩信号矩阵更新所述稀疏噪声矩阵：
Figure BDA0003573345110000153
其中，Nk+1为更新后的稀疏噪声矩阵，
Figure BDA0003573345110000154
λ为入射信号波长，μk为惩罚参数，X为接收信号矩阵，Pk+1为更新后的低秩信号矩阵，Yk为拉格朗日乘子。
S36、利用更新后的稀疏噪声矩阵更新拉格朗日乘子，如公式(13)所示。
S37、利用更新后的拉格朗日乘子更新惩罚参数，如公式(14)所示。
S38、判断是否满足迭代终止条件，若否，则令迭代次数加1次即k＝k+1，返回步骤S33继续迭代，若是，则循环结束，输出所述最优低秩信号矩阵。
在一个具体实施例中，迭代终止条件为：
Figure BDA0003573345110000155
其中，X为接收信号矩阵，Pk+1为更新后的低秩信号矩阵，Nk+1为更新后的稀疏噪声矩阵。
在另一个实施例中，迭代终止条件为：迭代次数达到预设迭代次数。具体的，当通过迭代终止条件无法得到最优低秩信号矩阵P*时，可以设置最大迭代次数来结束更新过程，得到最终无噪声的矩阵P*。
S4、利用所述最优低秩信号矩阵建立协方差矩阵，并结合所述协方差矩阵和合成导向矢量构建直达波仰角波达方向的最大似然估计，得到目标仰角估计值。具体包括步骤：
S41、利用所述最优低秩信号矩阵建立协方差矩阵：
Figure BDA0003573345110000161
其中，RP为协方差矩阵，N为阵元个数，E为，P*为最优低秩信号矩阵。
S42、确定角度搜索范围，利用先验信息将直达波强相干信号与反射波强相干信号合成为一个信号，构建所述合成导向矢量。
具体的，确定角度搜索范围并设置合适的角度步进值，利用衰减系数ρ、雷达架高ha等先验信息将直达波与反射波这两个强相干信号矢量合成为一个信号矢量，从而构建得合成导向矢量asum(θ)。通常当ha较小时，可近似认为直达角和反射角互为相反数，反射系数ρ＝-1，此时有合成导向矢量asum(θ)为：
asum(θ)＝a(θ)-exp(-j4πhasinθ/λ)a(-θ) (20)
其中，asum(θ)为合成导向矢量，a(θ)为直达波导向矢量，a(-θ)为对应的反射波导向矢量，ha为雷达架高，θ为入射角度，λ为入射信号波长。
S43、利用所述合成导向矢量构建最大似然投影矩阵：
Figure BDA0003573345110000162
其中，Psum(θ)为最大似然投影矩阵，asum(θ)为合成导向矢量。
S44、利用所述协方差矩阵和所述最大似然投影矩阵构建所述直达波仰角波达方向的最大似然估计，得到所述目标仰角估计值。
Figure BDA0003573345110000171
其中，θd为目标仰角估计值，Psum(θ)为最大似然投影矩阵，RP为协方差矩阵。
S5、利用所述目标仰角估计值计算目标高度。
在一个具体实施例中，根据谱峰搜索得到目标仰角估计值θd和已知的目标距离Rd、雷达架高ha等先验信息，根据米波雷达低仰角测高几何模型得到不考虑地球曲率时目标的高度：
ht＝Rdsinθd+ha (23)
其中，ht为目标高度，Rd为已知的目标距离，θd为目标仰角估计值，ha为雷达架高。
在一个具体实施例中，对于远距离目标，地球曲率不能忽略，Re为地球的曲率半径，考虑大气折射Re＝4R0/3，R0为地球平均半径，R0＝6370km，此时的目标高度为：
Figure BDA0003573345110000172
其中，ht为目标高度，Re为地球的曲率半径，Re＝4R0/3，R0为地球平均半径，R0＝6370km，Rd为已知的目标距离，θd为目标仰角估计值。
本实施例的米波雷达低仰角测高方法，采用鲁棒主成分分析降噪凸优化模型不断地优化阵列接收信号矩阵，降低噪声分量的同时又使得信号矩阵最优，从而逐渐消除阵列信号处理过程中噪声的影响，特别是非均匀高斯白噪声的影响；然后利用无噪声干扰的最优低秩信号协方差矩阵结合合成导向矢量最大似然进行DOA估计，在提高测角稳定性的同时又实现了超分辨，克服现有技术在非均匀噪声、低信噪比、少快拍数情形下测角性能下降、角度分辨力不足等问题，具有较高的测角稳定性以及良好的角度分辨力，在相干条件下具有良好的DOA估计性能。
实施例二
在实施例一的基础上，为了验证本实施例基于鲁棒主成分分析降噪的米波雷达低仰角测高方法的有效性，本实施例通过以下仿真实验进行说明。
仿真条件：考虑8阵元的均匀线阵，阵元间距为半波长，2个远场窄带信号入射到该天线阵列，信噪比均为10dB，采样快拍数为100，分别考虑噪声为高斯白噪声和非均匀高斯噪声情况，高斯噪声背景下噪声为0均值，方差为1的高斯白噪声；非均匀高斯噪声背景下噪声的协方差为Q＝diag{2,7,11,5,3.5,13,5.5,1.5}。
仿真实验1：
对比本实施例算法与SSMUSIC算法、SVML算法的空间谱在高斯噪声和非均匀高斯噪声背景下的DOA估计。
假设2个远场窄带信号的入射角度分别为-3.4°和3.4°，仿真结果如图4a-图4b所示。图4a-图4b为本发明实施例提供的几种算法空间谱在高斯噪声和非均匀高斯噪声背景下DOA估计示意图，其中，图4a为高斯背景下各算法的空间谱，图4b为非均匀高斯噪声背景下各算法的空间谱图。
从图中可看出当目标信号处于高斯白噪声背景下时，各算法均具有良好的DOA估计效果，但SVML算法和本实施例算法的分辨力更高，另外本实施例算法相较SVML算法的主峰宽度略小，旁瓣电平略低；当目标信号处于非均匀高斯噪声背景下时，SSMUSIC算法已经不能对这两个相干信号进行有效分辨，本实施例算法依然可以准确地区分出两个相干信号，且比SVML算法的谱峰更为尖锐，旁瓣电平更低。因此，本实施例所提算法在非均匀噪声背景下具有更高的角度分辨力。
仿真实验2：
对比本实施例算法与SSMUSIC算法、SVML算法的测角性能在非均匀高斯噪声背景下随信噪比的变化情况。
假设2个远场窄带信号的入射角度分别为-3.4°和3.4°，信噪比由-10dB按步长5dB向30dB变化，采样快拍数为100。仿真进行100次蒙特卡洛实验，结果取平均，仿真结果请参见图5，图5为本发明实施例提供的几种算法测角性能在非均匀高斯噪声背景下随信噪比的变化示意图。
从图5可看出，在非均匀高斯噪声背景下，各算法的测角精度均随信噪比SNR的增加逐渐降低，其中SSMUSIC算法测角误差最大，SVML次之，而本文算法测角误差最小。当信噪比SNR大于0dB小于10dB时，SSMUSIC和SVML算法测角误差较大导致其测角失效，而本文所提算法依然具有较高的测角精度。由此说明本文所提算法在低信噪比条件下具备更优的测角性能。
仿真实验3：
对比本实施例算法与SSMUSIC算法、SVML算法的测角性能(均方根误差)在非均匀高斯噪声背景下随快拍数的变化情况。
假设2个远场窄带信号的入射角度分别为-3.4°和3.4°，快拍数由50按步长50向300变化，信噪比均为15dB。仿真进行100次蒙特卡洛实验，结果取平均。仿真结果请参见图6，图6为本发明实施例提供的几种算法测角性能在非均匀高斯噪声背景下随快拍数的变化示意图。
从图6可看出在在非均匀高斯噪声背景下，本实施例所提算法的测角均方根误差最小，即使在较少快拍数条件下依然能具有较高的测角精度，测角性能稳定。
综上所述，本实施例提出的基于鲁棒主成分分析降噪的米波雷达低仰角测高方法能够有效估计非均匀噪声环境下的米波雷达低仰角目标的DOA，对非均匀噪声背景下的低信噪比、少快拍数条件更具鲁棒性，且可以实现相干信号的超分辨和高精度的测角能力，进一步证明了该算法的测角有效性。
以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单推演或替换，都应当视为属于本发明的保护范围。

Claims (10)

1.一种基于鲁棒主成分分析降噪的米波雷达低仰角测高方法，其特征在于，包括步骤：
S1、将米波雷达的阵列接收信号建模为二维接收信号矩阵，得到阵列接收信号模型；
S2、采用基于鲁棒主成分分析降噪方法对所述阵列接收信号模型分解为低秩信号矩阵和稀疏噪声矩阵，并建立所述低秩信号矩阵和所述稀疏噪声矩阵的凸优化模型；
S3、使用改进的非精确增广拉格朗日乘子法求解所述凸优化模型，交替迭代更新所述低秩信号矩阵和所述稀疏噪声矩阵，得到无噪声污染的最优低秩信号矩阵；
S4、利用所述最优低秩信号矩阵建立协方差矩阵，并结合所述协方差矩阵和合成导向矢量构建直达波仰角波达方向的最大似然估计，得到目标仰角估计值；
S5、利用所述目标仰角估计值计算目标高度。
2.根据权利要求1所述的基于鲁棒主成分分析降噪的米波雷达低仰角测高方法，其特征在于，所述阵列接收信号模型为：
X＝A(θ)S+N
其中，X＝[x(1),…x(L)]为M×L维的接收信号矩阵，N＝[n(1),…,n(L)]为M×L维的非均匀高斯噪声矩阵，S＝[s1(t),…,sK(t)]T为K×L维的信号复包络矩阵，A(θ)＝[a(θ1),…a(θK)]为M×K维的阵列流型矩阵，
Figure FDA0003573345100000011
为导向矢量矩阵，θk为入射角范围，θk∈[-90°,90°]，M为阵元个数，K＜M＜L，L表示采样快拍数，d为阵元间距，λ0为入射信号波长。
3.根据权利要求1所述的基于鲁棒主成分分析降噪的米波雷达低仰角测高方法，其特征在于，步骤S2包括：
S21、采用基于鲁棒主成分分析降噪方法将所述阵列接收信号模型分解为所述低秩信号矩阵和所述稀疏噪声矩阵，并建立所述低秩信号矩阵和所述系稀疏噪声矩阵的低秩优化模型：
Figure FDA0003573345100000021
其中，rank表示矩阵求秩，P为低秩信号矩阵，
Figure FDA0003573345100000022
S22、引入折中因子，将所述低秩优化模型由双目标优化问题转化为单目标优化问题，得到转化后的优化模型：
Figure FDA0003573345100000023
其中，λ为折中因子，λ＞0，||·||0为矩阵的l0范数；
S23、将所述转化后的优化模型松弛为所述凸优化模型：
Figure FDA0003573345100000024
其中，||·||*为矩阵核范数，||·||1为矩阵的l1范数。
4.根据权利要求3所述的基于鲁棒主成分分析降噪的米波雷达低仰角测高方法，其特征在于，所述折中因子的取值范围为：
Figure FDA0003573345100000025
其中，M为阵元个数，L为采样快拍数。
5.根据权利要求1所述的基于鲁棒主成分分析降噪的米波雷达低仰角测高方法，其特征在于，步骤S3包括：
S31、获取阵列接收信号、阵元数、快拍数和折中因子；
S32、初始化设置拉格朗日乘子、噪声矩阵、收敛因子、放大因子、收敛判断因子和迭代次数；
S33、对目标矩阵进行奇异值分解：
(U,∑,V)＝svd(X-Nk+Yk/μk)
其中，svd表示矩阵X-Nk+1+Yk/μk的奇异值分解，U为矩阵X-Nk+1+Yk/μk的左奇异矩阵，∑为非零对角矩阵，其对角线上的元素为矩阵X-Nk+1+Yk/μk的奇异值，V为矩阵X-Nk+1+Yk/μk对应的右奇异矩阵，X为接收信号矩阵，Nk为稀疏噪声矩阵，Yk为拉格朗日乘子，μk为惩罚参数；
S34、根据所述奇异值分解结果更新所述低秩信号矩阵：
Figure FDA0003573345100000031
其中，Pk+1为更新后的低秩信号矩阵，S*(·)为奇异值收缩因子或软阈值算子，
Figure FDA0003573345100000032
σij为矩阵∑上第i行第j列的奇异值，μk为惩罚参数；
S35、利用更新后的低秩信号矩阵更新所述稀疏噪声矩阵：
Figure FDA0003573345100000033
其中，Nk+1为更新后的稀疏噪声矩阵，
Figure FDA0003573345100000034
λ为入射信号波长，μk为惩罚参数，X为接收信号矩阵，Pk+1为更新后的低秩信号矩阵，Yk为拉格朗日乘子；
S36、利用更新后的稀疏噪声矩阵更新拉格朗日乘子：
Figure FDA0003573345100000035
其中，Yk+1为更新后的拉格朗日乘子，Yk为拉格朗日乘子，μk为惩罚参数，X为接收信号矩阵，
Figure FDA0003573345100000036
为低秩信号矩阵的最优解，
Figure FDA0003573345100000037
为稀疏噪声矩阵的最优解；
S37、利用更新后的拉格朗日乘子更新惩罚参数：
μk+1＝min(ρμk,μk.max)
其中，μk+1为更新后的惩罚参数，ρ＞1为常数，μk为惩罚参数，μk.max为k次迭代中最大的惩罚参数；
S38、判断是否满足迭代终止条件，若否，则令迭代次数加1次，返回步骤S33继续迭代，若是，则循环结束，输出所述最优低秩信号矩阵。
6.根据权利要求5所述的基于鲁棒主成分分析降噪的米波雷达低仰角测高方法，其特征在于，所述迭代终止条件为：
Figure FDA0003573345100000041
其中，X为接收信号矩阵，Pk+1为更新后的低秩信号矩阵，Nk+1为更新后的稀疏噪声矩阵。
7.根据权利要求5所述的基于鲁棒主成分分析降噪的米波雷达低仰角测高方法，其特征在于，所述迭代终止条件为：迭代次数达到预设迭代次数。
8.根据权利要求1所述的基于鲁棒主成分分析降噪的米波雷达低仰角测高方法，其特征在于，步骤S4包括：
S41、利用所述最优低秩信号矩阵建立协方差矩阵：
Figure FDA0003573345100000042
其中，RP为协方差矩阵，N为阵元个数，E为，P*为最优低秩信号矩阵；
S42、确定角度搜索范围，利用先验信息将直达波强相干信号与反射波强相干信号合成为一个信号，构建所述合成导向矢量：
asum(θ)＝a(θ)-exp(-j4πhasinθ/λ)a(-θ)
其中，反射系数设置为-1，asum(θ)为合成导向矢量，a(θ)为直达波导向矢量，a(-θ)为对应的反射波导向矢量，ha为雷达架高，θ为入射角度，λ为入射信号波长；
S43、利用所述合成导向矢量构建最大似然投影矩阵：
Figure FDA0003573345100000051
其中，Psum(θ)为最大似然投影矩阵，asum(θ)为合成导向矢量；
S44、利用所述协方差矩阵和所述最大似然投影矩阵构建所述直达波仰角波达方向的最大似然估计，得到所述目标仰角估计值：
Figure FDA0003573345100000052
其中，θd为目标仰角估计值，Psum(θ)为最大似然投影矩阵，RP为协方差矩阵。
9.根据权利要求1所述的基于鲁棒主成分分析降噪的米波雷达低仰角测高方法，其特征在于，当不考虑地球曲率半径时，所述目标高度为：
ht＝Rdsinθd+ha
其中，ht为目标高度，Rd为已知的目标距离，θd为目标仰角估计值，ha为雷达架高。
10.根据权利要求1所述的基于鲁棒主成分分析降噪的米波雷达低仰角测高方法，其特征在于，当考虑地球曲率半径时，所述目标高度为：
Figure FDA0003573345100000053
其中，ht为目标高度，Re为地球的曲率半径，Re＝4R0/3，R0为地球平均半径，R0＝6370km，Rd为已知的目标距离，θd为目标仰角估计值。
CN202210325727.1A 2022-03-30 2022-03-30 一种基于鲁棒主成分分析降噪的米波雷达低仰角测高方法 Pending CN114814830A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210325727.1A CN114814830A (zh) 2022-03-30 2022-03-30 一种基于鲁棒主成分分析降噪的米波雷达低仰角测高方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210325727.1A CN114814830A (zh) 2022-03-30 2022-03-30 一种基于鲁棒主成分分析降噪的米波雷达低仰角测高方法

Publications (1)

Publication Number Publication Date
CN114814830A true CN114814830A (zh) 2022-07-29

Family

ID=82532730

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210325727.1A Pending CN114814830A (zh) 2022-03-30 2022-03-30 一种基于鲁棒主成分分析降噪的米波雷达低仰角测高方法

Country Status (1)

Country Link
CN (1) CN114814830A (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117289262A (zh) * 2023-11-27 2023-12-26 中南大学 穿墙雷达目标检测方法及***
CN117991179A (zh) * 2024-04-03 2024-05-07 深圳大学 一种基于导向矢量矩阵重构的波达方向估计方法及装置

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117289262A (zh) * 2023-11-27 2023-12-26 中南大学 穿墙雷达目标检测方法及***
CN117289262B (zh) * 2023-11-27 2024-02-06 中南大学 穿墙雷达目标检测方法及***
CN117991179A (zh) * 2024-04-03 2024-05-07 深圳大学 一种基于导向矢量矩阵重构的波达方向估计方法及装置

Similar Documents

Publication Publication Date Title
CN108872926B (zh) 一种基于凸优化的幅相误差校正及doa估计方法
CN107315162B (zh) 基于内插变换和波束形成的远场相干信号doa估计方法
WO2018049595A1 (zh) 一种基于交替方向乘子法的稳健稀疏恢复stap方法及其***
CN110197112B (zh) 一种基于协方差修正的波束域Root-MUSIC方法
CN111046591B (zh) 传感器幅相误差与目标到达角度的联合估计方法
CN112305495B (zh) 一种基于原子范数最小的互质阵列协方差矩阵重构方法
CN107703478B (zh) 基于互相关矩阵的扩展孔径二维doa估计方法
CN113567913B (zh) 基于迭代重加权可降维的二维平面doa估计方法
CN115356678B (zh) 基于dpnalm算法的稀疏阵列doa估计方法
CN111257845A (zh) 一种基于近似消息传递的不在网格目标角度估计方法
CN110673119A (zh) 基于压缩感知的非正则化方位估计方法及***
CN110749855B (zh) 一种基于协方差域零化的均匀线阵波达方向估计方法
CN114814830A (zh) 一种基于鲁棒主成分分析降噪的米波雷达低仰角测高方法
CN115453528A (zh) 基于快速sbl算法实现分段观测isar高分辨成像方法及装置
CN113671485B (zh) 基于admm的米波面阵雷达二维doa估计方法
CN109696651B (zh) 一种基于m估计的低快拍数下波达方向估计方法
CN113625220A (zh) 一种多径信号波达方向和扩散角快速估计新方法
CN110531310B (zh) 基于子空间和内插变换的远场相干信号波达方向估计方法
CN109507634B (zh) 一种任意传感器阵列下的基于传播算子的盲远场信号波达方向估计方法
CN114167346B (zh) 基于协方差矩阵拟合阵元扩展的doa估计方法及***
CN113820653B (zh) 基于动态和差波束的米波雷达低仰角目标doa估计方法
CN114755628A (zh) 非均匀噪声下声矢量传感器阵列波达方向估计方法
CN113093098B (zh) 基于lp范数补偿的轴向不一致矢量水听器阵列测向方法
CN114996653A (zh) 一种基于原子范数最小化的二维鲁棒自适应波束形成方法
CN109633563B (zh) 基于多径信息的自适应相干波束形成方法

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