A method for reconstructing overpressure field of explosion shock wave based on finite element and tomographic inversion

Abstract

The present invention provides a method for reconstructing the overpressure field of an explosion shock wave based on finite element and tomographic inversion. First, the finite element method is used to calculate the overpressure of the explosion shock wave. The numerical calculation result is used as the initial model for inversion and reconstruction. The weighted generalized reverse travel time tomography inversion algorithm is used to perform tomographic inversion and reconstruction of the shock wave overpressure in the test area. At the same time, the numerical calculation result is used to constrain the abnormal values ​​appearing in the inversion process. The experimental results show that the method for reconstructing the overpressure field of an explosion shock wave based on finite element and tomographic inversion proposed by the present invention can greatly improve the accuracy of inversion and reconstruction, and also improve the speed of iterative convergence.

Description

Technical Field

[0001] The invention relates to a method for reconstructing an explosion shock wave overpressure field based on finite element and tomographic inversion, and belongs to the field of array signal processing and reconstruction. Background Art

[0002] When explosives explode in the air, high-temperature, high-pressure, and high-speed explosion products will be generated in a very short time. The products will diffuse to the surroundings at a very high speed, and the surrounding medium will be directly affected by the products and strongly compressed. Therefore, the pressure, density, and temperature of the surrounding medium will suddenly increase, thus forming an initial shock wave. Shock wave overpressure is one of the main factors that cause damage and destructive effects of ammunition explosions on personnel, equipment, and protective structures. Therefore, the analysis and testing of shock wave overpressure is of great significance in the engineering field, especially in the military field. Exploring the transmission laws of explosion shock waves can effectively predict the destructive effects caused by such explosions, and can also be used to develop powerful weapons and ammunition; in civilian use, it can also be used for earthquake prevention and disaster reduction.

[0003] Existing explosion shock wave overpressure testing methods are based on a limited number of devices for cost and technical considerations. For uneven damage caused by non-spherical, asymmetric ammunition or fuel-air explosives, local testing methods cannot obtain all-round explosion information.

[0004] Since the current task of studying explosion resistance of structures or components is very arduous, and explosion tests are very expensive, the finite element numerical calculation method, as a supplement to theoretical research and experiments, can intuitively understand the propagation of explosion shock waves and the process of their effect on the target, and can also solve the overpressure peak, becoming an effective means of explosion resistance research. Although numerical calculation can describe and analyze the explosion process, in the application process, different values ​​of state equation parameters, material model parameters and mesh density often have a greater impact on the simulation results, resulting in the research results not being consistent with the actual test. Summary of the invention

[0005] In order to overcome the shortcomings of the prior art, the present invention provides a method for reconstructing the overpressure field of an explosion shock wave based on finite element and tomographic inversion. The finite element method is used to establish a numerical calculation model of the explosion, and the explosion shock wave overpressure is calculated. Based on the principle of tomography, a travel-time tomography algorithm is used to invert and reconstruct the shock wave overpressure in the test area. During the inversion and reconstruction, the result of the finite element numerical calculation is used as the initial value of the tomographic inversion, and the numerical calculation result is used to constrain the inversion process, so as to achieve the purpose of greatly improving the inversion and reconstruction accuracy.

[0006] The technical solution adopted by the present invention to solve the technical problem is: a method for reconstructing the overpressure field of explosion shock wave based on finite element and tomographic inversion is provided, comprising the following steps:

[0007] S1. Grid division of the explosion shock wave overpressure field test area;

[0008] S2. Use the finite element modeling method to establish an explosion model and obtain the shock wave overpressure time history curve. The shock wave peak overpressure of each grid in the explosion shock wave overpressure field test area is directly obtained through the shock wave overpressure time history curve. The shock wave peak overpressure is converted into shock wave velocity and the inverse of the velocity is obtained as the slowness value. The initial model S of the tomography inversion model is established according to the slowness value of each grid. 0 ;

[0009] S3. Use the weighted generalized reverse travel time tomography method to invert the peak overpressure of the explosion shock wave, output a tomography inversion model that meets the convergence criterion, and complete the reconstruction of the explosion shock wave overpressure field.

