CN118130513A - 一种基于ct扫描与数值模拟的仿真gpr数据获取方法 - Google Patents

Abstract

本发明公开了一种基于CT扫描与数值模拟的仿真GPR数据获取方法，属于地下病害检测技术领域，针对需检测的不同材料、结构组合的道路进行钻心取样并通过X‑Ray ICT设备扫描构建具有待检测道路特征的多个剖面图，基于道路剖面图及MATLAB构建待检测道路模型并嵌入多类型病害模型，构建出1：1还原的病害滋生的道路模型，最终基于gprMax完成探地雷达正演并获取虚拟探地雷达数据进行分析，实现病害特征的鉴定。

Description

一种基于CT扫描与数值模拟的仿真GPR数据获取方法

技术领域

本发明属于地下病害检测技术领域，具体涉及一种基于CT扫描与数值模拟的仿真GPR数据获取方法。

背景技术

随着运营时间的增加，在荷载及环境因素的综合作用下，道路内部会滋生各类隐蔽病害，倘若不能及时对其进行检测与处治，则会诱发各类路表病害如表面裂缝、车辙、坑槽等，因此道路隐蔽病害的检测重要性极高。当前常用的道路隐蔽病害检测方法包括钻心取样与探地雷达，其中前者属破坏性检测方法，会对路面结构造成不可逆的永久损害且采样点数有限、代表性差。探地雷达凭借其无损性、连续性等优点近年来被逐渐应用于道路隐蔽病害检测，但其在道路领域的应用也存在着以下瓶颈：

1、由于道路材料、结构组合及内部病害的形状、尺寸各不相同，其在探地雷达图像上产生的特征也各不相同，为专业人员识别病害增大了难度。

2、专业人员往往通过钻心取样与雷达分析的经验做出检测判断，主观性较强，且不同的检测人员水平不同，导致了内部病害检测效果较差，只能识别典型病害，而对于部分异常特征难以确定其对应的病害类型。

3、当前的相关标准、规范给出的病害雷达图像特征均是通过大量钻心验证得到的病害特征，该方法首先不能概括全部的病害特征，代表性较差；其次对于不同的路面结构其普适性较差，最后该方法属于有损检测，对路面结构造成损害较大。

4、通过仿真模拟的方法确定雷达图像特征是一种有效的方法，其可以根据实际道路结构组合构建相应的道路及病害模型，通过模拟GPR数据采集过程获得带有病害特征的典型雷达图像特征。但当前雷达正演的方法大部分构建的是均匀路面介质与规则病害，即认为各结构层各处的介电常数值相同，这与实际探地雷达检测环境严重不符。

5、当前基于深度学习的地下病害检测技术逐渐兴起，但训练普适性强的模型需要海量数据，当前由于取芯数量有限，难以获得海量高质量雷达数据训练模型，阻碍了地下病害自动检测技术的发展。

因此，基于探地雷达的道路隐蔽病害检测领域存在着复杂病害特征分析难、高质量模型训练数据获取难的瓶颈。

发明内容

本发明的目的在于克服上述现有技术的缺点，提供一种基于CT扫描与数值模拟的仿真GPR数据获取方法，以解决现有技术难以确定复杂真实道路内部病害雷达特征及深度学习数据量难获取的瓶颈。

为达到上述目的，本发明采用以下技术方案予以实现：

一种基于CT扫描与数值模拟的仿真GPR数据获取方法，包括以下步骤：

步骤1，在待检测路段上钻取芯样，所述芯样由N个结构层组成，获得每一个结构层的厚度，N为≥1的自然数；

步骤2，扫描芯样的若干个纵向剖面，获得若干个剖面扫描图，所述芯样的轴线在纵向剖面上；按照结构层分割剖面扫描图，获得各个结构层的若干幅纵向剖面扫描图，阈值分割每一个结构层的所有纵向剖面扫描图，获得每一个结构层的剖面扫描图的组分图像，所述组分图像包括粘结剂灰度图、空隙灰度图和集料灰度图；将所有的组分图像转化为二值图像，组合后获得所述芯样每一个组分的二值图集合；

步骤3，在MATLAB中构建初始路面结构模型，设置空间步长，对所有的二值图集合进行尺寸变换，获得变换后的组分矩阵，变换后组分矩阵中像素点代表的实际尺寸与路面结构模型中的空间步长相等；根据组分矩阵中各个元素的材料类型赋值介电常数，获得路面结构模型；

步骤4，在路面结构模型中嵌入路面结构病害，获得病害嵌入的路面结构模型；

步骤5，运行gprMax正演程序，通过探地雷达对病害体嵌入的路面结构模型扫描，获得单道波正演数据，进一步的合成探地雷达A-scan图像及B-scan图像。

本发明的进一步改进在于：

优选的，步骤1中，钻取芯样的钻头直径大于100mm，钻头长度为400mm。

