Loading [MathJax]/jax/output/SVG/jax.js

基于扰动的结合Off-grid目标的层析SAR三维成像方法

杜邦 仇晓兰 张柘 雷斌 丁赤飚

李维新, 李明, 陈洪猛, 等. 针对回波数据异常时的雷达前视超分辨快速成像方法[J]. 雷达学报(中英文), 2024, 13(3): 667–681. doi: 10.12000/JR23209
引用本文: 杜邦, 仇晓兰, 张柘, 等. 基于扰动的结合Off-grid目标的层析SAR三维成像方法[J]. 雷达学报, 2022, 11(1): 62–70. doi: 10.12000/JR21093
LI Weixin, LI Ming, CHEN Hongmeng, et al. Fast radar forward-looking super-resolution imaging for abnormal echo data[J]. Journal of Radars, 2024, 13(3): 667–681. doi: 10.12000/JR23209
Citation: DU Bang, QIU Xiaolan, ZHANG Zhe, et al. L1 minimization with perturbation for off-grid tomographic SAR imaging[J]. Journal of Radars, 2022, 11(1): 62–70. doi: 10.12000/JR21093

基于扰动的结合Off-grid目标的层析SAR三维成像方法

DOI: 10.12000/JR21093
基金项目: 国家自然科学基金(61991421, 61991420)
详细信息
    作者简介:

    杜 邦(1995–),男,中国科学院空天信息创新研究院在读博士,研究方向为阵列信号处理、SAR三维成像

    仇晓兰(1982–),女,中国科学院空天信息创新研究院研究员,博士生导师,IEEE高级会员、IEEE地球科学与遥感快报副主编、雷达学报青年编委。主要研究方向为SAR成像处理、SAR图像理解

    张 柘(1988–),男,博士,中国科学院空天信息研究院、苏州空天信息研究院副研究员,硕士生导师,研究方向为稀疏信号处理与合成孔径雷达成像

    雷 斌(1978–),男,中国科学院空天信息创新研究院研究员,博士生导师,主要研究方向为合成孔径雷达及多源遥感信息处理与应用系统等

    丁赤飚(1969–),男,研究员,博士生导师,中国科学院院士,先后主持多项国家重点项目和国家级遥感卫星地面系统工程建设等项目,曾获国家科技进步奖一等奖、二等奖,国家技术发明二等奖等奖励。研究方向为合成孔径雷达、遥感信息处理和应用系统等领域

    通讯作者:

    仇晓兰 xlqiu@mail.ie.ac.cn

  • 责任主编:廖明生 Corresponding Editor: LIAO Mingsheng
  • 中图分类号: TN957.52

L1 Minimization with Perturbation for Off-grid Tomographic SAR Imaging

