CN113267287B

CN113267287B - 冲击波超压三维时空场重建方法

Info

Publication number: CN113267287B
Application number: CN202110725708.3A
Authority: CN
Inventors: 李剑; 孔慧华; 韩焱; 王黎明; 刘宾
Original assignee: North University of China
Current assignee: North University of China
Priority date: 2021-06-29
Filing date: 2021-06-29
Publication date: 2022-09-23
Anticipated expiration: 2041-06-29
Also published as: CN113267287A

Abstract

本发明涉及一种冲击波超压三维时空场重建方法，利用空中冲击波的速度信息以及冲击波射线到达时间信息相融合重建速度场，并通过速度和超压的线性关系间接计算超压场。与现有插值直接测量方法相比，能够在有限测试点的条件下，提高超压信号特征参数利用率，提高重建分辨率；采用压缩感知原理，将全变差TV最小化和字典学习相结合，实现对重建速度场的约束及优化，从而提高重建过程的稳定性，提升重建精度。

Description

冲击波超压三维时空场重建方法

技术领域

本发明属于爆破震动测试技术及重建领域，具体涉及一种冲击波超压三维时空场重建方法，尤其是一种空中爆炸冲击波超压威力场三维重建方法。

背景技术

近些年，随着工业、农业、及军事等相关安全领域的需求导向及学科的快速发展，爆炸冲击波动力学理论、超压测试理论、及冲击波场重建等方面的研究在国防及国民经济中起着至关重要的作用。目前的爆炸效能参数测试技术只能获得有限位置的参数值，不能获取全方位的爆炸信息，如何在有限点数据的基础上，重建出全方位的冲击波场信息，是近年来相关领域的研究者们研究的关键技术。

爆炸威力场分布是利用爆炸效能参数如冲击波超压、应力波能、动能等爆炸物理参量，以反演的方式重建其在爆炸后的强度分布，来表征爆炸后的威力场分布。目前主要有两种方法：(1) “点”测法，如生物实验法、等效罐法和等效靶板法，该方法把某种物品或者生物放到爆炸现场，在试验结束后，观察毁伤程度来评估毁伤威力，该方法虽然能够直观的看到地面上的毁伤效果，但其不能得到爆炸物理参数的具体数值及规律；(2)“场”测法，在不同的爆心距处布设传感器，将传感器阵列布设呈一定的形状(如一字阵、交叉阵等)，利用爆炸场物理参量(如冲击波超压、应力波能等)的衰减规律，通过拟合插值、数值模拟等方法重建爆炸威力场分布。与“点”测法相比，“场”测法通过反演爆炸参数的时空分布，能更全面了解爆炸的毁伤情况，准确评估爆炸的效能和威力，因此，该方法是实现对爆炸场毁伤威力有效评估的主要手段。

但目前“场”测法在针对空中爆炸后的冲击波威力场分布测试存在如下问题：

目前在威力场分布测试中，主要采用基于试验测试的离散点插值拟合法、基于理论分析的数值模拟法重建爆炸效能参数的分布。其中，离散点插值拟合法是利用有限测试节点的数据，以拟合插值方法得到物理参数的衰减变化规律，不能获取全方位爆炸信息；数值模拟法是利用有限元的方式对爆炸场进行数值模拟，虽然该方法能够对爆炸过程进行描述和分析，但在应用过程中，状态方程等参数的不同取值往往对模拟结果产生较大影响，导致研究结果与实际试验误差较大，时空场重建精度不高。

发明内容

本发明提供一种冲击波超压三维时空场重建方法，解决上述缺陷。

为了解决以上技术问题，本发明提供了一种击波超压三维时空场重建方法，其特征在于，具体包括以下步骤：

S1、将测试区域以炸点为中心设为正方形，并将其划分为多个区域，每个区域均以炸点为中心，并将每个区域划分为若干正方形网格单元；

S2、在划分好的每个区域边界处均布传感器；

S3、根据传感器获取的冲击波信号，提取冲击波超压时空场重建特征参数，包括每个传感器对应的超压峰值、到达时间；

S4、确定真实的冲击波超压峰值随距离的衰减规律模型；

