CN102967189B - Explosive blast overpressure space-time field reconstruction method - Google Patents

Abstract

A method for reconstructing the space-time field of an explosion shock wave. The overpressure sensor array in the shock wave test system is used to collect shock wave signals, and the overpressure space-time field reconstruction module of the shock wave is used to reconstruct the space-time field of the explosion shock wave. Then, the explosion point is located, and then the shock wave velocity field and the shock wave peak overpressure field are reconstructed to obtain the reconstruction result of the shock wave peak overpressure field; combined with the reconstruction result of the shock wave peak overpressure field at this point, the " The parameters of the "modified Friedlander equation" are calculated to obtain the shock wave overpressure space-time data; finally, the reconstruction data is converted into an image file to realize the visualization of the shock wave overpressure space-time field. The invention provides a basis for evaluating the damage efficiency of the shock wave and the structure and performance of the ammunition; meanwhile, it is also an important means for appraising the safety of personnel, the shock resistance of equipment, and evaluating the impact of the shock wave on the surrounding environment.

Description

Reconstruction Method of Spatio-temporal Field of Explosion Shock Wave Overpressure

technical field

The invention belongs to the field of array signal processing and reconstruction, and in particular relates to a reconstruction method of explosion shock wave overpressure space-time field.

Background technique

With the advancement of society and science and technology, and the needs of national economic development, explosive technology is increasingly used in various fields such as transportation, construction of water conservancy and hydropower facilities, and the process of moving mountains and reclamation. While the explosion technology is widely used, the impact and harm of the explosion operation on the surrounding environment and building facilities have become the focus of attention. How to assess the power of the explosion, control the hazards of the explosion, and reduce the impact of the explosion operation on the surrounding environment are important contents of the explosion safety technology research.

When the charge explodes in the air, a high-temperature and high-pressure gas is formed, which violently pushes the surrounding still air, and at the same time generates a series of compression waves that propagate to the surroundings, and finally form an air shock wave. The amount by which the peak pressure of the air shock wave front exceeds the normal air pressure is called shock wave overpressure. At the moment when the shock wave arrives, the air pressure suddenly rises from the initial pressure to the peak overpressure, and then the expansion of the detonation products gradually slows down, and the movement speed of the high-pressure gas decreases, so that the pressure of the shock wave rises continuously, and the overpressure drops to zero. Then there is a lower pressure than the surrounding gas. Shock wave overpressure is the main factor that produces killing and destructive effects, and is an important part of the energy that damages the target. The damage of the shock wave to the target is the result of the joint action of the shock wave overpressure field in space and time. As long as the attenuation law of shock wave overpressure with space position and time is known, the destructive effect of shock wave at different locations and at different times can be determined. Shock wave overpressure is an important characteristic quantity to measure the damage of shock wave to people and objects. For personnel engaged in the production and use of ammunition, the destructive effect of shock wave must be fully considered, especially when using explosion shock waves to operate and destroy a large number of ammunition The overpressure of the shock wave must be fully estimated to determine the safe distance without danger. At present, for attempted bombing cases, the procuratorate often requires the public security organs to make an appraisal of the power and destructive effect of the explosion as a basis for sentencing the suspect. Therefore, it is necessary to reconstruct the space-time field of explosion shock wave overpressure, that is, to obtain the distribution of overpressure field and the attenuation law of shock wave overpressure with space and time, so as to meet the needs of safe production and use.

At present, although relevant scholars at home and abroad have summarized some formulas for calculating shock wave overpressure, these formulas are all empirical formulas. Due to the fact that the test environment, test conditions, charge shape and test data results are not complete when conducting actual experimental research The same, therefore, makes direct use of these data and research results difficult.

Explosion of explosives is completed in a very short period of time, and the resulting shock wave propagates in all directions. For symmetrical explosives such as spherical explosives, the air shock wave generated by the explosion propagates outward in the form of spherical waves, and the shock wave distribution after the explosion of asymmetric explosives is different. Therefore, in order to evaluate the damage of shock wave, it is necessary to study the law of shock wave overpressure transmission and its distribution. The number of sensors used in the test is limited, and the space where the sensors are arranged is also limited. During a static explosion, if the sensor is very close to the explosion point, the sensor will be damaged; The layout is also subject to certain restrictions; therefore, only shock wave parameters at limited locations can be measured, and these measured values are far from enough for evaluating shock wave damage and analyzing ammunition performance.

Contents of the invention

The purpose of the present invention is to provide a reconstruction method for explosion shock wave overpressure space-time field in order to solve the deficiency that the existing local distribution test and intensity test cannot fully understand the shock wave propagation process. The distribution of the explosion shock wave peak overpressure field in the region and the attenuation law of the shock wave overpressure with distance and time have realized the reconstruction of the space-time field of the explosion shock wave overpressure, and the shock wave overpressure generated by the explosion with time and space can be clearly imaged Shown, the visualization of the space-time field of shock wave overpressure is realized.

