CN106772365B - 一种基于贝叶斯压缩感知的多径利用穿墙雷达成像方法 - Google Patents

Abstract

本发明公开了一种基于贝叶斯压缩感知的多径利用穿墙雷达成像方法。本发明采用单站单收发天线频率步进雷达系统在平行于墙面多个位置发射脉冲串并获取回波信号，通过随机测量矩阵对回波信号进行压缩采样获得相对少量的测量信号，然后联合依据墙体多径反射特性设计好的字典构建压缩感知模型，使用贝叶斯压缩感知方法将墙后场景信息从测量信号中重构出来，最后利用重构出的信息完成穿墙雷达成像。本发明降低了成本，缩减了记录时间，同时贝叶斯压缩感知技术的使用使得计算速度更快、所需观测数据量更少，结合多径利用方法得到的图像精度更高且有效地抑制了虚影的出现。

Description

一种基于贝叶斯压缩感知的多径利用穿墙雷达成像方法

技术领域

本发明属于穿墙雷达成像技术领域，具体涉及一种通过运用贝叶斯压缩感知(Bayesian Compressive Sensing,BCS)原理并利用墙体多径反射现象处理雷达接收信号从而达到穿墙成像目的的方法。

背景技术

近年来，穿墙雷达成像技术在军、民用领域都收到了越来越广泛的关注。在建筑拆除、灾后救援和城市反恐作战等诸多场景中，都存在着对建筑物内或障碍物后的隐蔽目标探测的难题。传统的穿墙雷达系统往往需要采用多个观测天线作为物理大孔径或者通过移动单个观测天线至多个观测位置合成一个大孔径来获得较大的距离分辨率。这样，系统成本高且记录数据时间长。同时，巨大的原始数据采样量会极大地增加数据传输和存储负担。

针对上述问题，引入近些年兴起的压缩感知(Compressive Sensing,CS)技术是有效的解决办法。根据压缩感知原理，如果一个信号是稀疏的，在远小于奈奎斯特(Nyquist)采样率的条件下，用随机采样获取信号的离散样本，然后通过非线性重建算法就可以重构原信号。一般情况下，墙后场景信息在成像领域通常可被视为绝对稀疏的，因此CS技术能够将墙后场景信息从少量的测量样本中恢复出来。

尽管基于CS的成像方法能够有效地减少数据采样时间、数据量以及节省信号带宽，然而穿墙雷达成像技术的应用场景往往存在低信噪比的特点，并且在墙体多径反射现象的影响下，传统CS技术的成像质量不高，会出现虚假目标(或称“虚影”)。

发明内容

本发明的目的在于提供一种基于贝叶斯压缩感知的多径利用穿墙雷达成像方法，通过少量的测量数据得到精确、高分辨率、低噪声和虚影少的墙后场景图像。

实现本发明目的的技术解决方案为：一种基于贝叶斯压缩感知的多径利用穿墙雷达成像方法，采用单站单收发天线频率步进雷达系统发射脉冲串并获取回波信号，然后通过基于贝叶斯压缩感知的多径利用方法处理回波信号并完成成像，具体实施步骤如下：

步骤1、将穿墙雷达置于初始位置即距墙体外墙面垂直距离d_r处，整个方法共需要M个穿墙雷达的观测位置，该M个观测位置距墙体外墙面垂直距离均为d_r；

步骤2、穿墙雷达的发射天线发射一串频率步进信号，该串频率步进信号包含N个频率步进脉冲，接收天线分别接收各频率下的回波信号z_n,m，n∈{1,2,...,N},m∈{1,2,...,M}；接收天线在第m个穿墙雷达观测位置获取的总回波信号表示为Z_m＝[z_1,m,z_2,m,...,z_N,m]^T；

步骤3、通过高斯随机测量矩阵Φ_m对回波信号Z_m做压缩采样，得到测量向量Y_m进行存储；

