2018, Vol. 39
针对频率步进连续波探地雷达[1-2](stepped-frequency continuous wave ground penetrating radar, SFCW-GPR)探测的块目标物体, 在探测区域不能满足稀疏特性的问题, 本文提出利用正交稀疏基[3-4]的特性来针对性地对块状目标进行稀疏化处理, 使其满足稀疏性条件; 然后通过字典矩阵和稀疏矩阵结合形成适用于块目标物体的新的观测矩阵; 通过压缩感知凸优化算法求解稀疏化系数, 接着把求解得到的稀疏化系数通过稀疏变换得到块目标的反射系数向量; 最后把块目标的反射系数按列重新进行排列得到目标的二维重构图像.
1 SFCW-GPR工作原理探地雷达在地面上沿着目标的一个切面水平移动, 首先在一个测量位置发射频率从f0到fH步进的电磁波信号(f0为起始频率, fH为最高频率), 接收天线接收，经过相应目标反射后的回波信号; 然后收发天线移动到下一个测量位置重新发射频率从f0到fH步进的电磁波信号, 接收天线再次接收，直到探测完整个探测区域.频率步进连续波探地雷达探测目标的工作示意图见图 1.
|
图 1 雷达目标探测示意图 Fig.1 Working schematic diagram of the GPR |
通过对图像进行稀疏变换, 可以使用少量的数据来准确地表示原图像的信息.本文利用离散余弦变换(DCT: discrete cosine transform)和离散小波变换[5](DWT: discrete wavelet transform)这两种稀疏正交基来实现对块目标物体的稀疏化处理.
2.2 块目标的稀疏化表示在雷达探测块目标的过程中, 由于目标物体不满足稀疏特性, 所以对其探测区域的网格进行拉伸延长划分后得到的反射系数向量也不可能满足稀疏特性, 如图 2所示.
|
图 2 块目标网格区域拉伸图 Fig.2 Mesh area graph of the block objects |
首先对网格区域内的块目标按列进行拉伸后组成反射系数向量, 然后通过正交稀疏基对反射系数进行正交投影变换, 在此设定稀疏正交基为ρ, 投影后的系数为α, 具体稀疏变换式为
|
(5) |
频率步进雷达工作频率是从f0开始, 以步进频率Δf呈阶梯式上升, 频率步进点数为N(N为在一个发射周期中电磁波发射的总频率点数), 则系统的第n个发射频率为fn=f0+nΔf, 其中n=0, 1, …, N-1;模拟测量过程中共有M个天线位置, 忽略地下媒质的衰减效应和地表的强反射效应, 则系统在第m(m=0, 1, …, M-1)个天线位置第n个工作频点的接收信号数学表达式为
|
(6) |
其中:σp是第p个散射目标对应的反射系数;fn是第n个发射频率;τj, m为电磁波在第m个天线测量位置与第j个点目标的双程时延.然后将目标成像区域分成K×L个均匀大小的网格.每个网格对应一个目标反射系数γ(k, l), 将这些反射系数纵向拉伸可以得到一个KL×1维的向量, 记为γ.系统接收到的回波信号可以表示为一个字典矩阵和反射系数相乘的形式:
|
(7) |
其中:sm=[sm(0),sm(1)，…，sm(N－1)]为N×1维测量数据向量;sm(n)为当工作频率为fn时, 第m个天线测量位置采集到的回波数据; Ψm为N×KL维的字典矩阵, 其第j列为
|
(8) |
其中τj, m为第j个网格目标与第m个天线测量位置的双程时延.双程时延采用一种近似算法来确定电磁波在空气与地下媒介间的折射点位置[6-7].
模拟雷达系统对应的过完备字典矩阵Ψ相应为MN×KL维的矩阵.系统在全部的M个天线测量位置接收到的总数据S和目标区域划分网格的反射率向量γ之间的关系如式(9)所示:
|
(9) |
结合式(5), 由于稀疏系数α具有稀疏特性, 相应的块目标压缩感知成像模型为
|
(10) |
其中
|
(11) |
然后依据凸优化工具箱[9]就能转化为可以解决的二次规划问题，进而求解该成像模型的稀疏系数向量α; 接着对稀疏系数向量α进行正交基逆变换就可以得到相应块目标的反射系数向量γ; 最后将目标的反射系数向量γ按网格区域大小进行逐列排布即可得到最终的块目标二维图像.
3 数据仿真及分析为了验证上述算法的有效性, 本节通过Matlab对块目标压缩感知的频率步进连续波探地雷达图像重构算法进行模拟仿真.系统实验参数设置如下:模拟系统发射电磁波的起始频率为f0=1 GHz, 最高频率为fH=3 GHz, 频率步进点数N=101, 频率步进间隔Δf=(fH－f0)(N－1)= 20 MHz, 天线测量位置总数M=30;雷达距离地面的距离0.1 m, 暂时不考虑地面强回波的影响, 成像区域沿方位向从1 cm到32 cm; 沿距离向从21 cm到52 cm.成像区域被均匀划分成32×32个网格目标.
首先针对单个块目标进行验证, 设定块目标的水平位置从4 cm到24 cm, 垂直距离从30 cm到46 cm, 采样矩阵的设计是从雷达天线水平移动的30个位置随机选取15个, 对应位置的101个发射频率中选取20个频率数据, 采样数据约占总数据量的(15×20)/(30×101)=9.9%.为使成像效果更明显, 对其重构成像图进行归一化处理, 动态显示的范围为[-15, 0], 下面运用不同的方法针对单个块目标的压缩感知重构算法进行仿真对比分析.
从图 3可以看出, 未稀疏化处理的压缩感知算法对块目标进行重建缺失严重, 误差大, 分辨率低; 但是通过系数稀疏变换后, 运用块目标的压缩感知成像模型来对其进行稀疏化重构, 基本上能准确重构出目标物体的位置及形状等信息.
|
图 3 单个块目标的压缩感知重构算法成像对比 Fig.3 Contrast of compressed sensing reconstruction algorithm for single block object (a)—原始单个块目标; (b)—未稀疏化CS算法; (c)—DCT稀疏化的CS算法; (d)—DCT稀疏化的CS算法+高通滤波; (e)—DWT稀疏化的CS算法; (f)—DWT稀疏化的CS算法+高通滤波. |
相对重构误差(relative reconstruction error)：
|
其中:x为原始目标物体散射系数向量; 
从表 1可以看出, 未稀疏化的L-OPT压缩感知算法相对重构误差最大, 分辨率最低, 而基于稀疏化处理的L-OPT压缩感知算法对块目标的重构误差最小, 精度最高.
|
表 1 单个块目标不同稀疏化下重建图像相对重构误差对比 Table 1 Comparison of MSE under different sparsity for single block object |
为了进一步验证该算法的可行性, 针对多个块目标物体进行仿真分析, 成像区域不变, 设定多个大小不同的块目标.目标物体1水平位置从4 cm到28 cm, 垂直位置从24 cm到31 cm; 目标物体2水平位置从4 cm到16 cm, 垂直位置从36 cm到47 cm; 目标物体3水平位置从21 cm到30 cm, 垂直位置从40 cm到47 cm.采样矩阵的选取同单个块目标实验, 采样数据约占总数据量的(15×20)/(30×101)=9.9%.为使成像效果更明显, 对其重构成像图进行归一化处理, 动态显示的范围为[-15, 0], 图 4为运用不同的方法针对块目标的压缩感知重构算法进行仿真对比分析.
|
图 4 多个块目标的压缩感知重构算法成像对比 Fig.4 Contrast of compressed sensing reconstruction algorithm for multiple block objects (a)—原始多个块目标; (b)—未稀疏化的L-OPT; (c)—DCT稀疏化+L-OPT; (d)—haar小波稀疏化+L-OPT+高通滤波; (e)—haar小波稀疏化+L-OPT; (f)—haar小波稀疏化+L-OPT+高通滤波. |
从图 4可以看出, 本算法对于多个目标物体同样适用, 都可以较为准确地分辨多个物体的位置大小信息, 下面通过重构相对误差来对其进行定量分析.
从表 2可知, 基于稀疏化的块目标压缩感知成像算法拥有比常规压缩感知算法更高的重构分辨率, 更快的重构速度, 相对重构误差更小的优势.
|
|
表 2 多个块目标不同稀疏化下重建图像相对重构误差对比 Table 2 Comparison of MSE under different sparsity for multiple block objects |
频率步进连续波技术具有很多优点, 所以在探地雷达领域的应用也越来越重视.针对块目标物体在成像区域不满足稀疏特性现象, 引入了正交稀疏基来对其进行稀疏化处理, 通过正交稀疏基变换能把目标物体的大部分信息投影到某些离散的点上, 即可满足稀疏性的条件.鉴于此构建一种新的适用于块目标压缩感知重构模型, 在该模型中把目标物体稀疏化后的系数作为待重构的压缩感知系数, 并通过压缩感知的凸优化重构算法求得稀疏系数；然后再由正交稀疏变换来得到重构的块目标图像.通过仿真验证该模型比未稀疏化处理压缩感知重构图像具有更高的精度和分辨率.
| [1] |
Couts T, Gurbuz A C.
Multistatic ground-penetrating radar experiments[J]. IEEE Transactions on Geoscience and Remote Sensing, 2007, 45(8): 2544–2553.
DOI:10.1109/TGRS.2007.900677 |
| [2] |
Zhuge X, Yarovoy A G.
Modified Kirchhoff migration for UWB MIMO array-based radar imaging[J]. IEEE Transactions on Geoscience and Remote Sensing, 2010, 48(6): 2692–2703.
DOI:10.1109/TGRS.2010.2040747 |
| [3] |
Donoho D L.
Compressed sensing[J]. IEEE International Conference on Information Theory, 2006, 52(4): 1289–1306.
|
| [4] |
Candes E J, Wakin M B.
An introduction to compressive sampling[J]. IEEE Signal Processing Magazine, 2008, 25(2): 21–30.
DOI:10.1109/MSP.2007.914731 |
| [5] |
Al-Haj A.
Combined DWT-DCT digital image watermarking[J]. Journal of Computer Science, 2007, 3(9): 740–746.
DOI:10.3844/jcssp.2007.740.746 |
| [6] |
Johansson E M, Mast J E. Three-dimensional ground-penetrating radar imaging using synthetic aperture time-domain focusing[C]//Spies International Symposium on Optics. San Diego, 1994: 205-214.
|
| [7] |
屈乐乐, 方广有, 杨天虹.
压缩感知理论在频率步进探地雷达偏移成像中的应用[J]. 电子与信息学报, 2011, 33(1): 21–26.
( Qu Le-le, Fang Guang-you, Yang Tian-hong. The application of compressed sensing to stepped-frequency ground penetrating radar migration imaging[J]. Journal of Electronics & Information Technology, 2011, 33(1): 21–26. ) |
| [8] |
刘福来, 彭泸, 汪晋宽, 等.
基于加权L1范数的CS-DOA算法[J]. 东北大学学报(自然科学版), 2013, 34(5): 655–657.
( Liu Fu-lai, Peng Lu, Wang Jin-kuan, et al. CS-DOA algorithm based on weighted L1 norm[J]. Journal of Northeastern University (Natural Science), 2013, 34(5): 655–657. ) |
| [9] |
Boyd S, Vandenberghe L.
Convex optimization[M]. Oxford: Cambridge University Press, 2004.
|