Loading [MathJax]/jax/element/mml/optable/MathOperators.js

基于对数域加性信号分解的时序SAR图像相干斑抑制方法

康健 童风雨 白雨松 丁翔 冀腾宇 张柘

杨悦, 叶星炜, 张方正, 等. 基于微波光子I/Q去斜接收的宽带线性调频雷达成像系统[J]. 雷达学报, 2019, 8(2): 224–231. doi: 10.12000/JR19002
引用本文: 康健, 童风雨, 白雨松, 等. 基于对数域加性信号分解的时序SAR图像相干斑抑制方法[J]. 雷达学报, 2023, 12(5): 1031–1043. doi: 10.12000/JR22242
YANG Yue, YE Xingwei, ZHANG Fangzheng, et al. Broadband LFM radar imaging system based on microwave photonic I/Q de-chirping[J]. Journal of Radars, 2019, 8(2): 224–231. doi: 10.12000/JR19002
Citation: KANG Jian, TONG Fengyu, BAI Yusong, et al. SAR time series despeckling based on additive signal component decomposition in logarithm domain[J]. Journal of Radars, 2023, 12(5): 1031–1043. doi: 10.12000/JR22242

基于对数域加性信号分解的时序SAR图像相干斑抑制方法

DOI: 10.12000/JR22242
基金项目: 国家自然科学基金(62101371),江苏省青年基金(BK20210707)
详细信息
    作者简介:

    康 健,副教授,硕士生导师,主要研究方向为SAR/InSAR信号处理、遥感图像解译等

    童风雨,硕士生,主要研究方向为目标检测技术

    白雨松,硕士生,主要研究方向为SAR干涉技术

    丁 翔,硕士生,主要研究方向为SAR图像去噪

    冀腾宇,讲师,硕士生导师,主要研究方向为张量分解等

    张 柘,研究员,博士生导师,主要研究方向为稀疏信号处理、稀疏微波成像等

    通讯作者:

    康健 jiankang@suda.edu.cn

    张柘 zhangzhe01@aircas.ac.cn

  • 责任主编:陈溅来 Corresponding Editor: CHEN Jianlai
  • 中图分类号: TP753

SAR Time Series Despeckling Based on Additive Signal Component Decomposition in Logarithm Domain

