非均质土壤中直埋式电缆的精确定位技术研究

于虹 李昊 周帅 沈锋

摘  要：近年来，探地雷达在地下电缆成像中得到了广泛的应用。然而，由于非均质土壤介电参数呈现随机性和复杂性等特征，使用单个相对介电常数来表征整个成像区域内土壤电学特性将会向成像结果中引入较大误差。如何精确表征土壤的介质特性，进而实现对直埋式电缆的精准定位仍然是亟待解决的问题。为此提出了一种改进的方法对地下管线类目标进行重建。该方法基于列聚类算法以及改进的Hough变换技术，提取了各管线目标的双曲线特征，针对每个目标所在区域适应性地构建土壤介电模型，并结合背向投影成像算法分别对单个目标进行重构；最后将各重构目标进行融合，实现对地下多根电缆的精细化重构。与传统方法相比，改进的方法具有更好的准确性与定位精度。

关键词：探地雷达；输电线路；非均质土壤

DOI：10.15938/j.jhust.2023.05.005

中图分类号： TP751

文献标志码： A

文章编号： 1007-2683（2023）05-0034-08

Direct Burial Cable Reconstruction Framework for

Accurate Positioning in Heterogeneous Soil

YU Hong1，  LI Hao1，  ZHOU Shuai1，  SHEN Feng2

（1.Electric Power Research Institute， Yunnan Power Grid Co.， Ltd.， Kunming 650217， China;

2.School of Instrumentation Science and Engineering， Harbin Institute of Technology， Harbin 150006， China）

Abstract：In recent years， GPR has been widely used in underground cable imaging. However， due to the randomness and complexity of dielectric parameters of heterogeneous soil， using a single relative dielectric constant to characterize the electrical characteristics of soil in the whole imaging area will introduce large errors into the imaging results. How to accurately characterize the medium characteristics of soil and realize the accurate positioning of directly buried cables is still an urgent problem to be solved. Therefore， an improved method is proposed to reconstruct underground pipeline targets. Based on the column clustering algorithm and the improved Hough transform technology， the hyperbolic characteristics of each pipeline target are extracted， the soil dielectric model is adaptively constructed for each target area， and the single target is reconstructed combined with the back projection imaging algorithm. Finally， the reconstruction targets are fused to realize the fine reconstruction of multiple underground cables. Compared with the traditional method， the improved method has better accuracy and positioning accuracy.

Keywords：GPR; transmission line; heterogeneous soil

收稿日期： 2022-05-24

基金项目： 国家自然科学基金（61673128）.

作者简介：

于  虹（1978—），女，博士，正高级工程师；

李  昊（1986—），男，硕士.

通信作者：

沈  锋（1981—），男，教授，博士研究生导师，E-mail：sf407@126.com.

0  引  言

電力电缆已经逐步取代架空输电线路成为城市电网体系中的主要传输方式，然而随着国内旧城区改造工作的不断推进，出现了大量相关老旧线路敷设信息丢失的问题，现有探测手段仍无法对线路进行可视化无损检测，导致政策推进过程中出现大量重复施工，极大增加了线路改造工作难度，严重阻碍了现代化城市的快速发展。近年来，基于探地雷达法的地下电缆无损检测技术受到了广泛关注。该方法利用电磁（EM）场检测介质相对介电常数的变化，已被广泛应用于地下目标结构参数估计，例如目标直径［1-2］、分布方向［3］和目标数量等［4-6］。

为了获得地下管网的精确位置信息，通常需要对探地雷达获取的原始数据进行精细化处理。具体来说包含以下几个方面：在数据预处理方面，常用的方法有时间零点矫正、带通滤波以及增益补偿等。其中，时间零点校正用于将空地界面反射的位置固定为公共时间零点位置［7］；带通滤波通过滤除频率超出工作带宽的信号分量来提高信噪比［8］；增益补偿功能用于辅助数据显示和解释［9］。对于收发天线间的直接耦合干扰，可以通过背景去除方法，如均值减法［10］、主成分分析［11-12］和奇异值分解（SVD）［13］，有效抑制由天线和地面反射之间的直接耦合而产生的水平杂波。在土壤的相对介电常数估计方面，通常采用双曲线拟合方法［14］或使用电介质探针计算电磁波在土壤中的波速来获得。

