PROCEEDINGS of the
23rd International Congress on Acoustics
9 to 13 September 2019 in Aachen, Germany

Study on large scale projectile impact point positioning method
based on trajectory shock wave
Difeng SUN1; Xubin LIANG1; Tianqing ZHAO1; Houlin FANG1; Zhang CHENG1; Hui ZHENG1;
Liangyong ZHANG1; Fang ZHANG1; Deyu SUN1; Yang LIU1
1

Northwest Institute of Nuclear Technology, China

ABSTRACTWith multiple nodes of single acoustic sensor layed around the target, the ballistic shock wave signal of the
projectile generated before hitting the ground can be detected at different locations. The trajectory can be
estimated by measuring the time difference between the signal of different nodes, and the impact point can becoordinated by crossing the trajectory with the target. In this article, the projectile impact point positioning
model is parameterised and fomulated by nonlinear equations with six patameters, and it can be solved bymeans of nonlinear least square (NLS) method. By means of Cramer-Rao lower bound (CRLB) analysis,
measurement error distribution by different incident angle and different impact location is calculated, andalso, parameter of exterior ballistics model and wind vector is under consideration.
Keywords: Impact Point Positioning, Projectory Trajectory, CRLB

1. INTRODUCTION
When the projectile flies in the air at supersonic speed, its front end sharply compresses the air infront, forming shock wave(1, 2). As it diffuses, the shock wave gradually degenerates into weakcompression wave, and the propagation speed drops from supersonic to sound, finally forming thecone wavefront(3).
At the sensor, the shock wave is measured in the form of N -wave(4). The trajectory of the
projectile can be estimated according to the arrival time of the N -wave measured at different
positions. The intersection point between the trajectory and the ground or obstacles is the projectile
impact point(5). In early work, the ballastic equation of motion was studied(6). Ballistic const Cb
was introduced to describe the air resistance effect(7). If the gravitational acceleration g is ignored,
the ballistic trajectory will be linear. Researches show that with ballistic constant C b known and dragcoefficient exponent η=0.5, the trajectory parameter and impact point can be well estimated(7).
In this paper, the impact point is estimated in the case that sensor coordinates and sound speed are
not exactly measured or presumed. Also, the influence of wind on sound speed vector is coarsely
modeled. This paper focuses on the study on error analysis method based on Cramer-Rao lower
bound (CRLB)(8). Through numerical simulation and CRLB method, the error distribution rule ofimpact point for large scale target is obtained, and then other scale targets are also taken into
account.

2. PROBLEM MODELING
2.1 Basic assumptions
In order to simplify the model of supersonic flight projectile ballistic shock wave, someassumptions are made. 1) The near-field nonlinear region of projectile is ignored and the velocityof shock wave is considered as sound speed. 2) The curvilinear induced by gravity acceleration g is
ignored since its acting time is too short. 3) After the projectile hits the ground, its motion equation
changes and some sensors’ signal data is affected. But here, this process is omitted since it can besolved in other ways such as excluding the abnormal sensors.
2.2 Geometric Model
Geometric model settings are similar to that in literature(9), as shown in Figure 1(a). N sensors S n
13346

sundifeng@nint.ac.cn

