2. 沈阳航空航天大学 电子信息工程学院, 辽宁 沈阳 110136
2. School of Electronic and Information Engineering, Shenyang Aerospace University, Shenyang 110136, China
由于频率步进探地雷达(SFGPR, stepped frequency ground penetrating radar)具有较高的动态范围和较强的射频抗干扰能力等优点, 近年来在探地雷达领域受到越来越多的关注[1-3].目前, 探地雷达不断朝着多通道、多极化、多波段和高分辨率方向发展, 使得SFGPR存在数据采集量大和数据采集时间长等缺点.近年来, 美国学者Donoho等提出了压缩感知(CS, compressive sensing)理论[4].目前, 国内外学者将CS理论应用到SFGPR的系统构建和成像重建等领域[5-7].
但在SFGPR的实际测量过程中, 由于强地面杂波的存在, CS成像算法的重建性能会急剧下降[8].所以需采用合适的杂波抑制方法对原始频域采样数据进行预处理以去除较强的地面反射波.传统的SFGPR杂波抑制算法大都是基于均匀采样数据[9], 而在CS成像过程中, 考虑到成像数据的大幅度减少且在空频域表现为非均匀采样数据, 经典杂波抑制算法往往都不适用.针对上述问题, 本文提出了一种基于子空间投影杂波抑制的压缩感知SFGPR成像重建算法, 该算法首先在每个天线测量位置基于接收信号的稀疏性重建所有频点的信号, 然后通过经典的子空间投影杂波抑制算法去除地面强反射波, 最后基于CS成像模型对地下目标进行高分辨成像重建, 从而实现在强杂波环境对地下目标进行准确成像定位.
1 SFGPR信号模型与压缩采样对于传统SFGPR系统而言, 设系统的初始频率为f0, 频率步长为Δf, 带宽B=(N-1)Δf.在发射带宽B内, 第n个频点的发射频率为
(1) |
设测量过程中共有M个天线测量位置, 如果忽略地下媒质的衰减效应, 则系统在第m(m=0, 1, …, M-1) 个天线位置第n个工作频点的接收信号可表示为
(2) |
其中:σs是对应空气与地面分界面的反射率; τs为系统到地面的双程传输延时; P为地下点目标的个数; σp(σp≪σs)是第p个点目标对应的反射率; τp, m为在第m个测量位置系统到第p个点目标的双程传输延时.
为得到地下待成像区域的2维反射率分布r(k, l), 其中k=0, 1, …, Nx-1;l=0, 1, …, Ny-1, Nx和Ny分别表示水平和深度方向划分的网格数, 首先可将r(k, l)通过列堆叠转换成NxNy×1维反射率向量r, 然后根据r的估计结果通过重排即可得到r(k, l).这样, 式(2) 可表示成矩阵形式:
(3) |
其中:
(4) |
根据M个测量孔径可得到MN×1维复合测量数据向量
(5) |
由于P往往远远小于网格数NxNy, 所以r具有较强的稀疏性.根据CS理论, 可以用Q1Q2
(6) |
其中Φ为Q1Q2×MN维测量矩阵, 考虑到在空频域实现联合压缩采样, 则Φ可表示为
(7) |
其中:⊗表示Kronecker积;IQ1为Q1×Q1维单位矩阵;Φs表示Q2×M维空域测量矩阵, 其构造可从M×M维单位矩阵中随机选取Q2行得到.Φm(m=0, 1, …, M-1) 表示Q1×N维频域测量矩阵, 其构造可从N×N维单位矩阵中随机选取Q1行得到.
图 1给出了传统SFGPR和压缩感知SFGPR的空频域采样示意图, 可以看出, 与传统SFGPR相比, 压缩感知SFGPR的空频域采样数据将大幅减少.对于给定的压缩采样测量数据向量
(8) |
式(8) 表示的稀疏最优化问题可利用正交匹配追踪算法得到求解[10].
2 压缩采样子空间投影杂波抑制 2.1 单个天线测量位置的频域数据重建在CS测量模型下, 由于测量数据大幅度减少且表现为空频域随机采样数据, 使得传统杂波抑制算法不再适用.为了能够使用经典杂波抑制算法, 本文提出在每个天线测量位置基于CS测量数据重建原始发射频点的均匀采样数据.首先将双程探测时窗τmax均匀划分为L个时延网格, 则可以得到L×1维双程传输时延向量τ=[τ1, τ2, …, τL]T.则在第m个测量位置, 压缩采样频域数据向量
(9) |
其中A是N×L维字典矩阵, 其第l列可表示为
(10) |
bm是L×1维回波幅度向量, 其第l个元素对应第l个时延网格的幅度, bm表示为
(11) |
在SFGPR实际探测过程中, 由于目标回波个数通常远小于划分的时延网格数, 所以bm是稀疏向量. bm可通过求解下面最优化问题得到:
(12) |
式(12) 可通过正交匹配追踪算法求解得到bm, 从而得到:
(13) |
当所有的天线位置处频域均匀采样数据得到重建后, 即可使用合适的杂波抑制方法去除地面强反射波.
2.2 子空间投影杂波抑制当在每个测量位置恢复N个频点测量信号后, 假设共有Q2个天线测量位置, 那么对应Q2个天线位置N个发射频点的所有测量数据可构成N×Q2维矩阵E.由于测量数据可分解为地面回波子空间和地下目标回波子空间[11], 所以可对矩阵E进行奇异值分解:
(14) |
其中H代表共轭转置, 正交矩阵U和V分别包含左右奇异向量, Λ为对角矩阵, 可表示为
(15) |
其中:λ1, λ2, …, λQ2为奇异值, 奇异值按λ1≥λ2≥…≥λQ2降序排列.
由于地面反射波幅度比目标回波幅度强很多, 通常可以用E的前K个较大的奇异值对应的奇异向量来构造地面回波子空间Aground.
(16) |
其中ui和vi分别代表前K个较大奇异值对应的左、右奇异向量.
地面回波子空间的正交子空间Aground⊥为
(17) |
其中I是单位矩阵.将E投影到Aground⊥上即可得到杂波抑制后的空频域采样数据:
(18) |
至此就能得到地面杂波抑制后的目标回波空频域测量数据.然后再采用CS成像算法对地下目标进行成像重建.
3 实测数据处理利用本文所提算法对佐治亚理工大学公布的SFGPR系统实测数据[1]进行处理, 系统工作频段为60 MHz~8.06 GHz, 频率步进间隔为20 MHz, 在每个天线测量位置共测得401个频点采样数据.实验在沙箱里进行, 收发天线相位中心距离地面的高度为27.8 cm, 收发天线间隔为12 cm, 移动步长为2 cm, 两个直径为12.7 cm的金属球埋入沙坑中, 沙的介电常数为4.两个金属球在地下的真实位置如图 2所示.选取天线位置从-60 cm到60 cm之间的61道数据, 对应61×401=24 461个空频域均匀采样数据.
把待成像区域划分成61×31个空间网格, 每个网格尺寸为1 cm×1 cm.对应CS算法, 随机抽取30个天线位置, 每个天线位置抽取80个频点采样数据, 共有2 400个空频域随机采样数据, 只占传统SFGPR系统数据采集量的10.1%.图 3给出了未采用杂波抑制技术预处理的后向投影(back projection, BP)成像重建结果和CS成像重建结果.为了便于比较, 成像结果均采用归一化dB值显示, 显示动态范围为40 dB.可以看出在有地面强反射波情况下, CS成像算法已无法重建准确的目标像.
图 4给出了经子空间投影杂波抑制预处理后的BP成像结果和CS成像结果.BP成像算法中直接采用子空间投影杂波抑制算法去除地面强反射波, 而CS成像算法采用本文方法去除地面强反射波.采用的计算环境为Intel(R) Core(TM) i7-4510U CPU @ 2.0 GHz, 4 GB内存, BP成像结果的重建时间为203.8 s, CS成像结果的重建时间为23.8 s.从图 4可以看出, 尽管两种成像算法都能得到准确的目标像, 但这时CS成像算法不仅所需数据量少, 而且成像结果的目标旁瓣较小, 分辨率高, 更利于对地下目标的辨识.
由于CS成像重建算法在强杂波测量环境中会失效, 而在CS测量模式下测量数据表现为空频域随机采样数据, 使得传统的杂波抑制算法不再适用.本文提出一种基于子空间投影杂波抑制的压缩感知SFGPR成像重建算法, 该算法首先在每个天线测量位置重建所有频域均匀采样数据, 然后基于正交子空间投影技术得到目标回波子空间, 最后利用CS成像模型重建高质量的地下目标像.实测数据处理结果表明所提成像重建算法能够实现在强杂波环境对地下目标进行准确成像定位, 而且成像所需数据量少, 成像结果具有旁瓣低和分辨率高等优点.
[1] | Counts T, Gurbuz A C, Scott W R. 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] |
屈乐乐, 方广有, 杨天虹.
频率步进探地雷达系统小型化的设计与实现[J]. 仪器仪表学报, 2010, 31(3): 699–703.
( Qu Le-le, Fang Guang-you, Yang Tian-hong. Design and implementation of miniature stepped-frequency ground penetrating radar system[J]. Chinese Journal of Scientific Instrument, 2010, 31(3): 699–703. ) |
[3] | Nicolaescu I. Improvement of stepped-frequency continuous wave ground-penetrating radar cross-range resolution[J]. IEEE Transactions on Geoscience and Remote Sensing, 2013, 51(1): 85–92. DOI:10.1109/TGRS.2012.2198069 |
[4] | Donoho D. Compressed sensing[J]. IEEE Transactions on Information Theory, 2006, 52(4): 1289–1306. DOI:10.1109/TIT.2006.871582 |
[5] | Yang J, Jin T, Huang X, et al. Sparse MIMO array forward-looking GPR imaging based on compressed sensing in clutter environment[J]. IEEE Transactions on Geoscience and Remote Sensing, 2014, 52(7): 4480–4494. DOI:10.1109/TGRS.2013.2282308 |
[6] | Krueger K R, McClellan J H, Scott W R. Efficient algorithm design for GPR imaging of landmines[J]. IEEE Transactions on Geoscience and Remote Sensing, 2015, 53(7): 4010–4021. DOI:10.1109/TGRS.2015.2388786 |
[7] | Ambrosanio M, Pascazio V. A compressive-sensing-based approach for the detection and characterization of buried objects[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2015, 8(7): 3386–3395. DOI:10.1109/JSTARS.2015.2421812 |
[8] | Qu L L, Yang T H. Investigation of air/ground reflection and antenna beam width for compressive sensing SFCW GPR migration imaging[J]. IEEE Transactions on Geoscience and Remote Sensing, 2012, 50(8): 3143–3149. DOI:10.1109/TGRS.2011.2179049 |
[9] | Wei X, Zhang Y. Interference removal for autofocusing of GPR data from RC bridge decks[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2015, 8(3): 1145–1151. DOI:10.1109/JSTARS.2015.2402211 |
[10] | Tropp J A, Gilbert A C. Signal recovery from random measurements via orthogonal matching pursuit[J]. IEEE Transactions on Information Theory, 2007, 53(12): 4655–4666. DOI:10.1109/TIT.2007.909108 |
[11] | Solimene R, Cuccaro A, Aversano A, et al. Ground clutter removal in GPR surveys[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2014, 7(3): 792–1151. DOI:10.1109/JSTARS.2013.2287016 |