CN102253371B - 一种用于探地雷达成像的散射强度加权处理方法 - Google Patents
一种用于探地雷达成像的散射强度加权处理方法 Download PDFInfo
- Publication number
- CN102253371B CN102253371B CN 201110097017 CN201110097017A CN102253371B CN 102253371 B CN102253371 B CN 102253371B CN 201110097017 CN201110097017 CN 201110097017 CN 201110097017 A CN201110097017 A CN 201110097017A CN 102253371 B CN102253371 B CN 102253371B
- Authority
- CN
- China
- Prior art keywords
- imaging
- grid
- time lag
- scattering
- ground penetrating
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
Images
Landscapes
- Radar Systems Or Details Thereof (AREA)
Abstract
Description
技术领域
本发明属于探地雷达成像技术领域,涉及一种用于探地雷达成像的散射强度加权处理方法。
背景技术
探地雷达是一种有效的无损探测技术。它通过空域扫描向探测区域发射电磁波并接收散射回波,可实现对未知区域内部的成像处理,获得未知区域中的隐蔽目标参数,即目标分布信息和散射强度信息,有效应用于市政工程、考古、地雷探测、反恐等多种场合。探地雷达的空域扫描有沿线的一维扫描和在表面的二维扫描。一维扫描时,发射天线和接收天线分别以一定的间隔沿线移动。在每个位置处,发射天线向探测区域发射电磁波,接收天线接收探测区域的散射回波。然后移动发射天线和接收天线到下一个位置,又可以获得一道散射回波。通过在整个测线上移动发射天线和接收天线,便可以获得多道散射回波。发射天线和接收天线可以装配在一起同时移动,也可以分别移动。接收天线还可以选择为阵列天线的形式。这些配置方式分别对应于探地雷达应用中不同的扫描方式。本专利适用于各种探地雷达天线一维沿线扫描方式下的散射强度加权处理。探地雷达成像的目的便是从这多道散射回波,即原始记录剖面中恢复出地下区域的散射强度分布信息,即成像结果。一维扫描可以获得二维成像结果,其中一维为横向扫描维,另一维为深度维。本专利针对一维扫描下探地雷达成像中的散射强度加权处理。设一维测线沿地表布置,测线方向设为x方向,测线范围为[A,B],该测线上共有L个测点,坐标分别为xi,i=1,…,L。发射天线和接收天线装配在一起同时移动。在测点xi处,发射天线向探测区域发射电磁波,接收天线接收地下探测区域的散射回波,该点处的一维散射回波记为si(t)=[si(t0)…si(tk)…si(tK-1)]T,其中K表示时间维采样点数,上标T表示转置。其采样时窗为W=tK-1-t0。则整个记录剖面数据E0(x,t)可表示为E0(x,t)=s1(t)…si(t)…sL(t)],即E0(x,t)为一个二维矩阵,其尺寸为K×L。
成像处理前,需设定成像区域并预知探测区域的背景媒质的电磁参数。对于下视探地雷达***而言,成像区域的横向维矢量的取值区间一般取为原始扫描的测线范围。对于前视或斜前视的探地雷达***而言,该取值区间需要根据具体的探测场景加以确定,此处统一记为[ha,hb]。深度维矢量需根据探地雷达的探测深度进行选取,与时窗W有关,记为[ga,gb]。设定成像区域后,需要对该成像区域进行二维离散化处理,将该区域分别沿深度维和横向维等间隔地划分为Lz和Lx个网格,则整个成像区域划分为Lz×Lx个网格。成像的目的便是获得该Lz×Lx个网格处的散射强度值,即O(zm,xn) m=1,…,Lz;n=1,…,Lx。
探地雷达的成像方法有多种,基于“延时-累加”处理的后向投影成像算法适用于非等间距采样下对复杂有耗媒质中点散射型目标的成像处理,广泛应用于探地雷达信号处理中。成像处理前需要对原始数据进行均衡、解振荡、去直达波、零点校正等预处理,设预处理后的记录剖面数据为E1(x,t)=[s′1(t)…s′i(t)…s′L(t)],仍为L列。对成像区域中的每一网格位置(zm,xn),根据探测扫描场景计算各测点xi处对应的电磁波双程传播时延τm,n,i,然后提取该测点处的一维回波信号s′i(t)在时刻τm,n,i处的值qm,n,i,从而生成一维信号wm,n=[qm,n,1,…,qm,n,L]。传统的成像方法是将该一维信号求和作为网格(zm,xn)处的成像结果,即(zm,xn)处的散射强度为遍历成像区域中所有的网格,分别计算各点的散射强度O(zm,xn)即可获得整个成像区域的成像结果O(zm,xn) m=1,…,Lz;n=1,…,Lx。
传统的“延时-累加”成像方法直接将成像网格处对应的时延曲线上的散射数据进行求和处理作为该网格处的成像结果,并没有充分利用各时延曲线上散射数据的统计特征,没有对各时延曲线上的散射数据进行不同的加权处理。因此,传统方法的成像结果分辨率差,旁瓣干扰大。
发明内容
本发明所要解决的技术问题是提出一种用于探地雷达成像的散射强度加权处理方法,该用于探地雷达成像的散射强度加权处理方法能提高探地雷达的成像质量。
本发明的技术解决方案如下:
一种用于探地雷达成像的散射强度加权处理方法,所述的用于探地雷达成像的散射强度加权处理方法包括以下步骤:
步骤1:提取散射数据:
提取成像区域中某网格(zm,xn)对应的时延曲线处的散射数据, 即wm,n=[qm,n,1,…,qm,n,L],m和n分别为成像区域中深度维网格下标和横向维网格下标,zm,xn分别表示纵向和横线的网格序号,L指接收阵列天线中阵元的数量;
步骤2:基于加权因子获取成像结果:
遍历成像区域中所有的网格,分别计算各网格(zm,xn)处的加权因子α(zm,xn)和各网格(zm,xn)处对应的时延曲线处散射数据之和进而获得整个成像区域的成像结果其中m=1,…,Lz;n=1,…,Lx;Lz和Lx分别为深度维和横向维的成像网格的数目。
Lz和Lx是根据探地雷达的有效探测深度和天线的有效覆盖范围确定的。
加权因子α(zm,xn)的计算方法为以下两种方法中的任意一种方法:
方法1: 其中m1和s分别为成像点(zm,xn)对应的时延曲线处的散射回波的均值和标准差,计算公式为:
方法2: 其中s、m1和e分别为成像点(zm,xn)对应的时延曲线处的散射回波的标准差、均值和能量值,能量值计算公式为
本发明的构思是:先将各成像单元对应的时延曲线上的散射回波数据在累加后,再乘以一个加权因子,作为各成像单元的散射强度值。加权因子是根据时延曲线上的散射回波数据计算得到的。本发明适用于各种扫描方式下的探地雷达成像处理。
当某个成像单元处有真实的点目标存在时,时延曲线处的散射回波体现出很好的一致性。而当某个成像单元处没有真实的点目标存在时,时延曲线处的散射回波体现出较强的波动性。本发明通过对时延曲线处散射回波的均值、标准差和能量值进行运算获得加权因子,定量描述了这种一致性,运用加权因子对各时延曲线处散射回波之和进行加权处理,获得了更好的成像结果。
有益效果:
本发明提出了一种用于探地雷达成像的散射强度加权处理方法。其特点在于:在后向投影成像过程中,对各成像单元对应的时延曲线上的散射数据进行累加处理后,再乘以一个加权因子作为各成像单元的散射强度值。加权因子是通过时延曲线上的散射数据计算得到的。现有的技术是直接将时延曲线上的散射数据进行累加处理获得成像结果的,没有进行加权处理。实验结果显示,本发明提出的方法与现有的成像方法相比,成像质量显著提高。
附图说明
图1示出了实孔径天线对探测区域的扫描示意图;
图2示出了图1所示的探测区域的空域散射回波;
图3示出了图2所示的原始散射回波经预处理后的雷达记录剖面;
图4示出了基于“延迟-累加”运算的雷达成像算法处理示意图;【图4.1为预处理后的雷达记录剖面,图4.2为成像区域,图4.3为实验曲线处的散射回波,图4.4为各分辨单元处的加权因子值(每个分辨单元对应一个值,在图中该值是用灰度来表示的),图4.3的散射回波经过计算后得到的加权因子值对应图4.4中方格所示的分辨单元的灰度值。
图5示出了未经加权处理时对图3所示的雷达记录剖面进行处理后的成像结果;
图6示出了经过第一种加权算法处理后对图3所示的雷达记录剖面进行处理后的成像结果;
图7示出了实孔径天线对探测区域的扫描示意图;
图8示出了图7所示的探测区域的空域散射回波;
图9示出了图8所示的原始散射回波经预处理后的雷达记录剖面;
图10示出了未经加权处理时对图9所示的雷达记录剖面进行处理后的成像结果;
图11示出了经过第一种加权算法处理后对图9所示的雷达记录剖面进行处理后的成像结果;
图12示出了经过第二种加权算法处理后对图9所示的雷达记录剖面进行处理后的成像结果。
标号说明:1-发射天线,2-接收阵列天线。
具体实施方式
以下将结合附图和具体实施例对本发明做进一步详细说明:
实施例1:
本实例是针对均匀背景媒质中目标的实孔径探测成像,但本技术并不局限于均匀背景媒质和实孔径扫描,对复杂背景媒质下合成孔径扫描和多发多收扫描方式也适用。
首先运用探地雷达对探测区域进行一维扫描,如图1所示。发射天线固定在x=1m、z=0m处,接收阵列放置在z=0m处,其横向扩展为[0.2 1.8]m,两相邻阵元间距0.04m,一共有41个阵元。探测区域如图1所示,在x=1m、z=0.7m处放置一个半径为0.02m的细长型铁棒,长度远大于其半径,导体棒的取向垂直于纸面。发射天线为线源形式,向地下发射电磁波。散射回波被阵列天线接收,共41道数据,如图2所示。图2中的原始散射数据经过预处理后的数据如图3所示。图4为基于“延迟-累加”运算的成像算法的示意图。成像处理前需要设定成像区域,本实例中的成像区域如图4.2所示。整个区域经过空域离散形成了水平向和竖直向的空间网格,每个网格便是一个分辨单元,此处的网格横向和纵向尺寸均为0.02m。成像的目的便是获得各空间网格处的散射强度值。根据探地雷达的扫描场景和电磁波传播规律,成像区域中的每个空间网格都可以计算获得一条时延曲线。对每个空间网格,根据费马原理计算电磁波从发射天线到成像网格坐标点的传播路径和时延,然后再计算电磁波从该成像网格坐标点到每个接收天线处的传播路径和时延,通过计算获得一条时延曲线。参考文献为【雷文太,刘立业,粟毅.平面分层媒质中目标的反向投影成像算法.信号处理.2007,23(5):680-685】。图4.2中的黑点对应的时延曲线在图4.1中以白线形式显示出来。该时延曲线的长度与雷达记录剖面的道数相等,也为41点。提取该时延曲线与雷达记录剖面相交处的值便可以获得一个长度为41点的一维信号,如图4.3所示。不同的空间网格对应于不同的时延曲线,也就可以获得不同的一维信号。原始成像处理是对该一维信号进行求和,然后将和值作为该空间网格处的散射强度值。遍历所有的空间网格,获得所有网格处的散射强度值,也就完成了成像运算。本专利通过对该一维信号的求和值进行加权处理以提高成像质量。本实例中运用第一种加权处理方法,即 计算获得该一维信号的加权因子后,将其显示在图4.4中对应的成像点处。遍历所有的成像点,可以获得一个二维加权因子矩阵,如图4.4所示。
原始成像方法的处理结果如图5所示。运用第一种加权处理方法后的成像结果如图6所示。从加权处理流程上讲,图6便是图5所示的原始成像结果与图4.4所示的加权因子矩阵对应点相乘的结果。从图中可见,成像质量显著提高。
实施例2:
探地雷达的天线配置同实施例1。探测区域如图7所示,在五个位置处放置细长型铁棒,长度远大于其半径,导体棒的取向垂直于纸面。发射天线为线源形式,向地下发射电磁波。散射回波被阵列天线接收,共41道数据,如图8所示。图8中的原始散射数据经过预处理后的数据如图9所示。
原始成像方法的处理结果如图10所示。运用第一种加权处理方法后的成像结果如图11所示。运用第二种加权处理方法后的成像结果如图12所示。两种加权方法都显著提高了成像质量,五个点目标的位置估计精度和分辨力都有较大程度的提高。
Claims (1)
1.一种用于探地雷达成像的散射强度加权处理方法,包括以下步骤:
步骤1:提取散射数据:
提取成像区域中某网格(zm,xn)对应的时延曲线处的散射数据,即wm,n=[qm,n,1,…,qm,n,L],m和n分别为成像区域中深度维网格下标和横向维网格下标,zm,xn分别表示纵向和横向的网格序号,L指接收阵列天线中阵元的数量;
步骤2:基于加权因子获取成像结果:
遍历成像区域中所有的网格,分别计算各网格(zm,xn)处的加权因子α(zm,xn)和各网格(zm,xn)处对应的时延曲线处散射数据之和 进而获得整个成像区域的成像结果 其中m=1,…,Lz;n=1,…,Lx;Lz和Lx分别为深度维和横向维的成像网格的数目;
加权因子α(zm,xn)的计算方法为以下两种方法中的任意一种方法:
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN 201110097017 CN102253371B (zh) | 2011-04-18 | 2011-04-18 | 一种用于探地雷达成像的散射强度加权处理方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN 201110097017 CN102253371B (zh) | 2011-04-18 | 2011-04-18 | 一种用于探地雷达成像的散射强度加权处理方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN102253371A CN102253371A (zh) | 2011-11-23 |
CN102253371B true CN102253371B (zh) | 2013-01-30 |
Family
ID=44980737
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN 201110097017 Expired - Fee Related CN102253371B (zh) | 2011-04-18 | 2011-04-18 | 一种用于探地雷达成像的散射强度加权处理方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN102253371B (zh) |
Families Citing this family (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102636501B (zh) * | 2012-03-26 | 2014-09-03 | 中南大学 | 一种剔除表层钢筋对高频电磁波影响的滤波方法 |
CN102621548B (zh) * | 2012-04-17 | 2013-11-27 | 中南大学 | 一种探地雷达多尺度后向投影成像方法 |
CN102830401B (zh) * | 2012-08-27 | 2014-09-17 | 中南大学 | 一种探地雷达加窗加权后向投影成像方法 |
CN105974405B (zh) * | 2016-05-04 | 2018-07-06 | 哈尔滨工业大学 | 基于幅度加权的探地雷达后向投影成像方法 |
CN107390213B (zh) * | 2017-07-14 | 2019-12-03 | 中南大学 | 一种基于滑动时窗的探地雷达记录剖面的时延曲线提取方法 |
Family Cites Families (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
NO322089B1 (no) * | 2003-04-09 | 2006-08-14 | Norsar V Daglig Leder | Fremgangsmate for simulering av lokale prestakk dypmigrerte seismiske bilder |
CN101285884A (zh) * | 2007-04-11 | 2008-10-15 | 汤子跃 | 一种逆合成孔径雷达成像时间选择新方法 |
CN101221239B (zh) * | 2008-01-25 | 2010-06-02 | 电子科技大学 | 一种基于水平集的合成孔径雷达图像分割方法 |
-
2011
- 2011-04-18 CN CN 201110097017 patent/CN102253371B/zh not_active Expired - Fee Related
Also Published As
Publication number | Publication date |
---|---|
CN102253371A (zh) | 2011-11-23 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN102830401B (zh) | 一种探地雷达加窗加权后向投影成像方法 | |
CN102253371B (zh) | 一种用于探地雷达成像的散射强度加权处理方法 | |
Cui et al. | Estimating tree-root biomass in different depths using ground-penetrating radar: Evidence from a controlled experiment | |
Longoni et al. | Surface and subsurface non-invasive investigations to improve the characterization of a fractured rock mass | |
Comite et al. | Multiview imaging for low-signature target detection in rough-surface clutter environment | |
CN105158808A (zh) | 一种浅海瞬变电磁海空探测及其解释方法 | |
CN102621548B (zh) | 一种探地雷达多尺度后向投影成像方法 | |
CN101915943A (zh) | 均匀背景介质的介电常数和隐蔽目标参数的联合反演方法 | |
Pérez-Gracia et al. | Horizontal resolution in a non-destructive shallow GPR survey: An experimental evaluation | |
CN103439707A (zh) | 一种探地雷达加窗距离偏移成像方法 | |
CN110458129A (zh) | 基于深度卷积神经网络的非金属地雷识别方法 | |
Ren et al. | 3D Imaging Algorithm for Down‐Looking MIMO Array SAR Based on Bayesian Compressive Sensing | |
CN103913733A (zh) | 极地冰川厚度探测方法 | |
Gloaguen et al. | Stochastic borehole radar velocity and attenuation tomographies using cokriging and cosimulation | |
Catapano et al. | Microwave tomography enhanced GPR surveys in Centaur’s Domus, Regio VI of Pompeii, Italy | |
Tajdini et al. | Real-time modeling of forward-looking synthetic aperture ground penetrating radar scattering from rough terrain | |
Lantini et al. | Application of Ground Penetrating Radar for mapping tree root system architecture and mass density of street trees. | |
Wang et al. | Velocity analysis of CMP gathers acquired by an array GPR system ‘Yakumo’: Results from field application to tsunami deposits | |
Zhou et al. | Reconstruction from antenna-transformed radar data using a time-domain reconstruction method | |
CN110297237A (zh) | 考虑天线方向图的探地雷达绕射叠加成像方法及*** | |
Wang et al. | Characteristic parameters extraction method of hidden Karst Cave from Borehole radar signal | |
Verdonck et al. | The urban survey: methodology | |
Atef et al. | Adaptive boxcar background filtering for real-time GPR utility detection | |
Li et al. | Research on 3D ground penetrating radar deep underground cavity identification algorithm in urban roads using multi-dimensional time-frequency features | |
Cao et al. | Subspace Projection Attention Network for GPR Heterogeneous Clutter Removal |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
C14 | Grant of patent or utility model | ||
GR01 | Patent grant | ||
CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20130130 |
|
CF01 | Termination of patent right due to non-payment of annual fee |