Technical scheme of the present invention is:

A reconstruction method of explosion shock wave overpressure space-time field, using a shock wave test system to realize the reconstruction of explosion shock wave overpressure space-time field, the shock wave test system consists of a shock wave overpressure sensor array, a signal conditioning module, a synchronous A/D conversion module, and a high-speed data transmission module , a data preprocessing module, a shock wave overpressure spatiotemporal field reconstruction module, and a computer, wherein the shock wave overpressure sensor array converts the shock wave overpressure signal into an electrical signal, which is conditioned into a standard voltage signal by a signal conditioning module and a synchronous A/D conversion module And convert it into a digital signal, and then transmit the digital signal to the computer through the high-speed data transmission module. The data preprocessing module and the overpressure space-time field reconstruction module are installed on the computer, and the shock wave overpressure signal after A/D conversion is input into the data preprocessing The processing module is input to the shock wave overpressure space-time field reconstruction module after data preprocessing, and the shock wave overpressure space-time field reconstruction module completes the reconstruction of the explosion shock wave overpressure space-time field; it is characterized in that the space-time field reconstruction method is: first The area is divided into grids, and the sensor array elements are arranged; secondly, the time difference method in the passive positioning algorithm is combined with the iterative algorithm to locate the explosion point, and then the shock wave peak value is determined based on the relationship between the travel time and the ray path and the relationship between the shock wave velocity and the peak overpressure. Reconstruct the overpressure field; then study the law and characteristics of shock wave overpressure attenuation with time at a fixed location in space, and combine the shock wave peak overpressure at this point to obtain the parameters of the "modified Friedlander equation" to obtain the shock wave overpressure space-time data; at last, transform the reconstruction data into an image file and display it on the monitor of the computer to realize the visualization of the shock wave overpressure space-time field; wherein: the reconstruction of the shock wave peak overpressure field consists of velocity field reconstruction and peak overpressure field reconstruction. Partial composition.

The specific method for reconstructing the shock wave peak overpressure field: extracting the measured travel time of the shock wave, optimizing the initial slowness model, calculating the ray path, setting the maximum number of iterations, calculating the theoretical travel time, inversion calculation of the slowness field, and conversion of the shock wave peak overpressure, The reconstruction result p _m (x, y) of the shock wave peak overpressure field is obtained.

The specific steps of the space-time field reconstruction method are as follows:

(1) During the grid division of the test area and the layout of the sensor array elements, the measurement area is divided into several regular grid units, and the sensor array elements are distributed on the nodes on the edge of the grid and on individual nodes in the grid according to the number of test use. The explosion point is located in the center of the test area;

(2) According to the relationship between the signal time delay and the signal speed, the distance difference is obtained, combined with the sensor position coordinates, an equation is established for each sensor, and then the iterative algorithm is used to calculate the spatial position of the explosion point; by measuring the shock wave reaching the sensor array element The relationship between the time and the ray path is reversed to obtain the shock wave velocity of each grid unit. According to the relationship between the shock wave velocity and the peak overpressure, the reconstruction result of the shock wave peak overpressure field is obtained, that is, the shock wave peak overpressure p _m (x,y );

(3) According to the time-varying curve of the shock wave overpressure at a fixed point obtained by the test sensor during the test, select the appropriate shock wave overpressure-time transmission model, and select the "modified Friedlander equation" by fitting comparison Combined with the reconstruction results p _m (x, y) of shock wave peak overpressure field, the relationship of shock wave overpressure with time and space is obtained, and the reconstruction result of shock wave overpressure space-time field is obtained, that is, the shock wave overpressure value p at any position and time (x,y,t);

(4) Establish the mapping between overpressure reconstruction data and color or geometric shapes, complete the transformation of application data to geometric element data such as points, lines, surfaces, and volumes; convert geometric data into image data, and generate image files from image data , and then smooth the generated image to eliminate individual abnormal points that do not conform to the law of shock wave transmission, and finally dynamically display the shock wave overpressure generated by the explosion with the change process of time and space position.

Outstanding feature and remarkable effect of the present invention are:

By using the shock wave overpressure space-time field reconstruction module in the shock wave test system to reconstruct various space-time parameters in the explosion process (such as the spatial position of the explosive point of the ammunition, the shock wave transmission speed, the shock wave overpressure distribution, the shock wave overpressure decay law with time and space, etc.) , to provide a basis for evaluating the damage effectiveness of the shock wave and the structure and performance of the ammunition; at the same time, it is also an important means to evaluate the safety of personnel, the shock resistance of equipment, and the impact of shock waves on the surrounding environment. Evaluation, tunnel explosion overpressure field reconstruction and power assessment provide a technical basis.

