CN108226928A - 基于期望传播算法的逆合成孔径雷达成像方法 - Google Patents

Abstract

本发明公开了一种基于期望传播算法的逆合成孔径雷达成像方法，本发明的实现步骤是：(1)稀疏贝叶斯建模；(2)利用期望传播算法求解稀疏贝叶斯模型中每个距离单元的实数权向量；(3)重构逆合成孔径雷达ISAR图像散射系数矩阵；(4)转置并画出矩阵，得到回波缺损及低信噪比情况下的目标逆合成孔径雷达ISAR成像结果。本发明基于稀疏信号重构理论实现了逆合成孔径雷达ISAR成像，可用于在回波缺损及低信噪比情况下对空间与空中非合作目标的二维成像。

Description

基于期望传播算法的逆合成孔径雷达成像方法

技术领域

本发明属于通信技术领域，更进一步涉及雷达信号处理技术领域中的一种基于期望传播算法的逆合成孔径雷达ISAR(Inverse Synthetic Aperture Radar)成像方法，可用于逆合成孔径雷达接收到目标回波信噪比低、存在缺损等情况下的空间目标二维成像。

背景技术

由于具有全天时、全天候、高分辨率和远距离等特点，逆合成孔径雷达ISAR在航空与航天目标观测中发挥着重要作用。当受雷达工作模式限制而无法对目标进行连续观测时，会产生方位缺损回波，此时传统基于脉冲压缩的成像方法会产生高旁瓣。此外，当逆合成孔径雷达ISAR对远距离小目标进行探测时，回波信噪比较低，从而很难得到高质量成像结果。

对于逆合成孔径雷达，目标强散射点数目远小于观测样本数，因此具有稀疏性。目前，充分利用其稀疏性，在回波缺损、低信噪比情况下实现空间目标高质量、高分辨成像是提高现有雷达对空间、空中目标探测与监视能力的关键技术。

西安电子科技大学在其申请的发明专利文献“基于稀疏孔径的机动目标逆合成孔径雷达成像方法”(公开号：103901429A，申请号：201410140123.5)中公开了一种基于稀疏重构的成像方法。该方法具体步骤为：对回波数据进行距离压缩和运动补偿，利用正交匹配追踪算法得到稀疏孔径重构回波信号，接着对重构回波信号进行快速傅里叶变换，最终实现逆合成孔径雷达ISAR的距离-多普勒成像。该方法虽然能够在回波缺损情况下实现高分辨成像，但是，该方法仍然存在的不足之处是，正交匹配追踪算法稀疏表征权向量能力不足，在回波缺损及低信噪比情况下容易产生虚假点，无法获得聚焦良好的ISAR图像。

吴称光，邓彬，苏伍各，王宏强，秦玉亮在其发表的论文“基于块稀疏贝叶斯模型的ISAR成像方法”(电子与信息学报2015，37(12)：2941-2947)中提出一种基于块稀疏贝叶斯模型的ISAR成像方法。该方法基于稀疏信号表示理论，将ISAR高分辨成像问题转化为稀疏信号表示问题，利用变分推断算法求解模型参数，最终实现目标逆合成孔径雷达ISAR成像。该方法虽然能够在回波缺损情况下实现二维ISAR成像，但该方法存在的不足之处是，在低信噪比情况下该方法参数估计误差较大，无法获得准确良好的ISAR像。

发明内容

本发明的目的是针对上述稀疏成像技术中数值优化方法稀疏表征能力不足，低信噪比条件下误差较高等局限性，提出一种基于期望传播算法的逆合成孔径雷达ISAR成像方法。该方法充分利用回波信号的稀疏性，通过回波信号的稀疏建模及稀疏系数求解最终实现目标回波缺损与低信噪比情况下的空间目标逆合成孔径雷达ISAR成像。

为实现本发明的目的的具体思路是：基于稀疏信号重构理论，将ISAR成像问题转化为稀疏线性回归问题，利用伯努利高斯模型对目标回波建模，进而采用期望传播算法求解权值向量，最终实现回波缺损及低信噪比情况下的空间目标二维成像。

为实现上述目的，本发明的主要步骤如下：

(1)稀疏贝叶斯建模：

(1a)计算距离向脉冲压缩后每个距离单元有效复数回波向量对应的实数向量；

