CN105117503B - 一种探地雷达数据可视化方法 - Google Patents

Abstract

本发明公开了一种探地雷达数据可视化方法，包括以下步骤：选定GPR数据集，分别提取GPR数据集不同尺度的无监督局部特征，得到GPR数据的局部特征集U；设定待匹配子序列长度范围Mlen，利用动态时间归整方法将待测的GPR数据集D与局部特征集U进行模式匹配；运用可视化方法将进行模式匹配的结果进行展示，得到GPR数据的直观图像表示。本发明有益效果：本发明对于GPR数据应用性强，并且能够更好的提取有效GPR数据特征，更为客观的呈现了GPR数据所蕴含的信息，用可视化方法直观的体现了探测结果。

Description

一种探地雷达数据可视化方法

技术领域

本发明属于计算机数据处理可视化领域，具体涉及一种基于多尺度局部特征和动态时间规整的探地雷达数据可视化方法。

背景技术

近年来，探地雷达(GroundPenetratingRadar，简称GPR)作为一种快速、连续、非接触电磁波探测技术，以其采集速度快、分辨率高的特点，成为物探领域研究的热点之一。在实际的应用中，GPR数据的解译需要较高的专业知识及经验，而且环境以及系统的干扰使得GPR数据掺杂各种杂波和噪声以及个人经验不同使得结果存在差异性，制约了从GPR技术的进一步发展和实用化。其中GPR数据解译过程中对其特征的提取传统方法本质都是基于线性预测、Fourier分析或小波分析这样的传统信号分析技术，对于地下介质类别数目就只得依赖解译人员专业知识及经验，过度或者不当的处理，都会使数据失去原有真实性。

数据可视化是描述和表达各种半结构化甚至非结构化问题的关系和模型的最佳方法和手段，可以成为解译GPR数据中蕴含信息的新的解决途径。通过可视化的方法对GPR数据进行图像显示，提供新的洞察力以挖掘传统方法提取不到的信息。

然而由于GPR数据采用声波反射测速技术所得，反射回波的不确定性和地质的多样性，使得GPR数据富含杂波、噪音和存在不连续性，以及获取的数据分辨率相对低的这些特点，给其可视化带来了极其大的困难与挑战。

发明内容

本发明的目的就是为了解决上述问题，提出了一种探地雷达数据可视化方法，该方法基于多尺度局部特征和动态时间规整，从而客观的挖掘更多有效特征信息，解决了GPR数据解译过程中过于依赖人工因素的问题。

为了解决上述问题，本发明采用如下技术方案：

一种探地雷达数据可视化方法，包括以下步骤：

(1)选定GPR数据集，分别提取GPR数据集不同尺度的无监督局部特征，得到GPR数据的局部特征集U＝{U(1),U(2),…,U(m)}；

(2)设定待匹配子序列长度范围[minlen,maxlen]，利用动态时间归整方法将待测的GPR数据集D＝{T(1),T(2),…,T(n)}与步骤(1)中所得局部特征集U进行模式匹配分类，U中的m即为此处分类数目；

(3)用可视化方法中的颜色映射方法将步骤(2)中进行模式匹配的结果进行展示，得到GPR数据的直观图像表示。

所述步骤(1)中提取GPR数据的不同尺度的无监督局部特征的方法为：

步骤S110：选取待测GPR数据集中某一序列T(ds)＝(t₁,t₂,…t_len),其中1≤ds≤n，初始化无监督局部特征数据集U＝Φ；

步骤S120：设定滑动窗口范围ql，移动步长l，从序列T(ds)中生成多个尺度的候选无监督局部特征集C；

步骤S130：计算候选集C中每一个候选序列的质量度量gap，并对其排序；

步骤S140：取质量度量gap最大的候选序列，更新到无监督局部特征数据集U中作为U(i)；

步骤S150：计算U(i)与GPR数据集中各序列的距离，设定阈值θ，移除距离小于阈值θ的序列，判断DIS_A是否稳定，若稳定，则输出特征集U；否则取距离最大的序列，返回步骤S110。

所述步骤S150中，

判断DIS_A稳定的条件具体为：

|DIS_A|＝1；

即：当小于分割点dt的距离的个数等于1达到稳定。

计算候选集C中每一个候选序列的质量度量gap的方法为：

S131：输入候选序列C(c)和数据集D＝{T(1),T(2),…,T(n)}，初始化maxGap＝0；

S132：计算候选序列C(c)与数据集D中各序列的距离，并将其排序，得到候选序列C(c)与数据集D中各序列距离的距离向量DIS＝{dis₁,dis₂,…dis_n}；