Description of drawings

Fig. 1 is the schematic diagram of shock wave testing system of the present invention;

Fig. 2 Flow chart of explosion shock wave overpressure space-time field reconstruction;

Fig. 3 is a schematic diagram of cell division in the present invention;

Fig. 4 is a schematic diagram of sensor layout of the present invention;

Fig. 5 is a flowchart of velocity field reconstruction in the present invention;

Fig. 6 is the shock wave pressure-time curve of the present invention;

Fig. 7 is a three-dimensional display diagram of the reconstruction result of the air explosion shock wave overpressure field.

Detailed ways

Taking the two-dimensional reconstruction of the overpressure field of the air explosion shock wave as an example, the present invention will be further described in detail in conjunction with the accompanying drawings.

As shown in Figure 1, the shock wave test system consists of a shock wave overpressure sensor array, a signal conditioning module, a synchronous A/D conversion module, a high-speed data transmission module, a data preprocessing module, a shock wave overpressure space-time field reconstruction module, and a computer. The shock wave overpressure sensor array converts the shock wave overpressure signal into an electrical signal, and through the signal conditioning module and the synchronous A/D conversion module, it is conditioned into a standard voltage signal and converted into a digital signal, and then the digital signal is transmitted to the computer through a high-speed data transmission module , the data preprocessing module and the overpressure space-time field reconstruction module are installed on the computer, and the shock wave overpressure signal after A/D conversion is input into the data preprocessing module. The pressure space-time field reconstruction module completes the reconstruction of the explosion shock wave overpressure space-time field.

As shown in Figure 2, the two-dimensional reconstruction of the air explosion shock wave overpressure field first divides the test area into a grid, and arranges the sensor array elements, and then uses the time difference method in the passive positioning algorithm combined with the iterative algorithm to locate the explosion point, and then based on The relationship between the travel time and the ray path and the relationship between the shock wave velocity and the peak overpressure are used to reconstruct the shock wave velocity field and the shock wave peak overpressure field; the law and characteristics of the shock wave overpressure decay with time are studied, and combined with the reconstruction results of the shock wave peak overpressure at this point, the The parameters of the "modified Friedlander equation" are calculated to obtain the shock wave overpressure space-time data; finally, the reconstruction data is converted into an image file to realize the visualization of the shock wave overpressure space-time field.

1) Mesh division of test area and layout of sensor array elements

In the two-dimensional reconstruction of the air explosion shock wave overpressure field, the explosion point is located in the center of the test area, and the sensor array elements and the explosion point are arranged on the same plane. In order to ensure the validity of the reconstruction data, the grid nodes at each edge of the reconstruction area The test sensors are arranged as much as possible so that the layout of the sensors can cover the measured area in as many directions as possible.

The grid coordinate system is shown in Figure 3. The first grid unit vertex on the lower left is taken as the origin O(0,0) of the coordinate system, and the coordinate of the explosion point is K(x ₀ ,y ₀ ), and the i-th ray The corresponding sensor coordinates are s _i ( _xi , y _i ), and 50KgTNT explosives are used to explode in the test. The sensors are arranged within the range of 32m×32m, and the explosion point is located in the center of the test area. 30 sensors are used. According to the total number of effective rays According to the size of the detection area and the ray density, the test area is divided into 64 × 64 grids of the same size, and the sensor array elements are distributed on the nodes on the edge of the grid and on individual nodes in the grid, as shown in Figure 4 .

2) Explosion point positioning

The reconstruction of the shock wave overpressure space-time field is carried out in a certain area centered on the explosion point, and the location of the explosion point will change at the moment of the explosion, so the precise location of the explosion point must be carried out first.

The sensor array elements arranged in space are used to receive target information in real time. Because the positions of the sensor array elements are different, the paths that the shock wave passes when reaching different sensor array elements are also different. There is a certain time difference between the sensor array elements receiving signals, that is, delay . According to the signal delay obtained by laying out the sensors in Figure 4, and the relationship between the sensor signal delay and the signal speed, the distance difference can be obtained. Combined with the sensor position coordinates, each sensor can establish an equation, and then use an iterative algorithm to calculate the explosion point. Spatial location.

3) Reconstruction of shock wave peak overpressure field

The reconstruction of shock wave peak overpressure field consists of two parts: velocity field reconstruction and peak overpressure field reconstruction. By measuring the relationship between the shock wave arrival time of the sensor array element and the ray path, the shock wave velocity of each grid unit can be reversed. According to the relationship between the shock wave velocity and the peak overpressure, the distribution of the shock wave peak overpressure field can be obtained.

(1) Velocity field reconstruction

