CN110726992B - 基于结构稀疏和熵联合约束的sa-isar自聚焦法 - Google Patents

Abstract

本发明属于雷达信号处理领域，具体涉及一种基于结构稀疏和熵联合约束的SA‑ISAR自聚焦法，包括以下步骤：S1：对包络对齐后的雷达回波进行回波建模；S2：对ISAR图像施加分层的结构化稀疏先验；S3：通过松弛的变分贝叶斯方法对ISAR图像和高层变量进行更新；S4：通过基于不动点的最小熵法对相位误差进行更新；S5：判断是否达到终止条件，如果达到终止条件则停止迭代循环，否则返回到S3，直到达到终止条件后输出自聚焦之后的图像。本发明可以提升稀疏孔径下ISAR图像的自聚焦精度，使得所成ISAR图像更为清晰；本发明的计算复杂度较低、迭代收敛速度较快，对抗稀疏孔径的能力强。

Description

基于结构稀疏和熵联合约束的SA-ISAR自聚焦法

技术领域

本发明属于雷达信号处理领域，具体涉及一种基于结构稀疏和熵联合约束的稀疏孔径逆合成孔径雷达(Sparse Aperture-Inverse Synthetic Aperture Radar，SA-ISAR)自聚焦法。

背景技术

成像技术可获取运动目标的高分辨率二维图像，从而捕获目标二维尺寸、结构特征，是空间目标识别的重要技术手段，已在空间目标检测、导弹防御、航空管制、雷达天文学等领域获得广泛应用。

在雷达成像领域，成像区间内脉冲完整的信号称为全孔径信号，而如果回波信号存在随机或者成段脉冲缺失，则称之为稀疏孔径(SA)信号，稀疏孔径下的ISAR成像简称为SA-ISAR。在ISAR系统中，许多因素均可造成回波信号孔径稀疏：其一，由目标距离较远、目标尺寸较小、以及复杂的空间电磁环境等导致的低信噪比将导致部分回波脉冲缺失；其二，目标的复杂运动会导致可用于成像的相干积累时间较短，造成孔径稀疏；其三，在激烈的空间攻防对抗条件下，越来越多样的干扰措施也将导致部分回波脉冲不可用；其四，随着雷达技术的不断提高，多功能雷达的广泛应用也会造成一定程度的孔径稀疏。为同时实现对目标的搜索、跟踪与成像，多功能雷达系统多采用“宽-窄”交替的工作模式，即交替发射窄带与宽带信号，通过窄带信号对目标进行跟踪，测量目标位置与速度，而通过宽带信号对目标进行成像，以获取目标尺寸结构信息。为降低硬件要求，这种工作模式下的雷达多间歇性发射宽带信号，从而导致回波的孔径稀疏。

对于全孔径信号，传统距离-多普勒成像算法(RD成像算法)可获得理想的ISAR图像，然而对于稀疏孔径信号，RD成像算法所获得的图像将受到严重的旁瓣、杂波干扰与主瓣展宽，导致分辨率降低，难以满足工程需求。另外，稀疏孔径还严重影响ISAR自聚焦性能，导致ISAR图像严重散焦。稀疏孔径条件下ISAR自聚焦技术，对改善稀疏孔径条件下ISAR成像质量具有重要的工程应用价值。

发明内容

本发明要解决的技术问题是，在稀疏孔径条件下，ISAR自聚焦性能下降，导致ISAR图像质量降低，难以满足工程实际需求。

本发明的思路是针对稀疏孔径条件下ISAR图像质量下降、分辨率降低的问题，提出一种基于结构稀疏和熵联合约束的SA-ISAR自聚焦法。该方法首先对ISAR图像进行先验建模，提出一种结构化的分层稀疏先验模型。然后采用变分贝叶斯推断对隐变量的后验参数进行学习，并采用一定的松弛方法对算法进行优化加速；对于初相误差，采用最小熵准则下的不动点法在后验参数迭代过程中进行估计和补偿。

为了公式表达的方便，对全文所使用的符号作出以下统一说明：