(1b)将待测目标周围环境中复数噪声向量对应的实数噪声向量的先验分布设置为高斯分布；

(1c)将每个距离单元的实数权向量的稀疏先验分布设置为伯努利高斯分布；

(2)利用期望传播算法计算每个距离单元的实数权向量：

(2a)设置初始距离单元的序号q＝1；

(2b)设置初始迭代次数k＝1，在第q个距离单元中，将先验分布R中所有高斯分布的方差值均初始化为p₀v_s，其中p₀表示方位存在散射点的概率，初始化为5×10^-3，v_s表示散射点幅度变化的方差值，初始化为0.1，将距离单元实数权向量初始化为2N_a×1维零向量，将噪声精度值初始化为4×10^-4；

(2c)利用迭代公式，依次计算当前迭代次数中的五个参数值，用当前迭代次数中的五个参数值分别替换上一次迭代次数中的五个参数值；

(2d)将当前迭代与上次迭代所求得的距离单元实数权向量相减得到差向量，判断差向量中所有元素最大值是否小于10^-5，若是，则停止迭代，得到距离单元对应的实数权向量，执行步骤(2e)，否则，在当前迭代次数上加1后执行步骤(2c)；

(2e)更新距离单元序号q＝q+1，判断更新后的距离单元序号是否大于距离单元总数，若是，则停止对距离单元的搜索，执行步骤(3)；否则，执行步骤(2b)；

(3)重构逆合成孔径雷达ISAR图像散射系数矩阵：

(3a)产生一个维数为2N_a×N_r的矩阵W，其中N_a表示方位向采样总数，N_r表示距离向采样总数，矩阵W的第q列为第q个距离单元的实数权向量；

(3b)利用公式W′＝W₁+jW₂，得到逆合成孔径雷达ISAR图像的散射系数矩阵W′，其中，W₁表示矩阵W第1行到第N_a行，第1列到第N_r列所有元素构成的矩阵，W₂表示矩阵W第N_a+1行到第2N_a行，第1列到第N_r列所有元素构成的矩阵，j表示虚数单位符号；

(4)转置并二维成像：

对逆合成孔径雷达ISAR图像散射系数矩阵W′进行转置操作，将得到的N_r×N_a维矩阵画出来，得到目标的二维逆合成孔径雷达ISAR图像。

本发明与现有的技术相比，具有以下优点：

第一，本发明在稀疏贝叶斯建模时，将待测目标周围环境中复数噪声向量对应的实数噪声向量的先验分布设置为高斯分布，以及将每个距离单元的实数权向量的稀疏先验分布设置为伯努利高斯分布，充分利用了噪声的先验信息和目标散射点分布的稀疏性，克服了现有技术在低信噪比条件下参数估计方法误差较大的问题，使得本发明在低信噪比条件下有效提高了成像质量。

第二，本发明利用期望传播算法，计算每个距离单元的实数权向量，克服了现有技术采用数值优化方法与贝叶斯推断方法稀疏表征权向量能力不足、容易产生虚假点等问题，使得本发明在回波缺损等条件下，能够获得空间目标聚焦良好的逆合成孔径雷达ISAR图像。

附图说明

图1为本发明的流程图；

图2为本发明的仿真图。

具体实施方式

下面结合附图对本发明做进一步的描述。

参照附图1，本发明的具体实施步骤如下：

步骤1，稀疏贝叶斯建模。

A.按照下式，计算距离向脉冲压缩后每个距离单元有效复数回波向量对应的实数向量：

y_f＝Φw_q+n

其中，y_f表示距离向脉冲压缩后第f个距离单元有效复数回波向量对应的实数向量，Φ表示与雷达方位向有效复数傅里叶字典对应的实数傅里叶字典矩阵，w_q表示第q个距离单元与待重构逆合成孔径雷达ISAR图像复数散射系数向量对应的实数权向量，n表示与待测目标周围环境中复数噪声向量对应的实数噪声向量，f与q取值相等。

所述距离向脉冲压缩后每个距离单元有效复数回波向量对应的实数向量由以下步骤得到：

(1)逆合成孔径雷达发射的电磁波在传播过程中遇到目标后，所反射的回波在传播过程中受到遮挡、稀疏观测等影响，最终被雷达接收机接收并记录为目标缺损回波矩阵。