are located in the target perimeter with position vector Rn =[Xn ,Y n,0] T.The projectile hit the ground ata certain elevation angle θ and azimuth angle ϕ, on point Q with vector RQ =[XQ ,Y Q,0] T. The unit
direction vector of trajectory is e v=-[cosθcosϕ, cosθsinϕ, sinθ] T. In the plane πn which determined bylinear trajectory and sensor Sn , a triangle is formed by the hitting point Q, the sensor point S n and thedetach point P n. The parameter set λ ^ X Q , YQ ,T , I ,Vr ` , is what we focus on, especially impact
point coordiantes X Q and YQ in it.
Z

x'n
1.24

Q
ev

ev

RQ
Q
α'n ©
Rn
ϕ

X

©n
sn

r'n

O

Pn

α'n

Y

Sn

SN

Sn

(a)
(b)
Figure 1 – (a) Geometric distribution of the sensors. The linear trajectary intersects the ground at Q(b)Ballistic shock wave reach sensor Sn
2.3 Mathematical model
In the ballistic model, g is ignored, the ballistic constant C b is known and drag coefficient
exponent η=0.5. The instantaneous velocity v '( x ') at x ' and the flight time t '( x ') from areference point to P n are calculated as in literature(9).
The time stamp of shock wave arriving at Sn is given by Equation (1) and sn is calculated withsine theory in triangle P nQS n.
W n t '( x 'n )  sn / c(1)
The Differencial Time of Arrival (DTOA)(9) between Sensors S m and S n (1 d m, n d N ) :W nm W n  W mt '( x 'n )  t '( x 'm )  ( sn  sm ) / c(2)

2.4 NLS Method
The estimated value of DTOA can be expressed as:
ˆ , cˆ, Cˆ )  n(3)bˆ ( τˆ ) b(λ; RbWhere bˆ [Wˆ21 ,Wˆ31 , ,WˆN 1 ]T , b [W 21 ,W 31, ,W N 1 ]T and n [n1 , n2 , nN 1 ]T are respectively theobservation vector, the model value vector and the total observation error vector of DTOA. The caret“^” on a certain variable represents an observation, presumption or estimation of the actual value orˆ [ Xˆ ,Yˆ , Xˆ ,Yˆ , Xˆ ,Yˆ ]T , ĉ , Cˆ are respectively thepresetting value. τˆ [Wˆ ,Wˆ , ,Wˆ ]T , R
1

2

N

1

1

2

2

N

N

b

observation or the presumption value (vector) of arrival time stamp vector τ [W1 ,W 2 ,coordination vector R [ X1 , Y1 , X 2 , Y2 ,observation or presumption error n τ,W N ]T , sensor

X N , YN ]T , sound speed c and ballistic constant Cb , with[nW 1 , nW 2 ,nW N ] , nR

[nX 1 , nY 1 , nX 2 , nY 2 ,nXN , nYN ], nc , nCb .

Besides, the observation error vector of DTOA caused by the sampling procedure can be described asnb [nτ 2  nτ1 , nτ 2  nτ1 , nτ N  nτ1 ]T . λ ^ X Q , YQ ,T , I ,Vr ` is the set of unknown parameter toresolve. Use a (N-1) demensional column vector F to describe the total error between observation
vector and model vector.
ˆ , cˆ, Cˆ ) bˆ ( τˆ )  b(λ; Rˆ , cˆ, Cˆ )(4)F(λ; τˆ , R
b
b

3347

Non-linear Least Square (NLS) method(9) can be used to estimate λ with Equation (11):ˆ , cˆ, Cˆ )T N 1F(λ '; τˆ , R
ˆ , cˆ, Cˆ )λˆ arg min F(λ '; τˆ , R
b
b
λ'

(5)

Since the observation error covariance matrix of DTOA is usually unknown, an assumed errorˆ V 2 (I  E) is alternative. In this formula, V 2 is the arrival time stamp
convariance matrix N
error caused by the sampling procedure, I is the (N-1)×(N-1) identity matrix and E is the unit
matrix. Equation (11) can be solved by numerical iterative algorithm.

3. ERROR ANALYSIS AND CRLB METHOD
3.1 Error Source Analysis
Traditional methods generally assume that the error of DTOA of two sensors is only caused by
the signal acquisition process. But in reality, it’s also indirectly influenced by various factors during
the flight.
The position error of sensors will affect the total distance of shock wave propagation in theair(10).
The different choice of trajectory model has a great influence on the impact point positioning
problem. The ballistic constant C b is one of the parameters of the model. In this paper, Cb is supposed
to be already known.
The shock wave speed is mainly affected by wind, if tempearture factor is ignored. The influenceof the wind should be the key to consider when the test is done outside (11), while most articles
choose to ignore its impact because it’s too complex.
3.2 Coarse One Dimensional Wind ModelingIt’s difficult to use a general method for wind modeling. Ignoring its influence on the trajectory of
projectile, wind mainly affects the velocity and direction of shock wave propagation relative to theearth. The wind is momentarily uniform for small area, but has strong inhomogeneity, instantaneity
and randomness for large area. As shown in Figure 2, the effect of wind speed of the same amplitude
on sound vector is equivalent to adding a circle whose center point is at the end of sound vector andradius is equal to wind speed amplitude.
c' c  w
(6)

r=|w|

Δθ
c
θ

c'

v

Figure 2 – Influence of wind on sound speed
The shock wave speed vector, namely the revised sound speed vector c ' is influenced by windspeed vector w . Since the amplitude of wind w is far less than sound speed c, the change of shock
wave angle Δθ is also far less than original angle θ. The maximum and minimum influence of wind
is obtained when revised angle T ' arcsin(c '/ v) meets the extreme value. Since w c , c ' and c
is approximately parallel and Equation (13) can be concluded.
c 'max | c  w
c 'min | c  w(7)
Between the maximum and minimun influence of wind, appropriate intersection of wind
distribution can be used to simply construct the wind model, in this way wind can be modeled in thesame (or opposite) direction of sound speed. A gaussion distribution of wind w is finally used inconvenience.
In this paper, the error of c from temperature is omitted, then we see that cˆ c ' and nc w , i.e.
the distribution of sound speed error is the same as the distribution of wind w. Although thisone-dimensional model leads to linear error distribution of impact point(as shown below in this3348

paper), which may not fit actual situation, the amplitude of error yet has reference value.
3.3 CRLB Error Analysis
Assume that the observation error vector n follows the zero-mean Gaussian distribution, withconvariance matrix N. When the observation (or presumption) of arrival time stamp τ̂ , sensor
coordinate R̂ , sound speed ĉ and ballistic constant Cˆ is known, the CRLB on conviance matrixb