用矩阵的下脚标来表示对矩阵的索引，即用A_·,n表示矩阵A第n列的列向量，用A_m,·表示矩阵A第m行的行向量，用A_m,n表示矩阵A的第m行第n列元素，用A_·,a表示由矩阵A标号为a的所有列向量组成的子矩阵，用a_i表示向量a的第i个元素(为了表示方便，标号均从0开始)；用diag(a)表示以n维向量a为对角线元素且其它元素为0的n×n方阵，用diag(A)表示方阵对角线元素组成的列向量；用A^T表示矩阵A的转置，用A^H表示矩阵A的Hermite转置；用tr(A)表示方阵A的迹；用vec(A)表示将A按列向量展开；用I_N表示N×N的单位矩阵；用||A||_F表示矩阵A的Frobenius范数(简称F-范数)；用||a||₂表示向量a的2范数；用

表示矩阵A和B的Kronecker乘积；用

表示同维矩阵A和B的Hadamard乘积；用

表示矩阵A和B(不含零元素)元素对应相除。以上概念在任何一本矩阵分析的教材中都有详细介绍。

用p(x₁,x₂,...,x_M|z₁,z₂,...,z_N；θ₁,θ₂,...,θ_K)表示x₁,x₂,...,x_M在随机变量z₁,z₂,...,z_N和参量θ₁,θ₂,...,θ_K下的联合概率密度函数。用

表示关于随机变量x，均值为u，协方差矩阵为Σ的多元复高斯分布；用

表示关于随机变量x，两参数分别为a,b的伽马分布；用

表示关于随机变量x，三个参数分别为λ,a,b的广义逆高斯分布。用<x>表示对x求期望。以上概念在任何一本概率论的教材中都有详细介绍。

本发明解决其技术问题所采取的技术方案是：一种基于结构稀疏和熵联合约束的SA-ISAR自聚焦法，包括以下步骤：

S1对包络对齐以后的雷达回波进行信号建模；

在稀疏孔径条件下，传统的最小熵/互相关包络对齐方法依然有效，本发明直接对已经完成包络对齐后的目标回波进行建模。在忽略目标的越分辨单元走动时，包络对齐后的目标回波可以近似离散表示为：

其中，s(m,n)表示回波信号，m(m＝0,1,...,M-1)为脉冲序号，n(n＝0,1,...,N-1)为快时间采样序号，M为全孔径的脉冲数，N为一个脉冲内总的采样数；σ(x,y)是目标的复散射分布，表示目标在成像平面坐标(x,y)处的复散射强度，C_T表示目标在成像平面的投影分布区域，λ_c是雷达载波的波长，K是发射线性调频信号的调频斜率，T_s是快时间的采样间隔，C是真空中光的传播速度，T_r是雷达发射信号的脉冲重复间隔，ω表示目标的等效转动角速度，

是第m个脉冲的初相误差；

将成像区域离散化，令x＝u·λ_c/(2ωT_rM)，y＝v·C/(KT_sN)，其中u(u＝0,1,...,M-1)是图像横向单元标号，v(v＝0,1,...,N-1)是图像距离向单元标号，假定目标分布区域满足

(1)进一步离散化(忽略了幅值的线性变化)为：

其中

将(2)矩阵化表示为：

Y^FA＝F^NX(F^M)^TE (3)

其中Y^FA为回波组成的N×M维矩阵，且Y^FA _n,m＝s(m,n)；X是离散化的目标二维散射分布矩阵，其维数为N×M，且

E是M×M的初相误差矩阵，且

F^M和F^N分别为M阶和N阶的傅里叶变换矩阵，即F^M _m,u＝exp(-j2πmu/M)，F^N _n,v＝exp(-j2πnv/N)；

将(3)向量化表示为：

y^FA＝Φx (4)

其中y^FA＝vec(Y^FA)，x＝vec(X)，

在有脉冲缺失时，即稀疏孔径条件下，(4)中M个脉冲中只有部分可以获得，假设可获得L(L＜M)个脉冲，其在全孔径中的位置下标为D，D是一个L维的列向量，每一个元素在0到M-1之间，代表了相应稀疏孔径脉冲的序号；这样可以把稀疏孔径回波采样表示为：