优选的，步骤1中，若取出的芯样结构层分离，步骤2中针对分离的每一个分层芯样均获得若干个剖面扫描图；

将同一个芯样的分层芯样同一个剖面的剖面扫描图组合后，获得所述芯样完整的剖面扫描图。

优选的，步骤2中，按照结构层从下到上顺序，组合获得每一个组分的二值图集合。

优选的，步骤3中，所述初始路面结构模型包括空气层和路面结构层；

所述空气层为单个矩阵，所述路面结构层包括i个矩阵，i为路面结构层的层数；所述i个矩阵的垂直网格数由空间垂直方向步长和实际取芯高度获得。

优选的，步骤3中，二值图集合进行尺寸变换的公式为：

其中，l_ori代表二值图对应的真实长度，h_ori代表二值图对应的真实宽度，dx代表长度方向的空间步长，dy代表宽度方向的空间步长。

优选的，步骤4中，所述路面结构病害包括不规则层间粘结不良病害、不规则脱空病害和裂缝病害。

优选的，所述不规则层间粘结不良病害或所述不规则脱空病害的构建过程为：

(1)确定病害的水平起始位置、病害垂直高度和病害水平延伸长度；

(2)在MATLAB中构建尺寸为Lc*1的列向量及循环函数，循环函数根据列向量中索引的值移动，重复遍历Lc个点后获得不规则层间粘结不良或不规则脱空病害。

优选的，步骤4中，所述裂缝病害的构建过程为：

(1)构建垂直裂缝，确定裂缝发生的起始位置、裂缝宽度和裂缝垂直长度；

(2)构建尺寸为H1*1的列向量和循环函数，循环函数根据列向量中索引的值移动，重复遍历H1个点后获得裂缝病害。

优选的，步骤5中，运行单道波堆积图生成程序，根据单道波正演数据获得A-scan图像；基于单道波正演数据，获得反射振幅数组，填充灰度伪彩图，获得B-scan图像。

与现有技术相比，本发明具有以下有益效果：

本发明公开了一种基于CT扫描与数值模拟的仿真GPR数据获取方法，针对需检测的不同材料、结构组合的道路进行钻芯取样并通过CT扫描构建具有待检测道路特征的多个剖面图，基于道路剖面图构建待检测道路模型并嵌入多类型病害模型，构建出1：1还原的病害滋生的道路模型；在整个过程中，基于X-Ray ICT设备扫描及MATLAB建模获取嵌入真实病害的道路仿真模型，最终基于gprMax完成探地雷达正演并获取虚拟探地雷达数据进行分析，最终基于探地雷达正演获得多维探地雷达图像，实现病害特征的鉴定。该方法减少了道路取芯的数量，从一个芯样上能够获得多个剖面，并且通过构建逼近真实的合成雷达图像进行先行病害分析，有助于帮助检测人员获取各类型病害的典型特征，构建定向道路的典型隐蔽病害雷达图谱，提高病害检测准确率；其次该方法也可以为训练深度学习模型提供海量高质量的数据，克服高质量数据获取难的瓶颈，减少取芯验证的频次；通过正演数值模拟方法可以用于确定病害雷达特征，模型能够无限量随机构建，并生成海量数据，解决了当前数值模拟建模缺乏真实性，该方法能够获取逼近真实GPR数据以分析复杂病害的雷达特征、扩充自动检测模型训练数据库。

进一步的，该方法能够布置于道路内部检测工作前期，用于确定待检测道路的内部病害的典型雷达图像特征，提高隐蔽病害可识度，减少隐蔽病害的漏识度，提高道路内部质量评估的全面性。同时，该方法只需在检测前期取芯几处，相比传统单公里取芯1～2个的模式大大减少了取芯成本，提高了检测效率。同时，通过基于真实道路剖面扫描图构建的三维道路模型在嵌入各类型病害后获得的正演图更加接近真实探地雷达图像，可用于指导专业人员获悉病害多类型复杂特征的表现形式，减少取芯验证成本。

附图说明

图1是本发明基于CT扫描与数值模拟的仿真GPR数据获取方法流程示意图；

图2是芯样扫描剖面及组分分割图；

图3是内嵌病害的道路仿真模型；

图4是探地雷达虚拟A-scan图像；

图5是探地雷达虚拟B-scan图像。

具体实施方式

下面结合附图对本发明做进一步详细描述：

在本发明的描述中，需要说明的是，术语“中心”、“上”、“下”、“左”、“右”、“竖直”、“水平”、“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制；术语“第一”、“第二”、“第三”仅用于描述目的，而不能理解为指示或暗示相对重要性；此外，除非另有明确的规定和限定，术语“安装”、“相连”、“连接”应做广义理解，例如，可以是固定连接，也可以是可拆卸连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通。对于本领域的普通技术人员而言，可以具体情况理解上述术语在本发明中的具体含义。

