CN106646466A - 一种基于主成分分析的加权后向投影算法的成像方法 - Google Patents

Abstract

本发明提供一种基于主成分分析的加权后向投影算法的成像方法，其方法包括：探地雷达信号预处理得到噪声和直达波得到抑制的数据；将预处理后的结果通过时延估计得到理想的时延坐标；将时延坐标位置取窗得到时延位置的幅值区域；利用主成分分析法求各道时延位置区域的能量和作为加权系数；最后利用加权系数和后向投影算法叠加得到成像结果。利用本发明可有效抑制噪声和干扰，提高信噪比和目标检测概率。

Description

一种基于主成分分析的加权后向投影算法的成像方法

技术领域

本发明涉及的是探地雷达信号处理技术、合成孔径成像技术、灾害地质调查、考古调查、公路工程质量检测以及地质勘探等领域，特别是涉及一种基于主成分分析的加权后向投影算法的成像方法。

背景技术

探地雷达(Ground Penetrating Radar,GPR)是用来对地下目标或场景进行探测的雷达系统。利用电磁波在地下介质电磁特性不连续之处产生的反射和散射现象，GPR可以获取地下目标或场景的几何和物理信息。与传统的对地探测方式相比，GPR具有不破坏探测场景、可重复性强、采样密集、安全可靠、高分辨率、操作灵活、经济性强等优势。直接使用一维或二维回波来完成检测和识别工作不方便直观，因此GPR成像技术应运而生，但是受到地下干扰严重，成像效果不佳。

在探地雷达成像方面，文献《基于Stolt偏移的探地雷达合成孔径成像研究》主要是通过一种新的stolt偏移插值实现方法，用于探地雷达合成孔径成像，克服传统方法中偏移能量不集中的缺点并保持处理速度快的优点，文献《实现浅地层探地雷达快速合成孔径成像的一种有效方法》主要是通过基于A-scan能量的特征来确定目标回波数据在整个回波数据中的分布区域，并只用该区域数据参与探地雷达合成孔径成像运算，从而排除非目标回波数据对合成孔径成像的影响，减少合成孔径运算量，实现快速合成孔径成像，文献《基于因式分解BP的车载前视地表穿透SAR快速成像算法》主要是针对双站SAR模型分析FFBP的距离误差，给出了极角采样率限制条件，提出了双站FFBP算法，通过合理设置子孔径中心提高聚焦效果，并结合距离相位补偿和线性差值技术，进一步提高聚焦精度。文献《基于幅度加权的探地雷达后向投影算法》主要是利用边缘提取获得双曲线回波上的点对应的响应幅值构造窗函数，然后对各通道的散射响应幅值进行加权成像，文献《一种基于Hilbert变换的探地雷达成像方法》主要是通过EEMD方法分解出本征模量，然后进行Hilbert-Huang变换，求出顺时振幅进行叠加成像，文献《一种用于探地雷达成像的散射强度加权处理方法》利用成像区域中对应的时延曲线处的散射数据求加权系数，而本发明利用主成分分析法求加权系数去提高成像效果，起到较好的旁瓣和杂波抑制效果，便于对地下目标进行检测和识别。

发明内容

本发明的目的在于克服现有技术的缺点与不足，提供了一种基于主成分分析的加权后向投影算法的成像方法，起到较好的旁瓣和杂波抑制效果。

本发明的目的通过以下技术方案实现：一种基于主成分分析的加权后向投影算法的成像方法，包括如下步骤：

步骤1、探地雷达采集数据，对采集到的数据信号进行预处理；

步骤2、将预处理后的结果通过时延估计得到时延坐标；将时延坐标位置取窗得到时延位置的幅值区域；

步骤3、利用主成分分析法求各道时延位置幅值区域的能量和作为加权系数；最后利用加权系数和后向投影算法叠加得到成像结果。

进一步的，所述预处理包括噪声抑制、天线耦合波、直达波抑制和射频干扰抑制。

