A method for identifying a muzzle blast

[0001] The present invention relates to law enforcement and security enhancing technologies, and more specifically to methods and systems for determining the starting point and direction of supersonic missiles. Methods and systems are capable of determining and disambiguating the position of the shooter even at large distances between the shooter and the sensor, and when there is no signal or when only a weak shot sound is received.

Background of the invention [0002] Systems and methods are known that can determine the overall direction and trajectory of supersonic bullets, such as small arms and artillery shells, by measuring parameters associated with the shock wave produced by the projectile. One such system, described in US Pat. No. 5,241,518, comprises at least three spaced apart sensors, each sensor having three acoustic transducers arranged in a plane. Sensors generate signals in response to the shockwave, which are related to the azimuth and angle of aiming, as well as to the start point of the shockwave. Shockwave-only measurements are unable to determine the distance between the sensor (s) and the origin of the shockwave. Distance information is usually obtained from the muzzle flame and muzzle blast.

[0003] The azimuth and angle of the shooter with respect to the position of the sensor is typically determined by measuring the arrival time information (TOA) from the muzzle signal and the shockwave signal for each sensor. Each sensor encounters signals at different times and generates a signal in response to outlet pressure and shock wave. The signals from individual sensors are processed and the direction (azimuth and aiming angle) from the sensor (sensors) relative to the starting point of the barrel and shock wave, and thus the projectile trajectory can be determined.

[0004] Conventional systems use microphones that can be spaced relatively close (e.g. at a distance of 1 m) or very dispersed (e.g. mounted on a vehicle or carried by soldiers on the battlefield) and can measure the discharge pressure and the shockwave all around in their appropriate positions. However, if the sensors are not relatively spaced apart and / or the trajectory is not within the antenna, very high accuracy of time measurements is needed to obtain accurate solutions based only on the shock wave and special techniques are required. [0005] The large size of the antenna can be a serious disadvantage of systems mounted e.g. on a vehicle. In addition, systems with minimal time resolution can generate ambiguous solutions in which information about the time of arrival of the shock wave by a given set of sensors is almost identical for the two positions of the shooter being their mirror image.

[0006] Conventional algorithms require at least 4 shockwave and muzzlewave detection points so that a 4x4 matrix can be converted to mapping a flatwave based on a shockwave TOA. Small errors in TOA and shockwave wave determination can cause significant errors in distance estimation. In addition, conventional algorithms assume a constant projectile speed along the projectile trajectory, which gives inaccurate distance estimates for distant shots fired from a distance of more than about 300 m.

[0007] Accordingly, there is a need for rapidly converging algorithms capable of accurately estimating distant shooter distances. There is also a need for disambiguation of solutions based only on the shock wave in order to determine the direction of the shooter. There is also a need to extract muzzle signals that may be unclear by acoustic features not associated with muzzle blast.

Summary of the Invention [0008] The present invention provides a method of muzzle blast identification according to claim 1. In a further aspect there is provided a muzzle blast identification system according to claim 8.

[0009] The invention addresses the disadvantages of known solutions and in various embodiments, provides methods and systems for estimating shooter range in the case of distant shots, especially when the muzzle signals are either weak or are detected in an insufficient number of detection channels. The disclosed methods and systems also improve the disambiguation of shot trajectory solutions based only on the shock wave, along with additional improvements obtained by including weak and / or uncertainly detected muzzle sound in the optimization process.

[0010] To prevent false signals from being interpreted as shockwave waveforms, the projectile trajectory can be eliminated as false if the acoustic energy of the measured shockwave waveform is less than a predetermined threshold in a fixed frequency band, e.g. in the frequency range between about 700 Hz and 10 kHz. Alternatively or additionally, the projectile trajectory can be eliminated as false if the time interval in which the measured shock waveform is positive is less than the minimum time or is greater than the maximum time, e.g. less than about 70 μs or greater than about 300 microseconds.

[0011] In preferred embodiments, the total energy can be determined by integrating the measured energy in the range, preferably omitting the fragments in the detected signal caused by the shockwave echoes. Preferably, the peak signal value can be determined in the range generating the maximum total energy, and if the peak signal value is greater than the measured total energy range by a fixed aspect ratio, the peak signal value can be identified as associated with the muzzle signal.

[0012] Embodiments of the invention may include one or more of the following features. The error distribution of the time calculation for the antenna and / or acoustic sensors can be associated with gain changes, sampling changes, and antenna sensor position changes. The confidence level of disambiguation depends on the size of the antenna, as a result of which the smaller antenna requires greater measurement accuracy. If there are two ambiguous solutions, the disambiguated projectile trajectory is selected based on the residual ratio for the two ambiguous solutions.

[0013] In other preferred embodiments, the arrival time difference (TDOA) for the sensor pairs can be determined by designing the sensor that first encounters the shockwave as a reference sensor, and setting the first flip-flop of the timing system when the amplitude of e.g. the beginning of the signal only the shock wave for the reference sensor exceeds the threshold value. This first flip-flop activates the start counters for each of the other sensors, with the counter in each of the other sensors running until the corresponding sensor encounters the shockwave. When one of the other sensors encounters e.g. the initial fragment of the shockwave only signal, it sets a second flip-flop for this sensor, which stops the start counter for that sensor. The TDOA values ​​for the other sensors relative to the reference sensor are then recorded.

[0014] Further features and advantages of the present invention will become apparent in the light of the following description of preferred embodiments and claims.

Brief Description of the Drawings [0015] These and other features and advantages of the invention will become more apparent from the following description with reference to the accompanying drawings, in which the elements are marked with similar reference numerals and which may not be drawn to scale.

Fig. 1 schematically shows a cross section of a Mach cone intersecting with an antenna;