参见图1，本发明公开了一种基于CT扫描与数值模拟的仿真GPR数据获取方法，该方法包括以下步骤：

步骤1，使用混凝土钻孔取芯机在选定点处进行钻芯取样，并进行道路路面结构的结构层信息计算及统计。

步骤1.1，针对待检测道路，根据道路建设资料，选定具有合适孔径及钻芯深度的钻头并选择几处代表性路表无病害处进行取芯，以获取完整的路面结构组合信息。取出的芯样尽量保证各结构层粘结、整体性好，便于后续的CT扫描。

该步骤中，为了更多地捕获道路材料组成、结构组合信息以及方便后期CT扫描成像，混凝土钻孔取芯机钻头直径应大于等于100mm。为了完整获取道路路面结构层信息，钻头可进入长度应超过路面结构总厚度至少50mm，以保证可将路面结构全部取出且不掉钻。

优选地，使用适配场景更广的直径为100mm，长度为400mm钻孔取芯机钻头，对于路面结构层厚度较厚的道路，可使用400mm长钻头分两次取样；而对于路面结构层较薄的道路，可使用400mm长钻头单次取样。

步骤1.2，对于取出的芯样使用直尺测量芯样不同的边线，所述边线为沿芯样厚度方向的直线；每次测量获得所述边线上每一个结构层的厚度，所述结构层包括面层、基层、底基层和垫层等结构层，将测量多个边线每一个结构层的厚度取平均值，记录各结构层材料、厚度等信息，便于后续送入X-Ray ICT平台进行数字化成像；对于不同的路段，结构层包括的类型数量能够不同。

优选的，测量的边线数量为3个，且三个边线沿芯样的周向等分设置。

步骤2，对芯样使用X-Ray ICT设备进行扫描并获取剖面扫描图，对其进行裁剪获取单结构层剖面扫描图；将单结构层剖面扫描图导入Mimics软件进行阈值分割(图2)，获得单组分二值图；

步骤2.1，对芯样进行标记，在芯样的顶端画若干个方向线，每一个方向线为过顶端中心点的方向线；

优选的，使用记号笔以沿72°为一个间隔在芯样顶端进行标记，共作5个方向线，代表沿5个方向线分别进行成像。

步骤2.2，将多个芯样表面清理干净放至X-Ray ICT试验平台进行扫描，设置扫描参数，开启试验平台对试件进行扫描。需要说明的是，对于路面结构厚度＜400mm的道路，其取出的芯样多数粘结良好，呈整体状，将其完整芯样进行成像；对于路面结构厚度≥400mm的芯样，因难以一次性取出厚度方向的所有试样，其取出的芯样多呈松散状，即面层、基层、底基层等结构层各呈一个芯样或部分结构层粘结、部分结构层脱离，定义为分层芯样，将每一个分层芯样的顶端画若干个方向线，各个分层芯样分别单独成像后组合，成为一个完整的芯样剖面图；扫描结束后，导出扫描图像，所获取的图像全部为灰度图，且获取的是沿扫描方向的多个图像。

优选的，每一个分层芯样的方向线数量及方向相同，方便后续同一个剖面不同分层的编号，能够反映同一个剖面的原材料分布变化。

步骤2.3，对于超过两层结构层粘结的芯样，按照结构层将芯样的扫描图像进行分割，按照结构层厚度信息将计算各结构层的分界像素点，依据像素点裁剪图像，最终获得面层、基层等各单结构层方向线所在平面的剖面扫描图。

步骤2.4，将各单结构层剖面扫描图导入Mimics软件，根据常用的阈值进行参数设定，进行阈值分割。不同的结构层，在分割时，阈值设定的不同；如对于路面结构第一层(面层)，设定三种阈值分别提取，分割完成后分别获得粘结剂灰度图、空隙灰度图、集料灰度图三种组分图像，并将其转化为二值图像；对于路面结构的第二层，会设定另外三种阈值，最终形成第二层的三种组分图像，其具体步骤如下：

对于第一层第一个方向得到的剖面扫描图，经过阈值分割及二值化处理最终得到的图像包括三张，分别为粘结剂二值图(A1-1)、空隙二值图(B1-1)、集料二值图(C1-1)，其中组分灰度为255(白)，图像其余部分均为0(黑)。对于其余各层(基层、底基层等)采用同样的方法得到其五个方向的组分二值图，并按照结构层出现顺序依次标号，如A2-3为第二层粘结剂的第三个方向获得的二值图，将每一个组分的，各个方向各个层的该组分二值图组合，最终得到的为粘结剂二值图集合A＝{A1-1，A1-2……，Am-n，……}，空隙二值图集合B＝{B1-1，B1-2……，Bm-n，……}，集料二值图集合C＝{C1-1，C1-2……，Cm-n，……}，其中m表示自上而下的结构层顺序，n表示第几次扫描获得(m≤结构层层数，n≤5)。