传统的目标检测方法通常将一个特定位置的测量或估计介电常数视为整个勘测现场中土壤的常数，使得电磁信号在整个探测区域内具有相同的波速［15-16］。一些广泛使用的商用探地雷达系统，如IDS Geo雷达［17］和GSSI系统［18］等会要求用户输入一个恒定的介质介电常数值。但在实际工况条件下，尤其是在湿度高的热带地区，土壤自身的非均质性导致电磁波速在土壤不同部位具有显著差异，从而严重影响后期地下目标定位精度。Gonzalez-Huici和G. Gennarelli等研究了土壤结构中不同方向上介电常数变化规律［18-19］，其中研究人员引入了一个等效介电常数函数，通过大量的解析计算初步得到了热带土壤的介电常数分布規律。在此基础上，A. Bitri等中进一步完成了水平层状非均匀介质中二维（2-D）探地雷达数据预处理，并实现了两条管道的二维成像定位［20］。然而，上述方法只能对地下空间中的稀疏分布目标进行精确定位，随着目标间距离逐渐缩小，传统方法的精度也急剧下降。根据N. R. Peplinski等的研究成果可以发现，由于含水量、成分和密度的不同，同一区域内不同位置的土壤相对介电常数可能会发生剧烈变化［21-22］，使用不匹配的介电常数会将导致对地下空间的错误重建。因此，采用传统方法使用相同的相对介电常数处理探地雷达数据可能会导致对目标水平和垂直位置的不准确估计。此外，土壤的异质性会造成干扰，从而掩盖紧密排布的地下管线反射信号。因此，探索一种有效、准确的方法来提高地下管线检测精度，最大限度地减少杂波噪声，在非均匀土壤环境下精确定位紧密排布的地下管线位置信息。

为此，本文提出了一种改进的方法对地下管线类目标进行重建，并使用天线和直埋式线缆之间的介质的垂直等效介电常数来计算电磁波波速。首先，通过列聚类（C3）算法对二维距离剖面（B-Scan）进行预处理，以提取双曲线［23-26］。在这个过程中，代表不同管线的双曲线被分离并提取到新的B-Scan中，同时去除随机杂波噪声。接下来，使用改进的Hough变换（IHT）方法来估计每个提取出来的B-Scan中双曲线所代表的目标位置以及垂直方向上等效土壤介电常数。在此基础上，通过优化后的介电常数可以获得每根管线的重建图像。最后，将上述重建结果融合，即可得到非均匀土壤环境下各管线目标的精确位置信息。

1  地下电缆高精度定位方法

传统的地下目标检测方法和本文提出的地下目标检测方法流程分别如图1（a）和（b）所示。可以发现，两种方法均以预处理后的B-Scan信号作为输入，其中预处理包括时间零点校正、带通滤波、SVD去除背景和自适应阈值。然而两种方法的不同点在于，传统方法仅任选一个双曲线作为参考来提取土壤介电常数，且不考虑勘探测区域内不同位置的土壤介电常数变化，提取到的介电常数将被直接用于处理整个探测区域中的B-Scan信号。

而在本文提出的方法中，首先结合聚类算法将各双曲线特征进行提取，并分别估计每一根双曲线所在区域土壤的等效介电常数，用于独立对每个目标进行重建。需要强调的是，由于土壤自身介电特性在不同的区域内具有显著差异，每个目标所对应的等效介电常数也不相同。这一步考虑了土壤的非均质性，从而能够显著提高对于直埋式电缆的定位精度。最后，将所有重建的目标进行合成，并最终得到优化后的重建目标信息。

1.1  基于聚类算法的双曲线分离提取技术

为了从B-Scan中提取双曲线特征，从而单独分析提取的双曲线，以减少双曲线相交区域和环境噪声的干扰，本文采用了列聚类（C3）算法对B-Scan中的双曲线信息进行提取。以图2所示的B-Scan为例，采用C3算法能够在提取双曲线1和2的同时忽略噪声1和2。

为了便于描述，将每列中的一组相邻像素点称为一个“区块”。对于像素值小于阈值的像素点所构成的区块被称为“噪声”，例如噪声1（包含2个像素点）和2（包含5像素点）。通常来说，阈值s取决于传感器的噪声级、雷达中心频率和采样频率。s的最大值与采样频率fs成正比，与中心频率fc成反比。为例抑制大部分噪声，一个理想的阈值s应大于大部分噪声中的像素数值，但应小于kfs/fc（k为常数）。

当收发天线从第一列C1扫描到最后一列C32。相邻列中的相邻区块按照下列步骤进行分类：

1）每一个新区域起始于前一列中没有相邻区块的区块，即C2中的红色区块和C15中的蓝色区块。