Fig. 2 schematically illustrates an example sensor array with 7 omnidirectional acoustic sensors;

Fig. 3 schematically shows the ambiguity inherent in determining trajectories only on the basis of a shock wave;

Fig. 4 schematically shows the probability density for measuring arrival time difference to determine Mach cone curvature;

Fig. 5 schematically shows the probability of correct disambiguation between shooter trajectories;

Fig. 6 is a flow chart of the correlation process;

Fig. 7 is a block diagram of the genetic algorithm used for correct disambiguation between shooter trajectories;

Fig. 8 is a block diagram of the process of discriminating between signals not related to a shock wave;

Fig. 9 schematically shows the shock wave arrival time (TOA) model;

Fig. 10 shows a block diagram of a distance estimation process;

Fig. 11 shows a block diagram of a genetic algorithm process for distance estimation.

Detailed Description of Some Embodiments Shown [0016] As described above in the Summary, the invention relates, in various embodiments, to methods and systems for estimating shooter range and disambiguating projectile trajectories. These systems and methods are particularly useful and advantageous when an insufficient number of parameters is detected to obtain an accurate solution or when such parameters cannot be reliably detected.

[0017] Trajectories of supersonic missiles are estimated only on the basis of the arrival time of the projectile shock wave, measured by a series of closely spaced sensors arranged in a "small" measurement space, referred to as the antenna. The measuring space is considered small if the sensor distance is 2 m or less. Once the projectile trajectory is determined, the shooter position is known except for the distance back along the trajectory. This distance can be calculated if the antenna also obtains the arrival time of the muzzle blast. However, the muzzle blast is not always detectable, so that an accurate solution based only on a shock wave is essential for determining the trajectory.

[0018] Referring now to Fig. 1, the shockwave surface is considered to be an expanding conical surface having an axis coinciding with the projectile's trajectory. The shockwave surface is also called Mach cone. To obtain a shock-wave-only solution, three properties must be determined based on arrival times measured by five or more antenna sensors: angle of incidence, radius of curvature, and spatial gradient of the radius of curvature of the expanding conical surface.

[0019] The angle of incidence of the conical surface that first reaches the antenna defines two possible relative angles (often called "ambiguous" angles) of the projectile trajectory relative to the angle of incidence on the antenna. The "ambiguous" angles will be described in more detail below with reference to Fig. 3. The radius of curvature of the conical surface at the antenna determines both the distance and the direction of the trajectory. The gradient of the radius of curvature along the path forming the surface determines the direction in which the projectile is moving, thus eliminating the "ambiguity" between two possible directions. Determining these three shockwave properties accurately and correctly deciding between two possible "ambiguous" trajectory angles requires very accurate measurements. For example, random errors should not be greater than about 1 μs to make the right choice between two alternative angular shooter positions.

[0020] The required accuracy can be estimated by considering the shockwave propagation characteristics of Fig. 1. Referring now to Fig. 2, the antenna 20 includes N sensors (N = 7) capable of determining the arrival times of the progressive conic shock wave. Since the trajectories of the incoming projectile can indeed be expected to originate anywhere, antenna elements 23 to 28 can be advantageously distributed evenly at C (Cxj, Cyj, Czj) on the spherical surface, with one element 22 located in the center on (Cx0 , Cy0, Cz0) so that homogeneous sensor apertures occur regardless of the angle of incidence. The moment when the first sensor, designated as the reference sensor, detects the progressing conical surface, is designated as t0. Other sensors detect the progressing conical surface at subsequent moments marked as ti. Sound propagation distances towards a conical surface are obtained by multiplying each of the time differences by the local speed of sound c, i.e. di = c^.(ti -to). If there are no measurement errors, then the conical surface passing through the reference sensor is also determined by other (N-1) sensors, where the three-dimensional coordinates of the N points ideally determine all the parameters of the shock wave cone. However, as mentioned above, errors in arrival time measurements and sensor coordinates can result in erroneous shock wave parameters and thus also projectile trajectories. The accuracy of the arrival time difference required to make the correct decision regarding both ambiguous trajectory angles will be further described.

[0021] The system preferably has features that ensure that it will not misinterpret non-ballistic signals, such as car noise, vibration, wind noise and electromagnetic interference, as a shooter. For example, the sensor mast can be mounted on a vehicle (not shown) with elastomer sleeves in mating connectors to prevent rattling. The sensors can be attached to the ends of the support rods with elastomer connectors having low resonance frequencies around 1 Hz to isolate them from the support rod vibration. The sensor support rods can be attached to a common hub that houses analogue electronics, which can also be attached to the sensor mast with elastomeric shock mounts to isolate it from mast vibration.

[0022] In addition, the following decision algorithm can be used to filter out signals that lack the signatures typically found in shockwave-derived signals. All values ​​are parameterized, i.e. relative, and can be tuned externally. The values ​​listed are for illustration only.

[0023] Referring now to Fig. 8, the process 800 determines if the detected signal is from a shockwave. The process 800 starts with step 802 and checks in step 804 whether the signal is a loud enough event to be considered an impact, e.g. if the peak signal value exceeds a given parameterized threshold of e.g. 500. If so, process 800 proceeds to step 806 and checks if there is a sharp transition from zero to peak signal value, ensuring that the transition to this peak value is not preceded by another signal of significant magnitude, e.g. 1/16 peak signal values.

[0024] If this is the case, the process 800 proceeds to step 808 and checks if the time between the shockwave minima and maxima is large enough, e.g. 200-400 μs. If so, process 800 proceeds to step 810 and checks if the magnitudes of the minima and maxima of the signal amplitudes are close, i.e. within 35% of each other. If this is the case, process 800 goes to step 812 and checks whether the transition of the pressure peak from the minimum signal peak to zero is sharp using essentially the same criteria as in step 806. If so, process 800 goes to step 814 and checks, whether the times between the maximum signal value and the zero crossing and between the zero crossing and the minimum signal value are comparable, e.g. within about 180 μs. If all steps give an affirmative answer, process 800 decides that the signal can be a shockwave and the signal is processed, step 816. On the contrary, if one of the six decision steps gives a negative response, the detected signal does not originate from the shockwave, step 818.

