一种探地雷达数据可视化方法
【技术领域】
[0001] 本发明属于计算机数据处理可视化领域，具体涉及一种基于多尺度局部特征和动 态时间规整的探地雷达数据可视化方法。
【背景技术】
[0002] 近年来，探地雷达（Ground Penetrating Radar，简称GPR)作为一种快速、连续、非 接触电磁波探测技术，以其采集速度快、分辨率高的特点，成为物探领域研究的热点之一。 在实际的应用中，GPR数据的解译需要较高的专业知识及经验，而且环境以及系统的干扰使 得GPR数据掺杂各种杂波和噪声以及个人经验不同使得结果存在差异性，制约了从GPR技 术的进一步发展和实用化。其中GPR数据解译过程中对其特征的提取传统方法本质都是基 于线性预测、Fourier分析或小波分析这样的传统信号分析技术，对于地下介质类别数目就 只得依赖解译人员专业知识及经验，过度或者不当的处理，都会使数据失去原有真实性。
[0003] 数据可视化是描述和表达各种半结构化甚至非结构化问题的关系和模型的最佳 方法和手段，可以成为解译GPR数据中蕴含信息的新的解决途径。通过可视化的方法对GPR 数据进行图像显示，提供新的洞察力以挖掘传统方法提取不到的信息。
[0004] 然而由于GPR数据采用声波反射测速技术所得，反射回波的不确定性和地质的多 样性，使得GPR数据富含杂波、噪音和存在不连续性，以及获取的数据分辨率相对低的这些 特点，给其可视化带来了极其大的困难与挑战。
【发明内容】
[0005] 本发明的目的就是为了解决上述问题，提出了一种探地雷达数据可视化方法，该 方法基于多尺度局部特征和动态时间规整，从而客观的挖掘更多有效特征信息，解决了 GPR 数据解译过程中过于依赖人工因素的问题。
[0006] 为了解决上述问题，本发明采用如下技术方案：
[0007] -种探地雷达数据可视化方法，包括以下步骤：
[0008] (1)选定GPR数据集，分别提取GPR数据集不同尺度的无监督局部特征，得到GPR 数据的局部特征集U= {U(l)，U(2)，"、U(m)};
[0009] (2)设定待匹配子序列长度范围[minlen，maxlen]，利用动态时间归整方法将待 测的GPR数据集D= {T(l)，T(2)，…，T(n)}与步骤（1)中所得局部特征集U进行模式匹配 分类，U中的m即为此处分类数目；
[0010] (3)用可视化方法中的颜色映射方法将步骤（2)中进行模式匹配的结果进行展 示，得到GPR数据的直观图像表示。
[0011] 所述步骤（1)中提取GPR数据的不同尺度的无监督局部特征的方法为：
[0012] 步骤SllO :选取待测GPR数据集中某一序列T(ds) = (^t2，…t J，其中 I < ds < n，初始化无监督局部特征数据集U = Φ ;
[0013] 步骤S120 :设定滑动窗口范围ql，移动步长1，从序列T(ds)中生成多个尺度的候 选无监督局部特征集C ;
[0014] 步骤S130 :计算候选集C中每一个候选序列的质量度量gap，并对其排序；
[0015] 步骤S140 :取质量度量gap最大的候选序列，更新到无监督局部特征数据集U中 作为U(i);
[0016] 步骤S150 :计算U(i)与GPR数据集中各序列的距离，设定阈值Θ，移除距离小于 阈值Θ的序列，判断DIS^否稳定，若稳定，则输出特征集U ;否则取距离最大的序列，返回 步骤Sl 10。
[0017] 所述步骤S150中，
[0018] 判断DISa稳定的条件具体为：
[0019] |disa|=i;
[0020] 即：当小于分割点dt的距离的个数等于1达到稳定。
[0021] 计算候选集C中每一个候选序列的质量度量gap的方法为：
[0022] S131 :输入候选序列C(c)和数据集D = {T(l), T(2),…，T(n)}，初始化maxGap = 〇 ；
[0023] S132 :计算候选序列C(c)与数据集D中各序列的距离，并将其排序，得到候选序列 C(c)与数据集D中各序列距离的距离向量DIS= Wis1, dis2, ."(IisrJ ;
[0024] S133 :根据距尚向量可得n-1个分割点dt，对每一个分割点dt分别计算分割点dt 两端的距离个数的比值R ;
[0025] S134 :判断R是否满足设定条件，若是则计算gap，否则计算下一分割点的R ;
[0026] S135 :判断计算的gap是否大于maxGap，是则更新maxGap，否则计算下一分割点 gap 〇
[0027] 所述步骤S133中，计算R的方法为：
[0029] 其中，DISa表示距离向量DIS中小于分割点dt的各距离集合，DIS B表示距离向量 DIS中大于分割点dt的各距离集合，丨DISa丨表示小于分割点dt的距离个数，丨DISb丨 表示大于分割点dt的距离个数。
[0030] 所述步骤S134中R的设定条件为：
[0031] R e (〇· 2,5) 〇
[0032] 所述步骤S134中，计算gap的方法具体为：
[0033] gap = μ Β- σ Β- ( μ Α+ σ B)
[0034] 其中，以8表示大于分割点dt的所有距离的平均值，σ 8表示大于分割点dt的所 有距离的方差，以&表示小于分割点dt的所有距离的平均值，σ a表示小于分割点dt的所 有距离的方差。
[0035] 所述步骤S132中，
[0036] 计算候选序列C(c)与数据集D中各序列距离的距离向量的方法为：
[0037] 31321:输入一个候选序列0=((：1，(^"(^")和数据集0={1(1)，1'(2)，··· ，T(n)};
[0038] S1322 :初始化距离向量DIS = Φ ;
[0039] S1323 :利用滑动窗口 ql = clen，clen为当前候选序列C = (C1, c2, "YtllfJ的长 度，依次取出T (ds)的所有子序列Z= (Z1, Z2，…Zv} ,Zi= (Z^Z2,…，ZclfJ ;
[0040] 其中：
[0041] T(ds) = U1, t2,…tlen)，Ien 为数据集 T(ds) = U1, t2,…tlen)的长度；
[0043] S1324 :计算候选序列C(C)与子序列集Z中各子序列之间的距离dis，取最小值为 候选序列C(C)与数据序列T (ds)的距离diSds，更新到DIS中；
[0044] S1325 :判断ds是否为n，若是，则对DIS排序得DIS = {diSl，Clis2，…disn}，否则， 转向步骤S1323。
[0045] 所述步骤（2)的具体方法为：
[0046] 步骤S201 :设定待匹配子序列长度范围[minlen，maxlen]，输入待测的GPR数据集 D = {T(1)，T(2)，…，T(n)}与局部特征集 U = {im)，U(2)，…，U(m)};
[0047] 步骤S202 :对于T(ds)，分别初始化起始点st = 1、长度值Mien = minlen以及距 离向量sdist = Φ ;
[0048] 步骤S203 :初始化i = 1 ;
[0049] 步骤S204 :计算[minlen，maxlen]范围内，与局部特征U(i)距离(Ii最小的子序列 M(st, Mlen) ' 更新到距离向量sdist ;
[0050] 步骤S205 :判断i是否为m，若是，则排序距离向量sdist，否则i = i+Ι，转向步骤 S204 ；
[0051] 步骤S206 :取子序列与局部特征集U中各特征的距离的距离向量中最小 的i值，作为M03UM的类别序号；
[0052] 步骤S207 :判断st+Mlen-Ι是否为len，若是，则下一步；否则更新起始点st = st+Mlen，并转向步骤S204 ;len为T(ds)的长度；
[0053] 步骤S208 :判断ds是否为n+1，若是，则下一步；否则返回步骤S202 ;
[0054] 步骤S209 :输出类别矩阵。
[0055] 所述步骤S204中，山距离计算公式为：
[0056] Cl1= DTff(M (stiMlen),U(i));
[0057] 其中，st为待匹配序列的起始点，Mlen为[minlen, maxlen]范围内某一长度，U(i) 为局部特征，为与局部特征U(i)距离d 1最小的子序列。
[0058] 本发明的有益效果：
[0059] 本发明对于GPR数据应用性强，并且能够更好的提取有效GPR数据特征，更为客观 的呈现了 GPR数据所蕴含的信息，用可视化方法直观的体现了探测结果。
[0060] 通过对数据集的多尺度局部特征的提取，解决了未知探测区域地下介质类目问 题，区分出了不同的地下介质，克服了传统GPR数据介意过程中特征提取方法所带来的多 解性的缺点。
[0061] 通过采用DTW距离计算能够对不等长序列进行相似的度量；并且DTW对于序列的 突变或者异常点并不敏感，这对于GPR数据中异常的噪点可以起到忽略的作用。
【附图说明】
[0062] 图1所示为本发明探地雷达数据可视化方法流程图。
【具体实施方式】
[0063] 下文将结合具体实施例详细描述本发明。应当注意的是，下述实施例中描述的技 术特征或者技术特征的组合不应当被认为是孤立的，它们可以被相互组合从而达到更好的 技术效果。
[0064] 如图1所示，本发明提供的一种基于多尺度局部特征和动态时间规整的探地雷达 数据可视化方法包括如下步骤：
[0065] 步骤SlOO:对GPR数据提取多尺度的局部特征；
[0066] 步骤S200:利用DTW距离对SlOO所述特征模型与GPR数据进行模式匹配；
[0067] 步骤S300:对S200所述模式匹配结果运用可视化方法得到GPR数据的直观图像 表不。
[0068] 上述步骤SlOO:对GPR数据提取多尺度的局部特征，具体地包括如下步骤：
[0069] 步骤Sl 10:输入GPR数据集中某一序列T (ds)，初始化无监督局部特征数据集U = Φ 〇
[0070] 步骤S120:设定滑动窗口范围ql，移动步长1，从T (ds)中生成多个尺度的候选无 监督局部特征集C。
[0071] 步骤S130:计算候选集C中每一个候选序列的质量度量gap，并对其排序。其中质 量度量gap计算公式如下：
[0072] gap = μ Β- σ Β- ( μ Α+ σ Β)
[0073] 步骤S140:取最大gap的候选序列，更新到无监督局部特征数据集U中作为U (i)。
[0074] 这样提出去的是某一个尺度的无监督局部特征的提取后，我们便提取到某一道 GPR数据序列T (ds)中的某一位置的一个无监督局部特征，表明此道GPR数据序列中的此位 置上的介质特征。然而GPR数据中可能包含更多的介质信息特征，因而我们对GPR数据中 含有此类介质特征相似的数据序列去除处理后，对数据集中剩余数据序列进行下一局部特 征的查找，其尺度取决于此特征的gap度量大小，因此最后我们会提取到多个尺度的局部 特征。
[0075] 步骤S150:计算U(i)与GPR数据集D中各序列距离，移除距离小于阈值Θ的序 列，判断DIS^否稳定，若稳定，则输出特征集U，否则取距离最大的序列，重复步骤S110。
[0076] 该稳定条件为：
[0077] DIS A| = 1
[0078] 即：当小于分割点dt的距离的个数等于1达到稳定。
[0079] 上述步骤S130:计算候选集C中每一个候选序列的质量度量gap，具体地包括如下 步骤：
[0080] S131:输入候选序列C(c)和数据集D = {T(l), T(2),…，T(n)}，初始化maxGap = 〇 ；
[0081] S1