2）当一个区块与任一区域中的区块相邻时，将该区块添加进此区域中，例如C10中的红色区块需要被添加到从C2到C9所属的红色区域中。

3）若一列中的某个区块被多个区域共享（例如C14中的蓝色和红色区域共享了C13中的绿色区块），则该区块将被分成两个区块以保证分类过程的连续性。

4）当下一列中没有像素相邻时，区域停止延伸（例如，C15中的红色区块和C22中的蓝色区块）

此外，在算法中还应再设置一个额外阈值，当某一区块和现有区域之间的像素值小于阈值时，该区域应停止扩展。综上，在分类完成后，记录每个区域中像素的坐标，通过提取原始B-Scan中的像素坐标和相应值，实现各双曲线的分离并提取。

1.2  各目标土壤等效介电常数的估算

在提取到代表直埋式线缆反射的双曲线后，论文使用改进的Hough变换（IHT）估计每个目标上方的土壤等效介电常数。如图3所示，对于提取到的单目标B-Scan，使用改进的Hough变换基于三参数抛物线方程对双曲线进行拟和，从而估算出土壤的相对介电常数以及物体的深度和水平位置。三参数抛物线方程可以表示为：

t2a2－（x－x0）2b2=1（1）

其中：

a=t0

b=t02cεr（2）

式中：t為信号的双向传播时间；x为沿扫描轨迹的水平位置。（x0，t0）为双曲线的中心；c为自由空间中的光速；εr为土壤的相对介电常数。

当波传播至空气和土壤以及土壤和直埋式线缆之间的分界面时，会反射回一个脉冲信号。反映在每一个A-Scan信号剖面中，则是会在远处产生一个如图4所示的一个小信号峰。为了保证计算过程的一致性，文章统一规定将每个回波信号的峰值所在时刻作为目标回波的参考时延，而将第一个回波信号的波谷作为参考时间的零点。为了求解式（1）中的未知参数a、b和x0，在图3（b）中的区域矩阵中随机选择三个不同的列，并在选择对应第k列中最大像素值所在行tk，从而能够得到该列最大像素值所在点坐标（xk， tk）。将三列中的坐标带入式（1）可求得一组双曲线参数（a，b，x0）。在此基础上，多次重复上述过程，可以计算出多组参数。其中，每重复一次上述过程，则需要求取三个双曲线参数中任意一参数当前的平均值。当平均值的变化小于0.1%时，停止重复实验，使用参数的平均值（am，bm，x0m）对双曲线进行重构，并结合目标的深度和水平位置来提取土壤等效介电常数。

1.3  基于背向投影算法的目标重建

虽然通过双曲线拟合的输出可用于确定目标的中心，但为了更准确的实现地下目标结构特征重建，需要引入背向投影成像方法（BPA）。GPR标准BP 算法的基本原理为时延-求和，然而在各道时延位置处的散射响应包含了所有成像点在此道时延位置处的响应和，即不仅包含了待成像点的响应，还有其他成像点的干扰。为了解决这一问题，本文基于互相关法对传统成像算法进行改进，利用各道回波延迟位置处信号与中心道延迟位置处信号的相关系数作为BP算法的加权系数，以实现对噪声和干扰的有效抑制。

具体来说，如图5所示，在基于互相关的后向投影成像算法（CBP）中，若合成孔径位置数目为Np，再对地下某一成像点进行成像时，需要先计算该点到各个合成孔径位置的时延：{τ1，τ2，…，τNp}，在各个时延处位置为中心，上下可以取一段有效信号作为发射窄脉冲段，即主要能量集中的部分，以合成孔径中间位置对应时延位置处的回波信号段为参考信号部分，将各段信号与其做互相关处理，得到互相关信号系数，其中此段信号的长度可以取：

L=2Fs2f0+1=2S+1（3）

其中：L为所取信号段长度，且为奇数，在延时位置处上下各取S个数据，再加上延时位置处数据组成回波段向量，符号“［·］”表示取整符号；Fs表示等效采样频率；f0表示发射信号中心频率。相应回波段的散射响应幅值向量可以表示为

xk=sk（［τk·Fs］－S∶［τk·Fs］+S），k=1，2，…，Np（4）

式中：xk表示第k个天线位置处回波段的散射响应幅值向量；sk表示第k个天线位置处录取的回波向量。若中间取样值超出了采样点数范围M，或者小于0，那么用0来补齐。通过该步骤生成一个关于成像点的L×Np维的矩阵［x1， x2，…，xNp］为参考回波段向量，即中间位置的散射响应值。若S=0，则得到的矩阵即为1×Np维的矩阵，与标准的BP 算法相同。

