基于变分贝叶斯学习算法的ISAR稀疏频带成像方法与流程

本发明属于雷达技术领域，更进一步涉及稀疏频带高分辨二维isar成像方法，可用于目标形状特征提取与识别。

背景技术：

随着逆合成孔径雷达isar的飞速发展，尽管现有的成像雷达能够提供较高的距离分辨率，但是在对太空碎片、小卫星、航天器等航天目标进行观测时，需要得到更高分辨率的二维雷达图像以对其特征进行准确地描述。目标距离-多普勒图像的距离分辨率由发射信号的带宽决定。获得高距离分辨率的途径有两种：一种是采用超宽带雷达，但是它对硬件系统的要求较高，需要很高的成本；另一种是采用现有的多部雷达在不同频段对目标同时进行观测，之后采用信号处理方法得到较大的等效带宽，最终得到目标的高分辨距离像。通过第二种途径，可以较小的硬件成本获得高质量、距离高分辨的目标图像。因此，研究稀疏频带条件下的高分辨isar成像具有重要的意义。

西北工业大学在其申请的发明专利“机动目标压缩感知isar成像方法”(公开号：cn102841350a，申请号：201210347782.7)中公开了一种机动目标压缩感知isar成像方法。该方法的具体步骤为：对回波数据进行距离压缩、运动补偿及徙动校正得到复矩阵sd，生成高斯随机矩阵φ对sd做降维观测，得到矩阵y；对矩阵y的每一列求解1范数凸优化方程，得到ti时刻的isar成像结果；遍历各个成像时刻，实现对机动目标各时段的isar成像。但是该方法采用的是传统的1范数,稀疏表示能力不足，无法实现低信噪比情况下的稀疏频带聚焦成像。

kevinm.cuomo在其发表的论文“ultrawide-bandcoherentprocessing”(ieeetransactionsonantennasandpropagation,1999,47(6):1094-1107)中提出了一种基于参数化谱估计的稀疏频带雷达目标成像方法。该方法的具体步骤为：对各个子频带雷达回波信号进行全极点模型建模；分别估计各个子频带全极点模型中的模型阶数、频率和振幅；利用得到的全极点模型内插和外推空带数据；采用脉冲压缩方法实现目标的超分辨距离成像。该方法可以获得较大的等效带宽从而实现距离高分辨成像。但是该方法采用的是传统的谱估计方法，对噪声比较敏感，且算法性能依赖于模型阶数估计的准确度，而模型阶数很难准确估计。

技术实现要素：

本发明的目的在于提出一种基于变分贝叶斯学习算法的稀疏频带高分辨成像方法，以实现在稀疏频带和低信噪比情况下对目标的精确成像，获得聚焦良好的二维isar图像。

本发明的基本思路是：基于压缩感知理论，将isar超分辨成像问题转化为稀疏表述问题，采用变分贝叶斯学习算法求解系数向量并进行全频带回波的精确重构，最终实现isar目标的超分辨二维成像，同时，采用多层字典降低算法的复杂度。其实现方案包括如下：

1.基于变分贝叶斯学习的逆合成孔径雷达稀疏频带成像方法，包括：

(1)通过逆合成孔径雷达录取低频带回波s1和高频带回波s2，s1的维数为n1×p，s2的维数为n2×p，其中n1为低频带采样点数，n2为高频带采样点数，p为方位采样点数；

(2)对低频带回波s1和高频带回波s2作预处理，得到低频带方位脉压后的雷达回波s′1和高频带方位脉压后的雷达回波s′2；

(3)将低频带方位脉压后的雷达回波s′1和高频带方位脉压后的雷达回波s′2按行相连，得到已知观测数据矩阵s＝[s′1；s′2]，其维数为n×p，其中n＝n1+n2；

(4)随机生成频点字典矩阵ψ，取ψ中的第n11行到第n12行作为低频带数据对应的字典矩阵ψ′1，取ψ中的第n21行到第n22行作为高频带数据对应的字典矩阵ψ′2，将ψ′1和ψ′2按行相连，得到维数为n×m的观测数据对应的字典矩阵ψ′＝[ψ′1；ψ′2]，其中n11和n12分别为低频带数据在全频带数据中的起始序号和终止序号，n21和n22分别为高频带数据在全频带数据中的起始序号和终止序号，m为字典矩阵ψ的列数；