(2)令逆合成孔径雷达到场景中心的距离为参考距离，令与逆合成孔径雷达发射信号载频、调频率相同，距离为参考距离的线性调频信号为参考信号S_ref，将参考信号S_ref取共轭后与有效回波矩阵S_e相乘，得到解线频调后信号S_rd：

其中，为距离快时间，t_m为方位慢时间，S_ref(·)为参考信号，S_rd为解线频调后的信号，*表示共轭操作。

(3)对解线频调后的信号S_rd沿快时间维作傅里叶变换，得到缺损回波距离脉冲压缩后的矩阵。剔除距离脉冲压缩后矩阵中缺损的列向量，对剩余矩阵元素进行转置操作，得到脉冲压缩后的有效复数回波矩阵。

(4)每个距离单元有效复数回波向量实部在上，虚部在下，得到对应的实数向量。

所述的与雷达方位向有效复数傅里叶字典对应的实数傅里叶字典矩阵是由以下步骤得到的：

(1)按照下式，构造一个雷达方位向傅里叶字典：

Φ′＝exp(-j2π·[1:N_a]^TΔt·[-f_r/2:Δf_d:f_r/2])

其中，Φ′表示雷达方位向傅里叶字典，exp表示以自然常数为底的指数操作，j表示虚数单位符号，π表示圆周率，1:N_a表示从1到N_a间隔为1的整数序列，N_a表示方位向采样总数，T表示转置操作，Δt表示时间间隔，定义Δt＝1/f_r，f_r表示雷达脉冲重复频率，-f_r/2:Δf_d:f_r/2表示从-f_r/2到f_r/2间隔为Δf_d的实数序列，Δf_d表示雷达多普勒分辨率，定义Δf_d＝f_r/N_a；

(2)删除雷达方位向傅里叶字典中与待测目标回波矩阵中缺损部分对应的字典向量，对剩余字典向量取转置，得到雷达方位向有效傅里叶字典矩阵Φ_e；

(3)按照构造实数傅里叶字典矩阵Φ，其中Re表示取实部操作，Im表示取虚部操作。

B.将待测目标周围环境中复数噪声向量对应的实数噪声向量的先验分布设置为高斯分布。

实数噪声向量的先验分布由下式给出：

p(n)＝N(n|0,α^-1I)

其中，p(n)表示实数噪声向量的概率密度，α表示噪声精度。

C.将每个距离单元的实数权向量的先验分布设置为伯努利高斯分布。

每个距离单元实数权向量的每个元素的先验分布由下式给出：

p(w_qi)＝z_iN(w_qi|0,v_s)+(1-z_i)δ(w_qi)

其中，p(w_qi)表示第q个距离单元的第i个元素的概率密度，z_i表示非零概率值，v_s表示散射点幅度变化的方差值，N(·)表示高斯分布，δ表示冲击函数。

非零概率值z_i的先验由下式给出：

p(z_i)＝Bernoulli(p₀)

其中，p(z_i)表示非零概率值z_i的概率密度，Bernoulli(·)表示伯努利分布，p₀表示方位单元存在散射点的概率。

步骤2，利用期望传播算法，计算每个距离单元的实数权向量。

A.设置初始距离单元序号q＝1；

B.设置初始迭代次数k＝1，在第q个距离单元中，将先验分布R中所有高斯分布的方差值初始化为p₀v_s，其中p₀表示方位单元存在散射点的概率，初始化为5×10^-3，v_s表示散射点幅度变化的方差值，初始化为0.1，将距离单元实数权向量初始化为2N_a×1维零向量，将噪声精度值初始化为4×10^-4；

C.利用迭代公式，依次计算当前迭代次数中的五个参数值，用当前迭代次数中的五个参数值分别替换上一次迭代次数中的五个参数值；

(1)按照下式，计算似然函数L的每个高斯分布的方差值：

其中，表示第k次迭代时似然函数L的第i₂个高斯分布的方差值，diag表示向量与对角矩阵的转换操作，表示第k-1次迭代时先验分布R所有伯努利-高斯分布方差值构成的方差向量，α表示取值为4×10^-4的噪声精度，Φ表示实数傅里叶字典，T表示转置操作，i₅表示向量中元素的序号，表示的第i₄个元素，向量中元素的序号i₂，i₄，i₅的取值对应相等；

