技术领域
[0001] 本
发明涉及地球物理探测技术领域，特别涉及一种基于PyTorch的探地雷达层析成像加速方法。
背景技术
[0002] 灾害资源事故是
煤矿安全面临的重大地质问题之一。地质安全保险制度建设不容忽视。采煤前查明采煤工作面或矿区的地质条件，是保证高产、高效益矿井安全生产的重要前提。
[0003] 探地雷达(GPR)是一种高
分辨率、高效率的非破坏性地球物理勘探技术。在地面雷达测量中，共偏移反射类型常用于地下结构的确定，而在煤矿安全探测中，传输式探地雷达的探测距离较大，其发射天线和接收天线分别位于探测剖面两侧。
[0004] 层析成像分为两种:速度层析成像和衰减层析成像。速度层析成像是根据初至的旅行时间来反演被测剖面的速度，而衰减层析成像则使用振幅或
频率数据。一些研究人员还开发了旅行时间和振幅的联合反演方法。雷达探测到的介质主要是电导率为零的有耗介质，
介电常数与衰减无关。目前，成像的速度还有待提高和优化。
发明内容
[0005] 本发明提出了一种基于PyTorch的探地雷达层析成像加速方法，该方法为利用PyTorch加速GPR速度层析成像反演的并行方法，可加速成像求解。
[0006] 本发明
实施例提供一种基于PyTorch的探地雷达层析成像加速方法，包括以下步骤：
[0007] S1、利用GPR进行层析探测时，探测区域两侧沿巷道布置两条平行线；发射天线端位于其中一行发射高频电磁脉冲，接收天线端位于另一行记录在探测介质内传播后的脉冲；
[0008] S2、使用GPR旅行时间层析成像进行地雷探测时，记录接收天线端的
电压作为时间的函数，称为轨迹；利用初至时间进行层析反演，得到发射天线端和接收天线端的速度；
[0009] S3、采用神经网络方法求解函数，采用PyTorch
框架调用CUDA加速GPU的运算；
[0010] S4、使用GPU来计算数据流图；
[0011] S5、使用Adam来实现更新规则，计算梯度的一阶矩估计和二阶矩估计，设计不同参数的独立自适应学习率。
[0012] 在一个实施例中，所述步骤S2还包括：
[0013] 基于射线的速度反演需要获得旅行时间，应用式(1)建立旅行时间与慢度的关系，即速度的反演；
[0014] t＝∫rays(l)dl
[0015] (1)
[0016] 其中t表示旅行时间，s表示代表慢度场，l表示射线路径长；
[0017] 对式(1)进行网格离散化和射线追踪，得到离散化的方程(2)；
[0018] T＝RS
[0019] (2)
[0020] 其中T表示旅行时间矩阵，R表示根据网格离散化结果的射线段长度矩阵，S表示慢度值矩阵。
[0021] 在一个实施例中，所述步骤S2还包括：
[0022] 构造一个基于正则化表达的最小二乘
算法的目标函数(3)；所述目标函数为损失函数，基于合成旅行时与观测旅行时的一致性来测量模型参数的
质量；
[0023] J(S)＝||RΔS-(Tobs-Tsyn)||2+λ||DS||2
[0024] (3)
[0025] λ表示拉格朗日乘数，D表示二阶空间导数矩阵，Tobs表示观测旅行时，Tsyn表示合成旅行时，△S表示慢度更新值；
[0026] 将模型参数m的最佳估计定义为最小化方程(4)中的损失函数的模型；对反演结果进行成像，以识别地质异常体；
[0027] S＝argmin(J(S))
[0028] (4)
[0029] 方程(4)应用PyTorch的神经网络框架求解。
[0030] 在一个实施例中，所述方法还包括：
[0031] S6、设计两个数值模型，用于测试和分析所提出的层析成像反演算法；根据到介质分布的不均匀性，将每种介质的相对介电常数和
波速设置为预定范围内的随机值；
[0032] S7、应用两次射线追踪程序进行反演，在每个
光线跟踪过程中，分别对基于神经网络NN和基于最小二乘QR分解法的方法应用预设次数
迭代，获得反演结果；
[0033] S8、根据公式(5)计算的时间差评估反演结果，基于神经网络NN的层析成像计算速度优于最小二乘QR分解法；
[0034]
[0035] 其中，N表示迹线数，DT表示时间差。
[0036] 本发明实施例提供的一种基于PyTorch的探地雷达层析成像加速方法，该方法采用PyTorch提供的神经网络框架作为并行方案，有效地利用了GPU的优势，利用PyTorch的自动求导机制和Adam算法对参数进行更新。与传统的最小二乘QR分解法(LSQR)反演方法相比，基于PyTorch神经网络框架的反演方法具有计算速度快、性能稳定等优点。利用PyTorch神经网络框架提供的自动求导和更新算法的特点，用户更容易以高效的并行方式实现自己的算法，可加速成像求解。
[0037] 本发明的其它特征和优点将在随后的
说明书中阐述，并且，部分地从说明书中变得显而易见，或者通过实施本发明而了解。本发明的目的和其他优点可通过在所写的说明书、
权利要求书、以及
附图中所特别指出的结构来实现和获得。
[0038] 下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。
附图说明
[0039] 附图用来提供对本发明的进一步理解，并且构成说明书的一部分，与本发明的实施例一起用于解释本发明，并不构成对本发明的限制。在附图中：
[0040] 图1为本发明实施例提供的雷达层析探测现场的工作示意图；
[0041] 图2为本发明实施例提供的煤矿探地雷达层析成像方法简要
流程图；
[0042] 图3为本发明实施例提供的层析成像反演的两层神经网络结构示意图；
[0043] 图4为本发明实施例提供的公式(2)的数据流图描述；
[0044] 图5a为本发明实施例提供的以干煤为背景的矩形气孔的模型；
[0045] 图5b为本发明实施例提供的以干煤为背景的两个矩形气孔的模型；
[0046] 图6a为本发明实施例提供的具有一个矩形气隙的基于NN的方法的反演结果示意图；
[0047] 图6b为本发明实施例提供的具有一个矩形气隙的基于LSQR的方法的反演结果示意图；
[0048] 图7a为本发明实施例提供的具有两个气隙的基于NN的方法的反演结果示意图。
[0049] 图7b为本发明实施例提供的具有两个气隙的基于LSQR的方法的反演结果示意图。
具体实施方式
[0050] 下面将参照附图更详细地描述本公开的示例性实施例。虽然附图中显示了本公开的示例性实施例，然而应当理解，可以以各种形式实现本公开而不应被这里阐述的实施例所限制。相反，提供这些实施例是为了能够更透彻地理解本公开，并且能够将本公开的范围完整的传达给本领域的技术人员。
[0051] 首次到达时间的速度层析反演通常采用最小二乘QR分解法(LSQR)方法。通过引入张量流，在
地震反演区引入了一种新的求解方法。在本发明中，提出了一种利用PyTorch加速GPR速度层析成像反演的并行方法。并与传统的LSQR方法进行了比较，验证了该方法的有效性。
[0052] 其中涉及的PyTorch是一个由Facebook开发的开源
深度学习框架，它封装了自动求导机制和更新算法，可以方便地应用。本发明实施例提供的基于PyTorch的探地雷达层析成像加速方法，通过引入张量流，在地震反演区引入了一种新的求解方法，包括：
[0053] 步骤一：如图1所示，利用GPR进行层析探测时，探测区域两侧沿巷道布置两条平行线。发射天线端位于其中一行发射高频电磁脉冲，接收天线端位于另一行记录在探测介质内传播后的脉冲。
[0054] 步骤二：如图2所示，使用GPR旅行时间层析成像技术进行地雷探测时，记录接收端的电压作为时间的函数，称为轨迹。然后利用初至时间进行层析反演，考虑到波沿射线路径传播，得到发射天线端和接收天线端的速度。
[0055] 基于射线的速度反演需要获得旅行时间，应用式(1)建立旅行时间与慢度的关系，即速度的反演。
[0056] t＝∫rays(l)dl
[0057] (1)
[0058] 其中t表示旅行时间，s表示代表慢度场，l表示射线路径长。
[0059] 然后对式(1)进行网格离散化和射线追踪，得到离散化的方程(2)；
[0060] T＝RS
[0061] (2)
[0062] 其中T表示旅行时间矩阵，R表示根据网格离散化结果的射线段长度矩阵，而S表示慢度值矩阵。
[0063] 为了稳定反演，需要构造一个基于正则化技术的最小二乘算法的目标方程(3)。目标函数也称为损失函数，它基于合成旅行时与观测旅行时的一致性来测量模型参数的质量。
[0064] J(S)＝||RΔS-(Tobs-Tsyn)||2+λ||DS||2
[0065] (3)
[0066] λ表示拉格朗日乘数，D表示二阶空间导数矩阵，Tobs表示观测旅行时，Tsyn表示合成旅行时，△S表示慢度更新值。
[0067] 将模型参数m的最佳估计定义为最小化方程(4)中的损失函数的模型。最后，对反演结果进行成像，以识别地质异常体。
[0068] S＝argmin(J(S))
[0069] (4)
[0070] 方程(4)应用PyTorch的神经网络框架求解。
[0071] 步骤三：采用神经网络方法求解函数，采用PyTorch框架调用CUDA加速GPU的运算。比如图3所示的两层神经网络可用于求解目标函数。
[0072] 使用PyTorch求解神经网络具有如下优点:用于计算的PyTorch中的张量类似于一个numpy数组，但可以在GPU上运行。利用张量构建计算图和自动求解梯度的自动封装。模型A神经网络层可以存储状态或可学习权值。
[0073] 可以用图4来描述PyTorch中的方程(2)。需要注意的是，在神经网络中通常存在一个偏置参数，但在GPR层析成像中不需要偏置参数，因为收发天线之间的时延和同步已经进行了调整。
[0074] 步骤四：使用GPU来计算数据流图，只需要使用一个不同的设备:(1)device＝torch.device('cuda:0')(2)x＝torch.randn(M,D,device＝device)。使用requires_grad＝True创建张量可以启用自动求导，使用requires_grad＝True创建张量的操作可以生成计算图。
[0075] 步骤五：使用Adam来实现更新规则，Adam是一种可以替代传统随机梯度下降过程的一阶
优化算法；它可以根据
训练数据对神经网络的权值进行迭代更新，Adam算法不同于传统的随机梯度下降算法；随机梯度下降保持一个单一的学习速率(alpha)来更新所有的权值，并且学习速率在训练过程中不会改变；Adam通过计算梯度的一阶矩估计和二阶矩估计，设计了不同参数的独立自适应学习率。
[0076] 步骤六：设计了两个数值模型，用于测试和分析所提出的层析成像反演算法。图5a显示了一个以干煤为背景的矩形气孔的模型，图5b显示了以干煤为背景的两个矩形气孔的另一个模型。考虑到介质分布的不均匀性，将每种介质的相对介电常数和波速设置为给定范围内的随机值。
[0077] 步骤七：应用两次射线追踪程序进行反演。在每个
光线跟踪过程中，分别对基于神经网络NN和基于LSQR的方法应用比如1000次迭代，然后获得反演结果。具有一个矩形气隙的第一个模型的反演结果如图6a～6b所示。图6a是基于NN的方法得出的，而图6b是基于LSQR的方法得出的。具有两个气隙的第二个模型的反演结果如图7a～7b所示。图7a是基于NN的方法得出的，图7b是基于LSQR的方法得出的。
[0078] 步骤八：为了评估这两种算法的性能，给出了反演结果与表1中的输入之间由公式(5)计算的时间差。应注意，表中的NN-1表示模型1的基于NN的方法，表1中的LSQR-1表示模型1的基于LSQR的方法。
[0079]
[0080] 其中N表示迹线数，DT表示时间差。
[0081] 表一：
[0082] Content NN-1 LSQR-1 NN-2 LSQR-2DT 4.70 5.24 4.33 16.83
[0083] 通过与传统的LSQR方法进行了比较，验证了该基于PyTorch的探地雷达层析成像加速方法的有效性。
[0084] 本实施例中，该方法采用PyTorch提供的神经网络框架作为并行方案，有效地利用了GPU的优势。利用PyTorch的自动求导机制和Adam算法对参数进行更新。数值实验表明，与传统的LSQR反演方法相比，基于PyTorch神经网络框架的反演方法具有计算速度快、性能稳定等优点。利用PyTorch神经网络框架提供的自动区分和更新算法的特点，用户更容易以高效的并行方式实现自己的算法。
[0085] 在GPR速度层析成像中，走时速度层析成像反演的第一步是剖面的离散化。剖面被离散为一系列单元。对于每个发射点-接收点对，计算通过单元的每个射线路径段的长度。所有的线段长度都组织在一个稀疏矩阵中，称为系数矩阵，它描述了射线的几何形状。因为每一条射线都穿过少数元素，所以系数矩阵中有很多零元素。这些零元素浪费了大量的存储内存和计算时间。本发明实施例设计了层析成像矩阵的稀疏存储和求解方法。在传统的稀疏矩阵结构中，每一行的非零元素个数和每一行的第一个非零元素个数相加。改进后的结构大大降低了层析成像反演算法的时间和空间复杂度，提高了反演效率。通过上述三次数值模拟试验表明，采用稀疏存储和求解方法，计算时间明显缩短。
[0086] 显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些
修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。