步骤4、沿平行于墙壁方向移动穿墙雷达至下一个观测位置，重复步骤2、3，继续移动穿墙雷达，直至获得M组测量信号为止；

步骤5、对成像区域进行网格化分，基于接收信号的相关域信息以及各网格像素与各观测位置间的传播时延构建M个字典Ψ_m；

步骤6、基于M组测量向量、M个字典和M个测量矩阵建立压缩感知模型；

步骤7、根据压缩感知模型，通过贝叶斯压缩感知算法将墙后图像信息矢量s重构，从而完成墙后场景成像。

本发明与现有技术相比，其显著优点为：(1)所需观测天线数量少、数据存储空间小、传输时间短，从而降低了成本，提高了效率；(2)BCS技术的使用使得该方法计算速度更快、所需观测数据量更少；(3)结合多径利用算法，得到的图像精度更高且有效地抑制了虚影的出现。

附图说明

图1为本发明基于贝叶斯压缩感知的多径利用穿墙雷达成像方法示意图。

图2为本发明基于贝叶斯压缩感知的多径利用穿墙雷达成像方法流程图。

图3为本发明基于贝叶斯压缩感知的多径利用穿墙雷达成像方法构成字典时成像区域网格划分示意图。

图4为本发明基于贝叶斯压缩感知的多径利用穿墙雷达成像方法成像多径几何关系模型图，其中，(a)是目标回波直接反射模型图，(b)是目标回波墙体内反射模型图，(c)是目标回波侧边内墙反射模型图。

图5为本发明基于贝叶斯压缩感知的多径利用穿墙雷达成像方法的实例效果图。

具体实施方式

下面结合附图和具体实施例进一步说明本发明基于贝叶斯压缩感知的多径利用穿墙雷达成像方法。

贝叶斯压缩感知算法(Bayesian Compressive Sensing,BCS)是用贝叶斯方法来处理压缩感知的信号重构问题，从统计的角度分析恢复问题，给原始信号赋予一个先验分布，由观测信号、稀疏基字典和观测矩阵恢复出原始稀疏信号的一个后验分布。与传统CS方法只能提供成像场景散射系数的点估计不同，该方法能够提供散射系数完整的后验密度估计，更加适合杂波环境；由于考虑了压缩采样过程中遇到的加性噪声，因而具有更好的鲁棒性；基于BCS方法得到的图像比基于传统CS方法得到的图像相比更加稀疏，且精度更高，计算速度更快，所需观测数量更少。因此，BCS方法相较于传统CS更加适用于穿墙雷达成像实际应用环境。

本发明是基于贝叶斯压缩感知的多径利用穿墙雷达成像方法，采用单站单收发天线频率步进雷达系统发射脉冲串并获取回波信号，然后通过基于贝叶斯压缩感知的多径利用方法处理回波信号并完成成像。本发明应用场景构造如图1所示，按图2所示，具体实施步骤如下：

步骤1、按图1所示，将穿墙雷达置于初始位置即距墙体外墙面垂直距离d_r处，整个方法共需要M个穿墙雷达的观测位置，该M个观测位置距墙体外墙面垂直距离均为d_r；

步骤2、穿墙雷达的发射天线发射一串频率步进信号，该串频率步进信号包含N个频率步进脉冲，接收天线分别接收各频率下的回波信号z_n,m，n∈{1,2,...,N},m∈{1,2,...,M}；接收天线在第m个穿墙雷达观测位置获取的总回波信号表示为Z_m＝[z_1,m,z_2,m,...,z_N,m]^T。

所述频率步进信号的脉冲周期为T_r，载频起始频率为f₁，频率步长为Δf，频率步进数为N，第n个脉冲信号频率为f_n，信号带宽f_B＝f_N-f₁，雷达分辨率即为c/2f_B，c为光在空气中的传播速度。这些参数可根据场景的实际需求进行设置。

所述接收天线接受的回波信号z_n,m包含目标反射墙面多径反射及噪声v_n,m：

