Sensor Distribution Design of Travel Time Tomography in Explosion

Optimal sensor distribution in explosion testing is important in saving test costs and improving experiment efficiency. Aiming at travel time tomography in an explosion, an optimizing method in sensor distribution is proposed to improve the inversion stability. The influence factors of inversion stability are analyzed and the evaluating function on optimizing sensor distribution is proposed. This paper presents a sub-region and multi-scale cell partition method, according to the characteristics of a shock wave in an explosion. An adaptive escaping particle swarm optimization algorithm is employed to achieve the optimal sensor distribution. The experimental results demonstrate that optimal sensor distribution has improved both indexes and inversion stability.


Introduction
Explosion technology is applied increasingly [1,2]. As explosion experiments are complex and the costs are high, optimizing sensor distribution is necessary to save test costs and improve experiment performance. The number of sensors to be used and their locations should be optimized [3][4][5].
When a shock wave propagates in air, the shock wave velocity is very fast. Supposing that a shock wave propagates by direct rays without propagation by grid boundary, in the course of shock wave propagation, the relationship of travel time is: (1) The wave front normal is defined as rays, i.e., r. Each sensor corresponds to a ray. v is velocity and s is slowness. The travel time t is the integral of slowness along with rays. The test region is divided into some regular or irregular cells as shown in Figure 1, the discretization expression of Equation (1) (2) where t i is the travel time that a shock wave travels from the explosive to the ith sensor. d ij is the length of the ith ray in the jth cell. S j is the slowness in jth cell. M is the sensor number and N is the cell number. Equation (2) can be written in a matrix form as: where 1 2 ( , , , )'  is the m-dimension column vector of travel time; 1 2 ( , , , )' N S s s s   is an unknown N dimensional column vector; D is the distance matrix of N M  and its element is d ij [6,7]. The elements of T can be obtained by the tested data from sensors. The matrix D can be calculated by the position of the explosive and sensors. The unknown slowness S can be obtained by solving Equation (3). Based on the above-mentioned principle, the velocity distribution of shock wave can be obtained. For the travel time tomography, the driving source is single and the sensors are few. The tomography rays are distributed sparsely. The Equation (3) is ill-posed and underdetermined [8]. For this tomography modality, optimizing sensor distribution is necessary for reducing costs and increasing acquired information. Optimizing sensor distribution is to solve ill-posed tomography problems from the perspective of mathematics and then the inversion stability can be ensured.
The ill-posed degree of Equation (3) is related to the structure of matrix D. Improving the ill-posed degree of Equation (3) is to improve the structure of matrix D. The structure of matrix D rests with a parameterized model of the test region and sensor distribution [9].

Effect of Eigenvalue and Rank on Inversion
In Equation (3), with a given data vector 0 t , the model vector 0 s ∈S and then minimized. This is accomplished by pre-multiplying Equation (3) by T D and taking a matrix inverse: is often near-singular, it results in instability in the solution. That is, some of its eigenvectors have extremely small eigenvalues   : 1, , Measurement errors in data space T propagate into the solution S in parallel to each eigenvector with an amplification i  / 1 . Hence, when small eigenvalues exist, the solution becomes unstable and the inverse problem is ill-posed [10]. When an eigenvalue is zero, the corresponding eigenvector in the data space can not be mapped into the model space [11]. The greater the rank is and the larger the eigenvalues are, the more stable the inversion problem is and the more independent pieces of information may be gained from the data. Hence, optimal sensor distribution can be obtained by maximizing the magnitude of eigenvalues of L and the rank of D. The evaluating function is: The evaluating function E 1 has two components, one represents the relative distribution of eigenvalues and the other indicates the relative size of null space. If the value of E 1 is minimal, the inversion results are optimal and stable.

Effect of Condition Number on Inversion
The Equation (3) is ill-posed if small perturbations of matrix D or T can result in a large change in solutions [12].
Assuming that the observed data T has a minor perturbation of T  , the perturbation of solution is S  . Equation (3) becomes: Then, According to the property of subordinate norm, is the condition number of matrix D. When the observed data T has a minor perturbation, the relative error of solution is determined by the condition number. Supposing that T is accurate and the matrix D has a minor perturbation of D  , the corresponding is the solution of perturbation equation, Similarly the following condition can be obtained: When the matrix D has a minor perturbation, the error of solution is also determined by the condition number.
Therefore, the condition number is the second judging index of matrix D's quality. The smaller condition number means the more well-posed equation.
The evaluating function expressed by the condition number can be written as: Optimal sensor distribution can reduce the condition number and the more stable inversion solutions can be obtained.

