一种粮堆中异物的识别方法

摘要

本发明涉及一种粮堆中异物的识别方法，用以解决现有技术干扰大、分类不精确、无法进行多种目标分类的问题。本方法首先利用图像处理技术和霍夫变换法在探地雷达图中对异物目标进行定位，然后在异物目标所在道上基于折射和反射定律对介电常数进行多层单道反演，最后根据反演的介电常数数值，依据介电常数与常见异物的对应关系确定异物的种类，实现异物的识别。该方法只需回波幅度信息，无需回波相位信息，因而算法更易于实现。同时，由于不同种类物体介电常数的差异，该方法对异物的分类和识别更加精确，并且可对同一区域中不同种类异物同时进行分类。

权利要求

1.一种粮堆中异物的识别方法，其特征在于，步骤如下：

a)用探地雷达探测，对粮堆中异物进行定位；

b)根据所测粮堆深度与异物深度，将粮堆沿垂直方向从上到下依次分成至少两个虚拟层；基于反演算法和反射折射定律，得出逐层反演公式：第一层介电常数：

第二层及之后层介电常数：

公式（12）中，εi是各虚拟层的介电常数，i表示从上到下层数的序号，ε0是空气介电常数，Ai是各虚拟层的反射回波幅度，ki是各虚拟层的线性误差校准因子，fi是各虚拟层幅度衰减因子；公式（11）中，Am是发射电磁波的幅度，A0是第一虚拟层反射回波幅度，k1是第一虚拟层的线性误差校准因子，f1是第一虚拟层幅度衰减因子；

c)将探地雷达天线放置在粮堆表面，测量探地雷达天线对应位置垂直扫描线上各虚拟层的反射回波幅度Ai；

d)采用测井法进行两次测量确定公式（11）、（12）中的参数f1，k1，fi，ki；采用一个发射天线与两个接收天线，第一次测量时，将发射天线和接收天线分别插入第一虚拟层，根据测量的电磁波幅度衰减和相移，并结合经几何衰减校正的电磁波传导模型，确定第一虚拟层介电常数ε1，并根据探地雷达测得的发射电磁波的幅度Am，空气介电常数ε0，第一层反射回波幅度A0，结合公式（11），设f1=1，确定k1；第二次测量时，将发射天线和接收天线分别插入第二虚拟层，根据测量的电磁波幅度衰减和相移，并结合经几何衰减校正的电磁波传导模型，确定第二虚拟层介电常数ε2，考虑到是同一介质，各层衰减相同，可把f1作为各层幅度衰减因子fk，即f1=fk，k=2,…,n；并结合公式（12），设ki=k1，确定幅度衰减校准因子fi；

e)利用公式（12）迭代计算探地雷达天线对应位置垂直扫描线上各虚拟层的介电常数；

f)根据各虚拟层的介电常数，判断异物的性质和种类。

2.根据权利要求1所述的一种粮堆中异物的识别方法，其特征在于，发射电磁波的幅度Am是通过将探地雷达天线紧贴一个用做反射镜面的金属板，进行测量得到的。

3.根据权利要求1所述的一种粮堆中异物的识别方法，其特征在于，步骤c）所述探地雷达天线的发射天线连接扫频仪，接收天线连接示波器。

4.根据权利要求1所述的一种粮堆中异物的识别方法，其特征在于，步骤a）中对异物定位包括：录取原始探地雷达剖面图；对原始探地雷达剖面图进行预处理；采用加窗统计法提取异物所在区域；采用二维滤波平滑异物所在区域图像进行分割；利用霍夫变换进行目标检测与定位。

5.根据权利要求1-4中任一项所述的一种粮堆中异物的识别方法，其特征在于，步骤b）中分层是将粮堆每5cm分为一层，共128层。

说明书

一种粮堆中异物的识别方法

技术领域

[0001] 本发明涉及一种属于电磁波探测领域，特别是一种利用电磁波探测技术对大型粮堆内部多个异常目标（异物）进行识别与分类的技术。

背景技术

[0002] 目前在电磁波探测领域，对媒质内部目标进行识别的主要处理手段包括成像识别和特征变量识别。成像处理通过对多个剖面的回波信号的处理，获取目标物体的几何特征，从而根据几何特征（主要是外形）对目标加以判别。即从所呈的二维和三维图像上判别有没有目标体，如果存在目标体，从图像上确定目标体的位置及形状。在比较理想的情况下，采用有效的探地雷达成像算法，可以获得较好的成像效果，从所成图像上能够清楚的知道是否存在异物，并可以完成对异物的定位以及异物形状的识别。但成像识别算法有赖于精确的系统建模，包括激励源信号的形式、天线表面电流的计算、并矢格林函数的计算，在成像过程中，需要进行多次积分方程的求解和散射场计算，计算量巨大，不适合工程应用。