Y＝Y^FA _·,D (5)

根据孔径的稀疏采样情况可以构造矩阵

其每一行有且仅有一个元素为1，其余元素为0，代表相应位置的回波采样可以获得；

表示LN×MN的实数矩阵；这样稀疏孔径的回波可以写为：

y＝vec(Y)＝PΦx (6)

其中，

表示稀疏孔径回波的向量化展开，

表示LN×1的复数矩阵；

通常，将回波中混有的噪声建模为复高斯噪声，即

其中β表示回波噪声方差的倒数；

噪声方差的倒数通常加入伽马分布先验，如下式(8)，其中参数c,d通常取一个小值c＝d＝10^-7：

以上就完成了对稀疏孔径回波的信号建模；

S2对ISAR图像施加分层的结构化稀疏先验；

为了充分挖掘目标的先验信息，本步骤设计针对ISAR图像的结构化稀疏先验模型，具体过程如下：

S2.1对ISAR图像施加零均值高斯先验分布：

其中，α是一个N×M的矩阵，其每一个元素都是X对应位置元素的先验方差；

S2.2对α施加伽马先验分布：

其中，γ是一个(N+2)×(M+2)的矩阵，其每一个元素都会以一定权重影响α对应相邻位置的先验均值和先验方差；

表示结构化先验的依赖模式，θ是相邻元素相互依赖的权重因子，通常取θ＝0.5，每一个Pattern和θ值可以对应一个平滑矩阵

其形式如下：

例如可取Pattern＝{(-1,0),(1,0),(0,-1),(0,1)}，这时Q取值为

S2.3对γ施加伽马先验分布：

其中n′n′＝0,1,...,N+1和m′m′＝0,1,...,M+1分别为矩阵γ在两个维度的索引变量；参数a通常取一个比1略大的值，例如a＝1+10^-7，b通常取一个小值，例如b＝10^-7

S3通过松弛的变分贝叶斯方法对ISAR图像和高层变量进行更新；

本步骤包括对ISAR图像和图像后验方差进行更新、对高层变量α、γ、β进行更新、对模型参数Z进行更新；

对于变分贝叶斯推断，通常采用证据下界(ELBO)作为优化目标，其表达式如下：

其中，p(y,x,α,γ,β)是模型中所有随机变量的联合概率密度函数，q(x,α,γ,β)是所有隐变量待优化的联合后验分布函数，其通常假设为可因式分解的形式：

式中q_x(x)、q_α(α)、q_γ(γ)、q_β(β)分别为x、α、γ、β的待优化后验分布函数；

分别是α_n,m、γ_n′,m′的待优化后验分布函数；

直接对(13)进行迭代寻优会出现大矩阵求逆使得计算量大大增加，下面对(13)进行松弛，以降低计算复杂度，考虑不等式：

其中z是新引入的与x同维数的中间参量，T是与缩放程度相关的一个参量，取值为2MN；对回波的似然函数p(y|x,β)进行缩放，得到松弛的似然函数G(y,x,β；z)：

将(13)中的p(y|x,β)用G(y,x,β；z)替换，得到松弛后的优化目标为reELBO(q；z)，其表达式为：

为了方便表示，定义一个维数为N×M的矩阵v：

根据变分贝叶斯推断的相关理论，最优的q_x(x)、

q_β(β)求解为q_x ^*(x)，

q_β ^*(β)，其分别为多元高斯分布、多元广义逆高斯分布、多元伽马分布，其表达形式如下：

其中期望值

的计算公式为：

期望值<|X_n,m|²>的计算公式为：

式(19)中u是x的后验均值，Σ是x的后验协方差矩阵，计算公式如下：

其中，1_N×M是维数为N×M的全1矩阵，I_MN是MN维的全1列向量，

是一个N×M的矩阵，它的第n行第m列的元素可以由下式计算：