[0025] Referring again to Figure 1, it is assumed that the projectile trajectory coincides with the 5th axis x. The Mach angle is given by the equation θ = arcsin (1 / M), where M is the Mach number determined as the projectile speed V divided by the speed c. L is the characteristic length of the antenna. The radii of curvature of the cone at the two ends of the antenna are r1 and r2, respectively. The front view in the left half of the drawing shows how the curvature of r is measured₁. The distance d is equal to d = rpcos ^). The angle φ is defined by sin ^) = L / 2 r₁ so that for small angles φ we get φ ~ L / 2 r₁. The time difference between the points on the antenna surface measuring the curvature, cutting the conical surface with a radius of r₁ is equal to dt₁ = Δd / c = (r₁ - d) / c ~ r₁ φ / 2 c = L / (8-yyyy).

Difference in curvature measurement time at₂ = year₁ - I / sinf)) is given the same expression, with r2 replacing r1. Accordingly, dt2 = dt1 + L sin ^) / 8r₁ c.

[0026] Assuming unloaded measurement errors, i.e. assuming that the differences in time measurements dt1 and dt2 are randomly distributed values ​​having different average values ​​dt1 and dt2 but the same statistically determined standard deviation σ, average measurement values ​​at both ends of the matrix they correctly designate the local curvature that exists there. Sample distribution of measurement values ​​for time differences dt1 and dt2 are shown in Fig.

4.

[0027] The test measurement taken at end 2 is designated X. The radius of curvature at end 2 (radius r2) is smaller than at end 1 (radius r1). As a result, all measurements taken at the end of one have values ​​greater than X and will result in the correct decision that the curvature at the end of 1 is greater than at the end of 2. The probability that the correct decision will be made when the measurement result at the end of 2 is X data is dependencies:

^{P (r}1 < ^r2 ^{1 x)} = P2^(X) JP1^{(X) dX} x

where and p1 (x) = p2 (x) =

<img file="PL2204665T3_D0001.tif" />

(x-dt₂)² e ^σ