is as shown(12, 13):

ˆ , cˆ, Cˆ ) [F T N 1 (F T )T ]1
Q(λ; τˆ , Rb(8)

Where  [w / wX Q , w / wYQ , w / wT , w / wI , w / wVr ]T is the gradient operator, and FT can beˆ , cˆ, Cˆ } and presetted value λ .
estimated by the observation or presumption value {τˆ , Rb
The total observation error vector n is determined by many aspects. Ignoring second and higher
order small terms, use first order Taylor expansion to describe the propagation law of error, we get:wFˆ , cˆ, Cˆ )  F( λ; τ, R, c, C )F (λ; τˆ , R¦ wT  (Tˆ  T )bb
T {τ , R , c ,Cb }
wF

wF

wF

wF

¦ wW  (Wˆ  W )  ¦ wX  ( Xˆ  X )  wc (cˆ  c)  wC (Cˆb  Cb )W {τ}

X {R}

(9)

b

wF ˆ
wF
wF ˆ
 (X  X ) (bˆ  b)  ¦(cˆ  c) 
(Cb  Cb )wc
wCb
X {R} wX
n

wF
wF
wF
 nX 
 nc 
 nCb
X
c
Cb
w
w
w
X {R}

nb  ¦(10){τ, R, c, Cb } is the set of all elements in τ, R, c, Cb , {τ} , {R} have the similar meaning. The
ˆ , cˆ, Cˆ and τ, R, c, C represents the total
function value of F of different paramaters τˆ , Rbbobservation error calculated by parameters with certain error distribution or not. F (Tˆ)
T {τ, R, c, Cb } represents the function value of all exact parameters except θ. F(λ ; τ , R ,c ,Cb ) 0
from definition, and F is the abbreviation of F(λ; τ, R, c, Cb ) . Suppose all n {n τ , nR , nc , nCb } areindependent identically distributed and follow zero-mean Gaussian distribution, we have:
Nb E(nb  nbT )E([nτ ,2  nτ ,1 , nτ ,2  nτ ,1 ,VNnτ ,NN  nτ ,11 ]T  [nτ ,22  nτ ,11 , nτ ,22  nτ ,11 ,nτ,NN  nτ,1 ])(11)

2
τ ( I  E)
T

E(n  n )

w
w
w
w
w
w
 nX   nc  nCb )  ( ¦ nX   nc  nCb )F T ]
X
c
C
Y
c
C
w
w
w
w
w
w
X {R}
Y {R}bb