根据(19)～(21)和概率学知识可求出x、α_n,m、β的后验均值如(27)～(29)所示，它们将作为变量迭代更新的规则：

<x>＝u (27)

κ_-1/2(·)、κ_3/2(·)和κ_1/2(·)分别表示阶数为-1/2、3/2和1/2的第二类修正贝塞尔函数；

对于γ，由于结构化依赖关系的引入，无法实现γ的完全变分贝叶斯推断，而采用变分贝叶斯最大期望法(VBEM)得到它的更新规则为：

其中χ为(N+4)×(M+4)的矩阵：

对于z，类似地，可以得到变量z基于VBEM的更新规则为：

z＝u (32)

在采用松弛的目标函数之后，Σ变成一个对角矩阵，对其的更新可以只看作对角元素的更新，这MN个对角元素与x的方差一一对应；将Σ的对角元素用矩阵表示为diagΣ＝vecΓ，那么Γ的元素就与X的方差一一对应，此时可以用维数较低的N×M矩阵实现迭代更新，而避免大维数矩阵(NM×NM)的求逆；将(24)(28)(29)(30)(32)的迭代规则用维数较低的N×M矩阵实现，描述如下(忽略中间的数学推导过程)：

S3.1：更新ISAR图像的后验方差Γ，后验均值X：

式(33)由(25)改写而来，矩阵

各个元素的更新可由式(26)实现；式(34)由(24)改写而来，且利用了傅里叶变换的快速算法加速实现，FFT2(·)表示N×M的二维快速傅里叶变换，IFFT2(·)表示N×M的二维快速傅里叶逆变换，其定义可见于任何一本《数字信号处理》教材；(·)_·,D表示按稀疏采样模式D，将N×M矩阵转化为N×L的矩阵，

表示按照稀疏采样模式D将缺失的脉冲位置垫零，从而将N×L的矩阵转化为N×M的矩阵，Z是z的矩阵形式，z＝vec(Z)；

S3.2：更新高层变量γ：首先通过式(31)对矩阵<α>进行延拓得到χ；再利用式(30)对γ更新，用N×M矩阵重新表达如下：

其中，χ*Q表示以Q为模板计算χ的valid模式的二维卷积，这一概念可见于任何一本《数字图像处理》的教材；

S3.3：更新高层变量α：首先根据(18)和(23)计算中间量，用N×M矩阵重新表达为v＝γ*Q，<|X|²>＝|<X>|²+Γ，再由(28)和(26)得到<α>和

用N×M矩阵重新表达为：

式中的(·)^1/2、和κ_λ(·),λ＝-1/2,1/2,3/2都是指对应到元素的运算；

S3.4：更新噪声方差β：根据式(29)即可得到，写成N×M矩阵形式为：

其中

其中sum(·)表示对矩阵的所有元素求和；

S3.5：更新参数Z：将(32)写成N×M的矩阵形式为Z＝<X>； (42)

S4：通过基于不动点的最小熵法对初相误差

进行更新；

众多研究文献中已经表明，最小熵准则能够较好地补偿ISAR回波的初始相位。本步骤在前述S3的贝叶斯稀疏恢复的框架下设计基于不动点的最小熵相位误差更新方法。

ISAR图像熵的表达式为：

其中的(·)²、|·|和log(·)分别表示对应到元素的取平方、取模和以自然常数e为底取对数；

Entropy(<X>)关于

的偏导数计算为：

其中对矩阵的偏导数表示对应到元素的偏导数，对复数的偏导数指对实部和虚部分别求导，Re(·)表示对应到元素取实部，(·)^*表示对应到元素取共轭；

由(34)可以求得：

其中IFFT(·)表示对矩阵的每一列进行一维快速傅里叶变换；

令(44)等于零可以求得初相误差的更新规则为：

S5：判断是否达到终止条件，如果达到终止条件则停止迭代循环，否则返回到S3，直到达到终止条件后输出自聚焦之后的图像。

终止条件的判断可以采用指定迭代次数法，例如指定迭代次数为30次，在达到30次循环后就终止迭代；也可以采用门限法，当ISAR图像的更新变化小于某个门限值时就停止迭代，用公式表示如下