Since the shock wave propagates in the air, the propagation medium can be regarded as a homogeneous medium, and the shock wave propagation speed is relatively large, so it can be approximately considered that the shock wave propagates along a straight line (wave front normal family). It is assumed that the shock wave does not propagate along the boundary of the grid unit, and the shock wave is approximately considered to propagate along a certain straight line in each grid unit, and each sensor corresponds to a ray. The time when the shock wave arrives at the sensor element is the travel time. During the transmission of the shock wave, its travel time is a function of velocity and geometric path, as shown in formula (1):

$\mathrm{<span\; class="notranslate">\; t</span>}<span\; class="notranslate">\; =</span>{<span\; class="notranslate">\; \&Integral;</span>}_{\mathrm{<span\; class="notranslate">\; L</span>}}\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; v</span>}}<span\; class="notranslate">\; \&Center\; Dot;</span>\mathrm{<span\; class="notranslate">\; dr</span>}<span\; class="notranslate">\; =</span>{<span\; class="notranslate">\; \&Integral;</span>}_{\mathrm{<span\; class="notranslate">\; L</span>}}\mathrm{<span\; class="notranslate">\; the\; s</span>}<span\; class="notranslate">\; \&Center\; Dot;</span>\mathrm{<span\; class="notranslate">\; dr</span>}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>$

In the formula, r is the ray; v is the velocity; s is the slowness; L is the integration path.

Discretize the above formula, according to the grid division of the test area, for the i-th ray:

$t_i=\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{\mathrm{the\; s}}_j{d}_{\mathrm{ij}}$ (i=1,2,3,…,M;j=1,2,3,…,N) (2)

In the formula, t _i : travel time of the i-th ray, that is, the time when the shock wave front arrives at the sensor;

d _ij : the ray length of the i-th ray passing through the j-th grid;

s _j : slowness in the jth grid;

M: number of rays;

N: Number of grids.

Write the above formula in matrix form as:

DS=T (3)

Among them, T=(t ₁ , t ₂ ... t _M )' is the M-dimensional column vector of the travel time of each ray; S=(s ₁ , s ₂ ...s _N )' is the slowness value of the discrete unit to be calculated, is an N-dimensional unknown column vector; D is an M×N order sparse matrix, and its elements are d _ij .

In the matrix equation (3), because the ray only passes through a small part of the grid cells, the matrix D is a large sparse matrix, and the non-zero elements account for a small proportion; the incomplete projection determines the underdetermination of the problem (non-unique nature); the existence of measurement error determines the incompatibility of the problem. In addition, due to noise and inaccurate measurement, it is very difficult to solve the exact solution of equation (3). Therefore, the equation is not directly solved here, but the difference between the observed travel time and the theoretical travel time is calculated first, and then the iterative method is used to solve .

The velocity field reconstruction process is shown in Figure 5.

a) Extract the measured travel time of the shock wave (that is, the time when the shock wave reaches the sensor element) data tc

By analyzing the signals of each sensor array element and extracting features, the time when the shock wave reaches the sensor array element is obtained.

b) Optimizing the initial slowness model

The initial slowness model is optimized according to the shock wave peak overpressure empirical formula:

When a spherical TNT charge with a charge density of 1.6 g/ ^cm3 explodes in an infinite medium, the peak overpressure p _m of the air shock wave can be calculated according to the following empirical formula:

$p_m=0.84\frac{\sqrt[3]{W}}{r}+2.7{\left(\frac{\sqrt[3]{W}}{r}\right)}^2+7{\left(\frac{\sqrt[3]{W}}{r}\right)}^3$ (for $1\&le;\stackrel{\&OverBar;}{r}\&le;10~15$ ) (4)

In the formula, W: charge weight (kg);

r: the distance from the measuring point to the center of the drug package (meters);

Scale distance (m/kg ^1/3 );

$p_m=\frac{20.06}{\stackrel{\&OverBar;}{r}}+\frac{1.94}{{\stackrel{\&OverBar;}{r}}^2}-\frac{0.04}{{\stackrel{\&OverBar;}{r}}^3}$ (for $0.05\&le;\stackrel{\&OverBar;}{r}\&le;0.50$ ) (5)

$p_m=\frac{0.67}{\stackrel{\&OverBar;}{r}}+\frac{3.01}{{\stackrel{\&OverBar;}{r}}^2}+\frac{4.31}{{\stackrel{\&OverBar;}{r}}^3}$ (for $0.50\&le;\stackrel{\&OverBar;}{r}\&le;70.9$ ) (6)

When other explosives are used or the charge density is different, the charge can be converted into an equivalent TNT equivalent according to the energy similarity principle and then calculated according to formulas (4) to (6).

The equivalent TNT equivalent calculation formula is as follows:

${\mathrm{<span\; class="notranslate">\; W</span>}}_{\mathrm{<span\; class="notranslate">\; T</span>}}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; W</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}\frac{{\mathrm{<span\; class="notranslate">\; Q</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}}{{\mathrm{<span\; class="notranslate">\; Q</span>}}_{\mathrm{<span\; class="notranslate">\; T</span>}}}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 7</span>}<span\; class="notranslate">\; )</span>$

In the formula, W _i : the weight of the explosive used (kg);

Q _i : heat of detonation of the explosive used (kcal/kg);

W _T : TNT equivalent converted from W _i (kg);

Q _T : Detonation heat of TNT (kcal/kg);

The relationship between shock wave velocity and peak overpressure is:

${\mathrm{<span\; class="notranslate">\; c</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; C</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; (</span>\frac{\mathrm{<span\; class="notranslate">\; 6</span>}{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; m</span>}}}{\mathrm{<span\; class="notranslate">\; 7</span>}{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 8</span>}<span\; class="notranslate">\; )</span>$

In the formula, p _m is the peak overpressure of the shock wave, c is the shock wave velocity, p ₀ is the initial pressure of the undisturbed air, C ₀ is the sound velocity of the undisturbed air, and for different temperatures: (m/s), T ₀ is the initial temperature of the undisturbed air (K).

The initial value of shock wave peak overpressure in each grid in the test area is calculated according to the above empirical formula, and then converted into the initial value of slowness according to formula (8).

c) Calculate the ray path d _ij

Because the path of the straight ray is the same in each iteration process of modifying the ray path, only the travel time caused by the change of wave velocity is different, so it only needs to be calculated once in the actual iteration process.

Take the vertex of the first grid unit at the bottom left as the origin O(0,0) of the coordinate system, set the coordinates of the explosion point as K(x ₀ ,y ₀ ), and the coordinates of the sensor corresponding to the i-th ray are s _i (x _i , y _i ), let the length and height of the rectangular test area be Length and Height respectively, and the grid is divided into H rows and L columns, then the length and height of each grid unit are Length/L, Height/H, i The ray equation is:

$\frac{\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}}{{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; the\; y</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}}<span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}}{{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 9</span>}<span\; class="notranslate">\; )</span>$

The horizontal dividing line equation is: y=mHeight/H (m is an integer, and 0≤m≤H),

The equation of the longitudinal dividing line is: x=nLength/L (n is an integer, and 0≤n≤L),

Simultaneously combining the ray equation and the transverse dividing line equation, and the ray equation and the longitudinal dividing line equation respectively, the intersection point of the straight ray and the transverse dividing line and the longitudinal dividing line can be obtained. For the i-th ray, according to the abscissa value of the intersection point on the ray path from the explosion point to the sensor s _i gradually decreases to sort all the intersection points.

Determine which grid unit the ray path passes through according to the horizontal and vertical coordinate values of two adjacent intersection points: Let the distribution of any two adjacent intersection points be (x _d , y _d ) and (x _d+1 , y _d+1 ) , then the coordinates of its midpoint are: $b_1=\frac{{\mathrm{the\; y}}_d+{\mathrm{the\; y}}_{d+1}}{2},$ Pick $a_2=\mathrm{INT}\left(\frac{a_1}{\mathrm{Length}/L}\right),$ $b_2=\mathrm{INT}\left(\frac{b_1}{\mathrm{Height}/h}\right),$ Then there is a grid unit serial number G=b ₂ ·L+a ₂ +1. Among them, INT represents rounding operation. Then calculate the straight-line distance between two adjacent intersection points in turn to obtain the path length d _ij .

d) Set the maximum number of iterations Maxq, and set the number of iterations to 100 in this calculation.

e) Calculate the theoretical travel time of each ray in turn according to the following formula

$t_{\mathrm{il}}^{\left(q\right)}=\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{d}_{\mathrm{ij}}{\mathrm{the\; s}}_j^{(q-1)}$ (i=1,2,3,...,M;j=1,2,3,...,N,q=1,2,3,...,Maxq) (10)

In the formula, q is the number of iterations (q=1,2,3,...,Maxq), is the theoretical travel time of the i-th ray at the qth iteration, (i=1,2,3,...,M; j=1,2,3,...,N); The resulting slowness value computed for the q-1th iteration. When q=1, the slowness value According to the initial model selection in b).

f) Calculate the difference between the measured travel time and the theoretical travel time of each ray in turn according to the following formula.

$\mathrm{\&delta;}{t}_i^{\left(q\right)}=t_{\mathrm{ic}}-t_{\mathrm{il}}^{\left(q\right)}$ (i=1,2,3,...,M,q=1,2,3,...,Maxq) (11)

g) The joint iterative reconstruction algorithm (SIRT) is used to sequentially calculate the average correction value of the slowness in each grid unit

Assuming that a total of K _j rays pass through in the jth grid cell, then:

$\mathrm{\&delta;}{\mathrm{the\; s}}_j^{\left(q\right)}=\frac{1}{K_j}\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}\mathrm{\&delta;}{t}_i^{\left(q\right)}\&Center\; Dot;d_{\mathrm{ij}}/\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{d}_{\mathrm{ij}}^2$ (i=1,2,3,...,M;j=1,2,3,...,N,q=1,2,3,...,Maxq) (12)

h) Correct the slowness value s _j of each grid cell

${\mathrm{the\; s}}_j^{\left(q\right)}={\mathrm{the\; s}}_j^{(q-1)}+\mathrm{\&delta;}{\mathrm{the\; s}}_j^{\left(q\right)}$ (j=1,2,3…N,q=1,2,3,…,Maxq) (13)

i) According to the measured velocity of the shock wave, the constraints are given.

${\mathrm{<span\; class="notranslate">\; the\; s</span>}}_{\mathrm{<span\; class="notranslate">\; min</span>}}<span\; class="notranslate">\; \&le;</span>{\mathrm{<span\; class="notranslate">\; the\; s</span>}}_{\mathrm{<span\; class="notranslate">\; j</span>}}^{<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; q</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; \&le;</span>{\mathrm{<span\; class="notranslate">\; the\; s</span>}}_{\mathrm{<span\; class="notranslate">\; max</span>}}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 14</span>}<span\; class="notranslate">\; )</span>$

In the formula, the values of s _min and s _max are determined according to the speed value obtained from the test, if but if ${\mathrm{the\; s}}_j^{\left(q\right)}<{\mathrm{the\; s}}_{\min },$ but ${\mathrm{the\; s}}_j^{\left(q\right)}={\mathrm{the\; s}}_{\min }.$

j) Set e=0.00001m/s, stop iteration when q≤Maxq satisfies formula (15), enter next round of iteration when q≤Maxq does not satisfy formula (15), repeat steps e) to i). When q>Maxq, but still does not satisfy formula (15), then stop the iteration.

$|{\mathrm{the\; s}}_j^{\left(q\right)}-{\mathrm{the\; s}}_j^{(q-1)}|\&le;e,$ (j=1,2,3…N,q=1,2,3,…,Maxq) (15)

Among them, c) and e) are the forward modeling process. f)-h) are slowness inversion process.

(2) Reconstruction of peak overpressure field

By comparing and fitting the shock wave velocity measured by the test with the peak overpressure value obtained by the sensor, the relationship between the shock wave peak overpressure at the test point and the shock wave velocity is obtained:

${\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; m</span>}}<span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; 7</span>}{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}}{\mathrm{<span\; class="notranslate">\; 6</span>}}<span\; class="notranslate">\; (</span>\frac{{\mathrm{<span\; class="notranslate">\; c</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}{{\mathrm{<span\; class="notranslate">\; C</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 16</span>}<span\; class="notranslate">\; )</span>$

where p _m is the shock wave peak overpressure, c is the shock wave velocity, p ₀ is the initial pressure of undisturbed air; C ₀ is the sound velocity of undisturbed air, for different temperatures: (m/s) ( _T0 is the initial temperature of the undisturbed air (K)).

Convert the reconstruction result of the slowness field into the velocity field, and then transform it into the reconstruction result p _m (x, y) of the shock wave peak overpressure field according to formula (16). The pressure value describes the spatial field distribution of the shock wave overpressure.

4) Shock wave overpressure space-time field reconstruction

According to the time-varying curve of shock wave overpressure at a certain point obtained by the test sensor during the test, an appropriate shock wave overpressure-time transmission model is selected, and through fitting and comparison, combined with the above shock wave peak overpressure field reconstruction results p _m (x,y ), choose the "modified Friedlander equation" to get the relationship of shock wave overpressure with time and space:

$\mathrm{<span\; class="notranslate">\; p</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; t</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; m</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; t</span>}<span\; class="notranslate">\; /</span>{\mathrm{<span\; class="notranslate">\; T</span>}}^{<span\; class="notranslate">\; +</span>}<span\; class="notranslate">\; )</span>{\mathrm{<span\; class="notranslate">\; e</span>}}^{<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; ct</span>}<span\; class="notranslate">\; /</span>{\mathrm{<span\; class="notranslate">\; T</span>}}^{<span\; class="notranslate">\; +</span>}}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 17</span>}<span\; class="notranslate">\; )</span>$

In the formula, p(x, y, t) is the overpressure value of the shock wave at the spatial position (x, y) at time t, p ₀ is the ambient pressure; p _m (x, y) is the Peak overpressure; T ⁺ is positive phase duration; c is attenuation parameter.

(1) p ₀ can be measured before the experiment;

(2) p _m (x, y) has been calculated and obtained through the reconstruction of the shock wave peak overpressure field;

(3) The following relationship model can be established between T ⁺ and the test charge W, and the distance R between the test sensor and the explosion point:

$\frac{{\mathrm{<span\; class="notranslate">\; T</span>}}^{<span\; class="notranslate">\; +</span>}}{\sqrt[\mathrm{<span\; class="notranslate">\; 3</span>}]{\mathrm{<span\; class="notranslate">\; W</span>}}}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; a</span>}{<span\; class="notranslate">\; (</span>\frac{\mathrm{<span\; class="notranslate">\; R</span>}}{\sqrt[\mathrm{<span\; class="notranslate">\; 3</span>}]{\mathrm{<span\; class="notranslate">\; W</span>}}}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; no</span>}}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 18</span>}<span\; class="notranslate">\; )</span>$

The above formula can be transformed into:

n(lnR-lnB)+lna=lnT ⁺ -lnB (19)

in Treat n and lna as unknown quantities x, y,

Then this equation can be transformed into a linear equation

z ₁ x+y=z ₂ (20)

Where z ₁ =lnR-lnB, z ² =lnT ⁺ -lnB.

The distance from the i-th test sensor to the explosion point is R _i (i=1,2,3,...,M), and the charge is W, then z _1i ＝lnR _i -lnB(i=1,2,3, ..., M); put the positive phase duration T _i ⁺ of each sensor test shock wave curve obtained during the test into z ₂ , and z _2i = lnT _i ⁺ -lnB(i=1,2,3,...,M) can be obtained , establish a linear equation system z _1i x+y=z _2i (i=1,2,3,…,M), and use the least squares method to fit it, you can get the parameters x, y, that is, get the parameters n, a. Put it into formula (18) to get the functional relationship between the parameter T ⁺ and the distance R and the dose W; determine the distance R between the test point and the explosion point according to the spatial position (x, y) in formula (17), and combine the dose W to obtain The positive phase duration T ⁺ in formula (17) can be determined.

(4) The parameter c is determined

The parameter c is not a constant, but exhibits complex properties. Through the analysis of the measured signal, it is found that the parameter c is not only related to the peak overpressure, but is actually a multivariate function, which is related to both the time t and the distance r between the explosion source and the measuring point, using the formula (17) The expression of the parameter c can be deduced backwards:

$\mathrm{<span\; class="notranslate">\; c</span>}<span\; class="notranslate">\; =</span><span\; class="notranslate">\; (</span>\frac{{\mathrm{<span\; class="notranslate">\; T</span>}}^{<span\; class="notranslate">\; +</span>}}{\mathrm{<span\; class="notranslate">\; t</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; [</span>\mathrm{<span\; class="notranslate">\; ln</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; m</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; T</span>}}^{<span\; class="notranslate">\; +</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; t</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; )</span><span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; ln</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; T</span>}}^{<span\; class="notranslate">\; +</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; p</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; the\; y</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; t</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; p</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; )</span><span\; class="notranslate">\; ]</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; twenty\; one</span>}<span\; class="notranslate">\; )</span>$

In the formula, p(x, y, t) is the overpressure value corresponding to different time t and different spatial positions (x, y) on the shock wave curve measured by the test sensor, which is obtained from the output data of the test sensor, and p(x, y, t ) and the corresponding T ⁺ into equation (21) can get the time history curve of parameter c.

Through the research and analysis of the time history curve of the parameter c of the measured signal, it is found that the logarithm of the logarithm of c (lnlnc) has a linear relationship with time, that is:

kt+m=lnlnc(0≤t≤T ⁺ ) (22)

Among them, k and m are undetermined parameters.

set matrix $A={\left[\begin{array}{cc}1& 1\\ 2& 1\\ 3& 1\\ ...& ...\\ \mathrm{data}& 1\end{array}\right]}_{\mathrm{data}\&times;2},$ $x=\left(\begin{array}{c}k\\ m\end{array}\right),$ $B={\left[\begin{array}{c}\ln \ln c_1\\ \ln \ln c_2\\ \ln \ln c_3\\ ...\\ \ln \ln c_{\mathrm{data}}\end{array}\right]}_{\mathrm{data}\&times;1}$

Among them, the first column of matrix A indicates the different time t corresponding to the sampling points (take data sampling points) on the shock wave curve measured by the test sensor, and matrix B indicates the lnlnc value at the corresponding moment of the sampling point, and the logarithm is obtained by formula (21) The logarithm is obtained. Equation (22) is a linear equation about k and m, which is brought into the measured signal data to establish a linear equation system as follows:

AX=B (23)

After transformation (23) becomes:

A ^{T A} X = A ^T B (24)

Where ^AT is the transpose matrix of matrix A. but

X＝(A ^T A) ^-1 A ^T B (25)

The parameters k and m can be determined by solving equation (25), and the parameter c can be determined by equation (22).

5) Visualization of shock wave overpressure space-time field