[0010] Step S1 uses a multi-scale grid to divide the explosion shock wave overpressure field test area into grids.

[0011] Step S2 uses a point detonation method to establish a finite element model, and the finite element model includes an explosive column part and an air system part.

[0012] Step S2 specifically includes the following process:

[0013] S2.1. Use the finite element modeling method to establish a finite element model, which includes an explosive column and an air system; use the JWL state equation to describe the explosive column:

[0014]

[0015] Where P is the detonation pressure, V is the relative volume, E0 is the initial internal energy per unit volume, ω is the first material constant, A is the second material constant, B is the third material constant, R1 is the fourth material constant, and R2 is the fifth material constant;

[0016] The air system uses the MAT-NULL material model and is described by the following state equation:

[0017] P=C0+C1μ+C2μ 2 +C3μ 3 +(C4+C5μ+C6μ 2 )E (2)

[0018] Where C0 is the first gas constant, C1 is the second gas constant, C2 is the third gas constant, C3 is the fourth gas constant, C4 is the fifth gas constant, C5 is the sixth gas constant, and C6 is the seventh gas constant; C0 = C1 = C2 = C3 = C6 = 0, C4 = C5 = γ-1, is the dynamic viscosity coefficient, ρ is the current density, ρ0 is the initial density, ρ0=1.3×10-3 g / cm 3 , γ=1.4, is the adiabatic index, E=2.5×10 5 Pa, is the internal energy per unit volume;

[0019] S2.2, numerical simulation calculation of explosion process is performed according to finite element model, equation (1) and equation (2);

[0020] S2.3. Extract the simulation calculation results to obtain the overpressure time history curves at different proportional distances, and read the shock wave peak overpressure of each grid in the test area from the curves;

[0021] S2.4, converting the shock wave peak overpressure of each grid into shock wave velocity according to the relationship between the shock wave peak overpressure and velocity;

[0022] S2.5. The inverse of the shock wave velocity of each grid is taken as the slowness value, and the slowness value of each grid is used as the initial model S of the tomography inversion model. 0 .

[0023] In step S2.4, the shock wave peak overpressure is converted to shock wave velocity by the following formula:

[0024]

[0025] Where c is the shock wave velocity, p m is the peak overpressure of the shock wave, p0 is the initial pressure of the undisturbed air, obtained by the sensor test during the test, C air is the speed of sound in undisturbed air, T0 is the initial temperature of the undisturbed air, which is obtained by the sensor test during the test.

[0026] Step S3 performs tomographic inversion of the peak overpressure of the shock wave in the explosion field, which specifically includes the following process:

[0027] S3.1. Initial assignment: The tomography inversion model S is initially assigned to the initial model S 0 , that is, S = S 0 ;

[0028] S3.2. According to the principle of travel time tomography,

[0029] DS=T (4)

[0030] Where T=(t1,t2…t M )' is the M-dimensional column vector of the travel time of each ray, obtained from experimental tests; S = (s1, s2…s N)' is the discrete unit slowness vector to be solved, which is an N-dimensional unknown column vector, that is, the final goal of inverting the tomography inversion model; D is the distance matrix, which is an M×N sparse matrix, where the elements are d ij , i.e., the ray length of the i-th ray passing through the j-th grid, obtained in the process of meshing the explosion shock wave overpressure field test area in step S1;

[0031] The data weighted matrix P is calculated, and its diagonal elements are:

[0032] diag{P}=T p -1 (5)

[0033] The model weight matrix QP is calculated, and its diagonal elements are:

[0034] diag{Q}=K (6)

[0035] The matrix T p The elements of are the travel times of each ray of the current model, that is, The matrix K is the total contribution of each ray on the jth grid, and its elements are the product of the length of the ray passing through the grid unit and the velocity of the corresponding unit, that is, Matrix T p Both and matrix K are diagonal matrices and symmetric positive definite;

[0036] Calculate the current tomography inversion model S h :