进一步的，所述步骤2具体为：

1)直线x＝0将场景分为两部分，上半部分为空气层，波速为光速c，下半部分为各向同性的均匀介质，波速其中ε_r代表媒质的相对介电常数；

2)收发天线的位置(x_k,-h)，A(x_A,z_A)为目标所在位置，(x_r,0)表示折射点位置，根据Snell折射定律，有如下关系式成立：

其中θ_i和θ_r分别表示入射角和折射角，利用几何关系，可知：

故可得到下面一元四次方程：

通过二分法求解此一元四次方程组，求得折射点的x_r值，则各道回波上对应于A点的时延分别为：

求解一组时间单元方向上的时延坐标，利用平移不变性原理便可完成整个场景中的时延计算；

3)利用发射信号的中心频率f₀和等效采样频率Fs，求得散射响应范围利用此范围结合时延位置，得到时延点对应的散射响应矩阵P：

K表示合成孔径的长度,其中表示第k个等效阵元位置到A点对应在时延位置的幅值，表示在第k道数据中距离第i个位置的幅值，若下坐标超过原数据的范围，则将其设为0。

进一步的，所述步骤3具体为：利用主成分分析法可得P矩阵中相关性强的能量为其最近秩1矩阵，计算最近秩1矩阵，并将该矩阵的能量和作为加权系数，具体为：

对P矩阵进行奇异值分解：

其中U＝{U_m,i}＝{u₁,u₂,…,u_S+1}∈R^(S+1)×(S+1),V＝{V_n,i}＝{v₁,v₂,…,v_K}∈R^K×K，U、V均为正交矩阵，U矩阵是属于大小为(S+1)*(S+1)的实数矩阵，V矩阵同理，T为转置符号，R为实数，D∈R^S+1×K为对角阵，对角元素即为奇异值，按由大到小排列，可以表示为D_1,1≥D_2,2≥…≥D_L,L；

令σ_i＝D_i,i，由奇异值分解的性质可知为矩阵P的“最近”秩l矩阵；故其最近秩1矩阵可表示为：

P⁽¹⁾＝σ₁W₁

利用最近秩1矩阵内所有值的能量和作为加权系数，即

利用加权系数和后向投影算法将点A在各道上的散射响应幅值进行叠加，完成对A点的成像E_A：

附图说明

图1本发明的方法流程图；

图2探地雷达原数据B-scan图；

图3探地雷达信号预处理数据B-scan图；

图4 GPR成像BP算法场景；

图5成像点位置关系的平移不变性；

图6时延曲线的平移不变性

图7探地雷达后向投影算法成像图；

图8主成分分析成像图；

图9基于主成分分析的加权后向投影算法成像图；

具体实施方式

下面将结合本发明实施例中的附图对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

结合图1，一种基于主成分分析的加权后向投影算法的成像方法，包括如下步骤：

步骤1、探地雷达采集数据，对采集到的数据信号进行预处理；图2给出了预处理前的探地雷达原数据，图3给出了预处理后的探地雷达数据；

步骤2、将预处理后的结果通过时延估计得到时延坐标；将时延坐标位置取窗得到时延位置的幅值区域；

步骤3、利用主成分分析法求各道时延位置幅值区域的能量和作为加权系数；最后利用加权系数和后向投影算法叠加得到成像结果。

其中预处理后通过计算各道数据的时延和计算中心频率对应时间单元数目，保存时延位置对应的幅值和幅值区域值，利用主成分分析法计算互相关系数作为后向投影算法加权系数，从而得到成像结果。

所述预处理包括噪声抑制、天线耦合波、直达波抑制和射频干扰抑制。

结合图4，所述步骤2具体为：

1)直线x＝0将场景分为两部分，上半部分为空气层，波速为光速c，下半部分为各向同性的均匀介质，波速其中ε_r代表媒质的相对介电常数；

