论文部分内容阅读
摘要：本文提出了一种基于原子集的压缩感知方法，将其应用到探地雷达成像信号处理上。首先依据理想目标点逆散射模型构造适用于探地雷达的原子集；其次利用构造原子的特性，降低图像的维数，优化重构算法参数，提高成像质量。实验结果表明，基于原子集的探地雷达成像算法具有数据采集时间短、成像速度快、算法快速稳定的特点。
关键词：探地雷达；压缩感知；成像算法；原子集
中图分类号：C35文献标识码： A
1、引言
探地雷达(GPR)通过发射电磁波，利用电磁波在媒质电磁特性不连续处产生的反射和散射实现浅表层成像进而辨识地表电磁特性的变化，实现对地下目标的无损探测。[1]近年来，探地雷达技术的快速发展得益于硬件电子工艺技术和数字信号处理技术的发展，而理论上没有质的飞跃。探地雷达的有效应用，故而离不开硬件系统和信号处理技术，但是如何引入新的理论、新的探地雷达成像算法和特征检测算法，是探地雷达技术发展的一个方向。现有的探地雷达遵循Nyquist-Shannon采样理论，常用成像算法有距离多普勒成像算法(RD)[2]、时域反向投影成像算法(Back Projection, BP)[3]、基于波场延拓的距离偏移成像算法[4]、逆时偏移成像算法[5]和Chirp Scaling成像算法[6]等。探地雷达具有成像速度快、空间分辨率高的特点，这就要求探地雷达具备复杂的匹配滤波器、高密度的采样孔径以及长时间的单孔径测量，最终得到庞大的数据量。
近年来，由Donoho, Candes, Romberg 和Tao等人[7-10]提出的压缩感知(CS)理论是一种充分利用信号稀疏性或可压缩性的全新信号获取和重构理论。压缩感知理论指出，只要信号是可压缩的或在某个变换域是稀疏的，就可以用一个与变换基不相关的观测基将变换所得的高维信号投影到一个低维空间上，最后通过求解一个优化问题就可以从低维空间的少量投影中以高概率重构原始信号。压缩感知理论是对传统信息理论的发展，突破了Nyquist-Shannon采样定理的极限，是科研领域新的研究热点。
本文将压缩感知理论应用到探地雷达成像上，根据成像区域目标稀疏性的先验性知识，构造合适的原子集对探测区域进行压缩测量从而降低采样孔径密度，缩短单孔径测量时间、减少获取数据量。仿真和实验结果表明，采用压缩感知探地雷达成像算法，可获得较传统探地雷达成像算法更好的成像效果。
2、压缩感知理论概述
若有一个k稀疏的长为N的实离散时间信号x(n)，n=1,2，…n。设x为一个N维列向量，即x∈RN×1，如果信号x是稀疏的或可压缩的，那么存在一组稀疏基可将x表示为
，
α和x均为N×1的矩阵，ψ为N×N的矩阵。
将此稀疏信号投影到一组测量基上，得到M(M 将稀疏基ψ与测量基φ结合起来，那么Y=Φx=Φψα=ʘα，其中ʘ=Φψ是M×N的矩阵，ʘ∈RM×N。那么由维的测量值Y来恢复维数为N(M （1）
来重构原始信号，然而零范數约束优化重构算法是NP（Non-deterministic Polynomial）问题，可通过求解l1范数约束次最优化来解决，即通过
（2）
来从M个线性测量中恢复出α。
3、压缩感知探地雷达成像算法
3.1理想点散射模型
从数据分析的角度看，探地雷达成像问题即信号的表示问题，由探地雷达成像理想点散射模型描述可知，对于含有P个理想点散射中心的目标而言，目标后向散射场表示为[12]
（3）
其中，f为频率，θ为方位角，c为光速，Ak，xk，yk分别表示第k个散射中心的复振幅和其在散射场中所处的位置坐标。如果将上述表达式（3）中的下面这一项作为一组基函数，估计目标散射幅度参数的问题，即可看成是有限的观测信号在特定基下的表示问题。
3.2 GPR原子集的构造
对于GPR图像而言，可基于目标散射模型来构造与之相适应的原子集。将成像区域划分为M×N大小的网格，目标后向散射场表示为:
（4）
在上式中，τ1=0,1，…M-1；τ2=0,1，…N-1。
设在M×N网格上每个点都是一个散射中心，若有目标，则其为真正的散射中心，系数A(k1，k2)不为零，否则，系数A(k1，k2)为零。
因此，此时的参数网格可取:
（5）
则原子集形式为
（6）
其中原子为：
（7）
原子集构造完毕后，要将压缩感知理论应用到探地雷达成像处理上，还应该考虑探地雷达具体的工作方式以及探测环境。在此，采用如图一所示的系统仿真模型，对探测区域进行沿X轴的B扫描，形成探地雷达二维图像。
4、算法仿真
为保证计算精度，要设定正演模拟的一系列参数。设土壤的极化、磁化性质都是均匀、线性各向同性的，其相对介电常数，电导率；采用中心频率为900M的振子天线，模型区域为0.6m×0.6m矩形区域，目标体为理想导体，如图一所示。
图一 系统仿真模型
在0.6m×0.6m的地下区域内，选择扫描道数为41道，、方向上的步长均为2.5mm，时间窗为S。设有两个1米长横截面为1mm×1mm大小的钢铁长方体，探地雷达沿着与长方体垂直的方向行进，这样可得到探地雷达二维图像。将0.6m×0.6m的探测区域划分为1357×41个网格，探地雷达在每一道扫描位置不移动位置时，可以得到每一道雷达回波幅值，如图三所示为第一道雷达回波幅值图，将41道雷达回波幅值都显示在一张图上得到如图四所示的雷达回波幅值图。鉴于二维图像的准确直观方便性，以及探地雷达一般为功率图像，二维图像每个像素点值的大小与地下的后向散射系数存在一一对应的关系，可以将扫描道数作为横轴，将回波延迟时间作为纵轴，将后向散射系数作为像素点灰度值的大小，从而得到如图五所示的探地雷达图像。
采用高斯随机矩阵作为投影矩阵对1357×41维的矩阵进行线性测量，采用公式4所示的原子集和高斯随机矩阵作为恢复基，在变换域中对线性测量值采用ROMP重构算法[13,14]进行恢复，压缩采样过程如图六所示，在图六中，采用随机矩阵对扫描的41道探地雷达信号随机测量30道，对1357个回波延迟用高斯随机矩阵测量256个值，得到如图七所示的256×30线性测量矩阵，将其作为恢复矩阵。为了选取更适合本仿真模型的重构算法参数，在matlab中，对算法的迭代次数和残差阈值设定了不同的值进行试验，所得仿真参数见表一：
迭代次数 残差阈值 PSNR（运行5次）
第1次 第2次 第3次 第4次 第5次
50 3 8.1675 7.8406 8.1943 8.3560 7.9053
9 7.9388 8.2143 7.9686 8.1770 8.0845
15 7.9591 8.4226 8.1461 8.0270 7.8954
30 7.9025 8.2905 7.8842 8.1506 8.0129
100 3 13.641 14.470 14.197 14.385 15.515
9 15.514 13.545 15.615 15.895 14.565
15 14.471 15.470 14.250 15.477 15.670
30 13.545 15.535 15.745 15.565 15.845
200 3 28.771 29.731 28.721 28.731 28.779
9 30.791 29.772 29.872 30.752 30.592
15 29.578 30.778 30.448 29.537 30.573
30 30.570 30.568 29.578 30.528 30.078
表一 ROMP算法参数表
由表一可见，残差阈值的选取较迭代次数而言对重构图像的峰值信噪比的影响较小，迭代次数的增加可以较为明显的提高峰值信噪比，但是当迭代次数增加时，程序运行时间会成指数规律上升，综合考虑重构图像的峰值信噪比和程序的运行时间，针对如图一所示的系统仿真模型，选取残差阈值为9，迭代次数为100作为算法参数。选取了9为残差阈值，100为迭代次数的算法参数，在基本不影响重构效果的前提下，可有效缩短程序运行时间。虽然峰值信噪比有所减少，但针对本次2个理想目标点的结构而言重构效果已经较好，重构效果见图八。
图二 原始目标点
图三 回波幅值
图四 41道回波幅值
图五 雷达回波
图六 随机测量
图七 恢复矩阵
图八 目标重构图
5、结束语
实验结果表明，基于原子集的探地雷达成像算法具有数据采集时间短、成像速度快、算法快速稳定的特点。但是，本文的模拟目标点少，测量矩阵只是随机矩阵，重构算法采用了凸优化方法，这些是有待提高的地方。下一步，将对多目标点乃至线目标进行仿真研究，在投影矩阵方面探索诸如Chirp测量矩阵、Alltop序列测量矩阵、随机卷积测量矩阵，在重构算法方面探索更适合于二维图像的最小全变分法等。
关键词：探地雷达；压缩感知；成像算法；原子集
中图分类号：C35文献标识码： A
1、引言
探地雷达(GPR)通过发射电磁波，利用电磁波在媒质电磁特性不连续处产生的反射和散射实现浅表层成像进而辨识地表电磁特性的变化，实现对地下目标的无损探测。[1]近年来，探地雷达技术的快速发展得益于硬件电子工艺技术和数字信号处理技术的发展，而理论上没有质的飞跃。探地雷达的有效应用，故而离不开硬件系统和信号处理技术，但是如何引入新的理论、新的探地雷达成像算法和特征检测算法，是探地雷达技术发展的一个方向。现有的探地雷达遵循Nyquist-Shannon采样理论，常用成像算法有距离多普勒成像算法(RD)[2]、时域反向投影成像算法(Back Projection, BP)[3]、基于波场延拓的距离偏移成像算法[4]、逆时偏移成像算法[5]和Chirp Scaling成像算法[6]等。探地雷达具有成像速度快、空间分辨率高的特点，这就要求探地雷达具备复杂的匹配滤波器、高密度的采样孔径以及长时间的单孔径测量，最终得到庞大的数据量。
近年来，由Donoho, Candes, Romberg 和Tao等人[7-10]提出的压缩感知(CS)理论是一种充分利用信号稀疏性或可压缩性的全新信号获取和重构理论。压缩感知理论指出，只要信号是可压缩的或在某个变换域是稀疏的，就可以用一个与变换基不相关的观测基将变换所得的高维信号投影到一个低维空间上，最后通过求解一个优化问题就可以从低维空间的少量投影中以高概率重构原始信号。压缩感知理论是对传统信息理论的发展，突破了Nyquist-Shannon采样定理的极限，是科研领域新的研究热点。
本文将压缩感知理论应用到探地雷达成像上，根据成像区域目标稀疏性的先验性知识，构造合适的原子集对探测区域进行压缩测量从而降低采样孔径密度，缩短单孔径测量时间、减少获取数据量。仿真和实验结果表明，采用压缩感知探地雷达成像算法，可获得较传统探地雷达成像算法更好的成像效果。
2、压缩感知理论概述
若有一个k稀疏的长为N的实离散时间信号x(n)，n=1,2，…n。设x为一个N维列向量，即x∈RN×1，如果信号x是稀疏的或可压缩的，那么存在一组稀疏基可将x表示为
，
α和x均为N×1的矩阵，ψ为N×N的矩阵。
将此稀疏信号投影到一组测量基上，得到M(M
来重构原始信号，然而零范數约束优化重构算法是NP（Non-deterministic Polynomial）问题，可通过求解l1范数约束次最优化来解决，即通过
（2）
来从M个线性测量中恢复出α。
3、压缩感知探地雷达成像算法
3.1理想点散射模型
从数据分析的角度看，探地雷达成像问题即信号的表示问题，由探地雷达成像理想点散射模型描述可知，对于含有P个理想点散射中心的目标而言，目标后向散射场表示为[12]
（3）
其中，f为频率，θ为方位角，c为光速，Ak，xk，yk分别表示第k个散射中心的复振幅和其在散射场中所处的位置坐标。如果将上述表达式（3）中的下面这一项作为一组基函数，估计目标散射幅度参数的问题，即可看成是有限的观测信号在特定基下的表示问题。
3.2 GPR原子集的构造
对于GPR图像而言，可基于目标散射模型来构造与之相适应的原子集。将成像区域划分为M×N大小的网格，目标后向散射场表示为:
（4）
在上式中，τ1=0,1，…M-1；τ2=0,1，…N-1。
设在M×N网格上每个点都是一个散射中心，若有目标，则其为真正的散射中心，系数A(k1，k2)不为零，否则，系数A(k1，k2)为零。
因此，此时的参数网格可取:
（5）
则原子集形式为
（6）
其中原子为：
（7）
原子集构造完毕后，要将压缩感知理论应用到探地雷达成像处理上，还应该考虑探地雷达具体的工作方式以及探测环境。在此，采用如图一所示的系统仿真模型，对探测区域进行沿X轴的B扫描，形成探地雷达二维图像。
4、算法仿真
为保证计算精度，要设定正演模拟的一系列参数。设土壤的极化、磁化性质都是均匀、线性各向同性的，其相对介电常数，电导率；采用中心频率为900M的振子天线，模型区域为0.6m×0.6m矩形区域，目标体为理想导体，如图一所示。
图一 系统仿真模型
在0.6m×0.6m的地下区域内，选择扫描道数为41道，、方向上的步长均为2.5mm，时间窗为S。设有两个1米长横截面为1mm×1mm大小的钢铁长方体，探地雷达沿着与长方体垂直的方向行进，这样可得到探地雷达二维图像。将0.6m×0.6m的探测区域划分为1357×41个网格，探地雷达在每一道扫描位置不移动位置时，可以得到每一道雷达回波幅值，如图三所示为第一道雷达回波幅值图，将41道雷达回波幅值都显示在一张图上得到如图四所示的雷达回波幅值图。鉴于二维图像的准确直观方便性，以及探地雷达一般为功率图像，二维图像每个像素点值的大小与地下的后向散射系数存在一一对应的关系，可以将扫描道数作为横轴，将回波延迟时间作为纵轴，将后向散射系数作为像素点灰度值的大小，从而得到如图五所示的探地雷达图像。
采用高斯随机矩阵作为投影矩阵对1357×41维的矩阵进行线性测量，采用公式4所示的原子集和高斯随机矩阵作为恢复基，在变换域中对线性测量值采用ROMP重构算法[13,14]进行恢复，压缩采样过程如图六所示，在图六中，采用随机矩阵对扫描的41道探地雷达信号随机测量30道，对1357个回波延迟用高斯随机矩阵测量256个值，得到如图七所示的256×30线性测量矩阵，将其作为恢复矩阵。为了选取更适合本仿真模型的重构算法参数，在matlab中，对算法的迭代次数和残差阈值设定了不同的值进行试验，所得仿真参数见表一：
迭代次数 残差阈值 PSNR（运行5次）
第1次 第2次 第3次 第4次 第5次
50 3 8.1675 7.8406 8.1943 8.3560 7.9053
9 7.9388 8.2143 7.9686 8.1770 8.0845
15 7.9591 8.4226 8.1461 8.0270 7.8954
30 7.9025 8.2905 7.8842 8.1506 8.0129
100 3 13.641 14.470 14.197 14.385 15.515
9 15.514 13.545 15.615 15.895 14.565
15 14.471 15.470 14.250 15.477 15.670
30 13.545 15.535 15.745 15.565 15.845
200 3 28.771 29.731 28.721 28.731 28.779
9 30.791 29.772 29.872 30.752 30.592
15 29.578 30.778 30.448 29.537 30.573
30 30.570 30.568 29.578 30.528 30.078
表一 ROMP算法参数表
由表一可见，残差阈值的选取较迭代次数而言对重构图像的峰值信噪比的影响较小，迭代次数的增加可以较为明显的提高峰值信噪比，但是当迭代次数增加时，程序运行时间会成指数规律上升，综合考虑重构图像的峰值信噪比和程序的运行时间，针对如图一所示的系统仿真模型，选取残差阈值为9，迭代次数为100作为算法参数。选取了9为残差阈值，100为迭代次数的算法参数，在基本不影响重构效果的前提下，可有效缩短程序运行时间。虽然峰值信噪比有所减少，但针对本次2个理想目标点的结构而言重构效果已经较好，重构效果见图八。
图二 原始目标点
图三 回波幅值
图四 41道回波幅值
图五 雷达回波
图六 随机测量
图七 恢复矩阵
图八 目标重构图
5、结束语
实验结果表明，基于原子集的探地雷达成像算法具有数据采集时间短、成像速度快、算法快速稳定的特点。但是，本文的模拟目标点少，测量矩阵只是随机矩阵，重构算法采用了凸优化方法，这些是有待提高的地方。下一步，将对多目标点乃至线目标进行仿真研究，在投影矩阵方面探索诸如Chirp测量矩阵、Alltop序列测量矩阵、随机卷积测量矩阵，在重构算法方面探索更适合于二维图像的最小全变分法等。