[0037] S h =Q -1 (PDQ -1 ) + P.T. (7)

[0038] The number of iterations h is initially assigned a value of 1;

[0039] S3.3, iterative calculation and finite element numerical calculation are used to constrain the inversion outliers, that is, to judge or Is it established? j h is the slowness value obtained at the hth iteration of the jth grid, s min and max are the minimum and maximum slowness values ​​set for the test, respectively;

[0040] If established, then The initial model S 0 The slowness value of the jth grid in , and then go to step S3.4; otherwise, directly go to step S3.4;

[0041] S3.4. Use the following convergence criteria to judge the convergence of iterative calculations:

[0042] Convergence criterion a, the mean square relative change between the h-th iteration and the h-1-th iteration solution estimate is less than a certain value ε1, ε1≤0.005, that is

[0043]

[0044] Convergence criterion b: travel time residual is less than ε2, ε2≤0.005, that is

[0045]

[0046] where t i is the measured travel time data of the i-th M-dimensional column vector T, d i is the i-th row vector of the distance matrix D, S h is the discrete unit slowness vector obtained at the hth iteration;

[0047] If the iteration satisfies the convergence criterion a or convergence criterion b, the iteration ends and the discrete unit slowness vector S = S is output. h , go to step S3.5; otherwise, assign h+1 to h, and return to step S3.2 for the next iterative calculation;

[0048] S3.5. The shock wave velocity is obtained based on the inverse relationship between slowness and speed, and then the shock wave peak overpressure in the test area is obtained based on the relationship between shock wave velocity and overpressure.

[0049] Step S3.5 uses formula (10) to convert the shock wave velocity into the shock wave peak overpressure in the test area:

[0050]

[0051] In the formula, p m is the peak overpressure of the shock wave, c is the shock wave velocity, and p0 is the initial pressure of the undisturbed air; c air is the undisturbed air sound speed, for different temperatures: T0 is the initial temperature of the undisturbed air.

[0052] The beneficial effects of the present invention based on its technical solution are:

[0053] The present invention provides a method for reconstructing an explosion shock wave overpressure field based on finite element and tomographic inversion. The method firstly applies the finite element method to calculate the explosion shock wave overpressure, uses the numerical calculation result as the initial model for inversion and reconstruction, and adopts the weighted generalized reverse traveltime tomography inversion algorithm to perform tomographic inversion and reconstruction on the shock wave overpressure in the test area. Meanwhile, the finite element numerical calculation result is used to constrain the abnormal values ​​appearing in the inversion process. Experimental results show that the method for reconstructing an explosion shock wave overpressure field based on finite element and tomographic inversion proposed by the present invention can greatly improve the accuracy of inversion and reconstruction, and also improve the speed of iterative convergence. BRIEF DESCRIPTION OF THE DRAWINGS

[0054] Figure 1 It is a schematic diagram of the grid division of the explosion shock wave overpressure field test area.

[0055] Figure 2 It is a schematic diagram of the finite element calculation model.

[0056] Figure 3 It is a schematic diagram of the shock wave overpressure time history curve of some points in the test area calculated by LS-DYNA, where the horizontal axis represents time / second and the vertical axis represents pressure / Pa.

[0057] Figure 4 It is the initial model of tomographic inversion shock wave velocity.

[0058] Figure 5 It is the inversion calculation flow chart.

[0059] Figure 6 is a schematic diagram of tomographic inversion (one quarter area).

[0060] Figure 7 This is a schematic diagram of the two-dimensional distribution of the inverted shock wave peak overpressure. DETAILED DESCRIPTION

[0061] The present invention will be further described below in conjunction with the accompanying drawings and embodiments.

[0062] The present invention provides a method for reconstructing an explosion shock wave overpressure field based on finite element and tomographic inversion, comprising the following steps:

[0063] S1. Grid the explosion shock wave overpressure field test area.

[0064] In this step, a multi-scale grid method can be used to grid the explosion shock wave overpressure field test area. For specific steps, please refer to the patent publication with publication number CN113094962A. The grid division result is as follows: Figure 1 shown.