(2)按照下式，计算似然函数L的每个高斯分布的均值：

其中，表示第k次迭代时似然函数L第i₁个高斯分布的均值，表示第k-1次迭代时先验分布R所有伯努利-高斯分布均值构成的均值向量，y表示距离向脉冲压缩后当前距离单元有效复数回波向量对应的实数向量，i₆表示向量中元素的序号，表示的第i₃个元素，向量中元素的序号i₁，i₂，i₃，i₄，i₆的取值对应相等；

(3)按照下式，计算先验分布R中每个伯努利-高斯分布的方差值：

其中，表示第k次迭代时先验分布R第i₄个伯努利-高斯分布的方差值；表示第k次迭代时先验分布R中第i₇个伯努利-高斯分布的非零概率值，v_s表示取值为0.1的散射点幅度变化的方差值，exp表示以自然常数为底的指数操作，向量中元素的序号i，i₁，i₂，i₇的取值对应相等；

所述第k次迭代时先验分布R中第i₇个伯努利-高斯分布的非零概率值由下式计算得到的：

其中，p₀表示取值为5×10^-3的方位单元存在散射点的概率，log表示以自然常数为底的对数操作，向量中元素的序号i₁，i₂，i₇的取值相等；

判断计算得到的是否小于零，如果小于零，则令否则，保持原值；

(4)按照下式，计算先验分布R的每个伯努利-高斯分布的均值：

其中，表示第k次迭代时先验分布R第i₃个伯努利-高斯分布的均值，向量中元素的序号i₂，i₃，i₄的取值对应相等；

(5)按照下式，计算当前距离单元的实数权向量的每个元素：

其中，表示第k次迭代时第q个距离单元实数权向量的第i个元素，向量中元素的序号i，i₁，i₂，i₃，i₄的取值对应相等。

D.将当前迭代与上次迭代所求得的距离单元实数权向量相减得到差向量，判断差向量的所有元素中最大值是否小于10^-5，若是，则停止迭代，得到距离单元对应的实数权向量，执行本步骤的第E步，否则，在当前迭代次数上加1后执行本步骤第C步；

E.更新距离单元序号q＝q+1，判断更新后的距离单元序号是否大于距离单元总数，若是，则停止对距离单元的搜索，执行步骤3；否则，执行本步骤B。

步骤3，重构逆合成孔径雷达ISAR图像散射系数矩阵。

A.产生一个维数为2N_a×N_r的矩阵W，其中N_a表示方位向采样总数，N_r表示距离向采样总数，矩阵W的第q列为第q个距离单元的实数权向量；

B.利用公式W′＝W₁+jW₂，得到逆合成孔径雷达ISAR图像的散射系数矩阵W′，其中，W₁表示矩阵W第1行到第N_a行，第1列到第N_r列所有元素构成的矩阵，W₂表示矩阵W第N_a+1行到第2N_a行，第1列到第N_r列所有元素构成的矩阵，j表示虚数单位符号。

步骤4，转置并二维成像。

对逆合成孔径雷达ISAR图像散射系数矩阵W′进行转置操作，将得到的N_r×N_a维矩阵画出来，得到目标的二维逆合成孔径雷达ISAR图像。

下面结合仿真实验对本发明的效果做进一步的说明。

1.仿真条件：

本发明的仿真实验采用工作在C波段的雷达实测得到的Yak-42飞机数据，对应载频为10GHZ，脉冲重复频率为1000。飞机方位向回波数为512次，距离向采样数为256个，回波数据的缺损率为50％，回波信噪比为-10dB。

2.仿真内容及其结果分析：

本发明的仿真实验有三个，其中：

仿真实验1：对第1～63列、192～319列及448～512列产生缺损的回波数据沿距离向做脉冲压缩，绘制其距离慢时间图，结果如图2(a)所示。

仿真实验2：利用现有技术中的变分推断方法，对距离向脉压后的缺损回波进行重构，绘制其成像结果，结果如图2(b)所示。

仿真实验3：利用本发明对距离脉压后缺损回波进行重构，绘制其成像结果，结果如图2(c)所示。