Funds: The National Natural Science Foundation of China (62101371), Jiangsu Province Science Foundation for Youths (BK20210707)
More Information
  • 摘要: 随着合成孔径雷达(SAR)在测绘带宽度、空间以及时间分辨率上的大幅提升,由不同时间获取的SAR图像配准得到的时间序列能更加精确地提供观测区域的动态变化信息。然而,相干斑噪声以及沿时间维度突变信号为后续的解译工作带来了严重挑战。尽管现有的主流方法可以对时序SAR图像的相干斑进行有效抑制,但沿时间维度突变信号会对去噪结果产生干扰。为更好地解决此问题,该文提出了一种基于对数域加性信号分解的方法,能同时抑制相干斑噪声并且对时序图像中的稳定信号和沿时间维度突变信号进行分离,从而消除突变信号对于去噪结果的影响。在仿真数据受到突变信号干扰的情况下,所提方法与其他主流方法相比,其去噪结果在峰值信噪比(PSNR)指标上取得了大约3 dB的提升。在哨兵1号数据上,所提方法能鲁棒地对时序图像中的相干斑噪声进行抑制,并且得到的突变信号成分也为后续的解译工作提供了参考数据。

     

  • 飞机尾流是飞机飞行时所产生的固有的流体力学现象,与周围大气相比,尾流具有独特的空气动力学特性和雷达特性,具有较强的雷达可探测性。尾流的雷达探测是航空安全、反隐身等领域的前沿研究课题。处于潮湿大气中的尾流具有较强的散射截面积(Radar Cross Section, RCS),但受到气象杂波的影响,需要根据功率谱特性来提高可探测性。飞机尾流回波经过傅里叶变换后,其功率谱是展宽的,其检测方法包括动目标检测(Moving Target Detection, MTD)方法、谱宽检测方法、神经网络方法等,这些方法有较多局限性。在传统的MTD技术中,首先估计观测数据的功率谱,然后对各频率通道进行恒虚警率(Constant False Alarm Rate, CFAR)检测,较难对展宽的尾流功率谱进行有效积累;而谱宽检测方法在云杂波环境下虚警较高;神经网络方法需要大量的训练样本。

    信息几何是在黎曼流形上采用微分几何方法研究信息领域和统计学问题的一门学科,其主要思想是把参数化的概率分布族构建为统计流形,并采用Fisher信息矩阵作为流形上的黎曼度量,以此为基础得到统计流形上的联络、曲率、测地线距离等概念,从而采用微分几何方法解决统计学中的问题。其研究对象是由概率分布族所构成的统计流形,矩阵流形就是一种重要的统计流形。与利用回波的功率谱相比,利用回波的协方差矩阵描述回波特性,其优点有:(1)不对回波进行多普勒处理,避免了多普勒处理所面临的功率谱展宽问题;(2)协方差矩阵所构成的空间在数学上为负曲率空间,适合用信息几何方法研究矩阵流形的几何结构,有利于更准确度量目标、噪声或杂波协方差矩阵间的差别。

    从信息几何的角度思考问题,能为许多问题提供更加符合实际的模型框架[1],并得出新的解决途径。1945年,印度统计学家Rao提出用Fisher信息矩阵来定义统计流形上的黎曼度量,Fisher度量与信息理论的相对熵具有紧密联系,都可以定量描述两个概率分布间可以区分的程度,从而开启了统计的几何学理论研究[2]。黎曼度量描述的是流形上切空间的性质,联络把流形上不同点处的切空间联系起来,从而使流形具有整体性。20世纪70年代,苏联数学家Chentsov引入了一族单参数的仿射联络,并证明Fisher信息与仿射联络在概率分布流形中的唯一性[3],随后,Amari用α-仿射联络对尽可能多的统计问题进行了几何解释,取得了很大的成功,并建立了统计流形的对偶几何结构来研究信息处理,从而丰富了传统的微分几何理论,由于能够较好地描述概率分布的内在特性,微分流形的理论和方法逐步进入统计和信息处理领域[4]。有了联络,就可以定义相对于该联络的曲率,曲率描述了流形中给定区域的体积与欧氏空间中与其相当区域体积的差异程度,曲率常用于研究Jacobi场的稳定性。测地线是欧氏空间上直线段在黎曼流形上的推广,测地线是连接流形上两点的最短距离的曲线,可用于计算多元概率分布间或协方差矩阵间的测地线距离,与Kulback-Leibler分离度相比,测地线距离的优势在于它满足距离的对称性和三角不等式等重要性质[5]。实际上,测量通常含有噪声或杂波,我们需要利用统计学方法减少测量的不确定性,为了求解噪声或杂波协方差矩阵的均值,需要使目标函数,即协方差矩阵均值的经验方差,取得最小值,文献[6]利用自然梯度下降算法给出了黎曼均值的迭代计算方法,与算术均值相比,黎曼均值具有更小的经验方差。统计流形上的曲率、联络、测地线和均值等概念大大扩展了统计流形的应用范围,促进了信息几何学的产生与发展。

    经过20余年的发展,信息几何学的理论迅速扩展,并在多个领域获得成功应用,已经开始应用到统计推断、控制理论、神经网络、医学成像、计算视觉、信号处理等领域[1,7]。由于较多的实际问题可以在矩阵流形上进行研究,如雷达信号处理[8,9]、矩阵方程求解[1]、医学图像处理[5,10]、流形学习、系统的稳定性与最优化等问题,将信息几何应用于矩阵流形上促进了矩阵信息几何的诞生。在用信息几何理论解决信息领域或统计学问题时,首先需要把所研究的问题构建成统计流形或矩阵流形,采用微分几何或李群方法研究流形的几何结构,从而用几何方法解决涉及的非线性问题。孙华飞等人将矩阵信息几何应用于求解矩阵方程问题,大多数矩阵方程难以得到解析解,通过将问题转化为矩阵流形上的优化问题,可以用自然梯度给出求解方案[1]。在雷达信号处理、目标检测和数据处理领域,法国Thales Air Systems研究员Barbaresco等人利用信息几何方法研究了近海小目标检测[8]、飞机尾流的检测和成像[8]、极化数据处理、STAP处理[11]等问题。文献[12]在信息几何框架下研究了雷达系统的信息分辨、信号检测、参数估计、目标跟踪等基础性和科学性问题,为雷达信号处理提供了一套全新的分析方法。总之,对于信号处理中涉及的非欧空间、非线性处理等问题,信息几何方法有优势。本文在前人研究基础上介绍了矩阵流形的几何结构,然后针对飞机尾流目标检测问题,提出了一种矩阵CFAR检测方法,分析了尾流目标的检测性能。

    飞机尾流是飞机产生升力的必然产物,主要包括尾喷流和翼尖尾流两类。翼尖尾流包括两个旋向相反的漩涡,其持续时间长达数十秒到几分钟,在空间上表现为飞机后一个长的柱状分布介质目标,一般到达飞机后部百倍翼展,长达数千米,属于一种较强的湍流。欧、美、俄等对飞机尾流的关注是从航空安全领域开始的。当飞机进入前机的尾流时,会因受力不均而出现机身抖动、下沉、翻转等现象,严重时会导致飞机坠毁。飞机尾流对民用航空安全构成了巨大威胁,近年来,欧盟和美国分别投入巨资启动了SESAR (2005-2020)计划和NGATS (2005-2025)计划以提高航空安全性能[13],飞机尾流的实时探测技术无一例外地成为各个研究计划的主要内容。从80年代开始至今,英国、法国和美国进行了多次尾流探测实验测量尾流的雷达反射率和多普勒特性,理论研究和实验结果表明,雷达能够探测到晴空尾流回波,并且能够在雨、雾等潮湿天气中有效工作,因而有望取代激光雷达成为重要的尾流探测手段。国防科技大学是国内较早从事飞机尾流探测技术研究的单位,已有周彬、李健兵、李军、刘俊凯、屈龙海、刘忠训、牛凤梁等多位博士研究尾流特性及其探测技术。

    处于云、雨、雾等潮湿大气中的飞机尾流具有较强的雷达散射截面,但受到气象杂波的影响,需要根据尾流功率谱特性来提高脉冲多普勒雷达对尾流的可探测性。在潮湿大气环境下,尾流的雷达后向散射主要是由大气中的水汽微粒产生的,水气微粒的运动规律可用尾流涡旋速度模型模拟[14],这里使用一群水汽微粒模拟尾流涡旋的运动特征,使用视频回波仿真方法对尾流回波进行建模仿真。

    首先,使用一群粒子表示飞机尾流内的水汽微粒,在时刻t把所有粒子的雷达回波矢量求和,得到一个时刻尾流的合成回波矢量。t时刻离散视频回波为:

    S(t)=Mm=1Amejφm(t) (1)

    其中,Am是第m个粒子的回波幅度,回波相位决定于粒子与雷达的距离,第m个粒子回波的相位与雷达距离的关系可表示为:

    φm(t)=4πRm(t)/λw (2)

    其中, Rm(t) 为第m个粒子t时刻距雷达的距离, λw 是雷达波长。

    接着,假设每个粒子以Burnham-Hallock模型描述的速度做匀速圆周运动,t+1时刻运动到新的位置后,得到新的合成回波矢量 S(t+1) 。其中,每个粒子运动的Burnham-Hallock速度模型表达式为[15]

    v(rm)=Γ02πrmrm2+r2c (3)

    其中, rm 表示尾涡上第m个粒子距尾涡中心的距离, rc 是尾涡的核半径, Γ0 是涡环量,表示点涡的强度,表达式为:

    Γ0=4MagπρVaBa (4)

    式中, Ma 为飞机重量,g为重力加速度, ρ 为室温空气密度(1.16 kg/m3), Va 为飞机的速度, Ba 为翼展。

    最后,得到一个相关处理间隔内若干个时刻的回波实部和虚部,对这一复回波序列进行傅里叶变换,得到功率谱。以波音747飞机为例进行仿真。波音747的典型参数为:飞机重量 Ma=350,000kg ,飞机的速度 Va=200m/s ,翼展 Ba=68m ,重力加速度 g=9.81m/s2 ,室温空气密度 ρ=1.16kg/m3 ,计算得到 Γ0=277 m2/s。一次仿真得到的飞机尾流回波如图1所示。进行多次仿真,可得到相似的归一化径向速度谱。

    图  1  仿真的飞机尾流雷达回波
    Figure  1.  Simulation of the aircraft wake radar echo

    1992年至1993年,美国通用电气公司马可尼研究中心的Sheppard使用X波段全相参脉冲多普勒雷达DX04进行了尾流雷达探测实验。在距离4.5 km处,经过相参处理,波音747飞机尾流回波接近接收机噪声电平,超过噪声尖峰2~3 dB,分析得到尾流RCS约为–73 dBsm, X波段雷达测量的功率谱如图2所示[16]。可见,视频回波建模方法仿真结果与马可尼研究中心测量的X波段功率谱具有相似的频域展宽特性,多普勒速度与尾涡旋转径向速度一致。

    图  2  马可尼研究中心X波段雷达测量的功率谱[16]
    Figure  2.  The Power Spectrum measured by the X-band radar in Marconi Research Center

    对于尾流目标来说,功率谱是展宽的且具有两个较大的峰值,在不同机型和不同飞行速度时,功率谱的峰值位置也不同,较难设计有效的频域CFAR检测器,特别是,需防止展宽的功率谱处于参考单元中,导致检测门限增大和检测概率降低。

    设观测样本 x=[x(0)x(1)···x(n1)]T ,传统的统计检测方法首先对观测样本提出假设:

    {H0:x=vH1:x=s+v (5)

    式中,s为目标信号向量, v 为噪声向量;然后从观测样本出发,制定检测准则,如奈曼-皮尔逊准则,通常可以得到似然比检测判决表达式:

    λ(x)=p(x|H1)p(x|H0)>η (6)

    若式(6)成立,则接受H1,否则接受H0,式中, η 为判决门限,由给定的先验条件确定。对左边的似然比进行化简,化为x的尽可能简单的表达式,称为检测统计量。

    假设信号和噪声都可看做零均值的高斯过程,则观测样本x服从零均值复多元高斯分布,即

    p(x|R)=1πn|R|exHR1x (7)

    其中,Rn阶协方差矩阵。从信息几何的角度看,零均值多元复高斯分布族可以构成一个统计流形[4]

    S={p(x|R)|RΘ} (8)

    式中, ΘCn×n n阶共轭对称正定矩阵空间。R为统计流形S的坐标,S又称为矩阵流形,其几何结构反应了概率密度函数族内在的本质属性。此时,H0H1条件下的似然函数 p(x|H0) p(x|H1) 分别对应统计流形S上的两点,将其坐标记为协方差矩阵R0R1

    在统计流形上以H0H1条件下似然函数间的距离,或者说,以统计流形S上两点R0R1间的距离 D(R1,R0) ,作为检测统计量,可以得到一类距离检测器,若

    D(R1,R0)>η (9)

    则接受H1,否则接受H0,式中, η 为根据虚警概率计算得到的判决门限。在检测器中,R0R1分别表示噪声或杂波协方差矩阵和检测单元观测样本协方差矩阵。

    由式(9),可得虚警概率和检测概率表达式

    Pfa=Pr[D(R1,R0)η|H0] (10)
    Pd=Pr[D(R1,R0)>η|H1] (11)

    实际中,为了满足检测器的恒虚警率特性,需要实时估计R0的值,R0的估计性能对检测性能具有重要影响。

    在信息论中应用较广的距离度量是Kulback-Leibler分离度,其定义为:

    KLD(R1,R0)=KLD[p(x|R1),p(x|R0)]=p(x|R1)lnp(x|R1)p(x|R0)dx

    Kulback-Leibler分离度容易计算,但也存在一定的局限性,即对于流形上相距较近的两点,它与测地线距离等价,但当两点相距较远时,由于未考虑流形的结构和连接两点的最短路径,它与测地线距离差异较大。另外,Kulback-Leibler分离度不满足距离定义的对称性、三角不等式等性质,并不是一个严格意义上的距离度量。程永强等人针对确定性信号检测和随机信号检测问题,将经典的似然比检测方法等效为统计流形上的“广义最小距离检测器”,并在信息几何框架下研究了广义似然比检测器和局部最大势检测器,从而将检测问题转变为统计流形上的几何问题来研究[12]。赵兴刚等人研究了低信杂比条件下基于Kulback-Leibler分离度的距离检测器和基于测地线距离的距离检测器,当信杂比较低时,统计流形上两假设分布间的距离很小,此时可以使用距离的近似形式来定义距离检测器,以使最终的检测器形式上更简单,更易实现[17]。这些检测器并不是严格意义上的恒虚警率(CFAR)检测器。

    在雷达目标检测中,要求检测器能实时估计检测单元的杂波功率、杂波概率分布类型或其参数,以满足恒虚警率特性,为此,前人提出了单元平均CFAR、有序统计类CFAR、自适应CFAR、双参数CFAR等检测器。在信息几何框架下,将观测数据协方差矩阵构成一个统计流形,它成为一个矩阵流形,类比CFAR检测的思想,首先计算参考单元中协方差矩阵的均值,以该均值作为检测单元中噪声或杂波协方差矩阵的估计值,然后利用检测单元协方差矩阵与参考单元协方差矩阵均值间定义的距离作为检测统计量,所获得的检测器称为矩阵CFAR检测器。文献[17]在矩阵流形上,研究了黎曼度量、测地线、测地线距离,以及矩阵黎曼均值的迭代计算方法,在此基础上对矩阵CFAR进行了讨论。这里仅给出主要概念和结论。

    在统计流形S上,研究流形上的几何性质。在统计流形上,Fisher信息矩阵是唯一的度量不变量[2],通常将Fisher信息矩阵作为黎曼度量张量,直接计算或通过仿射黎曼度量计算,得到黎曼度量分量为[1]

    gij=E[ilnp(x|R)jlnp(x|R)]=iR,jRR=tr(R1iRR1jR),i,j=1,2,···,n (12)

    式中,tr表示矩阵的迹。以R R+dR 为参数的相邻两个点 p(x|R) p(x|R+dR) 之间的微分距离 ds2 可以通过Fisher信息矩阵得到[9]

    ds2=ni=1nj=1gijdθidθj=tr[R1(ni=1iRdθi)R1(ni=1jRdθj)]=tr[(R1dR)2]=R1dR2F=R1/2dRR1/22F (13)

    其中, dR=ni=1iRdθi , F 为矩阵Frobenius范数。对于任意的 M=[mij]Cn×n ,则 M2F= ni=1nj=1|mij|2

    流形上两点间的距离与连接两点的曲线有关,在所有可选的曲线中,使距离最短的曲线称为测地线,相应的最小距离称为测地线距离。测地线距离是欧氏空间中直线距离在流形上的推广。对于矩阵流形上的任意点 AΘ ,设A点的切向量为 vTAΘ ( TAΘ A点的切空间)。利用指数映射的概念,则流形上以A为起点,沿 v 方向随参数t变化的测地线[18]为:

    γ(t)=A1/2etA1/2vA1/2A1/2,  0t1 (14)

    γ(0)=A ,令 γ(1)=B ,可得连接两点AB的测地线为:

    γ(t)=A1/2etlog(A1/2BA1/2)A1/2=A1/2(A1/2BA1/2)tA1/2,0t1 (15)

    沿测地线 γ(t) ds2 进行积分,可得AB之间的测地线距离[8]为:

    D(A,B)=minBAds2(R)=γ(t)[ds(R(t))dt]2dt=log(A1/2BA1/2)F=nk=1ln2(λk) (16)

    式中, R(t) 表示测地线 γ(t) 上对应点的矩阵, λk(k=1,···,n) 为对 A1/2BA1/2 进行谱分解得到的第k个特征值,利用矩阵的谱分解能够简化矩阵运算。与Kulback-Leibler分离度相比,测地线距离具有以下性质[5]

    ①非负性: D(A,B)0 , D(A,B)=0 A=B

    ②对称性: D(A,B)=D(B,A)

    ③三角不等式: D(A,C)D(A,B)+D(B,C)

    ④矩阵合同变换条件下不变性:D (A, B)= D(PAPH,PBPH) , Pn阶可逆矩阵;

    ⑤矩阵求逆条件下不变性: D(A,B)= D(A1,B1)

    这些性质使得测地线距离成为性能较好的检测统计量。

    对于已知的N个矩阵 Rk(k=1,2,···,N) ,其黎曼均值 ˉR 是使得目标函数,即下式中均方测地线距离表示的经验方差

    J(R)=1NNk=1D2(Rk,R) (17)

    取最小值时所对应的矩阵,即

    ˉR=argminRΘJ(R) (18)

    与算术均值相比,矩阵流形上的黎曼均值符合流形的几何性质,具有更高的估计精度。通常,利用目标函数 J(R) R的梯度计算局部最小值,该梯度

    J(R)=1NNk=1D2(Rk,R)=1NNk=1log(R1/2kRR1/2k)2F=RNNk=1log(R1/2kRR1/2k) (19)

    J(R) 获得局部最小值的充要条件为 J(R)=0 。当 N>2 时较难获得解析表达式,梯度下降算法是一种可行的计算方法。

    梯度下降方法的实现思想是:任意指定流形上的一个协方差矩阵作为黎曼均值的假定值,则该假定值与真实矩阵均值 ˉR 之间就确定了一条测地线 γˉR(t) , t0 ,该测地线上任一参数t对应的协方差矩阵记为 ˉR(t) ,则均值的假定值对应了 ˉR(0) , ˉR 对应了 ˉR(t) 在终点处的取值。利用式(14)所给的测地线方程,以 ˉR(t0) , t0=0 为初始值,切向量 ν 为目标函数 J(R) ˉR(t0) 处的负梯度 ν=J[ˉR(t0)] ,令 τ 为增量, τ<1 ,可得测地线 γˉR(t) t0+τ 处对应的点 ˉR(t0+τ) ,经多次迭代计算,直到目标函数 J(R) 小于一个较小值。

    文献[5]和文献[6]利用梯度下降算法给出了 ˉR 的迭代计算方法,将实矩阵情况推广到复矩阵情况下,令第i次迭代计算得到的黎曼均值估计值为 ˉRi ,第 i+1 次迭代计算得到的黎曼均值估计值为:

    ˉRi+1=ˉR1/2ieτNNk=1log(ˉR1/2iRkˉR1/2i)ˉR1/2i,0τ1   (20)

    选择迭代的初始值 ˉR0 和步长 τ ,就可以对矩阵的黎曼均值进行迭代运算。

    基于黎曼均值和测地线距离,矩阵CFAR检测器的框图如图3所示,其中,观测数据为根据复采样序列估计得到的样本协方差矩阵,当矩阵维数较大时,需要降低协方差矩阵的维数以提高矩阵处理的效率,滑窗中参考单元的长度为N,滑窗中的阴影部分为保护单元,RD为检测单元的协方差矩阵, ˉR 是参考单元中的N个协方差矩阵的黎曼均值,T是根据虚警概率计算得到的检测门限。

    图  3  基于信息几何的矩阵CFAR检测器框图
    Figure  3.  Block diagram of CA-CFAR detector based on information geometry

    检测过程中,首先利用式(20)计算参考单元的黎曼均值 ˉR ,用于估计背景噪声或杂波协方差矩阵;然后利用式(16)计算RD ˉR 之间的测地线距离平方作为检测统计量进行检测:

    D2(RD,ˉR) (21)

    由于虚警概率解析表达式难以得到,不易从理论上计算检测门限T,可利用Monte-Carlo仿真得到T。说明,为避免 D\left( {{{R}_{\rm{D}}},\bar{{R}}} \right) 中的根号运算和称呼简便,后续称 {D^2}\left( {{{R}_{\rm{D}}},\bar{{R}}} \right) 为测地线距离。

    矩阵CFAR检测方法直接利用观测样本序列的协方差矩阵描述回波特性,并进行目标检测,协方差矩阵中包括了目标及杂波时间相关性的原始信息,有利于检测飞机尾流目标这类的多普勒展宽目标,能避免常规的多普勒处理所面临的问题。由于利用矩阵流形的结构计算参考单元的协方差矩阵均值,即满足式(18)给出的条件,该均值具有较小的经验方差,从而使得CFAR检测器具有更稳健的恒虚警特性。

    需要说明的是,上面给出了一般的矩阵流形几何结构,当使用不同的统计流形的坐标、距离定义、均值定义和梯度下降算法等,能够获得不同的几何结构,同时也可以利用更多的几何结构,或采用更合适的检测策略或检测统计量,以便适应于不同的检测场景。

    当协方差矩阵的维数较大时,矩阵CFAR检测器的运算效率较低,不利于实时处理,因此有必要在检测性能下降较少的条件下,尽量降低协方差矩阵的维数。在对尾流回波协方差矩阵进行降维处理后,一种直观的评价标准是,降低维数之前与之后,根据相关系数估计得到的功率谱基本一致。

    由于矩阵CFAR检测方法利用的是观测数据的协方差矩阵,而常规MTD检测方法利用的是观测数据的功率谱,因此这里分别给出观测数据的自相关函数和功率谱,以对两种检测方法进行对比。自相关函数和功率谱是描述脉冲序列时间相关性的两种方式,利用傅里叶变换可以将自相关函数变换为功率谱。在脉冲数为16时,在噪声中加入SNR为6 dB的仿真尾流回波,噪声和尾流回波的自相关函数如图4(a)所示,功率谱如图4(b)所示。

    图  4  脉冲个数为16时的自相关函数与功率谱
    Figure  4.  Autocorrelation function and power spectrum when the pulse number is 16

    当尾流目标距离雷达较远时,需要进行长时间积累,实际使用的脉冲数较多。在脉冲数为64,SNR为6 dB时,利用协方差矩阵中延迟为0至15的相关系数估计功率谱,与直接利用FFT得到的功率谱进行比较,两种方法得到的功率谱如图5所示。可见,两种方法得到的功率谱较为一致,所以可以利用相关系数进行协方差矩阵的降维处理。该矩阵降维算法的缺点是功率谱的旁瓣较高,下一步将研究性能更优的矩阵降维算法。

    图  5  协方差矩阵降低维数之前与之后的功率谱
    Figure  5.  The power spectrum before and after the covariance matrix reduces the dimension

    下面以 64个脉冲为例,利用上节的方法降低协方差矩阵的维数为16,在矩阵流形上分析尾流目标的测地线距离随信噪比变化、黎曼均值估计性能和检测性能。

    在噪声背景下,计算仿真的尾流回波的协方差矩阵与噪声的协方差矩阵间的测地线距离,对100次仿真结果进行平均,图6给出了测地线距离 {D^2}\left( {{{R}_{\rm{D}}},\bar{{R}}} \right) 随信噪比变化的曲线。测地线距离随着信噪比增大而增大,所以取适当的检测门限,在信噪比大于某一确定值时,就能够检测到尾流目标。

    图  6  随SNR的变化测地线距离的变化
    Figure  6.  The variation of geodesic distance with the variation of SNR

    利用仿真噪声数据对矩阵均值迭代算法的估计性能进行仿真分析,其中已知矩阵的个数N=16,选择矩阵的算术均值作为迭代算法的初始值 {\bar{{R}}_0} ,步长 \tau = 0.1 ,可得协方差矩阵均值的估计性能如图7所示。由图可知,随着迭代次数的增加,均方测地线距离,即经验方法,不断减小,一般在20次迭代后,趋向于一个常数,利用梯度下降算法可以得到更符合物理意义的矩阵均值。

    图  7  矩阵均值迭代估计性能
    Figure  7.  Iterative estimation performance of matrix mean

    对于仿真的噪声和尾流回波,令检测统计量为尾流回波加噪声的协方差矩阵与噪声的协方差矩阵均值间的测地线距离,利用统计方法得到检测统计量的统计直方图,以及检测门限。令虚警概率为10–4,脉冲数目为16时,通过Monte-Carlo仿真得到检测门限为44.3999;脉冲数目为64时,检测门限为11.9567;脉冲数目为256时,检测门限为2.7531。脉冲数为16, SNR为6 dB时,噪声条件下和存在尾流条件下检测统计量的统计直方图如图8(a)所示,仿真1000次,检测概率为0.59。脉冲数为64, SNR为0 dB时,噪声条件下和存在尾流条件下检测统计量的统计直方图如图8(b)所示,仿真1000次,检测概率为0.79。

    图  8  有无尾流情况下检测统计量的统计直方图
    Figure  8.  The statistical histogram of detection statistics with and without the vortex target

    在脉冲数分别为16, 64和256时,利用矩阵CFAR检测方法进行尾流目标检测,其中,脉冲数大于16时,均降低协方差矩阵的维数为16。随着信噪比变化,尾流目标的检测概率如图9所示,可知,脉冲数为64时与脉冲数为16时相比,以及脉冲数为256时与脉冲数为64时相比,在相同的检测概率处检测增益均提高约6 dB,增加4倍脉冲数所获得的增益相当于相参积累所获得的增益,所以利用延迟为0至15的相关系数构成16维协方差矩阵,能够达到降低协方差矩阵的维数的目的,对矩阵CFAR检测器的性能不造成影响。

    图  9  基于矩阵CFAR的尾流目标的检测概率
    Figure  9.  Detection probability of the vortex target based on matrix CFAR

    图9可知,若用矩阵CFAR检测方法对尾流目标回波进行检测,在虚警概率10–4、检测概率80%、脉冲数为64时,所需信噪比约为0 dB。

    常规MTD检测方法通过MTD加CFAR检测器来实现尾流目标功率谱的频域峰值检测,设定如下的假设条件:

    (1) 在未知功率水平的白色高斯噪声及单脉冲平方律检测的假设下,背景功率水平由16个参考单元采样的均值估计得到,对非起伏目标,常规CA-CFAR检测器在虚警概率10–4、检测概率80%时,检测常规点目标所需的信噪比为11.6 dB;

    (2) 在信噪比为0 dB时,进行100次尾流谱仿真,统计出尾流功率谱的峰值与所有功率谱值之和的比值,由统计平均可知,尾流谱的峰值约占整个尾流谱20%的功率,理论上可做出假设,由于尾流谱的展宽,导致的检测损失约7 dB;

    (3) 在脉冲数为64时,经过相参积累,增益18 dB。

    综上,可以近似分析出,对于尾流目标功率谱的常规MTD检测方法,在虚警概率10–4、检测概率80%、脉冲数为64时,所需信噪比约为0.6 dB。

    由上面的分析可知,矩阵CFAR检测方法与常规方法相比,检测性能略有改善。另外,在信噪比为0 dB时,受噪声随机性的影响,100次仿真中有若干次,尾流目标功率谱的峰值与噪声相当,实际上常规方法的虚警概率将略高于10–4,检测性能还将低于上面的理论分析结果。更进一步,在使用常规MTD检测方法时,若参考单元选择不当,展宽的尾流目标功率谱落在参考单元中,还将抬高检测门限,导致目标遮蔽现象和检测概率下降。

    矩阵CFAR检测方法直接利用观测样本序列的协方差矩阵进行尾流目标检测,避免了常规多普勒处理所面临的问题,仿真实验证明了它在尾流目标检测中具有较大优势。矩阵CFAR检测方法的缺点是运算量较大,导致较多的处理时间,有待改善实时性,可能的解决思路是通过算术均值代替黎曼均值,利用协方差矩阵为Toeplitz Hermitian正定矩阵的特点,利用矩阵谱分解方法提高矩阵运算的效率等方面。

    飞机尾流是一种典型的角度、距离、多普勒展宽目标,同时也是一种空气介质目标,信噪比很低,针对某个特定的角度和距离单元上多普勒展宽目标的检测问题,本文研究了矩阵信息几何方法在飞机尾流目标检测中的应用。为保留原始观测信息,将协方差矩阵看作统计流形上的点,以测地线距离来定量描述目标和噪声杂波协方差矩阵间的差异,将统计检测问题转化为几何问题,提供了解决目标检测问题的一个新思路。对于功率谱展宽目标,该方法优于常规MTD检测方法。

    信息几何方法在雷达目标探测技术领域的应用仍处于起步阶段,除了在尾流目标检测、近海小目标检测、短脉冲串信号检测等方面的应用,在雷达信号处理和数据处理中也具有较广的应用前景,在信息几何应用中的理论问题有待深入探讨。信息几何涉及较深的数学知识,尤其是微分几何、统计流形等方面,由于黎曼度量、测地线、测地距离、黎曼均值等几何概念的物理意义较明确,在应用中发挥了重要作用,而信息几何中更多概念的物理意义还需要深入分析。在雷达信号处理和目标检测领域中,信息几何的优势与应用条件、非高斯概率分布流形的几何结构、基于功率谱的信息几何检测器、气象杂波环境下尾流目标的信息几何检测方法、矩阵CFAR的快速算法、适用于矩阵信息几何的矩阵降维算法等方面还需要进一步研究。

  • 图  1  含有沿时间维度突变信号的SAR幅度时间序列

    Figure  1.  SAR amplitude time series with outliers along the temporal dimension

    图  2  本文所提方法流程图

    Figure  2.  Flow chart of the proposed method

    图  3  本文采用的仿真数据

    Figure  3.  Simulation data used in this paper

    图  4  仿真数据相干斑抑制结果

    Figure  4.  Speckle suppression results of the simulated data

    图  5  本文方法的参数敏感性分析

    Figure  5.  Parameter sensitivity analyses of the proposed method

    图  6  所提方法对于图像数目的敏感性分析

    Figure  6.  Sensitivity analysis of the proposed method with respect to the number of images

    图  7  上海近海区域哨兵1号数据的相干斑抑制结果

    Figure  7.  Speckle suppression results for Sentinel-1 data in the Shanghai offshore region

    图  8  上海浦东机场区域哨兵1号数据的相干斑抑制结果

    Figure  8.  Speckle suppression results for Sentinel-1 data in Shanghai pudong airport region

    图  9  近海区域图像边缘提取结果对比

    Figure  9.  Comparison of image edge extraction results in offshore areas

    图  10  海上舰船区域图像边缘提取结果对比

    Figure  10.  Comparison of image edge extraction results of ship region

    表  1  仿真数据相干斑抑制结果的定量分析

    Table  1.   Quantitative analysis of speckle suppression results in simulation data

    方法PSNR (dB)MSSIM (dB)MCC
    0.5%1.0%0.5%1.0%0.5%1.0%
    SqueeSAR30.9529.8128.420.510.430.380.870.800.75
    MSAR-BM3D35.0118.0215.120.560.260.180.930.350.26
    RABASAR35.0732.9730.510.630.490.390.960.910.82
    DespecKS-NLLRTV29.5328.7028.410.540.310.280.900.880.88
    本文方法32.0732.8233.150.570.560.550.920.920.91
    下载: 导出CSV

    表  2  4块同质区域计算得到的等效视数

    Table  2.   The equivalent apparent number calculated by four homogeneous regions

    方法ENL
    A1A2A3A4
    原始图像0.8720.7940.9070.936
    SqueeSAR325.303.2564.3372.55
    MSAR-BM3D176.4015.4865.22117.07
    RABASAR26.3324.4024.1317.07
    DespecKS-NLLRTV54.8659.32110.7085.61
    本文方法107.3074.6880.11139.50
    下载: 导出CSV

    表  3  沿时间维度突变信号的熵值分析

    Table  3.   Entropy analysis of the outliers along the temporal dimension

    地点 测试日期DespecKS-NLLRTV本文方法
    近海2020-09-1514.6180.306
    2021-02-0615.7910.144
    机场2019-12-0218.6230.643
    2020-09-2716.0320.736
    下载: 导出CSV
  • [1] MOREIRA A, PRATS-IRAOLA P, YOUNIS M, et al. A tutorial on synthetic aperture radar[J]. IEEE Geoscience and Remote Sensing Magazine, 2013, 1(1): 6–43. doi: 10.1109/MGRS.2013.2248301
    [2] 邓云凯, 赵凤军, 王宇. 星载SAR技术的发展趋势及应用浅析[J]. 雷达学报, 2012, 1(1): 1–10. doi: 10.3724/SP.J.1300.2012.20015

    DENG Yunkai, ZHAO Fengjun, and WANG Yu. Brief analysis on the development and application of spaceborne SAR[J]. Journal of Radars, 2012, 1(1): 1–10. doi: 10.3724/SP.J.1300.2012.20015
    [3] 陈杰, 杨威, 王亚敏, 等. 高时相星载序贯SAR图像运动目标检测方法[J]. 雷达学报, 2022, 11(6): 1048–1060. doi: 10.12000/JR22184

    CHEN Jie, YANG Wei, WANG Yamin, et al. Moving target monitoring algorithm based on high-frame-rate SAR images[J]. Journal of Radars, 2022, 11(6): 1048–1060. doi: 10.12000/JR22184
    [4] XUE Feiyang, LV Xiaolei, DOU Fangjia, et al. A review of time-series interferometric SAR techniques: A tutorial for surface deformation analysis[J]. IEEE Geoscience and Remote Sensing Magazine, 2020, 8(1): 22–42. doi: 10.1109/MGRS.2019.2956165
    [5] 徐丰, 王海鹏, 金亚秋. 深度学习在SAR目标识别与地物分类中的应用[J]. 雷达学报, 2017, 6(2): 136–148. doi: 10.12000/JR16130

    XU Feng, WANG Haipeng, and JIN Yaqiu. Deep learning as applied in SAR target recognition and terrain classification[J]. Journal of Radars, 2017, 6(2): 136–148. doi: 10.12000/JR16130
    [6] ARGENTI F, LAPINI A, BIANCHI T, et al. A tutorial on speckle reduction in synthetic aperture radar images[J]. IEEE Geoscience and Remote Sensing Magazine, 2013, 1(3): 6–35. doi: 10.1109/MGRS.2013.2277512
    [7] VASILE G, TROUVE E, LEE J S, et al. Intensity-driven adaptive-neighborhood technique for polarimetric and interferometric SAR parameters estimation[J]. IEEE Transactions on Geoscience and Remote Sensing, 2006, 44(6): 1609–1621. doi: 10.1109/TGRS.2005.864142
    [8] LEE J S, GRUNES M R, and GRANDI G D. Polarimetric SAR speckle filtering and its implication for classification[J]. IEEE Transactions on Geoscience and Remote Sensing, 1999, 37(5): 2363–2373. doi: 10.1109/36.789635
    [9] BIOUCAS-DIAS J M and FIGUEIREDO M A T. Multiplicative noise removal using variable splitting and constrained optimization[J]. IEEE Transactions on Image Processing, 2010, 19(7): 1720–1730. doi: 10.1109/TIP.2010.2045029
    [10] DALSASSO E, YANG Xiangli, DENIS L, et al. SAR image despeckling by deep neural networks: From a pre-trained model to an end-to-end training strategy[J]. Remote Sensing, 2020, 12(16): 2636. doi: 10.3390/rs12162636
    [11] XIONG Kai, ZHAO Guanghui, WANG Yingbin, et al. SPB-Net: A deep network for SAR imaging and despeckling with downsampled data[J]. IEEE Transactions on Geoscience and Remote Sensing, 2021, 59(11): 9238–9256. doi: 10.1109/TGRS.2020.3034102
    [12] MOLINI A B, VALSESIA D, FRACASTORO G, et al. Speckle2Void: Deep self-supervised SAR despeckling with blind-spot convolutional neural networks[J]. IEEE Transactions on Geoscience and Remote Sensing, 2022, 60: 5204017. doi: 10.1109/TGRS.2021.3065461
    [13] DELEDALLE C A, DENIS L, POGGI G, et al. Exploiting patch similarity for SAR image processing: The nonlocal paradigm[J]. IEEE Signal Processing Magazine, 2014, 31(4): 69–78. doi: 10.1109/MSP.2014.2311305
    [14] XU Gang, GAO Yandong, LI Jinwei, et al. InSAR phase denoising: A review of current technologies and future directions[J]. IEEE Geoscience and Remote Sensing Magazine, 2020, 8(2): 64–82. doi: 10.1109/MGRS.2019.2955120
    [15] QUEGAN S, LE TOAN T, YU Jiongjiong, et al. Multitemporal ERS SAR analysis applied to forest mapping[J]. IEEE Transactions on Geoscience and Remote Sensing, 2000, 38(2): 741–753. doi: 10.1109/36.842003
    [16] TROUVE E, CHAMBENOIT Y, CLASSEAU N, et al. Statistical and operational performance assessment of multitemporal SAR image filtering[J]. IEEE Transactions on Geoscience and Remote Sensing, 2003, 41(11): 2519–2530. doi: 10.1109/TGRS.2003.817270
    [17] SU Xin, DELEDALLE C A, TUPIN F, et al. Two-step multitemporal nonlocal means for synthetic aperture radar images[J]. IEEE Transactions on Geoscience and Remote Sensing, 2014, 52(10): 6181–6196. doi: 10.1109/TGRS.2013.2295431
    [18] FERRETTI A, FUMAGALLI A, NOVALI F, et al. A new algorithm for processing interferometric data-stacks: SqueeSAR[J]. IEEE Transactions on Geoscience and Remote Sensing, 2011, 49(9): 3460–3470. doi: 10.1109/TGRS.2011.2124465
    [19] CHIERCHIA G, EL GHECHE M, SCARPA G, et al. Multitemporal SAR image despeckling based on block-matching and collaborative filtering[J]. IEEE Transactions on Geoscience and Remote Sensing, 2017, 55(10): 5467–5480. doi: 10.1109/TGRS.2017.2707806
    [20] PARRILLI S, PODERICO M, ANGELINO C V, et al. A nonlocal SAR image denoising algorithm based on LLMMSE wavelet shrinkage[J]. IEEE Transactions on Geoscience and Remote Sensing, 2012, 50(2): 606–616. doi: 10.1109/TGRS.2011.2161586
    [21] ZHAO Weiying, DELEDALLE C A, DENIS L, et al. Ratio-based multitemporal SAR images denoising: RABASAR[J]. IEEE Transactions on Geoscience and Remote Sensing, 2019, 57(6): 3552–3565. doi: 10.1109/TGRS.2018.2885683
    [22] BAIER G, HE Wei, and YOKOYA N. Robust nonlocal low-rank SAR time series despeckling considering speckle correlation by total variation regularization[J]. IEEE Transactions on Geoscience and Remote Sensing, 2020, 58(11): 7942–7954. doi: 10.1109/TGRS.2020.2985400
    [23] CHAMBOLLE A and POCK T. A first-order primal-dual algorithm for convex problems with applications to imaging[J]. Journal of Mathematical Imaging and Vision, 2011, 40(1/2): 120–145. doi: 10.1007/s10851-010-0251-1
    [24] BOYD S, PARIKH N, CHU E, et al. Distributed optimization and statistical learning via the alternating direction method of multipliers[J]. Foundations and Trends ® in Machine Learning, 2011, 3(1): 1–122. doi: 10.1561/2200000016
    [25] FERRETTI A, PRATI C, and ROCCA F. Permanent scatterers in SAR interferometry[J]. IEEE Transactions on Geoscience and Remote Sensing, 2001, 39(1): 8–20. doi: 10.1109/36.898661
    [26] DELEDALLE C A, DENIS L, TABTI S, et al. MuLoG, or how to apply Gaussian denoisers to multi-channel SAR speckle reduction?[J]. IEEE Transactions on Image Processing, 2017, 26(9): 4389–4403. doi: 10.1109/TIP.2017.2713946
    [27] WANG Yilun, YANG Junfeng, YIN Wotao, et al. A new alternating minimization algorithm for total variation image reconstruction[J]. SIAM Journal on Imaging Sciences, 2008, 1(3): 248–272. doi: 10.1137/080724265
    [28] PARIKH N and BOYD S. Proximal algorithms[J]. Foundations and Trends ® in Optimization, 2014, 1(3): 127–239. doi: 10.1561/2400000003
  • 加载中
图(10) / 表(3)
计量
  • 文章访问数: 916
  • HTML全文浏览量: 312
  • PDF下载量: 184
  • 被引次数: 0
出版历程
  • 收稿日期:  2022-12-25
  • 修回日期:  2023-02-28
  • 网络出版日期:  2023-03-14
  • 刊出日期:  2023-10-28

目录

/

返回文章
返回