Funds: The National Natural Science Foundation of China (61991421, 61991420)
More Information
  • 摘要: 层析合成孔径雷达(TomoSAR)通过组合在不同高度上获取的多基线二维SAR数据,实现合成孔径雷达的三维成像。TomoSAR的求解本质是一维谱估计问题,基于压缩感知的方法可以在非均匀分布的少量基线观测下实现求解,逐渐成为了主流的成像方式。在经典的压缩感知算法流程中,需要将连续的高程向划分成固定的网格,并且假定目标正好位于所划分的网格上。然而该假设通常难以成立,从而引起“基失配”问题,目前该问题在TomoSAR中很少被讨论。该文首先讨论了目标不在网格(Off-grid)上的TomoSAR观测模型,提出了采用加性扰动项来修正目标偏离网格所带来影响的求解模型。在此基础之上,引入局部优化算法与L1范数最小化结合的方法,求解所提出的Off-grid TomoSAR模型。最后,利用仿真数据和机载阵列干涉SAR实际飞行数据进行了验证,结果表明,对于Off-grid目标,该方法能够得到比基于L1范数最小化的经典方法更准确的位置、幅度和相位求解结果,证明了方法的优越性。

     

  • 机载雷达前视成像在态势感知、自主导航、地形跟踪和地形跟随等方面有着广泛的应用[13]。由于扫描雷达具有连续和广域观测的优势,因此常用于前视成像。然而与其他成像技术如双基前视合成孔径雷达成像方法相比,扫描雷达前视成像存在角度分辨率低的问题[47]

    近年来,研究学者提出了多种超分辨方法来获得高角度分辨率[8,9]。在超分辨方法中,扫描雷达回波信号可以看作天线方向图与目标散射系数之间的卷积。目标散射系数可以通过反卷积的方式得到,但是因为天线方向图的低通特性,直接进行反卷积是一个病态问题[10]。在矩阵分析中,天线方向图矩阵中小的奇异值求逆会提高噪声的幅值。根据这个特点,文献[11]提出了截断奇异值分解(Truncated Singular Value Decomposition, TSVD)方法,通过舍弃小的奇异值来抑制噪声。正则化方法也是一个很好地解决病态问题的方法,其核心是选择合适的正则化项[1214]。在文献[15]中,作者提出了Tikhonov正则化,其引入L2范数项来抑制噪声。除此之外,谱估计方法也可以用来提高角度分辨率[1618]。电子科技大学张永超等人[1921]将迭代自适应方法(Iterative Adaptive Approach, IAA)用于前视超分辨成像,并证明了IAA方法可以抑制噪声并提高角度分辨率。与整个成像场景相比,感兴趣成像目标是稀疏的,因此稀疏成像被引入到机载雷达前视超分辨成像中[2225]。在文献[26]中,作者提出了一个贝叶斯稀疏前视成像方法,该方法利用拉普拉斯分布来描述目标的稀疏特性。成像结果显示拉普拉斯分布可以很好地拟合目标,提高前视成像的角度分辨率。当飞机在船舶上自主着陆时,船舶的轮廓信息尤为重要。文献[27]将全变差(Total Variation, TV)和L1范数作为正则化约束项来实现前视成像。该方法在提高雷达前视成像角度分辨率的同时,很好地保留了目标的轮廓信息。

    相比于提高前视成像的角度分辨率,如何降低前视成像的计算复杂度也是当前研究的一个热点[2830]。文献[31]提出将雷达回波信号从时域转换到频域,利用特普利茨矩阵的低位移秩特点,采用Gohberg-Semencul分解实现了矩阵的快速求逆。因为矩阵的维度影响前视成像的实时性,文献[32]采用低秩近似思想对回波信号降维,并采用低秩迭代自适应方法(Low Rank Iterative Adaptive Approach, LRIAA)实现快速前视成像。文献[33]提出了一个二维前视成像处理方法,通过对整体的回波信号进行处理,避免了逐距离单元进行前视成像带来的计算负担,提高了前视成像的实时性。

    在前视成像中,雷达回波可能会受到电磁脉冲干扰,或者受设备性能不稳定影响,导致雷达回波数据出现异常[34,35]。传统的成像方法假设噪声的统计分布服从高斯分布。在面对异常值时,异常值可能会被重构为虚假目标,影响前视成像的准确性。在之前的工作中,我们提出采用学生t分布来拟合雷达回波数据异常时的噪声[35]。相比于高斯混合模型,基于学生t分布的超分辨方法避免了估计过多的系统参数。但是在估计过程中基于学生t分布的方法依然面临高计算复杂度问题。

    本文提出了一种机载雷达超分辨方法实现回波数据异常时的快速前视成像。首先,为了更好地拟合噪声和异常值,我们引入对异常值更加鲁棒的学生t分布。为了更好地重构目标,引入拉普拉斯分层先验来模拟目标的稀疏特征。然后基于期望最大化方法对成像参数进行估计。受TSVD方法的启发,将截断酉矩阵引入目标散射系数的估计公式中。通过矩阵变换降低了求逆矩阵的尺寸从而降低了参数估计的计算复杂度。仿真结果表明本文提出的快速方法可以通过更短的时间实现回波数据异常情况下的前视超分辨成像。

    图1为扫描雷达前视成像的几何示意图。图中机载平台距离地面高度为H,搭载着雷达以速度v向前飞行,天线以扫描速度ω逆时针扫描成像区域。初始时刻,载机位于点A所在位置。此时雷达与目标P1之间的初始距离记为R0,目标P1的方位角记为θ0,俯仰角记为φ0。经过t时刻,载机飞行至B,此时雷达与目标P1的方位角为θ(t),雷达与目标P1之间的瞬时距离记为R(t)

    图  1  扫描雷达几何示意图
    Figure  1.  Geometric model of scanning radar

    假设雷达发射线性调频信号

    s(τ)=rect(τ/T)exp(jπμτ2)exp(j2πfcτ) (1)

    其中,τ为快时间,T为脉冲持续时间,μ为脉冲调制频率,fc为脉冲载波频率,rect()为窗函数。

    对于目标P1,经过脉冲压缩和运动补偿后的回波信号可写为

    s(τ,t)=σh(tt0)sinc[Br(τ2R0c)]exp(j4πλR(t)) (2)

    其中,σ表示目标P1的散射系数。h()表示天线方向图调制函数。Br为信号带宽。λ为信号波长。式(2)中指数项exp(j4πR(t)/λ)为多普勒频移项。瞬时距离R(t)可写为

    R(t)=R20+(vt)22R0vtcosθ0cosφ0R0vtcosθ0cosφ0 (3)

    在某一距离单元,回波信号可写为

    y=(FH)σ+n (4)

    其中,y=[y1,y2,,yM]T为该距离单元上的回波向量。M为回波信号的长度。σ=[σ1,σ2,,σM]T为目标散射向量。n=[n1,n2,,nM]T为成像噪声。F为多普勒相位矩阵。H为由天线方向图向量h=[h1,h2,,hL]T组成的方向图矩阵。L为天线方向图长度。

    F=[111ej2πfd(θ1)PRFej2πfd(θ2)PRFej2πfd(θM)PRFej2πfd(θ1)M1PRFej2πfd(θ2)M1PRFej2πfd(θM)M1PRF]M×M (5)
    H=[hL/2hL/21h1hL/2+1hL/2h1hLh1hLh2hLhL1hL/2]M×M (6)

    fd(θm)为多普勒中心频率。PRF为脉冲重复频率。为向上取整操作。

    fd(θm)=2vcosθmcosφλ (7)

    φ为此距离单元上的俯仰角。为了获取目标散射向量的幅值和位置信息,本文仅考虑式(4)中的幅值信息,即

    |y|=|(FH)σ|+|n| (8)

    方便起见,将式(8)写为

    s=Hx+n (9)

    其中,s=[s1,s2,,sM]T为回波幅值向量,x=[x1,x2,,xM]T为目标散射幅值向量。

    本节提出了快速超分辨方法实现回波数据异常时的前视成像。首先,采用学生t分布来拟合回波数据异常时的噪声。接着引入拉普拉斯分层分布模拟目标的稀疏分布。然后采用期望最大化方法实现对成像参数的估计。最后将截断的酉矩阵引入到目标散射系数的估计公式中,并通过矩阵变换对目标估计公式进行重写,降低估计过程中的计算复杂度。

    在前视成像过程中,雷达边发射波束边扫描前方场景,通过接收目标的扫描信号对前视场景进行成像。然而在成像场景中可能存在无意的窄带电磁干扰,如非合作脉冲雷达发射的信号。在这种情况下,回波信号的某些方位角上将会出现异常值。此外,雷达受设备性能的影响,如应答器性能、通信链路性能、信号时间戳抖动/量化等多种因素。在这种情况下,整个回波信号中会在不同的方位向距离向随机出现异常值。

    受异常值影响,成像噪声变为系统固有高斯白噪声和异常值的组合。文献[36]指出,与高斯分布相比,学生t分布可以容忍回波信号中少数异常值的存在

    p(n|μ,η,v)=Γ(M+ν2)Γ(ν2)(ηπν)M2[1+(nμ)Tη(nμ)ν]M+ν2 (10)

    其中,μ为学生t分布的均值,η为学生t分布的精度。ν为自由度,当ν=1,学生t分布变为柯西分布,当ν,学生t分布变为高斯分布。Γ()为Gamma函数。

    图2给出了高斯分布和学生t分布的概率密度图比较。图中高斯分布的均值为0、方差为1.13,学生t分布的均值为0、精度为1、自由度为2。从图中可以看出,虽然高斯分布与学生t分布的峰值近似,但学生t分布拥有更高的拖尾,更高的拖尾意味着学生t分布可以容忍回波信号中少量的异常值存在,对异常值更加鲁棒。为了直观地展示学生t分布对异常值的鲁棒性,我们分别拟合了高斯数据和混合少量异常值的高斯数据。直方图拟合曲线如图3所示,其中图3(a)为高斯分布和学生t分布对高斯数据的直方图拟合曲线,图3(b)为高斯分布和学生t分布对混合少量异常值数据的直方图拟合曲线。相比于图3(a)中的拟合曲线,图3(b)中高斯分布的拟合曲线受异常值的影响变化较大,学生t分布的拟合曲线变化较小,说明学生t分布比高斯分布对异常值的鲁棒性更好。

    图  2  高斯分布和学生t分布的比较
    Figure  2.  Comparison of Gaussian distribution and student-t distribution
    图  3  直方图拟合曲线
    Figure  3.  Fitting curve of histogram distribution

    同时,学生t分布也可由高斯分布和Gamma分布得到

    p(n|μ,η,v)=0N(n|μ,(ηΣ)1)G(Σ|ν2,ν2)dΣ (11)

    其中,ηΣ=diag(ε1,ε2,,εM)共同表示高斯分布的逆协方差矩阵。N()为高斯分布,G()为Gamma分布。从式(11)可以看出,学生t分布由有限个均值相同、方差不同的高斯分布叠加组成。脉冲压缩后的回波信号已经实现了距离向高分辨,在超分辨过程中通过逐距离单元对方位向回波信号进行处理来提高分辨率。异常值不管出现在某些角度上或随机出现在整个回波信号中,在固定距离单元上异常值的形式是相同的。因此本文引入学生t分布来模拟成像噪声的统计特征。为了简化学生t分布,成像噪声服从

    p(n|η,Σ)=N(n|0,(ηΣ)1) (12)

    在式(12)中,假设成像噪声的均值为零。回波信号s服从

    p(s|x,η,Σ)=N(s|Hx,(ηΣ)1) (13)

    参数η服从Gamma分布

    p(η)=G(η|aη,bη) (14)

    其中,aηbη为超参数,共同控制η的分布。

    同时,参数矩阵Σ服从

    p(Σ|ν)=Mm=1G(εm|ν2,ν2) (15)

    自由度ν服从Gamma分布

    p(ν)=G(ν|aν,bν) (16)

    其中,超参数aνbν共同控制ν的分布。

    前视成像中感兴趣目标相对于成像场景是稀疏的,因此本文引入拉普拉斯分层先验分布模拟目标的稀疏特征

    p(x|λ)=Mm=1λm2exp(λm|xm|) (17)

    然后,λ=[λ1,λ2,,λM]T服从

    p(λ)=Mm=1p(λm)=Mm=1G(λm|aλ,bλ) (18)

    其中,超参数aλbλ共同控制λ的分布。

    上述成像参数的概率模型如图4所示,图中每个结点表示一个变量,链接表示变量之间存在概率关系,箭头表示生成方向,如图中噪声n和目标x共同作用于回波信号s

    图  4  概率图模型
    Figure  4.  Probabilistic graphical model

    3.1节给出了数据异常情况时机载前视成像中噪声和目标的概率模型,本节引入期望最大化方法对成像参数进行估计。期望最大化方法主要分为两步,E步骤和M步骤。E步骤的主要任务是确定目标散射系数xQ函数。M步骤是通过最大化Q函数来估计目标散射系数x

    E步骤:定义Q函数

    Q(ζ;ζ(q))=Eϕ|s,ζ(q)[lnp(s,ζ,ϕ)]=Eη,λ,ν|s,x(q),Σ(q)[lnp(x,η,Σ,λ,ν,s)] (19)

    其中,ϕ表示隐藏变量,ϕ={η,λ,ν}ζ表示模型参数,ζ={x,Σ}Eϕ|s,ζ(q)[]表示关于p(ϕ|s,ζ(q))的期望。ζ(q)表示第q次迭代中模型参数的值。

    M步骤:通过最大化E步骤中Q函数获得模型参数ζ

    ζ(q+1)=arg maxζ[Q( ζ;ζ(q))] (20)

    通过E步骤,可以得到ν的期望为

    ν=2aν+M2bνMm=1(lnεmεm)M (21)

    η的期望为

    η=M+2aη(sHx)TΣ(sHx)+2bη (22)

    λm的期望为

    λm=aλ+1|xm|+bλ (23)

    将式(21)、式(22)和式(23)代入式(19),Q函数可写为

    Q(ζ;ζ(q))=Eϕ|s,ζ(q)[lnp(s,ζ,ϕ)]=12lnΣη2(sHx)TΣ(sHx)+(ν21)lnΣν2ΣΛx1+const (24)

    其中,Λ=diag{λ1,λ2,,λM}

    在M步骤中,通过求解ζQ(ζ;ζ(q))=0找到最优模型参数ζ,可以得到

    εm=1+νη(smHmx)2+ν (25)
    x=η(ηHTΣH+W)1HTΣs (26)

    其中,Hm为天线方向图矩阵的第m行,W=diag{λ1(|x1|2+δ)12,λ2(|x2|2+δ)12,,λM(|xM|2+δ)12}δ为避免L1范数不可微引入的辅助参数。

    在目标散射系数x的估计过程中,需要计算矩阵(ηHTΣH+W)的逆,其计算复杂度为O(M3),本节引入辅助矩阵并通过矩阵变换对目标散射估计公式进行重写,降低了估计过程的计算复杂度。

    方向图矩阵H可以通过奇异值分解写为

    H=UAVT (27)

    其中,U=[u1,u2,,uM]V=[v1,v2,,vM]为酉矩阵,A=diag{α1,α2,,αM}H的奇异值矩阵。

    在TSVD方法中,目标散射系数x的估计方式为

    x=Hr1s=VrA1rUTrs=VrVTrx+ri=1viα1iuHin (28)

    其中,Hr=UrArVTr是截断的方向图矩阵,r为截断参数。Ur=[u1,u2,,ur],Vr=[v1,v2,,vr]为截断的酉矩阵。Ar=diag{α1,α2,,αr}为截断的奇异值矩阵。

    不考虑噪声的情况下,式(28)可以写为

    x=VrVTrx (29)

    从式(29)近似得到

    VrVTrI (30)

    利用式(30)的性质,可以将式(26)重写为

    x=η[ηHTΣH+W]1HTΣs=(HTVrVTrΣ12Σ12VrVTrH+Wη1)1HTVrVTrΣ12Σ12VrVTrs=(HTVrVTrΣ12Σ12VrVTrH+Wη1)1HTVrVTrΣ12(Σ12VrVTrHW1HTVrVTrΣ12+η1I)×(Σ12VrVTrHW1HTVrVTrΣ12+η1I)1Σ12VrVTrs=W1HTVrVTrΣ12(Σ12VrVTrHW1HTVrVTrΣ12+η1I)1Σ12VrVTrs=W1HTVr(VTrHW1HTVr+η1VTrΣ1Vr)1VTrs=W1˜HT(˜HW1˜HT+η1VTrΣ1Vr)1˜s (31)

    其中

    ˜H=VTrH (32)
    ˜s=VTrs (33)

    在原始基于学生t分布的方法中需要计算式(26)矩阵(ηHTΣH+W)的逆,矩阵的尺寸为M×M,其计算复杂度为O(M3)。在经过引入截断酉矩阵并通过矩阵变换后,式(31)中矩阵(˜HW1˜HT+η1VTrΣ1Vr)的尺寸为r×r,计算矩阵逆的计算复杂度变为O(r3)

    截断参数r的选择可以参考TSVD方法的思想。按照表1的雷达系统参数,图5给出了方向图矩阵(尺寸400×400)的奇异值曲线,从图中可以看出大奇异值数量为40左右。因为大的奇异值包含了更多的信息,因此本论文中在选择截断参数r时,可以根据奇异值曲线选择拐角处的截断参数,在尽可能保留大奇异值的同时,平衡加速方法的时间性能。

    表  1  仿真系统参数
    Table  1.  System parameters of simulation
    参数 数值 参数 数值
    扫描速度(/ss) 50 载波频率(GHz) 9.5
    扫描范围() ±10 信号带宽(MHz) 40
    脉冲重复频率(Hz) 1000 平台速度(m/mss) 30
    主瓣波束宽度() 3
    下载: 导出CSV 
    | 显示表格
    图  5  奇异值曲线
    Figure  5.  Singular value curve

    最后,所提方法的成像参数可以通过式(34)得到

    η=M+2aη(sHx)TΣ(sHx)+2bηλm=aλ+1|xm|+bλν=2aν+M2bνMm=1(lnεmεm)Mεm=1+νη(smHmx)2+νx=W1˜HT(˜HW1˜HT+η1VTrΣ1Vr)1˜s (34)

    在初始化参数中,本文将待估参数初始化为xm=HTms/(HTmHm)Σ设置为单位矩阵。超参数如aη, bη, aλ, bλ, aν, bν设置为较小的值,如104。手动选择合适的截断参数r,在保留大奇异值的同时兼顾算法的加速性能。

    本节对基于学生t分布的方法和其加速方法的计算复杂度进行了分析。假设固定距离单元上回波信号s的尺寸为M×1,迭代次数为q。在分析中,因为基于学生t分布的方法和加速方法中均对η, λ, ν, Σ进行了估计,因此本文忽略了这些参数的计算复杂度,仅考虑计算目标散射系数的计算复杂度。

    在基于学生t分布的方法中,计算矩阵HTΣH的复杂度为O(qM3)。计算矩阵(ηHTΣH+W)逆的复杂度为O(qM3)。计算η(ηHTΣH+W)1HTΣs的复杂度为O(qM3)。因此总的计算复杂度为O(3qM3)

    在其加速方法中,需要对方向图矩阵H进行奇异值分解,其计算复杂度为O(M3)。计算˜H的复杂度为O(rM2),计算˜s的复杂度为O(rM),计算矩阵˜HW1˜HT的复杂度为O(qMr2),计算VTrΣ1Vr的复杂度为O(qMr2)。计算矩阵(˜HW1˜HT+η1VTrΣ1Vr)逆的复杂度为O(qr3)。计算W1˜HT(˜HW1˜HT+η1VTrΣ1Vr)1˜s的复杂度为O(qMr2)。加速方法总的计算复杂度为O(3qMr2+qr3+M3+rM2+rM)

    加速方法相比于未加速方法需要进行奇异值分解以及其他附加运算,如计算VTrΣ1Vr, ˜H, ˜s。当截断参数与回波信号长度接近时,加速方法的计算复杂度会比未加速方法的计算复杂度高。根据前面截断参数的选择,rM,因此加速方法的计算复杂度低于基于学生t分布的计算复杂度。

    本节将采用仿真数据和半实测数据验证所提方法的有效性。本文将所提加速的基于学生t分布方法与Tikhonov正则化(Tikhonov regularization, REGU)方法,LRIAA方法,基于学生t分布的方法进行了比较。为方便起见,将基于学生t分布的方法简写为MBSD (Method Based on Student-t Distribution),加速的基于学生t分布的方法简写为AMBSD (Accelerated Method Based on Student-t Distribution)。在仿真中,假设信噪比(Signal-to-Noise Ratio, SNR)为信号与高斯白噪声的比值。仿真中的异常值为生成回波信号后随机加入。

    图6(a)为原始点目标场景,其中两个目标分别位于0.50.5,雷达仿真系统参数如表1所示。

    图  6  点目标仿真结果
    Figure  6.  Point target simulation results

    图6(b)为信噪比为15 dB的雷达回波信号,回波信号的尺寸为400×1。在回波信号中,我们加入了5个幅值随机的异常值。从图中可以看出,经过雷达波束的调制,两个目标被混合在一起无法分辨,并且有异常值融合在了两个目标之间。图6(c)图6(f)分别为不同方法的处理结果。在图6(c)中,REGU方法可以将两个目标大致区分开,但两目标之间的谷值还是较高。而且受异常值影响,很难判断目标的真实位置。图6中LRIAA方法可以很好地区分目标,但是同REGU方法一样,因为没有考虑异常值对前视成像的影响,一部分的异常值被重构为虚假目标。图6(e)图6(f)分别为MBSD方法以及AMBSD方法的处理结果,其中截断参数设为50。从图中可以看出,MBSD方法和AMBSD方法均可以将两相邻目标进行区分并且可以抑制回波中的异常值。

    为了评估所提方法的优越性,本文引入均方误差(Mean Square Error, MSE)来定量地观察不同方法的性能。MSE可以通过式(35)得到

    MSE=1Mc(1Mxˆx22) (35)

    其中,x是真实的目标值,ˆxx的估计值,Mc是蒙特卡罗实验的次数。仿真分别在信噪比为30 dB, 25 dB, 20 dB, 15 dB10 dB的环境下进行500次蒙特卡罗实验。MSE曲线如图7所示,可以观察到随着信噪比的提升,不同方法的性能均得到提升。相对于其他方法,MBSD方法和AMBSD方法可以获得更好的MSE结果。

    图  7  不同信噪比情况下的MSE曲线
    Figure  7.  MSE curves under different SNRs

    然后我们对不同方法的计算时间进行了对比。仿真时的CPU为 Intel Core i5-12500H 3.10 GHz,RAM为16 GB。仿真软件为MATLAB 2021a计算时间对比结果如图8所示。从图中可以看出REGU方法的计算时间相对于其他方法更低,但其角度分辨率较低。LRIAA方法的运行时间低于MBSD方法,但其MSE结果不如MBSD方法。反观AMBSD方法(r=50)在重构目标准确性和运算时间上都拥有一个良好的性能。

    图  8  不同方位维度情况下的运行时间曲线
    Figure  8.  Time curves under different azimuth dimensions

    基于点目标仿真,本节仿真了20 dB信噪比时不同截断参数情况下MSE和计算时间。截断参数对MSE影响的仿真结果如图9所示。从图中可以看出,随着截断参数r取值的逐渐增大,AMBSD方法的MSE逐渐下降并趋于MBSD方法的MSE。并且可以看出MSE曲线的趋势与图5中奇异值曲线的趋势类似,说明大奇异值保留了大部分的原始信息,对目标参数的重构误差影响较大。截断参数对计算时间的影响情况如图10所示。从图中可以看出当截断参数较大时,因为存在附加运算,因此加速方法的运行时间比未加速方法的运行时间长。当截断参数较小时,AMBSD方法的计算时间比MBSD方法的计算时间更短,并且截断参数与方位向回波尺寸的比值越小,AMBSD方法相比于MBSD方法的加速效果更加明显。

    图  9  截断参数对MSE影响
    Figure  9.  Influence of truncated parameter on MSE
    图  10  截断参数对运行时间影响
    Figure  10.  Influence of truncated parameter on running time

    本节我们通过面目标仿真验证所提方法的有效性。图11为原始目标场景,尺寸为300×400(距离向×方位向),图中8个合成目标分别分布于不同位置。AMBSD方法的截断参数设为50。表2给出了面目标仿真的雷达系统参数。受到电磁干扰或雷达设备性能影响,异常值可能会出现在回波中的某个方位角上或者随机出现在整个回波上,因此我们将面目标仿真分为两部分。

    图  11  原始目标场景
    Figure  11.  Original target scene
    表  2  面目标仿真系统参数
    Table  2.  System parameters of area target simulation
    参数 数值 参数 数值
    扫描速度(/ss) 50 载波频率(GHz) 9.5
    扫描范围() ±10 信号带宽(MHz) 40
    脉冲重复频率(Hz) 1000 平台速度(m/mss) 30
    主瓣波束宽度() 5
    下载: 导出CSV 
    | 显示表格

    (1) 受电磁干扰时前视成像:图12(a)为雷达回波信号,此时的信噪比为10 dB,有8个随机幅值的异常值随机分布在某些方位角上。从回波信号中可以看出目标被混合在一起,无法分辨。图12(b)为LRIAA方法的处理结果,因为没有考虑异常值,LRIAA方法受到异常值的影响,2500 m附近的重构目标被干扰严重。在图12(c)图12(d)中,MBSD方法和AMBSD方法可以很好地抑制异常值并且区分目标。为了更好地量化重构结果,表3给出了MSE结果对比和运行时间对比。从表中可以看出,所提AMBSD方法可以用更短的运行时间获得更好的重构结果。

    图  12  受电磁干扰时面目标仿真结果
    Figure  12.  Area target simulation results with electromagnetic interference
    表  3  受电磁干扰时面目标仿真的MSE和运行时间
    Table  3.  MSE and running time of area target simulation with electromagnetic interference
    方法MSE运行时间(s)
    LRIAA方法6.15×10–34.90
    MBSD方法0.65×10–323.14
    AMBSD方法0.70×10–33.82
    下载: 导出CSV 
    | 显示表格

    (2) 设备性能异常时前视成像:从表4的运行时间可以看出,所提AMBSD方法具有更快的运行时间。图13(a)为雷达回波信号,此时的信噪比为10 dB,异常值被随机地添加到回波信号中。图13(b)为LRIAA方法处理结果,图中目标可以被分开,但是异常值依旧零散地分布于整个重构场景中。相比于LRIAA方法的处理结果,MBSD方法和AMBSD方法获得了更好的重构结果。

    表  4  设备性能异常时面目标仿真的MSE和运行时间
    Table  4.  MSE and running time of area target simulation with abnormal equipment performance
    方法MSE运行时间(s)
    LRIAA方法1.09×10–34.80
    MBSD方法0.80×10–323.03
    AMBSD方法0.81×10–33.90
    下载: 导出CSV 
    | 显示表格
    图  13  设备性能异常时面目标仿真结果
    Figure  13.  Area target simulation results with abnormal equipment performance

    前面的仿真实验证明了所提AMBSD方法可以用更快的速度抑制回波信号中的异常值,提高前视成像的角度分辨率。本节采用半实测数据进一步验证所提方法的性能。图14(a)为机载雷达实测回波信号,其中红色矩形框内存在两个混合在一起的目标,回波信号尺寸为513×208(距离向×方位向)。AMBSD方法的截断参数设为50。图14(b)图14(c)为MBSD方法和AMBSD方法的处理结果。从处理结果中可以看出,各方法可以对红色矩形框中的目标进行分辨,提高前视成像的角度分辨率,证明了MBSD方法和AMBSD方法对实测数据的有效性。

    图  14  实测数据结果
    Figure  14.  Processed results of real data

    (1) 受电磁干扰时前视成像:图15(a)为半实测回波信号,我们将异常值随机添加到不同的方位角上。图15(b)图15(d)为不同方法的处理结果。从图15(b)中可以看出,LRIAA方法对异常值的抑制能力较差,并且红色矩形框中的目标没有被区分开。相比于LRIAA方法,MBSD方法和AMBSD方法均可以抑制异常值,并且红色矩形框中的目标也可以被区分出来。表5给出了不同方法的运行时间对比,可以看出AMBSD方法具有更短的运行时间。

    图  15  受电磁干扰时半实测数据仿真结果
    Figure  15.  Processed results of semi-real data with electromagnetic interference
    表  5  半实测数据运行时间(s)
    Table  5.  Running time of semi-real data (s)
    方法 受电磁干扰时运行时间 设备性能异常时运行时间
    LRIAA方法 1.98 1.89
    MBSD方法 4.57 4.58
    AMBSD方法 1.18 1.19
    下载: 导出CSV 
    | 显示表格

    (2) 设备性能异常时前视成像:图16(a)为回波信号,此时异常值随机分布于回波信号中。图16(b)图16(d)为不同方法的处理结果。从图16(c)图16(d)中可以看出MBSD方法和AMBSD方法可以抑制异常值并且区分红色矩形框中的目标。并且表5说明AMBSD方法可以用更短的运行时间获得与MBSD方法近似的分辨性能。相比于表3表4中AMBSD方法与MBSD方法的加速比值,因为在半实测仿真中截断参数与方位向尺寸的比值增大,因此表5中的加速比值更小。

    图  16  设备性能异常时半实测数据仿真结果
    Figure  16.  Processed results of semi-real data with abnormal equipment performance

    本文提出了一种机载雷达超分辨方法实现回波数据异常时快速前视成像。该方法通过利用学生t分布对异常值的鲁棒性,使噪声模型对存在异常值时的回波信号更加灵活。为了更好地重构目标,引入拉普拉斯分层先验来模拟目标的稀疏特征。然后基于期望最大化方法对成像参数进行估计。为了降低参数估计过程的运算量,本文将截断后的酉矩阵引入到目标散射系数的估计公式中,并通过矩阵变换的方法将原方法中需要对M×M维矩阵求逆转变为对r×r维矩阵求逆。从仿真数据和半实测数据中可以看出,本文所提出的快速方法可以用更短的时间实现在回波数据异常情况下的前视超分辨成像。

  • 图  1  TomoSAR成像几何关系

    Figure  1.  The geometry of TomoSAR imaging

    图  2  高程向目标重建结果的简单示意图

    Figure  2.  Sketch map of the reconstruct results in the elevation direction

    图  3  点目标仿真结果

    Figure  3.  Point target simulation results

    图  4  散射点位置估计值与真值距离的统计平均值

    Figure  4.  The average distance between the estimated scattering point position and the true value

    图  5  本文方法优于传统方法的概率与SNR的关系

    Figure  5.  Probability of the proposed algorithm providing better estimation versus SNR

    图  6  获得更精确估计的概率与观测数量的关系

    Figure  6.  Probability of the proposed algorithm providing better estimation versus number of acquisitions

    图  7  RMSE的计算结果

    Figure  7.  RMSE results

    图  8  从谷歌地球获取的实验地区于2015年的光学影像和对应的SAR图像

    Figure  8.  Optical image of the experiment area obtained from Google Earth in 2015 and SAR image of the corresponding area

    图  9  图8(a)中最后一排建筑光学影像

    Figure  9.  Optical image of the last row of buildings Fig. 8(a)

    图  10  图9中红框建筑的三维重建结果

    Figure  10.  TomoSAR results of the building marked in Fig. 9

    图  11  本文方法与传统L1范数最小化的方法对图8中红线标记的方位向切片的成像结果比较

    Figure  11.  Comparison of our method and the conventional process of L1 minimization in the slant range plane for the azimuth-bin marked in Fig. 8

    图  12  图11中区域的放大图

    Figure  12.  Zoom in of the circled area in Fig. 11

    表  1  BPLOT的流程

    Table  1.   The process of BPLOT

     算法 : BPLOT
     (1) 输入: A,A,g,ΔsK
     (2) L1范数求解:
       argmin{||gAgΓ||22+μ||Γ||1}
       Γ=[γγ]t,Ag=[AA]
     (3) For count t=1: K,
     (4) αt=argmaxp{|γp|+|γp|/|2πΔs|}, p(St1)
     (5) St=St1{αt}
     (6) 计算网格偏差 ΔkΔs=Re(γ/(j2πγ))
     (7) 计算散射值
       ˆAf=A(S+ΔkΔs)
       ˆγf=(ˆAfHˆAf)1ˆAfHg
     (8) 输出:ˆγf,S,Δk
    下载: 导出CSV

    表  2  仿真参数

    Table  2.   Simulation parameters

    参数参数值
    载波频率(GHz)14.5
    信号带宽(MHz)500
    通道数N8
    基线宽度(m)0.084
    脉冲重复频率(Hz)480
    平台高度(m)1200
    平台速度(m/s)70
    下视角(°)45
    下载: 导出CSV

    表  3  图3对应的位置计算结果

    Table  3.   Point calculation results corresponding to Fig. 3

    参数真值L1BPLOT
    K=2[0.1506, 0.3764][0.1400, 0.3750][0.1506, 0.3764]
    K=3[0.1501, 0.3764, 0.6120][0.1484, 0.3750, 0.6094][0.1501 0.3765 0.6119]
    下载: 导出CSV
  • [1] FORNARO G, SERAFINO F, and SOLDOVIERI F. Three-dimensional focusing with multipass SAR data[J]. IEEE Transactions on Geoscience and Remote Sensing, 2003, 41(3): 507–517. doi: 10.1109/TGRS.2003.809934
    [2] REIGBER A and MOREIRA A. First demonstration of airborne SAR tomography using multibaseline L-band data[J]. IEEE Transactions on Geoscience and Remote Sensing, 2000, 38(5): 2142–2152. doi: 10.1109/36.868873
    [3] ZHU Xiaoxiang and BAMLER R. Very high resolution spaceborne SAR tomography in urban environment[J]. IEEE Transactions on Geoscience and Remote Sensing, 2010, 48(12): 4296–4308. doi: 10.1109/TGRS.2010.2050487
    [4] FORNARO G, REALE D, and SERAFINO F. Four-dimensional SAR imaging for height estimation and monitoring of single and double scatterers[J]. IEEE Transactions on Geoscience and Remote Sensing, 2009, 47(1): 224–237. doi: 10.1109/TGRS.2008.2000837
    [5] CHAI Huiming, LV Xiaolei, YAO Jingchuan, et al. Off-grid differential tomographic SAR and its application to railway monitoring[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2019, 12(10): 3999–4013. doi: 10.1109/JSTARS.2019.2940730
    [6] AUSTIN C D, ERTIN E, and MOSES R L. Sparse signal methods for 3-D radar imaging[J]. IEEE Journal of Selected Topics in Signal Processing, 2011, 5(3): 408–423. doi: 10.1109/JSTSP.2010.2090128
    [7] FERRARA M, JACKSON J A, and AUSTIN C. Enhancement of multi-pass 3D circular SAR images using sparse reconstruction techniques[C]. Proceedings of SPIE 7337, Algorithms for Synthetic Aperture Radar Imagery XVI, Orlando, United States, 2009. doi: 10.1117/12.820256.
    [8] ZHU Xiaoxiang and BAMLER R. Super-resolution power and robustness of compressive sensing for spectral estimation with application to spaceborne tomographic SAR[J]. IEEE Transactions on Geoscience and Remote Sensing, 2012, 50(1): 247–258. doi: 10.1109/TGRS.2011.2160183
    [9] ZHU Xiaoxiang and BAMLER R. Tomographic SAR inversion by L1-norm regularization—The compressive sensing approach[J]. IEEE transactions on Geoscience and Remote Sensing, 2010, 48(10): 3839–3846. doi: 10.1109/TGRS.2010.2048117
    [10] 李杭, 梁兴东, 张福博, 等. 一种基于组稀疏的阵列干涉SAR三维重建方法[J]. 中国科学:信息科学, 2018, 48(8): 1051–1064. doi: 10.1360/N112017-00023

    LI Hang, LIANG Xingdong, ZHANG Fubo, et al. A novel 3-D reconstruction approach based on group sparsity of array InSAR[J]. Science in China Information, 2018, 48(8): 1051–1064. doi: 10.1360/N112017-00023
    [11] JIAO Zekun, DING Chibiao, QIU Xiaolan, et al. Urban 3D imaging using airborne TomoSAR: Contextual information-based approach in the statistical way[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2020, 170: 127–141. doi: 10.1016/j.isprsjprs.2020.10.013
    [12] MA Peifeng, LIN Hui, LAN Hengxing, et al. On the performance of reweighted L1 minimization for tomographic SAR imaging[J]. IEEE Geoscience and Remote Sensing Letters, 2015, 12(4): 895–899. doi: 10.1109/LGRS.2014.2365613
    [13] 廖明生, 魏恋欢, 汪紫芸, 等. 压缩感知在城区高分辨率SAR层析成像中的应用[J]. 雷达学报, 2015, 4(2): 123–129. doi: 10.12000/JR15031

    LIAO Mingsheng, WEI Lianhuan, WANG Ziyun, et al. Compressive sensing in high-resolution 3D SAR tomography of urban scenarios[J]. Journal of Radars, 2015, 4(2): 123–129. doi: 10.12000/JR15031
    [14] 王爱春, 向茂生. 基于块压缩感知的SAR层析成像方法[J]. 雷达学报, 2016, 5(1): 57–64. doi: 10.12000/JR16006

    WANG Aichun and XIANG Maosheng. SAR tomography based on block compressive sensing[J]. Journal of Radars, 2016, 5(1): 57–64. doi: 10.12000/JR16006
    [15] CHI Yuejie, SCHARF L L, PEZESHKI A, et al. Sensitivity to basis mismatch in compressed sensing[J]. IEEE Transactions on Signal Processing, 2011, 59(5): 2182–2195. doi: 10.1109/TSP.2011.2112650
    [16] HERMAN M A and NEEDELL D. Mixed operators in compressed sensing[C]. The 44th Annual Conference on Information Sciences and Systems (CISS), Princeton, USA, 2010. doi: 10.1109/CISS.2010.5464909.
    [17] CANDÈS E and ROMBERG J. Sparsity and incoherence in compressive sampling[J]. Inverse Problems, 2007, 23(3): 969–985. doi: 10.1088/0266-5611/23/3/008
    [18] STOICA P and BABU P. Sparse estimation of spectral lines: Grid selection problems and their solutions[J]. IEEE Transactions on Signal Processing, 2012, 60(2): 962–967. doi: 10.1109/TSP.2011.2175222
    [19] TANG Gongguo, BHASKAR B N, SHAH P, et al. Compressed sensing off the grid[J]. IEEE Transactions on Information Theory, 2013, 59(11): 7465–7490. doi: 10.1109/TIT.2013.2277451
    [20] YANG Zai and XIE Lihua. Enhancing sparsity and resolution via reweighted atomic norm minimization[J]. IEEE Transactions on Signal Processing, 2016, 64(4): 995–1006. doi: 10.1109/TSP.2015.2493987
    [21] 陈栩杉, 张雄伟, 杨吉斌, 等. 如何解决基不匹配问题: 从原子范数到无网格压缩感知[J]. 自动化学报, 2016, 42(3): 335–346. doi: 10.16383/j.aas.2016.c150539

    CHEN Xushan, ZHANG Xiongwei, YANG Jibin, et al. How to overcome basis mismatch: from atomic norm to gridless compressive sensing[J]. Acta Automatica Sinica, 2016, 42(3): 335–346. doi: 10.16383/j.aas.2016.c150539
    [22] FANNJIANG A and LIAO Wenjing. Super-resolution by compressive sensing algorithms[C]. 2012 Conference Record of the Forty Sixth Asilomar Conference on Signals, Systems and Computers (ASILOMAR), Pacific Grove, USA, 2012. doi: 10.1109/ACSSC.2012.6489036.
    [23] HERMAN M A and STROHMER T. General deviants: An analysis of perturbations in compressed sensing[J]. IEEE Journal of Selected Topics in Signal Processing, 2010, 4(2): 342–349. doi: 10.1109/JSTSP.2009.2039170
    [24] YANG Zai, XIE Lihua, and ZHANG Cishen. Off-grid direction of arrival estimation using sparse Bayesian inference[J]. IEEE Transactions on Signal Processing, 2013, 61(1): 38–43. doi: 10.1109/TSP.2012.2222378
    [25] FANNJIANG A and TSENG H C. Compressive radar with off-grid targets: A perturbation approach[J]. Inverse Problems, 2013, 29(5): 054008. doi: 10.1088/0266-5611/29/5/054008
    [26] CHEN S S, DONOHO D L, and SAUNDERS M A. Atomic decomposition by basis pursuit[J]. SIAM Review, 2001, 43(1): 129–159. doi: 10.1137/S003614450037906X
    [27] 仇晓兰, 焦泽坤, 彭凌霄, 等. SARMV3D-1.0: SAR微波视觉三维成像数据集[J]. 雷达学报, 2021, 10(4): 485–498. doi: 10.12000/JR21112

    QIU Xiaolan, JIAO Zekun, PENG Lingxiao, et al. SARMV3D-1.0: Synthetic aperture radar microwave vision 3D imaging dataset[J]. Journal of Radars, 2021, 10(4): 485–498. doi: 10.12000/JR21112
    [28] FANNJIANG A and LIAO Wenjing. Coherence pattern-guided compressive sensing with unresolved grids[J]. SIAM Journal on Imaging Sciences, 2012, 5(1): 179–202. doi: 10.1137/110838509
    [29] SCHMITT M. Reconstruction of urban surface models from multi-aspect and multi-baseline interferometric SAR[D]. [Ph. D. dissertation], Technische Universität München, 2014.
    [30] ZHU Xiaoxiang. Spectral estimation for synthetic aperture radar tomography[D]. [Master dissertation], Technische Universität München, 2008.
    [31] 李杭. 阵列干涉SAR高精度三维成像算法研究[D]. [博士论文], 中国科学院大学, 2017.

    LI Hang. Research on high-precision 3-D imaging algorithm of array interferometric synthetic aperture radar[D]. [Ph. D. dissertation], University of Chinese Academy of Sciences, 2017.
  • 期刊类型引用(1)

    1. 毛德庆,杨建宇,杨明杰,张永超,张寅,黄钰林. IAA-Net:一种实孔径扫描雷达迭代自适应角超分辨成像方法. 雷达学报. 2024(05): 1073-1091 . 本站查看

    其他类型引用(0)

  • 加载中
图(12) / 表(3)
计量
  • 文章访问数: 1666
  • HTML全文浏览量: 835
  • PDF下载量: 265
  • 被引次数: 1
出版历程
  • 收稿日期:  2021-07-01
  • 修回日期:  2021-09-17
  • 网络出版日期:  2021-10-15
  • 刊出日期:  2022-02-28

目录

/

返回文章
返回