(5)对字典矩阵ψ′进行删减，获得新字典矩阵φ¹；

(6)计算存在回波的方位单元的全频带回波数据：

(6a)设第一个存在回波的方位单元序号为h1，最后一个存在回波的方位单元序号为h2，令方位单元的初始序号q＝h1，令字典层数k的初始值为k＝1；

(6b)从观测数据矩阵s中取第q个方位单元对应的向量sq，其维数为n×1，设置迭代步数i的初始值设为i＝1，精度的初始值设为m＝1,…,m′，精度矩阵a的初始值设为噪声精度参数β的初始值设为β⁰＝0.01，四个参数初始值设为a⁰＝b⁰＝c⁰＝d⁰＝10-⁴，系数向量ω的初始值ω⁰设为m′×1阶零向量，迭代最大次数iter＝50，终止阈值wth＝10^-3；

(6c)依次计算协方差矩阵σⁱ，均值向量μⁱ，权向量均值ωⁱ，第一参数aⁱ，第二参数bⁱ，第三参数cⁱ，第四参数dⁱ，精度矩阵aⁱ和噪声精度参数βⁱ的估值，其中i为迭代次数；

(6d)判断是否满足终止条件max(ωⁱ-ω^i-1/ω^i-1)＜wth或者满足迭代次数i＞iter，若满足任一个条件，执行步骤(6e)，若不满足，令i＝i+1，返回步骤(6c)；

(6e)根据得到的ωⁱ对字典φ^k进行细化得到新的字典φ^k+1,其中k为字典层数；

(6f)判断是否满足k＝3，若满足，则终止迭代，执行步骤(6g)，若不满足，令k＝k+1，返回步骤(6b)；

(6g)利用式计算当前方位单元的全频带回波数据tq的实部re(tq)和虚部im(tq)，其中re()代表取实部，im()代表取虚部；

(6h)利用公式tq＝re(tq)+jim(tq)计算当前方位单元的全频带数据向量tq；

(6i)更新方位单元序号q＝q+1，当更新后的方位单元序号大于h2时，则停止对方位单元的搜索，获得所有目标对应的全频带回波数据，执行步骤(7)；否则，令字典层数k＝1，返回步骤(6b)；

(7)将所有方位单元的全频带数据向量按列组合，得到全频带数据矩阵s″＝[t1,t2,…,tp]，将全频带数据矩阵s″在距离维作一维傅里叶变换，得到全频带距离-多普勒图像。

本发明具有如下优点：

1.本发明充分利用数据的先验信息，采用变分贝叶斯学习算法求解稀疏表述模型的系数向量，解决了贝叶斯模型的精确求解问题，在低信噪比情况下，可以获得聚焦良好的高分辨二维isar图像。

2.本发明利用多层字典，在提高成像分辨率的同时避免了大字典所产生的高运算复杂度，提高了计算效率。

以下结合附图和具体实施方式对本发明的技术方案作进一步详细描述。

附图说明

图1是本发明的实现流程图；

图2是用本发明对低频带距离多普勒的成像图；

图3是用本发明对高频带距离多普勒的成像图；

图4是用本发明对全频带距离多普勒的成像图；

图5是利用本发明对全频带回波信号重构的距离多普勒成像图。

以下参照附图对本发明的实施例和效果作进一步详细描述。

具体实施方式

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

步骤1，逆合成孔径雷达录取目标的低频带回波s1和高频带回波s2。

逆合成孔径雷达录取目标的低频带和高频带回波，是指两部工作在不同频段的逆合成孔径雷达发射的电磁波在传播过程中遇到目标后，目标对电磁波发生的反射，所反射的回波被雷达接收机接收，并在雷达显示器上显示出目标的低频带回波s1和高频带回波s2。

步骤2，对低频带回波和高频带回波进行预处理，并获得低频带和高频带方位脉压后的信号。