具体的粘结剂二值图能够反映出粘结剂的形状及分布状况，比如小块状及在图中的位置；空隙灰度图反映的是空隙的形状及位置分布，如不规则块条状及在图中的位置；集料灰度图反映的是集料的形状及位置分布，如圆形、不规则多边形等形状及在图中的位置。

上述过程中，对于分层芯样，将分层芯样分别进行分割并分别得到三种组分图(多个扫描方向)，并进行对应的编号。

需要说明的是，扫描过程中所使用的无损X-Ray ICT设备允许扫描的工件尺寸一般为最大外直径170mm，选用扩展模式可达到260mm，允许工件最大高度为500mm，以保证各结构层均能被扫描到。

优选的，对试样在5个不同的方向分别进行成像，以获取5种差异化的数字化图像，保证后期建模的真实度。

步骤3，设置空间步长，在MATLAB中构建路面结构模型，获得路面结构结构层矩阵；对每一个组分的二值图集合进行尺寸变化，使其单像素代表尺寸对应空间步长，并导入MATLAB中获得组分矩阵；分别使用各组分各结构层的不同角度扫描获得的组分矩阵对相应的路面结构结构层进行填充，获得最终的路面结构模型。

步骤3.1，设置空间步长，在MATLAB中构建初始路面结构模型，其中模型分为空气层与路面结构层，空气层在路面结构层的上方。考虑到后续模型周边会设置完美匹配层(PML)，因此设置空气层以消除PML层对电磁波的吸收效应，并作为探地雷达收发天线的放置层。通过MATLAB中rand()函数构建随机数矩阵以表示空气层与路面结构层，空气层为单个矩阵，路面结构层包含i个矩阵，即S＝{M1，M2，……，Mi}，i为路面结构层层数。其中路面结构层矩阵的垂直网格数G_v应根据空间垂直方向步长dy与实际取芯高度h计算，如式(1)所示，其水平网格数G_h根据道路实际情况及计算机配置设置。

其中，Gy代表通过实际长度/空间步长得出的垂直方向上应该布置的网格数量，y代表垂直方向。

步骤3.2，分别对组分二值图集合A、B、C进行图像尺寸变换，使其单个像素点代表实际尺寸等于空间步长，新的图像尺寸由式(2)计算得到，通过自编写python程序将组分二值图A、B、C尺寸变换为新尺寸A_trans、B_trans、C_trans，并导入MATLAB中获得组分矩阵M-A_trans＝{M-A1-1，M-A1-2……，M-Am-n，……}，M-B_trans＝{M-B1-1，M-B1-2……，M-Bm-n，……}，M-C_trans＝{M-C1-1，M-C1-2……，M-Cm-n，……}，其中m表示自上而下的结构层顺序，n表示第几次扫描获得(m≤结构层层数，n≤5)。

其中，l_ori代表二值图对应的真实长度，h_ori代表二值图对应的真实宽度，dx代表长度方向的空间步长，dy代表宽度方向的空间步长。

通过二值图组分矩阵的转换，使得最后获得的道路模型更加精准。其中l_ori代表二值图对应的真实长度，h_ori代表二值图对应的真实宽度，dx代表长度方向的空间步长，dy代表宽度方向的空间步长。

步骤3.3，对M-A_trans、M-B_trans、M-C_transM中的每一个组分根据材料常见介电数值进行矩阵赋值，路面材料的常见介电常数如表1。

表1路面材料常见介电常数值

步骤3.4，使用M-A_trans、M-B_trans、M-C_trans中的各结构层各角度矩阵对S中对应路面结构结构层矩阵进行随机赋值填充，获得最终的路面结构模型。

优选的，为了方便路面及病害模型尺寸参数的设定，设定空间步长{dx,dy,dz}＝{0.002m,0.002m,0.002m}；优选的，考虑到病害模型的水平尺寸及模型参数量，路面结构层矩阵水平网格数设为1000，对应实际长度2m，空气层矩阵垂直网格数设为50，对应实际高度0.1m，水平网格数与路面结构层矩阵保持一致。

步骤4，对步骤3中构建成功的路面结构模型嵌入不规则隐蔽病害模型，包括不规则层间粘结不良、不规则脱空、裂缝等路面结构病害，构建病害嵌入的路面结构模型，具体包括以下内容：

(1)，嵌入不规则层间粘结不良或/与脱空病害。

a，以面层与基层不规则层间粘结不良建模为例(图3)，首先选定病害起始点P_s、病害垂直高度H_c、病害水平延伸长度L_c(均为网格数量)；

b，在MATLAB中使用rand()函数构建尺寸为L_c*1的列向量list，并对其中值选定范围分别进行-2、-1、0、1、2赋值，表示后续填充的移动方向，其中0代表不移动，正值代表向正方向移动，负值代表向负方向移动；