其中，目标反射表示为：

其中，P表示场景中存在的目标总数，表示第p个目标的反射系数，τ_m,p表示第m个观测位置到第p个目标的往返传播时延，f_n为频率步进信号中第n个脉冲信号频率；

墙面多径反射表示为：

其中，R表示墙反射路径总数，表示墙反射系数，β_r表示与第r条路径相关的衰变系数，表示第r条反射路径的往返传播时延。

步骤3、通过高斯随机测量矩阵Φ_m对回波信号Z_m做压缩采样，得到测量向量Y_m进行存储，表达式为：

Y_m＝Φ_mZ_m (4)

其中，Φ_m为L×N维的高斯随机矩阵，且L＜＜N，Z_m为N×1的一维列向量，因此Y_m是L×1的一维列向量。由于L＜＜N，Y_m相较于Z_m可以大大节省存储空间并降低传输所需时间。

步骤4、按图1所示，沿平行于墙壁方向移动穿墙雷达至下一个观测位置，重复步骤2、3，继续移动穿墙雷达，直至获得M组测量信号为止；

步骤5、对成像区域进行网格化分，基于接收信号的相关域信息以及各网格像素与各观测位置间的传播时延构建M个字典Ψ_m，具体为：

如图3所示，将墙后感兴趣的成像区域离散化为Q个网格的像素，网格距离单元的大小由成像系统决定，通常被设为雷达分辨率。

考虑到回波信号是目标直接反射和各路径多径反射的叠加，这样，所有图像信息归入一个RQ×1的一维列向量s。该向量中的元素s_q为一个加权指示函数，s_q代表了第q点的反射系数大小，定义为：

R表示墙反射路径总数，则s为RQ×1的一维列向量；

Z_m视为每个像素产生的回波信息叠加，利用每个像素点的特征设计字典Ψ_m，这样Z_m和s之间的关系如下式表示：

Z_m＝Ψ_ms (6)

其中，其中，字典Ψ_m由针对目标直接反射的子字典Ψ_m,t及针对各路径反射的子字典组成，即

Ψ_m,t＝[ψ_m,t(n,q)]是一个N×Q矩阵，其第(n,q)个元素表达式为：

其中，τ_m,q表示第m个位置天线到第q个像素的理论往返传播时延，f_n为频率步进信号中第n个脉冲信号频率；

是一个N×Q矩阵，表示第q个像素沿第r条路径反射的子字典，第(n,q)个元素表达式为：

其中，表示第m个观测位置到第q个像素沿第r条路径反射的理论往返传播时延；

τ_m,q、可根据图4所示穿墙雷达成像多径几何关系模型计算得出。对于图4(a)、(b)所示的模型，τ_m,q可计算为：

其中，l_air1是电磁波从发射天线到前墙面的传播距离，l_air2是电磁波从后墙面第q个像素的传播距离，l_wall是电磁波在墙体内部的传播距离，c和v分别是电磁波在空气中和在墙体内的传播速度。l_air1、l_air2、l_wall与入射角χ_mq、折射角以及第m个收发天线坐标(x_m,-d_r)、第q个像素坐标(x_q,y_q)之间的关系可写成如下表达式：

其中，式(10)中的d_w为墙壁厚度，式(11)中的W为电磁波在墙体内往返的次数。墙体介电常数ε已知，根据Snell定律，有如下表达式：

联立式(9)-(13)即可求得τ_m,q的值。

计算图4(c)所示侧边内墙反射为例传播时延时，可将第q个像素对应侧边内墙做轴对称等效至其虚影，即坐标(2x_w-x_q,y_q)处，然后用等效后的坐标替代原坐标(x_q,y_q)代入上述方法求解方程组即可得到

实际应用中，多径反射子字典的数量可以根据场景需求进行设置。一般情况下，只需考虑前墙体内一次反射和侧边内墙反射，其他路径反射往往由于衰减极大而可以忽略。