其中x^(new)和x^(old)分别表示本次迭代和上一次迭代计算的ISAR图像<X>，门限值Threshold根据实际情况选取，例如可以取10^-6。

本发明取得的有益效果为：本发明可以提升稀疏孔径下ISAR图像的自聚焦精度，使得所成ISAR图像更为清晰。本发明的计算复杂度较低、迭代收敛速度较快，在不进行初相误差粗估计的条件下，对256×256大小的回波数据能够在5s的计算时间内达到良好的图像聚焦效果；对抗稀疏孔径的能力强，即使在孔径采样比例仅为3/16的条件下依然能够有效重构出ISAR图像。

附图说明

图1本发明流程S2中的结构化稀疏先验模型示意图；

图2本发明的所提方法流程示意图；

图3C波段Yak-42数据目标光学图像和全孔径ISAR图像；

图4不同稀疏采样个数下的Yak-42自聚焦成像结果图；

图5不同方法在不同稀疏采样个数下的迭代时间对比图；

图6X波段民航飞机数据车载雷达示意图、目标光学图像、全孔径一维距离像、全孔径ISAR图像；

图7不同稀疏采样个数下的X波段雷达数据自聚焦成像结果图。

具体实施方式

下面结合附图对本发明进行进一步说明：

图1为S2中ISAR图像的结构化分层稀疏先验模型示意图。图中的X代表ISAR图像，α和γ代表高层变量，x代表向量化展开之后的ISAR图像，y代表包络对齐之后稀疏孔径回波的向量化展开。图中的实线箭头代表变量之间概率分布的依赖关系，α和γ之间的“Convolutional Connected”代表α_m,n的分布依赖于γ_m+1,n+1和相邻的若干元素，X和α之间的“Pointwise Connected”代表X_m,n的分布依赖于α_m,n，y和x之间的“Full Connected”代表y_i的分布依赖于x的所有元素，x与X之间标注为“Vectorize”的虚线箭头不代表概率分布的依赖关系，而代表x由X向量化展开而来，两者拥有相同的概率分布。

图2为本发明所提方法的流程示意图，包括以下步骤：

S1：对包络对齐后的雷达回波进行回波建模；

S2：对ISAR图像施加分层的结构化稀疏先验；

S3：通过松弛的变分贝叶斯方法对ISAR图像和高层变量进行更新；

S4：通过基于不动点的最小熵法对初相误差进行更新；

S5：判断是否达到终止条件，如果达到终止条件则停止迭代循环，否则返回到S3，直到达到终止条件后输出自聚焦之后的图像。

为了验证算法的有效性，先采用C波段实测数据进行验证，C波段雷达工作频率为5.52GHz，带宽为400MHz，宽带信号脉冲宽度为25.6us，雷达的脉冲重复频率为100Hz，去斜接收的中频信号快时间采样频率为10MHz。所照射的目标是型号为“Yak-42”的民航客机，快时间采样点数为306，全孔径脉冲数为256。图3左图是目标的光学图像，右图是目标全孔径数据在自聚焦之后的ISAR图像。

图4是不同稀疏孔径比(50％和25％)下的自聚焦成像结果图。对图3中的全孔径数据保留初相误差进行稀疏采样，左侧两幅图示出了稀疏孔径的采样模式，图中每一列代表一个一维距离像，一维像全为黑色代表孔径缺失；右侧两幅图是使用本发明对稀疏采样的回波数据自聚焦成像的结果，迭代次数为30，可以看到不论采样比为50％还是25％，本发明都能够有效补偿初相误差、重构出ISAR图像，且重构出的图像与全孔径得到的结果基本相同。

图5是使用该段“Yak-42”数据，对本发明所提方法的运算速度的测试结果。测试时，对全孔径数据进行不同比例的稀疏采样，进行30次迭代。“采用原始目标函数”是指不对ELBO进行松弛直接进行变分推断，“采用松弛目标函数”是指采用本发明对ELBO进行松弛再进行推断。可以看到采用本发明的松弛目标函数使得运算效率提高了一个数量级，对单幅尺寸为301×256的稀疏孔径ISAR图像自聚焦重构，本发明的消耗时间在4s左右。测试采用MATLAB2018a作为软件平台，处理器为Intel(R)Core(TM)i5-8250 CPU。