c，构建循环函数，对每一个点分别以此点为中心分别向上、向下填充H_c/2个单位。其中第一个点为P_s，第二个点在P_s的基础上沿正方向水平移动一个单位并根据list中相应索引对应的值进行垂直移动，在移动后继续以此点为中心分别向上、向下填充H_c/2个单位，随后继续重复上述操作，当遍历完L_c个点后完成不规则层间粘结不良建模(图3)。不规则脱空病害构建方法与层间粘结不良相同，只需根据实际情况调整P_s、H_c、L_c即可。

(2)，嵌入裂缝病害。

a，构建垂直裂缝。选定裂缝发生起始位置P_l、裂缝宽度L_l、裂缝垂直长度H_l，以P_l为左上角顶点构建长H_l、宽L_l的矩形。

b，构建不规则裂缝。①选定裂缝顶点位置P₁、裂缝宽度L_l、裂缝垂直长度H_l(均为网格数量)；②使用rand()函数构建尺寸为H₁*1的列向量，并对其中值选定范围分别进行-2、-1、0、1、2赋值，表示后续填充的移动方向，其中0代表不移动，正值代表向右移动，负值代表向左移动；③构建循环函数，对每一个点分别以此点为中心分别向左、向右填充L_l/2个单位。其中第一个点为P_l，第二个点在P_l的基础上垂直向下移动一个单位并根据list中相应索引对应的值进行移动，在移动后继续以此点为中心分别向左、向右填充L_l/2个单位，当遍历完H₁个点后完成不规则裂缝建模。

c，根据表1中的常用数值对病害模型进行介电常数赋值。

d，将构建完成的路面结构模型与病害模型导出为.h5格式的文件进行保存。

应当理解的是，上述内容给出了每一个病害嵌入的具体步骤，但是不同病害之间嵌入的顺序并没有特殊要求。

步骤5，撰写探地雷达正演模拟基础文件；调用计算机系统的控制台并输入程序运行指令，运行gprMax正演程序，完成探地雷达对所构建道路仿真模型的扫描；使用Python程序导入正演数据文件，获得正演反射振幅数组，分别绘制B-scan图像及A-scan堆积图像

步骤5.1，撰写探地雷达正演模拟基础文件，新建文本文档，根据实际待检测道路信息、所用探地雷达设备实际参数及前述步骤构建的路面结构及病害模型文件依次填写标题、建模域、网格大小、时窗、子波类型、源位置、接收位置、模型文件、显示模型、波场快照等信息。

步骤5.2，调用计算机系统的控制台并输入程序运行指令，运行gprMax 3D正演程序，完成探地雷达对所构建道路仿真模型的扫描，最终获得out后缀的数据文件，该文件合并了各单道波正演数据，即A-scan数据，记录着正演各单道波的反射振幅数组。

步骤5.3，获取合成探地雷达A-scan数据及B-scan数据；

步骤5.31，为了进一步分析病害特征的振幅、相位变化，使用Python运行单道波堆积图生成程序，根据正演设置的采样间隔及采集道数绘制相应数量的A-scan图像(图4)，其中单道波的正相用黑色填充，负相用白色填充。

步骤5.32，为了直观看到病害在目前最常用、最具代表性的雷达图像上的典型特征，使用python程序导入步骤5.2获得的out单道波合并文件，获得反射振幅数组，并以灰度伪彩图进行填充，获得合成B-scan图像(图5)，该探地雷达图像包含目标体病害反射特征及不均匀介质反射杂波特征。

步骤5中，优选的，模型的z方向长度设置为单个空间步长dz(0.002m)，以减少建模参数量，提高建模速度。

实施例1

实施例1描述了一种基于CT扫描与数值模拟的仿真GPR数据获取方法优选具体实施方式，具体实施情况如下：

步骤1，使用混凝土钻孔取芯机在选定点处进行钻芯取样，选择直径为100mm、长度为400mm的钻孔取芯机钻头，并进行道路路面结构结构层信息计算及统计；

步骤1.1，针对待检测道路，根据道路建设资料，选定具有合适孔径及钻芯深度的钻头并选择几处代表性路表无病害处进行取芯，以获取完整的路面结构组合信息。取出的芯样尽量保证各结构层粘结、整体性好，便于后续的CT扫描。

步骤1.2，对于取出的芯样使用直尺对其面层、基层、底基层、垫层等结构层进行测量，分别沿不同边线共进行三次测量并取平均值。及时记录各结构层材料、厚度等信息，便于后续送入X-Ray ICT平台进行数字化成像。本实施例以面层平均厚度10cm，基层平均厚度20cm的整体芯样为例进行说明。