4ϊπσ (4-dt1?

e ^2S [0028] Integration after x and substitution of variables results in the following probability of making the correct decision:

^{P (r}and &lt;U) = - and

2yp -_¥ e^-at erf (u - a) du where a = dt_- - dt _{2 =} L³ sin (q)

42s

428r-'cs [0029] Referring now to Fig. 5, the probability of a correct decision or disambiguation confidence level is plotted for two sample antenna sizes, L = 1 m and L = 2 m, as a function of the nearest approach point (CPA) r between the projectile trajectory and antenna 20. An audio velocity of c = 340 m / s was assumed. It is obvious that the larger antenna has a significantly extended range for unambiguous solutions based only on the shock wave. For large CPA values, the curvature difference at the two ends of the antenna (r1 and r2) is too small to be distinguishable so that the probability of a correct decision approaches 50%, or completely ambiguity. Accordingly, the confidence level depends on the size, i.e. diameter or spatial extent of the antenna.

[0030] As stated above, errors arise from timing errors and sensor coordinate uncertainty. The uncertainty of the sensor coordinates contributes to the occurrence of load errors, which are a highly variable function of the shock angle. However, for random incidence angles, sensor coordinate errors appear as random time difference errors.

[0031] Time measurement errors also arise due to both changes in gain and signal strength for individual channels. Arrival times are obtained when the sensor output values ​​rise to the set threshold value V0. Time measurement error dt due to the change in gain dg depends on the rate of voltage increase for the channel, where dg V0 dt dV / dt [0032] Time measurement errors also occur when the signal strength changes within the aperture. For aperture p of length L and a cylindrical sound source at a distance r, the maximum change in the signal level within the aperture is equal to p0 (L / 2r), where p0 is the sound pressure at the center of the aperture. The above equation for time measurement errors also applies to this type of error, with the expression

L / 2r replaces the relative gain change dg / g.

Amplitude errors are not random among the sensors, but they change evenly from the maximum through the entire aperture to zero in the middle. For ranges greater than 10 m, for an aperture of 1 m, the maximum amplitude factor is less than 0.05, which is smaller than the channel gain change parameter of 0.2 so that the effects of amplitude errors can be neglected. On the contrary, as described above, for ranges less than about 10 m, the Mach cone radius is small enough for an aperture length of 1 m so that the errors are not very significant.

[0033] Realistic estimates of time measurement errors caused by sensor uncertainty, assuming that the absolute values ​​of error vectors are statistically independent and evenly distributed between 0 and 1 mm, and that the error angles are statistically independent, the standard deviation of the equivalent evenly distributed random time difference errors will be equal

10^-3

340 · ΤΪ2 = 0.85ms.

The standard deviation of the binomial distributed random time sampling errors for system sampling at 1 MHz is equal to 0.25 μs. Time measurement errors due to gain changes are estimated to be approximately 0.75 μs for an example system with a channel bandwidth of about 18 kHz, which corresponds to a voltage change rate of about 0.02 V / gs. The acoustic sensors used in each matrix were selected to have sensitivities within ± 1.5 dB. Therefore, the relative changes in channel gain are approximately evenly distributed between 0.84 and 1.19 so that the standard deviation of the relative gain is approximately

1,19 - 0,84

ΤΪ2 = 0.10.

The threshold voltage is V0 = 0.15 V, which results in a standard deviation of time measurement errors of about 0.75 μs.

[0034] Total time measurement errors are estimated assuming that changes in channel gain, sampling changes, and sensor position changes are all statistically independent. Then the standard deviation of the time measurement error can be estimated as

-^0,85² + 0,75² + 0,25² = 1.1 ms.

[0035] It is difficult and expensive to achieve such accuracy using analogue-to-digital processing because high sampling frequencies with subsequent interpolation are necessary. The disclosed system uses two different systems to accurately measure the arrival time difference (TDOA).

[0036] In one embodiment, the exemplary system uses an analog time difference of arrival (TDOA) system that uses 1 MHz clocks in each channel. Clocks are activated when the sensor signal exceeds the threshold signal level in the reference sensor, which was defined above as the sensor that first encounters the shock wave. As discussed above, a clock frequency of 1 MHz is sufficient to eliminate in practice the importance of time sampling errors. The system works in analog mode, based on the detection of threshold levels, while the digital logic circuits perform the following functions:

1. The first flip-flop is set when the signal amplitude in the channel for the reference sensor that first encounters the shockwave exceeds the threshold.

2. The first flip-flop starts counters for each channel, which increases by one unit with each clock cycle. The processor is alerted.

3. The counter in each channel runs until the corresponding sensor encounters the shockwave. This sets a second flip-flop in the channel that stops counting on that channel. If the second flip-flop is not set, the corresponding counter runs until the upper limit value is reached.

4. The final state of each meter is saved in the TDOA digital register.

5. The processor reads the TDOA register.

6. The processor resets the counters to accept the next shockwave.

[0037] In another embodiment, the correlation for each channel is calculated with each other channel for a time period centered on the time of hardware detection.

TDOA. The correlation of two functions, designated as Corr (g, h), is defined by the relation:

Corr (g, h) ° J g (t +1) h (t) dt [0038] Correlation is a function of t, which is called "delay". Thus, it belongs to the time domain and has the following property:

Corr (g, h) G (f) H (- f) when g and h are real time functions. G (f) is the Fourier transform of the function g (t), while H (f) is the Fourier transform of the function h (t).

[0039] The total signal strength is:

+¥ +¥

Total power = J | h (t) |²dt = JH (fdf [0040] The arrival time signal has a finite length so that integration (or summation in the case of discrete data) should only be made after a specific time period, centered around the arrival time; the length of the data in one or both channels can be increased by supplementing zeros so that the duration of both signals is matched, as is known in the art.

[0041] In the discussion that follows, integrals of continuous functions have been used for simplicity, although the actual data are values ​​converted to digital and discrete values. Those skilled in the art will easily be able to replace integration by summation.

[0042] Referring now to Fig. 6, in process 60, shockwave time signal data gi (t), gj (t) are acquired on each channel i, j, steps 601, 602, and stored as a function of time. In steps 603, 604, the total signal strength in the channel is calculated and for further normalization of the correlation as:

Total channel power i ° J | g_t (T)² dt

Signal duration. [0043] Fourier transform Gi (f) shock time signal data gi (t) is calculated for the channel and a conjugate function transform Gi (-f) is created, step 605. Similarly, Fourier transform Gj (f) time data The shockwave signal gj (t) is calculated for all other channels j, step 606. Then, cross-correlation Gi (f) Gj (f) is created for each pair of channels (i, j), step 608, which is a function of fi, j (t) "delays" t.

TDOA for each pair of channels is the time tmax, where f (t) reaches its maximum value, step 610. Correlation between channels i and j can be defined as:

peak value f (t)

Corr (gg) = - ^, " ¹ yJ (Channel power i) * (Channel power j) [0044] The remainder for channel i is calculated by calculating the average value for the sensor and after all sensors j:

Rash (i) = average (Σ (1 - Corr (g_and, g_j)) ί &amp;

as given in step 612. TDOA and correlations for each channel with the best (ie, smallest) integer is then selected as the "best" solution, step 614.

[0045] As said above, channel data is typically sampled at discrete time intervals with a predetermined sampling frequency of, e.g., 41666.66 samples / sec. this corresponds to a width of 24 μs, reflecting the time resolution for the received signal. Correlation processing is performed with a time resolution that improves by a factor of 8 to 3 μs by collecting 333333 samples / 12 bek / sec.

[0046] Once the individual arrival time differences (TDOAs) between the sensors have been determined, the azimuth and angle of the shooter and the projectile trajectory can be determined based on the shockwave-only signals. Gunner positions, i.e. the distance of the gunner to the sensor array, can be determined if the muzzle blast signal is also known.

[0047] In the Cartesian coordinate system with the origin in the middle of the matrix, i.e. {(Cx0, Cy0, Cz0) = (0, 0, 0)}, the arrival time of the shock wave TOA at the given sensor (Cxj, Cyj, Cyj) (see Fig. 2) is given by:

^t Shock = ^t0 +~^{L (cos (b)} 4^{M - 1 sin (b))}

Mc where something (b) =

V_X (Xp - C) + V_s (X_p - C._s ) + V_from (X_p - Cz) LMc iV ^VUse

represents the speed of the supersonic missile

Mc = V 4V + V + V2, where c is the speed of sound and M is the Mach number. b represents the "box angle" between the shooter position and the projectile trajectory, which includes both azimuth and aiming angles. A direct hit will be b = 0. The Mach angle q defines the relationship

1 / M = sin (Q).

[0048] As said above and shown in Fig. 3, for a given shooter position and projectile trajectory, there is a different shooter position and projectile trajectory for which the shock wave TOA for a given set of sensors is almost identical. These two ambiguous solutions are actually identical if the simplified model assumes that the shock wave propagates through the sensor matrix like a flat wave. If the resolution

TDOA is high enough to solve the curvature of the shock wave, then these two almost identical solutions can be disambiguated. The fundamental ambiguity of TDOA solutions based only on shockwaves is shown in Figure 3.

[0049] Assuming sufficiently accurate TOA measurements, a true solution for shooter position and projectile trajectory can be obtained by calculating a guard / trajectory combination that minimizes the root mean square (RMS) of the remainder of the measured and calculated shockwave TDOA:

^dtmin = ^minJ Σ ^(toBL ^{- vol} measure when the sum is calculated from all sensors.

[0050] One approach to solve this problem is the Levenberg Marquardt L1 algorithm, described in detail in US Patent No. 5,930,202. The most classic point-by-point algorithms use the deterministic procedure of approaching the optimal solution, starting with a randomly guessed solution and determining the direction of the search based on a pre-determined transition rule, such as direct methods using the objective function and limit values, and gradient methods using first and second degree derivatives. However, these methods have disadvantages, e.g. such that the optimal solution depends on the initial solution chosen and such that the algorithm can "get stuck" on a suboptimal solution such as a local minimum, or when the surface of the cost function has a flat valley so that further iterations will not improve the result.

[0051] It has been found that the global minimum of shooter direction and projectile trajectory can be calculated faster and more reliably using an evolutionary genetic algorithm (GA). GA algorithms mimic natural evolutionary principles and apply them to search and optimization procedures.

[0052] A block diagram of the genetic algorithm (GA) is shown in Figure 7. Instead of starting with one conjecture solution, the GA process 70 begins its search by initiating a random population of solutions, step 71, and sets the generation counter to zero, condemning the initial set of solutions, step 72. After creating a random population of solutions, each of them is estimated in the context of a nonlinear programming problem, step 73, and each solution is assigned fitness (quality factor), step 74. Fitness can be represented as Euclidean distance At_mj_n between the calculated solution and the measured solution.

Dtmin = ^{min (vol}oBL ^{- τ}measure [0053] Intuitively, an algorithm having a low At value_min is better.

[0054] For example, when GA is used to disambiguate the solution for the shooter direction and projectile trajectory, the exemplary GA uses as the chromosome an initial population of 200 fours, with each four containing the following values:

[Azimuth Shot, Aiming Angle, Shot, Azimuth, Missing, Aiming Angle] [0055] [Azimuth, Shot, Aiming Shot] are defined by the angle (θ + β), while [Azimuth, Missing, Angle of Missing] are determined by the angle β (see Fig. 3). Since the muzzle blast is not used in the shock wave-based approach described above, a nominal distance of 100 m between the sensor array and the shooter is assumed.

[0056] The initial population is created by randomly selecting fours covering a significant and reasonable range of values ​​(all values ​​in degrees):

AzymutShooter = {0, ..., 360},

Aiming angle Shot = {-10, ..., 30},

AzymutChybienia = {-20, ..., 20},

Aiming angle: Missing = {-20, .... 20}.

[0057] In step 75, it is checked if the maximum number of iterations for GA, set e.g. at 25., has been reached. If the maximum number of iterations has been reached, the process 70 stops in step 80 and either the result can be accepted or further evaluated. Otherwise, step 76 checks if the adopted fitness criteria have been met. [0058] Fitness criteria may be, e.g., a calculated missed azimuth of &lt;15 ° and / or a ratio of the residues of two ambiguous solutions. If the fitness criteria are met, the process 70 stops at step 80; otherwise, a new population is created by crossing, step 77, and mutation, step 78, and the generation count increases by one, step 79.

[0059] In each generation, the "best" individual can survive unchanged, however

100 the best individuals, assessed according to their fitness criteria, are also surviving, but they are used to create the next 100 individuals from survivors pairs, using the crossover / mutation operators listed in Table 1.

[0060] The following exemplary crossover and mutation operators were used to demonstrate process 70:

Table 1

<td>Operator Name</td><td>Operator Type</td><td>Probability</td><td>Description</td>

<td>Azimuth crossing</td><td>crossing</td><td> 0,5</td><td>Exchanges shooter trajectory / azimuth between two chromosomes</td>

<td>Missing crosses</td><td>crossing</td><td> 0,5</td><td>Exchange miss azimuth / miss targeting angle between two chromosomes</td>

<td>Field mutation</td><td>mutations</td><td> 0,3</td><td>Replaces the given field (with a probability of 0.25 for the field) with a randomly selected new value in the range</td>

<td>Incremental mutation</td><td>mutations</td><td> 0,4</td><td>Introduces small mutations in all chromosome fields (within &lt;2 ° for shooter information; within &lt;0.5 ° for misinformation</td>

<td>Switching mutation</td><td>mutations</td><td> 0,1</td><td>Changes the solution to an ambiguous alternative solution</td>

<td>No mutation</td><td>mutations</td><td> 0,2</td><td>The chromosome stays intact</td>

[0061] Disambiguation is obtained and / or improved by performing gradient search at best and in the appropriate alternative. Residues and residue ratios are calculated for both ambiguous solutions. If the calculated missed azimuth is &lt;15 °, representing "close" arrows, and if the residual ratio is&gt; 2, then the solution with the smallest residue is chosen. Otherwise, no real choice is made and the solution with a smaller amount is designated as the "main" solution, while the other solution is designated as the "alternative" solution. [0062] In the case of shockwave only detection, the GA algorithm generates a solution using a computer with a 1 GHz clock operating under the operating system

Linux, within 0.15 s in a wide range of simulated shots. 97% of simulated shots were within 15 ° miss azimuth, and 86% of simulated shots were within 5 ° miss azimuth. Using the disambiguation algorithm described above, close shots, i.e. having missed azimuth &lt;15 °, were determined in 95% of cases. The disambiguation algorithm generated correct results for more distant shots in 70% of cases. It is expected that the accuracy of disambiguation will change depending on the sensor matrix geometry and the assumed distribution of shots, with shots having a small aiming angle easier to disambiguate.

[0063] The solutions described above for the projectile trajectory were obtained without detecting muzzle blast. However, it has been found that even a weak muzzle signal or muzzle signal received only on a limited number of channels can be advantageously used to improve distance determination and disambiguation.

[0064] Fig. 9 schematically shows the arrival time model (TOA), which is described in more detail in US Patent 6,178141. The TOA model can be used to estimate the projectile trajectory and shooter direction relative to the location of the sensor. TOA model is based on a ballistic model, taking into account certain physical parameters related to the projectile path, such as air density (which is related to temperature), P position (Px, Py, Pz) of the shooter, azimuth and angle of aiming of the rifle barrel, speed bullet muzzle (or equivalent Mach number) and sound velocity (which varies with air temperature / density). With this ballistic model, you can accurately calculate the time after which the shock wave and muzzle blast reach a specific point in space.

[0065] As shown in the diagram of Fig. 9, the shooter is located at the point P (PX, PY, PZ) relative to the beginning of the system (0, 0, 0), the individual sensors are located at the points S _j (Sxj, Syj, Szj), and the projectile trajectory is shown as running from the shooter towards

A. The vector from the shooter to the jth sensor is denoted by D, the nearest point of approach (CPA) of the projectile to the jth sensor is | R | = | D | sin (e), while the path from the point from which the shock wave is propagating from the trajectory to the jth sensor is marked as S (j sensor indicator is omitted). The mach angle of the projectile is equal to θ = sin^-1 (1 / M), M =

Vlc₀. M is the mach number of the projectile, V is the supersonic speed of the projectile, ac₀ is the speed of sound (depends on pressure and temperature). The "miss angle" between the trajectory and the jth sensor is designated β. The trajectory is characterized by its azimuth angle, measured counterclockwise from the x axis in the x - y plane, and the aiming angle measured up from the x - y plane. Equations defining the time of arrival of the shock wave, i.e. the unit vector at the j-th sensor recorded as a function of these geometrical quantities.

[0066] The arrival time is equal to the time the missile needs to travel the distance A to the point from which the sound propagated towards the j-sensor, plus the distance S from the point of propagation to the j-sensor.

WND. "^tj = ^t0 + v + - “^this + - sin (^ + 0),

V c0 c 0 where t0 is the time reference time (shot moment) and c0 is the speed of sound. The Mach angle q is also shown in Figure 9.

[0067] It can be safely assumed that the velocity V of the bullet remains constant in the path corresponding to the spacing of the sensors, such that there is a slight decrease in speed between the moments when the bullet sends a wave towards the individual sensors. However, it is known that at greater distances, the projectiles slow down due to air resistance. Air resistance can be expressed by the front drag coefficient Cb, which depends on the shape and caliber of the projectile. A mathematical ballistic model based on the laws of physics can predict the arrival time of a shockwave to any general point in space as a function of a full set of parameters describing a projectile (e.g., through its frontal resistance coefficient Cb), and its initial velocity and density of the surrounding air are known in advance.

[0068] The parameters required for accurate calculation are usually unknown in a real environment, such as e.g. a battlefield. However, distance estimation can be significantly improved by using the iterative process 200 shown in the block diagram in Figure 10, which takes into account the decrease in projectile speed along the trajectory. Process 200 begins with step 202 with the following assumptions:

c0 = sound speed modified for ambient temperature / air pressure (-340 m / s)

Cb = nominal frontal resistance coefficient averaged for the anticipated weapons

V0 = initial velocity of the bullet, at the time of firing, averaged for the anticipated weapons

M0 = V0 / c = initial Mach number of the bullet [0069] The first estimation of the shooter distance D0 is calculated in step 204 using the measured arrival time difference (TDOA) tms and the angle of incidence α between the shock wave and muzzle sound on the sensor matrix and assuming the initial constant V0 speed by Mach number M0, according to the equation D = - ^tms'^C

1- cos (a) [0070] Under these assumptions, the speed of the bullet at a distance a from the position of the shooter P can be calculated in step 206 from the equation:

Ma = M0

1V CJO,_} [0071] so that the time the bullet travels the distance along the trajectory, step 208, is expressed as:

T = a

V0^{and V}0 C [0072] The angle q is related to the Mach number Ma with the equation:

^{sin (Q)} = ~^~ ^Mand in which Mach number Ma initially receives the value M0. It should be noted that the instant bullet speed is taken as the speed of sound (i.e. Ma = 1) if the calculated bullet speed becomes less than the speed of sound. The corrected distance a = in step 210 then takes the value:

a = D- cos (b)

1tg (b) • Jm; [0073] Angles α, β and θ are related by the equation (α + β + θ) = 90 °. The process 200 then returns in a loop to step 206 inserting the calculated distance value a into the equation for Ma and Ta above, resulting in a correspondingly updated Mach number and updated projectile time of the projectile Ta, for the distance traveled. Measured TDOA tms and calculated updated values ​​of Ta and a is then used for the next update of the D value of the shooter distance:

D = c 0 - (tms + Ta) + s [0074] This process is repeated until either the maximum number of iterations is reached or the range of D values ​​begins to converge, as determined in step 212.

[0075] Process 200 also checks in step 214 whether the updated value D = distance between shooter and sensor matrix is ​​a "reasonable" value, in which case process 200 ends in step 216. For example, the value of D can be considered valid if the distance a traveled by the projectile and the distance s = a sin (b) between the sensor and the point whose sound wave is radiated from the projectile to the sensor are valid numerical values, i.e. not NAN. NAN is a special floating point value that represents the result of a calculation operation that cannot return a valid numeric value and is usually used to prevent error propagation in the calculation process. In addition, both angles α and β should be smaller than the set threshold indicating that the projectile was actually launched towards the sensor matrix.

[0076] As stated above, a pair of numbers (tms, α) are initially used to calculate the shooter distance D0 in zero approximation, disregarding the change in projectile velocity along the trajectory. If the iterative process 200 described above does not return consistent geometry supporting a pair of numbers (tms, α), then the solution is rejected.

[0077] Even if it may not be possible to obtain an accurate solution, the goal is to find values ​​for shooter distance D and misaligned azimuth and aiming angles (which are associated with β) that most closely match the measured shock wave TDOA and the measured exhaust sound TDOA. As has already been said, based only on the TDOA shock wave between individual sensors can be reliably measured in most situations. Gunner azimuth and gunner aiming angle, but not shooter distance, can be determined on the basis of TDOA only for the shock wave, using known sensor matrix coordinates (Sxj, Syj, Szj). It will be assumed that TDOA tms between the detected shock wave and muzzle sound can also be measured, so that muzzle sound may not be detected by all sensors.

[0078] If it is determined in step 214 that the iterative process 200 does not return a valid result, then the process 200 attempts to calculate the shooter distance by referring to the evolutionary genetic algorithm (GA) 300. GA algorithms mimic natural evolutionary principles and apply them to search procedures and optimization. GA begins its search with a random set of solutions instead of just one solution. When a random population of solutions is created, each solution is evaluated in the context of nonlinear programming, and each solution is assigned fitness (quality factor). In one embodiment, the fitness can be represented by the Euclidean distance between the calculated solution and the measured solution, e.g. by:

Dt =

Σ^(t j = 1 shock, calculated ^-tshock, measure ⁾ + abs ms, calculated

V ^M J

-t ms, measure [0079] Intuitively, an algorithm that generates smaller At values_mj_n is better.

[0080] Fig. 11 is a block diagram of a process 300 using GA 300. Process 300 uses the arrival time difference (TDOA) tms and angle of incidence α previously measured for process 200, step 302. An example of the number of triples having {RANGE values, MA, ME} is defined as the initial population, step 304, where

RANGE is the shooter distance D = D shown in Fig. 9, MA is the missed azimuth, and ME is the missed aim angle. MA and ME values ​​indicate how much the missile missed the target in azimuth and aiming angle. It is assumed that the target in the example shown is a sensor matrix. The initial population in step 304 is created by randomly selecting triples covering a significant and reasonable range of values:

Distance Shot = {1000, ..., 3000} [meters]

AzymutChbiebie = {-20, ..., 20} [degrees],

Aiming angle: Missing = {-20, 20} [degrees].

[0081] The calculations follow a similar process to the one outlined earlier for a solution based only on shockwaves. Initially, for the Gen = 0 generation, step 306, the shooter position vector P (Px, Py, Pz) is calculated for each three with a predetermined azimuth and angle of aim of the shooter, and with the assumed^r

RANGE for a specific three. Assuming the initial Mach number M0, the vector A (Ax, Ay, Az), i.e. the position from which the sound propagates, is calculated with the values ​​MA and

ME for each trio, step 308. Distance D = S j - P between the shooter and each sensor j that detects a shock wave is also calculated.

[0082] For each three, the angle β is calculated by the equation:

A • (P - S) something (b)

PS in which the "•" symbol represents the scalar product of two vectors. The updated values ​​of distance a, flight time Ta of the projectile at distance a, and Mach number are calculated by inserting the calculated value β and the initially assumed values ​​of Ma = M0 and a into the above equations for Ma, Ta, α and D, step 312. This process iterates a number of times for each of the three, e.g. three times, as determined in step 312, after which the At ™ residue as defined above which contains the muzzle signal, step 314, is calculated for each three.

[0083] In step 316, it is checked if the maximum number of iterations for GA, e.g. 25 iterations, has been reached. If the maximum number of iterations has been reached, then the process 300 stops at step 320, returning the three with the smallest residue. Otherwise, process 300 creates a new population through crossbreeding and mutation operations, step 318 and generation counter increases by one, step 322.

[0084] In each generation, the "best" individual can survive unchanged, however

100 the best individuals, assessed according to their fitness criteria, are also surviving, but they are used to create another 100 individuals from survivors pairs, using the crossover / mutation operators listed in Table 2 below. [0085] The following crossover and mutation operators were used to demonstrate process 300:

Table 2

<td>Operator Name</td><td>Operator Type</td><td>Probability</td><td>Description</td>

<td>Crossover missed azimuth</td><td>crossing</td><td> 0,5</td><td>Exchanges miss azimuth between two chromosomes</td>

<td>Crossing miss angle</td><td>crossing</td><td> 0,5</td><td>Lists the miss targeting angle between two chromosomes</td>

<td>Crossing miss distance</td><td>crossing</td><td> 0,5</td><td>Lists the miss distance between two chromosomes</td>

<td>Field mutation</td><td>mutations</td><td> 0,3</td><td>Replaces the given field (with a probability of 0.25 for the field) with a randomly selected new value in the range</td>

<td>Incremental mutation</td><td>mutations</td><td> 0,4</td><td>Introduces small interference in all chromosome fields (within ± 2 m for shooter range; within ± 0.1 ° for missed azimuth and aiming angle)</td>

[0086] The GA 300 process is performed with an initial population of 200 different triples, with a degree of replenishment of 50, for a total of 25 generations. GA is performed 5 times in parallel with different sets of initial triples, and the solution with the smallest remainder is chosen as the final solution for RANGE, miss azimuth and shooter miss angle, which allows calculating the vector D.

[0087] Recent experimental studies have shown a reduction in ambiguous shots from 95% to 8% with the same data set as a result of using at least one muzzle channel in addition to five or more shockwave channels, which is a significant improvement over solutions based only on on the shockwave.

[0088] Calculations that do not take into account the decrease in the projectile speed along its length due to air resistance tend to overstate. For some geometries and far enough shots, this overstatement may exceed 20%. The process described above eliminates this load from distance estimation when detecting distant shots.

[0089] As described above, ambiguous solutions based only on shockwaves can often be disambiguated by comparing residuals from two different trajectories and selecting a trajectory with a smaller residue.

[0090] If muzzle blast signals are detected in 4 or more sensor channels, then the shock wave - muzzle signal algorithms described above can be used to uniquely determine the position of the shooter, regardless of the number of channels associated with the shock wave. If muzzle blast signals are detected in less than 4 sensors, but the shockwave signals are detected in 5 or more shockwave channels, then the above GA can be used with the modified cost function or the rest, so that any available muzzle signals are "mixed in" to optimization function to disambiguate the solution based only on the shock wave and / or to refine the shooter distance estimation. However, if less than 3 muzzle and less than 5 blast channels are detected, then the alarm can be triggered without attempting to locate the shooter.

[0091] The muzzle signal may not be detected reliably on all channels because:

1. The detection level on one or more channels is too low to detect reliably.

2. You can't tell the difference between muzzle energy in a raw signal, which causes the system to associate it with "noise", giving uncertain TDOA estimates.

3. The shock wave echoes can be stronger than the muzzle blast and can arrive earlier than the muzzle blast, which causes the system to mistakenly detect the shock wave as an muzzle blast.

[0092] If only muzzle blast signal is detected on some channels, the rest in this case can be defined as

Dt =, Σ (D «muzzle, calculate ^dtmuzzle, measure)²+ ς ⁽dt

Shock calculated _dt^j

Shock, measure where the first component for the muzzle blast is added after the reduced number of sensors (&lt;4) that detect the muzzle blast, and j is added after the sensors detecting the shock wave (these are usually all sensors).

[0093] As demonstrated by the examples described above, the muzzle blast signal provides important information about the azimuth of the shooter and thus about the projectile trajectory, compared to solutions based only on the shock wave, so that the calculated trajectory solution approaches one of the ambiguous solutions more closely , i.e.

thus making the solutions more clear.

[0094] Without at least some certain muzzle blast signals, a significant number of ambiguous solutions based on shockwaves only can be generated, especially at large shooter distances, which is less desirable than a small number of unambiguous but less accurate solutions.

[0095] In the case of potentially uncertain muzzle blast detection, one can still initially attempt to detect muzzle signals, e.g., find the muzzle blast feature in the noisy signal and calculate the resulting TDOA. Detection of muzzle signals can be considered reliable if the muzzle signals are detected by a sufficient number of sensors with sufficient cross-correlation between channels u if there is a sufficiently strong correlation between the muzzle signal with an appropriate primary band on each channel (shifted by a certain number of frequency intervals to account for delays filters).

[0096] Otherwise, at least muzzle signals having inadequate correlation are deleted and reference is made to the following "coarse muzzle signal" logic:

- Search for peaks in shockwave energy by shockwave. Mark these peaks as likely "shock wave echoes", thus excluding them as muzzle blasts.

- Specify the maximum time the muzzle wave would need to pass through the sensor array and specify a "window" having the appropriate duration. Look for muzzle energy peaks by moving this window through essentially all detector channels following the detected shock wave, skipping sections in the detected signals that have been identified as shock wave echoes. Integrate energy around the window, i.e. look for the maximum:

f (i) = ΣΣ (mbi + j)² n = 0 j = 0 where square mb_andⁿ+ j represents the amount of energy, e.g. muzzle blast energy, measured by the nth sensor. (i + j) indicates the detection channel, wherein i denotes a discrete time interval between the moment the shock wave is detected and the beginning of the window, and j represents the time interval from the beginning of the window.

[0097] To distinguish from noise, the peak energy in the window generating the maximum of the function fmax (i) is checked to determine if the energy peak at this maximum is greater than the energy in the window area, in a given aspect ratio. If so, then the signal in the window is identified as muzzle blast detection and the cross-correlation is determined on all channels in the muzzle blast band to determine the muzzle wave TDOA.

[0098] The detected muzzle blast signal can then be used to determine the shooter distance and / or disambiguate the shockwave signal as described above.

In summary, the described system can accurately, quickly and often unambiguously provide shooter direction and projectile trajectory based only on shockwave measurements. Disambiguation can be improved and shooter distance can be estimated as soon as a weak muzzle blast is also detected. The system is relatively insensitive to false shooter signals in response to vehicle vibration and noise, wind noise, firecrackers or close shots fired away from the system.

[0100] It should be noted that the shockwave detection system performs two tests on the initial waveforms to determine if the signal can actually be attributed to the shockwaves. First, the measured total energy in the frequency band between about 700 Hz and 10 kHz is compared to an empirical threshold. Only if this threshold is exceeded can the signal waveform be derived from the shockwave. Secondly, the time interval of the detected initial positive pressure peak must be greater than about 70 μs and less than about 300 μs. These criteria ensure the system's resistance to impulsive noises such as firecrackers and harmless shooting. If these tests do not give positive results, the detected wave is not considered a shock wave and there is no attempt to obtain a solution for the shooter.

[0101] For example, although illustrative embodiments are illustrated with acoustic sensors such as microphones, this need not be the case. Instead, other types of pressure-sensitive mechanical or electrical sensors can be used. In addition, the values ​​given in Tables 1 and 2 for different operators are intended to be examples, and other values ​​may be selected depending on actual field conditions.