为了降低旁瓣和干扰的能量，本文根据各个合成孔径位置得到的回波段向量与参考向量的相关性，计算得到一组加权系数向量ρ={ρ1， ρ2，…， ρNp} ，其中相关系数的计算采用皮尔逊相关系数：

ρ=ρk=Cov（xk，xNp/2）Cov（xk，xk）·Cov（xNp/2，xNp/2）（5）

其中Cov（i，j）函数表示向量i与向量j的协方差。

最后通过系数加权叠加，可以得到成像点的幅值：

E=∑Nρk=1|ρk|·xk（S+1）（6）

2  仿真实验场景设计

为了更好地模拟真是环境下土壤结构的随机分布特性，论文基于时域有限差分法结合gprMax软件［27］对异质土壤环境中的直埋式管线的检测过程进行了仿真。环境中土壤结构单元的长度为0.002m。为了模拟直埋式电缆在地下的分布，如图6所示，设置了多个不同半径的线缆埋在非均质土壤的不同深度。5根电缆沿x轴从左到右，其直径分别为15cm（r1），10cm（r2），20cm（r3），10cm（r4），20cm（r5）；深度分别为0.5m（r1），0.5m（r2），0.7m（r3），0.5m（r4），0.7m（r5）。

在此基础上，本文基于Peplinski模型构造了非均质土壤仿真实验场景。在该场景中，土壤的砂粒分数（S）为0.3，黏土分数（C）为0.7，砂粒密度（ρs）为2.66g/cm3，容重（ρb）为2g/cm3，水体积分数（fu）为0.01至0.15。土壤介电常数数量和分型维度（β）分别设定为20和1.5。在此基础上，本文所提出的改进方法估算了每根线缆所处土壤空间的等效介电常数值，并列出和比较所提出方法和传统方法对管线深度和水平位置的检测误差，结果表明，改进的方法具有更好的性能。土壤场景的侧视图如图6所示。

上述仿真实验场景中采用了双站式探地雷达信号收发架构，其中收发天线的极化方向管线延伸方向相同，收发天线相距0.10m，位于土壤表面上方0.05m处。激励波形为中心频率为1GHz的Ricker脉冲。收发天线沿x方向的轨迹移动（见图6），步进为0.01m。其中r1到r5是代表5根管线的圆柱体。r4和r5附近土壤的相对介电常数与其他3根管线周围土壤的相对介电常数显著不同（用红色标识其轮廓）。收集180个A-Scan并合并成一个B-Scan。土壤相对介电常数分布在6到8的范围内。

3  仿真实验结果分析

图7～9展示了对图6所示目标的具体处理流程。整个流程共分为3大部分：A-Scan零点校准、B-Scan双曲线提取以及管线特征重构。

在A-Scan零点校准部分，如图7所示，通过设置地面反射的负最大峰值点作为时间零点位置，对A-Scan信号进行矫正。在图7（a）中，仅显示了在没有地下物体掩埋的情况下，由天线到天线分离和地面反射引起的振荡。在图7（b）中，由天线到天线分离和地面反射引起的振蕩与图7（a）中的相同点，还有另一个由地下目标反射引起的反射。

图8（a）和（b）分别显示了时间零点校正的原始B-Scan和使用基于SVD的背景去除方法和自适应阈值技术的预处理B-Scan。如图8（a）所示，经过时间零点校正后，两个偶极子之间的直接耦合和空气/土壤界面引起的反射在原始数据中具有较大的振幅，形成了明亮的水平带状。使目标反射难以区分。如图8（b）所示，去除背景后，收发天线间的直接耦合以及地表反射回波得到了抑制。在这种情况下，SVD获得了令人满意的背景去除结果。由于土壤的不均匀性和湿润环境中电磁波的高衰减，来自目标的反射与背景噪声融合在一起。

因此，本文使用自适应阈值技术进一步降低噪声。之后，显示了代表五根线缆反射的五条相交双曲线，但如图8（c）所示仍有一些噪声干扰。

在图6所示的非均质土壤中，目标r4和r5周围土壤的相对介电常数相比于其他三个目标具有显著差异。在传统方法中，使用单个相对介电常数值对整个测量区域内的对象进行成像，可能会导致不准确的检测结果。为了克服这一局限性，本文使用列聚类方法分别提取五条双曲线。从图8（c）中提取的单条双曲线如图9所示。此外，通过在双曲线提取过程中引入自适应二值化算法可以进一步降低背景噪声，产生了具有清晰物体反射的B-Scan。