S5、通过上述衰减规律模型计算出二维测试区域内所有网格处的超压峰值，根据超压峰值与冲击波速度以及慢度的关系，将超压场分布转换为慢度场分布，并将其作为走时层析成像反演的初始慢度值；

S6、建立走时层析成像模型；

S7、解算走时层析成像模型，得到二维慢度值；

S8、超压峰值场重建：将步骤S7中得到的二维慢度值转换为速度场，再转换为超压峰值Pm(x,y)，即得到冲击波超压峰值的空间分布；

S9、利用步骤S8中得到的超压峰值Pm(x,y)，得到具体的超压随时间和空间的关系式，重建冲击波超压时空场。

有益效果：

1、利用空中冲击波的速度信息以及冲击波射线到达时间信息相融合重建速度场，并通过速度和超压的线性关系间接计算超压场。与现有插值直接测量方法相比，能够在有限测试点的条件下，提高超压信号特征参数利用率，提高重建分辨率；

2、考虑到测试点数量有限，重建区域数据具有稀疏特性。采用压缩感知原理，将全变差TV(total variation)最小化和字典学习相结合，实现对重建速度场的约束及优化，从而提高重建过程的稳定性，提升重建精度。

附图说明

图1是网格划分示意图

图2是传感器布设示意图

图3是射线路径示意图。

具体实施方式

为使本发明的目的、内容和优点更加清楚，下面对本发明的具体实施方式作进一步详细描述。

本发明提出的一种冲击波超压三维时空场重建方法，其特征在于，具体包括以下步骤：

S1、划分网格

本发明采用均匀网格划分方法，

将测试区域以炸点为中心设为正方形，并将其划分为多个区域，每个区域均以炸点为中心，并将每个区域划分为若干正方形网格单元；具体为：

根据J.Henrych提出的超压峰值随距离衰减的经验公式中的比例距离

的区间，得到炸点到测点的距离r，根据该距离将 12m×12m的测试区域划分成对应的3个区域，并将每个区域划分为若干个1m×1m的正方形网格单元，网格划分越密，重建的精度就越高，具体划分如图1所示：

其中，Pm为冲击波超压峰值，

为比例距离，

r为炸点到测点的距离，W为装药重量。

S2、布设传感器

在爆炸走时层析成像中，当爆炸点位置确定时，一条射线对应着一个传感器的位置，因此，在网格划分确定的情况下，传感器的位置是影响投影矩阵结构的主要因素。传感器个数由实际情况决定，以四分之一为例，在划分好的3个区域边界处均布传感器，具体布设方案如图2所示；

S3、冲击波超压时空场重建特征参数提取

首先对传感器获取的冲击波信号进行预处理，采用小波变换对信号进行时频滤波，选择合适的阈值对小波系数进行阈值处理，实现信号去噪。

然后提取每个传感器对应的超压峰值、到达时间等效能参数。

S4、确定真实的冲击波超压峰值随距离的衰减规律模型

对于冲击波超压峰值Pm的计算，衰减规律模型都是以三项式来表达，公式如下：

其中，a、b、c为待求系数，Pm为冲击波超压峰值，

为比例距离。其中，根据炸点到测点的距离r以及

可以推倒出每个传感器对应的

结合试验数据与经验公式进行分析，将步骤3中每个传感器对应的冲击波超压峰值Pm和传感器的坐标信息代入式(2)，并通过最小二乘法进行拟合，得到参数a、b、c，建立更适合此次试验药量条件下的冲击波超压峰值随距离的衰减规律模型。

S5、计算走时层析成像反演的初始慢度值

通过上述衰减规律模型可计算出二维测试区域内所有网格处的超压峰值Pm(x,y)，已知冲击波超压峰值与冲击波速度的关系：

其中，Pm(x,y)为(x,y)位置处的超压峰值，P₀为环境大气压力，V(x,y) 为位置处(x,y)的冲击波速度，C₀为声速，

(T₀为未扰动的空气初始温度(K))。

冲击波速度V(x,y)与慢度S(x,y)的关系为：

由上式可将超压场分布转换为慢度场分布，并将其作为走时层析成像反演的初始慢度值。