步骤6、基于M组测量向量、M个字典和M个测量矩阵建立压缩感知模型，具体为：

按照式(14)方式用M个高斯随机测量矩阵Φ_m构成总测量矩阵Φ，用M个字典Ψ_m构成总字典Ψ，用M个测量向量Y_m构成总测量向量Y，这样根据压缩感知模型Y、Φ、Ψ及s的关系可由式(15)表示：

Y＝ΦΨs＝Θs (15)

其中，Θ为测量矩阵Φ和字典矩阵Ψ的积，称为传感矩阵；通过贝叶斯压缩感知方法求解Y＝Θs获得墙后图像信息s。

步骤7、根据压缩感知模型，通过贝叶斯压缩感知算法将墙后图像信息矢量s重构，从而完成墙后场景成像。

从贝叶斯方法角度考虑压缩感知信号重构问题，图像信息矢量s具有稀疏先验特性，其稀疏度用K表示。通过计算观测信号Y在传感矩阵Θ下得到的重构图像信息矢量权重的后验概率。

Y为总测量向量、Θ为测量矩阵Φ和字典矩阵Ψ的积、s为图像信息矢量、n_e为均值为0、方差为的高斯白噪声，贝叶斯压缩感知模型由式(15)修正重写为：

Y＝Θs+n_e (16)

Y的高斯似然模型为：

其中，K为图像信息矢量s的稀疏度；

至此，压缩感知问题转化为带有s的稀疏约束线性回归问题。Θ已知，则剩下工作是估计系数向量s和噪声方差

利用相对向量机(Relative Vector Machine,RVM)假设先验并估计参数。RVM采用层次先验模型促进稀疏性并满足共轭先验，对于s中的每个元素，假设均值为0的高斯先验：

其中，α_i表示高斯密度函数的精确度(方差的倒数)，α_i先验为

对超参数α进行边缘分布计算得到s的先验为：

其中，符合Student-t分布，调整置信区间(a，b)使该分布在s_i＝0获得峰值，所以此先验促进s的稀疏。

按照上述方法，为选择Γ(α₀|c,d)先验。在假设式(19)先验的情况下，对于给定的Y和Θ，根据贝叶斯公式，s的后验是多变量的高斯分布，s的均值和方差分别为：

μ＝α₀∑Θ^TY (20)

Σ＝(α₀Θ^TΘ+Α)^-1 (21)

其中，Α＝diag(α₁,α₂,…α_RQ)。参数α＝(α₁,α₂,…α_RQ)^T和参数α₀也要进行估计，通过第II类极大似然即Type-II ML(Type-II Maximum-Likelihood)逼近过程，计算α和α₀得到：

其中，μ_i是依据式(20)计算的第i次的s均值，Σ_ii是根据式(21)计算的第i次的s方差。

对s和α，α₀进交替行迭代计算，直到收敛为止。

照此方法计算出矢量s后，按照s对应每个网格以Q为周期进行R次赋值，就能得到墙后场景的成像结果。

实施例1

下面举一实例说明，该实例通过Matlab软件仿真。例子中，仿真参数设置如下：墙壁厚度d_w＝0.15m，介电常数ε＝6。雷达天线距离前墙面0.2m，以0.1m的间隔扫描从0m到5m间的感兴趣成像区域，即共有51个天线位置。发射信号为2～3GHz,1GHz带宽的步进频率信号(201个脉冲信号，步长5MHz)。1GHz带宽提供了0.15m的分辨率。将5m×5m的感兴趣成像区域划分为25×25的网格，网格数为625。考虑目标直接反射、墙体内一次反射和侧边内墙反射三条路径，分别就每条路径计算网格参数构建字典。测量数L设为10。三个目标分别位于(1.7m,0.9m)，(2.5m,4.5m)和(4.1m,2.7m)。