对低频带回波s1和高频带回波s2进行解线频调处理和方位脉压，具体步骤如下：

(2a)将逆合成孔径雷达到场景中心的距离作为参考距离，将与逆合成孔径雷达发射信号载频、调频率相同、距离为参考距离的线性调频信号作为参考信号，将参考信号取共轭后与接收的回波相乘，得到低频带和高频带解线频调处理后的信号：

其中，为距离快时间，tm为方位慢时间，sr1(·)为低频带参考信号，sr2(·)为高频带参考信号，s11为低频带解线频调处理后的信号，s22为高频带解线频调处理后的信号，*表示共轭操作；

(2b)将解线频调处理后的低频带信号s11和高频带信号s22分别在方位维作一维傅里叶变换，得到低频带方位脉压后的雷达回波s1′和高频带方位脉压后的雷达回波s′2。

步骤3，随机生成频点字典矩阵ψ′，并对ψ′进行删减。

(3a)随机生成频点字典矩阵ψ，其维数为q×m，其中q为全频带采样点数，m为字典矩阵ψ的列数，ψ中的每一个元素为q＝0,1,...,q-1为字典行序号，m＝0,1,2,...,m-1为字典列序号；

(3b)取ψ中的第n11行到第n12行作为低频带数据对应的观测矩阵ψ′1，取ψ中的第n21行到第n22行作为高频带数据对应的观测矩阵ψ′2，将ψ′1和ψ′2按行组起来，得到维数为n×m的观测数据对应的字典矩阵ψ′＝[ψ′1；ψ′2]，其中n11和n12分别为低频带数据在全频带数据中的起始序号和终止序号，n21和n22分别为高频带数据在全频带数据中的起始序号和终止序号；

(3c)对字典矩阵ψ′进行删减，获得新字典矩阵φ¹：

(3c1)分别将低频带方位脉压后的雷达回波s1′和高频带方位脉压后的雷达回波s′2在距离维作一维傅里叶变换，得到低频带距离多普勒图像i1和高频带距离多普勒图像i2；

(3c2)设低频带距离-多普勒图像i1中的距离向最低点和最高点对应的序号分别为m11和m12，得到字典矩阵ψ′中的最低点列序号和最高点列序号

(3c3)设高频带距离-多普勒图像i2中的距离向最低点和最高点对应的序号分别为m21和m22，得到字典矩阵ψ′中的最低点列序号和最高点列序号

(3c4)取m11和m21中的最小值记为m1，取m12和m22中的最大值记为m2，保留字典矩阵ψ′的m1列到m2列，删减其余列，得到维数为n×m′的新字典矩阵φ¹，其中m′＝m2-m1。

步骤4，计算存在回波的方位单元的全频带回波数据：

(4a)设第一个存在回波的方位单元序号为h1，最后一个存在回波的方位单元序号为h2，令方位单元的初始序号q＝h1，令字典层数k的初始值为k＝1；

(4b)从观测数据矩阵s中取第q个方位单元对应的向量sq，其维数为n×1，设置迭代步数i的初始值设为i＝1，精度的初始值设为m＝1,…,m′，精度矩阵a的初始值设为噪声精度参数β的初始值设为β⁰＝0.01，四个参数初始值设为a⁰＝b⁰＝c⁰＝d⁰＝10^-4，系数向量ω的初始值ω⁰设为m′×1阶零向量，迭代最大次数iter＝50，终止阈值wth＝10^-3；

(4c)依次计算协方差矩阵σⁱ，均值向量μⁱ，权向量均值ωⁱ，第一参数aⁱ，第二参数bⁱ，第三参数cⁱ，第四参数dⁱ，精度矩阵aⁱ和噪声精度参数βⁱ的估值，其中i为迭代次数，具体步骤如下：

(4c1)根据公式计算第i次迭代的第一参数aⁱ，根据公式计算第i次迭代的第三参数cⁱ；

(4c2)根据公式σⁱ＝(β^i-1[φ^k]^tφ^k+a^i-1)^-1计算第i次迭代的协方差矩阵σⁱ；