Transform the overpressure reconstruction data into geometric element data such as point, line, surface, body, etc.; convert the geometric data into image data, and transmit it to the display device to generate an image file; smooth the generated image, and finally The shock wave produced by the explosion is dynamically displayed with the change process of time and space position.

6) Reconstruction results

The patented method is used to carry out two-dimensional reconstruction of the explosion shock wave field within the range of 32m×32m, and the reconstruction result of the peak overpressure when the shock wave arrives is compared with the test result of the overpressure sensor during the test. The minimum difference is 0.000961Mp, and the two-dimensional explosion The results of 3D visualization of field overpressure reconstruction results are shown in Fig. 7.

Claims (3)

1. A reconstruction method for explosion shock wave overpressure space-time field, adopting shock wave test system to realize the reconstruction of explosion shock wave overpressure space-time field, shock wave test system is composed of shock wave overpressure sensor array, signal conditioning module, synchronous A/D conversion module, high-speed data It consists of a transmission module, a data preprocessing module, a shock wave overpressure space-time field reconstruction module, and a computer. Among them, the shock wave overpressure sensor array converts the shock wave overpressure signal into an electrical signal, which is conditioned into a standard by the signal conditioning module and the synchronous A/D conversion module. The voltage signal is converted into a digital signal, and then the digital signal is transmitted to the computer through a high-speed data transmission module. The data preprocessing module and the overpressure space-time field reconstruction module are installed on the computer, and the shock wave overpressure signal after A/D conversion is input The data preprocessing module is input into the shock wave overpressure space-time field reconstruction module after data preprocessing, and the shock wave overpressure space-time field reconstruction module completes the reconstruction of the explosion shock wave overpressure space-time field; it is characterized in that the space-time field reconstruction method is: first The test area is divided into grids, and the sensor array elements are arranged; secondly, the time difference method in the passive positioning algorithm is combined with the iterative algorithm to locate the explosion point, and then based on the relationship between the travel time and the ray path and the relationship between the shock wave velocity and the peak overpressure. Shock wave peak overpressure field reconstruction; then study the shock wave overpressure attenuation law and characteristics with time, combined with the grid unit shock wave peak overpressure, calculate the parameters of the "modified Friedlander equation", and obtain the shock wave overpressure space-time data; at last, transform the reconstruction data into an image file and display it on the monitor of the computer to realize the visualization of the shock wave overpressure space-time field; wherein: the reconstruction of the shock wave peak overpressure field consists of velocity field reconstruction and peak overpressure field reconstruction. Partial composition. the 2. A kind of explosion shock wave overpressure space-time field reconstruction method according to claim 1, is characterized in that: the concrete method of described shock wave peak overpressure field reconstruction: extract shock wave measured travel time, optimize slowness initial model, calculate ray path, set the maximum number of iterations, calculate the theoretical travel time, calculate the slowness field inversion, and convert the peak shock wave overpressure to obtain the reconstruction result p _m (x, y) of the shock wave peak overpressure field. 3. a kind of explosion shock wave overpressure space-time field reconstruction method according to claim 1, is characterized in that: the concrete steps of this space-time field reconstruction method are: (1) During the grid division of the test area and the layout of the sensor array elements, the measurement area is divided into several regular grid units, and the sensor array elements are distributed on the nodes on the edge of the grid and on individual nodes in the grid according to the number of test use. The explosion point is in the center of the test area; (2) According to the relationship between the signal time delay and the signal speed, the distance difference is obtained, combined with the sensor position coordinates, an equation is established for each sensor, and then the iterative algorithm is used to calculate the spatial position of the explosion point; by measuring the shock wave reaching the sensor array element The relationship between the time and the ray path is reversed to obtain the shock wave velocity of each grid unit. According to the relationship between the shock wave velocity and the peak overpressure, the reconstruction result of the shock wave peak overpressure field is obtained, that is, the shock wave peak overpressure p _m (x,y ); (3) According to the time-varying curve of the shock wave overpressure at a fixed point obtained by the test sensor during the test, select the appropriate shock wave overpressure-time transmission model, and select the "modified Friedlander equation" by fitting comparison Combined with the reconstruction results p _m (x, y) of shock wave peak overpressure field, the relationship of shock wave overpressure with time and space is obtained, and the reconstruction result of shock wave overpressure space-time field is obtained, that is, the shock wave overpressure value p at any position and time (x,y,t); (4) Establish the mapping between overpressure reconstruction data and color or geometric shape, and complete the transformation of application data to point, line, surface, and volume geometric element data; change geometric data into image data, and generate image files from image data, Then, the generated image is smoothed to eliminate individual abnormal points that do not conform to the law of shock wave transmission, and finally the shock wave overpressure generated by the explosion is dynamically displayed with the change process of time and space position. the