最后，对比如图10所示的重构结果，如表1所示，可以发现改进的方法信噪比分布在11dB左右，传统方法的信噪比分布在5dB左右，其中由于介质分布不均匀，目标R4和R5的信噪比要略低于目标R1～R3。可以发现改进的方法相比于传统方法具有更小的旁瓣以及更高的信噪比。

此外，在定位精度方面，如表2所示，相比于传统方法，改进方法具备更低的探测误差，对管线目标具备更高的定位精度。

4  结  论

本文提出了一种非均匀土壤环境下直埋式线缆探地雷达数据的信号处理方法，通过考虑土壤介电常数的变化来改善二维重建结果。具体来说，首先基于列聚类（C3）算法的双曲线分离和提取方法，以获得单独的双曲线，减少土壤异质性造成的杂波；随后，采用IHT技术估算了每条双曲线的土壤等效介电常数，并结合背向投影算法完成了对电缆截面形状的重构；最后，将所有聚焦图像合并成一幅图像，以恢复地下线缆的敷设场景。通过该方法能够有效改善地下管线检测精度。

参 考 文 献：

［1］  J. R. Butnor， J. A. Doolittle， L. Kress， et al. Use of Ground-penetrating Radar to Study Tree Roots in the Southeastern United States［J］. Tree Physiology， 2001， 21（17）： 1269.

［2］  DANNOURA M， HIRANO Y， IGARASHI T， et al. Detection of Cryptomeria Japonica Roots with Ground Penetrating Radar［J］. Plant Biosystems， 2008， 142（2）： 375.

［3］  SUN H H， LEE Y H， LUO W H， et al. Dual-Cross-Polarized GPR Measurement Method for Detection and Orientation Estimation of Shallowly Buried Elongated Object［J］. IEEE Transations on Instrumentation and Measurement， 2021： 70.

［4］  JAISWAL A K， THAKUR P K， KUMAR P， et al. Geospatial Modeling of Water Supply Distribution System： A Case Study of Dehradun City， India［J］. H2OPEN Journal， 2021， 4（1）： 393.

［5］  SHEN Y P， LIN Y R， LI P， et al. Simulation and Detection Leakage of Underground Water Pipeline by Ground Penetrating Radar［J］. Journal  of Testing and Evaluation， 2020， 48（3）： 2003.

［6］  JIN Y， DUAN Y. Wavelet Scattering Network-Based Machine Learning for Ground Penetrating Radar Imaging： Application in Pipeline Identification［J］. Remote Sensing， 2020， 12（21）：3655.

［7］  YELF R. Where is True Time Zero？［C］// Proceedings of the Tenth International Conference on Grounds Penetrating Radar， 2004： 279.

［8］  刘海，黄肇刚，岳云鹏，等. 地下管线渗漏环境下探地雷达信号特征分析［J］. 电子与信息学报，2022， 44（4）： 1257.

LIU Hai， HUANG Zhaogang， YUE Yunpeng， et al. Characteristics Analysis of Ground Penetrating Radar Signals for Groundwater Pipe Leakage Environment［J］. Journal of Electronics & Information Technology， 2022， 44（4）： 1257.

［9］  朱兆荣，赵守全，秦欣，等. 探地雷达在铁路隧道衬砌质量无损检测中的应用研究［J］. 工程地球物理学报. 2021， 18（5）： 703.

ZHU Zhaorong， ZHAO Shouquan， QING Xin， et al. Application of Ground Penetrating Radar in Nondestructive Testing of Railway Tunnel Lining Quality［J］. Chinese Journal of Engineering Geophysics， 2021， 18（5）： 703.

［10］PERSICO R， SOLDOVIERI F. Effects of Background Removal in Linear Inverse Scattering［J］. IEEE Transactions on Geoscience and Remote Sensing， 2008， 46（4）： 1104.

［11］CHEN G Y， QIAN S E. Denoising of Hyperspectral Imagery Using Principal Component Analysis and Wavelet Shrinkage［J］. IEEE Transactions on Geoscien and Remote Sensing， 2011， 49（3）： 973.

［12］薛策文，冯晅，李晓天，等. 全极化探地雷达多极化数据融合分析研究［J］. 雷达学报， 2021， 10（1）： 74.

XUE Cewen， FENG Xuan， LI Xiaotian， et al. Multi-polarization Data Fusion Analysis of Full-polarimetric Ground Penetrating Radar［J］. Journal of Radars， 2021， 10（1）： 74.