[0003] 对异物材质的特征变量识别主要是根据探地探地雷达的回波信号进行特征变量的提取而完成的。探地探地雷达探测目标是依靠发射电磁脉冲，接收目标体的回波信号来完成的，是基于电磁波的传播理论的，所以通过分析目标电磁波的波相可以进行目标的特征提取，进而实现目标识别。目前，基于特征量的异物识别法主要是回波振幅和方向分析法和目标散射回波频谱特征分析法两种方法。

[0004] 回波振幅和方向分析法根据电磁波传播理论中的反射、折射定律，电磁波反射的能量和相位取决于反射面上下层介质相对介电常数的差异，电介质界面处两侧的电磁学性质差异越大，反射波越强。当界面上层介质的介电常数大于界面下层介质的介电常数时，反射波与入射波的相位同相；反之，当界面上层介质的介电常数小于界面下层介质的介电常数时，反射波与入射波的相位反相。这是判定界面两侧介质性质与属性的依据，因而反射波的振幅和方向特征是判别目标介质的重要依据。

[0005] 目标散射回波的频谱特征分析法对介质内部目标进行超宽带窄脉冲扫描，不同目标的回波信号频谱相位具有不同的变化情况，因此，目前的一个比较可行的方法是通过比较接收信号与基准信号的相位谱，对目标材质属性进行分类。

[0006] 回波振幅和方向分析法简单实用，但由于电磁波在粮食介质中存在色散和衰减现象，如果对回波振幅不加以矫正，则分析结果可能与实际相差较大；回波频谱特征分析法对回波的信噪比，媒质的损耗要求较高，当媒质中存在多种异物时，会由于不同异物回波的相互干扰而使相位谱分析法失效。粮食是损耗较大、信噪比较低的色散媒质，目前可以采用的方法回波相位比较法干扰大、分类不精确、无法进行多种目标分类，难以取得令人满意的识别效果。

发明内容

[0007] 本发明的目的是提供一种粮堆中异物的识别方法，用以解决现有技术干扰大、分类不精确、无法进行多种目标分类的问题。

[0008] 为实现上述目的，本发明的方案步骤如下：

[0009] a) 用探地雷达探测，对粮堆中异物进行定位；

[0010] b) 根据所测粮堆深度与异物深度，将粮堆沿垂直方向从上到下依次分成至少两个虚拟层；基于反演算法和反射折射定律，得出逐层反演公式：

[0011] 第一层介电常数：

[0012] （11）

[0013] 第二层及之后层介电常数：

[0014] （12）

[0015] 公式（12）中， 是各虚拟层的介电常数，i表示从上到下层数的序号， 是空气介电常数，Ai是各虚拟层的反射回波幅度， ki是各虚拟层的线性误差校准因子，fi是各虚拟层幅度衰减因子；公式（11）中，Am是发射电磁波的幅度，A0是第一虚拟层反射回波幅度，k1是第一虚拟层的线性误差校准因子，f1是第一虚拟层幅度衰减因子；

[0016] c) 将探地雷达天线放置在粮堆表面，测量探地雷达天线对应垂直剖面上各虚拟层的反射回波幅度Ai；

[0017] d) 采用测井法进行两次测量确定公式（11）、（12）中的参数f1，k1，fi，ki；采用一个发射天线与两个接收天线，第一次测量时，将发射天线和接收天线分别插入粮堆第一虚拟层，根据测量的电磁波幅度衰减和相移，并结合经几何发散校正的电磁波传导模型，确定第一虚拟层介电常数 ，并根据探地雷达测得的发射电磁波的幅度Am，空气介电常数 ，第一层反射回波幅度A0，结合公式（11），设f1=1，确定k1；第二次测量时，将发射天线和接收天线分别插入粮堆第二虚拟层，根据测量的电磁波幅度衰减和相移，并结合经几何衰减校正的电磁波传导模型，确定第二虚拟层介电常数 ，并结合公式（12），设ki=k1，确定幅度衰减校准因子fi；

[0018] e) 利用公式（12）迭代计算探地雷达天线对应位置垂直扫描线上各虚拟层的介电常数；

[0019] f) 根据各虚拟层的介电常数，判断异物的性质和种类。

[0020] 该方法只需回波幅度信息，无需回波相位信息，因而算法更易于实现。同时，由于不同种类物体介电常数的差异，该方法对异物的分类和识别更加精确，并且可对同一区域中不同种类异物同时进行分类。具体优点如下：