S6、建立走时层析成像模型

走时层析成像的基础是Radon变化，假设爆炸波为连续速度场，爆炸波在传输的过程中，其走时是速度和几何路径的函数：

其中t为走时，即爆炸波传播到传感器节点的时间，v为速度，s为慢度，L为炸点到传感器的射线路径，dr为沿射线路径L的距离增量。

将(5)式离散化，以图3中的射线为例，即对第i条射线的投影数据有：

式中，t_i为第i条射线的走时，即爆炸波到达第i个传感器的时间； s_j为将二维网格区域中的慢度S(x,y)转换为1维列向量后，对应的第j个网格中的慢度；a_ij为第i条射线穿过第j个网格的射线长度；M为射线数，即传感器个数；N为划分的网格数，其中N是指将二维网格区域转换为1维列向量。s_j根据上述步骤S5中的 S(x,y)相对应。

因此可将反演问题这个非线性问题转化为线性问题，将(6) 式写成矩阵形式，得到走时层析成像模型：

AS＝T (7)

其中A为投影矩阵，S为慢度场重建图像，T为投影数据； T＝(t₁,t₂,…,t_M)'，即M维走时向量；S＝(S₁,S₂,…S_N)'，即N维慢度向量；A＝(a_ij)为M×N阶射线矩阵。

S7、解算走时层析成像模型

结合字典学习与TV(total variation)全变差，对走时层析成像模型构建目标函数：

其中，第一项为数据的保真项，第二项为TV(total variation)全变差正则项，第三项为字典学习正则项。μ为保真项系数，λ和β为正则项系数，E_i提取图像块的算子，v_j子字典学习正则化参数。上述系数在0-1之间，μ＝1，λ＝0.001，β＝0.01，v_j取常数0.002。

根据步骤5得到的初始慢度值，采用交替最小化的思想求解目标函数(8)，具体步骤如下：

步骤一：固定D和α_j，更新重建图像S，采用Split-Bregman 算法求解不可微正则化问题，首先引入辅助变量，将目标函数求解分为L₁范数和L₂范数两部分，然后采用Bregman算法求解。

目标函数变为：

为求解式(9)，引入中间变量

目标函数变为：

对于目标函数(10)，采用Bregman迭代求解，引入中间变量b_x，b_y，可得：

公式(11)的求解采用交替优化的思想可分为两个子问题，如下所示：

目标函数(13)，对于S的求解,利用梯度下降法，得：

S^k+1＝S^k-ξ·g (15)

其中ξ为下降步长，g为(13)式的导数，如下所示：

目标函数(14)对于d_x,d_y求解，可通过收缩算子求解：

步骤二：固定S，更新D和α_j，对应的目标函数如下：

字典学习分为字典更新与稀疏表示两个阶段。首先，固定α_j更新D，这里采用自适应字典，从重建的图像S中，使用K-SVD 方法学习字典。然后，固定D更新α_j，采用正交匹配追踪算法 (Orthogonal Matching Pursuit，OMP)用于更新稀疏表示系数α_j。

步骤三：当达到最大迭代次数k时，迭代停止，得到优化后的慢度值S^k，并将其转为二维慢度值S^k(x,y)。

S8、超压峰值场重建

将步骤S7中得到的二维慢度值S^k(x,y)根据式(4)转换为速度场V(x,y)，再根据式(3)转换为超压峰值Pm(x,y)，即冲击波超压峰值的空间分布。

S9、超压时空场重建

利用步骤S8中得到的超压峰值Pm(x,y)，结合“修正的弗里德兰德方程”作为超压-时间模型，得到冲击波超压随时间和空间的关系式：

式中，P(x,y,t)为位置处(x,y)冲击波在t时刻的超压值，P₀为环境压力，Pm(x,y)为超压峰值，T⁺为正相持续时间，c为衰减系数。下面对正相持续时间T⁺和衰减系数c进行求解。

(1)求解正相持续时间T⁺

已知正相持续时间T⁺和传感器到炸点的距离r、药量W有如下关系：

为了确定参数a和n，上式可转化为：

可进一步转化为：

xz₁+y＝z₂ (22)