S133：根据距离向量可得n-1个分割点dt，对每一个分割点dt分别计算分割点dt两端的距离个数的比值R；

S134：判断R是否满足设定条件，若是则计算gap，否则计算下一分割点的R；

S135：判断计算的gap是否大于maxGap，是则更新maxGap，否则计算下一分割点gap。

所述步骤S133中，计算R的方法为：

$R=\frac{\left|{\mathrm{DIS}}_A\right|}{\left|{\mathrm{DIS}}_B\right|};$

其中，DIS_A表示距离向量DIS中小于分割点dt的各距离集合，DIS_B表示距离向量DIS中大于分割点dt的各距离集合，〡DIS_A〡表示小于分割点dt的距离个数,〡DIS_B〡表示大于分割点dt的距离个数。

所述步骤S134中R的设定条件为：

R∈(0.2，5)。

所述步骤S134中，计算gap的方法具体为：

gap＝μ_B-σ_B-(μ_A+σ_B)

其中，μ_B表示大于分割点dt的所有距离的平均值，σ_B表示大于分割点dt的所有距离的方差，μ_A表示小于分割点dt的所有距离的平均值，σ_A表示小于分割点dt的所有距离的方差。

所述步骤S132中，

计算候选序列C(c)与数据集D中各序列距离的距离向量的方法为：

S1321：输入一个候选序列C＝(c₁,c₂,…c_clen)和数据集D＝{T(1),T(2),…,T(n)}；

S1322：初始化距离向量DIS＝Φ；

S1323：利用滑动窗口ql＝clen，clen为当前候选序列C＝(c₁,c₂,…c_clen)的长度，依次取出T(ds)的所有子序列Z＝{Z₁,Z₂,…Z_v}，Z_i＝(z₁,z₂,…,z_clen)；

其中：

T(ds)＝(t₁,t₂,…t_len)，len为数据集T(ds)＝(t₁,t₂,…t_len)的长度；

$1\≤ds\≤n,1\≤i\≤v,v=\frac{len(len+1)}{2};$

S1324：计算候选序列C(c)与子序列集Z中各子序列之间的距离dis，取最小值为候选序列C(c)与数据序列T(ds)的距离dis_ds，更新到DIS中；

S1325：判断ds是否为n，若是，则对DIS排序得DIS＝{dis₁,dis₂,…dis_n}，否则，转向步骤S1323。

所述步骤(2)的具体方法为：

步骤S201：设定待匹配子序列长度范围[minlen,maxlen]，输入待测的GPR数据集D＝{T(1),T(2),…,T(n)}与局部特征集U＝{U(1),U(2),…,U(m)}；

步骤S202：对于T(ds)，分别初始化起始点st＝1、长度值Mlen＝minlen以及距离向量sdist＝Φ；

步骤S203：初始化i＝1；

步骤S204：计算[minlen,maxlen]范围内，与局部特征U(i)距离d_i最小的子序列M_(st，Mlen)，更新到距离向量sdist；

步骤S205：判断i是否为m，若是，则排序距离向量sdist，否则i＝i+1，转向步骤S204；

步骤S206：取子序列M_(st，Mlen)与局部特征集U中各特征的距离的距离向量中最小的i值，作为M_(st，Mlen)的类别序号；

步骤S207：判断st+Mlen-1是否为len，若是，则下一步；否则更新起始点st＝st+Mlen，并转向步骤S204；len为T(ds)的长度；

步骤S208：判断ds是否为n+1，若是，则下一步；否则返回步骤S202；

步骤S209：输出类别矩阵。

所述步骤S204中，d_i距离计算公式为：

d_i＝DTW(M_(st,Mlen),U(i))；

其中，st为待匹配序列的起始点，Mlen为[minlen,maxlen]范围内某一长度，U(i)为局部特征，M_(st，Mlen)为与局部特征U(i)距离d_i最小的子序列。

本发明的有益效果：

本发明对于GPR数据应用性强，并且能够更好的提取有效GPR数据特征，更为客观的呈现了GPR数据所蕴含的信息，用可视化方法直观的体现了探测结果。

通过对数据集的多尺度局部特征的提取，解决了未知探测区域地下介质类目问题，区分出了不同的地下介质，克服了传统GPR数据介意过程中特征提取方法所带来的多解性的缺点。

通过采用DTW距离计算能够对不等长序列进行相似的度量；并且DTW对于序列的突变或者异常点并不敏感，这对于GPR数据中异常的噪点可以起到忽略的作用。