Effect of Ray Coverage on Inversion
When designing sensor location, rays coverage should be enlarged as much as possible, i.e., rays density and orthogonality should be maximized. The ray density represents the number of rays passing through each cell. The cells not being hit by any ray are the main factors giving rise to the null space. They make some column vectors of matrix D be zero (d j = 0) so that the equation can not be resolved. Therefore, tomography with sparse rays must ensure that any column vector is non-zero. Increasing the ray density can avoid zero vectors.
Ray orthogonality is measured by maximal sinusoidal quantity of angle between rays [13]. The small orthogonality makes some rows in D linearly dependent.
The greater the ray density is, the better the orthogonality is and the smaller inversion error can be achieved. The evaluating function expressed by rays density and orthogonality can be written as: where j  is the ray density in the jth cell. j O is the ray orthogonality in the jth cell.  There are many factors on the ill-posed and underdetermined equations except for the above analysis. However, the value of E can be used as a judging index with regard to optimizing sensor distribution and the experiment. The detailed process is as follows and a flow chart is illustrated in Figure 2. (1) Dividing cells according to model character.
(2) Giving a distribution model randomly according to the number of sensors and calculating matrix D and the rank of D.
If all column vectors of matrix D are non-zero and D is full rank, make this distribution model as an initial model; otherwise, go to Step (2). (4) When the number of initial model is equal to the given number K , the optimization algorithm is employed to obtain the optimum value of E and sensor distribution.

Sub-Region and Multi-Scale Cells Partition
The test region is divided into cells and each cell has the same velocity. The dividing pattern of cells affects matrix D. We propose a sub-region and multi-scale cell partition method. Cells are divided in terms of the solution distribution, i.e., the smaller size and the higher density of cells correspond to the bigger changing gradient of solutions, and the bigger size and the lower density correspond to the smaller changing of solutions.
When the explosive is exploding in air, the shock wave overpressure attenuates quickly in the near region to the explosive and the attenuation becomes slow with the increase in distance. According to the shock wave characteristic, the smaller cells are adopted in the near region to the explosive while the bigger cells are adopted in the far region. In order to avoid a row vector to be linearity dependent on matrix D aroused by symmetry, cells are divided into different size in the symmetrical region with consideration of the resolution of inversion.

Particle Swarm Optimization (PSO) Algorithm and Modification
Particle swarm optimization (PSO) algorithm is simple and easy to implement. However, PSO can fall into the local optimum [14,15]. The original algorithm is modified and an adaptive escaping particle swarm optimization algorithm (AEPSO) is proposed.
When P g does not change over M generations, all particles are close to P g . The flying speed of particles is very little and subsistence density is large. An escape strategy should be adopted to update the velocity of particles and enlarge the searching space. The escape strategy changes the velocity of particles and makes variation in order to increase the diversity of particles and jump out of local optimization. When escaping, the particles are divided into two components, one component updates their velocity according to Equation (16) and the other component updates their velocity according to Equation (17).
where 3 r and 4 r are the random numbers between (0, 1). A is a constant that controls the particle velocity. max v is the maximal velocity value that particles have.
The declining linearly inertia weight helps to swarm searching. The adjusting strategy is as follows: is the current inertia weight; max w and min w are the maximum and minimum inertia weights; max q is the maximum number of iterations; k is current iteration times.