2)收发天线的位置(x_k,-h)，A(x_A,z_A)为目标所在位置，(x_r,0)表示折射点位置，根据Snell折射定律，有如下关系式成立：

其中θ_i和θ_r分别表示入射角和折射角，利用几何关系，可知：

故可得到下面一元四次方程：

通过二分法求解此一元四次方程组，求得折射点的x_r值，则各道回波上对应于A点的时延分别为：

结合图5、图6，求解一组时间单元方向上的时延坐标，利用平移不变性原理便可完成整个场景中的时延计算；

3)利用发射信号的中心频率f₀和等效采样频率Fs，求得散射响应范围利用此范围结合时延位置，得到时延点对应的散射响应矩阵P：

K表示合成孔径的长度,其中表示第k个等效阵元位置到A点对应在时延位置的幅值，表示在第k道数据中距离第i个位置的幅值，若下坐标超过原数据的范围，则将其设为0。

结合图7、图8、图9，所述步骤3具体为：利用主成分分析法可得P矩阵中相关性强的能量为其最近秩1矩阵，计算最近秩1矩阵，并将该矩阵的能量和作为加权系数，具体为：

对P矩阵进行奇异值分解：

其中U＝{U_m,i}＝{u₁,u₂,…,u_S+1}∈R^(S+1)×(S+1),V＝{V_n,i}＝{v₁,v₂,…,v_K}∈R^K×K，U、V均为正交矩阵，U矩阵是属于大小为(S+1)*(S+1)的实数矩阵，V矩阵同理，^T为转置符号，R为实数，D∈R^S+1×K为对角阵，对角元素即为奇异值，按由大到小排列，可以表示为D_1,1≥D_2,2≥…≥D_L,L；

令σ_i＝D_i,i，由奇异值分解的性质可知为矩阵P的“最近”秩l矩阵；故其最近秩1矩阵可表示为：

P⁽¹⁾＝σ₁W₁

利用最近秩1矩阵内所有值的能量和作为加权系数，即

利用加权系数和后向投影算法将点A在各道上的散射响应幅值进行叠加，完成对A点的成像E_A：

利用主成分分析方法加权后向投影算法的成像结果如图9所示，图7为探地雷达后向投影算法成像图，图8为主成分分析方法成像图，比较图7和图9可以发现，通过加权的后向投影算法可以有效抑制干扰和噪声，提高信噪比和目标检测概率。

以上对本发明所提供的一种一种基于主成分分析的加权后向投影算法的成像方法，进行了详细介绍，本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本发明的限制。

Claims (4)

1.一种基于主成分分析的加权后向投影算法的成像方法，其特征在于，包括如下步骤：

步骤1、探地雷达采集数据，对采集到的数据信号进行预处理；

步骤2、将预处理后的结果通过时延估计得到时延坐标；将时延坐标位置取窗得到时延位置的幅值区域；

步骤3、利用主成分分析法求各道时延位置幅值区域的能量和作为加权系数；最后利用加权系数和后向投影算法叠加得到成像结果。

2.根据权利要求1所述的方法，其特征在于，所述预处理包括噪声抑制、天线耦合波、直达波抑制和射频干扰抑制。

3.根据权利要求1所述的方法，其特征在于，所述步骤2具体为：

1)直线x＝0将场景分为两部分，上半部分为空气层，波速为光速c，下半部分为各向同性的均匀介质，波速其中ε_r代表媒质的相对介电常数；

2)收发天线的位置(x_k,-h)，A(x_A,z_A)为目标所在位置，(x_r,0)表示折射点位置，根据Snell折射定律，有如下关系式成立：

$\frac{{\mathrm{sin\θ}}_i}{{\mathrm{sin\θ}}_r}=\sqrt{{\mathrm{\ϵ}}_r}$

其中θ_i和θ_r分别表示入射角和折射角，利用几何关系，可知：