附图说明

图1所示为本发明探地雷达数据可视化方法流程图。

具体实施方式

下文将结合具体实施例详细描述本发明。应当注意的是，下述实施例中描述的技术特征或者技术特征的组合不应当被认为是孤立的，它们可以被相互组合从而达到更好的技术效果。

如图1所示，本发明提供的一种基于多尺度局部特征和动态时间规整的探地雷达数据可视化方法包括如下步骤：

步骤S100:对GPR数据提取多尺度的局部特征；

步骤S200:利用DTW距离对S100所述特征模型与GPR数据进行模式匹配；

步骤S300:对S200所述模式匹配结果运用可视化方法得到GPR数据的直观图像表示。

上述步骤S100:对GPR数据提取多尺度的局部特征，具体地包括如下步骤：

步骤S110:输入GPR数据集中某一序列T(ds)，初始化无监督局部特征数据集U＝Φ。

步骤S120:设定滑动窗口范围ql，移动步长l，从T(ds)中生成多个尺度的候选无监督局部特征集C。

步骤S130:计算候选集C中每一个候选序列的质量度量gap，并对其排序。其中质量度量gap计算公式如下：

gap＝μ_B-σ_B-(μ_A+σ_B)

步骤S140:取最大gap的候选序列，更新到无监督局部特征数据集U中作为U(i)。

这样提出去的是某一个尺度的无监督局部特征的提取后，我们便提取到某一道GPR数据序列T(ds)中的某一位置的一个无监督局部特征，表明此道GPR数据序列中的此位置上的介质特征。然而GPR数据中可能包含更多的介质信息特征，因而我们对GPR数据中含有此类介质特征相似的数据序列去除处理后，对数据集中剩余数据序列进行下一局部特征的查找，其尺度取决于此特征的gap度量大小，因此最后我们会提取到多个尺度的局部特征。

步骤S150:计算U(i)与GPR数据集D中各序列距离，移除距离小于阈值θ的序列，判断DIS_A是否稳定，若稳定，则输出特征集U，否则取距离最大的序列，重复步骤S110。

该稳定条件为：

|DIS_A|＝1

即：当小于分割点dt的距离的个数等于1达到稳定。

上述步骤S130:计算候选集C中每一个候选序列的质量度量gap，具体地包括如下步骤：

S131:输入候选序列C(c)和数据集D＝{T(1),T(2),…,T(n)}，初始化maxGap＝0；

S132:计算候选序列C(c)与数据集D中各序列的距离，并将其排序，可得C(c)与D中各序列距离的距离向量DIS＝{dis₁,dis₂,…dis_n}；

S133:根据距离向量可得n-1个分割点dt，对每一个分割点dt计算分割点dt两端各距离个数的比值R；

R计算公式为：

$R=\frac{\left|{\mathrm{DIS}}_A\right|}{\left|{\mathrm{DIS}}_B\right|}$

S134:判断R是否满足条件，若是则计算gap，否则下一分割点；

其中R需要满足的条件是：

R∈(0.2,5)；

此条件用来控制分割点dt两边各距离数量的比例。

S135:判断是否大于maxGap，是则更新，否则计算下一分割点gap。

上述步骤S132：计算C(n)与D中各序列距离，具体地包括如下步骤：

S1321:输入一个候选序列C＝(c₁,c₂,…c_clen)和数据集D＝{T(1),T(2),…,T(n)}；

S1322:初始化距离向量DIS＝Φ；

S1323:利用滑动窗口ql＝clen，clen为当前候选序列C的长度，依次取出T(ds)＝(t₁,t₂,…t_len)的所有子序列；1≤ds≤n，len为数据集中T(ds)的长度；

S1324:计算候选序列C(c)与所有子序列距离dis，取最小值为C(c)与T(ds)的距离dis_ds，更新到DIS中；

S1325：判断ds是否为n，若是，则对DIS排序得DIS＝{dis₁,dis₂,…dis_n}，否则转向S1323。

GPR数据是通过分析GPR接收到的有效反射电磁波的特征来推断被测地下介质的空间分布状态。不同介质对于电磁波的具有不同的反射特性，表现在波形上就有不同的结果。经过多个尺度特征提取后，提取出GPR数据中蕴含的多个介质信息特征。传统的GPR数据处理方法多为去噪和抑制杂波，再加上较为依赖人工解译的主观判断，探测效果往往不尽人意。本文提出用DTW距离将提取出的局部特征与GPR数据进行模式匹配的方法，会得到更为客观的分析结果。