仿真结果如图5所示。由图5可见，三个目标都被准确重构，并且除目标周围有少量可见噪点外，无其他虚影产生。由此可以证明本发明一种基于贝叶斯压缩感知的多径利用穿墙雷达成像方法的可行性和准确性。

本发明有效地解决了传统穿墙雷达SAR成像方法所需观测天线数量多、数据存储空间大、传输时间长的问题，从而降低成本，缩减记录时间。同时，相比于传统压缩感知方法，贝叶斯压缩感知技术的使用使得本发明中的方法计算速度更快、所需观测数据量更少，结合多径利用方法得到的图像精度更高且有效地抑制了虚影的出现。

Claims (5)

1.一种基于贝叶斯压缩感知的多径利用穿墙雷达成像方法，其特征在于，采用单站单收发天线频率步进雷达系统发射脉冲串并获取回波信号，然后通过基于贝叶斯压缩感知的多径利用方法处理回波信号并完成成像，具体实施步骤如下：

步骤1、将穿墙雷达置于初始位置即距墙体外墙面垂直距离d_r处，整个方法共需要M个穿墙雷达的观测位置，该M个观测位置距墙体外墙面垂直距离均为d_r；

步骤2、穿墙雷达的发射天线发射一串频率步进信号，该串频率步进信号包含N个频率步进脉冲，接收天线分别接收各频率下的回波信号z_n,m，n∈{1,2,...,N},m∈{1,2,...,M}；接收天线在第m个穿墙雷达观测位置获取的总回波信号表示为Z_m＝[z_1,m,z_2,m,...,z_N,m]^T；

步骤3、通过高斯随机测量矩阵Φ_m对回波信号Z_m做压缩采样，得到测量向量Y_m进行存储；

步骤4、沿平行于墙壁方向移动穿墙雷达至下一个观测位置，重复步骤2、3，继续移动穿墙雷达，直至获得M组测量信号为止；

步骤5、对成像区域进行网格化分，基于接收信号的相关域信息以及各网格像素与各观测位置间的传播时延构建M个字典Ψ_m；

步骤6、基于M组测量向量、M个字典和M个测量矩阵建立压缩感知模型；

步骤7、根据压缩感知模型，通过贝叶斯压缩感知算法将墙后图像信息矢量s重构，从而完成墙后场景成像；

步骤2中所述频率步进信号的脉冲周期为T_r，载频起始频率为f₁，频率步长为Δf，频率步进数为N，第n个脉冲信号频率为f_n，信号带宽f_B＝f_N-f₁，雷达分辨率即为c/2f_B，c为光在空气中的传播速度；

步骤2中所述接收天线接受的回波信号z_n,m包含目标反射墙面多径反射及噪声v_n,m：

其中，目标反射表示为：

其中，P表示场景中存在的目标总数，表示第p个目标的反射系数，τ_m,p表示第m个观测位置到第p个目标的往返传播时延，f_n为频率步进信号中第n个脉冲信号频率；

墙面多径反射表示为：

其中，R表示墙反射路径总数，表示墙反射系数，β_r表示与第r条路径相关的衰变系数，表示第r条反射路径的往返传播时延。

2.根据权利要求1所述的基于贝叶斯压缩感知的多径利用穿墙雷达成像方法，其特征在于，步骤3中所述通过高斯随机测量矩阵Φ_m对回波信号Z_m做压缩采样，得到测量向量Y_m进行存储，表达式为：

Y_m＝Φ_mZ_m (4)

其中，Φ_m为L×N维的高斯随机矩阵，且L＜＜N，Z_m为N×1的一维列向量，因此Y_m是L×1的一维列向量。

3.根据权利要求1所述的基于贝叶斯压缩感知的多径利用穿墙雷达成像方法，其特征在于，步骤5所述对成像区域进行网格化分，基于接收信号的相关域信息以及各网格像素与各观测位置间的传播时延构建M个字典Ψ_m，具体为：