再使用X波段实测数据进行验证，图6左上图为采用的车载X波段雷达系统，该雷达发射LFM信号，中心频率为9GHz，带宽为1GHz，脉宽为100μs，PRF为100Hz。目标为某民航客机，其光学图像如图6右上。脉内采样点数为256，全孔径包含256个脉冲。同样，作为稀疏孔径条件下ISAR成像质量的参照标准，图6左下与图6右下分别展示了全孔径条件下目标包络对齐后的一维像序列与ISAR像。

图7是不同稀疏孔径比(50％和25％)下的自聚焦成像结果图。各图含义类比图4。可以看到不论采样比为50％还是25％，本发明都能够有效补偿初相误差、重构出ISAR图像。(注：为了显示出弱散射中心，图6右下和图7右上、图7右下都对图像数据取了对数)。

Claims (9)

1.一种基于结构稀疏和熵联合约束的SA-ISAR自聚焦法，其特征在于，该方法包括以下步骤：

S1对包络对齐以后的雷达回波进行信号建模；

在忽略目标的越分辨单元走动时，包络对齐后的目标回波可以近似离散表示为：

其中，s(m,n)表示回波信号，m为脉冲序号，n为快时间采样序号，m＝0,1,...,M-1，n＝0,1,...,N-1，M为全孔径的脉冲数，N为一个脉冲内总的采样数；σ(x,y)是目标的复散射分布，表示目标在成像平面坐标(x,y)处的复散射强度，C_T表示目标在成像平面的投影分布区域，λ_c是雷达载波的波长，K是发射线性调频信号的调频斜率，T_s是快时间的采样间隔，C是真空中光的传播速度，T_r是雷达发射信号的脉冲重复间隔，ω表示目标的等效转动角速度，

是第m个脉冲的初相误差；

将成像区域离散化，令x＝u·λ_c/(2ωT_rM)，y＝v·C/(KT_sN)，其中u是图像横向单元标号，v是图像距离向单元标号，u＝0,1,...,M-1，v＝0,1,...,N-1；假定目标分布区域满足

忽略幅值的线性变化后(1)进一步离散化为：

其中

将(2)矩阵化表示为：

Y^FA＝F^NX(F^M)^TE (3)

其中Y^FA为回波组成的N×M维矩阵，且Y^FA _n,m＝s(m,n)；X是离散化的目标二维散射分布矩阵，其维数为N×M，且

E是M×M的初相误差矩阵，且

F^M和F^N分别为M阶和N阶的傅里叶变换矩阵，即F^M _m,u＝exp(-j2πmu/M)，F^N _n,v＝exp(-j2πnv/N)；

将(3)向量化表示为：

y^FA＝Φx (4)

其中y^FA＝vec(Y^FA)，x＝vec(X)，

在有脉冲缺失时，即稀疏孔径条件下，(4)中M个脉冲中只有部分可以获得，假设可获得L个脉冲，L＜M，其在全孔径中的位置下标为D，D是一个L维的列向量，每一个元素在0到M-1之间，代表了相应稀疏孔径脉冲的序号；这样可以把稀疏孔径回波采样表示为

Y＝Y^FA _·,D (5)

根据孔径的稀疏采样情况可以构造矩阵

其每一行有且仅有一个元素为1，其余元素为0，代表相应位置的回波采样可以获得；

表示LN×MN的实数矩阵；这样稀疏孔径的回波可以写为：

y＝vec(Y)＝PΦx (6)

其中，

表示稀疏孔径回波的向量化展开，

表示LN×1的复数矩阵；

通常，将回波中混有的噪声建模为复高斯噪声，即

其中β表示回波噪声方差的倒数；

噪声方差的倒数通常加入伽马分布先验，如下式(8)：