[0021] 1，本方案采用图像处理技术和智能分析算法，只需借助扫描的探地雷达剖面图即可定位并识别异物种类，省去了频谱分析等复杂分析过程，使粮堆内部异物识别的过程变得相对简便，更具有智能化，因而提高了异物识别的效率。

[0022] 2，传统的通过相位比较法进行的异物识别，只能对识别结果进行粗粒度分类，如是否是异物，是金属还是非金属等，不能有效识别异物的种类。本方案通过对反演的介电常数分布进行数值分析，根据不同物质具有不同的介电常数，再根据探地雷达图上的异物的双曲线特征，能够精确估计异物的性质和种类。

[0023] 3， 传统的相位比较法在探测域有多种异物时，由于异物回波之间的相互影响，给散射回波相位与基准波相位的比较来带很大干扰，严重影响了相位比较的结果。本方案首先采用霍夫变换定位异物位置，然后在异物所在垂直扫描线上逐层进行介电常数反演，不同种类的异物在反演过程中互不影响，若采用并行算法还可加快反演运算过程，因此本方案可以适应有多种异物存在的情况。

[0024] 发射电磁波的幅度Am是通过将探地雷达天线紧贴一个用做反射镜面的金属板，进行测量得到的。步骤c）所述的测井法采用的天线，发射天线连接扫频仪，接收天线连接示波器。步骤a）中对异物定位包括：录取原始探地雷达剖面图；对原始探地雷达剖面图进行预处理；采用加窗统计法提取异物所在区域；采用二维滤波平滑异物所在区域图像进行分割；利用霍夫变换进行目标检测与定位。

附图说明

[0025] 图1是一幅原始探地雷达图；

[0026] 图2是经过预处理的探地雷达图；

[0027] 图3是经过霍夫变换的探地雷达图；

[0028] 图4是粮堆分层示意图；

[0029] 图5是电磁波在粮堆多个虚拟层的等效传播模型；

[0030] 图6是测井法测量示意图；

[0031] 图7是介电常数随粮堆深度分布示意图。

具体实施方式

[0032] 下面结合附图对本发明做进一步详细的说明，具体实施方式如下：

[0033] 1，对录取的探地雷达剖面图进行图像预处理。包括随机噪声抑制、直达波和射频干扰抑制、天线直耦波抑制、提高信噪比和图像增强等。原始探地雷达图如图1所示，可以看到图像中间靠近上方的双曲线表示的就是需要检测的异物，图像上方的亮条表示粮面反射波（称为直达波），底部的网状结构的亮条表示埋在地板中的钢筋网。可见，原始探地雷达图中除了所需检测的异物的探地雷达回波图像外，还存在着大量强弱不等的噪声。

[0034] 2， 采用加窗统计法提取异物所在区域。该步基于能量检测的统计方法进行初步检测，从大量数据中提取异物所在区域，用于快速地将绝大部分非目标区域排除。以减少后续的计算量。所加窗的窗函数可采用简单的矩形窗函数。如采用如下窗函数：

[0035] （20）

[0036] 其中M是时间窗内采样点数。窗宽度 。考虑到计算量，需要合理选择窗宽度。建立二元假设检验的模型如下：

[0037] (21)

[0038] 式中，N为扫描数， 为目标回波信号， 为噪声和干扰。设 服从正态分布，均值为0，方差为 。并设 各分量线性无关。

[0039] 对 加窗处理后为

[0040] (22)

[0041] 设 为 的二范数，加窗后二元假设的检验统计量为

[0042] (23)

[0043] 取阈值判断即可粗略划分异物所在区域，由统计曲线的峰值及选取的窗函数可以大致确定目标的位置范围，从而减少后续的计算量。

[0044] 3， 采用二维滤波平滑异物所在区域图像并进行图像分割。经过上述步骤提取出的异物所在二维区域中除目标双曲线外还存在大量随机起伏的干扰和噪声，会影响霍夫变换的效果，为此可采用二维滤波平滑图像，去除奇异亮点和非相关分量，提高异物所在区域的信噪比，以便提高参数的估计精度。如图图2是对原始的探地雷达图（图1）进行预处理的结果，处理结果显示已经滤除掉了大部分的噪声和干扰，只保留了需检测的异物的图像并进行了图像增强。