将墙后成像区域离散化为Q个网格的像素，网格距离单元的大小设为雷达分辨率，这样，所有图像像素信息归入一个一维列向量s，该向量中的元素s_q为一个加权指示函数，s_q代表了第q点的反射系数大小，定义为：

R表示墙反射路径总数，则s为RQ×1的一维列向量；

Z_m视为每个像素产生的回波信息叠加，利用每个像素点的特征设计字典Ψ_m，这样Z_m和s之间的关系如下式表示：

Z_m＝Ψ_ms (6)

其中，字典Ψ_m由针对目标直接反射的子字典Ψ_m,t及针对各路径反射的子字典组成，即

Ψ_m,t＝[ψ_m,t(n,q)]是一个N×Q矩阵，其第(n,q)个元素表达式为：

其中，τ_m,q表示第m个位置天线到第q个像素的理论往返传播时延，f_n为频率步进信号中第n个脉冲信号频率；

是一个N×Q矩阵，表示第q个像素沿第r条路径反射的子字典，第(n,q)个元素表达式为：

其中，表示第m个观测位置到第q个像素沿第r条路径反射的理论往返传播时延；

τ_m,q、根据墙壁特性及穿墙雷达成像多径几何关系计算得出，墙壁的介电常数、厚度分别为ε，d_w。

4.根据权利要求1所述的基于贝叶斯压缩感知的多径利用穿墙雷达成像方法，其特征在于，步骤6所述基于M组测量向量、M个字典和M个测量矩阵建立压缩感知模型，具体为：

按照式(9)方式用M个高斯随机测量矩阵Φ_m构成总测量矩阵Φ，用M个字典Ψ_m构成总字典Ψ，用M个测量向量Y_m构成总测量向量Y，这样压缩感知模型Y、Φ、Ψ及s的关系由式(10)表示：

Y＝ΦΨs＝Θs (10)

其中，Θ为测量矩阵Φ和字典矩阵Ψ的积，称为传感矩阵；通过贝叶斯压缩感知方法求解Y＝Θs获得墙后图像信息s。

5.根据权利要求1所述的基于贝叶斯压缩感知的多径利用穿墙雷达成像方法，其特征在于，步骤7所述根据压缩感知模型，通过贝叶斯压缩感知算法将墙后图像信息矢量s重构，从而完成墙后场景成像，具体为：

Y为总测量向量、Θ为测量矩阵Φ和字典矩阵Ψ的积、s为图像信息矢量、n_e为均值为0、方差为的高斯白噪声，贝叶斯压缩感知模型如下式所示：

Y＝Θs+n_e (11)

Y的高斯似然模型为：

其中，K为图像信息矢量s的稀疏度；

利用相对向量机假设先验并估计参数，相对向量机采用层次先验模型促进稀疏性并满足共轭先验，对于s中的每个元素s_i，假设均值为0的高斯先验：

其中，α_i表示高斯密度函数的精确度，α_i先验为

对超参数α进行边缘分布计算得到s的先验为：

其中，符合Student-t分布，调整置信区间(a，b)使该分布在s_i＝0获得峰值，所以此先验促进s的稀疏；

按照上述方法，为选择Γ(α₀|c,d)先验；在假设式(14)先验的情况下，对于给定的Y和Θ，根据贝叶斯公式，s的后验是多变量的高斯分布，s均值和方差分别为：

μ＝α₀∑Θ^TY (15)

Σ＝(α₀Θ^TΘ+Α)^-1 (16)

其中，Α＝diag(α₁,α₂,…α_RQ)，参数α＝(α₁,α₂,…α_RQ)^T和参数α₀也要进行估计，通过第II类极大似然即Type-II ML逼近过程，计算α和α₀得到：

其中，μ_i是依据式(15)计算的第i次的s均值，Σ_ii是根据式(16)计算的第i次的s方差；

对s和α、α₀进交替行迭代计算，直到收敛为止；

计算出矢量s后，按照s对应每个网格以Q为周期进行R次赋值，就能得到墙后场景的成像结果。