[0065] S2. Use the finite element modeling method to establish an explosion model and obtain the shock wave overpressure time history curve. The shock wave peak overpressure of each grid in the explosion shock wave overpressure field test area is directly obtained through the shock wave overpressure time history curve. The shock wave peak overpressure is converted into shock wave velocity and the inverse of the velocity is obtained as the slowness value. The initial model S of the tomography inversion model is established according to the slowness value of each grid. 0 The specific process includes the following:

[0066] S2.1. This step can use the 8-node entity unit of the finite element software LS-DYNA to establish the model, and use the point detonation method to establish the finite element model, which includes the explosive column part and the air system part. In order to avoid severe distortion of the grid, the explosive and air units use the Euler algorithm. For the experimental 50kgTNT charge explosion, the explosive radius is set to 10cm; the air system radius is 1600cm, and the calculation model is evenly divided into 19023 units, including 4 explosive units, and the units use the single-point ALE algorithm.

[0067] The JWL state equation is used to describe the explosive column:

[0068]

[0069] Where P is the detonation pressure, V is the relative volume, E0 is the initial internal energy per unit volume, ω is the first material constant, A is the second material constant, B is the third material constant, R1 is the fourth material constant, and R2 is the fifth material constant. The values ​​of each parameter refer to the following table:

[0070] A(GPa) B(GPa) ω <![CDATA[R1]]> <![CDATA[R2]]> <![CDATA[E0(GPa)]]> <![CDATA[V0]]> 373.8 3.747 0.35 4.15 0.9 6.0 1.0

[0071] Table 1 JWL state equation parameters

[0072] The air system uses the MAT-NULL material model and is described by the following state equation:

[0073] P=C0+C1μ+C2μ 2 +C3μ 3 +(C4+C5μ+C6μ 2 )E (2)

[0074] Where C0 is the first gas constant, C1 is the second gas constant, C2 is the third gas constant, C3 is the fourth gas constant, C4 is the fifth gas constant, C5 is the sixth gas constant, and C6 is the seventh gas constant; C0 = C1 = C2 = C3 = C6 = 0, C4 = C5 = γ-1, is the dynamic viscosity coefficient, ρ is the current density, ρ0 is the initial density, ρ0=1.3×10 -3 g / cm 3 , γ=1.4, is the adiabatic index, E=2.5×105 Pa is the internal energy per unit volume.

[0075] The schematic diagram of the finite element model is as follows Figure 2 shown.

[0076] S2.2, Combining the finite element model, equation (1) and equation (2) to perform numerical simulation of the explosion process, the results can be calculated in LS-DYNA finite element software as follows Figure 3 The shock wave overpressure time history curve is shown.

[0077] S2.3. Extract the simulation calculation results to obtain the overpressure time history curves at different proportional distances, and read the shock wave peak overpressure of each grid in the test area from the curves.

[0078] S2.4. According to the relationship between the peak overpressure and velocity of the shock wave, the peak overpressure of the shock wave of each grid is converted into the shock wave velocity by the following formula:

[0079]

[0080] Where c is the shock wave velocity, p m is the peak overpressure of the shock wave, p0 is the initial pressure of the undisturbed air, and p0=0.1Mp was measured during the experiment. air is the speed of sound in undisturbed air, T0 is the initial temperature of the undisturbed air (K), and the experimental result shows that T0 = 297 (K). The initial model of the tomographic inversion shock wave velocity is as follows: Figure 4 shown.

[0081] S2.5. The inverse of the shock wave velocity of each grid is taken as the slowness value, and the slowness value of each grid is used as the initial model S of the tomography inversion model. 0 .

[0082] S3. Use the weighted generalized reverse travel time tomography method to invert the peak overpressure of the explosion shock wave, output a tomography inversion model that meets the convergence criteria, and complete the reconstruction of the explosion shock wave overpressure field. Figure 5 , specifically including the following processes:

[0083] S3.1. Initial assignment: The tomography inversion model S is initially assigned to the initial model S 0 , that is, S = S 0 .