以上就完成了对稀疏孔径回波的信号建模；

S2对ISAR图像施加分层的结构化稀疏先验；

为了充分挖掘目标的先验信息，本步骤设计针对ISAR图像的结构化稀疏先验模型，具体过程如下：

S2.1对ISAR图像施加零均值高斯先验分布：

其中，α是一个N×M的矩阵，其每一个元素都是X对应位置元素的先验方差；

S2.2对α施加伽马先验分布：

其中，γ是一个(N+2)×(M+2)的矩阵，其每一个元素都会以一定权重影响α对应相邻位置的先验均值和先验方差；

表示结构化先验的依赖模式，θ是相邻元素相互依赖的权重因子，每一个Pattern和θ值可以对应一个平滑矩阵

其形式如下：

S2.3对γ施加伽马先验分布：

其中n′和m′分别为矩阵γ在两个维度的索引变量，n′＝0,1,...,N+1，m′＝0,1,...,M+1；

S3通过松弛的变分贝叶斯方法对ISAR图像和高层变量进行更新；

本步骤包括对ISAR图像和图像后验方差进行更新、对高层变量α、γ、β进行更新、对模型参数Z进行更新；

对于变分贝叶斯推断，通常采用证据下界作为优化目标，其表达式如下：

其中，p(y,x,α,γ,β)是模型中所有随机变量的联合概率密度函数，q(x,α,γ,β)是所有隐变量待优化的联合后验分布函数，其通常假设为可因式分解的形式：

式中q_x(x)、q_α(α)、q_γ(γ)、q_β(β)分别为x、α、γ、β的待优化后验分布函数；

分别是α_n,m、γ_n′,m′的待优化后验分布函数；

直接对(13)进行迭代寻优会出现大矩阵求逆使得计算量大大增加，下面对(13)进行松弛，以降低计算复杂度，考虑不等式：

其中z是新引入的与x同维数的中间参量，T是与缩放程度相关的一个参量，取值为2MN；对回波的似然函数p(y|x,β)进行缩放，得到松弛的似然函数G(y,x,β；z)：

将(13)中的p(y|x,β)用G(y,x,β；z)替换，得到松弛后的优化目标为reELBO(q；z)，其表达式为：

为了方便表示，定义一个维数为N×M的矩阵v：

根据变分贝叶斯推断的相关理论，最优的q_x(x)、

q_β(β)求解为q_x ^*(x)，

q_β ^*(β)，其分别为多元高斯分布、多元广义逆高斯分布、多元伽马分布，其表达形式如下：

其中期望值

的计算公式为：

期望值<|X_n,m|²>的计算公式为：

<|X_n,m|²>＝|u_m·N+n+1|²+Σ_{m·N+n+1,m·N+n+1} (23)

式(19)中u是x的后验均值，Σ是x的后验协方差矩阵，计算公式如下：

其中，1_N×M是维数为N×M的全1矩阵，I_MN是MN维的全1列向量，

是一个N×M的矩阵，它的第n行第m列的元素可以由下式计算：

根据(19)～(21)和概率学知识可求出x、α_n,m、β的后验均值如(27)～(29)所示，它们将作为变量迭代更新的规则：

<x>＝u (27)

κ_-1/2(·)、κ_3/2(·)和κ_1/2(·)分别表示阶数为-1/2、3/2和1/2的第二类修正贝塞尔函数；

对于γ，由于结构化依赖关系的引入，无法实现γ的完全变分贝叶斯推断，而采用变分贝叶斯最大期望法得到它的更新规则为：

其中χ为(N+4)×(M+4)的矩阵：

对于z，类似地，可以得到变量z基于VBEM的更新规则为：

z＝u (32)

将(24)(28)(29)(30)(32)的迭代规则用维数较低的N×M矩阵实现，描述如下：

S3.1：更新ISAR图像的后验方差Γ，后验均值X：

式(33)由(25)改写而来，矩阵