步骤2，对钻心试件使用X-Ray ICT设备进行扫描并获取剖面扫描图，对其进行裁剪获取单结构层剖面扫描图；将单结构层剖面扫描图导入Mimics软件进行阈值分割(图2)，获得单组分二值图；

步骤2.1，对各个芯样进行标记。优选的，使用红色记号笔以沿72°为一个间隔在芯样顶端进行标记，共作5个标记，代表沿5个方向分别进行成像。

步骤2.2，将试件表面清理干净放至X-Ray ICT试验平台进行扫描，设置扫描参数，开启试验平台对试件进行扫描。由于取出芯样面层、基层粘结好，呈一个整体芯样，且整体厚度为300mm小于400mm，因此将其完整芯样进行扫描；扫描结束后，导出扫描图像，所获取的图像全部为灰度图，且获取的是沿5个扫描方向的图像。

步骤2.3，对于芯样整体扫描图像进行分割，按照结构层厚度信息计算面层、基层的分界像素点，依据像素点裁剪图像，最终获得面层、基层的单结构层剖面扫描图。

步骤2.4，将各单结构层剖面扫描图导入Mimics软件，根据常用的阈值进行参数设定，进行阈值分割。以沥青混合料对于面层，设定三种阈值分别提取，根据Mimics灰度直方图处理及分析，本实施例中设定空隙、沥青砂浆、粗骨料的灰度值范围分别为0-60、60～120、120-255，分割完成后分别获得沥青粘结剂灰度图、空隙灰度图、集料灰度图三种组分图像，并将其转化为二值图像；对于基层，同样设定三种阈值分别提取，在本实施例中，根据Mimics灰度直方图处理及分析，确定空隙的灰度范围为50～120，水泥粘结剂的灰度范围为175～200，集料的灰度范围为210～255，分割完成后分别获得水泥粘结剂灰度图、空隙灰度图、集料灰度图三种组分图像，并将其转化为二值图像。优选的，步骤S2.4的具体步骤如下：

对于面层第一个方向得到的剖面扫描图，经过阈值分割及二值化处理最终得到的图像包括三张，分别为粘结剂二值图(A1-1)、空隙二值图(B1-1)、集料二值图(C1-1)，其中组分灰度为255(白)，图像其余部分均为0(黑)。对于基层采用同样的方法得到其五个方向的组分二值图，并按照结构层出现顺序依次标号，如A2-3为基层粘结剂的第三个方向获得的二值图。最终得到的为粘结剂二值图集合A＝{A1-1，A1-2……，A2-5}，空隙二值图集合B＝{B1-1，B1-2……，B2-5}，集料二值图集合C＝{C1-1，C1-2……，C2-5}，对于二值图的命名规则，Am-n，其中m表示自上而下的结构层顺序，n表示第几次扫描获得(m≤2，n≤5)，此处由于是两层结构层及5次扫描，因此m的最大值为2，n的最大值为5。

步骤3，设置空间步长，在MATLAB中构建路面结构模型，获得路面结构结构层矩阵；对组分二值图进行尺寸变化，使其单像素代表尺寸对应空间步长，并导入MATLAB中获得组分矩阵；分别使用各组分各结构层的不同角度扫描获得的组分矩阵对相应的路面结构结构层进行填充，获得最终的路面结构模型。

步骤3.1，设置空间步长，在MATLAB中构建路面结构模型，其中模型分为空气层与路面结构层，优选的，步骤3.1具体如下：

设定空间步长{dx,dy,dz}＝{0.002m,0.002m,0.002m}。在MATLAB中构建路面结构模型，其中模型分为空气层与路面结构层。考虑到后续模型周边会设置完美匹配层(PML)，因此设置空气层以消除PML层对电磁波的吸收效应，并作为探地雷达收发天线的放置层。通过MATLAB中rand()函数构建随机数矩阵以表示空气层与路面结构层，空气层包含单个矩阵，路面结构层包含2个矩阵，即S＝{M1，M2}。其中路面结构层矩阵的垂直网格数G_v应根据空间垂直方向步长dy与实际取芯高度h计算，如式(1)所示，此实施例中面层的垂直网格数为0.1/0.002＝50个，基层的垂直网格数为0.2/0.002＝100个，因此G_y等于150。其水平网格数G_h根据道路实际情况及计算机配置设置，优选的，考虑到病害模型的水平尺寸及模型参数量，路面结构层矩阵水平网格数G_h设为1000，对应实际长度2m，空气层矩阵垂直网格数设为50，对应实际高度0.1m，水平网格数与路面结构层矩阵保持一致。

步骤3.2，分别对组分二值图A、B、C进行图像尺寸变换，使其单个像素点代表实际尺寸等于空间步长，通过MATLAB获得各组分矩阵集合，优选的，步骤3.2具体如下：