模式匹配是度量数据序列相似程度的方法，在数据序列分析处理中具有基础性地位。其思想在于输入未知模式与已提取的特征模式相比较，具有相同或相似匹配的模式即为该位置模式的所属类型。两序列距离越小，越相似。在这里使用DTW距离直接匹配的优势在于：一方面在于特征提取的过程中，每一个提取出的局部特征可能会有不同的尺度，此时就可以充分体现DTW对于不等长序列进行相似度量的优势；而另一方面在于DTW对于序列的突变或者异常点并不敏感，这对于GPR数据中异常的噪点可以起到忽略的作用。

上述步骤S200:上述利用DTW距离对S100所述特征模型与GPR数据进行模式匹配，具体的包括如下步骤：

步骤S201:设定待匹配子序列长度范围[minlen,maxlen]，输入待测的GPR数据集D＝{T(1),T(2),…,T(n)}与局部特征集U＝{U(1),U(2),…,U(m)}；

步骤S202:对于T(ds)，分别初始化起始点st＝1、长度值Mlen＝minlen以及距离向量sdist＝Φ；

步骤S203:初始化i＝1；

步骤S204:计算[minlen,maxlen]范围内，与局部特征U(i)距离d_i最小的子序列M_(st，Mlen)，更新到距离向量sdist；

d_i距离计算公式为：

d_i＝DTW(M_(st,Mlen),U(i))；

其中，st为待匹配序列的起始点，Mlen为[minlen,maxlen]范围内某一长度，U(i)为局部特征，M_(st，Mlen)为与局部特征U(i)距离d_i最小的子序列。

步骤S205:判断i是否为m，若是，则排序距离向量sdist，否则i＝i+1，转向步骤S204；

步骤S206:取子序列M_(st，Mlen)与局部特征集U中各特征的距离的距离向量中最小的i值，作为M_(st，Mlen)的类别序号；

步骤S207:判断st+Mlen-1是否为len，若是，则下一步；否则更新起始点st＝st+Mlen，并转向步骤S204；len为T(ds)的长度；

步骤S208:判断ds是否为n+1，若是，则下一步；否则返回步骤S202；

步骤S209:输出类别矩阵。

上述步骤S300:对S200所述模式匹配结果运用可视化方法中颜色映射技术得到GPR数据的直观图像表示。

本文虽然已经给出了本发明的一些实施例，但是本领域的技术人员应当理解，在不脱离本发明精神的情况下，可以对本文的实施例进行改变。上述实施例只是示例性的，不应以本文的实施例作为本发明权利范围的限定。

Claims (9)

1.一种探地雷达数据可视化方法，其特征是，包括以下步骤：

(1)选定GPR数据集，分别提取GPR数据集不同尺度的无监督局部特征，得到GPR数据的局部特征集U＝{U(1),U(2),…,U(m)}；

(2)设定待匹配子序列长度范围[minlen,maxlen]，利用动态时间归整方法将待测的GPR数据集D＝{T(1),T(2),…,T(n)}与步骤(1)中所得局部特征集U进行模式匹配分类，U中的m即为此处分类数目；

(3)用可视化方法中的颜色映射方法将步骤(2)中进行模式匹配的结果进行展示，得到GPR数据的直观图像表示；

所述步骤(1)中提取GPR数据的不同尺度的无监督局部特征的方法为：

步骤S110:选取待测GPR数据集中某一序列T(ds)＝(t₁,t₂,…t_len),其中1≤ds≤n，初始化无监督局部特征数据集U＝Φ；

步骤S120:设定滑动窗口范围ql，移动步长l，从序列T(ds)中生成多个尺度的候选无监督局部特征集C；

步骤S130:计算候选集C中每一个候选序列的质量度量gap，并对其排序；

步骤S140:取质量度量gap最大的候选序列，更新到无监督局部特征数据集U中作为U(i)；

步骤S150:计算U(i)与GPR数据集中各序列的距离，设定阈值θ，移除距离小于阈值θ的序列，判断DIS_A是否稳定，若稳定，则输出特征集U；否则取距离最大的序列，返回步骤S110。

2.如权利要求1所述的一种探地雷达数据可视化方法，其特征是，所述步骤S150中，判断DIS_A稳定的条件具体为：

|DIS_A|＝1；

即：当小于分割点dt的距离的个数等于1达到稳定。

3.如权利要求1所述的一种探地雷达数据可视化方法，其特征是，计算候选集C中每一个候选序列的质量度量gap的方法为：

S131:输入候选序列C(c)和数据集D＝{T(1),T(2),…,T(n)}，初始化maxGap＝0；