各个元素的更新可由式(26)实现；式(34)由(24)改写而来，且利用了傅里叶变换的快速算法加速实现，FFT2(·)表示N×M的二维快速傅里叶变换，IFFT2(·)表示N×M的二维快速傅里叶逆变换；(·)·,_D表示按稀疏采样模式D，将N×M矩阵转化为N×L的矩阵，

表示按照稀疏采样模式D将缺失的脉冲位置垫零，从而将N×L的矩阵转化为N×M的矩阵，Z是z的矩阵形式，z＝vec(Z)；

S3.2：更新高层变量γ：首先通过式(31)对矩阵<α>进行延拓得到χ；再利用式(30)对γ更新，用N×M矩阵重新表达如下：

其中，χ*Q表示以Q为模板计算χ的valid模式的二维卷积；

S3.3：更新高层变量α：首先根据(18)和(23)计算中间量，用N×M矩阵重新表达为v＝γ*Q，<|X|²>＝|<X>|²+Γ，再由(28)和(26)得到<α>和

用N×M矩阵重新表达为：

式中的(·)^1/2、和κ_λ(·),λ＝-1/2,1/2,3/2都是指对应到元素的运算；

S3.4：更新噪声方差β：根据式(29)即可得到，写成N×M矩阵形式为：

其中

S₂＝-2·sum((Y-(FFT2(Z)·E)_·,D)⊙(FFT2(<X>-Z)·E)_·,D) (40)

其中sum(·)表示对矩阵的所有元素求和；

S3.5：更新参数Z：将(32)写成N×M的矩阵形式为

Z＝<X>； (42)

S4：通过基于不动点的最小熵法对初相误差

进行更新；

ISAR图像熵的表达式为：

其中的(·)²、|·|和log(·)分别表示对应到元素的取平方、取模和以自然常数e为底取对数；

Entropy(<X>)关于

的偏导数计算为：

其中对矩阵的偏导数表示对应到元素的偏导数，对复数的偏导数指对实部和虚部分别求导，Re(·)表示对应到元素取实部，(·)^*表示对应到元素取共轭；

由(34)可以求得：

其中IFFT(·)表示对矩阵的每一列进行一维快速傅里叶变换；

令(44)等于零可以求得初相误差的更新规则为：

S5：判断是否达到终止条件，如果达到终止条件则停止迭代循环，否则返回到S3，直到达到终止条件后输出自聚焦之后的图像。

2.一种根据权利要求1所述基于结构稀疏和熵联合约束的SA-ISAR自聚焦法，其特征在于：公式(8)中的参数c,d取值为c＝d＝10^-7。

3.一种根据权利要求1所述基于结构稀疏和熵联合约束的SA-ISAR自聚焦法，其特征在于：相邻元素相互依赖的权重因子θ取值为θ＝0.5。

4.一种根据权利要求1所述基于结构稀疏和熵联合约束的SA-ISAR自聚焦法，其特征在于：Pattern＝{(-1,0),(1,0),(0,-1),(0,1)}，此时Q取值为

5.一种根据权利要求1所述基于结构稀疏和熵联合约束的SA-ISAR自聚焦法，其特征在于：公式(12)中参数取值为a＝1+10^-7，b取值为b＝10^-7。

6.一种根据权利要求1所述基于结构稀疏和熵联合约束的SA-ISAR自聚焦法，其特征在于：S5中，终止条件的判断采用指定迭代次数法。

7.一种根据权利要求6所述基于结构稀疏和熵联合约束的SA-ISAR自聚焦法，其特征在于：指定迭代次数为30次，在达到30次循环后就终止迭代。

8.一种根据权利要求1所述基于结构稀疏和熵联合约束的SA-ISAR自聚焦法，其特征在于：S5中，终止条件的判断采用门限法，当ISAR图像的更新变化小于某个门限值时就停止迭代，用公式表示如下：

其中x^(new)和x^(old)分别表示本次迭代和上一次迭代计算的ISAR图像<X>，门限值Threshold根据实际情况选取。

9.一种根据权利要求8所述基于结构稀疏和熵联合约束的SA-ISAR自聚焦法，其特征在于：门限值Threshold取10^-6。

Images