[0084] S3.2. According to the principle of travel time tomography,

[0085] DS=T (4)

[0086] Where T=(t1,t2…t M)' is the M-dimensional column vector of the travel time of each ray, obtained from experimental tests; S = (s1, s2…s N )' is the discrete unit slowness vector to be solved, which is an N-dimensional unknown column vector, that is, the final goal of inverting the tomography inversion model; D is the distance matrix, which is an M×N sparse matrix, where the elements are d ij , that is, the length of the ray of the i-th ray passing through the j-th grid, is obtained in the process of meshing the explosion shock wave overpressure field test area in step S1; the schematic diagram of the tomographic inversion of the quarter area is shown in Figure 6 shown.

[0087] The data weighted matrix P is calculated, and its diagonal elements are:

[0088] diag{P}=T p -1 (5)

[0089] The model weight matrix QP is calculated, and its diagonal elements are:

[0090] diag{Q}=K (6)

[0091] The matrix T p The elements of are the travel times of each ray of the current model, that is, The matrix K is the total contribution of each ray on the jth grid, and its elements are the product of the length of the ray passing through the grid unit and the velocity of the corresponding unit, that is, Matrix T p Both and matrix K are diagonal matrices and symmetric positive definite;

[0092] Calculate the current tomography inversion model S h :

[0093] S h =Q -1 (PDQ -1 ) + P.T. (7)

[0094] The number of iterations h is initially assigned a value of 1, and the value of T is the shock wave arrival time measured by each sensor during the test.

[0095] S3.3, iterative calculation and finite element numerical calculation are used to constrain the inversion outliers, that is, to judge or Is it established? j h is the slowness value obtained at the hth iteration of the jth grid, s min and max are the minimum and maximum slowness values ​​set for the test, respectively;

[0096] If established, then The initial model S 0 The slowness value of the j-th grid in is obtained, and then the process goes to step S3.4; otherwise, the process goes directly to step S3.4.

[0097] S3.4. Use the following convergence criteria to judge the convergence of iterative calculations:

[0098] Convergence criterion a, the mean square relative change between the h-th iteration and the h-1-th iteration solution estimate is less than a certain value ε1, ε1≤0.005, that is

[0099]

[0100] in, is the slowness vector obtained at the hth iteration of the jth cell, is the slowness vector obtained at the h-1th iteration of the jth cell.

[0101] Convergence criterion b: travel time residual is less than ε2, ε2≤0.005, that is

[0102]

[0103] where t i is the measured travel time data of the i-th M-dimensional column vector T, d i is the i-th row vector of the distance matrix D, S h is the discrete unit slowness vector obtained at the hth iteration;

[0104] If the iteration satisfies the convergence criterion a or convergence criterion b, the iteration ends and the discrete unit slowness vector S = S is output. h , go to step S3.5; otherwise, assign h+1 to h, and return to step S3.2 for the next iterative calculation;

[0105] S3.5. The shock wave velocity is obtained based on the inverse relationship between slowness and speed, and then the shock wave peak overpressure in the test area is obtained based on the relationship between shock wave velocity and overpressure. The specific formula is as follows:

[0106]

[0107] In the formula, p m is the peak overpressure of the shock wave, c is the shock wave velocity, and p0 is the initial pressure of the undisturbed air; p0 = 0.1 Mp was measured in the experiment, c air is the undisturbed air sound speed, for different temperatures: T0 is the initial temperature of the undisturbed air. The experimental result shows that T0 = 297 (K).

[0108] The present invention can be verified by the following tests:

[0109] The test used 50kg of bulk TNT explosives, suspended 1.7 meters above the ground, and exploded in the air. The density of the explosives was 0.80g / cm 3 , the shape is close to spherical. Overpressure sensors are placed within 16 meters of the explosion center to conduct shock wave peak overpressure tests. The shock wave arrival time T measured by the sensor is substituted into formula (7) to perform overpressure field inversion and reconstruction. The test results are used as standard values. The explosion overpressure field numerical calculation method based on finite element software LS-DYNA, the weighted generalized reverse travel time tomography inversion algorithm based on the random initial model, and the inversion and reconstruction algorithm based on the combination of finite element numerical calculation and tomography inversion proposed in this patent are used to perform overpressure field inversion and reconstruction. The reconstruction results are shown in Table 2. Where ε1=ε2=0.001, s min =0.2s / km, s max =2.8s / km.