${\mathrm{sin\θ}}_i=\frac{x_r-x_k}{\sqrt{{(x_r-x_k)}^2+h^2}}$

${\mathrm{sin\θ}}_r=\frac{x_A-x_r}{\sqrt{{(x_A-x_r)}^2+{z_A}^2}}$

故可得到下面一元四次方程：

$\frac{{(x_r-x_k)}^2}{{(x_r-x_k)}^2+h^2}\·\frac{{(x_A-x_r)}^2+{z_0}^2}{{(x_A-x_r)}^2}={\mathrm{\ϵ}}_r$

通过二分法求解此一元四次方程组，求得折射点的x_r值，则各道回波上对应于A点的时延分别为：

${\mathrm{\τ}}_{A,k}=\frac{2\sqrt{{(x_r-x_k)}^2+h^2}}{c}+\frac{2\sqrt{{(x_A-x_r)}^2+{z_A}^2}}{v}$

求解一组时间单元方向上的时延坐标，利用平移不变性原理便可完成整个场景中的时延计算；

3)利用发射信号的中心频率f₀和等效采样频率Fs，求得散射响应范围利用此范围结合时延位置，得到时延点对应的散射响应矩阵P：

$P={\left[\begin{array}{ccccc}{x}_{(x_A,z_A)(-S/2,1)}& ...& x_{(x_A,z_A)(-S/2,k)}& ...& x_{(x_A,z_A)(-S/2,K)}\\ \·& & \·& & \·\\ \·& & \·& & \·\\ \·& & \·& & \·\\ x_{(x_A,z_A)(0,1)}& ...& x_{(x_A,z_A)(0,k)}& ...& x_{(x_A,z_A)(0,K)}\\ \·& & \·& & \·\\ \·& & \·& & \·\\ \·& & \·& & \·\\ x_{(x_A,z_A)(+S/2,1)}& ...& x_{(x_A,z_A)(+S/2,k)}& ...& x_{(x_A,z_A)(+S/2,K)}\end{array}\right]}_{(S+1)\×K}$

K表示合成孔径的长度,其中表示第k个等效阵元位置到A点对应在时延位置的幅值，表示在第k道数据中距离第i个位置的幅值，若下坐标超过原数据的范围，则将其设为0。

4.根据权利要求1所述的方法，其特征在于，所述步骤3具体为：利用主成分分析法可得P矩阵中相关性强的能量为其最近秩1矩阵，计算最近秩1矩阵，并将该矩阵的能量和作为加权系数，具体为：

对P矩阵进行奇异值分解：

$P={\mathrm{UDV}}^T=\underset{i=1}{\overset{L}{\Σ}}{u}_i{D}_{i,i}{v}_i^T,L=min\{S+1,K\}$

其中U＝{U_m,i}＝{u₁,u₂,…,u_S+1}∈R^(S+1)×(S+1),V＝{V_n,i}＝{v₁,v₂,…,v_K}∈R^K×K，U、V均为正交矩阵，U矩阵是属于大小为(S+1)*(S+1)的实数矩阵，V矩阵同理，^T为转置符号，R为实数，D∈R^S+1×K为对角阵，对角元素即为奇异值，按由大到小排列，可以表示为D_1,1≥D_2,2≥…≥D_L,L；

令σ_i＝D_i,i，由奇异值分解的性质可知为矩阵P的“最近”秩l矩阵；故其最近秩1矩阵可表示为：

P⁽¹⁾＝σ₁W₁

利用最近秩1矩阵内所有值的能量和作为加权系数，即

${\mathrm{\ρ}}_A=\underset{i=1}{\overset{S+1}{\Σ}}\underset{j=1}{\overset{K}{\Σ}}{P}^{\left(1\right)}$

利用加权系数和后向投影算法将点A在各道上的散射响应幅值进行叠加，完成对A点的成像E_A：

$E_A={\mathrm{\ρ}}_A\underset{k}{\Σ}{x}_{A,k}.$