On geometry of the acoustic sniper localization

(1)

On geometry of the acoustic sniper localization

E. Danicki^†, A. Kawalec^‡, J. Pietrasi`nski^‡

†Polish Academy of Sciences IPPT, Warsaw 00-049

‡Military University of Technology, Warsaw 00-908 POLAND

edanicki@ippt.gov.pl

Abstract

A supersonic bullet generating a cone of weak shock wave in a surrounding air is considered.

The characteristic shape of the acoustic signal is observed by two directional acoustic sensors of known spatial positions. The measured data include the wave arrival time and its propagation direction. The subject of this paper is to show that this set of data is necessary and sufficient for evaluation 1) the trajectory of the bullet (straight by assumption), and 2) its velocity. An analytical solution is presented and its sensitivity to measurement errors, depending on the distance to the bullet trajectory, its azimuth and elevation with respect to sensors.

1 Introduction

There is no much time left to destroy the sniper post after first and perhaps last detection of characteristic acoustic signal generated by a supersonic bullet. And there rarely are the other data that may help to pin-point the sniper.

Thus even approximate localization is better than none in a hostile environment. The localization problem should be solved using the data obtained from a single event measurement. And the measurement system should be as simple as possible to reduce its cost and rise its robustness in a hostile environment.¹

One may imagine a system of two or three microphones mounted on a typical military ve- hicle as its auxiliary equipment, perhaps com- plemented by one or more mounted on personal helmets to deliver auxiliary data. It could be a microphone thrown on the hostile field provided that its position is known from certain calibration measurements (i.e. acoustic). Such a remote sensor would transmit its raw detected signals to the system computer. The feasibility of such a system is studied below.

2 The localization problem

A supersonic bullet generates a weak shock wave in air.² For constant bullet velocity and straight its flight path, the shock wave has a shape of a cone with the bullet on its tip and conical angle ϑ depending on the bullet velocity v relative to the sound velocity in air c (Fig. 1) sin ϑ = c/v < 1. (1) The cone axis k coincides with the bullet path (and velocity direction) leading us to the sniper position (in the frame of the above assumption).

Finding the cone axis thus solves the simplest (canonical) case of the acoustic sniper localization problem. It is worthh tp note that:

• the shock wave propagates outward with the sound velocity c, in direction n normal to the shock wave conical surface,

(2)

2θ n

r r

1

n

2 1

2

weak shock wave

v k

c

sniper ^t

t

1 2

d t

2

+d/c

directional measurement

Fig. 1 Schematic measurement system (two directional microphones of known positions) of a weak shock wave generated by a supersonic bullet.

• the acoustic signal is detected when the shock wave hits the microphone (the detected is the waveform characteristic for a weak shock wave; this can improve the detection accuracy),

• it is detected at the same time (simulta- neously) if the microphones reside on the same shock wave cone, and

• the observation point shifted by the dis- tance d from the cone surface will detect the signal at different time, delayed by

∆t = d/c, (2)

• thus displacing the point by d = c∆t back along the local shock wave propagation direction n places the point on the same shock wave cone again (Fig. 1).

Let us take assumption that two directional microphones, placed at the known spatial positions r_1,2detect the passing shock wave at times t₁ and t₂ (of intensity and shape dependent on the distance to the bullet trajectory, assumed straight) and simultaneously measure the sound propagation directions (which are the outward normals to the shock wave cone) n₁, n₂ respec- tively. Our problem is to find the cone axis from these two measurements. It may be instructive to search it in the way presented below.

3 Fundamental solution

In the following analysis, the shock wave cone that includes the measurement point r₁ (a vector in cartesian space) is considered. The other measurement, r₂ resides on different cone, but

shifting it back along the normal to the cone surface n₂ by the distance c(t₂− t₁), places this point (now denoted by R₂) on the considered cone surface where the measured shock wave ar- rival time would be t⁰₂ = t₁ instead of t₂. Thus one has the reduced measurements:

R₁ = r₁, n₁; t₁,

R₂ = r₂− c(t₂− t₁)n₂, n₂; t⁰₂ = t₁. (3) As known, a plane can be determined by its normal n and a point x₀ residing on it;

(x − x₀) · n = 0 (4) is the plane equation satisfied by any point x = (x₁, x₂, x₃) residing on the plane p. The di- rectional measurements of n₁, n₂at points r₁, r₂ allows us to construct two independent planes p₁, p₂ tangential to the chosen shock wave cone along the yet unknown lines l_1,2:

(x − R_i) · n_i= 0, i = 1, 2 (5) are equations of these planes. They cross each other along the line l satisfying simultaneously both the above Eqs. (5). The third directional measurement (if exist) would solve the problem, as three tangential planes have a common point that is the cone tip O.

It is shown below that the point O can be found from only two directional measurements, as the cross section of the line l with the secant plane constructed from the available measurements r_i, n_i, i = 1, 2. First note that the vector l_i× n_i is tangential to the normal cross section of the cone that is a circle of radius c_i (Fig. 2).

In other words, it is coplanar with the circle plane normal to the cone axis k. Let’s consider the secant plane spanned on certain lines resid- ing on the cone surface, l₁ and l₂. The normal vector to this plane is l₁× l₂. It is shown below that this plane is parallel to n₁− n₂. It means that

(n₁− n₂) · l₁× l₂ = 0. (6) Proof. It is known from the vector calculus that a·b×c = c·a×b and a×b = −b×a, which allow us to rewrite the above in the form (k n_i k= 1)

l₂· n₁× l₁

| {z }

s1

= l₁· l₂× n₂

| {z }

s2

∼k l₁kk l₂k, (7)

(3)

where s_i are tangential to the corresponding normal cross section circles of the cone (having radius c_i, as shown in a perspective view from the point O, Fig. 2). Noticing that

l_i = c_i+ α_ik, (8) where α_i are certain scalars (k is normal to the cone cross section, coplanar with s_i), the last equation can be rewritten in the form

c₂· s₁= c₁· s₂ =k · kk · k cos(α − π/2), (9) where c_i, s_i are all coplanar vectors and α is the angle between c₁ and c₂ - the radiuses of the cone normal cross sections, Eq. (8). This concludes the proof.

The system of two (different by assumption) lines coplanar with the secant plane: R₁− R₂ and n₁− n₂, determine this plane by its normal n = (n₁− n₂) × (R₁ − R₂) and the point R₁ (or R₂) residing on it. Three different planes:

the above determined secant plane and two tangential planes (n₁, R₁) and (n₂, R₂) cross each other at the cone tip O (a vector in the cartesian space of measurements), satisfying

(O − R_i) · n_i = 0, i = 1, 2,

(O − R₁) · (n₁− n₂) × (R₁− r₂) = 0. (10) The cone axis k can be found as normal to two vectors: s₁ and s₂, where s_i = l_i× n_i and l_i = R_i−O, evaluated uniquely from the microphone positions in cartesian space (vectors R_i) and the above obtained cone tip position O:

k = ±[n₁× (R₁− O)] × [n₂× (R₂− O)], (11) again normalized and directed to obtain k·(R₁− O) > 0. Particular cases of n₁ − n₂ k R₁ − R₂ and R₁ − O k R₂− O (thus n₁ = n₂) are excluded from the analysis here. Alternatively, k is the cross-section of two planes normal to p_i and spanned over n_i, l_i, i = 1, 2.

The above results yield another important in- formation - the conical angle ϑ

sin ϑ = −n_1,2· k, 0 < ϑ < π/2, (12) and thus the bullet velocity v on the strength of Eq. 1 (k k k= 1 assumed).

n

n R

n -n r=R

1 2

1

2

l

2

1

O

k

shock wave cone

tangential planes

1 l

r

2

1 view from

O

p

₁

2

the

c

s

₁

2

c s

2

α

1

Fig. 2 Evaluation of the cone tip O, and sub- sequently the cone axis k, from two directional measurements: r1, n1and the reduced one R2, n2.

Fig. 3 shows the computed example of the above analysis, using the Matlab code presented in the Appendix below with short description.

Another measurement point r₃ is shown that does not reside on the r₁ cone; t₃ 6= t₁. This measurement is used to check Eq. (14) presented later below. Two blue lines are also shown in the figure besides k. They are k⁰ evaluated for changed (by measurement error, for instance) t₂ by 2% of k r₁− r₃ k /c.

4 Auxiliary measurements

Omnidirectional measurements are technically easier than directional ones. Thus it may be advantageous to set one or more auxiliary omnidirectional microphones in order to check or improve the above discussed sniper localization.

They deliver the measurements

r_i; t_i, i = 3, 4, ... (13) where r_i is the known microphone position in cartesian space, and t_i is the corresponding arrival time of the shock wave. Note that gener-

(4)

−0.5

0

0.5

−0.4

−0.2 0 0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4

n

1

r

1

P k

O r3

r

2

R

2

n

2

Fig.3 Evaluated sniper direction k from two direc- tional measurements; O is tip of the shock wave detected at r1. The cone normal cross section (the circle with r1 on it) has center P . The second measurement point r2 moved back by c(t2 − t1) along the normal n2 places R2 on the same cone.

ally, t_3,4,... 6= t₁ thus another shock wave cone is detected than the above considered one with r₁, O on its surface, Fig. 4. Consider the mea- surement No. 3. If correct, ξ₃= 0 results:

ξ₃ = k · (r₃− O) sin ϑ − b₃cos ϑ + c(t₃− t₁), b₃ =

q

k r₃− O k²−[k · (r₃− O)]² (14) (k k k= 1). First equation results from the below one multiplied by sin ϑ

k · (r₃− O) − b₃cot ϑ = v(t₁− t₃), v = c/ sin ϑ;

another auxiliary measurement data would satisfy analogous equations.

5 Directional accuracy

Let the shock wave arrival times (t_1,2,3,...) be measured with sufficient accuracy by all microphones residing at exactly known positions,

r_i, i = 1, 2, 3, ..., but there is some error in directional measurements resulting in slightly in- accurate n_1,2. This naturally results in certain error of the evaluated cone axis k. The problem considered here is whether the auxiliary (omnidirectional) measurements could correct k, making the localization more accurate.

k r O

3

b

3 3

1

v(t -t )

₁ ₃

shock wave front of

θ

r

₁

θ

c(t -t )

Fig. 4 Spatial position of the auxiliary omnidirectional measurement point with respect to the shock wave cone detected at r1(here, the case t1> t3).

First, the dependence of ξ_i on the n_1,2 error is analyzed. Assuming n₁ accurate, let n⁰₂ be varying with respect to its correct value in such a manner that the vector arrow draws a circle in the plane perpendicular to the correct n₂:

n⁰₂ = n₂+ ²(e₁cos γ + e₂sin γ),

e_1,2⊥ n₂, and e₁ ⊥ e₂, (15) again normalized to obtain k n⁰₂ k= 1; e_1,2 are also normalized vectors (we choose e₁ parallel to R₂− O and e₂ = e₁× n₂).

The parameter γ ∈ (0, 2π takes several values in the computed examples on Fig. 5 to present the resulting variation of the cone axis k⁰ (the cone tip O⁰ also varies with ², γ, but it is not shown in the figure; the vector k is attached to O to show only its directional variation that is more important for the sniper localization than the lateral shift O⁰− O). In general, k⁰ draws a closed curve around the correct k. It is interest- ing how the variation of n⁰₂ affects the equation ξ₃ = 0 for the third auxiliary measurement (accurate one in the current analysis). To illustrate this, different color marks n⁰₂ in Fig. 5 depend- ing whether the resulting ξ₃ > 0 or ξ₃< 0. It is evident that, for small ², there is certain γ = γ₂₃ for which ξ₃= 0 and invariant on n⁰₂(γ₂₃; ²)−n₂, being a function of ² (green dots in the figure).

Similar results can be obtained for variation of direction n⁰₁with respect to n₁(note however that variation of n⁰₂ changes R₂ into R⁰₂, while

(5)

−0.4

−0.2 0

0.2 0.4

0.6

−0.4

−0.2 0 0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4

k

O r₃

error 10 times exaggerated n_i and resulting k variations marked by green points satisfy ξ₃=0.

Fig. 5 The influence of the directional errors of n1,2(assumed 5% in directions perpendicular to ni) on the cone axis direction k (magenta dots at the end of k, in correct scale).

R₁ remains unchanged under n₁ variation, and that another e_1,2 are chosen), yielding its own invariant γ₁₃. Generally, for more omnidirectional measurements, several invariant γ_1i and γ_2i result. An important question arises: can omnidirectional measurements elim- inate the need of directional measurements or, at least, correct their inaccuracy?. This would be the case if two lines of invariant γ_2icross each other at certain point indicating the correct direction n₂, provided that n₁ is correct. Per- haps more omnidirectional measurements would allow to obtain both simultaneously corrected n_1,2, and thus the accurate k. This problem will be studied in the future, in the frame of pertur- bation analysis concerning infinitesimal ².

6 Conclusions

Presented above geometrical research of the acoustic localization for the case of supersonic bullet of straight path can be summarized as follows:

−0.6 −0.4 −0.2 0 0.2 0.4 0.6

−0.5

−0.4

−0.3

−0.2

−0.1 0 0.1 0.2 0.3

r₃

r₁ r₂

Fig. 5 cont. Top view. The assumed directional variation (simulated error) of ni is shown 10×

exaggerated. The sprigs at the variation circle of ni depict |ξ3|. Green dots mark variation of ² in direction resulting in ξ₃ = 0. The corresponding k variation (also green dots) indicates approximately the direction from the cone axis to ri.

1. two directional measurements are sufficient for evaluation of the direction to sniper. Mea- sured are: the arrival time of the weak shock wave at the measuring microphones known position, and the direction of the wave propagation (which is normal to the shock wave front having the form of a cone),

2. the bullet velocity is also evaluated simultaneously (this gives certain information of the arm used),

3. additional omnidirectional measurement of the shock wave arrival time at the known microphone position can improve the sniper localization accuracy.

Future investigations will include:

• generalization of the above theory to the case of curvilinear bullet path and its variable velocity,

• the algorithm optimization and the error anal- ysis of the solution depending on the positions of microphones,

• further study of the system feasibility, partic- ularly concerning the conditions on directional measurements.

Acknowledgements

Authors are indebted to Dr. R. Klemm for encour- agement. This work was supported by the State Sci- ence Committee under Grant O T00B 003 25 (2004).

(6)

References

1J. B`edard, S. Par`e,T. Johson, Ferret Small Arms’

Fire Detection System, this Symposium, 2004

2L.D. Landau, E.M. Lifshitz, Mechanics of continu- ous media, Mir Publ., Moscow, 1954 (Polish transl., PWN, 1958)

Appendix - Matlab code

r0=[.2 .4 .5];a0=pi/9.5;%observer’s space aa=[1 0 0;0 cos(a0) sin(a0);0 -sin(a0) cos(a0)];

r0=0*r0;aa=diag([1 1 1]);view([0,90]);

th=pi/6;c=1;v=c/sin(th);a=2*pi*(0:3:360)’/360;

kol=[sin(th)*[cos(a) sin(a)] ones(121,1)*cos(th)];

kol=[ones(121,1)*r0+kol]*aa’;pl(kol,1.5,’c’);hold on;

r0=r0*aa’;pp(r0+[0 0 cos(th)]*aa’,’+’,10,’b’);

%mic. No. 3 (omnidirectional)

a3=pi/2.5;z3=.5;t3=.5; r3=r0+[z3*[cos(a3) sin(a3)]*sin(th) ...

z3*cos(th)+v*t3]*aa’; pl([r0;r3],.5,’k’);pp(r3,’o’,5,’k’);

text(r3(1)+.02,r3(2)-.03,r3(3)-.02,’r_3’)

%Mic. No. 1 (directional)

a1=-pi/4.5;t1=0;z1=1; r1=r0+z1*[[cos(a1) sin(a1)]*sin(th) ...

cos(th)]*aa’; pp(r1,’.’,15,’b’); n1=z1*[[cos(a1) ...

sin(a1)]*cos(th)-sin(th)]; n1=n1/sqrt(n1*n1’)*aa’;

pl([r0;r1;r1+n1/4],1.5,’b’);

text(r1(1)+.01,r1(2)+.02,r1(3)+.08,’r_1’)

text(r1(1)+n1(1)/5-.1,r1(2)+n2(2)/5,r1(3)+n1(3)/5-.01,’n_1’)

%Mic. No. 2 (directional) at r2; reduced at R2

a2=6*pi/5;z2=.6;t2=.3; n2=z2*[[cos(a2) sin(a2)]*cos(th) -sin(th)];

n2=n2/sqrt(n2*n2’)*aa’; r2=r0+z2*[[cos(a2) sin(a2)]*sin(th) ...

cos(th)]*aa’+n2*(t2-t1);

pp(r2,’.’,5,’y’);pl([r0+n2*(t2-t1);r2;r2+n2/4],1,’y’);

R2=r2-n2*t2;pl([r0;R2;R2+n2/4],1.5,’c’);pp(R2,’.’,10,’c’);

text(r2(1)+.01,r2(2)-.04,r2(3),’r_2’) text(R2(1)-.03,R2(2)+.04,R2(3)+.06,’R_2’)

text(r2(1)+n2(1)/5,r2(2)+n2(2)/5,r2(3)+n2(3)/5-.05,’n_2’)

%Evaluation of the cone tip O and axis k

n=vp(n1-n2,r1-R2);O=[n1;n2;n]\[n1*r1’;n2*R2’;n*r1’];O=O’;

pp(O,’.’,10,’b’);text(O(1)-.01,O(2)+.02,O(3)-.08,’O’) k=vp(vp(R2-O,n2),vp(r1-O,n1));k=k/sqrt(k*k’); if ((r1-O)*k’)<0;k=-k;end;pl([O;O+k],2,’r’);

text(O(1)+k(1)+.01,O(2)+k(2)+.02,O(3)+k(3)+.06,’k’) P=[k ...

0;[diag([1 1 1]) -k’]]\[r1*k’;O’]; P=P(1:3)’;pp(P,’x’,5,’k’);

text(P(1)+.01,P(2)-.02,P(3),’P’)

x=O+(R2-O)*sqrt(((r1-O)*(r1-O)’)/((R2-O)*(R2-O)’));

pp(x,’.’,5,’g’);pl([R2;x;P;r1],1,’g’);

%does (r3,t3) satisfy equations?

r=r3-O;b=sqrt(r*r’-(r*k’)^2);

t=sqrt((P-O)*(P-O)’)/sqrt((r1-O)*(r1-O)’);

x=(r3*k’)*sqrt(1-t^2)-b*t-c*t3;display(x)%should be =0 ko=k;

%THE ERROR ANALYSIS

%ERROR in t2 = 2% of |r1-r2|/c

R2o=R2;t2o=t2;t2=t2o+.02*sqrt((r1-r2)*(r1-r2)’);R2=r2-t2*n2;

n=vp(n1-n2,r1-R2);O=([n1;n2;n]\[n1*r1’;n2*R2’;n*r1’])’;

pp(O,’.’,5,’m’);k=vp(vp(R2-O,n2),vp(r1-O,n1));k=k/sqrt(k*k’);

if ((r1-O)*k’)<0;k=-k;end;pl([r0;r0+k],.5,’b’);

t2=t2o-.01*sqrt((r1-r2)*(r1-r2)’);R2=r2-t2*n2;

pp(O,’.’,5,’m’);k=vp(vp(R2-O,n2),vp(r1-O,n1));k=k/sqrt(k*k’);

if ((r1-O)*k’)<0;k=-k;end;pl([r0;r0+k],.5,’b’);

t2=t2o;r2=r2o;k=ko;R2=R2o;

%end error in t2

%ERROR in n1

N1=n1;R1=r1;%error axis definitions:

e1=[-sin(a1) cos(a1) 0]*aa’;e2=[[cos(a1) sin(a1)]*sin(th) ...

cos(th)]*aa’; pl([R1+e1/4;R1;R1+e2/4],.5,’y’); ee=.05;

for i=1:46;eb=(i-1)*2*pi/45;

n1=N1+ee*(e1*cos(eb)+e2*sin(eb));n1=n1/sqrt(n1*n1’);

k=vp(vp(R2-O,n2),vp(r1-O,n1)); if ((r1-O)*k’)<0;k=-k;end;

k=k/sqrt(k*k’);pp(r0+k,’.’,4,’m’); P0=[k 0;[diag([1 1 1])...

-k’]]\[r1*k’;O’];P0=P0(1:3)’;

r=r3-O;b=sqrt(r*r’-(r*k’)^2);t=sqrt((P0-O)*(P0-O)’)...

/sqrt((r1-O)*(r1-O)’);

x=(r3*k’)*sqrt(1-t^2)-b*t-c*t3;r=ee/sqrt(1-(ko*k’)^2);

n=2.5*(n1-N1);%instead of (n1-N1)/4 --> 10*exagg. in the fig.

if x>0;pl([R1+N1/4;R1+N1/4]+[1;r]*n,.5,’k’);else;

pl([R1+N1/4;R1+N1/4]+[1;r]*n,.5,’r’);end;

if i==1;xs=x;ns=n1;x0=[];else;if x*xs>0;xs=x;ns=n1;else;

x0=[x0;(xs*n1-x*ns)/(xs-x)];xs=x;ns=n1;end;end;

end;%x0,

for i=1:6;n1=(i-1)*x0(2,:)+(6-i)*x0(1,:);n1=n1/sqrt(n1*n1’);

pp(R1+N1/4+2.5*(n1-N1),’.’,5,’g’);%10* exaggerated n=vp(n1-n2,r1-R2);O=([n1;n2;n]\[n1*r1’;n2*R2’;n*r1’])’;

k=vp(vp(R2-O,n2),vp(r1-O,n1));%pp(O,’.’,5,’b’) if ((r1-O)*k’)<0;k=-k;end;k=k/sqrt(k*k’);

pp(r0+k,’.’,5,’g’);end;n1=N1;r1=R1;k=ko;%endLOOPn1

%ERROR in n2

e1=[-sin(a2) cos(a2) 0]*aa’;e2=[[cos(a2) sin(a2)]*sin(th) ...

cos(th)]*aa’; pl([r2+e1/4;r2;r2+e2/4],.5,’y’);N2=n2;ee=.05;

for i=1:46;eb=(i-1)*2*pi/45;n2=N2+ee*(e1*cos(eb)+e2*sin(eb));

n2=n2/sqrt(n2*n2’);R2=r2-n2*t2;pp(r2,’.’,1,’g’);

k=vp(vp(R2-O,n2),vp(r1-O,n1)); if

((r1-O)*k’)<0;k=-k;end;k=k/sqrt(k*k’); pp(r0+k,’.’,4,’m’);

P0=[k 0;[diag([1 1 1]) -k’]]\[r1*k’;O’];P0=P0(1:3)’;

r=r3-O;b=sqrt(r*r’-(r*k’)^2);

t=sqrt((P0-O)*(P0-O)’)/sqrt((r1-O)*(r1-O)’);

x=(r3*k’)*sqrt(1-t^2)-b*t-c*t3;r=ee/sqrt(1-(ko*k’)^2);

n=2.5*(n2-N2);%10*exagger.

if x>0;pl([r2+N2/4;r2+N2/4]+[1;r]*n,.5,’k’);else;

pl([r2+N2/4;r2+N2/4]+[1;r]*n,.5,’r’);end;

if i==1;xs=x;ns=n2;x0=[];else;if x*xs>0;xs=x;ns=n2;else;

x0=[x0;(xs*n2-x*ns)/(xs-x)];xs=x;ns=n2;rs=r2;end;end;

end;for i=1:6;n2=(i-1)*x0(2,:)+(6-i)*x0(1,:);n2=n2/sqrt(n2*n2’);

pp(r2+N2/4+2.5*(n2-N2),’.’,5,’g’);R2=r2-n2*t2;

k=vp(vp(R2-O,n2),vp(r1-O,n1));%pp(O,’.’,5,’b’)

if ((r1-O)*k’)<0;k=-k;end;k=k/sqrt(k*k’);pp(r0+k,’.’,5,’g’);end;

n2=N2;k=ko;

%Support functions

function vectprd=vp(n1,n2) v=[n1(2)*n2(3)-n1(3)*n2(2) n1(3)*n2(1)-n1(1)*n2(3) n1(1)*n2(2)-n1(2)*n2(1)];

vectprd=v/sqrt(v*v’);

function pl3lin=pl(f,i,c)

set(plot3(f(:,1),f(:,2),f(:,3),’LineWidth’,i),{’Color’},{c});

function pl3pkt=pp(r,s,i,c)

set(plot3([r(:,1)],[r(:,2)],[r(:,3)],s,’Markersize’,i),{’Color’},{c});

On geometry of the acoustic sniper localization