[0110]

[0111] Table 2 Inversion of test shock wave peak overpressure and test results

[0112] The relative error of overpressure of the jth grid cell is defined as:

[0113]

[0114] Among them, P′ j is the overpressure value reconstructed (finite element calculation), P j It is the experimental test value.

[0115] The number of iterations of the inversion algorithm using the random initial model and the finite element numerical calculation results proposed in this patent as the initial model is shown in Table 3.

[0116]

[0117] Table 3 Comparison of the number of iterations of the inversion algorithm under different initial models

[0118] It can be seen that the relative error of the inversion algorithm proposed in the present invention is small, and the iterative convergence speed is faster. The two-dimensional distribution of the peak overpressure of the explosion shock wave reconstructed by the method of the present invention is as follows: Figure 7 shown.

[0119] The present invention provides a method for reconstructing an explosion shock wave overpressure field based on finite element and tomographic inversion. The method first applies the finite element method to calculate the explosion shock wave overpressure, uses the numerical calculation result as the initial model for inversion and reconstruction, and adopts a weighted generalized reverse traveltime tomography inversion algorithm to perform tomographic inversion and reconstruction on the shock wave overpressure in a test area. At the same time, the numerical calculation result is used to constrain abnormal values ​​appearing in the inversion process. Experimental results show that the method for reconstructing an explosion shock wave overpressure field based on finite element and tomographic inversion proposed by the present invention can greatly improve the accuracy of inversion and reconstruction, and also improve the speed of iterative convergence.

Claims

1. A method for reconstructing the overpressure field of explosion shock waves based on finite element and tomographic inversion, characterized in that The following steps are involved: S1. Grid division of the explosion shock wave overpressure field test area; S2. Use the finite element modeling method to establish an explosion model and obtain the shock wave overpressure time history curve. The shock wave peak overpressure of each grid in the explosion shock wave overpressure field test area is directly obtained through the shock wave overpressure time history curve. The shock wave peak overpressure is converted into shock wave velocity and the inverse of the velocity is obtained as the slowness value. The initial model S of the tomography inversion model is established according to the slowness value of each grid. 0 ; The finite element model is established by point detonation. The finite element model includes the explosive column part and the air system part, which specifically includes the following processes: S2.

1. Use the finite element modeling method to establish a finite element model, which includes an explosive column and an air system; use the JWL state equation to describe the explosive column: Where P is the detonation pressure, V is the relative volume, E0 is the initial internal energy per unit volume, ω is the first material constant, A is the second material constant, B is the third material constant, R1 is the fourth material constant, and R2 is the fifth material constant; The air system uses the MAT-NULL material model and is described by the following state equation: P=C0+C1μ+C2μ 2 +C3μ 3 +(C4+C5μ+C6μ 2 )E (2) Where C0 is the first gas constant, C1 is the second gas constant, C2 is the third gas constant, C3 is the fourth gas constant, C4 is the fifth gas constant, C5 is the sixth gas constant, and C6 is the seventh gas constant; C0 = C1 = C2 = C3 = C6 = 0, C4 = C5 = γ-1, is the dynamic viscosity coefficient, ρ is the current density, ρ0 is the initial density, ρ0=1.3×10 -3 g / cm 3 , γ=1.4, is the adiabatic index, E=2.5×10 5 Pa, is the internal energy per unit volume; S2.2, numerical simulation calculation of explosion process is performed according to finite element model, equation (1) and equation (2); S2.

3. Extract the simulation calculation results to obtain the overpressure time history curves at different proportional distances, and read the shock wave peak overpressure of each grid in the test area from the curves; S2.