[0045] 4，利用霍夫变换进行目标检测与定位。为兼顾检测精度和计算量，改进采用边缘检测的传统霍夫变换方法，直接对提取出的异物所在区域进行霍夫变换聚焦双曲线，不进行二值化等图像预处理。并根据探地雷达回波特点，将霍夫变换累加器单纯根据位置信息改为位置、幅度信息综合投票。

[0046] 根据目标回波的双曲线特征，采用霍夫变换公式：

[0047] (24)

[0048] 式中 为异物区域内 位置上的灰度值函数，即对应目标回波幅度值。表示经图像分割后由特征点 构成的双曲线，双曲线顶点位置 组成霍夫
变换的参数空间。

[0049] 根据 的分布情况可以判断是否有目标出现，有目标出现时 会出现陡峭的高峰，无目标时多为平缓的起伏。对 搜索最大峰值，就可以得到双曲线的顶点估计值 。如图3是对经处理过的探地雷达图（图2）进行了霍夫变换后的结果。
双曲线中的浅色线就是霍夫变换后对双曲线聚焦的结果。如果已知异物是层状目标，则无需进行霍夫变换，只需进行边缘检测，以确定层状目标的位置和深度。

[0050] 5，根据所测粮堆异物深度确定探测参数并划分迭代层数。在完成了上述四步后，一般已能确定异物的位置（即异物的深度和距离测量原点的水平距离）。根据异物的深度，确定需要使用电磁波探测的深度。如使用的天线频率、时窗和每扫采样点数。如估计的异物体积较大，可以采用每道较少的采样点数，以加快迭代时间，但会影响测量结果精度。根据时窗和每扫采样点数即可确定划分的每迭代层高度，原则上迭代层的高度不小于异物尺寸。结构层数目的划分依赖时窗和每扫抽样点数这两个参数。由于是纵深测量，采用发射天线和接收天线一体化的200M天线，或发射天线和接收天线分离的80M天线，将其直接置于粮面上，以减少传播衰减和几何衰减。测量时采用定时触发的工作模式，录取多次垂直方向上的单道回波信号。如采用40ns时窗和每扫128采样点的采集参数，假定电磁波在粮堆中传播的平均速度为6cm/ns（注意电磁波传播速率需采用已知目标深度法事先标定），则数据采集深度为40ns×6cm/ns=240cm。则粮堆中每一虚拟层的深度可定义为240cm/128≈1.875cm。一般在标准的高大平房仓中，按照2002年的国家粮食行业标准（LS/T 1203-2002），粮面距地面高度为6m，因此可将时窗调整为600cm/6cm/s=100ns，每一虚拟层的深度为600cm/128≈5cm。在实际测量中，考虑到介电常数的分辨率一般较低，可将相临虚拟层进行合并，以减少层次划分加快叠代运算时间，但测量精度会相应下降。

[0051] 6， 确定发射电磁波的幅度Am。由于电磁波在通常的介质表面会发生反射和折射，所以天线从第一层介质接收到的反射回波能量存在损失。为了确定发射电磁波的幅度，可以利用金属板作为反射镜面，并使天线紧贴金属板。这样，发射天线发出的电磁波能量几乎全部被反射回接收天线，而且不存在几何发散造成的能量衰减，所接收到的回波幅度可近似认为发射波的幅度Am。

[0052] 7，将探地雷达天线放置在粮堆表面，测量探地雷达天线对应垂直剖面上各虚拟层的反射回波幅度Ai；并且进行逐层反演法，详细过程如下，介质分层如图4所示。由于粮堆是散粒体，电磁波在粮堆中传播时会发生散射。接收天线在纵向上接收到的某采样点上的回波，为该采样点对应的虚拟层面上的散射回波的总和，这个散射回波总和可近似看作是该层面的反射回波。该层面上未到达接收天线的散射波和折射波将一并合成该层面的折射波，经过这样的抽象化后，电磁波在粮堆中的传播模型如图5所示。图中忽略了二次以上的多次反射波， 分别表示第 个媒质交界面的反射回波幅度和到达时间差。

[0053] 天线在位置 处发射信号 ，所接收信号可以近似看作各层面回波的叠加：

[0054] （8）

[0055] 其中L为媒质划分的层数， 和 为各层面回波的幅度和时延。由于粮堆媒质各层近似均匀，且媒质中一般无其他目标，且是仓内探测，通过预处理噪声和干扰已经得到了很好的处理和抑制，所以所接收回波可以应用该信息模型。

[0056] 只考虑一次反射波，基于平面波的折射和反射原理，可以得到介电常数与层面回波的关系式为：

[0057] （9）

[0058] 其中 代表空气与粮食的交界面。 为第 个层面回波的平均衰减， 为第 个媒质层的介电常数， 是空气的介电常数，通常取1。基于上面的关系，可以得到介电常数的递推公式为：