N b  E[F ( ¦

T

T

(12)

T

wF wF
wF wF
wF wF
)  V c2 ()  V C2b ()
wc wc
wCb wCbX {R} wX wX

V τ2 (I  E)  V R2 ( ¦Nb is the DTOA error covariance matrix caused by sampling procedure, and N is the total DOTA

observation error covariance matrix. V τ , V R , V c , V Cb are the corresponding standard deviation. From
Equation (10) we can conclude that N follows indeed the zero-mean Gaussion distribution if all
relevant error distributions follow zero-mean Gaussion distribution. n and N are both made up offour terms and influenced by them.
In special cases, if the last three term are zero or can be negelected, then
(13)n nb bˆ ( τˆ )  b(λ; R, c, Cb ), N Nb V τ2 (I  E)
This is in accordance with traditional methods(7, 9, 11).

3349

4. SIMULATION RESULTS
4.1 Simulation Setting
Suppose the square target is on the ground plane with side length of 4 a (a is the sensor spacing
and represents the scale of target), and N=16 sound detection nodes of single sensor are uniformly
distributed around it, as shown in Figure 1(a). The position vectors of N sensors are R 1=[0,0,0],R2=[0,a,0], etc. The standand deviation of X, Y and Z-coordinates is 3cm. Suppose the ballistic
constant C b=150 is exactly known in advance. Considering winds within gentle breeze (0~5.4m/s),
the standard derivation of shock wave speed, namely the standard derivation of wind is set to 2.7m/s.
Suppose the arrival time stamp error caused by signal acqusition procedure is 10 μs.
Computer simulations were carried out with MATLAB. To simulate a particular scenario, true
ˆ , cˆ and presetted value Cb was used to generate DTOA
value λ ^ X Q , YQ ,T , I ,Vr ` with erroneous R
observation value b̂ and then would be forgotten, see literature(9) for detailed calculation method.
ˆ , cˆ, C and b̂ known , λˆ
Then with R
could be resolved by NLS method withXˆ , Yˆ ,Tˆ, Iˆ,Vˆ

^

b

Q

Q

r

`

reference to Equation (4) and (5). For each scenario, simulations were carried out for 100 times.
CRLB on the DTOA error covariance matrix was calculated for every scenario with Equation ( 8),in which error covariance matrix N refers to Equation (12).
4.2 Influence of Error Sources on Simulation Results
This is Scenario 1, set a=50. In this 200×200m target area, the projectile hits the ground at point
Q(80,90,0)m with an elevation angle π/3 and azimuth angle 0, and the impact velocity is 700m/s.
Then the positioning error of the projectile impact point is simulated in the following cases.
Scenario 1.1 Only DTOA error caused by signal acquisition process is taken into account, 100 timesof simulations are carried out to calculate impact point position, its mean value and one standard
derivation (1-SD) ellipse, then actual value and CRLB ellipse are shown to compare in Figure 3(a).
Scenario 1.2 Only DTOA error caused by sensor position error is taken into account, and simulation
results are shown in Figure 3(b).
Scenario 1.3 Only DTOA error caused by wind, namely sound speed error is taken into account, andsimulation results are shown in Figure 3(c).
Scenario 1.4 All DTOA error sources in Scenerio 1.1~1.3 are taken into account, and results areshown in Figure 3(d).

3350

(a)

(b)

(c)

(d)

Figure 3 – Simulation results for different error source. (a)-(d) show the estimation results of
Scenerio 1.1-1.4. Dots donate estimates from 100 times of simulation value , crosses donate meanvalues of simulation, squares denote the actual values, solid lines denote CRLB error ellipses, anddashed lines denote 1-SD error ellipses.
As can be seen from the above figure, error discription method of CRLB fits the simulation ofprojectile impact point well. To targets of large area, the influence of arrival time acquisition error isthe smallest, the sensor coordinate error is the second, while the influence of wind speed is thegreatest. The total error is mainly affected by the wind speed.
In this scenario, the elevation angle θ=π/3 and the azimuth angle ϕ=0, so the projection of
trajectory parallels the X-axis, so the X-coordinate is mainly affected. The arrival time of Sensor 1 issubtracted in every equation of DTOA algorithm, resulting in the asymmetry of simulated impact
point coordinate distribution.
4.3 Positioning Error Vary with Trajectory Parameters
This is Scenario 2, and the basic setting is the same as Scenario 1. All error sources are
considered, and specific parameters vary in a certain range to obtain the distribution of CRLB standdeviation of DTOA.
Scenario 2.1 Scan X Q and Y Q, the CRLB standard deviation distribution of XQ , YQ and RQ is shownin Figure 4(a), (b) and (c).
Scenario 2.2 Scan θ and ϕ, the CRLB standard deviation distribution of XQ , Y Q and R Q is shown in
Figure 4(d), (e) and (f).

(a)

(b)

(c)

(d)

(e)

(f)

(g)
(h)
(i)
Figure 4 – Impact point error distribution versus hit position, angles and speed. (a) (d) (g) are3351

x-coordinate error distribution, (b) (e) (h) are y-coordinate error distribution and (c) (f) (i) are
position error distribution. (a)-(c) are error distribution versus XQ and YQ , (d)-(f) are error
distribution versus θ and ϕ, (g)-(i) are error distribution versus Vr and C b.
As is shown in the above figure (a)-(c), when the projectile hits the ground at different position,
the error of impact point is different. Near to Sensor 9 (at (200,200,0)m), the error of impact point islarge because that the DTOA algorithm is mainly based on time arrival difference between Sensor 1
and others. Points near Sensor 9 is relatively far from Sensor 1, forming a longer distance for shock
wave to travel, so the accumulated error is larger, and the impact point error is larger. When theimpact point is close to a sensor, the total error will reach the local minimum, because at thissituation, the arrival time acquisition error is minimal, and its contribution to the total error isminimal. The X-coordinate error of impact point is larger than the Y-coordinate error in general,
because the projection of the ballistic trajectory is parallel to the X-axis, and the spatial resolution of
the X-axis direction decreases(14) (Imaginne in extreme cases that θ=0 and ϕ=0, and thex-coordinate estimation error increases indefinitely).
Figure 4 (d-f) shows the error distribution versus elevation angle θ and azimuth angle ϕ, it can beseen that in general, the more the trajectory is prependicular to the ground, the less impact point
error is.
Figure 4 (g-i) are the error distribution versus ballistic const Cb and projectile hitting speed V r. As
the speed gets slower, the shock wave angle become greater, so the positioning error become greater.
4.4 Impact point positioning error versus sensor spacing
This is Scenaric 3, and the basic setting is the same as Scenario 1. Change the sensor spacing aand the impact point coordinate (1.6a, 1.8a, 0), the relationship between the error of impact pointand different error sources is obtained. The result is shown in Figure 5.

Figure 5 – Impact point error versus target scale for different error sources
Figure 5 shows that when the target size (or the sensor spacing) is greater than a certain threshold(here the threshold a=0.63, target size 2.5*2.5m), the influence of wind speed on the total error
become the greatest, and linearly increases with the increase of target size. When the target surface
amplitude is small, the sensor position error and the time acquisition error take the lead successively.
At this time, the sensor position is usually measured by rulers, with high accuracy, so the timeacquisition error has dominant influence. This conclusion is in according to the common set up inearly articles.

5. CONCLUSIONSIn this paper, NLS algorithm is applied to realize the localization of the impact point ofsupersonic flight projectile, and the CRLB method for the positioning error estimation is improved,
so that it is still applicable under the influence of the wind within level 3. Simulation results showthat the improved CRLB method can truly reflect the distribution law of simulation data and can beused as a basis for analyzing the error distribution characterist ics of impact points under different
circumstances. According to the improved method, analysis is conducted on which error sources are
dominant under different target scales, and the conclusion has certain guiding significance for theconfirmation and evaluation of field array scheme.

3352

REFERENCES
1. Rathakrishnan E. Applied gas dynamics. 2nd ed. Chichester, West Sussex, UK: John Wiley & Sons;
2019.
2. Anderson JD. Hypersonic and high-temperature gas dynamics. Reston, Virginia, USA: American
Institute of Aeronautics and Astronautics; 2006.
3. Cannella M, Cappa P, Sciuto SA. A novel approach for determining the trajectory and speed of asupersonic object. Measurement Science & Technology. 2003;14(8):1485.
4. Mäkinen T, Pertilä P. Shooter localization and bullet trajectory, caliber, and speed estimation based on
detected firing sounds. Applied Acoustics. 2010;71(10):902-13.
5. Levanon N. Acoustic hit indicator. IEEE Transactions on Aerospace & Electronic Systems.
2001;37(1):304-9.
6. Weinacht P, Cooper GR, Newill JF. Analytical prediction of trajectories for high-velocity direct-firemunitions. Aberdeen Proving Ground, MD: U.S. Army Research Laboratory; 2005.
7. Lo KW. Curvilinear trajectory estimation of a supersonic bullet using ballistic shock wave arrivals at
asynchronous acoustic sensor nodes. Journal of the Acoustical Society of America. 2017;141(6):4543.
8. Lindgren D, Wilsson O, Gustafsson F, Habberstad H, editors. Shooter localization in wireless sensor
networks. International Conference on Information Fusion; 2009.
9. Lo KW, editor Supersonic bullet trajectory estimation using ballistic shock wave arrivals at an acoustic
sensor array. Proceedings of ACOUSTICS 2016; 2016 9-11 November 2016; Brisbane, Australia.
10. Völgyesi P, Balogh G, Nádas A, Nash CB, Lédeczi Á, editors. Shooter localization and weapon
classification with soldier-wearable networked sensors. International Conference on Mobile Systems;
2007.
11. Lo KW, Ferguson BG. Localization of small arms fire using acoustic measurements of muzzle blast
and/or ballistic shock wave arrivals. Journal of the Acoustical Society of America.
2012;132(5):2997-3017.
12. Lindgren D, Wilsson O, Gustafsson F, Habberstad H. Shooter Localization in Wireless Microphone
Networks. Eurasip Journal on Advances in Signal Processing. 2010;2010(1):1-11.
13. Sengupta SK. Fundamentals of Statistical Signal Processing: Estimation Theory1993.
14. Grasing D, Ellwood B. Development of acoustic sniper localization methods and models. Proc Spie.
2010;7693:769312--13.

3353