Optimizing Sensor Distribution Based on AEPSO Algorithm
The adaptive escaping particle swarm optimization algorithm is adopted to obtain the optimum value of E and the sensor distribution. The detailed process is as follows and a flow chart is illustrated in Figure 3.
(1) Producing a particle. A d-dimension particle is produced by selecting a distribution model randomly according to the number of sensors, which is expressed as (2) Calculating matrix D and the rank of D.
(3) If all column vectors of matrix D are not zero and matrix D is full rank, put this particle into the initial particle swarm; otherwise, return to Step (1). (4) Calculating the optimal group fitness when the number of particles in initial swarm is equal to the given number. (5) Updating the velocity and position of each particle according to Equations (14) and (15) and calculating fitness. (6) Updating the individual optimal fitness. If the current fitness is superior to experienced optimal fitness i P , the experienced optimal fitness i P is replaced by the current fitness.
(7) Updating the group optimal fitness. If the current fitness is superior to the group optimal fitness, the group optimal fitness g P is replaced by current fitness.
(8) Recording group optimal fitness g P in each iteration. If the value of g P does not change over M generations, return to Step (9) and adopt escaping strategy; otherwise, return to Step (5). (9) Dividing the particles into two components and updating particles velocity according to Equations (16) and (17). (10) If the ending condition is met, terminate the iteration or return to Step (5).

Comparison of Cell Performance
To verify the sensor distribution design, a test region is 16 m × 16 m, explosive is placed in the lower-left and sensors are distributed on the region boundary. Two cell patterns are used in dividing the test region. The first one is a uniform rectangle with 49 cells as shown in Figure 4a. The second one is sub-region and multi-scale cells as shown in Figure 4b. The velocity of shock wave decreases exponentially with the distance according to the prior information. In the near region to the explosive, the attenuation amplitude of velocity is larger and the cells are smaller and denser. In the far region the attenuation amplitude of velocity is smaller and the cells are bigger and sparser. The region is divided into seven sub-regions according to the proportional distance and different size cells in the symmetry region are placed in order to avoid a row vector to be linearity dependent on matrix D . The total cells number is 58. The influence of cells on matrix D is compared with the same sensor distribution.
(1) To ensure each cell being passed through by at least one ray, the lowest number of sensor required by uniform rectangle cells is 9 as shown in Figure 4a. The multi-scale cells require at least five sensors as shown in Figure 4b. (2) Each index of matrix D in uniform rectangle cells and multi-scale cells is compared with 13 sensors and the same distribution. Simulation results are given in Table 1. It is clear that the indexes in multi-scale cells are superior to those of the uniform rectangle cells. The matrix D in multi-scale cells is full rank mostly and the condition number is far smaller than that of uniform rectangle cells. The ray density and orthogonality in multi-scale cells are larger than that of uniform rectangle cells generally.

Parameters Setting and Simulation of AEPSO Algorithm
To validate the AEPSO algorithm, two multi-modal functions are used to test the PSO and AEPSO algorithms.
Rastrigrin function: Griewank function: The two functions have many local minimum values, which has the global minimum value of 0 when x i = 0. The mean best fitness (MBF) and standard deviation (SD) are used to estimate the solutions. Parameters setting are given in Table 2. The simulation results are given in Table 3. From the table, it is concluded that AEPSO is superior to PSO algorithm in accuracy, convergence, and stability.

Optimizing Sensor Distribution
Each index based on multi-scale cells with symmetrical distribution and random distribution of sensors is shown in Figure 5a,b. According to the values of 1 E , 2 E ,  and O , the weight coefficients are selected as: . The adaptive escaping particle swarm optimization (PSO) algorithm is adopted to obtain the optimum value of E and sensor distribution as shown in Figure 5c. The values of 1 E , 2 E ,  and O with different sensor distributions are in Table 4. It is clear that the indexes in optimum distribution are superior to that of the others: smallest 1 E and 2 E , and the biggest 3 E . Figure 6 shows a comparison of the normalized eigenvalue spectra respectively in each case. The optimized model is again superior to that of the others. Numerical simulations are presented with the same initial model and inversion algorithm. From the results it can be concluded that the optimum distribution can result in the smallest inversion error.

Conclusions
For travel time tomography in explosion, optimizing sensor distribution is very important for saving the costs and increasing the acquired information. This paper shows the inversion stability can be improved by designing optimal sensor distribution. This paper analyzed the effect of eigenvalue, rank, condition number and rays coverage on improving matrix D and proposes the evaluating function on optimizing sensor distribution. An adaptive escaping particle swarm optimization algorithm is used to obtain the optimum sensor distribution based on sub-region and multi-scale cells. Simulation shows that the optimization method is feasible and the optimal sensor design achieves inversion stability.