［13］鄭俊，刘婉萍，马念茹，等. 探地雷达信号杂波抑制［J］. 电子测量技术， 2018， 41（21）： 63.

ZHENG Jun， LIU Wanping， MA Mianru， et al. Clutter Reduction for Ground Penetrating Radar Signal［J］. Electronic Measurement Technology， 2018， 41（21）： 63.

［14］MERTENS L， PERSICO R， MATERA L， et al. Automated Detection of Reflection Hyperbolas in Complex GPR Images With No A Priori Knowledge on the Medium［J］. IEEE Transactions on Geoscience and Remote Sensing， 2016， 54（1）： 580.

［15］LANTINI L， TOSTI F， GIANNAKIS I， et al. An Enhanced Data Processing Framework for Mapping Tree Root Systems Using Ground Penetrating Radar［J］. Remote Sensing， 2020， 12（20）： 15.

［16］ABOUDOURIB A， SERHIR M， LESSELIER D. A Processing Framework for Tree-Root Reconstruction Using Ground-Penetrating Radar Under Heterogeneous Soil Conditions［J］. IEEE Transactions on Geoscience and Remote Sensing， 2021， 59（1）： 208.

［17］LAMERI S， LOMBARDI F， BESTAGINI P， et al. Landmine Detection from GPR Data Using Convolutional Neural Networks［C］// 2017 25th European Signal Processing Conference （Eusipco）. 25th European Signal Processing Conference （Eusipco）， 2017： 508.

［18］GONZALEZ-HUICI M A， CATAPANO I， SOLDOVIERI F. A Comparative Study of GPR Reconstruction Approaches for Landmine Detection［J］. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing， 2014， 7（12）： 4869.

［19］GENNARELLI G， CATAPANO I， EEROBERT X， et al. A Ground Penetrating Radar Imaging Approach for a Heterogeneous Subsoil With a Vertical Permittivity Gradient［J］. IEEE Transactions on Geoscience and Remote Sensing， 2021， 59（7）： 5698.

［20］BITRI， GRANDJEAN G. Frequency-wavenumber Modelling and Migration of 2D GPR Data in Moderately Heterogeneous Dispersive Medial［J］. Geophysical Prospecting， 1998， 46（3）： 287.

［21］DOBSON M C， ULABY F T， HALLIKAINEN M T， et al. Microwave Dielectric Behavior of Wet Soil .2. Dielectric Mixing Models［J］. IEEE Transactions on Geoscience and Remote Sensing， 1985， 23（1）： 35.

［22］PEPLINSKI N R， ULABY F T， DOBSON M C. Dielectric Properties of Soils in the 0.3-1.3-GHz Range （Vol 33， pg 803， 1995）［J］. IEEE Transactions on Geoscience and Remote Sensing， 1995， 33（6）： 1340.

［23］DOU Q X， WEI L J， MAGEE D R， et al. Real-Time Hyperbola Recognition and Fitting in GPR Data［J］. IEEE Transactions on Geoscience and Remote Sensing， 2017， 55（1）： 51.

［24］何云斌，董 恒，萬 静.移动型数据与静态型数据的混合聚类算法［J］.哈尔滨理工大学学报，2021，26（2）： 26.

HE Yunbin， DONG Heng， WAN Jing. Hybrid Clustering Algorithm for Mobile Data and Static Data［J］.Journal of Harbin University of Science and Technology，2021，26（2）：26.

［25］万静，郑龙君，何云斌，等.高维数据的高密度子空间聚类算法［J］.哈尔滨理工大学学报，2020，25（4）： 79.

WAN Jing， ZHENG Longjun， HE Yunbin， et al. High Density Subspace Clustering Algorithm for High Dimensional Data［J］.Journal of Harbin University of Science and Technology，2020，25（4）：79.

［26］何云斌，孙暖，万静，等. 角度度量的动态增量聚类算法［J］.哈尔滨理工大学学报，2019，24（6）： 109.

HE Yunbin， SUN Nuan， WAN Jing， et al. A Kind of Dynamic Clustering Algorithm Based on Angle Measurement［J］. Journal of Harbin University of Science and Technology，2019，24（6）： 109.

［27］WARREN C， GIANNOPOULOS A， GIANNAKIS I. gprMax： Open Source Software to Simulate Electromagnetic Wave Propagation for Ground Penetrating Radar［J］. Computer Physics Communications， 2016， 209： 163.

（编辑：温泽宇）