4. According to the relationship between the shock wave peak overpressure and velocity, the shock wave peak overpressure of each grid is converted into shock wave velocity by the following formula: Where c is the shock wave velocity, p m is the peak overpressure of the shock wave, p0 is the initial pressure of the undisturbed air, obtained by the sensor test during the test, C air is the speed of sound in undisturbed air, T0 is the initial temperature of the undisturbed air, obtained by the sensor test during the test; S2.

5. The inverse of the shock wave velocity of each grid is taken as the slowness value, and the slowness value of each grid is used as the initial model S of the tomography inversion model. 0 ; S3, using the weighted generalized reverse travel time tomography method to invert the peak overpressure of the explosion shock wave, output a tomography inversion model that meets the convergence criterion, and complete the reconstruction of the explosion shock wave overpressure field, which specifically includes the following processes: S3.

1. Initial assignment: The tomography inversion model S is initially assigned to the initial model S 0 , that is, S = S 0 ; S3.

2. According to the principle of travel time tomography, DS=T (4) Where T=(t1,t2…t M )' is the M-dimensional column vector of the travel time of each ray, obtained from experimental tests; S = (s1, s2…s N )' is the discrete unit slowness vector to be solved, which is an N-dimensional unknown column vector, that is, the final goal of inverting the tomography inversion model; D is the distance matrix, which is an M×N sparse matrix, where the elements are d ij , i.e., the ray length of the i-th ray passing through the j-th grid, obtained in the process of meshing the explosion shock wave overpressure field test area in step S1; The data weighted matrix P is calculated, and its diagonal elements are: diag{P}=T p -1 (5) The model weight matrix QP is calculated, and its diagonal elements are: diag{Q}=K (6) The matrix T p The elements of are the travel times of each ray of the current model, that is, The matrix K is the total contribution of each ray on the jth grid, and its elements are the product of the length of the ray passing through the grid unit and the velocity of the corresponding unit, that is, Matrix T p The matrix K is diagonal and symmetric and positive definite; Calculate the current tomography inversion model S h : S h =Q -1 (PDQ -1 ) + P·T (7) The number of iterations h is initially assigned a value of 1; S3.3, iterative calculation and finite element numerical calculation are used to constrain the inversion outliers, that is, to judge or Is it established? j h is the slowness value obtained at the hth iteration of the jth grid, s min and max are the minimum and maximum slowness values ​​set for the test, respectively; If established, then The initial model S 0 The slowness value of the jth grid in , and then go to step S3.4; otherwise, directly go to step S3.4; S3.

4. Use the following convergence criteria to judge the convergence of iterative calculations: Convergence criterion a, the mean square relative change between the h-th iteration and the h-1-th iteration solution estimate is less than a certain value ε1, ε1≤0.005, that is Convergence criterion b: travel time residual is less than ε2, ε2≤0.005, that is where t i is the measured travel time data of the i-th M-dimensional column vector T, d i is the i-th row vector of the distance matrix D, S h is the discrete unit slowness vector obtained at the hth iteration; If the iteration satisfies the convergence criterion a or convergence criterion b, the iteration ends and the discrete unit slowness vector S = S is output. h , go to step S3.5; otherwise, assign h+1 to h, and return to step S3.2 for the next iterative calculation; S3.

5. The shock wave velocity is obtained according to the inverse relationship between slowness and speed, and then the shock wave peak overpressure in the test area is obtained according to the relationship between shock wave velocity and overpressure. Specifically, the shock wave velocity is converted into the shock wave peak overpressure in the test area using formula (10): In the formula, p m is the peak overpressure of the shock wave, c is the shock wave velocity, and p0 is the initial pressure of the undisturbed air; c air is the undisturbed air sound speed, for different temperatures: T0 is the initial temperature of the undisturbed air.

2. The method for reconstructing the overpressure field of explosion shock wave based on finite element and tomographic inversion according to claim 1 is characterized in that: Step S1 uses a multi-scale grid to divide the explosion shock wave overpressure field test area into grids.

Images