S132:计算候选序列C(c)与数据集D中各序列的距离，并将其排序，得到候选序列C(c)与数据集D中各序列距离的距离向量DIS＝{dis₁,dis₂,…dis_n}；

S133:根据距离向量可得n-1个分割点dt，对每一个分割点dt分别计算分割点dt两端的距离个数的比值R；

S134:判断R是否满足设定条件，若是则计算gap，否则计算下一分割点的R；

S135:判断计算的gap是否大于maxGap，是则更新maxGap，否则计算下一分割点gap。

4.如权利要求3所述的一种探地雷达数据可视化方法，其特征是，所述步骤S133中，计算R的方法为：

$R=\frac{\left|{\mathrm{DIS}}_A\right|}{\left|{\mathrm{DIS}}_B\right|};$

其中，DIS_A表示距离向量DIS中小于分割点dt的各距离集合，DIS_B表示距离向量DIS中大于分割点dt的各距离集合，〡DIS_A〡表示小于分割点dt的距离个数,〡DIS_B〡表示大于分割点dt的距离个数。

5.如权利要求3所述的一种探地雷达数据可视化方法，其特征是，所述步骤S134中R的设定条件为：

R∈(0.2，5)。

6.如权利要求3所述的一种探地雷达数据可视化方法，其特征是，所述步骤S134中，计算gap的方法具体为：

gap＝μ_B-σ_B-(μ_A+σ_B)

其中，μ_B表示大于分割点dt的所有距离的平均值，σ_B表示大于分割点dt的所有距离的方差，μ_A表示小于分割点dt的所有距离的平均值，σ_A表示小于分割点dt的所有距离的方差。

7.如权利要求3所述的一种探地雷达数据可视化方法，其特征是，所述步骤S132中，计算候选序列C(c)与数据集D中各序列距离的距离向量的方法为：

S1321：输入一个候选序列C＝(c₁,c₂,…c_clen)和数据集D＝{T(1),T(2),…,T(n)}；

S1322：初始化距离向量DIS＝Φ；

S1323：利用滑动窗口ql＝clen，clen为当前候选序列C＝(c₁,c₂,…c_clen)的长度，依次取出T(ds)的所有子序列Z＝{Z₁,Z₂,…Z_v}，Z_i＝(z₁,z₂,…,z_clen)；

其中：

T(ds)＝(t₁,t₂,…t_len)，len为数据集T(ds)＝(t₁,t₂,…t_len)的长度；

$1\≤ds\≤n,1\≤i\≤v,v=\frac{len(len+1)}{2};$

S1324：计算候选序列C(c)与子序列集Z中各子序列之间的距离dis，取最小值为候选序列C(c)与数据序列T(ds)的距离dis_ds，更新到DIS中；

S1325：判断ds是否为n，若是，则对DIS排序得DIS＝{dis₁,dis₂,…dis_n}，否则，转向步骤S1323。

8.如权利要求1所述的一种探地雷达数据可视化方法，其特征是，所述步骤(2)的具体方法为：

步骤S201：设定待匹配子序列长度范围[minlen,maxlen]，输入待测的GPR数据集D＝{T(1),T(2),…,T(n)}与局部特征集U＝{U(1),U(2),…,U(m)}；

步骤S202：对于T(ds)，分别初始化起始点st＝1、长度值Mlen＝minlen以及距离向量sdist＝Φ；

步骤S203：初始化i＝1；

步骤S204：计算[minlen,maxlen]范围内，与局部特征U(i)距离d_i最小的子序列M_(st，Mlen)，更新到距离向量sdist；

步骤S205：判断i是否为m，若是，则排序距离向量sdist，否则i＝i+1，转向步骤S204；

步骤S206：取子序列M_(st，Mlen)与局部特征集U中各特征的距离的距离向量中最小的i值，作为M_(st，Mlen)的类别序号；

步骤S207：判断st+Mlen-1是否为len，若是，则下一步；否则更新起始点st＝st+Mlen，并转向步骤S204；len为T(ds)的长度；

步骤S208：判断ds是否为n+1，若是，则下一步；否则返回步骤S202；

步骤S209：输出类别矩阵。

9.如权利要求8所述的一种探地雷达数据可视化方法，其特征是，所述步骤S204中，d_i距离计算公式为：

d_i＝DTW(M_(st,Mlen),U(i))；

其中，st为待匹配序列的起始点，Mlen为[minlen,maxlen]范围内某一长度，U(i)为局部特征，M_(st，Mlen)为与局部特征U(i)距离d_i最小的子序列。