其中，x＝n，

y＝lna，

将各个传感器到炸点的距离r_I、各个传感器测得的正相持续时间

及药量W代入式(22)，用最小二乘法进行拟合，得到x和 y，即参数a和n，得到T⁺与r和W的关系。

(2)求解衰减系数c

由于测试区域是依据J.Henrych提出的经验公式来划分为3 个区域，衰减系数c在整个场区不是恒定值，因此分区域计算该参数，下面以其中一个区域为例求解衰减系数c。

式(19)可写成参数c的表达式：

将此区域内各传感器S_I获取的在t时刻不同位置处(x,y)的超压值P(x,y,t)、超压峰值Pm(x,y)、正相持续时间T⁺等参数代入式 (23)，得到各个传感器对应得衰减系数c_I随时间t的曲线，且参数c和时间t存在如下关系：

kt+m＝ln ln c (24)

其中，k和m为待求系数。

在各条曲线上分别取若干个点，将各点对应的时间t及ln ln c 代入式(24)，采用最小二乘法进行拟合，得到参数k和m，即求得各条曲线对应得到衰减系数c_I。将3个区域对应的衰减系数c_I分别求均值，所得的c₁、c₂、c₃分别为三个区域的衰减系数。

在确定正相持续时间T⁺和衰减系数c后，并将其代入式(19)，得到具体的超压随时间和空间的关系式，即可重建冲击波超压时空场。

S10、超压时空场可视化

由当前爆炸区域各参数P₀、Pm(x,y)、T⁺、c，根据冲击波超压时空式(19)，可以得到测试区域冲击波超压的时空数据；将超压时空数据转换为线、面、体等的几何数据；即可得到测试区域冲击波超压的时空分布图。

基于python平台，通过matplotlib图形库函数，将几何数据进行归一化处理后转化为三维图像数据，对生成的图像进行平滑处理，选择一条合理的颜色带，进行三维可视化渲染。最后将爆炸冲击波超压随时间和空间的变化过程动态的显示出来，实现了爆炸冲击波超压时空场的动画可视化。

Claims

1.一种冲击波超压三维时空场重建方法，其特征在于，具体包括以下步骤：

S2、在划分好的每个区域边界处均布传感器；

S4、确定真实的冲击波超压峰值随距离的衰减规律模型；

S6、建立走时层析成像模型；

S7、结合字典学习与TV全变差，解算走时层析成像模型，得到二维慢度值；

2.根据权利要求1所述的一种冲击波超压三维时空场重建方法，其特征在于，S1具体为：根据冲击波超压峰值随距离衰减的经验公式中的比例距离

的区间，得到炸点到测点的距离r，根据该距离将测试区域划分成对应的3个区域，并将每个区域划分为若干个正方形网格单元，网格划分越密，重建的精度就越高：

其中，Pm为冲击波超压峰值，

为比例距离，

r为炸点到测点的距离，W为装药重量。

3.根据权利要求1所述的一种冲击波超压三维时空场重建方法，其特征在于，S4中，冲击波超压峰值Pm的计算的计算公式为：

其中，a、b、c为待求系数，根据炸点到测点的距离r以及

可以推倒出每个传感器对应的

将步骤S3中每个传感器对应的冲击波超压峰值Pm和传感器的坐标信息代入上述计算公式，并通过最小二乘法进行拟合，得到参数a、b、c，建立更适合此次试验药量条件下的冲击波超压峰值随距离的衰减规律模型。

4.根据权利要求1所述的一种冲击波超压三维时空场重建方法，其特征在于，S5中，冲击波超压峰值与冲击波速度的关系为：

其中，Pm(x,y)为(x,y)位置处的超压峰值，P₀为环境大气压力，V(x,y)为位置处(x,y)的冲击波速度，C₀为声速，

T₀为未扰动的空气初始温度；

冲击波速度V(x,y)与慢度S(x,y)的关系为：

由上式将超压场分布转换为慢度场分布。

5.根据权利要求1所述的一种冲击波超压三维时空场重建方法，其特征在于，S6中，假设爆炸波为连续速度场，爆炸波在传输的过程中，其走时是速度和几何路径的函数：