(4c3)根据公式μⁱ＝β^i-1σⁱ[φ^k]^tsq计算第i次迭代的均值向量μⁱ，并令第i次迭代的权向量均值ωⁱ＝μⁱ；

(4c4)根据公式计算第i次迭代的第二参数bⁱ，其中为权向量均值ωⁱ的第m个元素，为协方差矩阵σⁱ的第m行第m列对应的元素；根据公式计算第i次迭代的第四参数dⁱ；

(4c5)根据公式计算第i次迭代的精度矩阵aⁱ中的每个元素，根据公式计算第i次迭代的噪声精度参数βⁱ。

(4d)判断是否满足终止条件max(ωⁱ-ω^i-1/ω^i-1)＜wth或者满足迭代次数i＞iter，若满足任一个条件，执行步骤(4e)，若不满足，令i＝i+1，返回步骤(4c)；

(4e)根据得到的ωⁱ中每个元素的值对字典φ^k进行删减和细化，得到新的字典φ^k+1,其中k为字典层数，删减和细化字典的具体步骤如下：

(4e1)设置阈值wth＝wmax×0.01，其中wmax为系数向量ω中元素的最大值，令系数向量ωⁱ元素的第一个序号j＝1；

(4e2)从系数向量ωⁱ中取第j个元素wj，判断是否满足wj＜wth，若满足，将φ^k中的第j列φj进行删除，执行步骤(4e4)；若不满足，执行步骤(4e3)；

(4e3)将φ^k的第j-1列φj-1到第j+1列φj+1进行细化；

(4e4)更新j＝j+1，当更新后的j大于ωⁱ的元素个数时，则停止对j的搜索，获得更新后的字典矩阵φ^k+1；否则，返回步骤(4e2)；

(4f)判断是否满足k＝3，若满足，则终止迭代，执行步骤(4g)，若不满足，令k＝k+1，返回步骤(4b)；

(4g)利用式计算当前方位单元的全频带回波数据tq的实部re(tq)和虚部im(tq)，其中re()代表取实部，im()代表取虚部；

(4h)利用公式tq＝re(tq)+jim(tq)计算当前方位单元的全频带数据向量tq；

(4i)更新方位单元序号q＝q+1，当更新后的方位单元序号大于h2时，则停止对方位单元的搜索，获得所有目标对应的全频带回波数据；否则，令字典层数k＝1，返回步骤(4b)。

步骤5，将所有方位单元的全频带数据向量按列组合，得到全频带数据矩阵s″＝[t1,t2,…,tp]。

步骤6，将全频带数据矩阵s″在距离维作一维傅里叶变换，得到全频带距离多普勒图像。

本发明的效果可通过以下仿真进一步说明：

1.仿真参数

采用工作在c波段和x波段的雷达分别接收低频带和高频带的雷达回波信号，对应载频为6.8ghz和9.2ghz，带宽均为1.2ghz。目标的长度为8.75米，翼展宽度为6米,目标包含330个散射点。已知频带数据占全频带数据的40％，信噪比设为0db。

2.仿真内容

仿真1：对低频带回波信号进行距离多普勒成像，绘制其原始距离多普勒图像，结果如图2。

仿真2：对高频带回波信号进行距离多普勒成像，绘制其原始距离多普勒图像，结果如图3。

仿真3：对全频带回波信号进行距离多普勒成像，绘制其原始距离多普勒图像，结果如图4。

仿真4：利用本发明重构全频带回波信号，绘制其距离多普勒图像，结果如图5。

由图5与图2和图3对比可得，与高低频带距离多普勒图像相比，利用本发明得到的稀疏频带距离多普勒图像信噪比高、虚假点少，分辨率显著提高且图像聚焦良好。图5与图4对比可得，稀疏频带距离多普勒图像和全频带距离多普勒图像基本一致。

仿真结果表明，本发明利用压缩感知理论将isar超分辨成像问题转化为稀疏表述问题，采用变分贝叶斯学习方法求解系数向量并重构全频带回波，充分利用了雷达回波目标的先验信息，在低信噪比情况下得到高分辨率且聚焦效果好的距离多普勒图像。