新的图像尺寸由式(2)计算得到。对于面层，其图像新尺寸为(0.1/0.002)*(0.1/0.002)＝50*50；对于基层其图像新尺寸为(0.1/0.002)*(0.2/0.002)＝50*100。通过自编写python程序将组分二值图A、B、C尺寸变换为新尺寸A_trans、B_trans、C_trans，并导入MATLAB中获得组分矩阵M-A_trans＝{M-A1-1，M-A1-2……，M-A2-5}，M-B_trans＝{M-B1-1，M-B1-2……，M-B2-5}，＝{M-C1-1，M-C1-2……，M-C2-5}，其中各组分矩阵集合M-A_trans、M-B_trans、M-C_trans均各自包含2个结构层、5个方向共10个矩阵。通过二值图组分矩阵的转换，使得最后获得的道路模型更加精准。其中l_ori代表二值图对应的真实长度，h_ori代表二值图对应的真实宽度，dx代表长度方向的空间步长，dy代表宽度方向的空间步长。

步骤3.3，对M-A_trans、M-B_trans、M-C_transM根据材料常见介电数值进行矩阵赋值，路面材料的常见介电常数如表1。在本实施例中，面层与基层粘结剂成分不同，其他空气、集料两种组分性质相同，因此对M-A1-1、M-A1-2、M-A1-3、M-A1-4、M-A1-5均赋值沥青(3～5)，对于M-A2-1、M-A2-2、M-A2-3、M-A2-4、M-A2-5均赋值水泥(4-6)，对于M-B1-1、M-B1-2、M-B1-3、M-B1-4、M-B1-5、M-B2-1、M-B2-2、M-B2-3、M-B2-4、M-B2-5均赋值1，对于M-C1-1、M-C1-2、M-C1-3、M-C1-4、M-C1-5、M-C2-1、M-C2-2、M-C2-3、M-C2-4、M-C2-5均赋值石灰岩数值7。

表1路面材料常见介电常数值

步骤3.4，使用M-A_trans、M-B_trans、M-C_trans中的各结构层各角度矩阵对S中对应路面结构结构层矩阵进行随机赋值填充，获得最终的路面结构模型。

步骤4，对S3中构建成功的路面结构模型嵌入不规则隐蔽病害模型，包括不规则层间粘结不良、不规则脱空、裂缝等路面结构病害，构建病害体嵌入的路面结构模型，此处以层间粘结不良为例。

(1)，嵌入不规则层间粘结不良病害。

a，以面层与基层不规则层间粘结不良建模为例，选定病害起始位置P_s(x：0，y：100)、病害垂直高度H_c(8)、病害水平延伸长度L_c(1000)(均为网格数量)；

b，在MATLAB中使用rand()函数构建尺寸为1000*1的列向量list，并对其中值选定范围分别进行-2、-1、0、1、2赋值，表示后续填充的移动方向，其中0代表不移动，正值代表向正方向移动，负值代表向负方向移动；

c，构建循环函数，对每一个点分别以此点为中心分别向上、向下填充4(H_c/2)个单位。其中第一个点为P_s，第二个点在P_s的基础上沿正方向水平移动一个单位并根据list中相应索引对应的值进行垂直移动，在移动后继续以此点为中心分别向上、向下填充4(H_c/2)个单位，随后继续重复上述操作，当遍历完1000(L_c)个点后完成不规则层间粘接不良建模(图3)。

(2)，根据表1中的常用数值对病害模型进行介电常数赋值。

(3)，将构建完成的路面结构模型与病害模型导出为.h5格式的文件进行保存。

步骤5，撰写探地雷达正演模拟基础文件；调用计算机系统的控制台并输入程序运行指令，运行gprMax正演程序，完成探地雷达对所构建道路仿真模型的扫描；使用Python程序导入正演数据文件，获得正演反射振幅数组，分别绘制B-scan图像及A-scan堆积图像

步骤5.1，撰写探地雷达正演模拟基础文件，新建文本文档，根据实际待检测道路信息、所用探地雷达设备实际参数及前述步骤构建的路面结构及病害模型文件依次填写标题、建模域、网格大小、时窗、子波类型、源位置、接收位置、模型文件、显示模型、波场快照等信息。优选的，模型的z方向长度设置为单个空间步长dz，以减少建模参数量，提高建模速度。

步骤5.2，调用计算机系统的控制台并输入程序运行指令，运行gprMax 3D正演程序，完成探地雷达对所构建道路仿真模型的扫描，最终获得out后缀的数据文件，该文件合并了A-scan数据，记录着正演各单道波的反射振幅数组。

步骤5.3，获取合成探地雷达A-scan及B-scan数据，优选的，步骤5.3具体如下：

步骤5.31，为了进一步分析病害特征的振幅、相位变化，使用Python运行单道波堆积图生成程序，根据正演设置的采样间隔及采集道数绘制相应数量的A-scan图像(图4)，其中单道波的正相用黑色填充，负相用白色填充。