其中t为走时，即爆炸波传播到传感器节点的时间，v为速度，s为慢度，L为炸点到传感器的射线路径，dr为沿射线路径L的距离增量；

将(5)式离散化，对第i条射线的投影数据有：

式中，t_i为第i条射线的走时，即爆炸波到达第i个传感器的时间；s_j为将二维网格区域中的慢度S(x,y)转换为1维列向量后，对应的第j个网格中的慢度；a_ij为第i条射线穿过第j个网格的射线长度；M为射线数，即传感器个数；N为划分的网格数，其中N是指将二维网格区域转换为1维列向量；

将(6)式写成矩阵形式，得到走时层析成像模型为：

AS＝T

其中A为投影矩阵，S为重建图像，T为投影数据；T＝(t₁,t₂,…,t_M)'，即M维走时向量；S＝(S₁,S₂,…S_N)'，即N维慢度向量；A＝(a_ij)，即M×N阶射线矩阵。

6.根据权利要求1所述的一种冲击波超压三维时空场重建方法，其特征在于，S7具体为：

结合字典学习与TV全变差，对走时层析成像模型构建目标函数：

其中，第一项为数据的保真项，第二项为TV全变差正则项，第三项为字典学习正则项；μ为保真项系数，λ和β为正则项系数，E_i提取图像块的算子，v_j子字典学习正则化参数；上述系数在0-1之间。

7.根据权利要求6所述的一种冲击波超压三维时空场重建方法，其特征在于，μ＝1，λ＝0.001，β＝0.01，v_j取常数0.002。

8.根据权利要求6所述的一种冲击波超压三维时空场重建方法，其特征在于，根据步骤5得到的初始慢度值，求解目标函数(8)，具体步骤如下：

步骤一：固定D和α_j，更新重建图像S，首先引入辅助变量，将目标函数求解分为L1范数和L2范数两部分，然后采用Bregman算法求解；

步骤二：固定S，更新D和α_j，首先，固定α_j更新D，从重建的图像S中，使用K-SVD方法学习字典；然后，固定D更新α_j，采用正交匹配追踪算法用于更新稀疏表示系数α_j；

9.根据权利要求6所述的一种冲击波超压三维时空场重建方法，其特征在于，S9具体为：

利用步骤S8中得到的超压峰值Pm(x,y)，结合修正的弗里德兰德方程作为超压-时间模型，得到冲击波超压随时间和空间的关系式：

式中，P(x,y,t)为位置处(x,y)冲击波在t时刻的超压值，P₀为环境压力，Pm(x,y)为超压峰值，T⁺为正相持续时间，c为衰减系数；

正相持续时间T⁺和传感器到炸点的距离r、药量W的关系为：

xz₁+y＝z₂ (22)；

其中，x＝n，

y＝lna，

将各个传感器到炸点的距离r_I、各个传感器测得的正相持续时间T_I ⁺及药量W代入式(22)，用最小二乘法进行拟合，得到x和y，即参数a和n，即可得到T⁺与r和W的关系；

分区域计算参数c：

将此区域内各传感器S_I获取的在t时刻不同位置处(x,y)的超压值P(x,y,t)、超压峰值Pm(x,y)、正相持续时间T⁺代入式(23)，得到各个传感器对应得衰减系数c_I随时间t的曲线，且参数c和时间t存在如下关系：

kt+m＝lnlnc (24)

其中，k和m为待求系数；

在各条曲线上分别取若干个点，将各点对应的时间t及lnlnc代入式(24)，采用最小二乘法进行拟合，得到参数k和m，即求得各条曲线对应得到衰减系数c_I，将每个区域对应的衰减系数c_I分别求均值，分别得到每个区域的衰减系数；

10.根据权利要求9所述的一种冲击波超压三维时空场重建方法，其特征在于，还包括S10、将超压时空场可视化：由当前爆炸区域各参数P₀、Pm(x,y)、T⁺、c，根据超压随时间和空间的关系式，可以得到测试区域冲击波超压的时空数据；将超压时空数据转换为几何数据，即可得到测试区域冲击波超压的时空分布图。