[0059] （10）

[0060] 通过迭代可以逐层计算出各层的介电常数，但需要初始参数，且初始参数的精度将直接影响到递推精度。准确的初始参数通过校准获得。

[0061] 逐层反演方法一般假设天线发射平面波，并将粮堆假定为表面光滑、没有损耗和色散的分层媒质，但是在对粮仓这种近场目标测量时，天线所发射电磁波不能简单视为平面波，而且粮堆中的粮食是典型的色散媒质，为提高测量精度，将校准因子引入逐层反演公式：

[0062] （11）

[0063] （12）

[0064] 上式中，Am是反演回波幅度校准值，Ai是各层面的反射回波幅度， 是各层媒质的介电常数，i表示媒质的序号。考虑到计算出的是一维介电常数序列，故式中省略了位置x和深度z，二维连续介电常数 简化为 .fi和ki是对应不同层的校准因子，fi用于减少各层媒质的衰减，电磁波的几何扩散等因素的影响，ki是线性误差校准因子，同时又是一个总调节因子，可以减少时延和介电常数估计误差等带来的影响。

[0065] 8， 该方法是通过迭代计算各层的介电常数，故上面推导出的公式（11）、（12）需要有初始参数。第6步确定能够发射波幅度Am，上一步能够确定测量出第一层媒质的反射波幅度A0，而另外的初始参数需要采用测井法进行两次测量，采用一个发射天线与两个接收天线，第一次测量时，将发射天线和接收天线分别插入粮堆第一虚拟层，粮堆水分测量因为测量范围小，属近场工作方式，因此采用了先用平面波模型推导介电常数，然后再进行几何衰减校正的方法。

[0066] 采用经几何校正的平面波模型，采用测井法推导介电常数的方法原理如下：

[0067] 将电磁波在地层中的传播近似看成是平面波，即

[0068] （1）

[0069] 式中，k为波数。则对于由T-R1-R2组成的单发双收电极系，如图6所示，两个接收天线收到的电场强度的比值为：

[0070] （2）

[0071] 则幅度衰减为：

[0072] （3）

[0073] 相移为：

[0074] （4）

[0075] 由此可以推出：

[0076] （5）

[0077] 则介质的介电常数和电导率为：

[0078] （6）

[0079] 由于电磁波实际上是一种球面波，将它看作平面波只是一种近似，球面波具有几何衰减，因此需要进行几何衰减校正。定义：

[0080] （7）

[0081] 则几何发散校正因子为As=44.65+tp1(1.324+0.184tp1)。

[0082] 经几何发散校正后的幅度衰减为EATTc=EATT-As，校正后的 =（EATT-As）/8.686。

[0083] 将 代入（6）式，即可算出测点处的介电常数。

[0084] 运用以上方法进行该步测量时，注意井中天线不要插入过深（最好不要超过第一层），以免初始介电常数 计算不准确。

[0085] 在利用测井法取得了第一层的介电常数 ，发射电磁波幅度Am和反射波幅度A0后，就可利用（11）式计算介电常数误差校准因子k1。注意：这里假定了f1=1，这是因为：天线采用了屏蔽模式，且与粮面紧耦合，既无回波幅度的衰减，也无几何扩散等因素的影响，散射的回波几乎可以全部被接收天线接收。

[0086] 再次采用测井法，将发射和接收天线插入到划分的虚拟第二层中，然后测量介电常数，根据插入的深度读取探底探地雷达回波插入点处的幅度返回值。根据上步得到的线性误差校准因子k和探底探地雷达接收到的回波幅度值，代入递推式(12)，求幅度衰减校准因子f1,考虑到是同一介质，各层衰减相同，可把f1作为各层幅度衰减因子fk，即f1=fk(k=2,…,n)。这样公式（12）的初始值和参数都确定了。

[0087] 9，在确定了迭代式(17)的初始参数和迭代系数后，即可根据探底探地雷达采集的回波幅度值，反演出各层介电常数（其中包括了异物所在层的介电常数）。

[0088] 10，绘制介电常数单道扫描时的分布图，根据突变点处的数值，判断异物的性质和种类，各种常见异物的介电常数、电导率如表1所示，可以根据表1来判断异物种类。绘制的介电常数反演的结果分布图如图7所示。从图中可以看出，异物处的介电常数与媒质中其他位置的介电常数明显不同（几乎为1），据此可以根据异物介电常数的特征推断异物的种类（如介电常数为1处的异物是预先埋设的空瓶子）。

[0089]