步骤5.32，为了直观看到病害在目前最常用、最具代表性的雷达图像上的典型特征，使用python程序导入步骤5.2获得的out单道波合并文件，获得反射振幅数组，并以灰度伪彩图进行填充，获得合成B-scan探地雷达图像(图5)，该探地雷达图像包含目标体病害反射特征及不均匀介质反射杂波特征。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种基于CT扫描与数值模拟的仿真GPR数据获取方法，其特征在于，包括以下步骤：

步骤1，在待检测路段上钻取芯样，所述芯样由N个结构层组成，获得每一个结构层的厚度，N为≥1的自然数；

步骤2，扫描芯样的若干个纵向剖面，获得若干个剖面扫描图，所述芯样的轴线在纵向剖面上；按照结构层分割剖面扫描图，获得各个结构层的若干幅纵向剖面扫描图，阈值分割每一个结构层的所有纵向剖面扫描图，获得每一个结构层的剖面扫描图的组分图像，所述组分图像包括粘结剂灰度图、空隙灰度图和集料灰度图；将所有的组分图像转化为二值图像，组合后获得所述芯样每一个组分的二值图集合；

步骤3，在MATLAB中构建初始路面结构模型，设置空间步长，对所有的二值图集合进行尺寸变换，获得变换后的组分矩阵，变换后组分矩阵中像素点代表的实际尺寸与路面结构模型中的空间步长相等；根据组分矩阵中各个元素的材料类型赋值介电常数，获得路面结构模型；

步骤4，在路面结构模型中嵌入路面结构病害，获得病害嵌入的路面结构模型；

步骤5，运行gprMax正演程序，通过探地雷达对病害体嵌入的路面结构模型扫描，获得单道波正演数据，进一步的合成探地雷达A-scan图像及B-scan图像。

2.根据权利要求1所述的一种基于CT扫描与数值模拟的仿真GPR数据获取方法，其特征在于，步骤1中，钻取芯样的钻头直径大于100mm，钻头长度为400mm。

3.根据权利要求1所述的一种基于CT扫描与数值模拟的仿真GPR数据获取方法，其特征在于，步骤1中，若取出的芯样结构层分离，步骤2中针对分离的每一个分层芯样均获得若干个剖面扫描图；

将同一个芯样的分层芯样同一个剖面的剖面扫描图组合后，获得所述芯样完整的剖面扫描图。

4.根据权利要求1所述的一种基于CT扫描与数值模拟的仿真GPR数据获取方法，其特征在于，步骤2中，按照结构层从下到上顺序，组合获得每一个组分的二值图集合。

5.根据权利要求1所述的一种基于CT扫描与数值模拟的仿真GPR数据获取方法，其特征在于，步骤3中，所述初始路面结构模型包括空气层和路面结构层；

所述空气层为单个矩阵，所述路面结构层包括i个矩阵，i为路面结构层的层数；所述i个矩阵的垂直网格数由空间垂直方向步长和实际取芯高度获得。

6.根据权利要求1所述的一种基于CT扫描与数值模拟的仿真GPR数据获取方法，其特征在于，步骤3中，二值图集合进行尺寸变换的公式为：

其中，l_ori代表二值图对应的真实长度，h_ori代表二值图对应的真实宽度，dx代表长度方向的空间步长，dy代表宽度方向的空间步长。

7.根据权利要求1所述的一种基于CT扫描与数值模拟的仿真GPR数据获取方法，其特征在于，步骤4中，所述路面结构病害包括不规则层间粘结不良病害、不规则脱空病害和裂缝病害。

8.根据权利要求7所述的一种基于CT扫描与数值模拟的仿真GPR数据获取方法，其特征在于，所述不规则层间粘结不良病害或所述不规则脱空病害的构建过程为：

(1)确定病害的水平起始位置、病害垂直高度和病害水平延伸长度；

(2)在MATLAB中构建尺寸为Lc*1的列向量及循环函数，循环函数根据列向量中索引的值移动，重复遍历Lc个点后获得不规则层间粘结不良或不规则脱空病害。

9.根据权利要求7所述的一种基于CT扫描与数值模拟的仿真GPR数据获取方法，其特征在于，步骤4中，所述裂缝病害的构建过程为：

(1)构建垂直裂缝，确定裂缝发生的起始位置、裂缝宽度和裂缝垂直长度；

(2)构建尺寸为H1*1的列向量和循环函数，循环函数根据列向量中索引的值移动，重复遍历H1个点后获得裂缝病害。

10.根据权利要求1-9任意一项所述的基于CT扫描与数值模拟的仿真GPR数据获取方法，其特征在于，步骤5中，运行单道波堆积图生成程序，根据单道波正演数据获得A-scan图像；基于单道波正演数据，获得反射振幅数组，填充灰度伪彩图，获得B-scan图像。