图2(a)为第1～63列、192～319列及448～512列产生缺损的回波数据沿距离向脉压后的距离慢时间的分布图，图2(a)中的横坐标表示缺损的回波数据沿距离向脉压后的慢时间分布，纵坐标表示缺损的回波数据沿距离向脉压后的距离分布。由图2(a)可以看出，回波数据缺损及低信噪比情况下，图中所显示的距离慢时间图噪声点过多，影响直接进行距离多普勒算法成像的质量。

图2(b)为利用现有技术中的变分推断方法，对缺损回波数据沿距离向脉压后的数据进行重构的成像结果图，图2(b)中的横坐标表示成像结果的方位分布，纵坐标表示成像结果的距离分布。由图2(b)可以看出利用现有变分推断方法得到的成像结果虚假点较多，而且成像结果不准确。

图2(c)为利用本发明所提的方法，对缺损回波数据沿距离向脉压后的数据进行重构的成像结果图，图2(c)中的横坐标表示成像结果的方位分布，纵坐标表示成像结果的距离分布，图2(c)与图2(b)对比可得，利用本发明所得的成像结果能够清楚地呈现飞机目标的几何结构，虚假点少，去噪性能良好。

由上述仿真结果表明，本发明基于稀疏信号重构理论将逆合成孔径雷达ISAR高分辨成像问题转化为稀疏向量重构问题，采用期望传播算法求解实数权向量，充分利用了目标散射点分布的稀疏性与噪声数据的先验信息，在回波缺损及低信噪比情况下能够得到高质量、聚焦良好的逆合成孔径雷达ISAR图像。

Claims (4)

1.一种基于期望传播算法的逆合成孔径雷达成像方法，其特征在于，包括如下步骤：

(1)稀疏贝叶斯建模：

(1a)计算距离向脉冲压缩后每个距离单元有效复数回波向量对应的实数向量；

(1b)将待测目标周围环境中复数噪声向量对应的实数噪声向量的先验分布设置为高斯分布；

(1c)将每个距离单元的实数权向量的先验分布设置为伯努利高斯分布；

(2)利用期望传播算法，计算每个距离单元的实数权向量：

(2a)设置初始距离单元的序号q＝1；

(2b)设置初始迭代次数k＝1，在第q个距离单元中，将先验分布R中所有高斯分布的方差值均初始化为p₀v_s，其中p₀表示方位单元存在散射点的概率，初始化为5×10^-3，v_s表示散射点幅度变化的方差值，初始化为0.1，将距离单元实数权向量初始化为2N_a×1维零向量，将噪声精度值初始化为4×10^-4；

(2c)利用迭代公式，依次计算当前迭代次数中的五个参数值，用当前迭代次数中的五个参数值分别替换上一次迭代次数中的五个参数值；

(2d)将当前迭代与上次迭代所求得的距离单元实数权向量相减得到差向量，判断差向量中所有元素最大值是否小于10^-5，若是，则停止迭代，得到距离单元对应的实数权向量，执行步骤(2e)，否则，在当前迭代次数上加1后执行步骤(2c)；

(2e)更新距离单元序号q＝q+1，判断更新后的距离单元序号是否大于距离单元总数，若是，则停止对距离单元的搜索，执行步骤(3)；否则，执行步骤(2b)；

(3)重构逆合成孔径雷达ISAR图像散射系数矩阵：

(3a)产生一个维数为2N_a×N_r的矩阵W，其中N_a表示方位向采样总数，N_r表示距离向采样总数，矩阵W的第q列为第q个距离单元的实数权向量；

(3b)利用公式W′＝W₁+jW₂，得到逆合成孔径雷达ISAR图像的散射系数矩阵W′，其中，W₁表示矩阵W第1行到第N_a行，第1列到第N_r列所有元素构成的矩阵，W₂表示矩阵W第N_a+1行到第2N_a行，第1列到第N_r列所有元素构成的矩阵，j表示虚数单位符号；

(4)转置并二维成像：

对逆合成孔径雷达ISAR图像散射系数矩阵W′进行转置操作，将得到的N_r×N_a维矩阵画出来，得到目标的二维逆合成孔径雷达ISAR图像。

2.根据权利要求1所述的基于期望传播算法的逆合成孔径雷达成像方法，其特征在于，步骤(1a)中所述的距离向脉冲压缩后每个距离单元有效复数回波向量对应的实数向量是按照下述公式计算得到的：

y_f＝Φw_q+n

其中，y_f表示距离向脉冲压缩后第f个距离单元有效复数回波向量对应的实数向量，Φ表示与雷达方位向有效复数傅里叶字典对应的实数傅里叶字典矩阵，w_q表示第q个距离单元与待重构逆合成孔径雷达ISAR图像复数散射系数向量对应的实数权向量，n表示与待测目标周围环境中复数噪声向量对应的实数噪声向量，f与q取值相等。

3.根据权利要求2所述的基于期望传播算法的逆合成孔径雷达成像方法，其特征在于，所述的与雷达方位向有效复数傅里叶字典对应的实数傅里叶字典矩阵是由以下步骤得到的：

第一步，按照下式，构造一个雷达方位向傅里叶字典：

Φ′＝exp(-j2π·[1:N_a]^TΔt·[-f_r/2:Δf_d:f_r/2])

其中，Φ′表示雷达方位向傅里叶字典，exp表示以自然常数为底的指数操作，j表示虚数单位符号，π表示圆周率，1:N_a表示从1到N_a间隔为1的整数序列，N_a表示方位向采样总数，T表示转置操作，Δt表示时间间隔，定义Δt＝1/f_r，f_r表示雷达脉冲重复频率，-f_r/2:Δf_d:f_r/2表示从-f_r/2到f_r/2间隔为Δf_d的实数序列，Δf_d表示雷达多普勒分辨率，定义Δf_d＝f_r/N_a；

第二步，删除雷达方位向傅里叶字典中与待测目标回波矩阵中缺损部分对应的字典向量，对剩余字典向量取转置，得到雷达方位向有效傅里叶字典矩阵Φ_e；

第三步，按照构造实数傅里叶字典矩阵Φ，其中Re表示取实部操作，Im表示取虚部操作。

4.根据权利要求2所述的基于期望传播算法的逆合成孔径雷达成像方法，其特征在于，步骤(2c)中所述利用迭代公式，依次计算当前迭代次数中的五个参数值的具体步骤如下：

第一步，按照下式，计算似然函数L的每个高斯分布的方差值：

其中，表示第k次迭代时似然函数L的第i₂个高斯分布的方差值，diag表示向量与对角矩阵的转换操作，表示第k-1次迭代时先验分布R所有伯努利-高斯分布方差值构成的方差向量，α表示取值为4×10^-4的噪声精度，Φ表示实数傅里叶字典，T表示转置操作，i₅表示向量中元素的序号，表示的第i₄个元素，向量中元素的序号i₂，i₄，i₅的取值对应相等；

第二步，按照下式，计算似然函数L的每个高斯分布的均值：

其中，表示第k次迭代时似然函数L第i₁个高斯分布的均值，表示第k-1次迭代时先验分布R所有伯努利-高斯分布均值构成的均值向量，y表示距离向脉冲压缩后当前距离单元有效复数回波向量对应的实数向量，i₆表示向量中元素的序号，表示的第i₃个元素，向量中元素的序号i₁，i₂，i₃，i₄，i₆的取值对应相等；

第三步，按照下式，计算先验分布R中每个伯努利-高斯分布的方差值：

其中，表示第k次迭代时先验分布R第i₄个伯努利-高斯分布的方差值；表示第k次迭代时先验分布R中第i₇个伯努利-高斯分布的非零概率值，v_s表示取值为0.1的散射点幅度变化的方差值，exp表示以自然常数为底的指数操作，向量中元素的序号i，i₁，i₂，i₇的取值对应相等；

所述第k次迭代时先验分布R中第i₇个伯努利-高斯分布的非零概率值由下式计算得到的：

其中，p₀表示取值为5×10^-3的方位单元存在散射点的概率，log表示以自然常数为底的对数操作，向量中元素的序号i₁，i₂，i₇的取值相等；

判断计算得到的是否小于零，如果小于零，则令否则，保持原值；

第四步，按照下式，计算先验分布R的每个伯努利-高斯分布的均值：

其中，表示第k次迭代时先验分布R第i₃个伯努利-高斯分布的均值，向量中元素的序号i₂，i₃，i₄的取值对应相等；

第五步，按照下式，计算当前距离单元的实数权向量的每个元素：

其中，表示第k次迭代时第q个距离单元实数权向量的第i个元素，向量中元素的序号i，i₁，i₂，i₃，i₄的取值对应相等。